Как найти ошибку в параметрах для (хи-квадрат) кривой, подходящей для 5-го порядка? - PullRequest
Новая
       15

Как найти ошибку в параметрах для (хи-квадрат) кривой, подходящей для 5-го порядка?

0 голосов
/ 19 ноября 2018

Проясняю мой вопрос:

Я определил минимизированное значение хи-квадрат, назовем его def chisq (), для некоторых данных с учетом подбора, включающего 5 параметров: a, b, c, d и e.

Я хотел бы знать о функции или способе, где я могу минимизировать свою функцию хи-квадрат, сохраняя 1 параметр постоянным, в то же время изменяя другие так, чтобы стремиться к минимальному значению ... И

Способ, которым я могу сохранить 4 из 5 параметров постоянными, минимизируя до заданного значения (значение будет chisq () + 1, чтобы иметь тенденцию к экстремуму контура графика квадрата хи, но это просто некоторый добавленный контекст).

Так что я застрял ... В любом случае, я надеюсь, что это более логично, чем моя последняя попытка

Обновление: это то, где я сейчас нахожусь. Я позаимствовал основной код с другой страницы и пытаюсь адаптировать его под свои нужды. И я никогда не собирался делать участки уверенности, но они выглядят очень хорошо.

Однако я получаю сообщение об ошибке «Объект numpy.ndarray не может быть вызван». Я изучил это и думаю, что это связано с тем, что я ссылаюсь на массив, как на функцию. Это приводит меня к разделу минимизатора / мини, в котором, я думаю, и заключается проблема.

Так что любая помощь в том, как я могу исправить и адаптировать этот код, будет высоко ценится. Спасибо, пока @mikuszefski 

import numpy as np
import matplotlib.pyplot as plt
import lmfit


# Data

xval1 = np.array([0.11914471, 0.230193321, 0.346926427, 0.460501481, 0.575512013, 0.689201904, 0.804958885, 0.918017166, 1.034635434, 1.149990481, 1.266264234]) 
xval2 = np.array([0.116503429, 0.230078483, 0.348247067, 0.461420187, 0.577751359, 0.690120611, 0.80398276, 0.918878453, 1.033716728, 1.150794349, 1.266149395]) 
xval3 = np.array([0.115240208, 0.230710093, 0.346007721, 0.458032458, 0.576028785, 0.692359957, 0.80725565, 0.920888123, 1.035037368, 1.147521458, 1.267871969])
xval_av = (xval1 + xval2 + xval3 )/3

yval1 = np.array([2173.446136, 2400.969972, 2547.130603, 2658.022565, 2723.098769, 2760.481961, 2761.881166, 2734.638209, 2671.839502, 2559.251885, 2313.774753])
yval2 = np.array([2185.207051, 2441.547714, 2587.120886, 2701.697704, 2765.897692, 2792.190609, 2793.030024, 2761.191402, 2697.078931, 2583.572776, 2329.860358])
yval3 = np.array([2178.269249, 2438.81575, 2601.305985, 2704.228832, 2765.493933, 2802.801105, 2796.873214, 2760.426641, 2690.92674, 2575.961498, 2330.099396])
yval_av = (yval1 + yval2 + yval3 )/3

yerr1 = np.array([28.0041333, 17.22884875, 11.77504404, 12.0784658, 11.43116392, 10.17734322, 9.839420192, 7.879393332, 8.586145049, 8.129879332, 6.543378761])
yerr2 = np.array([28.84371257, 22.30391049, 16.3503791, 12.32503283, 10.6931908, 8.909895306, 9.310031644, 9.619524491, 7.725528299, 6.64584248, 5.483917187])
yerr3 = np.array([25.22103486, 17.58553765, 16.56186149, 14.93545788, 9.884470694, 11.15591573, 8.696907113, 9.602409632, 8.156617355, 6.966615987, 5.706426531])
yerr_av = (yerr1 + yerr2 + yerr3 )/3

# Parameters

a_original = 1772.73301389 
b_original = 4137.16888269 
c_original = -6903.58620042 
d_original = 5845.15206084 

#-----------------------------------------------------------------------------------------

x = xval_av
y = yval_av

p = lmfit.Parameters()
p.add_many(('a', a_original), ('b', b_original), ('c', c_original), ('d', d_original))

def residual(x,y,p):
    return p['a'] + p['b']*x + p['c']*x**2 + p['d']*x**3 - y 

# Minimiser
mini = lmfit.Minimizer(residual(x,y,p), p, nan_policy='omit')

# Nelder-Mead
out1 = mini.minimize(method='Nelder')

# Levenberg-Marquardt
# Nelder-Mead as initial guess
out2 = mini.minimize(method='leastsq', params=out1.params)

lmfit.report_fit(out2.params, min_correl=0.5)

ci, trace = lmfit.conf_interval(mini, out2, sigmas=[1, 2],
                                trace=True, verbose=False)
lmfit.printfuncs.report_ci(ci)

plot_type = 2

if plot_type == 0:
    plt.plot(x, y)
    plt.plot(x, residual(out2.params) + y)

elif plot_type == 1:
    cx, cy, grid = lmfit.conf_interval2d(mini, out2, 'b', 'd' )
    plt.contourf(cx, cy, grid, np.linspace(0, 1, 11))
    plt.xlabel('b')
    plt.colorbar()
    plt.ylabel('d')

elif plot_type == 2:
    cx, cy, grid = lmfit.conf_interval2d(mini, out2, 'a', 'd', 30, 30)
    plt.contourf(cx, cy, grid, np.linspace(0, 1, 11))
    plt.xlabel('a')
    plt.colorbar()
    plt.ylabel('d')

elif plot_type == 3:
    cx1, cy1, prob = trace['a']['a'], trace['a']['d'], trace['a']['prob']
    cx2, cy2, prob2 = trace['d']['d'], trace['d']['a'], trace['d']['prob']
    plt.scatter(cx1, cy1, c=prob, s=30)
    plt.scatter(cx2, cy2, c=prob2, s=30)
    plt.gca().set_xlim((2.5, 3.5))
    plt.gca().set_ylim((11, 13))
    plt.xlabel('a')
    plt.ylabel('d')

if plot_type > 0:
    plt.show()
...