12

だから私はいくつかのデータを2つのリストとして保存し、それらを使用してプロットしました

plot(datasetx, datasety)

次に、トレンドラインを設定します

trend = polyfit(datasetx, datasety)
trendx = []
trendy = []

for a in range(datasetx[0], (datasetx[-1]+1)):
    trendx.append(a)
    trendy.append(trend[0]*a**2 + trend[1]*a + trend[2])

plot(trendx, trendy)

しかし、元のデータセットのエラーである 3 番目のデータ リストがあります。エラーバーをプロットすることには問題ありませんが、これを使用して多項式トレンドラインの係数のエラーを見つける方法がわからないのです。

したがって、私のトレンドラインが 5x^2 + 3x + 4 = y であることが判明したとします。5、3、および 4 の値には何らかのエラーが必要です。

これを計算する NumPy を使用するツールはありますか?

4

2 に答える 2

14

ドキュメントcurve_fit)の機能が使えると思います。使用法の基本的な例:scipy.optimize

import numpy as np
from scipy.optimize import curve_fit

def func(x, a, b, c):
    return a*x**2 + b*x + c

x = np.linspace(0,4,50)
y = func(x, 5, 3, 4)
yn = y + 0.2*np.random.normal(size=len(x))

popt, pcov = curve_fit(func, x, yn)

ドキュメントに続いて、pcovは次のようになります。

poptの推定共分散。対角線は、パラメーター推定値の分散を提供します。

したがって、このようにして、係数の誤差推定を計算できます。標準偏差を得るには、分散の平方根を取ることができます。

これで係数にエラーが発生しましたが、これはydataと近似の間の偏差にのみ基づいています。ydata自体のエラーも考慮したい場合は、curve_fit関数は次のsigma引数を提供します。

シグマ:なしまたはN長シーケンス

Noneでない場合は、ydataの標準偏差を表します。このベクトルが与えられた場合、最小二乗問題の重みとして使用されます。

完全な例:

import numpy as np
from scipy.optimize import curve_fit

def func(x, a, b, c):
    return a*x**2 + b*x + c

x = np.linspace(0,4,20)
y = func(x, 5, 3, 4)
# generate noisy ydata
yn = y + 0.2 * y * np.random.normal(size=len(x))
# generate error on ydata
y_sigma = 0.2 * y * np.random.normal(size=len(x))

popt, pcov = curve_fit(func, x, yn, sigma = y_sigma)

# plot
import matplotlib.pyplot as plt

fig = plt.figure()
ax = fig.add_subplot(111)
ax.errorbar(x, yn, yerr = y_sigma, fmt = 'o')
ax.plot(x, np.polyval(popt, x), '-')
ax.text(0.5, 100, r"a = {0:.3f} +/- {1:.3f}".format(popt[0], pcov[0,0]**0.5))
ax.text(0.5, 90, r"b = {0:.3f} +/- {1:.3f}".format(popt[1], pcov[1,1]**0.5))
ax.text(0.5, 80, r"c = {0:.3f} +/- {1:.3f}".format(popt[2], pcov[2,2]**0.5))
ax.grid()
plt.show()

結果


次に、numpy配列の使用について別のことを説明します。numpyを使用する主な利点の1つは、配列に対する操作が要素ごとに適用されるため、forループを回避できることです。したがって、例のforループは次のように実行することもできます。

trendx = arange(datasetx[0], (datasetx[-1]+1))
trendy = trend[0]*trendx**2 + trend[1]*trendx + trend[2]

arangeリストの代わりにnumpy配列を返すため、範囲の代わりに使用します。この場合、numpy関数を使用することもできますpolyval

trendy = polyval(trend, trendx)
于 2011-08-25T09:03:55.433 に答える
1

numpy または python に組み込まれている係数でエラーを取得する方法を見つけることができませんでした。John Taylor のAn Introduction to Error Analysisのセクション 8.5 および 8.6 に基づいて作成した単純なツールがあります。あなたのタスクにはこれで十分かもしれません (デフォルトの戻り値は標準偏差ではなく分散であることに注意してください)。共分散が大きいため、(提供されている例のように) 大きなエラーが発生する可能性があります。

def leastSquares(xMat, yMat):
'''
Purpose
-------
Perform least squares using the procedure outlined in 8.5 and 8.6 of Taylor, solving
matrix equation X a = Y

Examples
--------
>>> from scipy import matrix
>>> xMat = matrix([[  1,   5,  25],
                   [  1,   7,  49],
                   [  1,   9,  81],
                   [  1,  11, 121]])
>>> # matrix has rows of format [constant, x, x^2]
>>> yMat = matrix([[142],
                   [168],
                   [211],
                   [251]])
>>> a, varCoef, yRes = leastSquares(xMat, yMat)
>>> # a is a column matrix, holding the three coefficients a, b, c, corresponding to
>>> # the equation a + b*x + c*x^2

Returns
-------
a: matrix
    best fit coefficients
varCoef: matrix
    variance of derived coefficents
yRes: matrix
    y-residuals of fit 
'''
xMatSize = xMat.shape
numMeas = xMatSize[0]
numVars = xMatSize[1]

xxMat = xMat.T * xMat
xyMat = xMat.T * yMat
xxMatI = xxMat.I

aMat = xxMatI * xyMat
yAvgMat = xMat * aMat
yRes = yMat - yAvgMat

var = (yRes.T * yRes) / (numMeas - numVars)
varCoef = xxMatI.diagonal() * var[0, 0]

return aMat, varCoef, yRes
于 2011-08-25T06:27:38.450 に答える