PyMC 中两个正态分布的拟合模型

机器算法验证 造型 Python pymc
2022-03-24 16:17:27

由于我是一名试图了解更多统计数据的软件工程师，所以在我开始之前你必须原谅我，这是一个严肃的新手领域......

我一直在学习PyMC并研究了一些非常（非常）简单的例子。我无法开始工作（也找不到任何相关示例）的一个问题是将模型拟合到从两个正态分布生成的数据。

假设我有 1000 个值；500 从 a 生成Normal(mean=100, stddev=20)，另外 500 从 a 生成Normal(mean=200, stddev=20)。

如果我想为它们拟合模型，即使用 PyMC 确定两个均值和单个标准差。我知道这有点像......

mean1 = Uniform('mean1', lower=0.0, upper=200.0)
mean2 = Uniform('mean2', lower=0.0, upper=200.0)
precision = Gamma('precision', alpha=0.1, beta=0.1)

data = read_data_from_file_or_whatever()

@deterministic(plot=False)
def mean(m1=mean1, m2=mean2):
    # but what goes here?

process = Normal('process', mu=mean, tau=precision, value=data, observed=True)

即，生成过程是正常的，但 mu 是两个值之一。我只是不知道如何表示一个值是否来自m1or之间的“决定” m2。

也许我只是完全采取了错误的方法来建模这个？谁能给我举个例子？我可以阅读 BUGS 和 JAGS，所以一切都很好。

2个回答

您绝对确定一半来自一个发行版，另一半来自另一个发行版吗？如果不是，我们可以将比例建模为随机变量（这是一个非常贝叶斯的事情）。

以下是我会做的，嵌入了一些提示。

from pymc import *

size = 10
p = Uniform( "p", 0 , 1) #this is the fraction that come from mean1 vs mean2

ber = Bernoulli( "ber", p = p, size = size) # produces 1 with proportion p.

precision = Gamma('precision', alpha=0.1, beta=0.1)

mean1 = Normal( "mean1", 0, 0.001 ) #better to use normals versus Uniforms (unless you are certain the value is  truncated at 0 and 200 
mean2 = Normal( "mean2", 0, 0.001 )

@deterministic
def mean( ber = ber, mean1 = mean1, mean2 = mean2):
    return ber*mean1 + (1-ber)*mean2


#generate some artificial data   
v = np.random.randint( 0, 2, size)
data = v*(10+ np.random.randn(size) ) + (1-v)*(-10 + np.random.randn(size ) )


obs = Normal( "obs", mean, precision, value = data, observed = True)

model = Model( {"p":p, "precision": precision, "mean1": mean1, "mean2":mean2, "obs":obs} )

与上述讨论有关的几点：

  1. 漫反射法线与均匀的选择非常学术，除非（a）你担心共轭，在这种情况下你会使用法线或（b）真实值可能在均匀的端点之外有一些合理的机会. 使用 PyMC，没有理由担心共轭，除非您特别想使用 Gibbs 采样器。

  2. 对于方差/精度参数之前的无信息者，伽玛实际上不是一个很好的选择。它最终可能会比您认为的提供更多信息。更好的选择是对标准偏差进行统一的先验，然后通过平方反比对其进行变换。有关详细信息，请参阅Gelman 2006。

其它你可能感兴趣的问题
上一篇构成协方差矩阵的变量之间的距离是多少？ 下一篇由于样本量非常大，如何纠正小 p 值