Demos
 Demos
Source Code of nlfit.num Document
Document 
 
` Nonlinear model fitting: the Levenberg-Marquardt method.
` This is an implementation of the algorithm as described
` in Numerical Recipes 15.5. The main modification is
` converting loops over data points to array operations.
` The algorithm is demonstrated using a simulated signal which
` is created as a sum of three Gaussians and a random noise.
` On each run, a new set of points is created.
`*********************************************************

` Functions

` model: computes the model function and its derivatives with
` respect to the parameters.
` Should be re-defined for each model.
` In this case: A sum of Gaussians each having a variable
` position, amplitude, and width (after example 15.5.16 NR, 'fgauss').
` Note that since we call the model function once for all the data
` points 'dyda' is a 2-dim array with dimensions [ndata x na] rather
` than a 1-dim array with 'na' elements. Similarly, 'y' is an array
` of size 'ndata' and not a scalar.
` Input: x,a
` Output: y,dyda

func 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
 ai = a[i] ` amplitude
 ai1 = a[i+1] ` position
 ai2 = a[i+2] ` width
 ar = (x-ai1)/ai2
 ex = exp(-ar*ar)
 fac = 2*ai*ex*ar
 ` the function and its derivatives
 y += ai*ex ` sum of Gaussians
 dyda[*,i] = ex ` column i : derivative - amplitude
 dyda[*,i+1] = fac/ai2 ` column i+1: derivative - position
 dyda[*,i+2] = fac*ar/ai2 ` column i+2: derivative - 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: y

func 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].
` The y's are defined as a sum of Gaussians + noise

func getdata(x,y,sig)
   ndata = 200 ` number of data points

 ` the simulated noise
 randomize ` initialize the random number generator
 n[ndata]:1 ` create array of 'ndata' points and set all
            ` elements to 1 (max. amplitude of noise)
 n = rand(n) - 0.5 ` create random noise around 0

 ` the x's
 x = 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's
 A = 1.2,1.8,1.5 ` amplitudes
 P = 2,5,8 ` positions
 W = 1.2,0.7,1.6 ` widths
 na = length(A)
 for k = 1 to na yt += A[k]*exp(-((x-P[k])/W[k])^2) ` exact values
 y = yt + n ` plus noise

 ` the standard deviation
 sig[ndata]:1 ` create array with uniform weight

`------------------------------------------------------------
` mrqcof: computes the matrix 'alpha', the vector 'beta', and 'chisq'.
` In this implementation of 'mrqcof' we changed the order of loops so
` 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 parameters
 clear alpha,beta ` initialize all elements to 0
 ymod[ndata]:0
 dyda[ndata,ma]:0

 ` call the model function once for all data points
 model(x,a,ymod,dyda) ` calculate the model function and its derivatives
 chisq = 0 ` initialize
 sig2 = 1/(sig*sig)
 dy = y - ymod

 ` define 'alpha' and 'beta'
 ` loop over the elements and sum over all the data points
 ` for each element

 for j = 1 to ma
    wt = 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)
   ` define linear system of equations: [new_alpha]*new_beta = beta
 new_alpha = alpha
 ma = length(a)
 j = 1 to ma` index of diagonal elements
 new_alpha[j,j] *= (1+lam)` modify diagonal elements
 ` Note: before we call 'linsol' to solve the linear system we disable
 ` Numerit's automatic error handling by the 'HideError' command. This
 ` allows us to catch the 'Singular Matrix' error (24) and handle it
 ` locally by restarting the process from the initial guess with a
 ` different '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 false
 ShowError` resume automatic error handling
 new_a = a + new_beta ` trial set of parameters
 new_chisq = chisq

 ` calculate new alpha, beta, chisq
 mrqcof(x,y,sig,new_a,new_alpha,new_beta,new_chisq)
 if new_chisq <= chisq ` improved
    randomize
 lam *= 0.2+rand(0.2) ` decrease lam (with some randomization)
 a = new_a
 alpha = new_alpha
 beta = new_beta
 chisq = 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 and
` also initialize 'alpha','beta','chisq', and 'lam'.
` 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 ` set all paremeters to 1
 chisq = 0
 lam = 0.01
 mrqcof(x,y,sig,a,alpha,beta,chisq)

`------------------------------------------------------------
` The main program


ma = 9  ` number of parameters
a[ma]:0 ` create the parameters array
alpha[ma,ma]:0 ` create the curvature matrix
beta[ma]:0 ` create the beta vector
chisq = 0
lam = 0

` start the process
x = 0; y = 0; sig = 0 ` are defined as arrays by 'getdata'
getdata(x,y,sig)
yfit = y ` initialize array of fitted values

` initialize a, alpha, beta, chisq, and lam
init(x,y,sig,a,alpha,beta,chisq,lam)
ochisq = chisq

` The following loop calls 'mrqmin' and checks for a termination condition
loop
   if mrqmin(x,y,sig,a,alpha,beta,chisq,lam) ` true = improvement
    ` define array of fitted values with the current set of parameters
 ` and display the current yfit

 modelf(x,a,yfit)
 refresh ` update viewers

 ` check termination condition
 if chisq = 0 break
 if chisq <= ochisq if (ochisq-chisq)/chisq < 1e-40 break
 ochisq = chisq

` end of process: calculate the covariance matrix
covar = inv(alpha)
j = 1 to ma ` index of diagonal elements
vari = covar[j,j] ` the variance elements
meanv = mean(vari) ` the mean variance
stdv = stdev(vari) ` standard deviation of variance
 
 Demos
Demos
Copyright © 2004 KEDMI Scientific Computing. All Rights Reserved. Document 
Document