python直方图上加正态分布线_python – 使用pymc将两个正态分布(直方图)与MCMC拟合？...

不是最简洁的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,反之亦然.