3

scipyのcurve_fit関数を使用して、歪んでシフトしたガウス曲線を近似しようとしていますが、特定の条件下では近似が非常に悪く、直線に近いか、正確に直線になることがよくあります。

以下のコードは、ドキュメントから派生していcurve_fitます。提供されているコードは、テスト目的の任意のデータセットですが、問題を非常によく示しています。

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

#def func(x, a, b, c):
#    return a*np.exp(-b*x) + c

def func(x, sigmag, mu, alpha, c,a):
    #normal distribution
    normpdf = (1/(sigmag*np.sqrt(2*math.pi)))*np.exp(-(np.power((x-mu),2)/(2*np.power(sigmag,2))))
    normcdf = (0.5*(1+sp.erf((alpha*((x-mu)/sigmag))/(np.sqrt(2)))))
    return 2*a*normpdf*normcdf + c

x = np.linspace(0,100,100)
y = func(x, 10,30, 0,0,1)
yn = y + 0.001*np.random.normal(size=len(x))

popt, pcov = curve_fit(func, x, yn,) #p0=(9,35,0,9,1))

y_fit= func(x,popt[0],popt[1],popt[2],popt[3],popt[4])

plt.plot(x,yn)
plt.plot(x,y_fit)

ガウス分布をゼロから大きくシフトしすぎると、この問題が発生するようです（を使用mu）。元の関数と同じものでも初期値を与えてみましたが、問題は解決しません。の値の場合mu=10、curve_fitは完全に機能しますが、使用するmu>=30とデータに適合しなくなります。

4

2 に答える 2

3

最小化の出発点を与えることは、しばしば不思議に働きます。最小値に最大値の位置と曲線の幅に関する情報を提供してみてください。

popt, pcov = curve_fit(func, x, yn, p0=(1./np.std(yn), np.argmax(yn) ,0,0,1))

コード内のこの1行をで変更するsigma=10とmu=50、 ここに画像の説明を入力してください

于 2013-03-14T13:23:54.380 に答える
3

curve_fitランダムな初期推測で何度も呼び出すことができ、エラーが最小のパラメーターを選択できます。

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

def func(x, sigmag, mu, alpha, c,a):
    #normal distribution
    normpdf = (1/(sigmag*np.sqrt(2*math.pi)))*np.exp(-(np.power((x-mu),2)/(2*np.power(sigmag,2))))
    normcdf = (0.5*(1+sp.erf((alpha*((x-mu)/sigmag))/(np.sqrt(2)))))
    return 2*a*normpdf*normcdf + c

x = np.linspace(0,100,100)
y = func(x, 10,30, 0,0,1)
yn = y + 0.001*np.random.normal(size=len(x))

results = []
for i in xrange(50):
    p = np.random.randn(5)*10
    try:
        popt, pcov = curve_fit(func, x, yn, p)
    except:
        pass
    err = np.sum(np.abs(func(x, *popt) - yn))
    results.append((err, popt))
    if err < 0.1:
        break

err, popt = min(results, key=lambda x:x[0])
y_fit= func(x, *popt)

plt.plot(x,yn)
plt.plot(x,y_fit)
print len(results)
于 2013-03-14T08:55:52.870 に答える