'*******************************************************
'*    Program to demonstrate the Parafit subroutine    *
'* --------------------------------------------------- *
'*  Reference: BASIC Scientific Subroutines, Vol. II   *
'*    By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981     *
'*    [BIBLI 01].                                      *                                
'* --------------------------------------------------- *
'* SAMPLE RUN:                                         *
'*                                                     *
'* Parametric least squares fit                        *
'*                                                     *
'* The input data are:                                 *
'*                                                     *
'* X( 1) =  1.000000   Y( 1) =  0.033702               *
'* X( 2) =  2.000000   Y( 2) =  0.249029               *
'* X( 3) =  3.000000   Y( 3) =  0.944733               *
'* X( 4) =  4.000000   Y( 4) =  1.840089               *
'* X( 5) =  5.000000   Y( 5) =  1.840089               *
'* X( 6) =  6.000000   Y( 6) =  0.944733               *
'* X( 7) =  7.000000   Y( 7) =  0.249029               *
'* X( 8) =  8.000000   Y( 8) =  0.033702               *
'* X( 9) =  9.000000   Y( 9) =  0.002342               *
'* X(10) = 10.000000   Y(10) =  0.000084               *
'*                                                     *
'* The coefficients are:                               *
'*  2.000000                                           *
'*  4.500000                                           *
'*  3.000000                                           *
'*                                                     *
'* The standard deviation of the fit is 0.000000       *
'*                                                     *
'* The number of iterations was 51                     *
'*******************************************************
DEFINT I-N
DEFDBL A-H, O-Z
DIM A(25), e1(25), X(25), Y(25)
n = 10: l = 3
PRINT " Parametric least squares fit"
PRINT
F$=" X(##) = ##.######   Y(##) = ##.######"
PRINT " The input data are:"
PRINT
FOR i = 1 TO n
  X(i) = i
  Y(i) = 2# * EXP(-(X(i) - 4.5#) * (X(i) - 4.5#) / 3#)
  PRINT using F$; i; X(i); i; Y(i)
NEXT i
PRINT
e = .1: e1 = .5: A(1) = 10#: A(2) = 10#: A(3) = 10#
GOSUB 1000
F$="  ##.######"
PRINT " The coefficients are:"
PRINT using F$; A(1)
PRINT using F$; A(2)
PRINT using F$; A(3)
PRINT
PRINT using " The standard deviation of the fit is ##.########"; d
PRINT
PRINT " The number of iterations was "; m
PRINT
END

'***************************************************************
'*      Parametric least squares curve fit subroutine          *
'* ----------------------------------------------------------- *
'* This program least squares fits a function to a set of data *
'* values by successively reducing the variance.  Convergence  *
'* depends on the initial values and is not assured.           *
'* n pairs of data values, X(i), Y(i), are given. There are l  *
'* parameters, A(j), to be optimized across.                   *
'* Required are initial values for the A(l) and e. Another     *
'* important parameter which affects stability is e1, which is *
'* initially converted to e1(l)for the first intervals.        *
'* The parameters are multiplied by (1 - e1(i)) on each pass.  *
'***************************************************************
1000 FOR i = 1 TO l
  e1(i) = e1
NEXT
m = 0
'Set up test residual
xl1 = 1000000#
'Make sweep through all parameters
1050 FOR i = 1 TO l
  a0 = A(i)
  'Get value of residual
  A(i) = a0
1100  GOSUB 1200
  'Store result in xm0
  xm0 = xl2
  'Repeat for xm1
  A(i) = a0 * (1# - e1(i))
  GOSUB 1200
  xm1 = xl2
  'Change interval size if called for
  'If variance was increased, halve E1(i)
  IF xm1 > xm0 THEN e1(i) = -e1(i) / 2#
  'If variance was reduced, increase step size by increasing E1(i)
  IF xm1 < xm0 THEN e1(i) = 1.2# * e1(i)
  'If variance was increased, try to reduce it
  IF xm1 > xm0 THEN A(i) = a0
  IF xm1 > xm0 THEN GOTO 1100
NEXT i
'End of a complete pass
'Test for convergence
m = m + 1
IF xl2 < 0 THEN RETURN
IF ABS((xl1 - xl2) / xl2) < e THEN RETURN
'If this point is reached, another pass is called for
xl1 = xl2
GOTO 1050
'Residual generation subroutine
1200 xl2 = 0
FOR j = 1 TO n
  X = X(j)
  'Obtain function
  GOSUB 1500
  xl2 = xl2 + (Y(j) - Y) * (Y(j) - Y)
NEXT j
d = SQR(xl2 / (n - l))
RETURN

1500 'Function subroutine
  Y = A(1) * EXP(-(X - A(2)) * (X - A(2)) / A(3))
RETURN

' End of file parafit.bas