R/exact_algorithm_mixG.R

In rchan26/mixGaussTempering: What the Package Does (One Line, Title Case)

##########################
# these following functions make up the Exact Algorithm to simulate Langevin diffusion for tempered mixture Gaussian
# the Exact Algorithm requires the first and second derivatives
# the first part of this script deals with functions to obtain these derivatives
# we obtain these easily by noticing that differentiating a mixture Gaussian is easy as it is a linear combination of Gaussians
# and noticing that we can differentiate a tempered mixture Gaussian easily using the chain rule
##########################

######################################## derivatives of tempered mixture Gaussian to get 'phi function' ########################################

#' Check vectors have equal dimension
#'
#' Function to check that a vectors in a list are of the same length
#'
#' @param list_of_vectors list with elements that are vectors
#'
#' @return boolean value: TRUE if all vectors in list are the same length, FALSE otherwise
#'
#' @examples
#' check_equal_dim(list(c(0,1,2), c(0,1))) # FALSE
#' check_equal_dim(list(c(0,1,2), c(0,1,2))) # TRUE
#'
#' @export
check_equal_dim <- function(list_of_vectors) {
  # function to check the dimensions of weights, means and sds are equal
  check_eq <- sapply(list_of_vectors[-1], function(x) length(x)==length(list_of_vectors[[1]]))
  if (prod(check_eq)==0) {
    # one of the dimensions don't match
    return(FALSE)
  } else {
    return(TRUE)
  }
}

#' Simulate exactly a mixture Gaussian distribution
#'
#' Simulate exactly a mixture Gaussian distribution using component sampling
#'
#' @param n number of samples required
#' @param weights vector: weights of mixture Gaussian
#' @param means vector: means of mixture Gassuan
#' @param sds vector: st.devs of mixture Gaussian
#'
#' @return samples: samples of a mixture Gaussian distribution
#'
#' @examples
#' rmix_norm(n = 10000, weights = c(0.4, 0.6), means = c(-3, 9), sds = c(2, 0.4))
#'
#' @export
rmix_norm <- function(n, weights, means, sds) {
  # function to sample directly from mixture Gaussian
  if (!(sum(weights)==1)) stop('Error: Weights do not add up to 1.')
  if (!check_equal_dim(list(weights, means, sds))) stop('Dimensions of input vectors do not match')

  w <- sample(1:length(weights), prob=weights, size=n, replace=T)
  samples <- rnorm(n, mean=means[w], sd=sds[w])
  return(samples)
}

#' Function of a un-normalised tempered mixture Gaussian
#'
#' Returns the value of a un-normalised tempered mixture Gaussian at given point x
#'
#' @param x real value
#' @param w vector: weights of mixture Gaussian
#' @param m vector: means of mixture Gassuan
#' @param s vector: st.devs of mixture Gaussian
#' @param b real value between 0 and 1
#'
#' @return value: value of a un-normalised tempered mixture Gaussian at given point x
#'
#' @examples
#' curve(dnorm_mix_tempered_unnorm(x, weights = c(0.4, 0.6), means = c(-3, 9), sds = c(1,1), beta = 1/4), -15, 20)
#'
#' @export
dnorm_mix_tempered_unnorm <- function(x, w, m, s, b) {
  value <- 0
  for (i in 1:length(w)) {value <- value + w[i]*dnorm(x, mean=m[i], sd=s[i])}
  value <- value^(b)
  return(value)
}

