#library(MASS)
library(compiler)

# returns rho/sqrt(1-rho^2) rather than rho
# a useful transformation as discussed in
# Johnson, Klotz and Balakrishnan
# Continuous Univariate Distributions
# x, which is the current value, is assumed
# to have come from the same transformation

corrMH <- cmpfun(function(a, b, N, x) {
  p <- NaN

  while (is.nan(p)) {
    y <- rcauchy(1, location=x, scale=a^(-0.5))

    p <- ((1+y^2)/(1+x^2))^(0.5*N-1.5) *
      exp(b*(y*sqrt(1+y^2)-x*sqrt(1+x^2))-0.5*a*(y^2-x^2))
  }
         
  if (runif(1)<p) return(y)
  else return(x)
})

