gistfile1.r

dHoyt <- function(x, qpar, omega) {
    fac1 <- x*(1+qpar^2) / (qpar*omega)
    fac2 <- exp( -x^2*(1+qpar^2)^2/(4*qpar^2*omega))
    besselArg <- (x^2*(1-qpar^4)  /(4*qpar^2*omega))
    fac3 <- exp(log(besselI(besselArg, nu=0, expon.scaled=TRUE)) + besselArg)
    return(fac1*fac2*fac3)
}

simR <- function(n, qpar, omega) {
    ev2  <- omega / ((1/qpar^2) + 1) # 1st eigenvalue
    ev1  <- omega - ev2              # 2nd eigenvalue
    eVal <- sort(c(ev1, ev2), decreasing=TRUE)  # sorted eigenvalues

    ## simulated 2D normal vectors with mean 0 with identity cov-mat
    X  <- matrix(rnorm(2*n), nrow=n)
    xy <- X %*% diag(sqrt(eVal), 2)
    return(sqrt(rowSums(xy^2)))      # distances to center
}

## simulate radii
r <- simR(10000, qpar=0.5, omega=10)

## distribution of simulated radii vs. Hoyt density
hist(r, freq=FALSE, breaks="FD")
curve(dHoyt(x, qpar=0.5, omega=10), add=TRUE, col="blue")

## distribution of squared simulated radii vs. Gamma density
hist(r^2, freq=FALSE, breaks="FD")
curve(dgamma(x, shape=0.5, scale=0.5/10), add=TRUE, col="blue")
## maybe parametrization is different?
curve(dgamma(x, shape=0.5, rate=0.5/10), add=TRUE, col="blue")