examples for Klemens paper in R

klemens.R

################################################################################
# Definition 4
################################################################################

norm_for_Lp <- function(x,p){
  dnorm(x = x,mean = p,sd = 1,log = FALSE)
}

Lp <- function(p){
  integrate(f = norm_for_Lp,lower = -Inf,upper = Inf,p=p)$value
}

Ldp <- function(d,p){
  dnorm(x = d,mean = p,sd = 1,log = FALSE) / Lp(p)
}

p <- 3
D <- seq(-10,10,by=0.01)
xx <- Ldp(D,p)

plot(D,xx,type="l",xlab="Data Space",ylab="L(d|p)")


################################################################################
# Definition 5
################################################################################

norm_for_Ld <- function(x,d){
  dnorm(x = d,mean = x,sd = 1,log = FALSE)
}

Ld <- function(d){
  integrate(f = norm_for_Ld,lower = -Inf,upper = Inf,d=d)$value
}

Lpd <- function(p,d){
  dnorm(x = d,mean = p,sd = 1,log = FALSE) / Ld(d)
}

d <- 10 # our datum

Lpd(d = d,p = 10)

p <- seq(0,20,by=0.001)

xx <- Lpd(p = p,d = d)
plot(p,xx,type="l",xlab = "Parameter Space",ylab = "L(p|d)")


################################################################################
# 4.5 Data Constraint
################################################################################

# example 1: truncation
L_prime_normalising_constant <- function(ub,p){
  integrate(f = dnorm,lower = -Inf,upper = ub,mean = p,sd = 1,log = FALSE)$value
}

L_prime <- function(d,p,ub){
  if(d < ub){
    dnorm(x = d,mean = p,sd = 1,log = FALSE) / L_prime_normalising_constant(ub = ub,p = p)
  } else {
    0
  }
}

x <- seq(-6,6,by=0.01)

# constrain our Gaussian at 2.135 (totally arbitrary) with a mean of 0.855
ub <- 2.135
p <- 1.855
l <- sapply(x,L_prime,p=p,ub=ub)

plot(x,l,type="l",xlab="Data Space",ylab="L'(d,p)")

# check it integrates to 1 over the truncated region
suppressWarnings(integrate(f=L_prime,lower = -Inf,upper = ub,p=0.5,ub=ub))
# voila, a proper distribution.

# example 2: censoring
L_prime_2 <- function(d,mu,sd,lb){
  if(d >= lb){
    dnorm(x = d,mean = mu,sd = sd) / integrate(f = dnorm,lower = lb,upper = Inf,mean = mu,sd = sd)$value
  } else {
    0
  }
}

L_prime_censoring <- function(d,mu,sd){
  (pnorm(q = 0,mean = mu,sd = sd)*1) + ((1 - pnorm(q = 0,mean = mu,sd = sd))*L_prime_2(d = d,mu=mu,sd=sd,lb=0))
}

data <- runif(n = 1e3,min = -1,max = 5)
data[data < 0] <- 0 # censor

mu <- 0.85
sd <- 2
like <- sapply(data,L_prime_censoring,mu=mu,sd=sd)

plot(data,like,pch=16,
     col=adjustcolor("firebrick3",alpha.f = 0.5))


################################################################################
# 4.6 Differentiable Transformation
################################################################################

library(numDeriv)
library(truncnorm)

# we want to lower bound our X at -8.25, because reasons
a <- -8.25

# our X is truncated normal with some parameters
mean <- -3.96
sd <- 2.3

# but we might be throwing this X into an optimization or MCMC routine that wants
# to sample from the whole real line. So we need to transform X.
# here's how:

# plug in x, get y
f <- function(x,a){
  log(x - a)
}

# plug in y, get x
f_inv <- function(y,a){
  exp(y) + a
}

# we can calculate the determinant of the jacobian by hand (its exp(y)) 
# but let's assume we can't (usually), we can get it numerically:
det_jacob <- function(y,a){
  abs(jacobian(func = f_inv,x = y,a=a))[1,1]
}

# density of x
# it was truncated at -8.25
p_x <- function(x,a,mean,sd){
  dtruncnorm(x = x,a = a,b = Inf,mean = mean,sd = sd)
}

# density of y on the unconstrained space
p_y <- function(y,a,mean,sd){
  p_x(f_inv(y,a),a,mean,sd) * det_jacob(y,a)
}

x <- seq(ceiling(a),5,by=0.01)
y <- f(x,a)

# evaluate densities
# computing p_y takes awhile, its doing lots of computation to save
# you the need to do math (always good)
x_like <- p_x(x = x,a = a,mean = mean,sd = sd)
y_like <- p_y(y = y,a = a,mean = mean,sd = sd)

par(mfrow=c(1,2))

plot(x,x_like,type="l",col="steelblue",
     main="Constrained X",lwd=2.5,
     xlab="X",ylab="P_{X}")

plot(y,y_like,type="l",col="firebrick3",
     main="Unconstrained Y",lwd=2.5,
     xlab="Y",ylab="P_{Y}")

par(mfrow=c(1,1))