draw sample from distribution in interval

[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?