# Nested Integrals

• MATLAB

## Main Question or Discussion Point

Hello,

How can I evaluate the following nested integral in MATLAB for a general value of $K$

$${\int\limits_{u_1=0}^{\gamma}\int\limits_{u_2=0}^{\gamma-U_1}\cdots \int\limits_{u_{K}=0}^{\gamma-\sum_{k=1}^{K-1}U_k}}f_{U_1}(u_1)f_{U_2}(u_2)\cdots f_{U_{K}}(u_{K})\,du_1du_2\cdots du_{K}$$

where $f_{U_k}(u_k)$ is the PDF of the random variable $U_k$, and $\gamma\geq 0$ is a constant.

Thanks in advance

## Answers and Replies

Related MATLAB, Maple, Mathematica, LaTeX News on Phys.org
This is for double integrals. I am looking for a K nested integrals.

Any hints?

kreil
Gold Member
I'm not familiar with the exact application you have, and you might be able to use properties of the PDF, or knowledge about whether the random variables are independent, to simplify things. That said...

Do you have the function values as vectors or matrices? If so you can use trapz and nest the integrals to your hearts content. For a given value of K the array of function values would have K dimensions, and each integral evaluates over one dimension, so the final integral reduces a vector to a scalar.

Otherwise, if you have the functional expressions for the different $f_{U_k}$, then you would need to use integral/integral2/integral3, which limits how many nested numerical integrations you could do.

The third option is to use int() in Symbolic Math Toolbox which lets you calculate the answer symbolically. You should be able to nest as many integrals as you want with that.

Unfortunately, I couldn't simply the results numerically, or write it as a product of one-dimensional integrals. The random variable $U_k$ are independent, but the integrations limits dependent on other random variables, so they are inseparable. I am looking to write a piece of code in MATLAB, that allows me to evaluate the K nested integrals without changing the code for each value of K. Could you elaborate on the last method by using int()?

kreil
Gold Member
Here's an example that does a 4-D integral

$$\int_0^5 \int_0^5 \int_0^5 \int_0^5 \left( w^2 + x^2 + y^2 + z^2 \right) \mbox{dw} \, \mbox{dx} \,\mbox{dy} \,\mbox{dz}$$

Code:
>> syms w x y z
>> f = w^2 + x^2 + y^2 + z^2;
>> dw = int(f,'w',0,5)
dw =

5*x^2 + 5*y^2 + 5*z^2 + 125/3

>> dx = int(dw,'x',0,5)
dx =

25*y^2 + 25*z^2 + 1250/3

>> dy = int(dx,'y',0,5)
dy =

125*z^2 + 3125

dz = int(dy,'z',0,5)
dz =

62500/3
You can do this for definite OR indefinite integrals, and you can nest as many calls as necessary. You can implement this as a loop that runs K times to do the K integrals, with the first pass integrating the defined function, and the subsequent passes integrating the result of the previous pass.

That's interesting. Yes, I am looking to write a for loop to evaluate the integrals to automate the process. I will check this way out. Thanks

I can write a for loop for the example you provided, because the limits are independent of each others. This takes only one for loop. But in my case, the limits are interdependent. How to handle this using this approach? I suspect we need $K-1$ for loops to iterate over all possible values of the limits, which is not efficient.

andrewkirk
Science Advisor
Homework Helper
Gold Member
I can write a for loop for the example you provided, because the limits are independent of each others. This takes only one for loop. But in my case, the limits are interdependent. How to handle this using this approach? I suspect we need $K-1$ for loops to iterate over all possible values of the limits, which is not efficient.
First note that the above integral is written incorrectly. The order of the $du_k$s on the right should be the reverse of the order of the integrals on the left. A last-in-first-out principle applies. I'll assume the intended order is that of the integrals on the left.

You can write your integral with independent limits, with every upper limit being $\gamma$, by defining the integrand to be
• what is written above if $u_j<\gamma - \sum_{k=1}^{j-1}u_k$ for all $j\in 1,...,K$, where we set $u_0=0$, and
• zero otherwise.
What that does is implicitly apply the limits by setting the integrand to zero whenever one or more of the limits is breached.

Prima facie, I would expect that to be less efficient than manually programming the integration by adding up the areas of incremental rectangles, because that can use the interdependent limits to dramatically reduce the number of iterations needed.

But it is possible that MATLAB has some super-efficient routines for numerical integration when you use its built-in functions (eg they would use compiled code rather than interpreted code) that would more than offset manually programming the integration. I'd be surprised if that were the case, but it may be.

For the manual programming, you need to define the nested integral as a recursive function. The mathematics of this is as follows:

EDIT: The following equations and the code that follows contain some minor, but important corrections, plus a simplification (reduction in number of arguments to the recursive function) that were made on 19/7/18.

We can write this using a recursive function $G:\mathbb R\times \mathbb Z_+\rightarrow\mathbb R$, defined as:
$$G(h,m) = \begin{cases} 1 &\mbox{if } m=0 \\ \int_0^{h} G(h - u_m,m-1)\,f_{U_m}(u_m)\,du_m & \mbox{otherwise } \end{cases}$$

