R/polar_two_stage.R

In jean997/bvpolar: Estimate Distribution of Angles from Bivariate Data

#'Get posterior distribution of theta and empirical distribution or r
#'@param rhat Vector of rhat
#'@param thetahat Vector of thetahat
#'@param orient indicates that data have been oriented so that thetahat is between -pi/2 and pi/2
#'@param gr Object of class unimix giving elements of r mixture
#'@param rprior Prior for elements of gr
#'@param tvals,tprior grid and prior values of theta
#'@export
polar_two_stage_uni <- function(rhat, thetahat, orient=FALSE,
                                gr=NULL, tvals=NULL, rprior=NULL, tprior=NULL){

  if(orient & (is.null(tvalls) | is.null(tprior))){
    tvals <- seq(-1/2, 1/2, lenght.out=21)*pi
    tprior <- rep(c(10, 1, 10, 1, 10), c(1, 9, 1, 9, 1))
  }else if(is.null(tvalls) | is.null(tprior)){
    tvals <- seq(-1, 1, lenght.out=21)*pi
    tprior <- rep(c(10, 1, 10, 1, 10), c(1, 4, 1, 4, 1, 4, 1, 4, 1))
  }
  #Stage 1. Fit r
  cat("Fitting r\n")
  t0 <- system.time(fit_r <- polar_est_r_mixture_uni(rhat, gr, rprior))

  cat("\n", t0, "\n")
  #Stage 2: Fit theta
  cat("Fitting theta\n")
  t2 <- system.time(fit_t <-polar_est_t(rhat, thetahat, fit_r, tvals, tprior, orient))

  cat("\n", t2, "\n")
  ret <- list("fit_r"=fit_r, "fit_t"=fit_t)
  return(ret)
}