2

I am trying to recreate the models in John Kruschke's 'Doing Bayesian Data Analysis' and am currently trying to model ordinal data (chapter 23 in the book. This is the JAGS model that I'm trying to recreate:

total = length(y) 
#Threshold 1 and nYlevels-1 are fixed; other thresholds are estimated.
#This allows all parameters to be interpretable on the response scale.
nYlevels = max(y)  
thresh = rep(NA,nYlevels-1)
thresh[1] = 1 + 0.5
thresh[nYlevels-1] = nYlevels-1 + 0.5

modelString = "
model {
  for ( i in 1:Ntotal ) {
    y[i] ~ dcat( pr[i,1:nYlevels] )
    pr[i,1] <- pnorm( thresh[1] , mu[x[i]] , 1/sigma[x[i]]^2 )
    for ( k in 2:(nYlevels-1) ) {
      pr[i,k] <- max( 0 ,  pnorm( thresh[ k ] , mu[x[i]] , 1/sigma[x[i]]^2 )
                       - pnorm( thresh[k-1] , mu[x[i]] , 1/sigma[x[i]]^2 ) )
    }
    pr[i,nYlevels] <- 1 - pnorm( thresh[nYlevels-1] , mu[x[i]] , 1/sigma[x[i]]^2 )
  }
  for ( j in 1:2 ) { # 2 groups
    mu[j] ~ dnorm( (1+nYlevels)/2 , 1/(nYlevels)^2 )
    sigma[j] ~ dunif( nYlevels/1000 , nYlevels*10 )
  }
  for ( k in 2:(nYlevels-2) ) {  # 1 and nYlevels-1 are fixed, not stochastic
    thresh[k] ~ dnorm( k+0.5 , 1/2^2 )
  }
}

This model involves using the cumulative normal function (pnorm in R). I tried to do it as follows, using Scipy's norm.cdf:

nYlevels = df.Y.cat.categories.size
thresh = np.arange(1.5, nYlevels)

with pm.Model() as ordinal_model:  

  #priors for mu and sigma

  mu = pm.Normal('mu', (1+ nYlevels)/2.0, 1.0/(nYlevels**2))
  sig = pm.Uniform('sig', nYlevels/1000.0, nYlevels*10.0)

  #priors for the other parameters (one for each ordinal level, except the fixed lowest and highest levels)

  theta2 = pm.Normal('theta2',mu=thresh[1], sigma = (1/2)**2)
  theta3 = pm.Normal('theta3',mu=thresh[2], sigma = (1/2)**2)
  theta4 = pm.Normal('theta4',mu=thresh[3], sigma = (1/2)**2)
  theta5 = pm.Normal('theta5',mu=thresh[4], sigma = (1/2)**2)

  #probabilities for the ordinal levels

  prob_theta1 = norm(loc=mu, scale = (1 / (sig ** 2))).cdf(thresh[0])
  prob_theta2 = np.max([0, (norm(mu=mu, sigma = (1 / (sig ** 2))).cdf(theta2)) - (NormalDist(mu=mu, sigma=(1 / (sig ** 2))).cdf(thresh[0]))])
  prob_theta3 = np.max([0, (norm(mu=mu, sigma = (1 / (sig ** 2))).cdf(theta3)) - (NormalDist(mu=mu, sigma=(1 / (sig ** 2))).cdf(theta2))])
  prob_theta4 = np.max([0, (norm(mu=mu, sigma = (1 / (sig ** 2))).cdf(theta4)) - (NormalDist(mu=mu, sigma=(1 / (sig ** 2))).cdf(theta3))])
  prob_theta5 = np.max([0, (norm(mu=mu, sigma = (1 / (sig ** 2))).cdf(theta5)) - (NormalDist(mu=mu, sigma=(1 / (sig ** 2))).cdf(theta4))])
  prob_theta6 = 1 - (norm(mu=mu, sigma=(1 / (sig ** 2))).cdf(theta5))

  pr = np.array([prob_theta1,prob_theta2,prob_theta3,prob_theta4,prob_theta5,prob_theta6])

  y = pm.Categorical('y', p = pr, observed=df.Y.cat.codes.values)

However, when I run this, I get an error: 'TypeError: Unsupported dtype for TensorType: object'. It occurs at the line I use norm.cdf. Perhaps stats.scipy.norm cannot work with pymc3 objects?

So I was wondering how to use the cumulative normal function in a PyMC3 model, or how to make norm.cdf work here.

NotSinJicx
  • 21
  • 3
  • Thanks for updating. Yeah, PyMC3 RandomVariable objects are not immediately compatible with NumPy because they don't readily cast to numerics - they are nodes in a Theano compute graph. Instead, try to rewrite only using [Theano Tensor operations](https://theano-pymc.readthedocs.io/en/latest/library/tensor/index.html) and PyMC3 methods. It's probably worth noting that all PyMC3 distributions implement a `logcdf`. [The tutorial on distributions](https://docs.pymc.io/en/v3/Probability_Distributions.html) may also be useful to look at. Sorry I don't have time work out the answer explicitly. – merv Dec 14 '21 at 17:35
  • thanks, this helps, I'll try to figure it out. – NotSinJicx Dec 15 '21 at 10:14
  • 2
    I ran into this same issue trying to implement an item-response model where the characteristic curve was a normal CDF. I was able to accomplish it by using the Theano tensor operator `erf` operation like @merv suggested: `import theano.tensor as tt; with pm.Model(): b = Normal("b", mu=0, sigma=10); phat = .5 * (1 + tt.erf((b - 0) / np.sqrt(2))); Y = Bernoulli("Y", p=phat, observed=1)` Hope that's helpful! – sjplural Dec 16 '21 at 22:41

1 Answers1

2

One way to do this is to use a Theano operator. Follow the example here. In that example, an ODE solution has been embeded in the model. The as_op decorator approach is very easy to use. It is also very flexible, and you can use SciPy's functions. But this approach does limit the type of samplers you may use to draw samples. Speedwise, it is a bit slower.

If you care about run-time speed, you can supply custom gradient functions. See here for an example. But programming time will be much longer, and this approach is prone to programming errors.

Personally, I have only used the decorator approach, and I recommend you to do the same for your model.

hyena
  • 318
  • 1
  • 9