- #1

- 195

- 52

Greetings,

I am desparately trying to solve a double integral via Monte Carlo integration. I integrate over two probability densities, the Beta distributions. Python can easily draw samples from these densities and calculate its function values. The integral can be seen here:

Now my idea was to sample n values ## \theta_1 ## from the first Beta function and insert it into the function.

For each ## \theta_1 ## I would then calculate the inner integral through sampling n ## \theta_2 ## from the second Beta function.

I would then calculate the entire integral via:

## \frac{1}{n} \sum_{i=1}^n ## Beta(##\theta_1^{i}## | 91,11) *## C_i##

with ##C_i=\frac{\theta_1^{i}}{m} \sum_{j=1}^m ## Beta(##\theta_2^{j}## | 3,1)

However, I could only consider those ## \theta_2 ## with ##\theta_2 < \theta_1 ##.

At first, I used an if-structure (commented part of my code) to store only the fitting m ##\theta_2## in a list for the calculation of the integral. But I figured that ##\theta_2 < \theta_1## is often not fulfilled, so the approximation of the inner integral would need much larger sample sizes. It did not lead to correct results.

Then I tried to approximate the inner integral through random sampling. For each ##\theta_1##, just sample n random ##\theta_2 \in (0, \theta_1) and approximate the integral with the described sum above.

But I get a result of 6.71...??

As a probability, the right result should be < 1. In fact, I know it to be around 0.7.

Can someone see my error?

Regards!

I am desparately trying to solve a double integral via Monte Carlo integration. I integrate over two probability densities, the Beta distributions. Python can easily draw samples from these densities and calculate its function values. The integral can be seen here:

Now my idea was to sample n values ## \theta_1 ## from the first Beta function and insert it into the function.

For each ## \theta_1 ## I would then calculate the inner integral through sampling n ## \theta_2 ## from the second Beta function.

I would then calculate the entire integral via:

## \frac{1}{n} \sum_{i=1}^n ## Beta(##\theta_1^{i}## | 91,11) *## C_i##

with ##C_i=\frac{\theta_1^{i}}{m} \sum_{j=1}^m ## Beta(##\theta_2^{j}## | 3,1)

However, I could only consider those ## \theta_2 ## with ##\theta_2 < \theta_1 ##.

At first, I used an if-structure (commented part of my code) to store only the fitting m ##\theta_2## in a list for the calculation of the integral. But I figured that ##\theta_2 < \theta_1## is often not fulfilled, so the approximation of the inner integral would need much larger sample sizes. It did not lead to correct results.

Then I tried to approximate the inner integral through random sampling. For each ##\theta_1##, just sample n random ##\theta_2 \in (0, \theta_1) and approximate the integral with the described sum above.

But I get a result of 6.71...??

As a probability, the right result should be < 1. In fact, I know it to be around 0.7.

Can someone see my error?

Regards!

Python:

```
from scipy.stats import beta
import random
n=3000 #sample sizes
sample_theta1=[beta.rvs(91,11) for x in range(n)]
beta_values_theta1=beta.pdf(sample_theta1,91,11)
summands_total_integral=[]
sample_theta2=[]
beta_values_theta2=[]
integrals_theta2=[]
#calculate integral over theta2 for each theta1
for x in sample_theta1:
# sample_theta2=[beta.rvs(3,1) for z in range(n)]
# for z in range(n):
# if sample_theta2[z]<x:
# beta_values_theta2.append(beta.pdf(sample_theta2[z],3,1))
# integrals_theta2.append(float(x) / len(beta_values_theta2) * sum(beta_values_theta2))
# beta_values_theta2.clear()
sample_theta2=[random.uniform(0, x) for z in range(n)]
beta_values_theta2=beta.pdf(sample_theta2, 3, 1)
integrals_theta2.append(float(x) / n * sum(beta_values_theta2))
#calculate summands for approximation of theta1 integral
for x in range(n):
summands_total_integral.append(beta_values_theta1[x] * integrals_theta2[x])
result=float(1) / n * sum(summands_total_integral)
print(result)
```