Speeding up scipy.integrate.dblquad for updating Bayesian priors (credit rating probability of default model)

[binomial-torrent]

Hi all,

I am working on a model to create a transition matrix for rating transitions of various credit portfolios. I assume a normal prior distribution (standardized for values only between 0 and 1, which is the range of the enumerated credit rating values) of credit transitions given parameters μ and σ, as well as a normal (standardized to be between 0 and 1) prior density over the parameters μ and σ themselves. I know normal isn't the best assumption, but I just want to get a good first order analysis of the data. Anyways, I am following the method described in Chapter 8.3 of The Elements of Statistical Learning (page 268 in the linked pdf). To integrate over the parameter space and update the prior for the transition density function of each rating, I have to compute the integral Pr(z_new|Z) = ∫ Pr(z_new|θ) · Pr(θ|Z) dθ where θ = (μ,σ). I am doing this using scipy.integrate.dblquad; however, it takes astronomically long to update given one z, let along many observations Z. Does anyone know of a way to speed it up? I am attempting to use Cython without much luck. Should I discretize the parameter space θ along with its associated probability density, by using a discrete distribution such that P(θ) ~ O(e^(|θ-θ_mean|^2))?

 

[Daniel Duffy]

page 268, BTW (not sure if your link is ethical.. have you checked copyright???).

Can you post the core code to give an idea how you are doing it?

You might be using QUADPACK wrong because Fortran is the fastest language on the planet.

 

[binomial-torrent]

page 268, BTW (not sure if your link is ethical.. have you checked copyright???).

Can you post the core code to give an idea how you are doing it?

You might be using QUADPACK wrong because Fortran is the fastest language on the planet.

The link is legal, the pdf is made available from the Stanford website.

Here is the code:[CODE lang="python"]from scipy.integrate import quad, dblquad

from math import isnan

prev_r2 = defaultdict(lambda: defaultdict(str))

posterior_probabilities_discretized = defaultdict(lambda: np.zeros((num_different_ratings,num_different_ratings)))

g_prior = defaultdict(dict)
g_given_prev_r = defaultdict(lambda: defaultdict(list))

prior_dist = defaultdict(dict)

r2s = defaultdict(lambda: defaultdict(list))

f_mu_sigma = defaultdict(lambda: defaultdict(list))

for boid in boids:
for r1 in universe_of_ratings:
g_prior[boid][r1] = g[boid][rating_states_enum[r1]]
prior_dist[boid][r1] = lambda r2: f(r2, mean_prior[rating_states_enum[r1]], sigma_prior)

g_given_prev_r_dict = defaultdict(lambda: defaultdict(list))
f_g_convolution_dict = defaultdict(lambda: defaultdict(list))
P_r2_new_given_prev_r2s_dict = defaultdict(lambda: defaultdict(list))

for boid in boids:
for r1 in universe_of_ratings:
if r1 in df_by_BOID[boid]["RATUSED"].unique():
if not df_by_BOID[boid][df_by_BOID[boid]["RATUSED"]==r1].empty:
transitions = df_by_BOID[boid][df_by_BOID[boid]["RATUSED"]==r1]["RATUSED_transition"]
for idx, transition in enumerate(transitions):
if type(transition)==tuple:
r2_ = transition[1]
r2s[boid][r1].append(r2_)

r2s_idx = max(0,len(r2s[boid][r1])-1)

# This O(n) recursion is to express g^<n>(theta) = P(theta|r2_0, r2_1,...,r2_n) as a function
# of g^<n-1>(theta) = P(theta|r2_0, r2_1,...,r2_(n-1)), more specifically,
# g^<n>(theta) = f(r2_n|theta)*g^<n-1>(theta) / ∫ f(r2_n|theta)*g^<n-1>(theta) d(theta)
# where f(r2_n|theta) = P(r2_n|theta).
# The base case is g^<0>(theta) = f(r2_0|theta)g^<prior>(theta)

def g_given_previous_n_r2s(mu, sigma, n):
if n==0:
print("recursion depth: " + str(n))

standardizing_const = dblquad(lambda mu, sigma: f(r2s[boid][r1][0], mu, sigma) * g_prior[boid][r1](mu, sigma), 0, 1, lambda sigma: 0, lambda sigma: 1)[0]
print(f(r2s[boid][r1][0], mu, sigma) * g_prior[boid][r1](mu, sigma) / standardizing_const)
return f(r2s[boid][r1][0], mu, sigma) * g_prior[boid][r1](mu, sigma) / standardizing_const

else:
print("recursion depth: " + str(n))
### The double integral part
standardizing_const = dblquad(lambda mu, sigma: f(r2s[boid][r1][n], mu, sigma) * g_given_previous_n_r2s(mu, sigma, n-1), 0, 1, lambda sigma: 0, lambda sigma: 1)[0]
print(standardizing_const)
return f(r2s[boid][r1][n], mu, sigma) * g_given_previous_n_r2s(mu, sigma, n-1) / standardizing_const

