knitr::opts_chunk$set(echo = TRUE, eval=FALSE) library(assertthat) library(numDeriv)
x <- seq(-10, 10, by=01) # log-concave example h <- function(x) { return(2*x - 10*log(1 + exp(x)) - 0.5*x^2 + 50) } # log-concave example h <- function(x) { return(log(dnorm(x))) } # NOT log-concave h_not <- function(x) { return(exp(x)) }
get_z <- function(j, x, h) { h_xj <- h(x[j]) h_prime_xj <- grad(h,x[j]) h_xjnext <- h(x[j+1]) h_prime_xjnext <- grad(h, x[j+1]) z_numerator <- h_xjnext - h_xj - ( x[j+1] * h_prime_xjnext ) + (x[j] * h_prime_xj) z_denominator <- h_prime_xj - h_prime_xjnext z_result <- z_numerator / z_denominator ## ASSERTION ## .. check that obtained z intersection is within x assert_that(z_result > x[1] && z_result < x[length(x)], msg = "z intersection is outside of X") return(z_result) }
get_z_all <- function(x, h, D) { return(unlist(lapply(1:(length(x) - 1), get_z, x, h))) }
input indices j+1, ..., l-1 avoid unbounded D
check_log_concave <- function(j,z,h) { assert_that(is.finite(z[j]), msg = "z_j not finite") assert_that(is.finite(z[j+1]), msg ="z_{j+1} not finite") cat('full z is: ', z, '\n') # integrate u(x) over x_{index} and x_{index + 1} get_u_integral <- function() { uj <- function(t) { return(u[[t]]$intercept + u[[t]]$slope * t) } fun_u <- function(t) { exp(uj(t)) } return(integrate(fun_u, z[j], z[j+1])$value) } u_integral <- get_u_integral() ## store result, I_u(x) # integrate h(x) over x_{index} and x_{index + 1} get_h_integral <- function() { fun_h <- function(t) { exp(h(t)) } return(integrate(fun_h, z[j], z[j+1])$value) } h_integral <- get_h_integral() ## store result, I_h(x) # OUTPUT I_u(x) and I_h(x) cat("u_integral: ", u_integral, "\nh_integral: ", h_integral, "\n") # CHECK if I_u(x) > I_h(x) and OUTPUT cat(ifelse(u_integral > h_integral, "upper bound segment > h(x)\n" , "ERROR upper bound segment < h(x)\n")) cat("\n") # integrate l(x) over x_{index} and x_{index + 1} get_l_integral <- function() { lj <- function(t) { return(l[[t]]$intercept + l[[t]]$slope*t) } fun_l <- function(t) { exp(lj(t)) } return(integrate(fun_l, z[j], z[j+1])$value) } l_integral <- get_l_integral() ## store result, I_l(x) # OUTPUT I_l(x) and I_h(x) cat("l_integral: ", l_integral, "\nh_integral: ", h_integral, "\n") # CHECK if I_l(x) < I_h(x) and OUTPUT cat(ifelse(l_integral < h_integral, "lower bound segment < h(x)\n" , "ERROR lower bound segment > h(x)\n")) assert_that(l_integral < h_integral, msg = "lower bound segment is not less than INPUT density") assert_that(u_integral > h_integral , msg = "upper bound segment is not greater than INPUT density") }
check_log_concave <- function(j,x,h) { # integrate u(x) over x_{index} and x_{index + 1} get_u_integral <- function() { uj <- function(t) { tmp_u <- get_u_segment(j,x,h) return(tmp_u$intercept + tmp_u$slope*t) } fun_u <- function(t) { exp(uj(t)) } return(integrate(fun_u, x[j], x[j+1])$value) } u_integral <- get_u_integral() ## store result, I_u(x) # integrate h(x) over x_{index} and x_{index + 1} get_h_integral <- function() { fun_h <- function(t) { exp(h(t)) } return(integrate(fun_h, x[j], x[j+1])$value) } h_integral <- get_h_integral() ## store result, I_h(x) # OUTPUT I_u(x) and I_h(x) cat("u_integral: ", u_integral, "\nh_integral: ", h_integral, "\n") # CHECK if I_u(x) > I_h(x) and OUTPUT cat(ifelse(u_integral > h_integral, "upper bound segment > h(x)\n" , "ERROR upper bound segment < h(x)\n")) cat("\n") # ---------------------- # # ---------------------- # # ---------------------- # # integrate l(x) over x_{index} and x_{index + 1} get_l_integral <- function() { lj <- function(t) { tmp_l <- get_l_segment(j,x,h) return(tmp_l$intercept + tmp_l$slope*t) } fun_l <- function(t) { exp(lj(t)) } return(integrate(fun_l, x[j], x[j+1])$value) } l_integral <- get_l_integral() ## store result, I_l(x) # OUTPUT I_l(x) and I_h(x) cat("l_integral: ", l_integral, "\nh_integral: ", h_integral, "\n") # CHECK if I_l(x) < I_h(x) and OUTPUT cat(ifelse(l_integral < h_integral, "lower bound segment < h(x)\n" , "ERROR lower bound segment > h(x)\n")) # ASSERT proper upper and lower bounds for u(x) and l(x), respectively # ... otherwise, error message and exit assert_that(l_integral < h_integral, msg = "lower bound segment is not less than INPUT density") assert_that(u_integral > h_integral , msg = "upper bound segment is not greater than INPUT density") }
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.