I'm working on a problem with PyMC3 that makes me think I need to better understand how it deals with random variables whose parameters are vector-valued.
Data description and problem setup
I have $N$ samples of binomial data, $(n_i, k_i)$ where $i = 1...N$, that is, $N$ experiments where for the $i$th experiment $n_i$ trials were performed and $k_i$ of those were successes (so $0 \leq k_i \leq n_i$).
For each of the $N$ data points, I also have a feature $x_i$, which happens to be a real number between 0 and 1. But it's known for each of the $N$ data points and fixed.
I'd like to model the binomial data with a Binomial distribution, and then model the distribution's probability parameter $p$ as Beta-distributed, which is the Binomial's conjugate prior.
I have reason to believe that the Beta prior's parameters, $\alpha = m_\alpha x + b_\alpha$ and $\beta = m_\beta x + b_\beta$, that is, a parameter of the Beta is a linear function of $x_i$.
So each of my binomial data points can be said to $$ k_i \sim Binomial(n_i, p_i), $$ where $$p_i \sim Beta(\alpha_i, \beta_i)$$ and where $$\alpha_i = m_\alpha x_i + b_\alpha$$ and $$\beta_i = m_\beta x_i + b_\beta$$
Synthetic data generation
The following generates some sample data:
import numpy as np from scipy.stats import beta, binom Ndata = 1000 feature = np.random.rand(Ndata) # x's trials = np.random.randint(1, 20, (Ndata,)) # n's alphaSlope = 10.0 # m_alpha alphaIntercept = 5.0 # b_alpha betaSlope = -3.0. # m_beta betaIntercept = 8.0 # b_beta obs = [binom(n, beta(al, be).rvs(1)).rvs(1) for al, be, n in zip( alphaSlope * feature + alphaIntercept, betaSlope * feature + betaIntercept, trials)] # k's
PyMC model 1
I model this in PyMC (version 3.8) as:
import theano.tensor as tt import pymc3 as pm with pm.Model() as model: # weights: slope, intercept walpha = pm.Uniform('walpha', lower=-20, upper=20, shape=2) wbeta = pm.Uniform('wbeta', lower=-20, upper=20, shape=2) # parameters of Beta a = (walpha * feature + walpha).clip(1, 5000) b = (wbeta * feature + wbeta).clip(1, 5000) # prior and likelihood p = pm.Beta('p', alpha=a, beta=b, testval=0.7) x = pm.Binomial('x', n=trials, p=p, observed=obs) trace = pm.sample(10000, tune=2000, cores=5)
(Before proceeding, I understand I could have modeled this using PyMC's
BetaBinomial random variable. I don't want to use this because in my actual problem, the Binomial isn't parameterized by probability $p$ but rather $p^\delta$, a nonlinear exponentiation.)
But the above PyMC model produces incorrect results: using Pandas to describe the MCMC trace with
pm.trace_to_dataframe(trace).describe() reveals both $\alpha$ and $\beta$'s slope and intercept parameters to be far from true values:
walpha__0 walpha__1 wbeta__0 wbeta__1 p count 50000.000000 50000.000000 50000.000000 50000.000000 50000.000000 mean 19.870350 19.942757 13.847049 14.402839 0.589303 std 0.128837 0.057044 0.725568 0.465387 0.005422 min 18.562779 19.464832 10.930696 12.495120 0.563810 25% 19.819790 19.920913 13.356722 14.090275 0.585660 50% 19.909900 19.960104 13.842200 14.396608 0.589340 75% 19.962834 19.983450 14.335798 14.713298 0.593011 max 19.999999 19.999999 17.108659 16.592652 0.612058
PyMC's trace plot shows how bad the results are (
What exactly is the line,
p = pm.Beta('p', alpha=a, beta=b, testval=0.7), doing? Note that
b here are vector-valued (1000 by 1, since I have 1000 points of data), each a linear combination of two random variables, so perhaps I need to specify the
p manually? This is my second PyMC model:
PyMC model 2
This is the exact same model as above except I add a
shape parameter to the Beta random variable, explicitly stating that it should be a 1000-long vector of a thousand different random variables:
with pm.Model() as model: # weights: slope, intercept walpha = pm.Uniform('walpha', lower=-20, upper=20, shape=2) wbeta = pm.Uniform('wbeta', lower=-20, upper=20, shape=2) # parameters of Beta a = (walpha * feature + walpha).clip(1, 5000) b = (wbeta * feature + wbeta).clip(1, 5000) # prior and likelihood p = pm.Beta('p', alpha=a, beta=b, testval=0.7, shape=len(obs)) # !!! x = pm.Binomial('x', n=trials, p=p, observed=obs) trace = pm.sample(10000, tune=2000, cores=5)
This model runs much slower (~200 draws per second, versus model 1's ~1500 draws per second), but produces reasonable results:
walpha__0 walpha__1 wbeta__0 wbeta__1 p__0 \ count 50000.000000 50000.000000 50000.000000 50000.000000 50000.000000 mean 13.071932 6.050287 -3.331585 9.634910 0.579830 std 2.793833 1.244391 1.643823 1.603906 0.090586 min 2.112455 2.620591 -11.579191 5.225143 0.208103 25% 11.097158 5.167857 -4.386282 8.506066 0.518849 50% 12.977088 5.926347 -3.280822 9.479486 0.581282 75% 14.968447 6.794506 -2.206138 10.582560 0.642689 max 19.995363 13.397899 2.877322 18.900355 0.880307
The true values for
walpha__0 walpha__1 wbeta__0 wbeta__1 are, recall, 10, 5, -3, and 8. (I'm a bit surprised that with 1000 data samples, and 50'000 MCMC iterations, it couldn't get closer to the true
walpha__0 so maybe there's something weird going on.)
Model 3: BetaBinomial
As mentioned above, I can't use PyMC's
BetaBinomial random in my actual application because I need a nonlinearly-transformed Beta prior on the Binomial parameter, but for just checking the above result, I tried it:
with pm.Model() as model: # weights: slope, intercept walpha = pm.Uniform('walpha', lower=-20, upper=20, shape=2) wbeta = pm.Uniform('wbeta', lower=-20, upper=20, shape=2) # parameters of Beta a = (walpha * feature + walpha).clip(1, 5000) b = (wbeta * feature + wbeta).clip(1, 5000) # prior and likelihood x = pm.BetaBinomial('x', observed=obs, n=trials, alpha=a, beta=b) trace = pm.sample(12000, tune=6000, cores=5)
BetaBinomial, PyMC runs at ~600 draws per second, faster than model 2 above but slower than model 1. I didn't need to specify any
shape parameter but it appears to have found the same results as model 2 above:
walpha__0 walpha__1 wbeta__0 wbeta__1 count 50000.000000 50000.000000 50000.000000 50000.000000 mean 13.049388 6.060485 -3.359202 9.649070 std 2.737753 1.250730 1.657210 1.600937 min 3.351191 2.292545 -11.563000 4.695184 25% 11.122850 5.171321 -4.418419 8.514347 50% 12.953331 5.943487 -3.293241 9.501761 75% 14.914511 6.806098 -2.238418 10.626369 max 19.999812 13.462526 3.996577 18.123882
My main question is, am I correctly specifying my model by including a
shape parameter in
p = pm.Beta('p', alpha=a, beta=b, testval=0.7, shape=len(obs))? I think this is doing the right thing because an explicit
BetaBinomial random variable returns extremely similar results, but I also think this could be improved because it runs so much slower than
BetaBinomial (which I cannot use in my actual example). Is there a better way of expressing this?
Secondarily, what exactly does PyMC do with
p = pm.Beta('p', alpha=a, beta=b, testval=0.7), that is, model 1's random variable where the
beta parameters are vectors but
p is apparently a scalar random variable? What kind of model is that statistically fitting? Is it perhaps just looking at the first elements of the
beta vectors, thereby giving weird results for the slope and intercept terms?
Any intuition about how PyMC handles random variables with vector parameters, and the impact of the
shape argument in those cases, would be helpful.