Having some problems with the ABC SMC algorithm. I'm trying to implement the methods taken from here: Simulation-based model selection for dynamical systems in systems and population biology
How do weights scale down with increasing parameters?
I am having problems with this part of the weight calculation, w(m ,theta) . My current implementation is giving a higher weighting to particles with more parameters, which should not be the case.
In the function S, KP is a perturbation kernel that describes the probability of particle having come from the previous population. A small KP makes the denominator small, and therefore the overall weighting high.
Is it not the case that the more parameters that are present, the smaller KP becomes? I can't see anywhere in the paper that describes how the weighting is adjusted according to the number of parameters.
I've written a small example to illustrate the problem, and removed parts involving model perturbation that are included in the figures. Model a has one parameter, model b has six parameters. We assume they are all accepted and given weights. What am I missing to make the weighting proportional to the number of parameters?
def get_pdf_uniform(lower_bound, upper_bound, x):
if (x > upper_bound) or (x < lower_bound):
print(x)
print(upper_bound)
print(lower_bound)
return 0.0
else:
return (1 / (upper_bound - lower_bound))
def weight_problem_example():
# Two priors, one small one large
prior_a = [[1, 5]]
model_a = "a"
prior_b = [[1, 5], [1, 5], [1, 5], [1, 5], [1, 5], [1, 5]]
model_b = "b"
# Previous population had equal weights
prev_weight_a = 0.5
prev_weight_b = 0.5
# Sampled parameters are within prior range
particle_a_params = [np.random.uniform(x[0], x[1]) for x in prior_a]
particle_b_params = [np.random.uniform(x[0], x[1]) for x in prior_b]
particle_a_prev_params = [np.random.uniform(x[0], x[1]) for x in prior_a]
particle_b_prev_params = [np.random.uniform(x[0], x[1]) for x in prior_b]
packed_a = [prior_a, particle_a_params, particle_a_prev_params, prev_weight_a, model_a]
packed_b = [prior_b, particle_b_params, particle_b_prev_params, prev_weight_b, model_b]
particles = []
particles = particles + [packed_a for x in range(1000)]
particles = particles + [packed_b for x in range(1000)]
model_a_weights = []
model_b_weights = []
for part in particles:
# unpack particle data
prior = part[0]
params = part[1]
prev_params = part[2]
prev_weight = part[3]
model = part[4]
# Numerator
particle_prior_prob = 1
for idx, p in enumerate(params):
particle_prior_prob = particle_prior_prob * get_pdf_uniform(prior[idx][0], prior[idx][1], p)
# Assume equal model probability
particle_prior_prob = particle_prior_prob * 1
# Denominator, assuming we only have two particles
s = 0
for theta in particles:
if theta[4] == model:
theta_prior = part[0]
theta_params = part[1]
theta_prev_params = part[2]
theta_prev_weight = part[3]
theta_model = part[4]
k_prob = 1
for theta_idx, theta_p in enumerate(theta_prev_params):
# get parameter kernel pdf
k_prob = k_prob * get_pdf_uniform(
theta_p - theta_prior[theta_idx][0]*2,
theta_p + theta_prior[theta_idx][1]*2,
params[theta_idx])
# print(s)
s += (theta_prev_weight * k_prob)
part_weight = particle_prior_prob/s
if model == "a":
model_a_weights.append(part_weight)
if model == "b":
model_b_weights.append(part_weight)
print(np.mean(model_a_weights))
print(np.mean(model_b_weights))
