我试图与在CCD光谱仪检测,以适应线轮廓。 为了便于审议,我已经包含了一个示范,如果解决了,是非常相似的一个其实我是想解决的问题。
我看这样的: https://stats.stackexchange.com/questions/46626/fitting-model-for-two-normal-distributions-in-pymc和其他各种问题和答案,但他们正在做的事情根本不同比我想做的事情。
import pymc as mc
import numpy as np
import pylab as pl
def GaussFunc(x, amplitude, centroid, sigma):
return amplitude * np.exp(-0.5 * ((x - centroid) / sigma)**2)
wavelength = np.arange(5000, 5050, 0.02)
# Profile 1
centroid_one = 5025.0
sigma_one = 2.2
height_one = 0.8
profile1 = GaussFunc(wavelength, height_one, centroid_one, sigma_one, )
# Profile 2
centroid_two = 5027.0
sigma_two = 1.2
height_two = 0.5
profile2 = GaussFunc(wavelength, height_two, centroid_two, sigma_two, )
# Measured values
noise = np.random.normal(0.0, 0.02, len(wavelength))
combined = profile1 + profile2 + noise
# If you want to plot what this looks like
pl.plot(wavelength, combined, label="Measured")
pl.plot(wavelength, profile1, color='red', linestyle='dashed', label="1")
pl.plot(wavelength, profile2, color='green', linestyle='dashed', label="2")
pl.title("Feature One and Two")
pl.legend()
我的问题是:在 PyMC(或其变型)给我的平均值,幅度和西格玛上面使用的两个组件?
请注意,我会真正适合我的真正的问题的功能并不高斯 - 所以请提供使用(在我的例子一样GaussFunc)的通用功能的例子,而不是一个“内置” pymc.Normal()型功能。
另外,我理解模型的选择是另外一个问题:所以用电流噪声,1个组件(配置文件)可能是所有在统计学上是有道理的。 不过,我想看看什么是最好的解决方案,1,2,3等成分会。
我还没有结婚使用PyMC的想法 - 如果scikit学习,astroML,或其它的什么包似乎完美的,请让我知道!
编辑:
我失败了多种方法,但我认为是正确的轨道上的事情之一是以下几点:
sigma_mc_one = mc.Uniform('sig', 0.01, 6.5)
height_mc_one = mc.Uniform('height', 0.1, 2.5)
centroid_mc_one = mc.Uniform('cen', 5015., 5040.)
但我不能构建工作了mc.model。
这不是最简洁的PyMC代码,但我做了决定,以帮助读者。 这应该运行,并给予(真的)准确的结果。
我做了使用统一的先验,与宽松的范围,因为我真的不知道我们正在建模的决定。 但是可能有一个围绕质心位置的想法,并且可以使用更好的先验存在。
### Suggested one runs the above code first.
### Unknowns we are interested in
est_centroid_one = mc.Uniform("est_centroid_one", 5000, 5050 )
est_centroid_two = mc.Uniform("est_centroid_two", 5000, 5050 )
est_sigma_one = mc.Uniform( "est_sigma_one", 0, 5 )
est_sigma_two = mc.Uniform( "est_sigma_two", 0, 5 )
est_height_one = mc.Uniform( "est_height_one", 0, 5 )
est_height_two = mc.Uniform( "est_height_two", 0, 5 )
#std deviation of the noise, converted to precision by tau = 1/sigma**2
precision= 1./mc.Uniform("std", 0, 1)**2
#Set up the model's relationships.
@mc.deterministic( trace = False)
def est_profile_1(x = wavelength, centroid = est_centroid_one, sigma = est_sigma_one, height= est_height_one):
return GaussFunc( x, height, centroid, sigma )
@mc.deterministic( trace = False)
def est_profile_2(x = wavelength, centroid = est_centroid_two, sigma = est_sigma_two, height= est_height_two):
return GaussFunc( x, height, centroid, sigma )
@mc.deterministic( trace = False )
def mean( profile_1 = est_profile_1, profile_2 = est_profile_2 ):
return profile_1 + profile_2
observations = mc.Normal("obs", mean, precision, value = combined, observed = True)
model = mc.Model([est_centroid_one,
est_centroid_two,
est_height_one,
est_height_two,
est_sigma_one,
est_sigma_two,
precision])
#always a good idea to MAP it prior to MCMC, so as to start with good initial values
map_ = mc.MAP( model )
map_.fit()
mcmc = mc.MCMC( model )
mcmc.sample( 50000,40000 ) #try running for longer if not happy with convergence.
重要
请记住,该算法是不可知的标签,所以结果可能会显示profile1
与所有从特性profile2
,反之亦然。