Monte Carlo error analysis using Genetic Optimisation.

 The following code requires that you have the genetic optimisation package installed (http://www.igorexchange.com/project/gencurvefit).
The snippet uses Monte Carlo resampling to estimate fit uncertainties and parameter correlations for a given set of data and a fit function. (see, for example http://www.casaxps.com/help_manual/error_analysis.htm)

It works by fitting datasets that are synthesized from the original input data. Each synthesized dataset has new ordinate values (yvalues) that are generated by adding random gaussian noise to each datapoint, with the standard deviation of the gaussian distribution being the same as the experimental uncertainty (error) for that given datapoint. To obtain sensible results many (e.g. 1000) fits are done using genetic optimisation.

Useage:
Moto_montecarlo("myfitfunction", myydata, myxdata, myedata, "1010", 1000)

Outputs:
M_correlation - Correlation Matrix
W_sigma954 - uncertainty in fitted parameter
M_montecarlo - fitted parameters for each fit iteration. (you can histogram all the values obtained for a particular parameter to see their distribution)

Function Moto_montecarlo(fn, w, yy, xx, ee, holdstring, Iters)
    String fn
    Wave w, yy, xx, ee
    string holdstring
    variable Iters
    //useage:
   
    string cDF = getdatafolder(1)
    variable ii,jj,kk, summ

    newdatafolder/o root:packages
    newdatafolder/o root:packages:motofit
    newdatafolder/o root:packages:motofit:old_genoptimise

    try
        //get initialisation parameters for genetic optimisation
        struct GEN_optimisation gen
        gen.GEN_Callfolder = cDF
        GEN_optimise#GEN_Searchparams(gen)
   
        //get limits
        GEN_setlimitsforGENcurvefit(w, holdstring, cDF)
        Wave limits = root:packages:motofit:old_genoptimise:GENcurvefitlimits
   
        //make a wave to put the montecarlo iterations in
        make/o/d/n=(Iters, dimsize(w, 0)) M_montecarlo

        //now lets do the montecarlo fitting
        variable timed = datetime
        for(ii=0 ; ii<Iters ; ii+=1)
            Gencurvefit/MC/q/n/X=xx/I=1/W=ee/K={gen.GEN_generations, gen.GEN_popsize,gen.k_m, gen.GEN_recombination}/TOL=(gen.GEN_V_fittol) $fn, yy, w, holdstring, limits
            M_montecarlo[ii][] = w[q]
            print "montecarlo ", ii, " done - total time = ", datetime-timed
        endfor
        print "overall time took: ", datetime - timed , " seconds"
   
        //now work out correlation matrix and errors.
        //see Heinrich et al., Langmuir, 25(7), 4219-4229
        make/n=(dimsize(w, 0))/o W_sigma954, means, stdevs
        make/n=(dimsize(w,0), dimsize(w, 0))/o M_correlation
        M_correlation = NaN
   
        for(ii = 0 ; ii<dimsize(w, 0) ; ii+=1)
            make/o/d/n=(Iters) goes
            goes = M_montecarlo[p][ii]
            Wavestats/alph=0.045501/M=2/q/w goes
            Wave M_wavestats
            W_sigma954[ii] = M_wavestats[25]- M_wavestats[24]
            means[ii] = M_wavestats[3]
            stdevs[ii] = M_wavestats[4]
            if(stringmatch(holdstring[ii], "1"))
                W_sigma954[ii] = NaN
            endif
        endfor
        for(ii=0 ; ii< dimsize(w, 0) ; ii+=1)
            for(jj= ii ; jj<dimsize(w,0) ; jj+=1)
                if(ii==jj || stringmatch(holdstring[ii], "1") || stringmatch(holdstring[jj], "1"))
                    M_correlation[ii][jj]=NaN
                else           
                    summ = 0
                    for(kk = 0 ; kk < Iters ; kk+=1)
                        summ += (M_montecarlo[kk][ii]-means[ii])*(M_montecarlo[kk][jj]-means[jj])
                    endfor
                    M_correlation[ii][jj] = summ / (Iters-1) / (stdevs[ii] * stdevs[jj])
                endif  
                M_correlation[jj][ii] = M_correlation[ii][jj]
            endfor
        endfor
    catch

    endtry
    killwaves/z M_wavestats, goes, means, stdevs
    setdatafolder $cDF
End
Hello,
I need help applying this snippet for MC estimation. I copied the code into the procedure window and then entered the usage line into the command window. I used "Motofit" for "myfitfunction", the names of my x, y and error datasets, and 1000 iterations. I got a command error: "Expected wave name, variable name, or operation." In the usage line in the command window, "Moto_montecarlo" is highlighted. I noticed that the usage line doesn't match the function-could this be the problem? I would appreciate any suggestions.

Thanks.

April 15, 2010 at 08:14 pm - Permalink

You are right the usage instructions aren't the same as the snippet. The usage should be:

Moto_montecarlo("myfitfunction", mycoefs, myydata, myxdata, myedata, "1010", 1000)

If you have the latest version of Motofit installed (see http://dav1-platypus.nbi.ansto.gov.au/ for links) this Moto_montecarlo procedure will already be there. You have to select the "Genetic + MC analysis" fitting method from the reflectivity panel. It produces the same output.

Please note:
1) The Monte Carlo error analysis does not increase the likelihood of the fit being correct, the only advantage is to show the probability density of each fitted parameter.
2) If the parameter uncertainty is normally distributed then doing the MC analysis only gives information that could have been obtained by Levenberg Marquardt anyway.
3) To produce a good fit the MC analysis relies on the correct choice of model, just as much as a normal fit method does.
4) MC analysis ignores systematic errors as well, etc.

April 16, 2010 at 12:40 am - Permalink

A small point- you might want to add the optional second parameter for gnoise(stddev, rng) and set it to 2 to use the Mersenne Twister RNG. It has better statistical characteristics than the standard one used if you don't supply that parameter. The standard RNG is pretty good, but MT is better.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

April 16, 2010 at 10:22 am - Permalink

I have been trying to implement the Monte Carlo routine with a custom fitting function. I found the updated Moto_montecarlo function in the Genetic Optimisation procedure (geneticOptimization.ipf). However, the updated function does not seem to call the genetic optimization procedure.

At the end of the run, the rows in the M_Montecarlo matrix are all the same and the figure with the correlations contain only a single point per box.

Do I need to call a different function other than Moto_montecarlo to initialize the Monte Carlo routine?

March 19, 2013 at 12:34 pm - Permalink

Forum

Support

Gallery

Igor Pro 9

Learn More

Igor XOP Toolkit

Learn More

Igor NIDAQ Tools MX

Learn More