【问题标题】:Data Fitting in Python for multiple peaksPython中的数据拟合多个峰
【发布时间】:2019-12-12 06:49:39
【问题描述】:
import numpy as np
from sympy.physics.wigner import wigner_6j
import matplotlib.pyplot as plt
xr=np.arange(0,33)
[Jo, Ju, I, Ao, Au]=[4.5, 4.5, 2.5,674.4,929]
Ao=Ao*0.00003335640954804
Au=Au*0.00003335640954804
xr1=100000000/np.array(xr)
positions=xr1
centroid = positions.mean(axis=0)
newo=0.005+100000000/4715.2274
Fo=[]
Fu=[]
new=[]
In=[]
Fomax=Jo+I
Fomin=abs(Jo-I)
Fumax=Ju+I
Fumin=abs(Ju-I)
no=int(2*min(Jo,I)+1)
for i in range(0,no):
    Fo.append(Fomax-i)
nu=int(2*min(Ju,I)+1)
for i in range(0,nu):
    Fu.append(Fumax-i)
for i in range(0,no):
    for j in range(0,nu):
        if abs(Fo[i]-Fu[j])<2:
            new.append(newo+(Ao/2)*(Fo[i]*(Fo[i]+1)-Jo*(Jo+1)-I*(I+1))-(Au/2)*(Fu[j]*(Fu[j]+1)-Ju*(Ju+1)-I*(I+1)))
            In.append((2*Fo[i]+1)*(2*Fu[j]+1)*(wigner_6j(Jo,Fo[i],I,Fu[j],Ju,1))**2/(2*I+1))
            max1=np.max(In)
for i in range(0,len(new)):
    for j in range(0,len(new)):
        if new[i]>new[j]:
            temp=new[j]
            new[j]=new[i]
            new[i]=temp
            temp=In[j]

            In[j]=In[i]
            In[i]=temp
zr=[]
sigma=0.031
x2=[]
y2=[]
y2r=[]
for i in range(0, len(new)):
    mew=new[i]
    for j in range(-100,100):
        c=mew+j/1000
        cc=In[i]*(1/(sigma*(44/7)**0.5))*np.exp(-1*((c-mew)/sigma)**2)
        y2.append(cc)
        x2.append(c)
max2=np.max(y2)
for i in range(0,len(new)):
    In[i]=In[i]/max1
for i in range(0,len(y2)):
    y2[i]=y2[i]/max2
for i in range(0,len(y2)):
    y2r.append(y2[i])
for i   in range(0,15):
    a=5
print(centroid)
fig, ax = plt.subplots()
plt.plot(x2, y2r,label="fitted data")
plt.legend()
plt.show()

我有这段代码,它产生了多个峰值。图 1 显示了多个峰相互重叠的数据,但我试图通过使用这些重叠的峰来仅实现一条曲线,如图 2 中的“红色”线所示。

但问题是必须适合最后一张图片所示的线条

【问题讨论】:

  • 一个问题是“使用什么类型的方程或模型?”我认为,如果您在遍历 X 轴时提取 Y 的最大值,您将需要拟合数据,但适合什么?我从您的绘图图像中手动提取了这些数据点,即使是高阶多项式也不适合。我唯一的建议是使用可以拟合此类数据的样条,但下一个问题是“样条是否足以满足您的要求?”样条曲线是一系列与数据部分相匹配的重叠多项式 - 实际上,是一个分段模型,分段之间具有平滑连接。
  • 您好,请问:哪里合适? ...或者说:要拟合什么数据?
  • @JamesPhillips Gaussian 模型是我想要的。但在这个阶段,任何适合都是非常可观的。我不知道样条,但如果你能做样条发送给我,我会检查它,看看它是否对我有好处..
  • @mikuszefski 我们必须拟合数据.. x 轴的数据在 "x2" 中给出,y 轴是 "x2r" 。蓝线是原始数据,在图 2 中需要红线。我用手做了红线。 ;)
  • ...但是蓝线是构建的,对吧?我猜红线应该是所有蓝线的总和。那么要装什么呢?或者那应该是一些通用数据?抱歉,我还是没听懂。

