I’m trying to formulate a simple example in PyMC3 of using a Multinomial distribution to infer the probabilities of a die throw, except where the side probabilities are linked by a specific model with some lesser number of parameters. In the example below it is just a linearly decreasing probability with some characteristic slope m, which is the parameter to infer. I’ve got it working using the code below by setting up a Uniform random variable for the prior and then making a Deterministic variable for each side probability. But I don’t really like having all the terms written out individually, especially if one would extend the number of possible states it would get very ugly. I was thinking of using a custom DensityDist to specify the probabilities, but I can’t quite figure out if that is feasible or not. Or possibly I could write a custom distribution, but also it is not readily apparent how to do that. Does anyone have any insight on how to improve the code below? Thanks!
import numpy as np import pymc3 as mc def randomDataLinearDiceModel(p, n): # Calculate the cumulative probabilities. c = np.cumsum(p) # Generate a random sequence of dice rolls. r = np.random.rand(n) x = np.zeros((n,), dtype=np.uint8) for i in range(n): x[i] = np.argmax(r[i]<c) # Get the counts for each die side. (counts,edges) = np.histogram(x,np.arange(0-0.5,len(p)+0.5,1.0)) centers = (edges[:-1]+edges[1:])/2 return centers,counts def bayesianInferenceLinearDiceModel(p, n, counts): with mc.Model() as model: # Prior distributions. m = mc.Uniform('m', lower=0, upper=0.05) p0 = mc.Deterministic('p0', (1/6) + 3*m) p1 = mc.Deterministic('p1', (1/6) + 2*m) p2 = mc.Deterministic('p2', (1/6) + 1*m) p3 = mc.Deterministic('p3', (1/6) - 1*m) p4 = mc.Deterministic('p4', (1/6) - 2*m) p5 = mc.Deterministic('p5', (1/6) - 3*m) # Likelihood. counts_likelihood = mc.Multinomial('counts', n=n, p=[p0,p1,p2,p3,p4,p5], observed=counts) # Perform the sampling. trace = mc.sample(2000, tune=1000) return trace m=0.045 p = np.array([1/6+3*m,1/6+2*m,1/6+1*m,1/6-1*m,1/6-2*m,1/6-3*m]) n=10000 sides,rolls = randomDataLinearDiceModel(p, n) trace = bayesianInferenceLinearDiceModel(p, n, rolls) print("Posterior distribution summary:") print(mc.summary(trace)) print("Finished.")