Установка гауссова, получение прямой линии. Python 2.7

Question

Как мой заголовок предлагает, я пытаюсь поместить гауссов к некоторым данным, и я просто получаю прямую линию. Я смотрел на эти другие обсуждения Gaussian fit for Python и Fitting a gaussian to a curve in Python, которые, кажется, предлагают в основном то же самое. Я могу заставить код в этих обсуждениях работать нормально для данных, которые они предоставляют, но он не будет делать это для моих данных.

Мой код выглядит следующим образом:

import pylab as plb 
import matplotlib.pyplot as plt 
from scipy.optimize import curve_fit 
from scipy import asarray as ar,exp 

y = y - y[0] # to make it go to zero on both sides 
x = range(len(y)) 

max_y = max(y) 

n = len(y) 
mean = sum(x*y)/n 
sigma = np.sqrt(sum(y*(x-mean)**2)/n) 
# Someone on a previous post seemed to think this needed to have the sqrt. 
# Tried it without as well, made no difference. 

def gaus(x,a,x0,sigma): 
    return a*exp(-(x-x0)**2/(2*sigma**2)) 

popt,pcov = curve_fit(gaus,x,y,p0=[max_y,mean,sigma])  
# It was suggested in one of the other posts I looked at to make the 
# first element of p0 be the maximum value of y. 
# I also tried it as 1, but that did not work either 

plt.plot(x,y,'b:',label='data') 
plt.plot(x,gaus(x,*popt),'r:',label='fit') 
plt.legend() 
plt.title('Fig. 3 - Fit for Time Constant') 
plt.xlabel('Time (s)') 
plt.ylabel('Voltage (V)') 
plt.show() 

Данные, которые я пытаюсь соответствовать следующим образом:

y = array([ 6.95301373e+12, 9.62971320e+12, 1.32501876e+13, 
    1.81150568e+13, 2.46111132e+13, 3.32321345e+13, 
    4.45978682e+13, 5.94819771e+13, 7.88394616e+13, 
    1.03837779e+14, 1.35888594e+14, 1.76677210e+14, 
    2.28196006e+14, 2.92781632e+14, 3.73133045e+14, 
    4.72340762e+14, 5.93892782e+14, 7.41632194e+14, 
    9.19750269e+14, 1.13278296e+15, 1.38551838e+15, 
    1.68291212e+15, 2.02996957e+15, 2.43161742e+15, 
    2.89259207e+15, 3.41725793e+15, 4.00937676e+15, 
    4.67187762e+15, 5.40667931e+15, 6.21440313e+15, 
    7.09421973e+15, 8.04366842e+15, 9.05855930e+15, 
    1.01328502e+16, 1.12585509e+16, 1.24257598e+16, 
    1.36226443e+16, 1.48356404e+16, 1.60496345e+16, 
    1.72482199e+16, 1.84140400e+16, 1.95291969e+16, 
    2.05757166e+16, 2.15360187e+16, 2.23933053e+16, 
    2.31320228e+16, 2.37385276e+16, 2.42009864e+16, 
    2.45114362e+16, 2.46427484e+16, 2.45114362e+16, 
    2.42009864e+16, 2.37385276e+16, 2.31320228e+16, 
    2.23933053e+16, 2.15360187e+16, 2.05757166e+16, 
    1.95291969e+16, 1.84140400e+16, 1.72482199e+16, 
    1.60496345e+16, 1.48356404e+16, 1.36226443e+16, 
    1.24257598e+16, 1.12585509e+16, 1.01328502e+16, 
    9.05855930e+15, 8.04366842e+15, 7.09421973e+15, 
    6.21440313e+15, 5.40667931e+15, 4.67187762e+15, 
    4.00937676e+15, 3.41725793e+15, 2.89259207e+15, 
    2.43161742e+15, 2.02996957e+15, 1.68291212e+15, 
    1.38551838e+15, 1.13278296e+15, 9.19750269e+14, 
    7.41632194e+14, 5.93892782e+14, 4.72340762e+14, 
    3.73133045e+14, 2.92781632e+14, 2.28196006e+14, 
    1.76677210e+14, 1.35888594e+14, 1.03837779e+14, 
    7.88394616e+13, 5.94819771e+13, 4.45978682e+13, 
    3.32321345e+13, 2.46111132e+13, 1.81150568e+13, 
    1.32501876e+13, 9.62971320e+12, 6.95301373e+12, 
    4.98705540e+12]) 

я покажу вам, как это выглядит, но, видимо, я не имеют достаточные очки репутации ...

Кто-нибудь понял, почему он не подходит должным образом?

Спасибо за вашу помощь :)

Answer 1

Важность первоначального предположения, p0 в списке аргументов curve_fit «s по умолчанию, не может быть подчеркнуто достаточно.

Обратите внимание, что docstring упоминает, что

[p0] Если None, то начальные значения не будут все 1

Так что, если вы не предоставите, он будет использовать начальное предположение из 1 для всех параметров, которые вы пытаетесь оптимизировать. Выбор p0 влияет на скорость, с которой базовый алгоритм изменяет вектор угадывания p0 (ссылка на документацию least_squares).

Когда вы посмотрите на данные, которые у вас есть, вы заметите, что максимальное и среднее значение, mu_0, гауссова типа набора данных y, являются 2.4e16 и 49 соответственно. При таком большом значении пика алгоритм должен был бы внести коренные изменения в свое первоначальное предположение, чтобы достичь этого большого значения.

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

Используя данные, вы можете поставить хорошее начальное приближение для peak_value, в mean и sigma, записывая их так:

y = np.array([...]) # starting from the original dataset 
x = np.arange(len(y)) 
peak_value = y.max() 
mean = x[y.argmax()] # observation of the data shows that the peak is close to the center of the interval of the x-data 
sigma = mean - np.where(y > peak_value * np.exp(-.5))[0][0] # when x is sigma in the gaussian model, the function evaluates to a*exp(-.5) 
popt,pcov = curve_fit(gaus, x, y, p0=[peak_value, mean, sigma]) 
print(popt) # prints: [ 2.44402560e+16 4.90000000e+01 1.20588976e+01] 

Обратите внимание, что в вашем коде, для среднего взять sum(x*y)/n, что странно, потому что это будет модулировать гауссовский полином степени 1 (он умножает гауссовский с монотонно увеличивающейся линией постоянного наклона), прежде чем принимать среднее значение. Это компенсирует среднее значение y (в данном случае направо). Аналогичное замечание можно сделать для вашего расчета sigma.

Заключительное замечание: гистограмма y не будет похожа на гауссовскую, так как y уже является гауссовой.Гистограмма будет всего лишь bin (количество) значений в разные категории (отвечая на вопрос «сколько точек данных в y достигает значения от [a, b]?»).