R/polar_est_mixture_uni.R

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

#'Estimate the distribution or r from rhat
#'@param rhat Vector of rhat
#'@param g object of class "unimix"
#'@param prior Prior weights for each element of g
#'@param ... Args to be passed to get_grid_r
#'@export
polar_est_r_mixture_uni <- function(rhat, g=NULL, prior=NULL, ...){
  if(is.null(g)){
    g <- get_grid_r(rhat, ...)
    prior <- g$prior
    g <- g$g
  }

  stopifnot(class(g)=="unimix")
  stopifnot(all(g$a==0))
  K <- length(g$b)
  p <- length(rhat)
  llmat <- matrix(nrow=p, ncol=K)

  for(i in 1:K){
    cat(i, " ")
    if(g$b[i] > 0){
      lls <- sapply(rhat, function(rh){
        polar_log_lik_marg_r_uni(g$b[i], rh)
      })
      llmat[,i] <- lls - log(g$b[i])
    }else{
      lls <- sapply(rhat, function(rh){
        log(polar_lik_marg_r(g$b[i], rh))
      })
      llmat[,i] <- lls
    }
  }
  matrix_llik = llmat - apply(llmat, 1, max)
  matrix_lik = exp(matrix_llik)
  w_res = ashr:::mixIP(matrix_lik =matrix_lik, prior=prior)
  mypi <-  w_res$pihat
  g$pi <- mypi
  ret <- list("fitted_g"=g, "llmat"=llmat, "prior"=prior)
  return(ret)
}

#'Simulate observations from a mixture of uniforms
#'@param n number of data points
#'@param g object of class "unimix"
#'@export
runimix <- function(n, g){
  stopifnot(class(g)=="unimix")
  k <- length(g$pi)
  z <- sample(1:k, size=n, replace=TRUE, prob = g$pi)
  x <- rep(NA, n)
  for(i in 1:k){
    x[z==i] <- runif(n=sum(z==i), min=g$a[i], max=g$b[i])
  }
  return(x)
}

#'Choose a reasonable set of uniform distributions to make up mixture for r
#'@param rhat Vector of rhat
#'@param mode mode of distribution (probably 0)
#'@param gridmult,grange,prior,nullweight See ashr
#'@export
get_grid_r <- function(rhat, mode=0, gridmult=sqrt(2),
                       grange=c(0, Inf), prior="nullbiased", nullweight=10){
  p <- length(rhat)
  data <- list("x"=rhat, "s"=rep(1, p), "lik"=ashr::lik_normal())
  mixsd = ashr:::autoselect.mixsd(data,gridmult,mode,grange,"+uniform")
  mixsd <- c(0, mixsd)
  k = length(mixsd)
  prior = ashr:::setprior(prior,k,nullweight,1)
  pi = ashr:::initpi(k,length(rhat),1)
  g = ashr::unimix(pi,rep(mode,k),mode+mixsd)
  return(list("g"=g, "prior"=prior))
}