Estimating an integral

Want to compute the integral

\[ I(g) = \int_0^1 g(x) dx \,, \] for some function \(g:[0,1] \to \Re\).

Estimating an integral

Hence, if \(X_1, \ldots, X_n\) are \(\mathcal{U}(0,1)\), and \(E(|g(X_i)|) < \infty\), then \[ \lim_\limits{n \to \infty} \frac{1}{n} \sum_{i=1}^n g(X_i) \overset{a.s.}{=} E(g(X_i)) = \int_0^1 g(x) dx \,. \]

Monte Carlo approximation

Draw a large sample \(X_1,\ldots, X_n \sim \mathcal{U}(0,1)\) and approximate

\[ \int_0^1 g(x) dx \approx \frac{1}{n} \sum_{i=1}^n g(X_i) \,. \]

Uncertainty quantification

Our approximation \[ \bar{I}(g) = \frac{1}{n} \sum_{i=1}^n g(X_i) \] for \(I(g) = \int_0^1 g(x)dx\) is a random variable. To make statements about the precision of our approximation, we can appeal to the central limit theorem.

Uncertainty quantification

Hence, for \(Y_i = g(X_i)\), \(X_i \sim \mathcal{U}(0,1)\), we get that

\[ \Pr \left( a < \frac{\bar{I}(g) - I(g)}{\sigma / \sqrt{n}} < b \right) \longrightarrow \Psi(b) - \Psi(a) \quad n \to \infty \,, \] where \(\sigma = \int_0^1 (g(x) - I(g))^2 dx\).

Again using the Monte Carlo method, we can approximate \(\sigma\) as \[ \sigma^2 \approx \frac{1}{n} \sum_{i = 1}^n g(X_i)^2 - \bar{I}(g)^2 \,. \]

Integration by substitution

Want to compute an integral with integration bounds \(a,b \in [-\infty, \infty]\), i.e.  \[\int_a^b g(x) dx \,.\]

Recipe

