Вычисление доверительного интервала для двойного интеграла с гауссовым взвешенным ⇐ Python

[Anonymous]

Требование : предоставить 96% доверительный интервал для .>
. Проблема
** координатное преобразование
определить: u = x и v = y - 1, тогда: d x d y = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =. d u d v .
Интеграл становится: выборка) v ) - это плотность n (0, i 2 ).
** Переписывание j как ожидание
Мы можем написать: .
Вычисление важности веса:
altage "/> altavation"/> altacation "/> altavation"/> altavation " src = "https://i.sstatic.net/jucb562c.png" />.образное. v i ) от n (0, i 2 ). Тогда оценка Монте -Карло j : .> ** Оценка стандартного отклоненности Определение определения "src =" https://i.sstatic.net/zin5ykms.png "/>, Среднее значение выборки: , выборка: . 96% доверительный интервал
Предполагая z 0,98 ≈ 2,05 (с α = 0,04, поскольку 96% доверительный интервал соответствует /2 = 0,02 в каждом хвосте), достоверное интерв. alt = "доверительный интервал" src = "https://i.sstatic.net/vaoju2th.png" />.
* Реализация кода Python приведенного выше конструкции Monte Carlo
python Code приведенного выше конструкции Monte Carlo
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()
< /code>

Надеюсь, кто -то сможет подтвердить мое решение, а также улучшить его, если это возможно. < /p>

Подробнее здесь: https://stackoverflow.com/questions/794 ... e-integral