Гауссовская подгонка для Python

Я пытаюсь поместить гауссовские данные (что уже грубо гауссовское). Я уже проконсультировал их здесь и попробовал curve_fit и leastsq, но я думаю, что мне не хватает чего-то более фундаментального (в этом я понятия не имею, как использовать эту команду). Вот посмотрите на script Я до сих пор

import pylab as plb
import matplotlib.pyplot as plt

# Read in data -- first 2 rows are header in this example. 
data = plb.loadtxt('part 2.csv', skiprows=2, delimiter=',')

x = data[:,2]
y = data[:,3]
mean = sum(x*y)
sigma = sum(y*(x - mean)**2)

def gauss_function(x, a, x0, sigma):
    return a*np.exp(-(x-x0)**2/(2*sigma**2))
popt, pcov = curve_fit(gauss_function, x, y, p0 = [1, mean, sigma])
plt.plot(x, gauss_function(x, *popt), label='fit')

# plot data

plt.plot(x, y,'b')

# Add some axis labels

plt.legend()
plt.title('Fig. 3 - Fit for Time Constant')
plt.xlabel('Time (s)')
plt.ylabel('Voltage (V)')
plt.show()

То, что я получаю от этого, - гауссова форма, которая является моими исходными данными, и прямая горизонтальная линия.

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

Ответы

Ответ 1

Здесь скорректированный код:

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

x = ar(range(10))
y = ar([0,1,2,3,4,5,4,3,2,1])

n = len(x)                          #the number of data
mean = sum(x*y)/n                   #note this correction
sigma = sum(y*(x-mean)**2)/n        #note this correction

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

popt,pcov = curve_fit(gaus,x,y,p0=[1,mean,sigma])

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

результат:

Ответ 2

Вы получаете горизонтальную прямую, потому что она не сходится.

Лучшая сходимость достигается, если первый параметр фитинга (p0) помещен как max (y), 5 в примере вместо 1.

Ответ 3

Потеряв часы, пытаясь найти мою ошибку, проблема заключается в вашей формуле:

sigma = sum (y * (x-mean) ** 2)/n это неверно, правильная формула - это квадрат этого!!

SQRT (сумма (у * (х-среднее) ** 2)/п)

Надеюсь, что это поможет

Ответ 4

Существует еще один способ выполнить подгонку, используя пакет 'lmfit'. Он в основном использует cuve_fit, но намного лучше подходит для установки и предлагает комплексную установку. Подробные пошаговые инструкции приведены в приведенной ниже ссылке. http://cars9.uchicago.edu/software/python/lmfit/model.html#model.best_fit

Ответ 5

sigma = sum(y*(x - mean)**2)

должен быть

sigma = np.sqrt(sum(y*(x - mean)**2))

Ответ 6

Объяснение

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

Оценочные параметры должны быть близки к конечным значениям (используйте взвешенное среднее арифметическое - разделите на сумму всех значений):

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

x = np.arange(10)
y = np.array([0, 1, 2, 3, 4, 5, 4, 3, 2, 1])

# correction for weighted arithmetic mean
mean = sum(x * y) / sum(y)
sigma = np.sqrt(sum(y * (x - mean)**2) / sum(y))

def Gauss(x, a, x0, sigma):
    return a * np.exp(-(x - x0)**2 / (2 * sigma**2))

popt,pcov = curve_fit(Gauss, x, y, p0=[max(y), mean, sigma])

plt.plot(x, y, 'b+:', label='data')
plt.plot(x, Gauss(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()

Я предпочитаю использовать numpy.

Комментарий к определению среднего (включая ответ разработчика)

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

Ваше определение возвращает:

>>> sum(x * y)
125

Определение разработчика возвращается:

>>> sum(x * y) / len(x)
12.5 #for Python 3.x

Средневзвешенное среднее арифметическое:

>>> sum(x * y) / sum(y)
5.0

Аналогичным образом вы можете сравнить определения стандартного отклонения (sigma). Сравните с фигурой результирующей подгонки:

Комментарий для пользователей Python 2.x

В Python 2.x вы должны дополнительно использовать новое подразделение, чтобы не запускаться в странные результаты или явно преобразовывать числа перед делением:

from __future__ import division

или, например,

sum(x * y) * 1. / sum(y)