[SciPy-user] Puzzle with Monte Carlo error estimation

Iain Day iain at day-online.org.uk.invalid
Sun Nov 12 14:55:07 CST 2006


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



More information about the SciPy-user mailing list