使用pymc进行MCMC拟合两个正态分布（直方图）？

Question

使用pymc进行MCMC拟合两个正态分布（直方图）？

14

我正在尝试拟合通过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()

Plot of true values and measured values

我的问题：PyMC（或其变体）能否为上述使用的两个分量提供平均值、振幅和sigma？

请注意，我实际要拟合的函数并不是高斯函数——因此，请使用通用函数（例如我的示例中的GaussFunc），而不是“内置”的pymc.Normal()类型函数来提供示例。

此外，我了解模型选择是另一个问题：因此，在当前噪声下，可能只有1个分量（轮廓）是统计上证明的。但我想看看1、2、3等分量的最佳解决方案。

我也不固执于使用PyMC——如果scikit-learn、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。

- JBWhitmore

虽然我能理解你为什么选择 mc.Uniform('sig', 0.01, 6.5)，但是从安全角度考虑，最好将范围从0开始而不是0.01。数学上，0确实会破坏模型，但在MCMC算法中，0也永远不会被选中。包括0（到0.01）的原因是为了包含真实值在0到0.01之间的小概率。 - Cam.Davidson.Pilon

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Cam.Davidson.Pilon · Accepted Answer

这不是最简明的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的所有特征，反之亦然。