#' Density of tempered mixture Gaussian
#'
#' Returns the value of a normalised tempered mixture Gaussian at given point x
#'
#' @param x real value
#' @param weights vector: weights of mixture Gaussian
#' @param means vector: means of mixture Gassuan
#' @param sds vector: st.devs of mixture Gaussian
#' @param beta real value between 0 and 1
#'
#' @return value: value of a normalised tempered mixture Gaussian at given point x
#'
#' @examples
#' curve(dnorm_mix_tempered(x, weights = c(0.4, 0.6), means = c(-3, 9), sds = c(1,1), beta = 1/4), -15, 20)
#'
#' @export
dnorm_mix_tempered <- function(x, weights, means, sds, beta) {
  # check dimensions of input values
  if (!check_equal_dim(list(weights, means, sds))) stop('Dimensions of input vectors do not match')

  norm_constant <- integrate(function(x) dnorm_mix_tempered_unnorm(x=x, w=weights, m=means, s=sds, b=beta),
                             lower=-Inf, upper=Inf)$value
  value <- dnorm_mix_tempered_unnorm(x, w=weights, m=means, s=sds, b=beta) / norm_constant
  return(value)
}

#' Density of a mixture Gaussian
#'
#' Returns the value of a mixture Gaussian at given point x
#'
#' @param x real value
#' @param weights vector: weights of mixture Gaussian
#' @param means vector: means of mixture Gassuan
#' @param sds vector: st.devs of mixture Gaussian
#'
#' @return value: value of a mixture Gaussian at given point x
#'
#' @examples
#' curve(dnorm_mix(x, weights = c(0.4, 0.6), means = c(-3, 9), sds = c(1,1)), -10, 15)
#'
#' @export
dnorm_mix <- function(x, weights, means, sds) {
  # check dimensions of input values
  if (!check_equal_dim(list(weights, means, sds))) stop('Dimensions of input vectors do not match')

  value <- 0
  for (i in 1:length(weights)) {value <- value + weights[i]*dnorm(x, mean=means[i], sd=sds[i])}
  return(value)
}

#' First derivative of a normal distribution
#'
#' Returns the value of the first derivative of a normal distribution at a given point x
#'
#' @param x real value
#' @param mean real value (default 0)
#' @param sd real value (default 1)
#'
#' @return value: value of a first derivative of a normal distribution at a given point x
#'
#' @examples
#' curve(dnorm_deriv1(x), -5, 5)
#'
#' @export
dnorm_deriv1 <- function(x, mean=0, sd=1) {
  return(((mean-x)/(sd^2))*dnorm(x, mean, sd))
}

#' Second derivative of a normal distribution
#'
#' Returns the value of the second derivative of a normal distribution at a given point x
#'
#' @param x real value
#' @param mean real value (default 0)
#' @param sd real value (default 1)
#'
#' @return value: value of a second derivative of a normal distribution at a given point x
#'
#' @examples
#' curve(dnorm_deriv2(x), -5, 5)
#'
#' @export
dnorm_deriv2 <- function(x, mean=0, sd=1) {
  return(((((x-mean)^2)/(sd^4))*dnorm(x, mean, sd)) - (dnorm(x, mean, sd)/(sd^2)))
}

#' Function of a the first derivative of a un-normalised tempered mixture Gaussian
#'
#' Returns the value of the first derivative of a un-normalised tempered mixture Gaussian at given point x
#'
#' @param x real value
#' @param weights vector: weights of mixture Gaussian
#' @param means vector: means of mixture Gassuan
#' @param sds vector: st.devs of mixture Gaussian
#'
#' @return value: value of the first derivative of a un-normalised tempered mixture Gaussian at given point x
#'
#' @examples
#' curve(dnorm_mix_deriv1(x, weights = c(0.4, 0.6), means = c(-3, 9), sds = c(1,1)), -15, 20)
#'
#' @export
dnorm_mix_deriv1 <- function(x, weights, means, sds) {
  # check dimensions of input values
  if (!check_equal_dim(list(weights, means, sds))) stop('Dimensions of input vectors do not match')

  value <- 0
  for (i in 1:length(weights)) {value <- value + weights[i]*dnorm_deriv1(x, mean=means[i], sd=sds[i])}
  return(value)
}

