尝试使用任意曲线在Python中拟合多个Moyal函数

问题描述

我一直试图拟合具有多个Landau峰的分布（所以我正在使用scipy.stats moyal）。我曾尝试使用scipy curve_fit，但我一直收到错误消息。该函数将具有许多曲线，由“ nMIP”指定。功能如下：

x = np.linspace(0,len(ADC[0])-1,len(ADC[0]))

def multMoyal(x,nMIPs,*args):
    moy = 0
    for i in range(int(nMIPs)):
        r = 3*i
        moy = moy + args[r]*moyal.pdf(x,args[r+1],args[r+2])
    return moy

其中ADC是光谱。这是粒子碰撞数据；它实际上是带电粒子穿过检测器中的瓷砖，并且每条态曲线都表示带电粒子的数量（通常为2-5，这取决于实验参数和瓷砖位置）。

当我用nMIPs = 2调用它时，该函数可以正常工作：

args = [1,1,0.5,0.8,2,0.5]
y = multMoyal(x,*args)

但是当我尝试在curve_fit中输入它时，会不断出现错误（“ IndexError：元组索引超出范围”）。 curve_fit线如下：

    param,param_cov = curve_fit(multMoyal(x=x,nMIPs=nMIPs),x,y[0],maxfev=5000)

其中y [0]是ADC频谱之一，但在x上进行了平滑处理以消除噪声。

那么，如何获得任意数量的曲线却仍可以在scipy的curve_fit中工作？还是我应该使用curve_fit以外的东西？

相关：只要我可以使用它，如何设置某些参数，但其他参数不确定？（我想对第一条曲线和权重的峰值进行猜测，但将其他所有内容留空。）

谢谢！

解决方法

我找到了一个解决方案，但我不知道这是否是正确的方法（但它确实有效）。我将功能更新为：

def multMoyal(x,nMIPs,*args):
    moy = 0
    for i in range(int(nMIPs)):
        l = 3*i
        moy = moy + args[l]*moyal.pdf(x,args[l+1],args[l+2])
    return moy

然后，我添加了界限和猜测。最初对scipy.optimise.curve_fit的猜测都是1，所以您输入的几乎所有内容都会比这更好。我选择了一些似乎适合该问题的范围，然后从那里解决了。要设置单个参数（或您喜欢的多个参数），我使用了界限。将边界设置为+-np.inf将使该参数自由。

例如，我使用了从nMIP到nMIPs + 1的nMIP的绑定，这实际上将其设置为nMIP（当然，是对nMIP的猜测）。看起来如下：

fitStart = []
startVal = []
weightGuess = []
guess = []
boundsInit = []
boundsEnd = []
for q in range(int(len(ADC))):
    MIPguess = exMax[q][np.where(exMax[q] > exMin[q][0])[0][0]]-10
    starter = int(exMin[q][0]+abs(exMin[q][0]-MIPguess)/2)
    fitStart.append(starter)
    startVal.append(MIPguess)
    weightGuess.append(y[q][MIPguess])
    guess.append([nMIPs,weightGuess[q],startVal[q],0.2*startVal[q]])
    boundsInit.append([nMIPs,fitStart[q],0])
    boundsEnd.append([nMIPs+1,np.inf,0.3*startVal[q]])
    for i in range(1,nMIPs):
        l = i+1
        guess[q].extend([weightGuess[q]/l,l*startVal[q],0.2*startVal[q]])
        boundsInit[q].extend([0,fitStart[q]*l*0.75,0])
        boundsEnd[q].extend([np.inf,np.inf])

如果我只想设置nMIP，则可以将所有开始的猜测都设为1，并将界限（除了nMIP之外）设置为+-np.inf。然后按照以下步骤进行拟合：

param,param_cov = curve_fit(
            multMoyal,x[fitStart[q]:],ADC[q][fitStart[q]:],maxfev=5000,p0=guess[q],bounds=(boundsInit[q],boundsEnd[q]))

它可以工作，但是我非常怀疑这是最严格的代码。我仍然觉得我缺少一种更简单的方法，但是现在就可以了。我很想知道专家们对此可能会回答什么，但是我将其作为答案，因为还没有任何问题，至少可以解决问题。

curve-fitting python scipy scipy