import numpy as np
import matplotlib.pyplot as pl
%matplotlib inline

data = 5.0 + 2.0 * np.random.randn(10)
meanD = np.mean(data)
stdD = np.std(data)

def abcAcceptance(pars):
    if (pars[1] < 0):
        return False
    samples = pars[0] + pars[1]*np.random.randn(10)
    
    diffMean = np.abs(meanD - np.mean(samples))
    diffStd = np.abs(stdD - np.std(samples))
    if ((diffMean < 0.1) and (diffStd < 0.2)):
        return True
    else:
        return False

start = np.asarray([4.0,2.3])
chain = np.zeros((2,25001))
chain[:,0] = start
for i in range(25000):
    proposal = chain[:,i] + 0.7 * np.random.randn(2)
    
    if (abcAcceptance(proposal)):
        chain[:,i+1] = proposal
    else:
        chain[:,i+1] = chain[:,i]

pl.plot(chain[0,:],'.')

[<matplotlib.lines.Line2D at 0x562d210>]

pl.plot(chain[1,:],'.')

[<matplotlib.lines.Line2D at 0x5a25b90>]

print "{0} - {1}".format(np.mean(chain[0,:]),np.std(chain[0,:]))

5.71802874288 - 0.681588282195

print "{0} - {1}".format(np.mean(chain[1,:]),np.std(chain[1,:]))

2.25793396108 - 0.553244176481