[Python] Computing a Confidence Interval for a Gaussian-Weighted Double Integral

Requirement: Provide a 96% confidence interval for .


Here is my solution to this problem:

* Transforming the expectation problem

** Coordinate transformation

Define: u = x and v = y - 1, then: dxdy = dudv.

The integral becomes: .

** Choosing a proposal distribution (importance sampling)

A convenient choice is to use the independent bivariate normal distribution with density: .

Thus, g(u, v) is the density of N(0, I2).

** Rewriting J as an expectation

We can write: .

Computing the importance weight:

Thus: .

* Monte Carlo algorithm

Suppose we generate N independent samples (ui, vi) from N(0, I2). Then the Monte Carlo estimate of J is: .

** Estimating the standard deviation

Define: , Sample mean: , Sample variance: .

The standard error of the estimate is: .

** 96% confidence interval

Assuming z0.98 ≈ 2.05 (with α = 0.04, since a 96% confidence interval corresponds to α / 2 = 0.02 in each tail), the 96% confidence interval for J is: .

* Python code implementation of the above Monte Carlo algorithm

import random
import math

def estimate_J(N):
"""
Estimate J using the Monte Carlo method with the proposed distribution N(0, I_2).

Pseudocode:
sum_Y = 0
sum_Y2 = 0
Loop i from 1 to N:
Generate (u, v) from the 2D normal distribution N(0, I_2)
Compute Y = 2 * π * exp( (u * v) / 8 )
sum_Y = sum_Y + Y
sum_Y2 = sum_Y2 + Y^2
mean_Y = sum_Y / N
variance_Y = (sum_Y2 - N * mean_Y^2) / (N - 1)
s_Y = sqrt(variance_Y)

I_hat = mean_Y // Estimate of J
SE = s_Y / sqrt(N)
z = 2.05 // 0.98 quantile of the normal distribution
CI_lower = I_hat - z * SE
CI_upper = I_hat + z * SE

Return (I_hat, CI_lower, CI_upper)
"""

sum_Y = 0.0
sum_Y2 = 0.0

for i in range(N):
# Generate u, v from the normal distribution N(0,1)
u = random.gauss(0, 1)
v = random.gauss(0, 1)
# Compute Y = 2 * π * exp((u * v) / 8)
Y = 2 * math.pi * math.exp((u * v) / 8)
sum_Y += Y
sum_Y2 += Y ** 2

mean_Y = sum_Y / N
variance_Y = (sum_Y2 - N * (mean_Y ** 2)) / (N - 1)
s_Y = math.sqrt(variance_Y)

# Estimate J
I_hat = mean_Y
SE = s_Y / math.sqrt(N)
z = 2.05 # 0.98 quantile of the normal distribution
CI_lower = I_hat - z * SE
CI_upper = I_hat + z * SE

return I_hat, CI_lower, CI_upper

def main():
N = 1000000 # Number of samples, can be adjusted for desired accuracy
I_hat, CI_lower, CI_upper = estimate_J(N)
print("Estimated J =", I_hat)
print("96% confidence interval: [{}, {}]".format(CI_lower, CI_upper))

if __name__ == "__main__":
main()


I hope anyone can validate my solution, and also improve it if possible.

Continue reading...