标签: python scipy curve-fitting least-squares data-fitting


【解决方案1】:

我希望我现在得到这个问题。如果我看对了,问题基本上是x 值不同。最重要的是,所有内容都合并在一个列表中。为了解决这个问题,我将print(centroid) 之后的所有内容都更改为

from scipy.interpolate import interp1d

def partition( inList, n ):
    return zip( *[ iter( inList ) ] * n )

xSplit = partition( x2, 200 ) ###manually set to 200 as data is created with range(-100,100)
ySplit = partition( y2r, 200 )
allx = sorted( x2 )
ally = np.zeros( len( allx ), np.float )
funcDict = dict()
for i in range( len( xSplit ) ):
    funcDict[i] = interp1d( xSplit[i], ySplit[i], kind='linear', bounds_error=False, fill_value=0 ) 
for i in range( len( xSplit ) ):
    ally += funcDict[i]( allx )

fig, ax = plt.subplots()
ax.plot( allx, ally, linewidth=2 )
for col1, col2 in zip( xSplit, ySplit ):
    plt.plot( col1, col2, linestyle='--' )
plt.legend()
plt.show()

给你

这是总和,但使用数据的插值。是这个主意吗?

编辑 似乎 OP 需要更多的信封而不是总和。关于解决方案由 James Phillips 给出。甚至可以使用 numpy 更改来缩短此时间

ally += funcDict[i]( allx )

到

ally = np.maximum(ally, funcDict[i]( allx ) )

然后给出

【讨论】:

  • 它看起来不错..但是当我运行这个程序时它给了我错误.... for i in range(len(xSplit)): TypeError: object of type 'zip' has no len ()
  • 我尝试了所有方法,但错误仍然存​​在。 ''for i in range(len(xSplit)): TypeError: object of type 'zip' has no len()'' 你如何运行代码????
  • @Zewo 你的红线情节让我认为红线适合最大值,而不是总和值。哪个更适合您的工作?
  • 要摆脱类型错误,请将 zip 对象转换为列表。我通过使用 Python 的分号操作符将第 76 行从 "funcDict = dict()" 更改为 "funcDict = dict();xSplit = list(xSplit);ySplit = list(ySplit)" 来链接命令,这样我就可以把在这里的评论中。
  • @JamesPhillips 我只是想获得红线。总和值有点偏差,但如果你能适应最大线,那就太棒了..
【解决方案2】:

将用于求和的@mikuszefski 代码与此代码相结合以获得最大值可能是您正在寻找的,使用它将“print(centroid)”之后的所有内容替换为:

from scipy.interpolate import interp1d

def partition( inList, n ):
    return zip( *[ iter( inList ) ] * n )

xSplit = partition( x2, 200 ) ###manually set to 200 as data is created with range(-100,100)
ySplit = partition( y2r, 200 )
allx = sorted( x2 )
ally = np.zeros( len( allx ), np.float )
funcDict = dict();xSplit = list(xSplit);ySplit = list(ySplit)
for i in range( len( xSplit ) ):
    funcDict[i] = interp1d( xSplit[i], ySplit[i], kind='linear', bounds_error=False, fill_value=0 ) 
for i in range( len( xSplit ) ):
    for j in range(len(allx)):
        ally[j] = max(ally[j],funcDict[i](allx[j])) 

fig, ax = plt.subplots()
ax.plot( allx, ally, linewidth=2 )
for col1, col2 in zip( xSplit, ySplit ):
    plt.plot( col1, col2, linestyle='--' )
plt.legend()
plt.show()

【讨论】:

  • 关于 OP 的最大值和zip 的帮助很好。猜猜这是 Python 2.7 Python 3 的东西。关于最大值,我使用 numpy 的maximum 进行了更新。 ...干杯。
猜你喜欢
  • 2013-10-26
  • 2020-07-15
  • 2021-06-22
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 2016-01-14
  • 1970-01-01
  • 2021-08-09
相关资源
最近更新 更多