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

Questions
1

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.

2

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