创建(正常)平均值和标准差的两层估计
创建(正常)平均值和标准差的两层估计
提问于 2021-05-03 21:06:32
我有一个正态分布的变量x(如产品需求)，一个索引id_1 (如产品编号)和第二个索引id_2 (如产品组)。我的目标是按层次估计x的平均值和标准差(全部>产品组>产品)。这是我的数据：
import numpy as np
import pymc3 as pm
import arviz as az
# data
my_step_1 = 0.4
my_step_2 = 4.1
sd_step_1 = 0.1
sd_step_2 = 0.2
my = 10
sd = .1
grp_n = 8
grps_1 = 5
grps_2 = 4
x = np.round(np.concatenate([np.random.normal(my + i * my_step_1 + j * my_step_2, \
sd + i * sd_step_1 + j * sd_step_2, grp_n) \
for i in range(grps_1) for j in range(grps_2)]), 1) # demand
id_1 = np.repeat(np.arange(grps_1 * grps_2), grp_n) # group, product number
id_2 = np.tile(np.repeat(np.arange(grps_2), grp_n), grps_1) # super-group, product group
shape_1 = len(np.unique(id_1))
shape_2 = len(np.unique(id_2))我只管理了一个层次结构：
with pm.Model() as model_h1:
#
mu_mu_hyper = pm.Normal('mu_mu_hyper', mu = 0, sd = 10)
mu_sd_hyper = pm.HalfNormal('mu_sd_hyper', 10)
sd_hyper = pm.HalfNormal('sd_hyper', 10)
#
mu = pm.Normal('mu', mu = mu_mu_hyper, sd = mu_sd_hyper, shape = shape_1)
sd = pm.HalfNormal('sd', sd = sd_hyper, shape = shape_1)
y = pm.Normal('y', mu = mu[id_1], sd = sd[id_1], observed = x)
trace_h1 = pm.sample(1000)
#az.plot_forest(trace_h1, var_names=['mu', 'sd'], combined = True)但是，我如何编写2个层次结构呢？
# 2 hierarchies .. doesn't work
with pm.Model() as model_h2:
#
mu_mu_hyper2 = pm.Normal('mu_mu_hyper2', mu = 0, sd = 10)
mu_sd_hyper2 = pm.HalfNormal('mu_sd_hyper2', sd = 10)
sd_mu_sd_hyper2 = pm.HalfNormal('sd_mu_sd_hyper2', sd = 10)
sd_hyper2 = pm.HalfNormal('sd_hyper2', sd = 10)
#
mu_mu_hyper1 = pm.Normal('mu_hyper1', mu = mu_mu_hyper2, sd = mu_sd_hyper2, shape = shape_2)
mu_sd_hyper1 = pm.HalfNormal('mu_sd_hyper1', sd = sd_mu_sd_hyper2, shape = shape_2)
sd_hyper1 = pm.HalfNormal('sd_hyper1', sd = sd_hyper2, shape = shape_2)
#sd_hyper1 = pm.HalfNormal('sd_hyper1', sd = sd_hyper2[id_2], shape = shape_2)??
#
mu = pm.Normal('mu', mu = mu_mu_hyper1, sd = mu_sd_hyper1, shape = shape_1)
sd = pm.HalfNormal('sd', sd = sd_hyper1, shape = shape_1)
y = pm.Normal('y', mu = mu[id_1], sd = sd[id_1], observed = x)
trace_h2 = pm.sample(1000)回答 1
Stack Overflow用户
发布于 2021-08-06 07:40:09
您可以尝试遍历产品组，并使用组的均值、std作为属于此特定组的产品的约束。
# sample product group to product mapping
group_product_mapping = {0: [1, 2, 3], 1: [4, 5, 6]}
total_groups = len(group_product_mapping.keys())
with pm.model() as model_h:
mu_all = pm.Normal('mu_all', 0, 10)
sd_all = pm.HalfNormal('sd_all', 10)
sd_mu_group = pm.HalfNormal('sd_mu_group', 10)
# group parameters constrained to mu, sd from all
mu_group = pm.Normal('mu_group', mu_all, sd_all, shape=total_groups)
sd_group = pm.HalfNormal('sd_group', sd_mu_group, shape=total_groups)
mu_products = dict()
# iterate through groups and constrain product parameters to the product group they belong to
for idx, group in enumerate(group_product_mapping.keys()):
mu_products[group] = pm.Normal(f'mu_products_{group}', mu_group[idx], sd_group[idx], shape=len(group_product_mapping[group]))
sd_produtcs[group] = pm.HalfNormal(f'sd_mu_products_{group}', 10)页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/67369450
复制相关文章
相似问题
腾讯云开发者
