Bad initial energy, check any log probabilities that are inf or -inf, nan or very small:

I am trying to implement the following research paper https://people.ok.ubc.ca/pgill/research/cricketsim.pdf which is a simulation model for cricket matches. See section 2.1 for model specification (“a Bayesian latent variable model which is
related to the classical cumulative-logit models for ordinal responses”).

I first implemented the model in WinBUGS like the authors of the paper, but it didn’t seem to run nearly fast enough. Here is how I specified the model in WinBUGS. I am not sure if this is correct, but it seems to be a direct implementation of the paper.

> model
> {
> 	for(l in 1:9){
> 		z[l,1] <- - 999999999
> 		for(m in 2:7){
> 			z[l,m] ~ dnorm(0,s)
> 		}
> 		z[l,8] <- 999999999
> 
> 		for(o in 1:8){
> 			a[l,o] <- ranked(z[l,],o)
> 		}
> 		
> 		d[l] <- 0 + step(l-4) * delta1 + step(l-7) * delta2
> 	}	
> 		for(k in 1:7){
> 			p[,k] <- 1/(1+exp(-(a[L,k+1] - u1[batsmanid[n]] + u2[bowlerid[n]] - d[L]))) - 1/(1+exp(-(a[L,k] - u1[batsmanid] + u2[bowlerid]-d[L])))
> 		}
> 		X ~ dcat(p)
> 	
> 	for(i in 1:I){
> 		u1[i] ~ dnorm(0,t)
> 	}
> 	
> 	for(j in 1:J){
> 		u2[j] ~ dnorm(0,t)
> 	}
> 	delta1 ~ dunif(0,1)
> 	delta2 ~ dunif(0,1)	
> 	t ~ dgamma(1,1)
>  s ~ dgamma(1,1)
> }

Now I am trying to implement it with PyMC3, but I am getting the following error:

Bad initial energy, check any log probabilities that are inf or -inf, nan or very small:

I faced problems in removing the for-loops and vectorising the implementation, which seemed challenging and I am not sure if the way I have done it is okay. Here is my implementation.

> INF = 5
> testval = [[-INF + x * (2 * INF)/7.0 for x in range(0, 8)] for i in range(0, 9)]
> l = [i for i in range(1, 10)]

> model = pm.Model()
> with model:
>     delta_1 = pm.Uniform("delta_1", lower=0, upper=1)
>     delta_2 = pm.Uniform("delta_2", lower=0, upper=1)
>     inv_sigma_sqr = pm.Gamma("sigma^-2", alpha=1.0, beta=1.0)
>     inv_tau_sqr = pm.Gamma("tau^-2", alpha=1.0, beta=1.0)
>     a = pm.Normal("a", mu=0, sigma=1/pm.math.sqrt(inv_sigma_sqr), transform=tr.Chain([tr.LogOdds(), tr.Ordered()]), shape=(9,8), testval=testval)
>     mu_1 = pm.Normal("mu_1", mu=0, sigma=1/pm.math.sqrt(inv_tau_sqr), shape=N)
>     mu_2 = pm.Normal("mu_2", mu=0, sigma=1/pm.math.sqrt(inv_tau_sqr), shape=N)
>     delta = pm.math.ge(l, 4) * delta_1 + pm.math.ge(l, 4) * delta_2
>     p = invlogit(a[L-1] - stack([mu_1[id1]] * 8, axis=1) + stack([mu_2[id2]] * 8, axis=1) - stack([delta[L-1]] * 8, axis=1)) - \
>         invlogit(concatenate([np.zeros(shape=(N,1)), a[L-1][:,:-1]], axis=1) - stack([mu_1[id1]] * 8, axis=1) + stack([mu_2[id2]] * 8, axis=1) - stack([delta[L-1]] * 8, axis=1))
>     X_obs = pm.Categorical("X_obs", p, observed=X-1)

I am a complete newbie to Bayesian inference and pymc3 so please help me out :frowning:
Thank you!

EDIT: I just discovered the OrderedLogistic distribution which I think is what I should be using. Any ideas how to implement the model using that since my predictors can have different sets of “cutpoints” depending on their value.

Yep you should definitely use OrderedLogistic. There is a simple example in the docstring : https://docs.pymc.io/api/distributions/discrete.html#pymc3.distributions.discrete.OrderedLogistic, otherwise, the discussion below might also helpful to you:
Ordered probit model for ordinal data

1 Like