#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))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.