The full nested integral is
$$G(\gamma,K)$$

This reverses the order of nesting of the integrals, but that doesn't change the result, since the random variables are all independent.

The Python code would be:
Code:
granularity = 1000
gamma = # set value of gamma ....
K = # set value of K ....

# assume that all the f_{U_k} are the same, so we only need one definition:
def f_U(u)
# define function f_U here

def integ(h, m)
if m == 0
res = 1
else
du = h / granularity
res, u = 0
# integrate across this u_m
while u < h
u += du
res += integ(h - u, m - 1) * f(u) * du
return res
# end of definition of function integ

# call the function to find the nested integral
prob = integ(gamma, K)
print(prob)

Last edited:
First note that the above integral is written incorrectly. The order of the $du_k$s on the right should be the reverse of the order of the integrals on the left. A last-in-first-out principle applies. ...
I think that the limits of the integrals determine the order. But you are right, it is more accurate to reverse the order of the integrands, because the derivation shows them that way (as illustrated below).

For the implementation, I still don't get it. I think the reason is that I don't get the math behind it. Let us step back a bit. The following integration can be approximated using the rectangular method as

$$\int_a^bf(x)\,dx = \lim_{N\to \infty}\sum_{i=1}^{N-1}f(x_i)(x_{i+1}-x_i)$$

where $N$ is the number of points, $x_1 = a$, and $x_N = b$.

Now let's take the simple case of 2 random variables in our case. We want to evaluate
$$\text{Pr}\left[U_1+U_2\leq \gamma\right]=\text{E}_{U_1}\left\{\text{Pr}\left[U_2\leq \gamma-U_1/ U_1=u_1\right]\right\}\\=\int_{u_1=0}^{\gamma}\text{Pr}\left[U_2\leq \gamma-u_1\right]f_{U_1}(u_1)\,du_1\\=\int_{u_1=0}^{\gamma}\underbrace{F_{U_2}(\gamma-u_1)f_{U_1}(u_1)}_{G_1(u_1)}\,du_1=\lim_{N\to \infty}\sum_{i=1}^{N-1}G_1(u_{1,i})(u_{1, i+1}-x_{1, i})$$

where $u_{1,1}=0$, and $u_{1,N}=\gamma$.

For the case of 3 random variables, we want to evaluate
$$\text{Pr}\left[U_1+U_2+U_3\leq \gamma\right]=\text{E}_{U_1}\left\{\text{Pr}\left[U_2+U_3\leq \gamma-U_1/ U_1=u_1\right]\right\}\\=\int_{u_1=0}^{\gamma}\text{Pr}\left[U_2+U_3\leq \gamma-u_1\right]f_{U_1}(u_1)\,du_1=\\=\int_{u_1=0}^{\gamma}\text{E}_{U_2}\left\{\text{Pr}\left[U_3\leq \gamma-u_1-u_2/U_2=u_2\right]\right\}f_{U_1}(u_1)\,du_1\\=\int_{u_1=0}^{\gamma}\left[\int_{u_2=0}^{\gamma-u_1}\text{Pr}\left[U_3\leq\gamma-u_1-u_2\right]f_{U_2}(u_2)\,du_2\right]f_{U_1}(u_1)\,du_1= \int_{u_1=0}^{\gamma}\left[\int_{u_2=0}^{\gamma-u_1}\underbrace{{F_{U_3}(\gamma-u_1-u_2)f_{U_2}(u_2)\,du_2}}_{G_2(u_1+u_2)}\right]f_{U_1}(u_1)\,du_1$$

I get lost here on how to approximate the nested integrals as summations. That is why I don't get the logic behind the code.

EDIT (to give it a try): for a given $u_1$, the inner integral can be approximated as
$$\int_{u_2=0}^{\gamma-u_1}G_2(u_1+u_2)\,du_2=\lim_{N\to \infty}\sum_{i=1}^{N-1}G_2(u_{1}+u_{2,i})(u_{2, i+1}-u_{2, i})$$

where $u_{2,1} = 0$, and $u_{2,1}=\gamma-u_1$. But $u_1$ changes from $0$ to $\gamma$ as $u_{1,1}=0,\,u_{1,2},\ldots,\,u_{1,N}=\gamma$. So, at each point of $u_{1,j}$ for $j=1,\,2,\,\ldots,\,N$, we have

$$\int_{u_2=0}^{\gamma-u_1}G_2(u_{1,j}+u_2)\,du_2=\lim_{N\to \infty}\sum_{i=1}^{N-1}G_2(u_{1,j}+u^{(j)}_{2,i})(u^{(j)}_{2, i+1}-u^{(j)}_{2, i})$$

where $u^{(j)}_{2, 1} = 0$, and $u^{(j)}_{2, N} = \gamma-u_{1,j}$. Right? So, now to compute the whole integral as a summation approximation, we can write

