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

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.