#' Function of a the second derivative of a un-normalised tempered mixture Gaussian
#'
#' Returns the value of the second derivative of a un-normalised tempered mixture Gaussian at given point x
#'
#' @param x real value
#' @param weights vector: weights of mixture Gaussian
#' @param means vector: means of mixture Gassuan
#' @param sds vector: st.devs of mixture Gaussian
#'
#' @return value: value of the second derivative of a un-normalised tempered mixture Gaussian at given point x
#'
#' @examples
#' curve(dnorm_mix_deriv2(x, weights = c(0.4, 0.6), means = c(-3, 9), sds = c(1,1)), -15, 20)
#'
#' @export
dnorm_mix_deriv2 <- function(x, weights, means, sds) {
  # check dimensions of input values
  if (!check_equal_dim(list(weights, means, sds))) stop('Dimensions of input vectors do not match')

  value <- 0
  for (i in 1:length(weights)) {value <- value + weights[i]*dnorm_deriv2(x, mean=means[i], sd=sds[i])}
  return(value)
}

#' Function of a the first derivative of a normalised tempered mixture Gaussian
#'
#' Returns the value of the first derivative of a normalised tempered mixture Gaussian at given point x
#'
#' @param x real value
#' @param weights vector: weights of mixture Gaussian
#' @param means vector: means of mixture Gassuan
#' @param sds vector: st.devs of mixture Gaussian
#' @param beta real value between 0 and 1
#'
#' @return value: value of the first derivative of a normalised tempered mixture Gaussian at given point x
#'
#' @examples
#' curve(dnorm_mix_tempered_deriv1(x, weights = c(0.4, 0.6), means = c(-3, 9), sds = c(1,1), beta = 1/4), -15, 20)
#'
#' @export
dnorm_mix_tempered_deriv1 <- function(x, weights, means, sds, beta) {
  # check dimensions of input values
  if (!check_equal_dim(list(weights, means, sds))) stop('Dimensions of input vectors do not match')

  q <- dnorm_mix(x, weights, means, sds)
  q_prime <- dnorm_mix_deriv1(x, weights, means, sds)
  norm_constant <- integrate(function(x) dnorm_mix_tempered_unnorm(x=x, w=weights, m=means, s=sds, b=beta),
                             lower=-Inf, upper=Inf)$value

  return((beta*(q^(beta-1))*q_prime)/norm_constant)
}

#' Function of a the second derivative of a normalised tempered mixture Gaussian
#'
#' Returns the value of the second derivative of a normalised tempered mixture Gaussian at given point x
#'
#' @param x real value
#' @param weights vector: weights of mixture Gaussian
#' @param means vector: means of mixture Gassuan
#' @param sds vector: st.devs of mixture Gaussian
#' @param beta real value between 0 and 1
#'
#' @return value: value of the second derivative of a normalised tempered mixture Gaussian at given point x
#'
#' @examples
#' curve(dnorm_mix_tempered_deriv2(x, weights = c(0.4, 0.6), means = c(-3, 9), sds = c(1,1), beta = 1/4), -15, 20)
#'
#' @export
dnorm_mix_tempered_deriv2 <- function(x, weights, means, sds, beta) {
  # check dimensions of input values
  if (!check_equal_dim(list(weights, means, sds))) stop('Dimensions of input vectors do not match')

  q <- dnorm_mix(x, weights, means, sds)
  q_prime <- dnorm_mix_deriv1(x, weights, means, sds)
  q_double_prime <- dnorm_mix_deriv2(x, weights, means, sds)

  term1 <- beta*(beta-1)*(q^(beta-2))*(q_prime^2)
  term2 <- beta*(q^(beta-1))*(q_double_prime)
  norm_constant <- integrate(function(x) dnorm_mix_tempered_unnorm(x=x, w=weights, m=means, s=sds, b=beta),
                             lower=-Inf, upper=Inf)$value

  return((term1 + term2)/norm_constant)
}



######################################## 'phi function' for the Exact Algorithm for tempered mixture Gaussians ########################################