$$\int_{u_1=0}^{\gamma}\left[\int_{u_2=0}^{\gamma-u_1}G_2(u_1+u_2)\,du_2\right]f_{U_1}(u_1)\,du_1\\=\lim_{N\to \infty}\sum_{j=1}^{N-1}\left[\sum_{i=1}^{N-1}G_2(u_{1,j}+u^{(j)}_{2,i})(u^{(j)}_{2, i+1}-u^{(j)}_{2, i})\right]f_{U_1}(u_{1,j})(u_{1,j+1}-u_{1,j})$$

It this correct?

Last edited:
andrewkirk
Science Advisor
Homework Helper
Gold Member
@EngWiPy : I spotted some silly errors in my post. The last argument to the function G, the vector, should not be there. It is not used in the Python code and should not be used in the function definition. We do need to pass a vector along, and use it, when the central integrand is a joint density function of all components of the vector and the components are interdependent. But in this case the components are independent, so we can deal with them one at a time as we emerge outwards through the nesting levels.

Also, I forgot to put the recursive call into the function in the code so, as written, it never dived down through the nesting levels.

I apologise for the errors, which may have made things harder than necessary to follow.
I've corrected my post above - both the code and the equations - and I have also made a simplification that reduces to two the number of arguments to the recursive function, which we can do by reversing the order of nesting of the integrals (which has no mathematical effect).

Regarding your post:

Your approach has the right general idea, but differs from (the corrected version of) what I wrote in the following ways:

1. I have a single function G (labelled 'integ' in the code), whereas you have a series of functions $G_1,G_2,....$. The differences between different Gs in your case is covered in mine by the 2nd argument to my (corrected) G function. When the 2nd argument to my function is $K-m$, it corresponds to your $G_m$, but with the following difference:

2. My G evaluates as an integral that gives a value of the CDF $F_m$, except in the terminal case where $m=K+1$ where it evaluates to 1. Your $G_m$ is the whole integrand, which is the CDF integral multiplied by a value of the PDF $f_m$.
Perhaps with my corrections and simplification it will be clearer how the recursive process works.

Last edited:
Thanks for the corrections, and pointing out the differences. Having a single function $G$ with two arguments, makes it simpler later in the code when the function calls itself recursively. Could you write the mathematical definition of your $G$ function?

I think it is more convenient for mathematical notations to average over $U_K, U_{K-1},\,\ldots,U_2$ instead of averaging over $U_1, U_2, \,\ldots,\,U_{K-1}$. In this case we can write $G(\gamma, K) = \text{Pr}\left\{\sum_{k=1}^KU_k\leq \gamma\right\}$ as

$$G(h, k) = \int_{u_k=0}^{h}G(h-u_k, k-1)f_{U_k}(u_k)\,du_k$$

