Skip to content

Instantly share code, notes, and snippets.

/gist:9094590

Created Feb 19, 2014
Embed
What would you like to do?
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")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment