# This code is maintained at: # # https://github.com/AlexanderLyNL/jasp-desktop/blob/development/JASP-Engine/JASP/R/correlationbayesian.R # # and is part of the JASP software package: # # https://jasp-stats.org/ # # ----------------------------------- # Function definitions for the Marsman IMH sampler # An example of its usages is given below # .stretchedBeta <- function(rho, alpha, beta){ # Beta density stretched over -1, 1 # Use: Fit for the prior and posterior result <- 1/2*dbeta((rho+1)/2, alpha, beta) return(result) } .betaParameterEstimates <- function(someMean, someVar){ # Use: To fit a stretched beta to the posterior samples given some mean and variance # Note: someMean \in (0, 1) # some.a <- someMean*(someMean*(1-someMean)/someVar-1) some.b <- (1-someMean)*(someMean*(1-someMean)/someVar-1) result <- list(betaA=some.a, betaB=some.b) return(result) } .logTarget <- function(z, n, r, alpha){ # Auxilary definition: For the IMH sampler # z is Fisher's transformation for r, but also use it for rho # The Fisher z transform and the log (likelihood*prior*Jacobian) of the tranformation (0.5*(n - 1))*log(1 - tanh(z)^2) - (n - 1 - 0.5)*log(1 - tanh(z)*r)+log(1-tanh(z)^2)*alpha } .logProposal <- function(z, n, r){ # Auxilary definition: For the IMH sampler # z is Fisher's transformation for r, but also use it for rho -(n-3)/2*(z-atanh(r))^2 } .independentOneStep <- function (rhoCurrent, n, r, alpha=1) { # One step sampler zCurrent <- atanh(rhoCurrent) zCandidate <- rnorm(1, mean=atanh(r), sd=1/sqrt(n-3)) candidateAcceptance <- .logTarget(zCandidate, n, r, alpha)+ .logProposal(zCurrent, n, r)- (.logTarget(zCurrent, n, r, alpha)+ .logProposal(zCandidate, n, r)) if (log(runif (1)) <= candidateAcceptance) { rhoCandidate <- tanh(zCandidate) return(rhoCandidate) } else { return (rhoCurrent) } } .marsmanIMHSampler <- function(n, r, alpha=1, nIters=50000){ # Use: Posterior samples rhoMetropolisChain <- NULL yTemp <- r for (iter in 1:nIters) { yTemp <- .independentOneStep(yTemp, n, r, alpha) rhoMetropolisChain[iter] <- yTemp } acceptanceRate <- length (unique (rhoMetropolisChain)) / nIters metropolisVar <- var(rhoMetropolisChain)/2^2 metropolisMean <- mean((rhoMetropolisChain+1)/2) mhFit <- .betaParameterEstimates(metropolisMean, metropolisVar) mhFit$acceptanceRate <- acceptanceRate mhFit$samples <- rhoMetropolisChain return(mhFit) } # ----------------------------------- # Example: # somePosterior <- .marsmanIMHSampler(24, 0.6) somePosterior$acceptanceRate plot(density(somePosterior$samples), xlim=c(-1, 1), type="l", lwd=2) rhoDomain <- seq(-1, 1, by=0.001) lines(rhoDomain, .stretchedBeta(rhoDomain, somePosterior$betaA, somePosterior$betaB), col="red")