I have this Python code that implements a rectangular numerical integration. It evaluates the (K-1)-dimensional integral for arbitrary integer \$K \geq 1\$
$$\int_{u_K = 0}^{\gamma}\int_{u_{K-1} = 0}^{\gamma-u_K}\cdots\int_{u_2}^{\gamma-u_K-\cdots-u_3}F_U(\gamma-\sum_{k=2}^Ku_k)f_U(u_2)\cdots f_U(u_K)\,du_2\cdots du_K$$
where \$F_U\$ corresponds to the cdf function in the code, and \$f_U\$ to the pdf. I implemented it using recursion as follows:
#************************** Import necessary libraries******************************************
import numpy as np
import matplotlib.pyplot as plt
import time
#******************************Set the constant scalars and vectors***************************
start_time = time.time()
KU = 3
eta_dB = 3
eta = 10**(eta_dB/10)
ExpanF = 50
tdB = np.arange(-5,11,4)
tVec = 10**(tdB/10)
thVec = (tVec/eta)*(ExpanF-eta*(KU-1))
N = 10000 # For simulation
du = 0.01
#******************************Define functions to be used***************************************
#Define the CDF of U
def CDF(u):
return 1-1/(u+1)
#Define the PDF of U
def pdf(u):
return 1/((1+u))**2
def FK(h, k):
#print(f'h is {h}, and k is {k}')
if k == 1:
res = CDF(h)
else:
#n_iter = int(h/du)
res = 0
u = 0
while u < h:
res += FK(h-u, k-1)*pdf(u)*du
u += du
return res
#*******************Find the numerical and simulation values of the integral******************
ResultNum = []
ResultSim = []
if (ExpanF-eta*(KU-1)) > 0:
for t in thVec:
# Numerical
ResultNum.append(1-FK(t, KU))
# Simulation
count = 0
for n in range(0, N):
if np.sum(np.random.exponential(1, (1, KU))/np.random.exponential(1, (1, KU))) <= t:
count += 1
ResultSim.append(1-count/N)
The parameter du
is the resolution of the numerical integral, and the smaller, the more accurate the numerical approximation is. However, decreasing du
comes at the expense of more computational time using the above code. For example I tried to run the above code on my machine Intel i5 @2.5 GH with 6 GB of RAM for du=0.0001
, and it took 15 hours and didn't finish, at which case I had to abort it.
Is there anyway the above code can be optimized to run significantly faster? I tried C++, and I got a 30x speed factor using the same code and algorithm, but I was wondering if I can get a similar or better speed up factor in Python.
scipy.integrate.nquad
? \$\endgroup\$func(x0, x1, ..., xn, t0, t1, ..., tm)
. Integration is carried out in order. That is, integration overx0
is the innermost integral, andxn
is the outermost." \$\endgroup\$