Nothing
#' MCMC for rotation data
#'
#' Use non-informative Bayesian methods to infer about the central orientation and concentration
#' parameter for a sample of rotations.
#'
#' The procedures detailed in \cite{bingham2009b} and \cite{bingham2010} are implemented to obtain
#' draws from the posterior distribution for the central orientation and concentration parameters for
#' a sample of 3D rotations. A uniform prior on SO(3) is used for the central orientation and the
#' Jeffreys prior determined by \code{type} is used for the concentration parameter.
#'
#' @param x \eqn{n\times p}{n-by-p} matrix where each row corresponds to a random rotation in matrix (\eqn{p=9}) or quaternion (\eqn{p=4}) form.
#' @param type Angular distribution assumed on R. Options are \code{\link{Cayley}}, \code{\link{Fisher}} or \code{\link{Mises}}
#' @param S0 initial estimate of central orientation
#' @param kappa0 initial estimate of concentration parameter
#' @param tuneS central orientation tuning parameter, concentration of proposal distribution
#' @param tuneK concentration tuning parameter, standard deviation of proposal distribution
#' @param burn_in number of draws to use as burn-in
#' @param m number of draws to keep from posterior distribution
#' @return list of
#' \itemize{
#' \item \code{S} Draws from the posterior distribution for central orientation S
#' \item \code{kappa} Draws from the posterior distribution for concentration parameter kappa
#' \item \code{Saccept} Acceptance rate for central orientation draws
#' \item \code{Kaccept} Acceptance rate for concentration draws
#' }
#' @details bingham2009b bingham2010
#' @export
#' @examples
#' #Not run due to time constraints
#' \donttest{
#' Rs <- ruars(20, rfisher, kappa = 10)
#' draws <- MCMCSO3(Rs, type = "Fisher", S0 = mean(Rs), kappa0 = 10, tuneS = 5000,
#' tuneK = 1,burn_in = 1000, m = 5000)}
MCMCSO3<-function(x,type,S0,kappa0,tuneS,tuneK,burn_in,m=5000){
UseMethod("MCMCSO3")
}
#' @rdname MCMCSO3
#' @export
MCMCSO3.SO3<-function(x,type,S0,kappa0,tuneS,tuneK,burn_in,m=5000){
if(type %in% c("Cayley","cayley")){
lpangle <- 1
}else if(type %in% c("Fisher","fisher")){
lpangle <- 2
}else if(type %in% c("Mises","mises")){
lpangle <- 3
}else{
stop("Invalid choise of type: please choose Cayley, Fisher or Mises.")
}
listRes<-both_MCMC_CPP(x,S0, kappa0,tuneS,tuneK,burn_in,m, lpangle)
class(listRes$S)<-"SO3"
#listRes$S<-as.SO3(listRes$S)
return(listRes)
}
#' @rdname MCMCSO3
#' @export
MCMCSO3.Q4<-function(x,type,S0,kappa0,tuneS,tuneK,burn_in,m=5000){
Rs<-as.SO3(x)
S0<-as.SO3(matrix(S0,3,3))
SO3Res<-MCMCSO3(Rs,type,S0,kappa0,tuneS,tuneK,burn_in,m)
Q4Res<-list(Q=as.Q4(SO3Res$S),kappa=SO3Res$kappa,Qaccept=SO3Res$Saccept,Kaccept=SO3Res$Kaccept)
return(Q4Res)
}
#' Bayes credible regions
#'
#' Find the radius of a \eqn{100(1-\alpha)}\% credible region for the central orientation and concentration parameter using
#' non-informative Bayesian methods.
#'
#' Compute the radius of a \eqn{100(1-\alpha)}\% credible region for the central orientation and concentration parameter
#' as described in \cite{bingham2009b} and \cite{bingham2010}. The posterior mode is returned along with the radius
#' of the credible region centered at the posterior mode.
#'
#' @param x \eqn{n\times p}{n-by-p} matrix where each row corresponds to a random rotation in matrix (\eqn{p=9}) or quaternion (\eqn{p=4}) form.
#' @param type Angular distribution assumed on R. Options are \code{\link{Cayley}}, \code{\link{Fisher}} or \code{\link{Mises}}
#' @param S0 initial estimate of central orientation
#' @param kappa0 initial estimate of concentration parameter
#' @param tuneS central orientation tuning parameter, concentration of proposal distribution
#' @param tuneK concentration tuning parameter, standard deviation of proposal distribution
#' @param burn_in number of draws to use as burn-in
#' @param m number of draws to keep from posterior distribution
#' @param alp alpha level desired, e.g. 0.05 or 0.10.
#' @return list of
#' \itemize{
#' \item \code{Shat,Qhat} Mode of the posterior distribution for the central orientation S
#' \item \code{Radius} Radius of the credible region centered at the posterior mode
#' }
#' @seealso \code{\link{fisheretal}}, \code{\link{prentice}}, \code{\link{chang}}, \code{\link{zhang}}
#' @details bingham2009b bingham2010
#' @export
#' @examples
#' #Not run due to time constraints
#' \donttest{
#' Rs <- ruars(20, rvmises, kappa = 10)
#'
#' #Compare the region size of the moment based theory mean estimator to the
#' #Bayes region.
#'
#' region(Rs, method = "direct", type = "theory", estimator = "mean", alp=0.1, m = 100)
#' bayesCR <- region(Rs, type = "Mises", method = "Bayes", estimator = "mean", S0 = mean(Rs),
#' kappa0 = 10, tuneS = 5000, tuneK = 1, burn_in = 1000, alp = .01, m = 5000)
#'
#' bayesCR$Radius #Region size is give by "Radius"
#' bayesCR$Shat #The Bayes region is centered around the posterior mode: "Shat"}
bayesCR<-function(x,type,S0,kappa0,tuneS,tuneK,burn_in,m=5000,alp=0.1){
UseMethod("bayesCR")
}
#' @rdname bayesCR
#' @export
bayesCR.SO3<-function(x,type,S0,kappa0,tuneS,tuneK,burn_in,m=5000,alp=0.1){
if(type %in% c("Cayley","cayley")){
#lpangle <- lpcayley
lpangle <- 1
}else if(type %in% c("Fisher","fisher")){
#lpangle <- lpfisher
lpangle <- 2
}else if(type %in% c("Mises","mises")){
#lpangle <- lpvmises
lpangle <- 3
}else{
stop("Invalid choise of type: please choose Cayley, Fisher or Mises.")
}
listRes<-both_MCMC_CPP(x,S0, kappa0,tuneS,tuneK,burn_in,m, lpangle)
Sdraws<-as.SO3(listRes$S)
Shat<-mean(Sdraws)
rs<-rot.dist(Sdraws,Shat)
return(list(Shat=Shat,Radius=stats::quantile(rs,1-alp)))
}
#' @rdname bayesCR
#' @export
bayesCR.Q4<-function(x,type,S0,kappa0,tuneS,tuneK,burn_in,m=5000,alp=0.1){
Rs<-as.SO3(x)
S0<-as.SO3(matrix(S0,3,3))
SO3Res<-bayesCR(Rs,type,S0,kappa0,tuneS,tuneK,burn_in,m,alp)
Q4Res<-list(Qhat=as.Q4(SO3Res$Shat),Radius=SO3Res$Radius)
return(Q4Res)
}
#' Parameter estimates based on non-informative Bayes
#'
#' Use non-informative Bayes to estimate the central orientation and concentration parameter of a sample of rotations.
#'
#' The procedures detailed in \cite{bingham2009b} and \cite{bingham2010} are implemented to obtain
#' draws from the posterior distribution for the central orientation and concentration parameters for
#' a sample of 3D rotations. A uniform prior on SO(3) is used for the central orientation and the
#' Jeffreys prior determined by \code{type} is used for the concentration parameter.
#'
#' @param x \eqn{n\times p}{n-by-p} matrix where each row corresponds to a random rotation in matrix (\eqn{p=9}) or quaternion (\eqn{p=4}) form.
#' @param type Angular distribution assumed on R. Options are \code{\link{Cayley}}, \code{\link{Fisher}} or \code{\link{Mises}}
#' @param S0 initial estimate of central orientation
#' @param kappa0 initial estimate of concentration parameter
#' @param tuneS central orientation tuning parameter, concentration of proposal distribution
#' @param tuneK concentration tuning parameter, standard deviation of proposal distribution
#' @param burn_in number of draws to use as burn-in
#' @param m number of draws to keep from posterior distribution
#' @return list of
#' \itemize{
#' \item \code{Shat} Mode of the posterior distribution for the central orientation S
#' \item \code{kappa} Mean of the posterior distribution for the concentration kappa
#' }
#' @seealso \code{\link{mean.SO3}}, \code{\link{median.SO3}}
#' @details bingham2009b bingham2010
#' @export
#' @examples
#' Rs <- ruars(20, rvmises, kappa = 10)
#'
#' Shat <- mean(Rs) #Estimate the central orientation using the projected mean
#' rotdist.sum(Rs, Shat, p = 2) #The projected mean minimizes the sum of squared Euclidean
#' rot.dist(Shat) #distances, compute the minimized sum and estimator bias
#'
#' #Estimate the central orientation using the posterior mode (not run due to time constraints)
#' #Compare it to the projected mean in terms of the squared Euclidean distance and bias
#' \donttest{
#' ests <- bayes.mean(Rs, type = "Mises", S0 = mean(Rs), kappa0 = 10, tuneS = 5000,
#' tuneK = 1, burn_in = 1000, m = 5000)
#'
#' Shat2 <- ests$Shat #The posterior mode is the 'Shat' object
#' rotdist.sum(Rs, Shat2, p = 2) #Compute sum of squared Euclidean distances
#' rot.dist(Shat2) #Bayes estimator bias}
bayes.mean<-function(x,type,S0,kappa0,tuneS,tuneK,burn_in,m=5000){
UseMethod("bayes.mean")
}
#' @rdname bayes.mean
#' @export
bayes.mean.SO3<-function(x,type,S0,kappa0,tuneS,tuneK,burn_in,m=5000){
if(type %in% c("Cayley","cayley")){
#lpangle <- lpcayley
lpangle <- 1
}else if(type %in% c("Fisher","fisher")){
#lpangle <- lpfisher
lpangle <- 2
}else if(type %in% c("Mises","mises")){
#lpangle <- lpvmises
lpangle <- 3
}else{
stop("Invalid choise of type: please choose Cayley, Fisher or Mises.")
}
listRes<-both_MCMC_CPP(x,S0, kappa0,tuneS,tuneK,burn_in,m, lpangle)
postRes<-list(Shat=mean.SO3(listRes$S),kappa=mean(listRes$kappa))
return(postRes)
}
#' @rdname bayes.mean
#' @export
bayes.mean.Q4<-function(x,type,S0,kappa0,tuneS,tuneK,burn_in,m=5000){
Rs<-as.SO3(x)
S0<-as.SO3(matrix(S0,3,3))
SO3Res<-MCMCSO3(Rs,type,S0,kappa0,tuneS,tuneK,burn_in,m)
Q4Res<-list(Qhat=as.Q4(mean(SO3Res$S)),kappa=mean(SO3Res$kappa))
return(Q4Res)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.