具有抛硬币问题的简单贝叶斯网络

Question

我正在尝试实现贝叶斯网络并使用 PYMC3 解决回归问题。在我的模型中，我有一个公平的硬币作为父节点。如果父节点为H，则子节点选择正态分布N(5,0.2)；如果为 T，则孩子选择 N(0,0.5)。这是我的网络的插图。

为了模拟这个网络，我生成了一个示例数据集并尝试使用下面的代码进行贝叶斯回归。目前，模型只对子节点进行回归，就好像父节点不存在一样。如果有人能让我知道如何实现条件概率 P(D|C)，我将不胜感激。最终，我有兴趣找到 mu1 和 mu2 的概率分布。谢谢！

# Generate data for coin flip P(C) and store in c1
theta_real = 0.5 # unkown value in a real experiment
n_sample = 10
c1 = bernoulli.rvs(p=theta_real, size=n_sample)

# Generate data for normal distribution P(D|C) and store in d1
np.random.seed(123)
mu1 = 0
sigma1 = 0.5
mu2 = 5
sigma2 = 0.2

d1 = []
for index, item in enumerate(c1):
    if item == 0:
        d1.extend(normal(mu1, sigma1, 1))
    else:
        d1.extend(normal(mu2, sigma2, 1))

# I start building PYMC3 model here
c1_tensor = theano.shared(np.array(c1))
d1_tensor = theano.shared(np.array(d1))
with pm.Model() as model:
   # define prior for c1.  I am not sure how to do this.
   #c1_present = pm.Categorical('c1',observed=c1_tensor)

   # how do I incorporate P(D | C)
   mu_prior = pm.Normal('mu', mu=2, sd=2, shape=1)  
   sigma_prior = pm.HalfNormal('sigma', sd=2, shape=1)
   y_likelihood = pm.Normal('y', mu=mu_prior, sd=sigma_prior, observed=d1_tensor)

Answer 1

此答案是对@balleveryday's answer 的补充，它建议使用高斯混合模型，但在使对称性破坏起作用时遇到了一些麻烦。诚然，the official example 中的对称性破坏是在 Metropolis-Hastings 采样的背景下完成的，而我认为 NUTS 可能对遇到不可能的值更敏感（不确定）。这对我有用：

import numpy as np
import pymc3 as pm
from scipy.stats import bernoulli
import theano.tensor as tt

# everything should reproduce
np.random.seed(123)
n_sample = 2000

# Generate data for coin flip P(C) and store in c1
theta_real = 0.2 # unknown value in a real experiment

c1 = bernoulli.rvs(p=theta_real, size=n_sample)

# Generate data for normal distribution P(D|C) and store in d1
mu1, mu2 = 0, 5
sigma1, sigma2 = 0.5, 0.2

d1 = np.empty_like(c1, dtype=np.float64)
d1[c1 == 0] = np.random.normal(mu1, sigma1, np.sum(c1 == 0))
d1[c1 == 1] = np.random.normal(mu2, sigma2, np.sum(c1 == 1))

with pm.Model() as gmm_asym:
    # mixture vector
    w = pm.Dirichlet('p', a=np.ones(2))

    # Gaussian parameters (testval helps start off ordered)
    mu = pm.Normal('mu', 0, 20, shape=2, testval=[-10, 10])
    sigma = pm.HalfNormal('sigma', sd=2, shape=2)

    # break symmetry, forcing mu[0] < mu[1]
    order_means_potential = pm.Potential('order_means_potential',
                                         tt.switch(mu[1] - mu[0] < 0, -np.inf, 0))

    # observed
    pm.NormalMixture('like', w=w, mu=mu, sigma=sigma, observed=d1)

    # reproducible sampling
    tr_gmm_asym = pm.sample(tune=2000, target_accept=0.9, random_seed=20191121)

这会生成带有统计信息的样本

          mean      sd        mc_error  hpd_2.5   hpd_97.5  n_eff       Rhat
