Как включить ошибки для моих данных в минимизацию методом наименьших квадратов lmfit и что это за ошибка для функции conf_interval2d в lmfit?

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

    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: ДАЛЕЕ, когда я пытаюсь использовать следующее, чтобы построить доверительные интервалы, xs, ys, grid = conf_interval2d (mi, 'a', 'b', 20, 20)

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

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

or

Параметр 4 подпрограммы DGESV был неправильным. Ошибка параметра BLAS в Mac OS в DGESV, параметр № 0 (недоступен) равен 0


python uncertainty confidence-interval least-squares
person Nicole Goebel    schedule 06.03.2013    source источник

Ответы (1)


arrow_upward
5
arrow_downward

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

Из lmfit документов (хотя найти их не так-то просто):

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

Однако сделать это не так уж и сложно. Из статьи Википедии о взвешенном подборе методом наименьших квадратов:

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

Однако 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

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

person Thriveth    schedule 12.04.2013
comment
Что здесь означает ошибка? это отклонение? стандартное отклонение ?? - person Mehdi; 16.10.2015
comment
Извините, ошибка в данном случае означает стандартное отклонение. - person Thriveth; 16.10.2015
comment
Я не думаю, что это правильно, поскольку в документации прямо говорится, что остаточная функция должна возвращать только остатки. - person rhody; 08.08.2016
comment
@rhody Остатки возвращаются только в целях минимизации. Вы можете вернуть любой массив для минимизации - и описанная выше процедура дает статистически правильный вектор минимизации для получения взвешенного соответствия. - person Thriveth; 08.08.2016
comment
Вы правы, я нашел еще одно предложение, похороненное в документации, в котором говорится, что мы можем взвесить остатки, что я сделал, и это работает. Однако я не уверен, следует ли вычислять полный хи-квадрат в функции остатков. Документы здесь неясны, но я предполагаю, что после вызова остатков lmfit сам вычисляет хи-квадрат? - person rhody; 09.08.2016
comment
Это хороший момент, и именно поэтому я позволил коду в примере возвращать взвешенный хи-неквадрат. - person Thriveth; 10.08.2016