curve_fitから信頼区間を取得する方法

Question

私の質問には統計とpythonが関係しており、私は両方の初心者です。私はシミュレーションを実行しており、独立変数（X）の各値に対して、従属変数（Y）の1000個の値を生成します。私がしたことは、Xの各値に対してYの平均を計算し、これらの平均をscipy.optimize.curve_fitを使って当てはめたことです。曲線はきれいにフィットしますが、私も信頼区間を描きたいと思います。私がやっていることが正しいのか、私がやりたいことができるのか分かりませんが、私の質問は、curve_fitによって作られた共分散行列から信頼区間を得る方法です。コードはまずファイルから平均値を読み込み、その後単にcurve_fitを使用します。

import numpy as np 
import matplotlib.pyplot as plt 
from scipy.optimize import curve_fit 


def readTDvsTx(L, B, P, fileformat): 
    # L should be '_Fixed_' or '_' 
    TD = [] 
    infile = open(fileformat.format(L, B, P), 'r') 
    infile.readline() # To remove header 
    for line in infile: 
     l = line.split() # each line contains TxR followed by CD followed by TD 
     if eval(l[0]) >= 70 and eval(l[0]) <=190: 
      td = eval(l[2]) 
      TD.append(td) 
    infile.close() 
    tdArray = np.array(TD) 

    return tdArray 


def rec(x, a, b): 
    return a * (1/(x**2)) + b 



fileformat = 'Densities_file{}BS{}_PRNTS{}.txt' 
txR = np.array(range(70, 200, 20)) 
parents = np.array(range(1,6)) 
disc_p1 = readTDvsTx('_Fixed_', 5, 1, fileformat) 


popt, pcov = curve_fit(rec, txR, disc_p1) 


plt.plot(txR, rec(txR, popt[0], popt[1]), 'r-') 
plt.plot(txR, disc_p1, '.') 

print(popt) 
plt.show() 

そして、ここで結果フィットさ：

Answer 1

は、ここでは簡単と間違った答えです：あなたは、その対角線の平方根としてごaとbのパラメータの共分散行列からエラーを近似することができます。 np.sqrt(np.diagonal(pcov))。パラメータの不確かさを使用して信頼区間を描くことができます。

データをモデルに適合させる前に、平均したdisc_p1ポイントの誤差の見積もりが必要になるため、答えが間違っています。平均化すると、母集団の散布に関する情報が失われ、curve_fitには、それを供給しているy点が絶対的ではないと信じています。これは、パラメータエラーの過小評価を引き起こす可能性があります。

平均Y値の不確実性を推定するには、誤差の絶対値が「curve_fit」になるまで分散尺度を推定し、それを渡す必要があります。以下は、各点が正規分布から描かれた1000個のサンプルからなるランダムなデータセットに対してこれを行う方法の例です。

from scipy.optimize import curve_fit 
import matplotlib.pylab as plt 
import numpy as np 

# model function 
func = lambda x, a, b: a * (1/(x**2)) + b 

# approximating OP points 
n_ypoints = 7 
x_data = np.linspace(70, 190, n_ypoints) 

# approximating the original scatter in Y-data 
n_nested_points = 1000 
point_errors = 50 
y_data = [func(x, 4e6, -100) + np.random.normal(x, point_errors, 
      n_nested_points) for x in x_data] 

# averages and dispersion of data 
y_means = np.array(y_data).mean(axis = 1) 
y_spread = np.array(y_data).std(axis = 1) 

best_fit_ab, covar = curve_fit(func, x_data, y_means, 
           sigma = y_spread, 
           absolute_sigma = True) 
sigma_ab = np.sqrt(np.diagonal(covar)) 

from uncertainties import ufloat 
a = ufloat(best_fit_ab[0], sigma_ab[0]) 
b = ufloat(best_fit_ab[1], sigma_ab[1]) 
text_res = "Best fit parameters:\na = {}\nb = {}".format(a, b) 
print(text_res) 

# plotting the unaveraged data 
flier_kwargs = dict(marker = 'o', markerfacecolor = 'silver', 
        markersize = 3, alpha=0.7) 
line_kwargs = dict(color = 'k', linewidth = 1) 
bp = plt.boxplot(y_data, positions = x_data, 
       capprops = line_kwargs, 
       boxprops = line_kwargs, 
       whiskerprops = line_kwargs, 
       medianprops = line_kwargs, 
       flierprops = flier_kwargs, 
       widths = 5, 
       manage_xticks = False) 
# plotting the averaged data with calculated dispersion 
#plt.scatter(x_data, y_means, facecolor = 'silver', alpha = 1) 
#plt.errorbar(x_data, y_means, y_spread, fmt = 'none', ecolor = 'black') 

# plotting the model 
hires_x = np.linspace(50, 190, 100) 
plt.plot(hires_x, func(hires_x, *best_fit_ab), 'black') 
bound_upper = func(hires_x, *(best_fit_ab + sigma_ab)) 
bound_lower = func(hires_x, *(best_fit_ab - sigma_ab)) 
# plotting the confidence intervals 
plt.fill_between(hires_x, bound_lower, bound_upper, 
       color = 'black', alpha = 0.15) 
plt.text(140, 800, text_res) 
plt.xlim(40, 200) 
plt.ylim(0, 1000) 
plt.show() 

編集： あなたがデータポイント上の本質的な誤差を考慮していない場合、あなたはおそらく、私は前に述べた「qiuckと間違った」ケースを使用して大丈夫です。次に、共分散行列の対角成分の平方根を使用して、信頼区間を計算することができます。

from scipy.optimize import curve_fit 
import matplotlib.pylab as plt 
import numpy as np 

func = lambda x, a, b: a * (1/(x**2)) + b 

n_ypoints = 7 
x_data = np.linspace(70, 190, n_ypoints) 

y_data = np.array([786.31, 487.27, 341.78, 265.49, 
        224.76, 208.04, 200.22]) 
best_fit_ab, covar = curve_fit(func, x_data, y_data) 
sigma_ab = np.sqrt(np.diagonal(covar)) 

# an easy way to properly format parameter errors 
from uncertainties import ufloat 
a = ufloat(best_fit_ab[0], sigma_ab[0]) 
b = ufloat(best_fit_ab[1], sigma_ab[1]) 
text_res = "Best fit parameters:\na = {}\nb = {}".format(a, b) 
print(text_res) 

plt.scatter(x_data, y_data, facecolor = 'silver', 
      edgecolor = 'k', s = 10, alpha = 1) 

# plotting the model 
hires_x = np.linspace(50, 200, 100) 
plt.plot(hires_x, func(hires_x, *best_fit_ab), 'black') 
bound_upper = func(hires_x, *(best_fit_ab + sigma_ab)) 
bound_lower = func(hires_x, *(best_fit_ab - sigma_ab)) 
# plotting the confidence intervals 
plt.fill_between(hires_x, bound_lower, bound_upper, 
       color = 'black', alpha = 0.15) 
plt.text(140, 630, text_res) 
plt.xlim(60, 200) 
plt.ylim(0, 800) 
plt.show() 

あなたは絶対誤差を含めるかどうか、あなたにそれらを推定する方法がわからない場合：しかし、信頼区間は、我々は不確実性を落としたことになりましたを縮小していることに注意してください場合は、Cross Validatedでアドバイスを求める方がよいでしょう。スタックオーバーフローは、主に回帰方法の実装に関する議論であり、基礎となる統計に関する議論ではないからです。