#' phi-function for a tempered mixture Gaussian
#'
#' phi-function for the Exact Algorithm for a tempered mixture Gaussian
#'
#' @param x real value
#' @param weights vector: weights of mixture Gaussian
#' @param means vector: means of mixture Gassuan
#' @param sds vector: st.devs of mixture Gaussian
#' @param beta real value between 0 and 1
#'
#' @return real value
#'
#' @examples
#' curve(phi_function_mixG(x, weights = c(0.4, 0.6), means = c(-3, 9), sds = c(1,1), 1), from = -15, to = 20)
#'
#' @export
phi_function_mixG <- function(x, weights, means, sds, beta) {
  # check dimensions of input values
  if (!check_equal_dim(list(weights, means, sds))) stop('Dimensions of input vectors do not match')

  fc <- dnorm_mix_tempered(x, weights, means, sds, beta)
  fc_prime <- dnorm_mix_tempered_deriv1(x, weights, means, sds, beta)
  fc_double_prime <- dnorm_mix_tempered_deriv2(x, weights, means, sds, beta)

  term1 <- (fc_prime/fc)^2
  term2 <- ((fc)*(fc_double_prime)-(fc_prime)^2) / (fc^2)

  return((term1+term2)/2)
}

#' Obtain bounds for phi function
#'
#' Finds the lower and upper bounds of the phi function between an interval
#'
#' @param weights vector: weights of mixture Gaussian
#' @param means vector: means of mixture Gassuan
#' @param sds vector: st.devs of mixture Gaussian
#' @param beta real value between 0 and 1
#' @param lower lower end of interval
#' @param upper upper end of interval
#'
#' @return vector of lower and upper bound of phi-function between [lower, upper]
#'
#' @examples
#' curve(phi_function_mixG(x, weights = c(0.4, 0.6), means = c(-3, 9), sds = c(1,1), 1), from = -15, to = 20)
#' abline(v=c(-10, -2))
#' bounds <- bounds_phi_function_mixG(weights = c(0.4, 0.6), means = c(-3, 9), sds = c(1,1), beta = 1, lower = -10, upper = -2)
#' abline(h=c(bounds$fun_ubound, bounds$fun_lbound))
#'
#' @export
bounds_phi_function_mixG <- function(weights, means, sds, beta, lower = min(interval), upper = max(interval)) {
  # finding maxima and minimma in the interval
  maxi <- optimise(function(x) phi_function_mixG(x, weights, means, sds, beta),
                   interval = c(lower, upper), maximum = TRUE)$objective
  mini <- optimise(function(x) phi_function_mixG(x, weights, means, sds, beta),
                   interval = c(lower, upper), maximum = FALSE)$objective

  # checking end points
  lower_end <- phi_function_mixG(lower, weights, means, sds, beta)
  upper_end <- phi_function_mixG(upper, weights, means, sds, beta)
  # if the max and min of the end points are larger or smaller then they are replaced

  up_bound <- ifelse(test = max(lower_end, upper_end) > maxi, yes = max(lower_end, upper_end), no = maxi)
  low_bound <- ifelse(test = min(lower_end, upper_end) < mini, yes = min(lower_end, upper_end), no = mini)

  return(list('fun_ubound' = up_bound, 'fun_lbound' = low_bound))
}

#' Lower bound finder for tempered mixture Gaussian
#'
#' Finds the lower bound of a tempered mixture Gaussian
#'
#' @param weights vector: weights of mixture Gaussian
#' @param means vector: means of mixture Gassuan
#' @param sds vector: st.devs of mixture Gaussian
#' @param beta real value between 0 and 1
#'
#' @return value: lower bound of a tempered mixture Gaussian
#'
#' @examples
#' curve(phi_function_mixG(x, weights = c(0.4, 0.6), means = c(-3, 9), sds = c(1,1), 1), from = -15, to = 20)
#' abline(h = lower_bound_phi_function_mixG(weights = c(0.4, 0.6), means = c(-3, 9), sds = c(1,1), 1))
#'
#' @export
lower_bound_phi_function_mixG <- function(weights, means, sds, beta) {
  # using 'optimise' to find the minima around the means
  # local minima of phi occur around the means - so the global maxima is typically the smallest of these
  mins <- sapply(1:length(means), function(mode) optimise(function(x) phi_function_mixG(x, weights, means, sds, beta),
                                                          interval = means[mode]+c(-2*sds[mode], 2*sds[mode]), maximum = FALSE)$objective)

  return(min(mins))
}

