19

私の数学の知識は限られているので、おそらく行き詰まっています。2つのガウスピークを適合させようとしているスペクトルがあります。最大のピークにはフィットできますが、最小のピークにはフィットできません。2つのピークのガウス関数を合計する必要があることは理解していますが、どこが間違っているのかわかりません。現在の出力の画像が表示されます：

電流出力

青い線は私のデータであり、緑の線は私の現在の適合です。次のコードを使用して、現在フィットしようとしているデータのメインピークの左側に肩があります。

import matplotlib.pyplot as pt
import numpy as np
from scipy.optimize import leastsq
from pylab import *

time = []
counts = []


for i in open('/some/folder/to/file.txt', 'r'):
    segs = i.split()
    time.append(float(segs[0]))
    counts.append(segs[1])

time_array = arange(len(time), dtype=float)
counts_array = arange(len(counts))
time_array[0:] = time
counts_array[0:] = counts


def model(time_array0, coeffs0):
    a = coeffs0[0] + coeffs0[1] * np.exp( - ((time_array0-coeffs0[2])/coeffs0[3])**2 )
    b = coeffs0[4] + coeffs0[5] * np.exp( - ((time_array0-coeffs0[6])/coeffs0[7])**2 ) 
    c = a+b
    return c


def residuals(coeffs, counts_array, time_array):
    return counts_array - model(time_array, coeffs)

# 0 = baseline, 1 = amplitude, 2 = centre, 3 = width
peak1 = np.array([0,6337,16.2,4.47,0,2300,13.5,2], dtype=float)
#peak2 = np.array([0,2300,13.5,2], dtype=float)

x, flag = leastsq(residuals, peak1, args=(counts_array, time_array))
#z, flag = leastsq(residuals, peak2, args=(counts_array, time_array))

plt.plot(time_array, counts_array)
plt.plot(time_array, model(time_array, x), color = 'g') 
#plt.plot(time_array, model(time_array, z), color = 'r')
plt.show()
4

3 に答える 3

18

このコードは、2つのガウス分布の組み合わせである関数のみを適合させるという条件で私のために機能しました。

2つのガウス関数を加算し、実際のデータからそれらを減算する残差関数を作成しました。

Numpyの最小二乗関数に渡したパラメーター（p）には、最初のガウス関数の平均（m）、1番目と2番目のガウス関数との平均の差（dm、つまり水平シフト）、標準偏差が含まれます。最初の（sd1）の、および2番目の（sd2）の標準偏差。

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

######################################
# Setting up test data
def norm(x, mean, sd):
  norm = []
  for i in range(x.size):
    norm += [1.0/(sd*np.sqrt(2*np.pi))*np.exp(-(x[i] - mean)**2/(2*sd**2))]
  return np.array(norm)

mean1, mean2 = 0, -2
std1, std2 = 0.5, 1 

x = np.linspace(-20, 20, 500)
y_real = norm(x, mean1, std1) + norm(x, mean2, std2)

######################################
# Solving
m, dm, sd1, sd2 = [5, 10, 1, 1]
p = [m, dm, sd1, sd2] # Initial guesses for leastsq
y_init = norm(x, m, sd1) + norm(x, m + dm, sd2) # For final comparison plot

def res(p, y, x):
  m, dm, sd1, sd2 = p
  m1 = m
  m2 = m1 + dm
  y_fit = norm(x, m1, sd1) + norm(x, m2, sd2)
  err = y - y_fit
  return err

plsq = leastsq(res, p, args = (y_real, x))

y_est = norm(x, plsq[0][0], plsq[0][2]) + norm(x, plsq[0][0] + plsq[0][1], plsq[0][3])

plt.plot(x, y_real, label='Real Data')
plt.plot(x, y_init, 'r.', label='Starting Guess')
plt.plot(x, y_est, 'g.', label='Fitted')
plt.legend()
plt.show()

コードの結果。

于 2012-04-13T23:36:22.480 に答える
17

scikit-learnのガウス混合モデルを使用できます。

from sklearn import mixture
import matplotlib.pyplot
import matplotlib.mlab
import numpy as np
clf = mixture.GMM(n_components=2, covariance_type='full')
clf.fit(yourdata)
m1, m2 = clf.means_
w1, w2 = clf.weights_
c1, c2 = clf.covars_
histdist = matplotlib.pyplot.hist(yourdata, 100, normed=True)
plotgauss1 = lambda x: plot(x,w1*matplotlib.mlab.normpdf(x,m1,np.sqrt(c1))[0], linewidth=3)
plotgauss2 = lambda x: plot(x,w2*matplotlib.mlab.normpdf(x,m2,np.sqrt(c2))[0], linewidth=3)
plotgauss1(histdist[1])
plotgauss2(histdist[1])

ここに画像の説明を入力してください

以下の関数を使用して、必要なガウス数をncompパラメーターに適合させることもできます。

from sklearn import mixture
%pylab

def fit_mixture(data, ncomp=2, doplot=False):
    clf = mixture.GMM(n_components=ncomp, covariance_type='full')
    clf.fit(data)
    ml = clf.means_
    wl = clf.weights_
    cl = clf.covars_
    ms = [m[0] for m in ml]
    cs = [numpy.sqrt(c[0][0]) for c in cl]
    ws = [w for w in wl]
    if doplot == True:
        histo = hist(data, 200, normed=True)
        for w, m, c in zip(ws, ms, cs):
            plot(histo[1],w*matplotlib.mlab.normpdf(histo[1],m,np.sqrt(c)), linewidth=3)
    return ms, cs, ws
于 2013-10-04T13:48:09.153 に答える
4

係数0と4は縮退しています。データには、これらの間で決定できるものはまったくありません。2つではなく1つのゼロレベルパラメータを使用する必要があります（つまり、コードから1つを削除します）。これがおそらくあなたの適合を妨げているものです（これは不可能であるというここでのコメントを無視してください-そのデータには明らかに少なくとも2つのピークがあり、確かにそれに適合できるはずです）。

（なぜ私がこれを提案しているのかは明確ではないかもしれませんが、何が起こっているのかというと、係数0と4は互いに打ち消し合う可能性があります。両方ともゼロにすることも、一方を100にして、もう一方を-100にすることもできます。これは同じように良いです。これはフィッティングルーチンを「混乱」させます。フィッティングルーチンは、正しい答えが1つもない場合に、それらがどうあるべきかを理解するために時間を費やします。フィット感は同じになります）。

実際、プロットからは、ゼロレベルはまったく必要ないように見えます。私はそれらの両方をドロップして、フィット感がどのように見えるかを見てみます。

また、係数1と5（またはゼロ点）を最小二乗法に適合させる必要はありません。代わりに、モデルは線形であるため、ループごとに値を計算できます。これにより処理が速くなりますが、重要ではありません。私はあなたがあなたの数学があまり良くないと言っているのに気づいたので、おそらくこれを無視してください。

于 2012-04-14T01:05:34.507 に答える