MATLAB Evaluating Nested Integral in MATLAB for General K

[andrewkirk]

Where did I make a mistake?

Best not to use u inside the G function to denote the limit for the previous integral, because u is used for the integration variable. Your call of the function G correctly passes that previous limit as the 3rd argument, but when you pass it on as you recurse down the stack, we need to pass on not u but that previous limit minus u. The function should read as follows:

def G(h, k, h_prev):
    #print(f'h is {h}, and k is {k}')
    if k == 1:
        res = F(h)-F(h_prev)
    else:
        #du = h/granularity
        res =  0
        u = 0
        while u < h:
            res += G(h-u, k-1, max(0,h_prev-u))*f(u)*du
            u += du
    return res

Note:

• u=0 has been uncommented, as 0 is the correct starting point for the integrating variable in all cases except the base case
• in the base case k==1 we deduct the CDF at h_prev.
  
• the use of max in the recursing call to prevent the lower limit going below zero.

 

[EngWiPy]

Yes, you are right about all three points, but interestingly, the evaluation time has increased for this approach about 1.3 times relative to the original approach! Does this make sense?

 

[andrewkirk]

Yes, you are right about all three points, but interestingly, the evaluation time has increased for this approach about 1.3 times relative to the original approach! Does this make sense?

I think that can happen, and I think I realize why. Fortunately, thinking over that has given me a ripper of an idea! I've got a new program that runs in under 10 seconds with very high accuracy. It uses iteration rather than recursion, yet is completely general about depth of nesting. The trick is to calculate the CDFs from the bottom up, thereby avoiding all the repetition that goes on in existing versions. Let $Z_k$ be the sum of $k$ iid versions of RV $X$ that has PDF $f_U$ and CDF $F_U$. Then the PDF $f_{Z_1}$ and CDF $F_{Z_1}$ are the same as $f_U$ and $F_U$. Then for each level $k$ we can calculate the CDF of $Z_k$ in terms of $f_U$ and $F_{Z_{k-1}}$, as the following convolution:

$$F_{Z_k}(u) = \int_0^\infty f_U(v) F_{Z_{k-1}}(u-v)\,dv$$

Then, by working from the bottom up we can calculate the CDFs for sums of increasing numbers of $U$s.

