# Numerical Integration

## Main Question or Discussion Point

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.