Population Monte Carlo implementation

Question

I am trying to implement the Population Monte Carlo algorithm as described in this paper (see page 78 Fig.3) for a simple model (see function model()) with one parameter using Python. Unfortunately, the algorithm doesn't work and I can't figure out what's wrong. See my implementation below. The actual function is called abc(). All other functions can be seen as helper-functions and seem to work fine.

To check whether the algorithm workds, I first generate observed data with the only parameter of the model set to param = 8. Therefore, the posterior resulting from the ABC algorithm should be centered around 8. This is not the case and I'm wondering why.

I would appreciate any help or comments.

# imports

from math import exp
from math import log
from math import sqrt
import numpy as np
import random
from scipy.stats import norm


# globals
N = 300              # sample size
N_PARTICLE = 300      # number of particles
ITERS = 5            # number of decreasing thresholds
M = 10               # number of words to remember
MEAN = 7             # prior mean of parameter
SD = 2               # prior sd of parameter


def model(param):
  recall_prob_all = 1/(1 + np.exp(M - param))
  recall_prob_one_item = np.exp(np.log(recall_prob_all) / float(M))
  return sum([1 if random.random() < recall_prob_one_item else 0 for item in range(M)])

## example
print "Output of model function: \n" + str(model(10)) + "\n"

# generate data from model
def generate(param):
  out = np.empty(N)
  for i in range(N):
    out[i] = model(param)
  return out

## example

print "Output of generate function: \n" + str(generate(10)) + "\n"


# distance function (sum of squared error)
def distance(obsData,simData):
  out = 0.0
  for i in range(len(obsData)):
    out += (obsData[i] - simData[i]) * (obsData[i] - simData[i])
  return out

## example

print "Output of distance function: \n" + str(distance([1,2,3],[4,5,6])) + "\n"


# sample new particles based on weights
def sample(particles, weights):
  return np.random.choice(particles, 1, p=weights)

## example

print "Output of sample function: \n" + str(sample([1,2,3],[0.1,0.1,0.8])) + "\n"


# perturbance function
def perturb(variance):
  return np.random.normal(0,sqrt(variance),1)[0]

## example 

print "Output of perturb function: \n" + str(perturb(1)) + "\n"

# compute new weight
def computeWeight(prevWeights,prevParticles,prevVariance,currentParticle):
  denom = 0.0
  proposal = norm(currentParticle, sqrt(prevVariance))
  prior = norm(MEAN,SD)
  for i in range(len(prevParticles)):
    denom += prevWeights[i] * proposal.pdf(prevParticles[i])
  return prior.pdf(currentParticle)/denom


## example 

prevWeights = [0.2,0.3,0.5]
prevParticles = [1,2,3]
prevVariance = 1
currentParticle = 2.5
print "Output of computeWeight function: \n" + str(computeWeight(prevWeights,prevParticles,prevVariance,currentParticle)) + "\n"


# normalize weights
def normalize(weights):
  return weights/np.sum(weights)


## example 

print "Output of normalize function: \n" + str(normalize([3.,5.,9.])) + "\n"


# sampling from prior distribution
def rprior():
  return np.random.normal(MEAN,SD,1)[0]

## example 

print "Output of rprior function: \n" + str(rprior()) + "\n"


# ABC using Population Monte Carlo sampling
def abc(obsData,eps):
  draw = 0
  Distance = 1e9
  variance = np.empty(ITERS)
  simData = np.empty(N)
  particles = np.empty([ITERS,N_PARTICLE])
  weights = np.empty([ITERS,N_PARTICLE])

  for t in range(ITERS):
    if t == 0:
      for i in range(N_PARTICLE):
        while(Distance > eps[t]):
          draw = rprior()
          simData = generate(draw)
          Distance = distance(obsData,simData)

        Distance = 1e9
        particles[t][i] = draw
        weights[t][i] = 1./N_PARTICLE

      variance[t] = 2 * np.var(particles[t])
      continue


    for i in range(N_PARTICLE):
      while(Distance > eps[t]):
        draw = sample(particles[t-1],weights[t-1])
        draw += perturb(variance[t-1])
        simData = generate(draw)
        Distance = distance(obsData,simData)

      Distance = 1e9
      particles[t][i] = draw
      weights[t][i] = computeWeight(weights[t-1],particles[t-1],variance[t-1],particles[t][i])


    weights[t] = normalize(weights[t])  
    variance[t] = 2 * np.var(particles[t])

  return particles[ITERS-1]



true_param = 9
obsData = generate(true_param)
eps = [15000,10000,8000,6000,3000]
posterior = abc(obsData,eps)
#print posterior

Answer 1

I stumbled upon this question as I was looking for pythonic implementations of PMC algorithms, since, quite coincidentally, I'm currently in the process of applying the techniques in this exact paper to my own research.

Can you post the results you're getting? My guess is that 1) you're using a poor choice of distance function (and/or similarity thresholds), or 2) you're not using enough particles. I may be wrong here (I'm not very well-versed in sample statistics), but your distance function implicitly suggests to me that the ordering of your random draws matters. I'd have to think about this more to determine whether it actually has any effect on the convergence properties (it may not), but why don't you simply use the mean or median as your sample statistic?

I ran your code with 1000 particles and a true parameter value of 8, while using the absolute difference between sample means as my distance function, for three iterations with epsilons of [0.5, 0.3, 0.1]; the peak of my estimated posterior distribution seems to be approaching 8 just like it should on each iteration, alongside a reduction in the population variance. Note that there is still a noticeable rightward bias, but this is because of the asymmetry of your model (parameter values of 8 or less can never result in more than 8 observed successes, while all parameters values greater than 8 can, leading to a rightward skewedness in the distribution).

Here's the plot of my results: