问题:我想将经验数据拟合到双峰正态分布中,从物理背景中我知道峰值的距离(固定),而且两个峰值必须具有相同的标准差。在
我试图用^{
细节:我避免了loc
和scale
参数,并将它们作为m
和{_pdf
-方法中,因为峰值距离{fit
-方法中将它们固定为floc=0
和fscale=1
,实际上需要m
,s
的拟合参数和w
的权重
我在示例数据中期望的是一个峰值在x=-450
和x=450
(=>;m=0
)附近的分布。stdev s
应该在100或200左右,但不是1.0,权重w
应该大约是0.5
from __future__ import division
from scipy.stats import rv_continuous
import numpy as np
class norm2_gen(rv_continuous):
def _argcheck(self, *args):
return True
def _pdf(self, x, m, s, w, delta):
return np.exp(-(x-m+delta/2)**2 / (2. * s**2)) / np.sqrt(2. * np.pi * s**2) * w + \
np.exp(-(x-m-delta/2)**2 / (2. * s**2)) / np.sqrt(2. * np.pi * s**2) * (1 - w)
norm2 = norm2_gen(name='norm2')
data = [487.0, -325.5, -159.0, 326.5, 538.0, 552.0, 563.0, -156.0, 545.5, 341.0, 530.0, -156.0, 473.0, 328.0, -319.5, -287.0, -294.5, 153.5, -512.0, 386.0, -129.0, -432.5, -382.0, -346.5, 349.0, 391.0, 299.0, 364.0, -283.0, 562.5, -42.0, 214.0, -389.0, 42.5, 259.5, -302.5, 330.5, -338.0, 508.5, 319.5, -356.5, 421.5, 543.0]
m, s, w, delta, loc, scale = norm2.fit(data, fdelta=900, floc=0, fscale=1)
print m, s, w, delta, loc, scale
>>> 1.0 1.0 1.0 900 0 1
我在做了一些调整后,使您的分布符合数据:
w
时,有一个隐含的约束:0<;=w
<;=1。fit()
方法使用的解算器没有意识到这个约束,因此w
可能以不合理的值结束。处理此类约束的一种方法是允许w
为任意实数,但在PDF的公式中,使用phi = 0.5 + arctan(w)/pi
将w
转换为0和1之间的一个分数{fit()
方法使用数值优化例程来寻找最大似然估计。像大多数这样的例程一样,它需要一个优化的起点。默认的起点都是1,但这并不总是有效的。您可以通过在数据后面提供值作为fit()
的位置参数来选择不同的起点。我在脚本中使用的值是有效的;我没有探究结果对这些起始值的敏感程度。在我做了两个估计。在第一个中,我让
delta
作为一个自由参数,在第二个中,我将delta
固定为900。在下面的脚本生成以下绘图:
脚本如下:
相关问题 更多 >
编程相关推荐