We specifying a maximum allowable increment of probability weight per u increment and use that to calculate the required u granularity. We use that granularity to set a partition of the interval from 0 to the maximum (you were using about 606, and I've kept that), then we start building up from $Z_1=U$ up to $Z_k$ for $k$ as large as we like. For each step the computational complexity is only $O(n^2)$, where $n$ is the number of segments used in the range of possible values. I think the complexity of the recursion approach is at least$O(n^K)$, which is why it takes so long.

EDITED the above para to change X to U, to be consistent with the OP.

I'll post the R code and two graphs - one log-scale and one not. You can see it's very accurate. The red line is from the sim and the black line from the numerical integration. Then I'll try to do a Python version for Python practice and post that later.

time.save <- Sys.time()
#*****************************Set the constant scalars and vectors**************************
max.dp <- 0.0001 # granularity will be set so that no interval has an incremental probability weight greater than this
numbins.sim <- 1000 # number of equal-frquency bins to use for graphing the simulation results
K <- 3 # number of iid variables being convoluted ( = nbr of levels of integral nesting)
N <- 100 * 1000 # Number of simulations
max.lim <- 606 # CDF will be calculated in interval from 0 to here
set.seed(123) # set seed for random number generator so that simulation results are reproducible
#*****************************Define functions to be used**************************************
# Define the CDF of U
F <- function(u){
      1-1/(u+1)
}
# Define the PDF of U
f <- function(u){
     1/((1+u))^2
}
#******************Find the numerical values of the integral*****************
# find smallest integration increment to satisfy max.dp constraint
x <- max(f(0), f(max.lim)) # start search loop with x > max.dp
y <- 0
while (x > max.dp){
  y <- y + 1
  z <- 2^-y
  du <- max.lim * z
  u <- du * 0:(2^y)  # divide up interval into 2^y equal segments
  x <- max(f(u) * 2^-y)
}
n <- 2^y + 1

convolution <- array(F(u), dim=c(n, K)) # initialise this array. It  will hold all the convolved CDFs

for (j in 2:K){
  # calc CDF for a convolution of j iid variables of PDF f and CDF F
  convolution[,j] <- 0
  for (i in 1:n){
    # calculate probability that sum of j iid variables is <= u[i]
    x <- du * f(u[i:1]) * convolution[1:i, j - 1]
    if (i %% 2 == 0){ # have extra segment to add to simpson integration
      y <- sum(x[1:2]) / 2
      x <- x[-1] # drop first element of integrand vector , so there are an odd number of points
    } else { # pure simpson integration, so make add-on zero
      y <- 0
    }
    m <- length(x) # this is the number of integration points for the Simpson calc, and it will be odd
    if (m == 1)
      y <- x[1] else
        y <- y + (1/3) * (2 * sum(x) - x[1] - x[m] + 2 * sum(x[2 * (1:((m - 1) / 2))]))
    convolution[i,j] <- min(1, y)
  }
}
Pout <- 1 - convolution[,3]
#******************Find the simulation values of the integral*****************
exp.inv.cdf <- function(x){ # inverse CDF of exponential with parameter 1
  -log(1 - x)
}
simx <- exp.inv.cdf(runif(K * N)) / exp.inv.cdf(runif(K * N))
dim(simx) <- c(N, K)
simx <- apply(simx, 1, sum)
simx <- simx[order(simx)]
u.sim <- simx[round((1:numbins.sim - 0.5) * N / numbins.sim)]
PoutSim <- 1 - (1:numbins.sim) / numbins.sim
#*************************Visualization*****************
time.end <- Sys.time()
time.gone <- difftime(time.end, time.save, units = "mins")
print(paste('time taken is', time.gone, "minutes"))

plotcurves <- function(x1, y1, x2, y2, log = FALSE){
  if (log) {
    y <- c(y1, y2)
    # y1.1 <- pmax(y1, ) # set min of the min positive Pout value
    # lims <- rep(max(c(y1.1, y2)), 2)
    # lims[1] <- min(c(y1.1, y2))
    y.min <- min(y[y > 0])
    y1 <- pmax(y.min, y1)
    y2 <- pmax(y.min, y2)
    lims <- c(y.min, 1)
    plot(x1, y1, type='l',log="y", ylim = lims)
  } else {
    lims <- 0:1
    plot(x1, y1, type='l', ylim = lims)
  }
  lines(x2, y2, type='l', col='red')
  grid()
}

png('plain.graph.png')
plotcurves(u, Pout, u.sim, PoutSim)
dev.off()
png('log.graph.png')
plotcurves(u, Pout, u.sim, PoutSim, log = T)
dev.off()

 

[EngWiPy]

Very interesting. I am waiting to see your Python code, as I cannot understand the R code, and how you implemented the algorithm using loops instead of recursions.

How long did the original code with recursion and using the same numerical approach and parameters take? You mentioned that the cubic spline was the most accurate, and took 3min 30sec, but I am sure what numerical approach you are using here.

 

[andrewkirk]

I am waiting to see your Python code

Here we are. I got it working in the end:

#************************** Import necessary libraries******************************************
import numpy as np
import math
import matplotlib.pyplot as plt
import datetime
#************************** Set global constants ******************************************
max_dp = 0.0001# 001 # granularity will be set so that no interval has an incremental probability weight greater than this
numbins_sim = 1000# 00 # number of equal-frquency bins to use for graphing the simulation results
K = 3 # number of iid variables being convoluted ( = nbr of levels of integral nesting)
N = 100 * 1000 # Number of simulations
max_lim = 606 # CDF will be calculated in interval from 0 to here
np.random.seed(123) # set seed for random number generator so that simulation results are reproducible

time_start = datetime.datetime.now()
#*****************************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+u)**-2
#******************Find the numerical values of the integral*****************
# find smallest integration increment to satisfy max.dp constraint
x = max(f(0), f(max_lim)) # start search loop with x > max.dp
y = 0
while x > max_dp:
  y = y + 1
  n = 2**y + 1
  z = 2**-y
  du = max_lim * z
  u = [x * du for x in range(n)] 
  x = max([f(u1) for u1 in u]) * z

# now start the convolution process, working from one RV up to K RVs
convolution = [[F(u[i]) for i in range(n)]] # base case of one RV
# now work up from two to K RVs
for j in range(1,K):
  # calc CDF for a convolution of j iid variables of PDF f and CDF F
  convolution.append([0] * n)
  for i in range(n):
    # calculate probability that sum of j iid variables is <= u[i]
    x = [du * f(u[i-k]) * convolution[j-1][k] for k in range(i + 1)]
    if (i % 2 == 1): # have extra segment to add to simpson integration
      y = sum(x[:2]) / 2
      x = x[1:] # drop first element of integrand vector , so there are an odd number of points
    else:  # pure simpson integration, so make add-on zero
      y = 0    
    m = len(x) # this is the number of integration points for the Simpson calc, and it will be odd
    if m == 1:
      y = x[0]
    elif m == 0:
      y = 0
    else:
      y = y + (1/3) * (2 * sum(x) - x[0] - x[m-1] + 2 * sum([x[2 * k + 1] for
      k in range(round((m - 1) / 2))]))
    convolution[j][i] = min(1, y) 
Pout = [1 - convolution[2][k] for k in range(n)]
#******************Find the simulation values of the integral*****************
exp1 = np.random.exponential(1, [N, K])
exp2 = np.random.exponential(1, [N, K])
simx = [[exp1[i][j] / exp2[i][j] for j in range(K)] for i in range(N)]
sumsimx = [sum([simx[i][j] for j in range(K)]) for i in range(N)]
sumsimx = sorted(sumsimx) # sort the list
u_sim = [sumsimx[round(N * (j + 0.5) / numbins_sim)] for j in range(numbins_sim)]
PoutSim = [1 - j / numbins_sim for j in range(numbins_sim)]
#*************************Visualization*****************
time_end = datetime.datetime.now()
print("Time elapsed for run is:")
print(time_end - time_start)

fig, ax = plt.subplots(figsize = (10, 8))
ax.plot(u, Pout, 'b-')
ax.plot(u_sim, PoutSim, 'r--')
ax.axis([0, 607, 0, 1])
plt.show()
fig.savefig("test_ak.png", bbox_inches='tight')

It produces results that match the sim and numeric results to each other very closely, as shown in the graph below, which is plain scale, as I don't know how to do log scale in Python. The numeric curve is unbroken blue and the sim surve is dashed red. If you look closely, you can see that they overlay one another as closely as the naked eye can make out. On a log scale one might be able to see a bit of separation.

The Python version runs in 2 mins 22 seconds, compared to only 11 seconds for the R version. That may be just because I'm a clumsy Python programmer, having only just started learning it. Perhaps someone more experienced with Python than me could implement the algorithm using operations that are more efficient than those I have used.

 

[EngWiPy]

Thanks. I think you use too many list comprehensions, which are in effect FOR loops. I am not an expert in Python, but I think we should be able to use vectorization. The fact that R runs in 11 secs, while the Python version runs in 2 min and 22 secs indicates an inefficiency in the coding. Do you use vectorization in your R code instead of the for loops or list comprehension in Python? I will see what I can do.

 

[andrewkirk]

Do you use vectorization in your R cod

Yes the R code is vectorised wherever possible, which is in most of the code. I look forward to learning how to vectorise in Python!

 

[EngWiPy]

Yes the R code is vectorised wherever possible, which is in most of the code. I look forward to learning how to vectorise in Python!

That's why the R code is faster. I am not sure how to vectorize the code in Python myself, but I am researching it. Will update if I find anything. Thanks

 

[EngWiPy]

Let's us take these two lines

...
u = [x * du for x in range(n)]
x = max([f(u1) for u1 in u]) * z
...

These can be replaced by

import numpy as np
#Vectorize the PDF function to be applied to a vector
f_vec = np.vectorize(f)

u = x*np.arange(n)
x = max(f_vec(u))*z

I will try to spot where vectorization can be used in your code.

 

[EngWiPy]

After some thought, I don't think I get your algorithm. I understand you go from $Z_1$ to $Z_1+Z_2$ up to $Z_1+Z_2+\cdots+Z_K$, but not sure how $\text{Pr}\left[Z_1+Z_2+Z_3\leq z\right]$ could be computed from $\text{Pr}\left[Z_1+Z_2\leq z\right]$. The integration you wrote is the same integration we are trying to evaluate recursively. How do you go forward instead of backward? Let take a simple example where $K=3$, could you list the procedures of how to evaluate the 2-dimensional integral forwardly?

 

[andrewkirk]

integration you wrote is the same integration we are trying to evaluate recursively.

here's the integral again:
$$
F_{Z_k}(u) = \int_0^\infty f_U(v) F_{Z_{k-1}}(u-v)\,dv
$$
In the recursive version we start at the top, with $k=3$ and in the formula we use a recursive call to evaluate $F_{Z_{k-1}}(u-v)$. Recall that $F_{Z_k}$ is the CDF of $Z_k$, which is the sum of k iid copies of the base random variable $U$ whose PDF is $f_U$.

The light-bulb moment was to realize that, if we start at the bottom (k=1) instead of the top (k=3), at each level k, we can save a vector that gives $F_{Z_k}(u)$ for $u$ covering the support interval at small increments. Then when we increment $k$ by 1 and apply the above formula to calculate the next level (one higher than the previous one), we don't need to recurse down to calculate $F_{Z_{k-1}}(u-v)$. Rather we just look up the value in the vector we created at the previous step.

If you look at the code you'll see that it builds up a K x n matrix called convolution, for which the $k$th column is a vector that saves $F_{Z_k}(u)$. We look up values in the $(k-1)$th column in the course of calculating each entry in the $k$th column. We get a vector of integrand values and we then use Simpson's rule to estimate the integral that gives us $F_{Z_k}(u)$. There's a bit of fiddling about to cope with the integration in cases where we have an even number of points, or only one point (Simpson requires an odd number greater than 2).

The vector of points to be integrated (multiplied by du) is set up as x in this line:

 x = [du * f(u[i-k]) * convolution[j-1][k] for k in range(i + 1)]

note how it looks up the previous column (j-1) of convolution to get values of $F_{j-1}$, rather than recursing down to calculate them.

EDITED: changed X to U for consistency with OP

 

[EngWiPy]

OK, right. Yesterday I understood the math, but somehow I forgot it today. You are right. We view $Z_k$ as $Z_{k-1}+U_k$ for $k=2,\,3,\,\ldots, K$, and $Z_1=U_1$. That is a neat idea. Now I need to get my head around the implementation. How do you deal with the $\infty$ limit in the convolution? I think the upper limit should be $u$, right?

I know I keep doing this, but you introduced a new notation here which is $X$. To eliminate any confusion, and to make the notations consistent, I think the integral we want to evaluate eventually for $K=3$ is

F_3(z) = F_{Z_3}(z)=\int_0^{z}F_{Z_2}(z-u_2)f_U(u_2)\,du_2

where $z$ here is $\gamma$ in the original formulation.

I need to read your post a couple more times to understand it. The code even in Python is not easy to follow.

Thanks

 

[EngWiPy]

@andrewkirk To be honest, I still don't understand the implementation. Let me explain what I do understand. We start with $Z_1 = U_1$, then we have $F_{Z_1}(z) = F_U(z) = 1-1/(1+z)$.

Then we can formulate $Z_2$ as $Z_2 = Z_1 + U_2$, and thus

F_{Z_2}(z)=\int_{u_2=0}^zF_{Z_1}(z-u_2)f_U(u_2)\,du_2

This integration is easy to evaluate, because we have $F_{Z_1}(z)$ in closed form. We can define a vector $\mathbf{u}_2= 0,\,du,\ldots,\,(N-1)du$, where $du$ is the step size, and $N = \lceil z/du \rceil$ (I assume a uniform step size for the simplicity of exposition). And using vectorization, the above integral can be evaluated as

\text{sum}\left(F_{Z_1}(z-\mathbf{u}_2)*f_U(\mathbf{u}_2)*du\right)

After than I get lost. I mean we write $Z_3 = Z_2 + U_3$, and

F_{Z_3}(z)=\int_{u_3=0}^zF_{Z_2}(z-u_3)f_U(u_3)\,du_3

but the implementation will be

\text{sum}\left(F_{Z_2}(z-\mathbf{u}_3)*f_U(\mathbf{u}_3)*du\right) = \text{sum}\left(\left[F_{Z_1}(z-\mathbf{u}_3-\mathbf{u}_2)*f_U(\mathbf{u}_3)*du\right]*f_U(\mathbf{u}_2)*du\right)

where $\mathbf{u}_3-\mathbf{u}_2$ is all possible values of $u_3$ minus all possible values of $u_2$. But we need o make sure that $z-u_3-u_2 > 0$, so, we can write the argument as $\min(0, z-u_3-u_2)$, or loop as

for m1 in len(u): #u = 0, du, 2*du, ... , (N-1)*du
      for m2 in len(u)-m1:
            z = z - u[m1] - u[m2]

Which I can evaluate. For $u_3 = 0$, $F_{Z_1}(z-u_3-u_2)*f_U(u_2)*du$ is already evaluated in the previous step. Am I in the correct direction? I know it is a messy presentation, but it is just because I don't understand the approach very well. So, bear with me, please.

When you implemented the new approach in R, what is the speed up factor compared to the original recursive code also in R?

 

[andrewkirk]

Where your exposition differs from my algorithm is that you are making $du_k$ dependent on $z$ whereas in my algorithm $the step size$du## is selected at the beginning by a loop ('while x > max_dp') that selects it to be small enough that $f(u)du<\textit{max_dp}$ for some pre-set tolerance max_dp. Perhaps the algorithm can be done with varying increments but it would make things unnecessarily complex.

Define a function simp that is given an odd number of consecutive values of a function at evenly spaced increments du, and returns the Simpson's Rule estimate of the integral of the function over the relevant interval.
Then define another function x_simp (Extended Simpson) that is given any finite number of values of the function and estimates the integral as:

• zero if there is only one value
• simp applied to those values if there is an odd number of values
• the trapezoidal integral for the first incremental volume, added to simp applied to the odd sequence of values obtained by dropping the first one, if there's a nonzero, even number of values

Let the practical support of $f$ be $[0,\textit{max_lim})$, so that we can ignore the probability of $U$ exceeding max_lim. Then we define M to be max_lim / du rounded up to the next integer. At each nesting level k, we will estimate $F_{Z_k}(u)$ for $u\in\{0, du, ...,M\,du\}$. The elements of this estimate make up a vector $(A_0^k, ..., A_M^k)$, and we can calculate this directly for $k=1$ as $A_r^k=F_U(r\, du)$.

Then for $k>1$ and natural number $j\le M$ we have:

$$F_{Z_k}(j\,du) =
\int_0^z F_{Z_{k-1}}(z-j\,du)f_U(j\,du)\,du
\approx
\textit{x_simp}(A_0^{k-1}f_U(j\,du), A_1^{k-1}f_U((j-1)du), ..., A_{j-1}^{k-1}f_U(du), A_j^{k-1}f_U(0))
$$

Note how we store in the vector of $A^{k-1}_r$s all the information we'll need later about $F_{Z_{k-1}}$ when we do level $k-1$. Then we use that, together with the values of $f_U$ at multiples of $du$, to calculate the next set of $A^k_r$s.

So no recursive expansions are ever needed. In your post you've done a recursive expansion over two levels, so that your last formula contains both $u_2$ and $u_3$. Rather than recursively expanding, we just access the A values from teh previous level and use them in the Simpson integration estimate.

I hope that's clearer.

As I recall, the recursive version in R ran in about 3 mins 30, while the new version runs in about seven seconds, so the speed-up factor is about thirty. That's because the recursive version, at levels below the top one, is always redoing calcs that have already been done, whereas this approach avoids that.

 

[EngWiPy]

Oh, really? I am doing recursion again I told you I am still not sure about the implementation. As usual, I will have to read and analyse your post to understand it. It is very complex for me, to be honest.

A speed-up factor of 30x compared to the original is great. I believe we can get another 15-30x speed up factor if we use C++ using the same algorithm. Imagine when I have to evaluate not (K-1)-nested integrals, but 2*(K-1)-nested integrals, when I consider the order statistics case as we were discussing earlier on another thread! Computational speed is a stumbling block here.