for $k=2,\,3,\,\ldots,\,K$, and $G(h, 1)$ is the CDF of a single variable evaluated at $h$ (I don't see this value returned in your code. Maybe you are doing it differently?).

I will try to write the code in another way using Python. For simplicity, let us assume that the random variables $\{U\}_{k=1}^K$ are i.i.d. exponential random variables with parameter 1. So, the PDF and CDF of each random variables are defined as $f_U(u)=e^{-u}$ and $F_U(u)=1-e^{-u}$, respectively.

Now we can write the code for this as following:

Python:
import math
#Define the CDF of U
def F(u):
return 1-math,exp(-u)
#Define the PDF of U
def f(u):
return math.exp(-u)
granularity = 1000
def G(h, k):
if k == 1:
res = F(h)
else:
du = h/granularity
res, u = 0, 0
while u < h:
res += G(h-u, k-1)*f(u)*du
u += du
return res

#Find the total probability
prob = G(gamma, K)
print(prob)
Is this correct?

Using Monte-Carlo simulations, the value of $G(1,4)$ calculated as

Python:
count = 0
N = 1000000
for n in range(0, N):
if np.sum(np.random.exponential(1, (1, 4))) <= 1:
count += 1
print(count/N)
which gives 0.017 as a result. The numerical evaluation of the same function for a granularity of 300 gives 0.019 (granularity of 1000 took too long). I tried for other values $G(0.5, 5)$, and Monte-Carlo simulations gives 0.00018, while the numerical evaluation with granularity 100 gives 0.00016 (as $K$ increases, increasing the granularity increases the time for the numerical integration to be evaluated. So, I had to reduce it for $K=5$). I think this seems to work, but the question is: Is this a good approximation?

Last edited:
andrewkirk
Science Advisor
Homework Helper
Gold Member
Could you write the mathematical definition of your $G$ function?
The mathematical definition of G is in the corrected post #10, immediately above the code.

I am away from home this weekend with only limited access to a computer. But I did manage to do (on paper) an analytic evaluation of the nested integral in the case where U is exponential with PDF $f_U(u) = e^{-u}$. We can show by induction that, if the nesting depth is $n$, the value of the $m$ innermost levels for $m\le n$ is:

$$1 - e^{-(\gamma - s_m)}\sum_{k=0}^{m-1} \frac{(\gamma - s_m)^k}{k!}$$

where $s_m\triangleq \sum_{j=1}^{n-m} u_j$.

Hence the value of the full nested set of $n$ integrals is:
$$1 - e^{-\gamma}\sum_{k=0}^{n-1} \frac{\gamma ^k}{k!}$$

That value can be evaluated directly, without the need for any computer programming.

Note that the summation part is the first $n$ terms of the Taylor series for $e^{\gamma}$ so that, for a fixed $\gamma$, the value asymptotically approaches zero as $n\to\infty$, which is what we would expect.

The exponential assumption is just to test if the numerical approach works. In my case, as we were discussing in the other thread, $U_k=X_k/Y_k$, where $X_k$, and $Y_k$ are exponential random variables with parameter 1. The summation of $K$ i.i.d. exponential random variables is $\chi^2$-random variable with $K$ degrees of freedom.

OK, now I put the code together for the actual PDF and CDF of the random variable, and did the numerical and Monte-Carlo simulation for a range of values of $\gamma$. The code is as follows:

Python:
#************************** Import necessary libraries******************************************
import numpy as np
import math
import matplotlib.pyplot as plt
%matplotlib inline

#******************************Set the constant scalars and vectors***************************
granularity = 700
K = 3
beta_dB = 3
beta = 10**(beta_dB/10)
EF = 100

gammadB = np.arange(-5,11)
gammaVec = 10**(gammadB/10)

N = 100000 # For simulation

#******************************Define functions to be used***************************************
#Define the CDF of U
def F(u):
return 1-1/(u+1)
#Define the PDF of U
def f(u):
return 1/((1+u))**2

def G(h, k):
#print(f'h is {h}, and k is {k}')
if k == 1:
res = F(h)
else:
du = h/granularity
res =  0
u = 0
while u < h:
res += G(h-u, k-1)*f(u)*du
u += du
return res
#*******************Find the numerical and simulation values of the integral******************
Pout = []
PoutSim = []
for gamma in gammaVec:
gammath = (gamma/beta)*(EF-beta*(K-1))
# Numerical
Pout.append(1-G(gamma, K))
# Simulation
count = 0
for n in range(0, N):
if np.sum(np.random.exponential(1, (1, K))/np.random.exponential(1, (1, K))) <= gamma:
count += 1
PoutSim.append(1-count/N)

#**************************Visualization******************

fig, ax = plt.subplots(figsize = (10, 8))
major_ticks = np.arange(0.1, 1.1, 0.1)
minor_ticks = np.arange(0.1, 1.1, 0.1)

ax.set_yticks(major_ticks)
ax.set_yticks(minor_ticks, minor=True)

ax.set_yscale('log')

ax.grid(which='minor', alpha=0.2, linestyle='--', linewidth=1, color='k')
ax.grid(which='major', alpha=0.5, linestyle='-', linewidth=1, color='k')

ax.plot(gammadB, Pout, 'r-*', label = 'Numerical')
ax.plot(gammadB, PoutSim, 'k-o', label = 'Sim')

ax.legend()
ax.grid(True)

plt.show()

fig.savefig("test.png", bbox_inches='tight')
The resulted figure is attached as test.png. It shows that both yield very close curves. Which is good.

What I actually want to do, is to replace gamma, by gammath. The changes will be

Python:
 ...
Pout.append(1-G(gammath, K))
...
if np.sum(np.random.exponential(1, (1, K))/np.random.exponential(1, (1, K))) <= gammath:
....
This resulted in test1.png, which is also attached. Here the curves are very different. Why? The simulation curve makes more sense than the numerical evaluation. This implies that there is something wrong in the numerical implementation of the integral. What is wrong in the numerical implementation?

#### Attachments

• 7.9 KB Views: 280
• 8.5 KB Views: 249
andrewkirk
Science Advisor
Homework Helper
Gold Member
@EngWiPy You are using a constant relative granularity as your integrating range increases. That means that the size of the integrating increment du gets bigger as the outer limit gamma or gammath gets bigger. Since the overstatement of the simple rectangular approach is proportional to the size of the increments, the integrating error gets larger when you get to very large ranges, as is the case with gammath, which goes to over 600.

You could fix this by setting the granularity to be a fixed absolute, rather than relative, amount. eg set du globally to be 0.01.

@EngWiPy You are using a constant relative granularity as your integrating range increases. That means that the size of the integrating increment du gets bigger as the outer limit gamma or gammath gets bigger. Since the overstatement of the simple rectangular approach is proportional to the size of the increments, the integrating error gets larger when you get to very large ranges, as is the case with gammath, which goes to over 600.

You could fix this by setting the granularity to be a fixed absolute, rather than relative, amount. eg set du globally to be 0.01.
Thanks. I am testing it. 0.01 wasn't small enough, so I decreased it to 0.001, and now it is taking too long. I will report back when I have some results.

andrewkirk
Science Advisor
Homework Helper
Gold Member
Thanks. I am testing it. 0.01 wasn't small enough, so I decreased it to 0.001, and now it is taking too long. I will report back when I have some results.
Re run-time, there are two simple strategies that, taken together, should greatly reduce it:

1. Use trapezoidal rather than rectangular incremental areas. To do that you'll need to insert another variable to hold the previous value in each iteration of the 'while u<h' loop, which you'll use together with the new value to calc the incremental trapezium's area.

2. Based on the PDF, generate a set of integrating increments that starts at say 0.001 and increases as the PDF becomes flatter. Then for each u in the central loop, select the du that corresponds to the current value of u. That way you'll be using fine increments in the early part of the curve where approximations are least accurate, and coarse increments in the right tail. eg you could set du to the greater of 0.001 and $0.0005\,e^{u/5}$. You could change the coefficient and the 5 in the exponent as needed to balance between runtime and accuracy.

Where to insert No. 2 in the code?

What I did is that I fixed du for each gammath as gammath/granularity. To do this, I added a new argument to the function G(h, K, du), such that du remains constant for each call of the funcion from outside the function. I got a little faster, but not by that much.

I ran the same code on MATLAB, and it was MUCH faster (I ran the code the whole night on Python, and it wasn't finished in the morning. When I got to the office, and use MATLAB for the same parameter values, I got the result in less 5 minutes!!).

Python has the option of parallel programming, but not sure if we can use it here, and how.

PS: OK, on my office computer, I have i7 with 8GB of RAM, while at home on my laptop I have i5 with 6GB of RAM. So, maybe the difference lies there. Unfortunately, I do't have Python on my office computer to test this hypothesis.

Last edited:
This implementation in MATLAB runs in less than 30 seconds on my laptop. Interesting problem. If I understand your math correctly, there may be a way to call the built-in integration routines for a comparison. If I get a chance, I will look at this again.

Code:
function [] = nestint()
% www.physicsforums.com/threads/nested-integrals.951242/#post-6030169
K = 3;
beta_dB = 3;
% beta = 10^(beta_dB/10);
% EF = 100;
gammadB = -5:10;
L = length(gammadB);
gammaVec = 10.^(gammadB/10);
N = 100000;
Pout = zeros(1,L);
PoutSim = zeros(1,L);

for gamma = 1:L
%     gammath = (gammaVec(gamma)/beta)*(EF-beta*(K-1));
Pout(gamma) = (1-G(gammaVec(gamma), K));
count = 0;

for n =0:(N-1)
if sum(exprnd(1,1,K)./exprnd(1,1,K)) <= gammaVec(gamma)
count = count +  1;
end
end

PoutSim(gamma) = 1-count/N;
end

semilogy(gammadB, Pout)
hold on
semilogy(gammadB, PoutSim)

function res = G(h, k)
du = .001;

if k == 1
res = F(h);
else
res =  0;
u = 0;
while u < h
res = res +  G(h-u, k-1)*f(u)*du;
u = u + du;
end
end

function x = F(u)
x = 1-1./(u+1);

function x = f(u)
x = 1./((1+u)).^2;

andrewkirk
Science Advisor
Homework Helper
Gold Member
I've done an R implementation (Python turned out to be too painful, with my relative inexperience in that language) below. I used two techniques to improve accuracy:
1. Setting the integration increment to be inversely proportional to the integrand, with a maximum increment size
2. Using cubic spline interpolation to reduce errors from nonlinearity
It starts to misbehave when the gamma integration limit is very large, I think because the integrand gets so close to zero, and it gets in a mess if it goes below zero, which can easily happen with a cubic spline interpolation. Putting a floor of zero creates big errors - I don't know why yet, cos I'd expect the errors from that to be negligible. So instead I got it to switch to trapezoidal integration if any of the cubic splines dip below zero.

Attached are the code and three plots. The most accurate used cubic splines with a maximum increment of 3, and took 3.5 mins to run. The next was cubic with max increment of 4, and took 1.5 mins to run. Last is one that used trapezoidal integration, which took 14 secs to run. You can see it loses accuracy for gamma (integration limit) beyond 100, but the probability is already below 3% at that stage. The log scale exaggerates the errors.

I think the key point is that for deeply nested integrals like this, Monte Carlo simulation is actually more accurate for a given runtime than numerical integration. The MC sim runs in a couple of seconds and is much smoother (which I take to be a sign of accuracy.
Code:
dev.new()
time.save <- Sys.time()
#****************************^Set the constant scalars and vectors************************^*
spline.on <- FALSE # TRUE
diag.on <- FALSE
do.big.ones <-  FALSE # TRUE #
granularity <- 700 # 700
max.increment <- 10 # set to 4 for high accuracy, with run taking 3.5 mins. 3 gives reasonable and runs in 1.5 mins
K <- 3
beta0_dB <- 3
beta0 <- 10^(beta0_dB/10)
EF <- 100

set.seed(123)

gamma.dB <- (-5):11
gamma.Vec <- 10^(gamma.dB/10)

# du <- max(gamma.Vec) / granularity
dV <- 1 / granularity

N <- 100 * 1000 # For simulation

#****************************^Define functions to be used************************************^*
print.diag <- function(t){ # print the argument t only if global flag diag.on is TRUE
if (diag.on) print (t)
}

# Define the CDF of U
F <- function(u){
1-1/(u+1)
}
# Define the PDF of U
f <- function(u){
1/((1+u))^2
}
# curve(f,from = 0, to = 10)

$a function for integration on a piece of the spline piecewise_int <- function (hi, yi, bi, ci, di) { yi * hi + bi * hi ^ 2 / 2 + ci * hi ^ 3 / 3 + di * hi ^ 4 / 4 } integ_fn <- function(x, y){ # do an integral of the cubic spline fitted to the points x, y sf <- splinefun(x, y, 'natural') # information for integration n <- length(x) int_info <- with(environment(sf)z, list(h = diff(x), y = y[-n], b = b[-n], c = c[-n], d = d[-n]) ) res <- do.call(piecewise_int, int_info) if (any(res < 0)){ # do trapezoidal integration instead sum(diff(x) * (y[-length(y)] + 0.5 * diff(y))) } else$ cubic spline integration on all pieces
sum(pmax(0,res))
}

big.up.gamma <- function(g) (g / beta0) * (EF - beta0 * (K-1))

G <- function(h, h0, k, diag){ # prob that the sum of k iid vbls lies between h0 and h
u.all <- I.all <- NULL
if (k == 1){
res <- F(h) - F(h0)
} else { # not base case. Set up to recurse
res <- res0 <- curv <- 0
u3 <- I3 <- rep(0,3)
while (u3[2] < h) {
I3[3] <- G(h - u3[3], max(0, h0 - u3[3]), k-1, diag1) * f(u3[3])
if (u3[3] > 0){
I.av <- 0.5 * sum(I3[2:3])
res0 <- res0 + du * I.av
if (u3[2] > 0){
du <- min(max.increment, dV / I.av) # else leave du at initial low value

curv <- ((I3[3] - I3[2]) / (u3[3] - u3[2]) - (I3[2] - I3[1]) / (u3[2] - u3[1])) / (0.5 * (u3[3] - u3[1]))
curvhist <<- c(curvhist, curv)

# currently disable the curvature adjustment, as it slows down run and doesn't change output
# du <- du / max(1, abs(curv) / 0.5)

} # end inner if
} else # end outer if, except for else on next line
du <- 1 / granularity # first point in integration
u.all <- c(u.all, u3[3])
I.all <- c(I.all, I3[3])
u3 <- c(u3[-1], min(h, u3[3] + du))
I3 <- c(I3[-1], -1)
} # end while
x <- c(u.all, I.all)
n <- length(u.all)
if (n > 0){
if (spline.on){
res <- integ_fn(u.all, I.all)
#res <- max(0, res)
} else {
res <- sum(diff(u.all) * (I.all[-n] + 0.5 * diff(I.all)))
}
} else res <- 0
} # end outermost if

res
}
#****************^*Find the numerical and simulation values of the integral****************^
curvhist <- NULL
Pout <- 1
PoutSim <- NULL

call.count <- 0

exp.inv.cdf <- function(x){ # inverse CDF of exponential with parameter 1
-log(1 - x)
}

gamma.Vec <- c(gamma.Vec, big.up.gamma(gamma.Vec))
# Do Numerical integration version
gamma.index <- 1
gamm.last <- 0
for (gamma0 in gamma.Vec){
gamm <- if (do.big.ones) big.up.gamma(gamma0) else gamma0
print(paste('gamma.index is', gamma.index,'-------------------------'))
x <- tail(Pout, 1) - G(gamm, gamm.last, K, c(gamm, rep(-1, K - 1)))
Pout <- c(Pout, x)

if (x < 0){
x=x
}

gamm.last <- gamm
gamma.index <- gamma.index + 1
}
Pout <- Pout[-1]

# now do simulation version
simx <- exp.inv.cdf(runif(K * N)) / exp.inv.cdf(runif(K * N))
dim(simx) <- c(N, K)
PoutSim <- outer(apply(simx, 1, sum), gamma.Vec, "<=")
PoutSim <- 1 - apply(PoutSim, 2, sum) / N

#************************^Visualization****************^
time.end <- Sys.time()
time.gone <- difftime(time.end, time.save, units = "mins")
print(paste('time taken is', time.gone, "minutes"))

plotcurves <- function(gamma.Vec.x, Pout.x, PoutSim.x){
Pout.x1 <- pmax(Pout.x, min(Pout.x[Pout.x > 0])) # set min of the min positive Pout value
lims <- rep(max(c(Pout.x1, PoutSim.x)), 2)
lims[1] <- min(c(Pout.x1, PoutSim.x))
plot(gamma.Vec.x, Pout.x1, type='b',log="y", ylim = lims)
lines(gamma.Vec.x, PoutSim.x, type='l')
grid()
}

dev.new()
plotcurves(gamma.Vec, Pout, PoutSim)

#### Attachments

• 24.3 KB Views: 369
• 23.2 KB Views: 389
• 25.9 KB Views: 405
Last edited:
This implementation in MATLAB runs in less than 30 seconds on my laptop. Interesting problem. If I understand your math correctly, there may be a way to call the built-in integration routines for a comparison. If I get a chance, I will look at this again.

Code:
function [] = nestint()
% www.physicsforums.com/threads/nested-integrals.951242/#post-6030169
K = 3;
beta_dB = 3;
% beta = 10^(beta_dB/10);
% EF = 100;
gammadB = -5:10;
L = length(gammadB);
gammaVec = 10.^(gammadB/10);
N = 100000;
Pout = zeros(1,L);
PoutSim = zeros(1,L);

for gamma = 1:L
%     gammath = (gammaVec(gamma)/beta)*(EF-beta*(K-1));
Pout(gamma) = (1-G(gammaVec(gamma), K));
count = 0;

for n =0:(N-1)
if sum(exprnd(1,1,K)./exprnd(1,1,K)) <= gammaVec(gamma)
count = count +  1;
end
end

PoutSim(gamma) = 1-count/N;
end

semilogy(gammadB, Pout)
hold on
semilogy(gammadB, PoutSim)

function res = G(h, k)
du = .001;

if k == 1
res = F(h);
else
res =  0;
u = 0;
while u < h
res = res +  G(h-u, k-1)*f(u)*du;
u = u + du;
end
end

function x = F(u)
x = 1-1./(u+1);

function x = f(u)
x = 1./((1+u)).^2;
May I ask what are your computer specifications? Because you don't seem to change anything in the code. I am not sure if Python is slower than MATLAB, or it is a machine speed issue, because I have Python on my laptop (i5 with 6 GB of RAM) but not MATLAB, and on my office computer (i7 and 8 GB of RAM) I have MATLAB but not Python (I am an administrator and cannot download it without asking it to be installed).

I've done an R implementation (Python turned out to be too painful, with my relative inexperience in that language) below. I used two techniques to improve accuracy:
1. Setting the integration increment to be inversely proportional to the integrand, with a maximum increment size
2. Using cubic spline interpolation to reduce errors from nonlinearity
It starts to misbehave when the gamma integration limit is very large, I think because the integrand gets so close to zero, and it gets in a mess if it goes below zero, which can easily happen with a cubic spline interpolation. Putting a floor of zero creates big errors - I don't know why yet, cos I'd expect the errors from that to be negligible. So instead I got it to switch to trapezoidal integration if any of the cubic splines dip below zero.

Attached are the code and three plots. The most accurate used cubic splines with a maximum increment of 3, and took 3.5 mins to run. The next was cubic with max increment of 4, and took 1.5 mins to run. Last is one that used trapezoidal integration, which took 14 secs to run. You can see it loses accuracy for gamma (integration limit) beyond 100, but the probability is already below 3% at that stage. The log scale exaggerates the errors.

I think the key point is that for deeply nested integrals like this, Monte Carlo simulation is actually more accurate for a given runtime than numerical integration. The MC sim runs in a couple of seconds and is much smoother (which I take to be a sign of accuracy.
View attachment 228436 View attachment 228437 View attachment 228438

Code:
dev.new()
time.save <- Sys.time()
#****************************^Set the constant scalars and vectors************************^*
spline.on <- FALSE # TRUE
diag.on <- FALSE
do.big.ones <-  FALSE # TRUE #
granularity <- 700 # 700
max.increment <- 10 # set to 4 for high accuracy, with run taking 3.5 mins. 3 gives reasonable and runs in 1.5 mins
K <- 3
beta0_dB <- 3
beta0 <- 10^(beta0_dB/10)
EF <- 100

set.seed(123)

gamma.dB <- (-5):11
gamma.Vec <- 10^(gamma.dB/10)

# du <- max(gamma.Vec) / granularity
dV <- 1 / granularity

N <- 100 * 1000 # For simulation

#****************************^Define functions to be used************************************^*
print.diag <- function(t){ # print the argument t only if global flag diag.on is TRUE
if (diag.on) print (t)
}

# Define the CDF of U
F <- function(u){
1-1/(u+1)
}
# Define the PDF of U
f <- function(u){
1/((1+u))^2
}
# curve(f,from = 0, to = 10)

$a function for integration on a piece of the spline piecewise_int <- function (hi, yi, bi, ci, di) { yi * hi + bi * hi ^ 2 / 2 + ci * hi ^ 3 / 3 + di * hi ^ 4 / 4 } integ_fn <- function(x, y){ # do an integral of the cubic spline fitted to the points x, y sf <- splinefun(x, y, 'natural') # information for integration n <- length(x) int_info <- with(environment(sf)z, list(h = diff(x), y = y[-n], b = b[-n], c = c[-n], d = d[-n]) ) res <- do.call(piecewise_int, int_info) if (any(res < 0)){ # do trapezoidal integration instead sum(diff(x) * (y[-length(y)] + 0.5 * diff(y))) } else$ cubic spline integration on all pieces
sum(pmax(0,res))
}

big.up.gamma <- function(g) (g / beta0) * (EF - beta0 * (K-1))

G <- function(h, h0, k, diag){ # prob that the sum of k iid vbls lies between h0 and h
u.all <- I.all <- NULL
if (k == 1){
res <- F(h) - F(h0)
} else { # not base case. Set up to recurse
res <- res0 <- curv <- 0
u3 <- I3 <- rep(0,3)
while (u3[2] < h) {
I3[3] <- G(h - u3[3], max(0, h0 - u3[3]), k-1, diag1) * f(u3[3])
if (u3[3] > 0){
I.av <- 0.5 * sum(I3[2:3])
res0 <- res0 + du * I.av
if (u3[2] > 0){
du <- min(max.increment, dV / I.av) # else leave du at initial low value

curv <- ((I3[3] - I3[2]) / (u3[3] - u3[2]) - (I3[2] - I3[1]) / (u3[2] - u3[1])) / (0.5 * (u3[3] - u3[1]))
curvhist <<- c(curvhist, curv)

# currently disable the curvature adjustment, as it slows down run and doesn't change output
# du <- du / max(1, abs(curv) / 0.5)

} # end inner if
} else # end outer if, except for else on next line
du <- 1 / granularity # first point in integration
u.all <- c(u.all, u3[3])
I.all <- c(I.all, I3[3])
u3 <- c(u3[-1], min(h, u3[3] + du))
I3 <- c(I3[-1], -1)
} # end while
x <- c(u.all, I.all)
n <- length(u.all)
if (n > 0){
if (spline.on){
res <- integ_fn(u.all, I.all)
#res <- max(0, res)
} else {
res <- sum(diff(u.all) * (I.all[-n] + 0.5 * diff(I.all)))
}
} else res <- 0
} # end outermost if

res
}
#****************^*Find the numerical and simulation values of the integral****************^
curvhist <- NULL
Pout <- 1
PoutSim <- NULL

call.count <- 0

exp.inv.cdf <- function(x){ # inverse CDF of exponential with parameter 1
-log(1 - x)
}

gamma.Vec <- c(gamma.Vec, big.up.gamma(gamma.Vec))
# Do Numerical integration version
gamma.index <- 1
gamm.last <- 0
for (gamma0 in gamma.Vec){
gamm <- if (do.big.ones) big.up.gamma(gamma0) else gamma0
print(paste('gamma.index is', gamma.index,'-------------------------'))
x <- tail(Pout, 1) - G(gamm, gamm.last, K, c(gamm, rep(-1, K - 1)))
Pout <- c(Pout, x)

if (x < 0){
x=x
}

gamm.last <- gamm
gamma.index <- gamma.index + 1
}
Pout <- Pout[-1]

# now do simulation version
simx <- exp.inv.cdf(runif(K * N)) / exp.inv.cdf(runif(K * N))
dim(simx) <- c(N, K)
PoutSim <- outer(apply(simx, 1, sum), gamma.Vec, "<=")
PoutSim <- 1 - apply(PoutSim, 2, sum) / N

#************************^Visualization****************^
time.end <- Sys.time()
time.gone <- difftime(time.end, time.save, units = "mins")
print(paste('time taken is', time.gone, "minutes"))

plotcurves <- function(gamma.Vec.x, Pout.x, PoutSim.x){
Pout.x1 <- pmax(Pout.x, min(Pout.x[Pout.x > 0])) # set min of the min positive Pout value
lims <- rep(max(c(Pout.x1, PoutSim.x)), 2)
lims[1] <- min(c(Pout.x1, PoutSim.x))
plot(gamma.Vec.x, Pout.x1, type='b',log="y", ylim = lims)
lines(gamma.Vec.x, PoutSim.x, type='l')
grid()
}

dev.new()
plotcurves(gamma.Vec, Pout, PoutSim)
I am familiar with MATLAB and Python, but not R. One trick is how to control du, I think. What do you mean you set the increment to be inversely proportional to the integrad with a maximum increment size? Could you write that mathematically, please?

Monte-Carlo simulations run fast and accurate, but can be used as a verification to the numerical simulations when you write an academic paper. Usually, it cannot be used alone as an evaluation method.

I had to consider C++ as another language to implement the integration. For the same code, I got 30x speed factor over that written in Python, for a constant increment size. There is the option of multi-threading because the values in the FOR loop are independent, and could be implemented in parallel. I tried this tentatively, and it took 1min 9sec mins for a constant increment size du=0.001, while the original code with the same resolution in C++ took 2min 18sec, which isn't that a significant speed up factor!!

Last edited: