[SciPy-user] Puzzle with Monte Carlo error estimation

Iain Day iain at day-online.org.uk.invalid
Sun Nov 12 15:52:42 CST 2006

```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

```