如何適應與SciPy的

Question

你會如何使用scipy.optimize.curve_fit在任何y以適應任意鍾/孔曲線或X抵消任何高斯曲線？

我發現了幾個例子，但他們通常是隻能夠適應非常特殊類型的鐘形曲線的，如果最初的猜測是不緊密，或曲線已經反轉或者甚至它們通常會失敗略微偏移。

該段代碼描述我目前的嘗試：

import os,sys 
import numpy as np 
import pylab 
import scipy 
from scipy.optimize import curve_fit 
from pylab import * 
from numpy import * 

def gauss1(x, *p): 
    A, mu, sigma = p 
    return A*np.exp(-(x-mu)**2/(2.*sigma**2)) 

def gauss2(x, *args): 
    m1, m2, m3, s1, s2, s3, k1, k2, k3 = args 
    ret = k1*scipy.stats.norm.pdf(x, loc=m1 ,scale=s1) 
    ret += k2*scipy.stats.norm.pdf(x, loc=m2 ,scale=s2) 
    ret += k3*scipy.stats.norm.pdf(x, loc=m3 ,scale=s3) 
    return ret 

def gauss3(x, a, b, c, d, x0): 
    return a*np.exp(-(x-x0)**2/(2*d**2)) + c 

ydata = np.array([0.001, 0.1, 0.5, 0.75, 0.9, 1.0, 0.9, 0.75, 0.5, 0.1, 0.001]) 
xdata = np.array(xrange(len(ydata))) 

# Test transformations. 
#ydata = 1 - ydata # flip => completely breaks fit 
#ydata = ydata + 100 # offset y => distorts fit 
#xdata = xdata + 10 # offset x => completely breaks fit 

# Equations and initial parameters. 
func, p0 = gauss1, [1., 0., 1.] #doesn't handle invert, of x/y offsets 
#func, p0 = gauss2, [-50, 3, 20, 1, 1, 1, 1, 1, 1] #handles x-offset, but nothing else 
#func, p0 = gauss3, [1,1,1,1,1] #doesn't handle anything 

coeff, var_matrix = curve_fit(func, xdata, ydata, p0=p0) 

yfit = func(xdata, *coeff) 

pylab.plot(xdata, ydata, 'o', label='data', markersize=5) 
plt.plot(xdata, yfit, label='Fitted data') 
pylab.ylim(ydata.min()-1, ydata.max()+1) 
pylab.legend(loc='best') 
fn = 'test5_gaussian.png' 
pylab.savefig(fn) 
pylab.show() 
print 'plot saved to %s' % fn 

我使用3種不同類型的高斯公式，沒有太大的成功。您可以通過取消註釋不同的轉換和/或方程來查看圖表，瞭解每個圖表是如何破解的。

是否有更好的高斯公式我可以使用，或以某種方式計算初始猜測最大化，它會適應可能性有多大？

Answer 1

作爲替代曲線擬合，可以使用method of moments這往往是當更多的穩定應用。爲此，您可以選擇一些統計數據，並根據分佈統計信息對統計數據進行抽樣。

高斯曲線完全由只有兩個參數確定：

1. 位置或平均值
2. 規模或標準差

爲了獲得樣本均值&標準差，則需要通過它們的和，這導致在每個採樣點的概率正常化 -data y：

from scipy.stats import norm 

ys = np.array([0.001, 0.1, 0.5, 0.75, 0.9, 1.0, 0.9, 0.75, 0.5, 0.1, 0.001]) 
xs = np.fromiter(range(len(ys)), dtype=ys.dtype) 

prob = ys/ys.sum()      # probabilities 
mu = prob.dot(xs)      # mean or location param. 
sigma = np.sqrt(prob.dot(xs**2) - mu**2) # standard deviation or scale param. 

xs1 = np.linspace(-2, 12, 100) 
ys1 = norm.pdf(xs1, mu, sigma) * ys.sum() 

fig, ax = plt.subplots() 
ax.plot(xs, ys, 'o') 
ax.plot(xs1, ys1)