Python -- PYMC beta-binomial with shrinkage

Question

I'm struggling with the implementation of a Bayesian Hierarchical model.

Basically, I'm trying to build a model with this structure...

I've got data for how many times individuals are able to hit a ball when it's thrown to them (modeled as a binomial). The individual chance to hit each ball (theta) is drawn from a beta distribution of all possible thetas. The mode of this individual-level beta distribution is itself drawn from a higher-level "player-position" beta distribution that represents the possible modes for all individuals that play in that position. The idea is to put this all into a hierarchical model, to allow shrinkage at the individual-level distributions.

I've tried implementing it using pymc like this, but it doesn't seem to be sampling properly (i.e. remains at starting values for everything). Any help appreciated!

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)

Python -- PYMC beta-binomial with shrinkage

0 Answers0