2013-03-06 7 views
3

Я новичок в python и пытаюсь использовать пакет lmfit для проверки моих собственных вычислений, однако я не уверен (1) относительно того, как включать ошибки для данных (sig) для следующего теста (и 2) Я получаю ошибку с conf_interval2d показано ниже):Как включить ошибки для моих данных в минимизацию минимальных квадратов lmfit, и какова эта ошибка для функции conf_interval2d в lmfit?

import numpy as np 
    from lmfit import Parameters, Minimizer, conf_interval, conf_interval2d, minimize, printfuncs 


    x=np.array([ 0.18, 0.26, 1.14, 0.63, 0.3 , 0.22, 1.16, 0.62, 0.84,0.44, 1.24, 0.89, 1.2 , 0.62, 0.86, 0.45, 1.17, 0.59, 0.85, 0.44]) 
    data=np.array([ 68.59, 71.83, 22.52,44.587,67.474 , 55.765, 20.9,41.33783784,45.79 , 47.88, 6.935, 34.15957447,44.175, 45.89230769, 57.29230769, 60.8,24.24335594, 34.09121287, 42.21504003, 26.61161674]) 
    sig=np.array([ 11.70309409, 11.70309409, 11.70309409, 11.70309409,11.70309409, 11.70309409, 11.70309409, 11.70309409,11.70309409, 11.70309409, 11.70309409, 11.70309409,11.70309409, 11.70309409, 11.70309409, 11.70309409,11.70309409, 11.70309409, 11.70309409, 11.70309409]) 

    def residual(pars, x, data=None): 
     a=pars['a'].value 
     b=pars['b'].value 
     model = a + (b*x) 
     if data is None: 
      return model 
     return model-data 

    params=Parameters() 
    params.add('a', value=70.0) 
    params.add('b', value=40.0) 

    mi=minimize(residual, params, args=(x, data)) 
    #mi=minimize(residual, params, args=(x,), kws={'data': data})#is this more correct? 
    ci, trace = conf_interval(mi, trace=True) 

Это прекрасно работает до сих пор, но, как вопрос выше, как я включаю ошибки для данных (sig_chla), так что я могу рассчитал взвешенные и снижение хи-квадрат ?

Часть 2: КРОМЕ, когда я пытаюсь использовать следующее, так что я могу построить доверительные интервалы, хз, Ю.С., сетки = conf_interval2d (мили, «а», «б», 20, 20)

я получаю следующее сообщение об ошибке:

* ValueError: не удалось создать намерение (кэш | скрыть) | необязательно array-- должен быть определен размеры, но есть (0,)

или

паритет ameter 4 рутинной DGESV неверен ошибка Mac OS BLAS параметр в DGESV, параметр # 0, (нет), 0

ответ

4

Вы должны взвешивать свои данные, ошибки самостоятельно в функции residual().

Из lmfit документации (не очень легко найти, хотя):

Note that the calculation of chi-square and reduced chi-square assume that the returned residual function is scaled properly to the uncertainties in the data. For these statistics to be meaningful, the person writing the function to be minimized must scale them properly.

Это не так трудно сделать, хотя. Из записи Википедии на weighted least-square fitting:

If, however, the measurements are uncorrelated but have different uncertainties, a modified approach might be adopted. Aitken showed that when a weighted sum of squared residuals is minimized, is BLUE if each weight is equal to the reciprocal of the variance of the measurement.

Однако lmfit занимает в остаточном, а не квадрат остатка, так что вместо того, чтобы просто идти

# This is what you do with no errorbars -> equal weights. 
    resids = model - data 
    return resids 

вы должны сделать что-то вроде этого (с scipy ввозимым в качестве sp):

# Do this to include errors as weights. 
    resids = model - data 
    weighted = sp.sqrt(resids ** 2/sig ** 2) 
    return weighted 

Это должно дать вам правильно взвешенное подгонку.

+0

Что здесь означает ошибка? это дисперсия? стандартное отклонение? – Mehdi

+0

Извините, «ошибка» в этом случае означает стандартное отклонение. – Thriveth

+0

Я не думаю, что это правильно, поскольку документы явно говорят, что остаточная функция должна возвращать только остатки. – rhody