Determine the hyper-parameters of an Approximate Bayesian computation

Question

I created my own Approximate Bayesian computation algorithm in Python 3 and I try to find the best hyper-parameters for it but it is really slow. This question is not about the math involved in it. For that, you can read my other question on Cross Validated.

The objective of the ABC algorithm is to find the best distribution for the parameters of a model by using a Bayesian approach. To do that, I create a dataset with the real parameters of my model and I try to find all the values of \$a\$ and \$b\$ that produce the same outcome. Then, I use those new values of \$a\$ and \$b\$ to create my posterior distribution.

The purpose of my piece of code is to find the hyper-parameters \$k\$ and \$n\$ who give me the best compromise between precision and time of execution. The \$k\$ parameter corresponds to the size of the dataset I generate at each step and the \$n\$ parameter corresponds to the number of points I generate from my prior distribution.

For those two hyper-parameters, the higher the better but I want to check if increasing one of them has a big enough impact on the standard deviation of the parameter of my model.

So, here is my question: what can I do to improve the speed of my algorithm?

import numpy as np
from scipy.stats import norm
from scipy.stats import uniform
from scipy.stats import lognorm
from scipy.stats import gaussian_kde
import matplotlib.pyplot as plt

# Because I want to find the best hyper-parameters for my model,
# I will create a dataset with fixed parameters
real_a = 0.9
real_b = -0.02

# Make m steps for each set of hyper-parameters
m = 100

# Return the probability that a Bernoulli trial will success
def gen_model(a, b, odd):
    prob = a / odd + b
    prob = 0.0 if prob < 0.0 else prob
    prob = 1.0 if prob > 1.0 else prob
    return prob

def bernoulli_trial(p):
    return np.random.rand() <= p

# Return the outcome of a Bernoulli trial with the fixed parameters
def real_exp(odd):
    return bernoulli_trial(gen_model(real_a, real_b, odd))

# k is one of my hyper-parameters. It represents the number of rows of my dataset
for k in range(5,11):
    # n = 10 ** p is the second hyper-parameter. It represents the number of tuples (a,b) I generate for each run
    for p in range(5,8):
        n = 10 ** p
        # Three results for each hyper-parameters tuple
        for l in range(3):
            # The first time, I pick the parameters from a normal distribution
            first = True

            posterior = None

            last_a_mean = 1.0
            last_a_std = 0.0

            last_b_mean = 0.0
            last_b_std = 0.0

            for run in range(m):
                selected_parameters = [[], []]

                # Accept the parameters only if at least 100 are non-rejected
                while(len(selected_parameters[0]) < 100):
                    current_loop_odd = []
                    current_loop_result = []

                    for idx in range(k): # Creation of my dataset
                        current_loop_odd.append(1 / uniform.rvs(0.01, 0.99))
                        current_loop_result.append(real_exp(current_loop_odd[idx]))

                    for loop in range(n):
                        add_selected_parameters = True

                        if first: # My prior is a and b are independent and follow a normal distribution (mean of 1 for a and mean of 0 for b)
                            a = norm.rvs(1.0,0.1)
                            b = norm.rvs(0.0,0.1)
                        else: # For the next steps, I use datapoints generated from my posterior distribution
                            a = points[0][loop]
                            b = points[1][loop]

                        for idx in range(k): 
                            if bernoulli_trial(gen_model(a, b, current_loop_odd[idx])) != current_loop_result[idx]:
                                add_selected_parameters = False # The tuple is rejected because it doesn't correspond to my dataset

                        if add_selected_parameters:                
                            selected_parameters[0].append(a)
                            selected_parameters[1].append(b)

                posterior = gaussian_kde(selected_parameters) # I generate a new posterior distribution from my non-rejected tuples
                points = posterior.resample(n) # And I create new datapoints for the next step
                first = False

                last_a_mean = np.mean(selected_parameters[0])
                last_a_std = np.std(selected_parameters[0])

                last_b_mean = np.mean(selected_parameters[1])
                last_b_std = np.std(selected_parameters[1])

            # Print the result
            print(str(k) + " - " + str(n) + " : (" + str(l) + ") -> " + str(last_a_mean) + " , " + str(last_a_std) + " / " + str(last_a_mean) + " , " + str(last_a_std))

Answer 1

what can I do to improve the speed of my algorithm?

You have an (impressive) seven-deep nested loop, where only one of the levels of the nested loop takes up to \$10^7\$ iterations. Taking a guess at iteration count order of magnitude,

$$ (11 - 5) (8 - 5) \cdot 3 \cdot 100 \cdot 10^7 \cdot 10 = 5.4 \cdot 10^{12} $$

It won't exactly be equal to this, because you have one while loop whose scale is unclear to me, and \$k\$ and \$n\$ are variable.

Even if this were well-written in Python using perfect vectorisation, it would be slow; you need to abandon it. If you can't redesign your algorithm to cut some of these iterations down, then your only probable hope is to drop down to C/C++, and do at least the bare minimum in multiprocessing though threads or ideally GPGPU.