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