R/polar_log_lik.R

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

#This file contains functions for calculating likelihoods


#Log likelihood of observing rhat and thetahat given the truth r and theta
#Assumes z1hat ~ N(z1, 1); z2hat ~ N(z2,1) independent (for now)
#rhat = sqrt(z1hat^2 + z2hat^2); thetahat = atan2(z2hat, z1hat)
# r= sqrt(z1^2 + z2^2); theta = atan2(z2,z1)
#'@title Calculate log p(rhat, thetahat | r, theta)
#'@description Assumes z1hat ~ N(z1, 1); z2hat ~ N(z2,1) independent
#'rhat = sqrt(z1hat^2 + z2hat^2); thetahat = atan2(z2hat, z1hat)
#'@param r r
#'@param theta theta
#'@param rhat rhat
#'@param thetahat thethat
#'@export
polar_log_lik_rt <- function(r, theta, rhat, thetahat){
  ll <- log(rhat)-log(2*pi) -0.5*(rhat^2 + r^2 - 2*rhat*r*cos(thetahat-theta))
  return(ll)
}


#Log likelihood of observing rhat and thetahat after the data have been oriented so that
#theta is between -pi/2 and pi/2 (i.e. so that beta_hat_1 is positive)
#given the truth r and theta
#Assumes z1hat ~ N(z1, 1); z2hat ~ N(z2,1) independent (for now)
#rhat = sqrt(z1hat^2 + z2hat^2); thetahat = atan2(z2hat, z1hat)
# r= sqrt(z1^2 + z2^2); theta = atan2(z2,z1)

#'@title Calculate log p(rhat, thetahat_o | r, theta)
#'@description Assumes z1hat ~ N(z1, 1); z2hat ~ N(z2,1) independent
#'z1hat_o = abs(z1hat); z2hat_o = sign(z1hat)*z2hat
#'This forces thetahat to be between -pi/2 and pi/2
#'rhat = sqrt(z1hat^2 + z2hat^2); thetahat = atan2(z2hat, z1hat)
#'@param r r
#'@param theta theta
#'@param rhat rhat
#'@param thetahat thethat_o
#'@export
polar_log_lik_rt2 <- function(r, theta, rhat, thetahat){
  stopifnot(-pi/2 <= thetahat & thetahat <= pi/2)
  ll <- log(rhat)-log(2*pi) -0.5*(rhat^2 + r^2) +
        log(exp(-rhat*r*cos(thetahat-theta)) + exp(rhat*r*cos(thetahat-theta)))
  return(ll)
}

#marginal likelihood of observing rhat given truth r (does not depend on theta)
#when r is away from 0, this can be approximated by N(r, 1)
#'@title Calculate the marginal log likelihood log p(rhat | r) (which does not depend on theta)
#'@param r r
#'@param rhat rhat
#'@export
polar_lik_marg_r <- function(r, rhat){
  lik <- rhat*exp(-0.5*(rhat^2 + r^2))*besselI(rhat*r, 0)
  return(lik)
}

#'@title Calculate b*p(rhat | r uniform between 0 and b)= int_0^b p(rhat | r)dr
#'@param b b
#'@param rhat rhat
#'@return Returns b*p(rhat | r uniform between 0 and b)
#'@export
polar_log_lik_marg_r_uni <- function(b, rhat){
  f <- function(x){
    rhat*exp(-0.5*(rhat^2 + x^2))*besselI(rhat*x, 0)
  }
  val <- integrate(f, lower=0, upper=b)$value
  return(log(val))
}

#'@title Calculate b*p(rhat, thetahat | theta (fixed), r uniform between 0 and b)
#'@param b b
#'@param theta theta
#'@param rhat rhat
#'@param thetahat thetahat
#'@return Returns b*p(rhat, thetahat | theta (fixed), r uniform between 0 and b)
#'@export
polar_log_lik_rt_runi <- function(b, theta, rhat, thetahat){
  if(r==0 & theta!=0) return(0)
  f <- function(x){
    exp(polar_log_lik_rt(x, theta, rhat, thetahat))
  }
  val <- integrate(f, lower=0,upper=b)$value
  return(log(val))
}

#'@title Calculate b*p(rhat, thetahat_o | theta (fixed), r uniform between 0 and b)
#'@param b b
#'@param theta theta
#'@param rhat rhat
#'@param thetahat thetahat_o
#'@return Returns b*p(rhat, thetahat_o | theta (fixed), r uniform between 0 and b)
#'@export
polar_log_lik_rt2_runi <- function(r, theta, rhat, thetahat){
  if(r==0 & theta!=0) return(0)
  f <- function(x){
    exp(polar_log_lik_rt2(x, theta, rhat, thetahat))
  }
  val <- integrate(f, lower=0,upper=r)$value
  return(log(val))
}



#likelihood of observing thetahat given truth r and theta
#when r is away from 0, this can be approximated by N(theta, 1/r^2)
#'@title Calculate the (approximat) marginal log likelihood log p(thetahat | r, theta) (does depend on r!)
#'@param r r
#'@param theta theta
#'@param thetahat thetahat
#'@export
polar_lik_marg_t <- function(r, theta, thetahat){
  if(r==0 & theta!=0) return(0)
  f <- function(x){
    exp(polar_log_lik_rt(r, theta, x, thetahat))
  }
  val <- integrate(f, lower=0,upper=100)$value
  return(log(val))
}