• \(y:[a, b] \to [0,1]\) with \(\lim_\limits{x \to a} y(x) = 0\), \(\lim_\limits{x \to b} y(x) = 1\)
• solve for \(x(y)\)
• compute \(x'(y) = \frac{dx}{dy}\)
• substitue \(x(y)\) for \(x\), \(x'(y) dy\) for \(dx\)
• adjust inegration bounds, e.g. \(a \mapsto \lim_\limits{x \to a} y(x) = 0, \, b \mapsto \lim_\limits{x \to b} y(x) = 1\)

Integration by substitution

Then \[ \int_a^b g(x) dx = \int_{0}^{1} g(x(y)) x'(y) dy \,. \]

Question 1

1.a

Create a vector of \(10^5\) random numbers \(X_i \sim \mathcal{U}(0,1)\)

# Sample size 
N <- 10^5 

# Sampling uniform random variables 
set.seed(123)
X <- runif(N) 

Question 1

1.b

# Histogram 
hist(X, freq = FALSE) 
abline(c(1,0), col = "red", lty = 2, lwd = 3)

# All but last X 
X1 <- X[-N] 
# All but first X 
X2 <- X[-1] 

# Scatter plot 
plot(X1, X2, pch = 20, cex = 0.1, col = "royalblue1") 
plot_title <- expression("Scatter plot " ~ X[i] * "," ~ X[i + 1])
title(main = plot_title)

# ACF function plots by default 
acf(X) 

Question 1

1.c

Calculate the mean, median and variance of the sample

# mean of X 
cat("mean: ", mean(X)) 

mean:  0.4992992

# median of X 
cat("median: ", median(X))

median:  0.4972005

# variance of X 
cat("variance: ", var(X)) 

variance:  0.08317841

1.d

Estimate the probability \(\Pr(X_i > 0.6)\).

# estimate pr(x > 0.6) 
pr <- mean(X > 0.6)
cat("Pr(X_i > 0.6) ≈ ", pr)

Pr(X_i > 0.6) ≈  0.39867

Question 2

If \(y = \exp\{-x\}\) calculate the integral \[ \int_0^\infty g(x) dx = \int_0^1 \frac{1}{y} g(-\log(y)) dy \]

where \(g(-y) = \exp \left\{-\log(y)^2 / 2 \right\}\). Run the algorithm for \(10^4\) simulations.

# setting seed 
set.seed(123) 
# sample size 
Nsim <- 10^4
# draw sample 
y <- runif(Nsim)
# calculate function within integral 
x <- (1/y)*exp((-log(y)^2)/2)
# take mean
cat("Integral approximation: ", mean(x))

Integral approximation:  1.256539

# approximation error 
cat("Approximation errror: ", mean(x) - sqrt(pi / 2)) 

Approximation errror:  0.003224685

Question 3

Estimate the value of the following integrals. In all cases you have to give a measure of how accurate your approximation is such as the standard error or a \(95\%\) confidence interval.

3.a

\[ \int_0^1 \exp\{-x^{1.5}\} dx \]

Question 3

3.a Integral approximation

\[ \begin{aligned} \int_0^1 \exp\{-x^{1.5}\} dx &= \mathrm{E}_{\mathcal{U}(0,1)} \left[ \exp\{-x^{1.5}\} \right] \\ & \overset{\textrm{LLN}}{\approx}\frac{1}{n} \sum_{i = 1}^n \exp\{-X_i^{1.5}\}\,, \end{aligned} \]
where \(X_i\) are i.i.d. \(\mathcal{U}(0,1)\).

Question 3

3.a 95% confidence interval

• Central limit theorem: As \(n \to \infty\): \[ \Pr\left( \left| \frac{\bar{I}(g) - I(g)}{\sigma / \sqrt{n}} \right| < z \right) \longrightarrow 2 \Psi(z) - 1\,, \] for \(\sigma = \int_0^1 (g(x) - I(g))^2 dx\).

• Finding \(z\): \[ 2 \Psi(z) - 1 = 0.95 \iff \Psi(z) = \frac{1 + 0.95}{2} \iff z = \Psi^{-1}(0.975) \]

Therefore \[ \Pr\left( \bar{I}(g) + z \frac{\sigma}{\sqrt{n}}< I(g) < \bar{I}(g) + z \frac{\sigma}{\sqrt{n}} \right) \approx 0.95 \,, \] for large \(n\).

Question 3

3.a

# set seed 
set.seed(1) 
# sample size 
N <- 10^5 
# uniform sample 
X <- runif(N) 
# function within integral 
y <- exp(-X^(1.5)) 
# integral approximation 
int_approx <- mean(y) 
# estimating standard error 
se <- sqrt(mean(y^2) - int_approx^2) / sqrt(N) 
# confidence interval 
conf_level <- 0.95 
z <- qnorm((1 + conf_level) / 2)  
uci <- int_approx + z * se 
lci <- int_approx - z * se 
# real integral value 
int_val = 0.699792 
# approximation error 
approx_error <- int_val - int_approx
print(rbind(int_val, int_approx, approx_error, se, lci, uci))

                      [,1]
int_val       0.6997920000
int_approx    0.6999995477
approx_error -0.0002075477
se            0.0006244706
lci           0.6987756079
uci           0.7012234875

Question 3

3.b.

\[ \int_0^\infty \exp\{-x^{1.5}\} dx\,, \] try various transformation for this to find out which one works best.

Question 3

3.b. Integration by substitution

• \(y = \frac{1}{1 + x}\)
• \(x = \frac{1}{y} - 1\)
• \(\frac{d}{dy} \frac{1}{y} - 1 = - \frac{1}{y^2}\)
• \(\lim_\limits{x \to 0} \frac{1}{1 + x} = 1\), \(\lim_\limits{x \to \infty} \frac{1}{1 + x} = 0\)

Therefore

\[ \begin{aligned} \int_0^\infty \exp\{-x^{1.5}\} dx &= - \int_1^0 \frac{1}{y^2} \exp\left\{-\left(\frac{1}{y} - 1\right)^{1.5}\right\} dy \\ &= \int_0^1 \frac{1}{y^2} \exp\left\{-\left(\frac{1}{y} - 1\right)^{1.5}\right\} dy \end{aligned} \]

Question 3

3.b. What makes a good transformation?

# plotting transformed integrand with y = 1 / (1 + x) on unit interval 
ys <- seq(from = 0.00001, to = 1, length.out = 1000) 
integrand <- 1 / (ys^2) * exp(-(1 / ys - 1)^(1.5)) 
par(mar = c(5, 4, 4, 12))
plot(ys, integrand, type = "l", lwd = 2, ylim = c(0, 2))

# plot another transoformation y = exp(-x^1.5)
other_integrand <- 2 / 3 * (-log(ys))^(-1/3)
lines(ys, other_integrand, type = "l", lwd = 2, col = "green") 

# plot another transoformation y = exp(-x)
other_integrand2 <- 1 / ys * exp(-(-log(ys))^(1.5)) 
lines(ys, other_integrand2, type = "l", lwd = 2, col = "cyan") 

# add hline at value of integral 
f <- function(x){(1/x^2)*exp(-(1/x-1)^(1.5))}
Integral = integrate(f,0,1)$value 
abline(Integral, 0, lwd = 2, lty = 2, col = "red")

# add legend 
legend("right", 
       legend = c(expression(y == frac(1, 1+x)), 
                expression(y == exp(-x^1.5)), 
                expression(y == exp(-x)),  
                "Integral value"), 
       col = c("black", "green", "cyan", "red"), 
       lty = c(1, 1, 1, 2),          
       lwd = 2, 
       inset = -0.35, 
       xpd = TRUE)

set.seed(123) 
N <- 10^5
X <- runif(N) 
y1 <- 1 / (X^2) * exp(-(1 / X - 1)^(1.5)) 
y2 <- 2 / 3 * (-log(X))^(-1/3)
y3 <- 1 / X * exp(-(-log(X))^(1.5)) 

se1 <- sqrt(mean(y1^2) - mean(y1)^2) / sqrt(N) 
se2 <- sqrt(mean(y2^2) - mean(y2)^2) / sqrt(N) 
se3 <- sqrt(mean(y3^2) - mean(y3)^2) / sqrt(N)

ci1 <- mean(y1) + c(-z * se1, z * se1)
ci2 <- mean(y2) + c(-z * se2, z * se2)
ci3 <- mean(y3) + c(-z * se3, z * se3)

cat("y = 1 / (1 + x) \nEstimate:", mean(y1), "\nError: ", mean(y1) - Integral, "\nCI: ", ci1) 

y = 1 / (1 + x) 
Estimate: 0.9015331 
Error:  -0.001212228 
CI:  0.8976504 0.9054157

cat("y = exp(-x^1.5) \nEstimate:", mean(y2), "\nError: ", mean(y2) - Integral, "\nCI: ", ci2) 

y = exp(-x^1.5) 
Estimate: 0.901386 
Error:  -0.001359323 
CI:  0.8977233 0.9050486

cat("y = exp(-x) \nEstimate:", mean(y3), "\nError: ", mean(y3) - Integral, "\nCI: ", ci3) 

y = exp(-x) 
Estimate: 0.9026241 
Error:  -0.0001211551 
CI:  0.9005703 0.9046779

Question 3

3.c.

\[ \int_{-\infty}^\infty \frac{1}{1 + x^2} dx \] This is known to be \(\pi\)

Question 3

3.c. Integration by substitution

Before thinking of transforms, let us simplify things

• Symmetry: \[ \int_{-\infty}^\infty \frac{1}{1 + x^2} dx = 2 \int_{0}^\infty \frac{1}{1 + x^2} dx \]

• Split into parts: \[ \int_{0}^\infty \frac{1}{1 + x^2} dx = \int_{0}^1 \frac{1}{1 + x^2} dx + \int_{1}^\infty \frac{1}{1 + x^2} dx \]

Question 3

3.c. Integration by substitution

Hence, we need to find a transformation for \[\int_{1}^\infty \frac{1}{1 + x^2} dx\]

Let \(y = \frac{1}{x}\):

• \(x = \frac{1}{y}\)
• \(\frac{d}{dy} x(y) = - \frac{1}{y^2}\)
• \(\lim_\limits{x \to 1} y(x) = 1\), \(\lim_\limits{x \to \infty} y(x) = 0\)

Question 3

3.c. Integration by substitution

Thus \[ \begin{aligned} \int_{1}^\infty \frac{1}{1 + x^2} dx &=- \int_{1}^0 \frac{1}{1 + (1/y)^2} \frac{1}{y^2} dy \\ &= \int_{0}^1 \frac{1}{1 + (1/y)^2} \frac{1}{y^2} dy \\ &= \int_{0}^1 \frac{1}{\frac{y^2 + 1}{y^2}}\frac{1}{y^2} dy \\ &= \int_{0}^1 \frac{y^2}{1 + y^2} \frac{1}{y^2} dy \\ &= \int_{0}^1 \frac{1}{1 + y^2} dy \\ \end{aligned} \]

Question 3

3.c. Integration by substitution

Hence, we have found that \[ \begin{aligned} \int_{-\infty}^\infty \frac{1}{1 + x^2} dx &= 2 \int_{0}^\infty \frac{1}{1 + x^2} dx \\ &= 2 \left( \int_{0}^1 \frac{1}{1 + x^2} dx + \int_{1}^\infty \frac{1}{1 + x^2} dx \right) \\ &= 2 \left( \int_{0}^1 \frac{1}{1 + x^2} dx + \int_{0}^1 \frac{1}{1 + y^2} dy \right) \\ &= 4 \int_{0}^1 \frac{1}{1 + y^2} dy \\ &= 4 \textrm{E}_{\mathcal{U}(0,1)} \left[\frac{1}{1 + X^2} \right] \\ & \approx 4 \frac{1}{n} \sum_{i=1}^n \frac{1}{1 + X_i^2}, \quad X_i \overset{i.i.d.}{\sim}\mathcal{U}(0,1) \end{aligned} \]

Question 3

3.c. Geometric solution

x <- seq(0, 1, length.out = 500)
f <- sqrt(1 - x^2)

plot(x, f, type = "l", xlab = "x", ylab = "", xlim = c(0, 1), ylim = c(0, 1), lwd = 2, col = "black")
polygon(c(x, rev(x)), c(f, rep(0, length(f))), col = "lightgrey")

legend("topright", legend = c(expression(sqrt(1 - x^2)), expression(frac(pi, 4))),
       lwd = c(2, NA), pch = c(NA, 15), col = c("black", "lightgrey"))

Question 3

3.c.

set.seed(123) 
N <- 10^5 
X <- runif(N) 

y1 <- 4 / (1 + X^2)
y2 <- 4 * sqrt(1 - X^2) 
integral <- pi 

se1 <- sqrt(mean(y1^2) - mean(y1)^2) / sqrt(N) 
se2 <- sqrt(mean(y2^2) - mean(y2)^2) / sqrt(N) 

ci1 <- mean(y1) + c(-z * se1, z * se1)
ci2 <- mean(y2) + c(-z * se2, z * se2)

cat("y = 1 / (1 + x^2) \nEstimate:", mean(y1), "\nError: ", mean(y1) - integral, "\nCI: ", ci1) 

y = 1 / (1 + x^2) 
Estimate: 3.143417 
Error:  0.001823983 
CI:  3.139433 3.1474

cat("y = sqrt(1 - x^2) \nEstimate:", mean(y2), "\nError: ", mean(y2) - integral, "\nCI: ", ci2) 

y = sqrt(1 - x^2) 
Estimate: 3.143952 
Error:  0.002359275 
CI:  3.138422 3.149481

Question 3

3.d.

\[ I(g) = \int_0^1 \int_0^1 \int_0^1 \int_0^1 \exp\{-\sqrt{xywz}\} dx dy dz dw \]

Note that \[ \begin{aligned} I(g) = \textrm{E}_{\mathcal{U}(0,1)^4} \left[ \exp\{-\sqrt{XYWZ}\} \right] \approx \frac{1}{n} \sum_{i = 1}^n \exp\{-\sqrt{X_i Y_i W_i Z_i}\}\,, \end{aligned} \] where \(X_i, Y_i, W_i, Z_i\) are independent, \(\mathcal{U}(0,1)\) random variables

set.seed(123) 
N <- 10^5
X <- runif(N)
Y <- runif(N) 
Z <- runif(N) 
W <- runif(N)  
xywz <- exp(-sqrt(X * Y * Z * W))
int_approx <- mean(xywz) 
se <-sqrt(mean(xywz^2) - int_approx^2) / sqrt(N) 
ci <- mean(int_approx) + c(-z * se, z * se)
cat("Estimate:", int_approx, "\nCI: ", ci) 

Estimate: 0.8299637 
CI:  0.8292357 0.8306917