draw sample from distribution in interval

[SchroedingersLion]

Greetings,

I am trying to numerically calculate a double integral over probability distributions (beta distributions), see picture.
I want to approximate the integrals as sums, taking samples according to the distributions.

I know I can sample from the beta function via beta.rvs(a,b), but since my second integration goes only from 0 to theta1, I can only use samples that are < theta1 there.

How can I solve this problem efficiently?



Regards!

[SchroedingersLion]

from scipy.stats import beta

n=10000 #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()

#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)

This is my current code. At first I sample n theta1 values from the first distribution.
I save the values of the beta function at these points in a list.

For each theta1 (first for loop), I sample n theta2 values from the second distribution.
However, I only consider those theta2 that are < theta1 (if structure).
These theta2 get inserted into the beta function and the values are saved in a list.
I sum over these function values (divided by number of samples and multiplied by x, namely theta1) in order to get the approximation to the "inner" integral for a given theta1.

In the end, I get my final result by summing over the products of the values of the beta function at the theta1 points with the respective inner integrals.

Any comments?

[SchroedingersLion]

Ah, I got my error. Aside from my inexperienced Python coding, my mathematical method was flawed.