# [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

```