доверительный интервал с наименьшим соответствием в питона - PullRequest
Новая
       15

доверительный интервал с наименьшим соответствием в питона

4 голосов
/ 28 апреля 2011

Как рассчитать доверительный интервал для наименьших квадратов (scipy.optimize.leastsq) в питоне?

Ответы [ 3 ]

8 голосов
/ 28 апреля 2011

Я бы использовал метод начальной загрузки.
См. Здесь: http://phe.rockefeller.edu/LogletLab/whitepaper/node17.html

Простой пример для шумного гауссиана: 

x = arange(-10, 10, 0.01)

# model function
def f(p):
    mu, s = p
    return exp(-(x-mu)**2/(2*s**2))

# create error function for dataset    
def fff(d):
    def ff(p):
        return d-f(p)
    return ff

# create noisy dataset from model
def noisy_data(p):
    return f(p)+normal(0,0.1,len(x))

# fit dataset to model with least squares    
def fit(d):
    ff = fff(d)
    p = leastsq(ff,[0,1])[0]
    return p

# bootstrap estimation        
def bootstrap(d):
    p0 = fit(d)
    residuals = f(p0)-d
    s_residuals = std(residuals)

    ps = []
    for i in range(1000):
        new_d = d+normal(0,s_residuals,len(d))
        ps.append(fit(new_d))

    ps = array(ps)
    mean_params = mean(ps,0)
    std_params = std(ps,0)

    return mean_params, std_params

data = noisy_data([0.5, 2.1])
mean_params, std_params = bootstrap(data)

print "95% confidence interval:"
print "mu: ", mean_params[0], " +/- ", std_params[0]*1.95996
print "sigma: ", mean_params[1], " +/- ", std_params[1]*1.95996
4 голосов
/ 28 апреля 2011

Я не уверен, что вы подразумеваете под доверительным интервалом.

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

Как указывалось в строке документации функции, вы могли бы использовать эту информацию вместе с остатками (разницамежду выбранным решением и фактическими данными), чтобы вычислить ковариацию оценок параметров, которая является локальной оценкой доверительного интервала.

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

2 голосов
/ 28 апреля 2011

Самый простой способ оценки доверительного интервала (CI) - умножить стандартные ошибки (стандартное отклонение) на константу. Чтобы вычислить постоянную, вам нужно знать количество степеней свободы (DOF) и уровень достоверности, для которого вы хотите рассчитать CI. Оцениваемый таким образом CI иногда называют асимптотическим CI. Подробнее об этом можно прочитать в статье «Подгонка моделей к биологическим данным с использованием линейной и нелинейной регрессии», автор Motulsky & Christopoulos ( google books ). Эта же книга (или очень похожая) доступна бесплатно как руководство для авторского программного обеспечения .

Вы также можете прочитать , как рассчитать CI с использованием библиотеки C ++ Boost.Math . В этом примере CI рассчитывается для распределения одной переменной. В случае подгонки наименьших квадратов DOF - это не N -1, а N-M , где M - количество параметров. Это должно быть легко сделать то же самое в Python.

Это самая простая оценка. Я не знаю метод начальной загрузки, предложенный Зефиром, но он может быть более надежным, чем метод, о котором я писал.

...