Implementing numerical integration in Python
I have this Python code that implements a rectangular numerical integration.


In what sense "rectangular"? Or is that as opposed to contour 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.


It would probably help to rewrite it in fully symmetric form: $$\underset{\substack{u_i \ge 0 \\ u_1 + \cdots + u_K \le \gamma}}{\int \cdots \int} f_U(u_1)\cdots f_U(u_K) \textrm{d}u_1 \cdots \textrm{d}u_K$$

Then you can exploit the symmetry to divide the polytope of interest into equivalent chunks and only integrate each chunk once.




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



So far I understand.


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



What do any of those mean? "constant scalars and vectors" is not an adequate description.


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



If \$f_U\$ is that simple, it might make more sense to tackle it with symbolic integration rather than quadrature.


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 &lt; h:
            res += FK(h-u, k-1)*pdf(u)*du
            u += du
    return res



This could use some comments. It seems to be first order quadrature. If you use a higher order, maybe a larger step size will give you equally accurate results much faster.


#*******************Find the numerical and simulation values of the integral******************
ResultNum = []
ResultSim = []

if (ExpanF-eta*(KU-1)) &gt; 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))) &lt;= t:
                count += 1
        ResultSim.append(1-count/N)



I'm missing a lot of context here.