Complicated nested integrals in terms of python scipy.integrate

I want to write the following complicated function using python scipy.integrate:

where \phi and \Phi are the pdf and cdf of the standard normal respectively and the sum i != 1,2 is over the thetas in *theta, and the product s != 1,i,2 is over the thetas in *theta that is not theta_i

To make things simpler, I break down the problem into a few easier pieces by defining a few more functions:

And here is what I have tried:

import scipy.integrate as integrate
from scipy.integrate import nquad
from scipy.stats import norm
import math
import numpy as np

def f(v, u, t1, ti, t2, *theta):
  prod = 1
  for t in theta:
    prod = prod * (1 - norm.cdf(u + t2 - t))
  return norm.pdf(v) * norm.cdf(v + ti - t1) * prod

def g(v, u, t1, t2, *theta):
  S = 0
  for ti in theta:
    S = S + integrate.quad(f, -np.inf, u + t2 - ti, args=(u, t1, ti, t2, *theta))[0]
  return S * norm.pdf(u)

def P(t1, t2, *theta):
  return integrate.quad(g, -np.inf, np.inf, args=(t1, t2, *theta))[0]

But it doens't work because I think in the definition of P the integral is integrating w.r.t. the first variable, namely v but we should be integrating u but I don't know how we can do that.

Is there any quick way to do it?

For checking, here is what the correct P would produce:

P(0.2, 0.1, 0.3, 0.4) = 0.08856347190679764

P(0.2, 0.1, 0.4, 0.3, 0.5) = 0.06094233268837703

The first thing I did was move prod outside of f since it isn't a function of v. I realized that your prod includes the ti term, so I added a check for ts (renamed from t) equal to ti.

You also have way too many variables for your functions. g is a function of u and the thetas and f is a function of v and the thetas, so they don't both need u and v.

This isn't the most efficient code, but it works.

from scipy.integrate import quad
from scipy.stats import norm
import numpy as np

def f(v, t1, ti):
  return norm.pdf(v) * norm.cdf(v + ti - t1)

def g(u, t1, t2, *theta):
  S = 0
  for ti in theta:
      prod = 1
      for ts in theta:
          if ts != ti:
              prod *= (1 - norm.cdf(u + t2 - ts))
      S += prod*quad(f, -np.inf, u + t2 - ti, args=(t1, ti))[0]
  return S * norm.pdf(u)

def P(t1, t2, *theta):
  return quad(g, -np.inf, np.inf, args=(t1, t2, *theta))[0]

print(P(0.2, 0.1, 0.3, 0.4))            # 0.08856347187084088
print(P(0.2, 0.1, 0.4, 0.3, 0.5))       # 0.060942332693157915