Numerit[WIN32][1700][1703]x 3#3 qffffff)@j@fffffvq@ffffff9@ffffff9@ffffff9@ffffff9@ffffff)@ffffff)@         Times New RomanArialSymbol Courier New G c@ffffffI@yfitx getdata:ytx>@>@>@>@>@>@>@ >@ >@ >@>@>@>@>@>@ >@  >==>>??02.50.40.010.0101020.010.01X 33333,@ @chisqG c@ffffffI@yfitxyx>@>?>@>@>@>@>@ >@ >@ >@>@>@>@>@>@ >@  >==>>??02.50.40.010.0101020.010.01X 33333,@ @meanvX 33333,@ @stdvvvvpvvvpvvvpvvvpvvvpvvvpvvv 3vvv.Nonlinear fitting to a sum of three Gaussians*vvv.using the Levenberg-Marquardt methodvvv vvv In this experiment we use a simulated data signal and fit to it a sum of three Gaussians. Each Gaussian may have any amplitude, position, and width. The data is created by taking a specific signal of three Gaussians and adding to it random noise. vvv Figure 1 displays the simulated data (in red) and the fitted curve (in black). When the program is running we can see the process of finding the best fit in action.vvv "# 8vvv ! . Data (red dots) and fitted curve (solid black) vvvpvvv It is interesting to see how close we get to the original signal which was used to create the simulated data, so in Fig. 2 we show this theoretical curve (red) together with the fitted curve (black). vvv ! Bvvv " . Theoretical (dashed red) and fitted (solid black) curves. vvv + vvv And finally, the quality of the fit is determined by the value of "c2 and the covariant matrix. These are presented below. The covariance matrix is represented by the mean value of its diagonal elements (which give the variance of each fitted parameter), and by their standard deviation:NvvvffB~B33B "c2  ! =" U Variance: 9 i Mean =$ B Std. Dev. =% ffvv " q8ffffff)@j@fffffvq@?@ffffff9@?@ffffff9@ffffff)@ffffff)@         Times New RomanArialSymbol Courier New ;` This sample program demonstrates an implementation of the=` Levenberg-Marquardt algorithm for fitting experimental data=` to a given model. The model we use in this example is a sum=` of three Gaussians with variable positions, amplitudes, and ` widths.9` This is an implementation of the algorithm as described` in Numerical Recipes 15.5.>` We tried to keep the algorithm as similar as possible to the>` one in NR while converting loops over data points into array ` operations.` Notes:=` 1. Here we assume that ALL the parameters should be fitted,C` unlike the case in NR where an array 'ia' defines which of theB` parameters are to be fitted and which are held fixed to their` input values.A` 2. The algorithm is demonstrated using a simulated signal whichD` is created as a sum of three Gaussians and a random noise. Each>` time you run the program, a new set of points is created.*` To fit a different model you need to:D` * Change the function 'model' which computes the model function` and its derivatives.E` * Change the function 'modelf' which computes the model functionF` only. This is used just for displaying the fitted signal and is7` the same as 'model' but without the derivatives.C` * Change the function 'getdata' to read data from a file or to(` create a proper simulated signal.<` * Define a proper initial guess in the function 'init'.:`********************************************************* ` Functions=` model: computes the model function and its derivatives with` respect to the parameters.&` Should be re-defined for each model.9` In this case: A sum of Gaussians each having a variableF` position, amplitude, and width (after example 15.5.16 NR, 'fgauss').I` Note that since we call the model function once for all the data pointsD` 'dyda' is a 2-dim array with dimensions [ndata x na] rather than aD` 1-dim array with 'na' elements. Similarly, 'y' is an array of size` 'ndata' and not a scalar. ` Input: x,a` Output: y,dydafunc model(x,a,y,dyda)na = length(a)!clear y ` set all elements to 0+for i = 1 to na by 3 ` loop over Gaussians*` some coefficients: to save computation` and array access timeai = a[i] ` amplitudeai1 = a[i+1] ` positionai2 = a[i+2] ` widthar = (x-ai1)/ai2ex = exp(-ar*ar)fac = 2*ai*ex*ar"` the function and its derivatives-y += ai*ex ` sum of GaussiansMdyda[*,i] = ex ` column i : derivative with respect to amplitudeLdyda[*,i+1] = fac/ai2 ` column i+1: derivative with respect to positionJdyda[*,i+2] = fac*ar/ai2 ` column i+2: derivative with respect to width=`------------------------------------------------------------>` modelf: computes the model function but without derivatives.&` Used for displaying the current fit.&` Should be re-defined for each model. ` Input: x,a ` Output: yfunc modelf(x,a,y)na = length(a)clear y;for i = 1 to na by 3 y += a[i]*exp(-((x-a[i+1])/a[i+2])^2)=`------------------------------------------------------------)` getdata: defines the experimental data.<` Normally data is read from a file but here we simulate it.*` The x's are defined in the range [0,10].3` The y's are defined as a sum of Gaussians + noisefunc getdata(x,y,sig)(ndata = 200 ` number of data points` the simulated noise3randomize ` initialize the random number generatorEn[ndata]:1 ` create array of 'ndata' points (max. amplitude of noise)0n = rand(n) - 0.5 ` create random noise around 0 ` the x's2x = 0 to 10 len ndata ` ndata points from 0 to 10 ` the y's)yt[ndata]:0 ` array for theoretical y's$y[ndata]:0 ` array for actual y'sA = 1.2,1.8,1.5 ` amplitudesP = 2,5,8 ` positionsW = 1.2,0.7,1.6 ` widthsna = length(A)Afor k = 1 to na yt += A[k]*exp(-((x-P[k])/W[k])^2) ` exact valuesy = yt + n ` plus noise` the standard deviationsig[ndata]:1 ` uniform array=`------------------------------------------------------------F` mrqcof: computes the matrix 'alpha', the vector 'beta', and 'chisq'.E` In this implementation of 'mrqcof' we changed the order of loops soC` that the loop over the data points is moved inside and is carried'` out by array operations and by 'sum'.` Input: x,y,sig,a` Output: alpha,beta,chisq'func mrqcof(x,y,sig,a,alpha,beta,chisq)ndata = length(x)+ma = length(a) ` number of parameters3clear alpha,beta ` initialize all elements to 0 ymod[ndata]:0dyda[ndata,ma]:02` call the model function once for all data pointsGmodel(x,a,ymod,dyda) ` calculate the model function and its derivatives!chisq = 0 ` initializesig2 = 1/(sig*sig) dy = y - ymod` define 'alpha' and 'beta'9` loop over the elements and sum over all the data points` for each elementfor j = 1 to mawt = dyda[*,j]*sig2.for k = 1 to j alpha[j,k] = sum(wt*dyda[*,k])beta[j] += sum(dy*wt)(` define the symmetric elements of alpha:for j = 2 to ma for k = 1 to j-1 alpha[k,j] = alpha[j,k]` define 'chisq'chisq = sum(dy*dy*sig2)=`------------------------------------------------------------,` mrqmin: the Levenberg-Marquardt algorithm.:` This is more or less similar to 'mrqmin' in NR with some` modifications.(` Input: x,y,sig,a,alpha,beta,chisq,lam ` Output: a,alpha,beta,chisq,lam+func mrqmin(x,y,sig,a,alpha,beta,chisq,lam)B` define linear system of equations: [new_alpha]*new_beta = betanew_alpha = alphama = length(a)Lj = 1 to ma ` index array to address diagonal elements<new_alpha[j,j] *= (1+lam) ` modify diagonal elementsE` Note: before we call 'linsol' to solve the linear system we disableE` Numerit's automatic error handling by the 'HideError' command. ThisC` allows us to catch the 'Singular Matrix' error (24) and handle itK` locally by restarting the process from the initial guess with a differentE` 'lam' parameter (the automatic error handler pauses the program and` shows the error). HideError;new_beta = linsol(new_alpha,beta) ` solve the linear system!if errval = 24 ` singular matrix"ResetError ` reset error flag/ShowError ` resume automatic error handling;beep ` make a beep to let us know that it happened%` initialize process and increase lam$init(x,y,sig,a,alpha,beta,chisq,lam) randomize=lam *= 3+rand(5) ` increase lam (with some randomization) return false8ShowError ` resume automatic error handling0new_a = a + new_beta ` trial set of parametersnew_chisq = chisqUmrqcof(x,y,sig,new_a,new_alpha,new_beta,new_chisq) ` calculate new alpha, beta, chisq!if new_chisq <= chisq ` improved randomize=lam *= 0.2+rand(0.2) ` decrease lam (with some randomization) a = new_aalpha = new_alphabeta = new_betachisq = new_chisq return true` worse randomize=lam *= 3+rand(5) ` increase lam (with some randomization) return false=`------------------------------------------------------------` init: Set initial values.?` Here we define the initial guess for the parameters array and4` also initialize 'alpha','beta','chisq', and 'lam'.B` In this case we use a simple initial guess (all parameters = 1);(` normally we would have a better guess.)func init(x,y,sig,a,alpha,beta,chisq,lam) a[*] = 1 chisq = 0 lam = 0.01#mrqcof(x,y,sig,a,alpha,beta,chisq)=`------------------------------------------------------------` The main program*ma = 9 ` number of parameters1a[ma]:0 ` create the parameters array1alpha[ma,ma]:0 ` create the curvature matrix,beta[ma]:0 ` create the beta vector chisq = 0lam = 0` start the process:x = 0; y = 0; sig = 0 ` are defined as arrays by 'getdata'getdata(x,y,sig)9yfit = y ` initialize array of fitted values+` initialize a, alpha, beta, chisq, and lam$init(x,y,sig,a,alpha,beta,chisq,lam)ochisq = chisqK` The following loop calls 'mrqmin' and checks for a termination condition.N` The termination condition used here may not be universal. It works fine withO` the current model function but we may have to use a different condition for a` different model.loopAif mrqmin(x,y,sig,a,alpha,beta,chisq,lam) ` true = improvementD` define array of fitted values with the current set of parameters` and display the current yfitmodelf(x,a,yfit)refresh ` update viewers` check termination conditionif chisq = 0 break8if chisq <= ochisq if (ochisq-chisq)/chisq < 1e-40 breakochisq = chisq1` end of process: calculate the covariance matrix HideErrorcovar = inv(alpha)8if errval <> 24 ` calculate only if not singularFShowError ` resume automatic error handling and continueKj = 1 to ma ` array index - used to extract the diagonal elements-vari = covar[j,j] ` the variance elements)meanv = mean(vari) ` the mean variance7stdv = stdev(vari) ` standard deviation of variance:ShowError ` show the 'Singular Matrix' errord:\num\num1.5\work\ maaalphabetachisq lamx y sigyfit ochisqcovarjvarimeanv stdv     4 , 4 , 4 ,        <?          <?          <7   <?  \76  [7 A C Y76  6   ]7 e  4    d8bQ__ model modelxaydydanaiaiai1ai2arexfac a*+ , -   X80 41  @42  @4 3  A C 4 G B 5  B B B 7  BN8 4 9  @4 C:  @4 B C N6<> \_ modelfmodelfxaynai ,DE F G   X8 4   @4A  @4C EGBN N6H> '* getdata getdataxysigndatanyt APWnak aOP ST 4 ,U  AX   f[ 4 ,\ 4 ,]   5^  5_  5 `  a  X8 4  4A 4 C EGBN N 6b  @e 4 ,f>;RU mrqcofmrqcofxysigaalphabetachisqndatamaymoddydasig2dyjwtk op q  r  s 4 , t  4 , v   <?w x  BC y  A ~   X8 4 B   X8  4  4 B N6 4 BS N6   X8    AX8  4  4 N6  N6 B B> 2]`AQTc|ix{ mrqminmrqminxysigaalphabetachisqlam new_alphamaj new_betanew_a new_chisq  |     e  4 @U    \7         <?  @P=  @      <?  [7  @P    =  @P=>?o init initxysigaalphabetachisqlam 4          <?> 1?3@2@200i@0.5?010$@1.2333333?1.8?1.5?5@8 @0.7ffffff?1.6?248@0.2?0.01{Gz?9"@1e-40Ww'&l7