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.
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, 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 the montecarlo waves that you will actually fit duplicate/o yy, root:packages:motofit:old_genoptimise:y_montecarlo duplicate/o xx, root:packages:motofit:old_genoptimise:x_montecarlo duplicate/o ee, root:packages:motofit:old_genoptimise:e_montecarlo Wave y_montecarlo = root:packages:motofit:old_genoptimise:y_montecarlo Wave x_montecarlo = root:packages:motofit:old_genoptimise:x_montecarlo Wave e_montecarlo = root:packages:motofit:old_genoptimise:e_montecarlo //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) y_montecarlo[] = yy[p] + gnoise(ee[p]) Gencurvefit/q/n/X=x_montecarlo/I=1/W=e_montecarlo/K={gen.GEN_generations, gen.GEN_popsize,gen.k_m, gen.GEN_recombination}/TOL=(gen.GEN_V_fittol) $fn, y_montecarlo, 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
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.
You are right the usage
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.
A small point- you might
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