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))

                    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:                

                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))
  • \$\begingroup\$ Explain at the top, specific requests at the top, code at the bottom. \$\endgroup\$ – FreezePhoenix Aug 2 at 12:54
  • \$\begingroup\$ A style note: Comments in code are generally incomplete - you could stick a noun in front of them, and it would work: # make m steps for each set of hyper-parameters however this suggestion does not qualify for an answer \$\endgroup\$ – FreezePhoenix Aug 2 at 12:56
  • \$\begingroup\$ @FreezePhoenix I put the code at the end of the question. \$\endgroup\$ – Pierre Aug 2 at 13:52
  • \$\begingroup\$ A general code-review.se rule: When someone suggests a change specific to the code, such as I did in my second comment, you don't edit your question. And when someone answers a question, you still don't edit it. This is to keep answers from having to change as comments are added to the post. \$\endgroup\$ – FreezePhoenix Aug 3 at 13:36
  • \$\begingroup\$ @FreezePhoenix I'm OK with what you say about changing the question after an answer but I didn't have any answer to my question. \$\endgroup\$ – Pierre Aug 4 at 8:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.