library(cubature)
context("Test known integration results")
test_that("Test a product of Cosine functions", {
expected <- 0.708073
tol <- 1e-4
## Product of cosines
testFn0 <- function(x) prod(cos(x))
result <- cubature::adaptIntegrate(f = testFn0,
lowerLimit = rep(0, 2),
upperLimit = rep(1, 2),
tol = tol)
testthat::expect_equal(0, result$returnCode, info = "Integration unsuccessful!")
testthat::expect_equal(expected, result$integral, tolerance = tol, scale = abs(result$integral),
info = "Relative error not reached")
testthat::expect_equal(expected, result$integral, tolerance = tol, scale = 1,
info = "Absolute error not reached")
testFn0_v <- function(x) {
r <- apply(x, 2, function(z) prod(cos(z)))
matrix(r, ncol = ncol(x))
}
result <- cubature::adaptIntegrate(f = testFn0_v,
lowerLimit = rep(0, 2),
upperLimit = rep(1, 2),
tol = tol,
vectorInterface = TRUE)
testthat::expect_equal(0, result$returnCode, info = "Vector mode Integration unsuccessful!")
testthat::expect_equal(expected, result$integral, tolerance = tol, scale = abs(result$integral),
info = "Relative error not reached in vector mode")
testthat::expect_equal(expected, result$integral, tolerance = tol, scale = 1,
info = "Absolute error not reached in vector mode")
})
test_that("Test a Gaussian integral remapped to [0, infinity] limits", {
expected <- 1.00001
tol <- 1e-4
## Gaussian function
testFn1 <- function(x) {
val <- sum(((1 - x) / x)^2)
scale <- prod((2 / sqrt(pi)) / x^2)
exp(-val) * scale
}
result <- cubature::adaptIntegrate(f = testFn1,
lowerLimit = rep(0, 3),
upperLimit = rep(1, 3),
tol = tol)
testthat::expect_equal(0, result$returnCode, info = "Integration unsuccessful!")
testthat::expect_equal(expected, result$integral, tolerance = tol, scale = abs(result$integral),
info = "Relative error not reached")
testthat::expect_equal(expected, result$integral, tolerance = tol, scale = 1,
info = "Absolute error not reached")
testFn1_v <- function(x) {
val <- matrix(apply(x, 2, function(z) sum(((1 - z) / z)^2)), ncol(x))
scale <- matrix(apply(x, 2, function(z) prod((2 / sqrt(pi)) / z^2)), ncol(x))
exp(-val) * scale
}
result <- cubature::adaptIntegrate(f = testFn1_v,
lowerLimit = rep(0, 3),
upperLimit = rep(1, 3),
tol = tol,
vectorInterface = TRUE)
testthat::expect_equal(0, result$returnCode, info = "Vector mode Integration unsuccessful!")
testthat::expect_equal(expected, result$integral, tolerance = tol, scale = abs(result$integral),
info = "Relative error not reached in vector mode")
testthat::expect_equal(expected, result$integral, tolerance = tol, scale = 1,
info = "Absolute error not reached in vector mode")
})
test_that("Test volume of a hypersphere (integrating a discountinuous function)", {
expected <- 0.19728
tol <- 1e-4
## discontinuous objective: volume of hypersphere
testFn2 <- function(x) {
radius <- 0.50124145262344534123412
ifelse(sum(x * x) < radius * radius, 1, 0)
}
result <- cubature::adaptIntegrate(f = testFn2,
lowerLimit = rep(0, 2),
upperLimit = rep(1, 2),
tol = tol)
testthat::expect_equal(0, result$returnCode, info = "Integration unsuccessful!")
testthat::expect_equal(expected, result$integral, tolerance = tol, scale = abs(result$integral),
info = "Relative error not reached")
testthat::expect_equal(expected, result$integral, tolerance = tol, scale = 1,
info = "Absolute error not reached")
testFn2_v <- function(x) {
radius <- 0.50124145262344534123412
matrix(apply(x, 2, function(z) ifelse(sum(z * z) < radius * radius, 1, 0)), ncol = ncol(x))
}
result <- cubature::adaptIntegrate(f = testFn2_v,
lowerLimit = rep(0, 2),
upperLimit = rep(1, 2),
tol = tol,
vectorInterface = TRUE)
testthat::expect_equal(0, result$returnCode, info = "Vector mode Integration unsuccessful!")
testthat::expect_equal(expected, result$integral, tolerance = tol, scale = abs(result$integral),
info = "Relative error not reached in vector mode")
testthat::expect_equal(expected, result$integral, tolerance = tol, scale = 1,
info = "Absolute error not reached in vector mode")
})
test_that("Test a simple polynomial (product of coordinates)", {
expected <- 1
tol <- 1e-4
## product of coordinates
testFn3 <- function(x) prod(2 * x)
result <- cubature::adaptIntegrate(f = testFn3,
lowerLimit = rep(0, 3),
upperLimit = rep(1, 3),
tol = tol)
testthat::expect_equal(0, result$returnCode, info = "Integration unsuccessful!")
testthat::expect_equal(expected, result$integral, tolerance = tol, scale = abs(result$integral),
info = "Relative error not reached")
testthat::expect_equal(expected, result$integral, tolerance = tol, scale = 1,
info = "Absolute error not reached")
testFn3_v <- function(x) matrix(apply(x, 2, function(z) prod(2 * z)), ncol = ncol(x))
result <- cubature::adaptIntegrate(f = testFn3_v,
lowerLimit = rep(0, 3),
upperLimit = rep(1, 3),
tol = tol,
vectorInterface = TRUE)
testthat::expect_equal(0, result$returnCode, info = "Vector mode Integration unsuccessful!")
testthat::expect_equal(expected, result$integral, tolerance = tol, scale = abs(result$integral),
info = "Relative error not reached in vector mode")
testthat::expect_equal(expected, result$integral, tolerance = tol, scale = 1,
info = "Absolute error not reached in vector mode")
})
test_that("Test Gaussian centered at 1/2", {
expected <- 1
tol <- 1e-4
## guassian centered at 1/2
testFn4 <- function(x) {
a <- 0.1
s <- sum((x - 0.5)^2)
((2 / sqrt(pi)) / (2. * a))^length(x) * exp (-s / (a * a))
}
result <- cubature::adaptIntegrate(f = testFn4,
lowerLimit = rep(0, 2),
upperLimit = rep(1, 2),
tol = tol)
testthat::expect_equal(0, result$returnCode, info = "Integration unsuccessful!")
testthat::expect_equal(expected, result$integral, tolerance = tol, scale = abs(result$integral),
info = "Relative error not reached")
testthat::expect_equal(expected, result$integral, tolerance = tol, scale = 1,
info = "Absolute error not reached")
testFn4_v <- function(x) {
a <- 0.1
r <- apply(x, 2, function(z) {
s <- sum((z - 0.5)^2)
((2 / sqrt(pi)) / (2. * a))^length(z) * exp (-s / (a * a))
})
matrix(r, ncol = ncol(x))
}
result <- cubature::adaptIntegrate(f = testFn4_v,
lowerLimit = rep(0, 2),
upperLimit = rep(1, 2),
tol = tol,
vectorInterface = TRUE)
testthat::expect_equal(0, result$returnCode, info = "Integration unsuccessful!")
testthat::expect_equal(expected, result$integral, tolerance = tol, scale = abs(result$integral),
info = "Relative error not reached")
testthat::expect_equal(expected, result$integral, tolerance = tol, scale = 1,
info = "Absolute error not reached")
})
test_that("Test double Gaussian", {
expected <- 0.999994
tol <- 1e-4
## double guassian
testFn5 <- function(x) {
a = 0.1
s1 = sum((x - 1 / 3)^2)
s2 = sum((x - 2 / 3)^2)
0.5 * ((2 / sqrt(pi)) / (2. * a))^length(x) * (exp(-s1 / (a * a)) + exp(-s2 / (a * a)))
}
result <- cubature::adaptIntegrate(f = testFn5,
lowerLimit = rep(0, 3),
upperLimit = rep(1, 3),
tol = tol)
testthat::expect_equal(0, result$returnCode, info = "Integration unsuccessful!")
testthat::expect_equal(expected, result$integral, tolerance = tol, scale = abs(result$integral),
info = "Relative error not reached")
testthat::expect_equal(expected, result$integral, tolerance = tol, scale = 1,
info = "Absolute error not reached")
testFn5_v <- function(x) {
a = 0.1
r <- apply(x, 2, function(z) {
s1 <- sum((z - 1 / 3)^2)
s2 <- sum((z - 2 / 3)^2)
0.5 * ((2 / sqrt(pi)) / (2. * a))^length(z) * (exp(-s1 / (a * a)) + exp(-s2 / (a * a)))
})
matrix(r, ncol = ncol(x))
}
result <- cubature::adaptIntegrate(f = testFn5_v,
lowerLimit = rep(0, 3),
upperLimit = rep(1, 3),
tol = tol,
vectorInterface = TRUE)
testthat::expect_equal(0, result$returnCode, info = "Integration unsuccessful!")
testthat::expect_equal(expected, result$integral, tolerance = tol, scale = abs(result$integral),
info = "Relative error not reached")
testthat::expect_equal(expected, result$integral, tolerance = tol, scale = 1,
info = "Absolute error not reached")
})
test_that("Test Tsuda's example", {
expected <- 0.999998
tol <- 1e-4
## Tsuda's example
testFn6 <- function(x) {
a <- (1 + sqrt(10.0)) / 9.0
prod( a / (a + 1) * ((a + 1) / (a + x))^2)
}
result <- cubature::adaptIntegrate(f = testFn6,
lowerLimit = rep(0, 3),
upperLimit = rep(1, 3),
tol = tol)
testthat::expect_equal(0, result$returnCode, info = "Integration unsuccessful!")
testthat::expect_equal(expected, result$integral, tolerance = tol, scale = abs(result$integral),
info = "Relative error not reached")
testthat::expect_equal(expected, result$integral, tolerance = tol, scale = 1,
info = "Absolute error not reached")
testFn6_v <- function(x) {
a <- (1 + sqrt(10.0)) / 9.0
r <- apply(x, 2, function(z) prod( a / (a + 1) * ((a + 1) / (a + z))^2))
matrix(r, ncol = ncol(x))
}
result <- cubature::adaptIntegrate(f = testFn6_v,
lowerLimit = rep(0, 3),
upperLimit = rep(1, 3),
tol = tol,
vectorInterface = TRUE)
testthat::expect_equal(0, result$returnCode, info = "Integration unsuccessful!")
testthat::expect_equal(expected, result$integral, tolerance = tol, scale = abs(result$integral),
info = "Relative error not reached")
testthat::expect_equal(expected, result$integral, tolerance = tol, scale = 1,
info = "Absolute error not reached")
})
test_that("Test Morokoff & Calflish example", {
expected <- 1.00001
tol <- 1e-4
## Morokoff & Calflish function
testFn7 <- function(x) {
n <- length(x)
p <- 1/n
(1 + p)^n * prod(x^p)
}
result <- cubature::adaptIntegrate(f = testFn7,
lowerLimit = rep(0, 3),
upperLimit = rep(1, 3),
tol = tol)
testthat::expect_equal(0, result$returnCode, info = "Integration unsuccessful!")
testthat::expect_equal(expected, result$integral, tolerance = tol, scale = abs(result$integral),
info = "Relative error not reached")
testthat::expect_equal(expected, result$integral, tolerance = tol, scale = 1,
info = "Absolute error not reached")
testFn7_v <- function(x) {
matrix(apply(x, 2, function(z) {
n <- length(z)
p <- 1/n
(1 + p)^n * prod(z^p)
}), ncol = ncol(x))
}
result <- cubature::adaptIntegrate(f = testFn7_v,
lowerLimit = rep(0, 3),
upperLimit = rep(1, 3),
tol = tol,
vectorInterface = TRUE)
testthat::expect_equal(0, result$returnCode, info = "Integration unsuccessful!")
testthat::expect_equal(expected, result$integral, tolerance = tol, scale = abs(result$integral),
info = "Relative error not reached")
testthat::expect_equal(expected, result$integral, tolerance = tol, scale = 1,
info = "Absolute error not reached")
})
test_that("Test cubature web example", {
expected <- 13.69609
tol <- 1e-4
result <- cubature::adaptIntegrate(f = function(x) exp(-0.5 * sum(x^2)),
lowerLimit = rep(-2, 3),
upperLimit = rep(2, 3),
tol = tol)
testthat::expect_equal(0, result$returnCode, info = "Integration unsuccessful!")
testthat::expect_equal(expected, result$integral, tolerance = tol, scale = abs(result$integral),
info = "Relative error not reached")
testthat::expect_equal(expected, result$integral, tolerance = tol, scale = 1,
info = "Absolute error not reached")
result <- cubature::adaptIntegrate(f = function(x) matrix(
apply(x, 2, function(z) exp(-0.5 * sum(z^2))),
ncol = ncol(x)),
lowerLimit = rep(-2, 3),
upperLimit = rep(2, 3),
tol = tol,
vectorInterface = TRUE)
testthat::expect_equal(0, result$returnCode, info = "Integration unsuccessful!")
testthat::expect_equal(expected, result$integral, tolerance = tol, scale = abs(result$integral),
info = "Relative error not reached")
testthat::expect_equal(expected, result$integral, tolerance = tol, scale = 1,
info = "Absolute error not reached")
})
test_that("Test Wang-Landau sampling 1d example", {
expected <- 1.63564436296
tol <- 1e-5
## Numerical integration using Wang-Landau sampling
## Y. W. Li, T. Wust, D. P. Landau, H. Q. Lin
## Computer Physics Communications, 2007, 524-529
## Compare with exact answer: 1.63564436296
##
I.1d <- function(x) {
sin(4 * x) *
x * ((x * ( x * (x * x - 4) + 1) - 1))
}
result <- cubature::adaptIntegrate(f = I.1d,
lowerLimit = -2,
upperLimit = 2,
tol = tol)
testthat::expect_equal(0, result$returnCode, info = "Integration unsuccessful!")
testthat::expect_equal(expected, result$integral, tolerance = tol, scale = abs(result$integral),
info = "Relative error not reached")
testthat::expect_equal(expected, result$integral, tolerance = tol, scale = 1,
info = "Absolute error not reached")
I.1d_v <- function(x) {
matrix(apply(x, 2, function(z)
sin(4 * z) *
z * ((z * ( z * (z * z - 4) + 1) - 1))),
ncol = ncol(x))
}
result <- cubature::adaptIntegrate(f = I.1d_v,
lowerLimit = -2,
upperLimit = 2,
tol = tol,
vectorInterface = TRUE)
testthat::expect_equal(0, result$returnCode, info = "Integration unsuccessful!")
testthat::expect_equal(expected, result$integral, tolerance = tol, scale = abs(result$integral),
info = "Relative error not reached")
testthat::expect_equal(expected, result$integral, tolerance = tol, scale = 1,
info = "Absolute error not reached")
})
test_that("Test Wang-Landau sampling 2d example", {
expected <- -0.01797992646
tol <- 1e-5
## Numerical integration using Wang-Landau sampling
## Y. W. Li, T. Wust, D. P. Landau, H. Q. Lin
## Computer Physics Communications, 2007, 524-529
## Compare with exact answer: -0.01797992646
##
I.2d <- function(x) {
x1 <- x[1]; x2 <- x[2]
sin(4 * x1 + 1) * cos(4 * x2) * x1 * (x1 * (x1 * x1)^2 - x2 * (x2 * x2 - x1) +2)
}
result <- cubature::adaptIntegrate(f = I.2d,
lowerLimit = rep(-1, 2),
upperLimit = rep(1, 2),
tol = tol)
testthat::expect_equal(0, result$returnCode, info = "Integration unsuccessful!")
testthat::expect_equal(expected, result$integral, tolerance = tol, scale = abs(result$integral),
info = "Relative error not reached")
testthat::expect_equal(expected, result$integral, tolerance = tol, scale = 1,
info = "Absolute error not reached")
I.2d_v <- function(x) {
matrix(apply(x, 2,
function(z) {
x1 <- z[1]; x2 <- z[2]
sin(4 * x1 + 1) * cos(4 * x2) * x1 * (x1 * (x1 * x1)^2 - x2 * (x2 * x2 - x1) +2)
}),
ncol = ncol(x))
}
result <- cubature::adaptIntegrate(f = I.2d_v,
lowerLimit = rep(-1, 2),
upperLimit = rep(1, 2),
tol = tol,
vectorInterface = TRUE)
testthat::expect_equal(0, result$returnCode, info = "Integration unsuccessful!")
testthat::expect_equal(expected, result$integral, tolerance = tol, scale = abs(result$integral),
info = "Relative error not reached")
testthat::expect_equal(expected, result$integral, tolerance = tol, scale = 1,
info = "Absolute error not reached")
})
test_that("Test Multivariate Normal", {
expected <- 0.3341125
tol <- 1e-5
dmvnorm <- function (x, mean, sigma) {
x <- matrix(x, ncol = length(x))
distval <- stats::mahalanobis(x, center = mean, cov = sigma)
logdet <- sum(log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values))
logretval <- -(ncol(x) * log(2 * pi) + logdet + distval)/2
exp(logretval)
}
m <- 3
sigma <- diag(3)
sigma[2,1] <- sigma[1, 2] <- 3/5 ; sigma[3,1] <- sigma[1, 3] <- 1/3
sigma[3,2] <- sigma[2, 3] <- 11/15
result <- cubature::adaptIntegrate(f = dmvnorm,
lowerLimit = rep(-0.5, 3),
upperLimit = c(1, 4, 2),
tol = tol,
mean=rep(0, m), sigma=sigma)
testthat::expect_equal(0, result$returnCode, info = "Integration unsuccessful!")
testthat::expect_equal(expected, result$integral, tolerance = tol, scale = abs(result$integral),
info = "Relative error not reached")
testthat::expect_equal(expected, result$integral, tolerance = tol, scale = 1,
info = "Absolute error not reached")
dmvnorm_v <- function (x, mean, sigma) {
x <- t(x)
distval <- stats::mahalanobis(x, center = mean, cov = sigma)
logdet <- sum(log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values))
logretval <- -(ncol(x) * log(2 * pi) + logdet + distval)/2
logretval <- matrix(logretval, ncol = nrow(x))
exp(logretval)
}
result <- cubature::adaptIntegrate(f = dmvnorm_v,
lowerLimit = rep(-0.5, 3),
upperLimit = c(1, 4, 2),
tol = tol,
vectorInterface = TRUE,
mean=rep(0, m), sigma=sigma)
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.