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) < 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[loop] b = points[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.append(a) selected_parameters.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) last_a_std = np.std(selected_parameters) last_b_mean = np.mean(selected_parameters) last_b_std = np.std(selected_parameters) # 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))