Bivariate cumulative probit function

Hi folks,

I am coding up a “non-standard” likelihood which at some point involves the bivariate CDF of the normal function. That is, I want to make use of the greta equivalent of, say, the pmvnorm function in the R package mvtnorm, which calculates the CDF for a multivariate normal distribution given a mean vector and covariance matrix.

Is greta able to do this straightforwardly?

Thanks heaps!

In general it’s pretty straightforward to add most new link functions or other functions - just write out the analytic expression as you would in R. The problem with the (multivariate) CDF of a bivariate normal is that it doesn’t have an analytic expression. mvtnorm::pmvnorm uses a Monte Carlo simulation to approximate the value of the CDF.

That means it’s not differentiable, so can’t be used with HMC. It might be possible to use it with some other MCMC algorithm that doesn’t use gradient information. But in any case, to use it in greta would require a Tensorflow (or C, linked to Tensorflow) version of the Genz-Bretz algorithm they use to do the efficient Monte Carlo simulation. That’s doable, but very hard work, and you still wouldn’t be able to use HMC, so probably not worth doing in greta.

An alternative might be too recast your model to explicitly consider (rather than marginalise) the latent bivariate Gaussian. Not sure if that’s feasible though.

Thanks a lot for the advice and background Nick!

Much appreciated.

Yours non-significantly,
FH