[SciPy-user] Puzzle with Monte Carlo error estimation

[Iain Day]

Iain Day wrote:
> 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


Okay, sorry to follow up my own post, but I've figured something out, 
and I'm not sure I understand it fully.

When I create p0, I create it as <type 'list'>, and optimize.leastsq 
returns p1 as <type 'numpy.ndarray'>.

If I pass p0 as <type 'numpy.ndarray'> then it gets changed to the best 
fit parameters, whereas if it is <type 'list'> it doesn't.

So, the fix for my code is to have:

p1 = list(p1)

just before the Monte Carlo run.

Have I missed (or misunderstood) something obvious here?

Thanks,

Iain