######################################## exact algorithm to simulate langevin diffusion ########################################

#' Exact Algorithm for langevin diffusion with pi as a tempered mixture Gaussian
#'
#' Simulate langevin diffusion using the Exact Algorithm with pi = exp(-(beta*x^4)/2)
#'
#' @param x0 start value
#' @param y end value
#' @param s start time
#' @param t end time
#' @param K lower bound of the phi function
#' @param weights vector: weights of mixture Gaussian
#' @param means vector: means of mixture Gassuan
#' @param sds vector: st.devs of mixture Gaussian
#' @param beta real value between 0 and 1
#' @param accept_prob boolean value: defaults to F, determined whether or not to return acceptance prob (T) or a skeleton of the diffusion (F)
#'
#' @return accept_prob real value if accept_prob = TRUE
#' @return layered_bb matrix holding skeleton of the diffusion
#'
#' @examples
#' lower_bound <- lower_bound_phi_function_mixG(weights = c(0.4, 0.6), means = c(-3, 9), sds = c(1,1), 1)
#' test <- simulate_langevin_diffusion_mixG(x0 = 0.6, y = 1.2, s = 0, t = 2, K = lower_bound, accept_prob = FALSE,
#'                                          weights = c(0.5, 0.5), means = c(-2, 2), sds = c(1,1), beta = 1)
#' plot(x = test['time',], y = test['X',], type = 'l')
#'
#' @export
simulate_langevin_diffusion_mixG <- function (x0, y, s, t, K, weights, means, sds, beta, accept_prob = FALSE) {
  # simulates langevin diffusion for f_{c} = mixture Gaussian
  # if accept_prob = TRUE, function returns the acceptance probability only

  repeat {
    # first simulate layer information to bound sample path
    a_seq <- seq(from = 0, to = ceiling(t-s), by = 0.1)
    bes_layer <- bessel_layer_simulation(x0, y, s, t, a_seq)
    lbound <- min(x0,y) - bes_layer$a[bes_layer$l]
    ubound <- max(x0,y) + bes_layer$a[bes_layer$l]

    # calculate upper and lower bounds of phi given the simulated sample layer information
    bounds <- bounds_phi_function_mixG(weights=weights, means=means, sds=sds, beta=beta, lower=lbound, upper=ubound)
    LX <- bounds$fun_lbound
    UX <- bounds$fun_ubound

    if (!(accept_prob)) {
      if (runif(1,0,1) < (1 - exp(-(LX-K)*(t-s)))) {
        next # go back to step 1: recall the function to start again
      }
    }

    # simulating skeletal points via Poisson thinning: simulating number of points, then simulating Poisson times
    kap <- rpois(1, (UX - LX)*(t-s))
    pois_times <- sort(runif(kap, s, t))

    # simulating sample path at Poisson times
    layered_bb <- layered_brownian_bridge(x0, y, s, t, bes_layer$a, bes_layer$l, pois_times)

    # layered_bb includes points for times s,t,tau, so only keeping those from pois_times
    pois_points <- layered_bb[1, (layered_bb[2,] %in% pois_times)]

    # calculating acceptance probability
    acc_prob <- ifelse(test = length(pois_times)==0, yes = 1,
                       no = prod((UX - phi_function_mixG(pois_points, weights, means, sds, beta))/(UX - LX)))

    if (accept_prob) {
      return(acc_prob * (exp(-(LX-K)*(t-s))))
    }

    if (runif(1,0,1) < acc_prob) {
      return(layered_bb) # accept sample path and return the entire path
    }
  }
}