mu__0     0.004549  0.011975  0.000226 -0.017398  0.029375  2425.487301 0.999916
mu__1     5.007663  0.008993  0.000166  4.989247  5.024692  2181.134002 0.999563
p__0      0.789983  0.009091  0.000188  0.773059  0.808062  2417.356539 0.999788
p__1      0.210017  0.009091  0.000188  0.191938  0.226941  2417.356539 0.999788
sigma__0  0.497322  0.009103  0.000186  0.480394  0.515867  2227.397854 0.999358
sigma__1  0.191310  0.006633  0.000141  0.178924  0.204859  2286.817037 0.999614

和痕迹

Answer 2

您可以使用狄利克雷分布作为抛硬币的先验，并使用NormalMixture 作为两个高斯分布的先验。在下面的 sn-p 中，我改变了硬币的公平性并增加了掷硬币的次数，但你可以随意调整这些：

import numpy as np
import pymc3 as pm
from scipy.stats import bernoulli

# Generate data for coin flip P(C) and store in c1
theta_real = 0.2 # unkown value in a real experiment

n_sample = 2000
c1 = bernoulli.rvs(p=theta_real, size=n_sample)

# Generate data for normal distribution P(D|C) and store in d1
np.random.seed(123)
mu1 = 0
sigma1 = 0.5
mu2 = 5
sigma2 = 0.2

d1 = []
for index, item in enumerate(c1):
    if item == 0:
        d1.extend(np.random.normal(mu1, sigma1, 1))
    else:
        d1.extend(np.random.normal(mu2, sigma2, 1))

with pm.Model() as model:
   w = pm.Dirichlet('p', a=np.ones(2))
   mu = pm.Normal('mu', 0, 20, shape=2)
   sigma = np.array([0.5,0.2]) 
   pm.NormalMixture('like',w=w,mu=mu,sigma=sigma,observed=np.array(d1))
   trace = pm.sample()

pm.summary(trace) 

这将为您提供以下信息：

           mean        sd  mc_error   hpd_2.5  hpd_97.5        n_eff      Rhat
mu__0  4.981222  0.023900  0.000491  4.935044  5.027420  2643.052184  0.999637
mu__1 -0.007660  0.004946  0.000095 -0.017388  0.001576  2481.146286  1.000312
p__0   0.213976  0.009393  0.000167  0.195602  0.231803  2245.905021  0.999302
p__1   0.786024  0.009393  0.000167  0.768197  0.804398  2245.905021  0.999302

从跟踪图中也可以看出，参数恢复得很好：

上面的实现会给你theta_real、mu1和mu2的后验，但是当我添加sigma1和sigma2作为数据估计的参数时我无法收敛（即使先验范围很窄）：

with pm.Model() as model:
   w = pm.Dirichlet('p', a=np.ones(2))
   mu = pm.Normal('mu', 0, 20, shape=2)
   sigma = pm.HalfNormal('sigma', sd=2, shape=2)
   pm.NormalMixture('like',w=w,mu=mu,sigma=sigma,observed=np.array(d1))
   trace = pm.sample()


print(pm.summary(trace))
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [sigma, mu, p]
Sampling 4 chains: 100%|██████████| 4000/4000 [00:10<00:00, 395.57draws/s] 
The acceptance probability does not match the target. It is 0.883057127209148, but should be close to 0.8. Try to increase the number of tuning steps.
The gelman-rubin statistic is larger than 1.4 for some parameters. The sampler did not converge.
The estimated number of effective samples is smaller than 200 for some parameters.
              mean        sd  mc_error  ...  hpd_97.5     n_eff        Rhat
mu__0     1.244021  2.165433  0.216540  ...  5.005507  2.002049  212.596596
mu__1     3.743879  2.165122  0.216510  ...  5.012067  2.002040  235.750129
p__0      0.643069  0.248630  0.024846  ...  0.803369  2.004185   30.966189
p__1      0.356931  0.248630  0.024846  ...  0.798632  2.004185   30.966189
sigma__0  0.416207  0.125435  0.012517  ...  0.504110  2.009031   17.333177
sigma__1  0.271763  0.125539  0.012533  ...  0.497208  2.007779   19.217223

[6 rows x 7 columns]

基于此，如果您还想估计此数据的两个标准差，您很可能需要重新参数化。