PyMC3 Normal,每列有方差
PyMC3 Normal with variance per column
我正在尝试定义一个 pymc3.Normal 变量,其中包含以下 mu:
import numpy as np
import pymc3 as pm
mx = np.array([[0.25 , 0.5 , 0.75 , 1. ],
[0.25 , 0.333, 0.25 , 0. ],
[0.25 , 0.167, 0. , 0. ],
[0.25 , 0. , 0. , 0. ]])
epsilon = pm.Gamma('epsilon', alpha=10, beta=10)
p_ = pm.Normal('p_', mu=mx, shape = mx.shape, sd = epsilon)
问题是 p_ 中的所有随机变量都获得相同的标准差 (epsilon)。我希望第一行使用 epsilon1,第二行使用 epsilon2 等
我该怎么做?
可以为 shape
参数传递参数来实现此目的。为了证明这一点,让我们制作一些假数据以观察到的方式传递,其中我们使用 epsilon 的固定值,我们可以将其与推断的值进行比较。
示例模型
import numpy as np
import pymc3 as pm
import arviz as az
# priors
mu = np.array([[0.25 , 0.5 , 0.75 , 1. ],
[0.25 , 0.333, 0.25 , 0. ],
[0.25 , 0.167, 0. , 0. ],
[0.25 , 0. , 0. , 0. ]])
alpha, beta = 10, 10
# fake data
np.random.seed(2019)
# row vector will use a different sd for each column
sd = np.random.gamma(alpha, 1.0/beta, size=(1,4))
# generate 100 fake observations of the (4,4) random variables
Y = np.random.normal(loc=mu, scale=sd, size=(100,4,4))
# true column sd's
print(sd)
# [[0.90055471 1.24522079 0.85846659 1.19588367]]
# mean sd's per column
print(np.mean(np.std(Y, 0), 0))
# [0.92028042 1.24437592 0.83383181 1.22717313]
# model
with pm.Model() as model:
# use a (1,4) matrix to pool variance by columns
epsilon = pm.Gamma('epsilon', alpha=10, beta=10, shape=(1, mu.shape[1]))
p_ = pm.Normal('p_', mu=mu, sd=epsilon, shape=mu.shape, observed=Y)
trace = pm.sample(random_seed=2019)
这个样本很好,并给出了以下总结
这清楚地限制了 HPD 内标准偏差的真实值。
我正在尝试定义一个 pymc3.Normal 变量,其中包含以下 mu:
import numpy as np
import pymc3 as pm
mx = np.array([[0.25 , 0.5 , 0.75 , 1. ],
[0.25 , 0.333, 0.25 , 0. ],
[0.25 , 0.167, 0. , 0. ],
[0.25 , 0. , 0. , 0. ]])
epsilon = pm.Gamma('epsilon', alpha=10, beta=10)
p_ = pm.Normal('p_', mu=mx, shape = mx.shape, sd = epsilon)
问题是 p_ 中的所有随机变量都获得相同的标准差 (epsilon)。我希望第一行使用 epsilon1,第二行使用 epsilon2 等
我该怎么做?
可以为 shape
参数传递参数来实现此目的。为了证明这一点,让我们制作一些假数据以观察到的方式传递,其中我们使用 epsilon 的固定值,我们可以将其与推断的值进行比较。
示例模型
import numpy as np
import pymc3 as pm
import arviz as az
# priors
mu = np.array([[0.25 , 0.5 , 0.75 , 1. ],
[0.25 , 0.333, 0.25 , 0. ],
[0.25 , 0.167, 0. , 0. ],
[0.25 , 0. , 0. , 0. ]])
alpha, beta = 10, 10
# fake data
np.random.seed(2019)
# row vector will use a different sd for each column
sd = np.random.gamma(alpha, 1.0/beta, size=(1,4))
# generate 100 fake observations of the (4,4) random variables
Y = np.random.normal(loc=mu, scale=sd, size=(100,4,4))
# true column sd's
print(sd)
# [[0.90055471 1.24522079 0.85846659 1.19588367]]
# mean sd's per column
print(np.mean(np.std(Y, 0), 0))
# [0.92028042 1.24437592 0.83383181 1.22717313]
# model
with pm.Model() as model:
# use a (1,4) matrix to pool variance by columns
epsilon = pm.Gamma('epsilon', alpha=10, beta=10, shape=(1, mu.shape[1]))
p_ = pm.Normal('p_', mu=mu, sd=epsilon, shape=mu.shape, observed=Y)
trace = pm.sample(random_seed=2019)
这个样本很好,并给出了以下总结
这清楚地限制了 HPD 内标准偏差的真实值。