for r2 in universe_of_ratings:
print("length of transitions: " + str(r2s_idx))
posterior_probabilities_discretized[boid][rating_states_enum[r1]][rating_states_enum[r2]] = dblquad(lambda mu, sigma, r, n: f(r, mu, sigma) * g_given_previous_n_r2s(mu, sigma, n), 0, 1, lambda sigma: 0, lambda sigma: 1, args=(r2,r2s_idx))

else:
for r2 in universe_of_ratings:
posterior_probabilities_discretized[boid][rating_states_enum[r1]][rating_states_enum[r2]] = prior_dist[boid][r1](r2)[/CODE]

I'm just implementing the base case right now, the recursive calls to update the priors will be even more complex.

 

[Daniel Duffy]

I copied and pasted the code but it does not run for me. Can you attach the .py file please?

 

[binomial-torrent]

I copied and pasted the code but it does not run for me. Can you attach the .py file please?

Just sent you the code and data.

 

[Daniel Duffy]

I had an initial look at the posted code. My feeling was QUAD is not the problem and I still think so. The code has a number of language-independent issues and in the Python case might make things worse because everything is done at run-time.

First impressions.

0, Undocumented code. Hard to understand/modify. Especially a) if it's someone else, b) yourself after 2-3 months gone from long-term memory..
1. Recursion is the darling of CS textbooks, but in numerical analysis it is not at all efficient. This is a bottleneck for sure. Why not use an iterative solution?
2. I think the loops can be optimized,

Loop optimization - Wikipedia

en.wikipedia.org

NB loop unswitching.
3. Lots of nested arrays and indexing, not sure if dict is a good idea (are its keys hashable). Plan B use standard numeric data structures.
AFAIR access in dict is O(n log n). Would a matrix be more appropriate and efficient?

NumPy - Matrix Library

NumPy - Matrix Library - NumPy package contains a Matrix library numpy.matlib. This module has functions that return matrices instead of ndarray objects.

www.tutorialspoint.com

4. ? use numba to parallelise the loops.

Numba - Wikipedia

en.wikipedia.org

Anyway, that would be my approach. hth

 

[binomial-torrent]

Hmm interesting. In order to debug I did put in print statements. Don’t think the recursion was a bottleneck as much as the fact that dblquad was calling the lowest recursive depth, 0 (evaluating the prior function many times as part of the integration process presumably). This is because when I discretized μ and σ and their distribution, the problem goes away and the recursion bubbles back up really quickly according to debug print statements.

I’ll see what I can do with numba/cython. Pretty sure dictionary access is O(1) more or less in Python because of hashing. Loop optimization I’ll look into as well, but the main bottleneck was huge chunks of time between successive calls to the integrand function when calling dblquad.

 

[Daniel Duffy]

I forgot to ask why the def g_given_previus is defined _inside_ the loop?

 

[binomial-torrent]

I forgot to ask why the def g_given_previus is defined _inside_ the loop?

it depends on the specific sequence of r2s that come after r1. basically, it is P_boid_r1(theta|r2_0, r2_1,...,r2_n), where the sequence being conditioned on is {r2 for r2 in (r1, r2)}, (r1,_) being a subset of the transitions with r1 fixed. so it depends on r1 and boid and their specific sequence of relevant transitions.

 

[Daniel Duffy]

it depends on the specific sequence of r2s that come after r1. basically, it is P_boid_r1(theta|r2_0, r2_1,...,r2_n), where the sequence being conditioned on is {r2 for r2 in (r1, r2)}, (r1,_) being a subset of the transitions with r1 fixed. so it depends on r1 and boid and their specific sequence of relevant transitions.

I'll take your word for it. Since the code is not documented in such a way that I can see he guts of the corresponding algorithm.
I would do a new design.

The problem is forward<->backward traceability (algorithm<->code). It's just difficult to follow things. This entails that profiling performance is also different.

Instead, call a 'stub' (empty) quad function and measure performance. Then you know where the bottleneck is.

 

[binomial-torrent]

I’m going to write the documentation today through to tomorrow actually.

 

[Daniel Duffy]

On a related, follow-on remark..

In general, academia pays little attention to designing maintainable software but in industry it is vital. We all know the bed-time Rocky Horror Software Show.

Languages such as C++ and C# force you to be disciplined but in Python you can do anything at any time, anywhere. This is a maintenance concern going forward if you want to support "production" Python code.

Plan B; find Python library (C++, Fortran) that does the job and use it with scripting commands.Horror

 

[Daniel Duffy]

"but the main bottleneck was huge chunks of time between successive calls to the integrand function when calling dblquad."

Have you tried profiling using e.g. cProfile?