1

我正在努力实现贝叶斯层次模型。

基本上，我正在尝试建立一个具有这种结构的模型......

分层模型图片

我有关于个人在将球扔给他们时能够击球多少次的数据（建模为二项式）。击中每个球的个人机会 (theta) 是从所有可能的 theta 的 beta 分布中得出的。这种个人级别的贝塔分布的模式本身是从更高级别的“玩家位置”贝塔分布中得出的，它代表了在该位置上玩的所有个人的可能模式。这个想法是将所有这些都放入一个层次模型中，以允许在单个级别的分布上进行收缩。

我已经尝试像这样使用 pymc 来实现它，但它似乎没有正确采样（即，所有东西都保持在起始值）。任何帮助表示赞赏！

import pandas as pd
import pymc as pm

### EXAMPLE DATA

initial_data=pd.DataFrame({'player_position':pd.Categorical([1,1,1,2,2,2]),
                       'player_ids':pd.Categorical([1, 2, 3, 4, 5, 6]),
                       'total_chances':pd.array([8, 7, 4, 4, 6, 13]), 
                       'hits':np.array([7, 7, 3, 3, 5, 11])})

### declarations
total_chances = initial_data.total_chances.values
hits = initial_data.hits.values
player_ids= initial_data.player_ids.values
position_ids = initial_data.player_position.values
num_players = len(initial_data.player_ids.unique())
num_positions = len(initial_data.player_position.unique())

### SPECIFYING THE MODEL
### overall-level distributions
overall_mode = pm.Beta('overall_mode', 1, 1)
overall_conc = pm.Gamma('overall_conc', 0.1, 0.1)

### position-level distributions
position_mode = pm.Beta('position_mode', (overall_mode*overall_conc)+1, 
                           ((1-overall_mode)*overall_conc)+1, 
                            size=num_positions)
position_conc = pm.Gamma('position_conc', 0.1, 0.1, 
                            size=num_positions)

### individual-level distribution 
individual_theta = pm.Beta('individual_theta',
     (position_mode[position_ids]*position_conc[position_ids])+1, 
                            ((1-
     position_mode[position_ids])*position_conc[position_ids])+1, size=num_players)

### likelihood scoring
hits_likelihood = pm.Binomial('hits_likelihood', n=total_chances, 
p=individual_theta, value=hits, observed=True)

### model construction
model = pm.Model([overall_mode, overall_conc, position_mode, 
position_conc, individual_theta, hits_likelihood])

### RUNNING THE MODEL
map_ = pm.MAP(model)
map_.fit()
mcmc = pm.MCMC(model)
mcmc.sample(2000, 400, 1)

pm.Matplot.plot(mcmc)
4

0 回答 0