[SciPy-user] Puzzle with Monte Carlo error estimation

[Iain Day]

Hi,

I've got a puzzle doing some Monte Carlo stuff to get standard 
deviations for some fitting parameters I've calculated. The code is 
below. The problem is that p1 is changing each time optimize.leastsq is 
called in the loop and I can't see why. Any ideas?

Thanks,

Iain



from scipy import *

pts = linspace(0, 10000.0, 1000)
expdata = 3500.0 * (1.0 - exp(-pts / 565.0)) + 1870.0 + rand(len(pts))
p0 = [4000.0, 600.0, 2000.0]
fix0 = []
num_mc = 1000 # This is the number of MC runs

p1, success = optimize.leastsq(residuals, p0, args=(fix0, pts, expdata))

sigma = std(residuals(p0, fix0, pts, expdata))
mocked_para = zeros((num_mc, len(p0)), dtype=float)

for i in range(num_mc):
     mock_data = buildup(p1, fix0, time_points) + sigma * 
rand(len(time_points))
     p2, success = optimize.leastsq(residuals, p1, args=(fix0, ts, expdata))
     mocked_parameters[i] = p2

mc_errors = std(mocked_parameters, axis=0)


def buildup(varpars, fixpars, time):
     amplitude, tau, offset = varpars
     buildup = amplitude * (1.0 - exp(- time / tau)) + offset
     return buildup

def residuals(varpars, fixpars, time, expdata):
    err = expdata - buildup(varpars, fixpars, time)
    return err