Фиксация до гистограммы Пуассона

Я пытаюсь подогнать кривую по гистограмме распределения Пуассона, которая выглядит так:

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

def histo(bsize):
    N = bsize
    #binwidth
    bw = (dt.max()-dt.min())/(N-1.)
    bin1 = dt.min()+ bw*np.arange(N)
    #define the array to hold the occurrence count
    bincount= np.array([])
    for bin in bin1:
        count = np.where((dt>=bin)&(dt<bin+bw))[0].size
        bincount = np.append(bincount,count)
    #bin center
    binc = bin1+0.5*bw
    plt.figure()
    plt.plot(binc,bincount,drawstyle= 'steps-mid')
    plt.xlabel("Interval[ticks]")
    plt.ylabel("Frequency")
histo(30)
plt.xlim(0,.5e8)
plt.ylim(0,25000)
import numpy as np
from scipy.optimize import curve_fit
delta_t = 1.42e7
def func(x, t):
    return t * np.exp(- delta_t/t) 
popt, pcov = curve_fit(func, np.arange(0,.5e8),histo(30))
plt.plot(popt)

Ответ 1

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

Привязанные наименьшие квадраты Fit

В общем, вы не можете получить много, гораздо проще:

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


# get poisson deviated random numbers
data = np.random.poisson(2, 1000)

# the bins should be of integer width, because poisson is an integer distribution
entries, bin_edges, patches = plt.hist(data, bins=11, range=[-0.5, 10.5], normed=True)

# calculate binmiddles
bin_middles = 0.5*(bin_edges[1:] + bin_edges[:-1])

# poisson function, parameter lamb is the fit parameter
def poisson(k, lamb):
    return (lamb**k/factorial(k)) * np.exp(-lamb)

# fit with curve_fit
parameters, cov_matrix = curve_fit(poisson, bin_middles, entries) 

# plot poisson-deviation with fitted parameter
x_plot = np.linspace(0, 20, 1000)

plt.plot(x_plot, poisson(x_plot, *parameters), 'r-', lw=2)
plt.show()

Это результат:

Нераспределенный максимальный размер правдоподобия

Еще лучшая возможность - не использовать гистограмму вообще и вместо этого сделать максимальную правдоподобие.

Но при ближайшем рассмотрении даже это не нужно, поскольку оценка максимального правдоподобия для параметра пуассоновского распределения - среднее арифметическое.

Однако, если у вас есть другие, более сложные pdf файлы, вы можете использовать это в качестве примера:

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
from scipy.misc import factorial


def poisson(k, lamb):
    """poisson pdf, parameter lamb is the fit parameter"""
    return (lamb**k/factorial(k)) * np.exp(-lamb)


def negLogLikelihood(params, data):
    """ the negative log-Likelohood-Function"""
    lnl = - np.sum(np.log(poisson(data, params[0])))
    return lnl


# get poisson deviated random numbers
data = np.random.poisson(2, 1000)

# minimize the negative log-Likelihood

result = minimize(negLogLikelihood,  # function to minimize
                  x0=np.ones(1),     # start value
                  args=(data,),      # additional arguments for function
                  method='Powell',   # minimization method, see docs
                  )
# result is a scipy optimize result object, the fit parameters 
# are stored in result.x
print(result)

# plot poisson-deviation with fitted parameter
x_plot = np.linspace(0, 20, 1000)

plt.hist(data, bins=np.arange(15) - 0.5, normed=True)
plt.plot(x_plot, poisson(x_plot, result.x), 'r-', lw=2)
plt.show()

Ответ 2

Спасибо за прекрасную дискуссию!

Вы можете рассмотреть следующие вопросы:

1) Вместо вычисления "Пуассона", вычислите "лог-Пуассона", для лучшего численного поведения

2) Вместо использования "lamb", используйте логарифм (позвольте мне назвать его "log_mu"), чтобы избежать подгонки "блуждающих" в отрицательные значения "mu". Так

log_poisson(k, log_mu): return k*log_mu - loggamma(k+1) - math.exp(log_mu)

Где "loggamma" - это функция scipy.special.loggamma.

На самом деле, в приведенном выше примере термин "loggamma" только добавляет постоянное смещение к минимизируемым функциям, поэтому можно просто сделать:

log_poisson_(k, log_mu): return k*log_mu - math.exp(log_mu)

ПРИМЕЧАНИЕ: log_poisson_() не то же самое, что log_poisson(), но при использовании для минимизации описанным выше способом даст тот же подобранный минимум (то же значение mu, вплоть до числовых проблем). Значение минимизируемой функции будет смещено, но обычно это все равно не волнует.