home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Interactive Guide / c-cplusplus-interactive-guide.iso / c_ref / csource4 / 209_01 / simplexr.doc < prev    next >
Encoding:
Text File  |  1990-03-05  |  9.8 KB  |  335 lines

  1. SIMPLEXR.DOC       VERS:- 01.00  DATE:- 09/26/86  TIME:- 10:02:52 PM
  2.  
  3. notes on the simplex method of function minimization, 
  4.     by use of the Nelder-Mead algorithm.
  5.  
  6.  
  7.  
  8.  
  9.  
  10.  
  11.  
  12.         By J. A. Rupley, Tucson, Arizona
  13.  
  14.  
  15.          NOTES ON THE ALGORITHM OF THE SIMPLEX FITTING PROGRAM
  16.  
  17.  
  18.         A. THE SIMPLEX METHOD OF FUNCTION MINIMIZATION
  19.  
  20.  
  21.              The fitting of a model with variable parameters to a set of
  22.         data points is handled by minimizing an error function, such as
  23.         the sum of squares of the differences between observed values and
  24.         values calculated according to the model.
  25.  
  26.              We describe here a simplex approximation procedure for
  27.         estimating the parameter values that give a function minimum.
  28.         The strategy is that of Nelder and Mead, and their article
  29.         (Computer Journal, vol 7, p 308, 1965) should be consulted for
  30.         more information than is given below.
  31.  
  32.              A simplex is constructed in parameter space, with vertices
  33.         that are arbitrary but not too far from the point at which the
  34.         function is a minimum and that describe a volume that includes
  35.         the minimum point.
  36.  
  37.              The position of each vertex in parameter space defines (or
  38.         equally, is defined by) a set of parameter values, for which one
  39.         calculates the function value for that vertex.  The simplex has
  40.         (M + 1) vertices, where M is the number of free parameters.
  41.  
  42.              In successive iterations the vertex with the highest
  43.         function value is moved, to obtain a new vertex position of
  44.         smaller function value.  The movement is directed with respect to
  45.         the center of the reduced simplex, which is the simplex less the
  46.         highest vertex.  The movements can be the following:
  47.              reflection;
  48.              expansion beyond the reflected position;
  49.              reflection after failed expansion;
  50.              contraction toward the center of the reduced simplex from
  51.                   the original position;
  52.              contraction from the reflected position;
  53.              if none of these operations gives a lower function value
  54.                   then the entire simplex is shrunk about the lowest
  55.                   vertex.
  56.  
  57.              Exit from the iterative minimization is if the fractional
  58.         RMS deviation of the function values at the vertices is less than
  59.         a test value and if the centroid of the simplex is within two RMS
  60.         deviations of the mean of the vertices.  The default test value
  61.  
  62.         is 1E-8 * the mean function value.
  63.  
  64.              Thus at exit from the minimization, the centroid of the
  65.         simplex gives an arbitrarily close approximation of the parameter
  66.         values at the function minimum.
  67.  
  68.  
  69.                                         1
  70.  
  71.  
  72.  
  73.  
  74.  
  75.  
  76.  
  77.  
  78.  
  79.  
  80.              The routine that calculates the function values must be
  81.         adapted to the model to be optimized.
  82.  
  83.              Entry of the data and starting parameter values can be at
  84.         compilation, by initializing the data and parameter arrays, or by
  85.         reading a text file.
  86.  
  87.              Results should be written to a file, for convenience of
  88.         review and for recovery from crashes.
  89.  
  90.  
  91.  
  92.  
  93.  
  94.  
  95.  
  96.  
  97.  
  98.  
  99.  
  100.  
  101.  
  102.  
  103.  
  104.  
  105.  
  106.  
  107.  
  108.  
  109.  
  110.  
  111.  
  112.  
  113.  
  114.  
  115.  
  116.  
  117.  
  118.  
  119.  
  120.  
  121.  
  122.  
  123.  
  124.  
  125.  
  126.  
  127.  
  128.  
  129.  
  130.  
  131.  
  132.  
  133.  
  134.  
  135.  
  136.                                         2
  137.  
  138.  
  139.  
  140.  
  141.  
  142.  
  143.  
  144.  
  145.  
  146.  
  147.         B. THE QUADRATIC FIT TO THE LEAST SQUARES FUNCTION SURFACE
  148.  
  149.  
  150.              Standard deviations of the parameters are calculated by
  151.         fitting a quadratic function to the surface about the minimum in
  152.         parameter space and then using the properties of this function to
  153.         calculate the variance-covariance matrix for the parameters.  The
  154.         strategy is a modification of that of Nelder and Mead (1965).
  155.  
  156.              The contracted simplex obtained through the previous
  157.         minimization process is used to define a new coordinate system of
  158.         M oblique axes in parameter space.  The origin of the axes is at
  159.         the centroid of the simplex.  Each axis corresponds to a free
  160.         parameter.  The unit value (scale) for each axis is set in the
  161.         new system at the average of the deviations of the vertices of
  162.         the simplex from the centroid value.
  163.  
  164.              In effect, one constructs in the new coordinate system a new
  165.         simplex, based on the one obtained by the previous minimization
  166.         process.  The (M + 1) vertices of the new construct are at the
  167.         centroid (the zero of the new system) and at the unit points on
  168.         each free parameter axis.  As before, the values of the least
  169.         squares function at these vertices describe a surface.  It is
  170.         this least squares surface to which a quadratic function is to be
  171.         fit.
  172.  
  173.  
  174.              If necessary, the scales of the axes (the unit values) in
  175.         the new coordinate system can be adjusted to give better behavior
  176.         of the surface in the vicinity of the minimum.
  177.  
  178.              The diagonal matrix Q transforms the new coordinate system
  179.         back to the original system.  The elements of Q are for each
  180.         parameter the (adjusted) average differences between the vertices
  181.         and the centroid of the original simplex.
  182.  
  183.              The function value, y, in the region near the function
  184.         minimum generally can be approximated in the new coordinate
  185.         system by the quadratic vector function:
  186.  
  187.               y  =  a(o)   +   2 a'.x   +   x'.B.x                   (1)
  188.  
  189.              The elements of the vector x are the values of the M free
  190.         parameters at a point in parameter space.  The scalar a(o), the
  191.         vector a and the matrix B are determined by the shape of the
  192.         surface y near the minimum and can be calculated with the use of
  193.         simple numerical approximations, evaluated at x = 0 (the position
  194.         of the centroid of the original simplex):
  195.  
  196.  
  197.  
  198.  
  199.  
  200.  
  201.  
  202.  
  203.                                         3
  204.  
  205.  
  206.  
  207.  
  208.  
  209.  
  210.  
  211.  
  212.  
  213.                 a(0)   =   y(0)
  214.  
  215.                 ai     =   0.5 * (dy/dxi)
  216.                        =   0.25 * (yi - y-i)
  217.  
  218.                 bij    =   0.5 * (d2y/dxi.dxj)
  219.                        =   0.25 * (yij + y-i-j + 2 * y(0)
  220.                                            - y-i - y-j - yi - yj)
  221.  
  222.                 bii    =   0.5 * (d2y/dxi2)
  223.                        =   0.5 * (yi + y-i - 2 * y(0))
  224.  
  225.              The point x = 0 (the centroid position) may not be the true
  226.         position of the function minimum, although it is assumed to be
  227.         close to it.  A refined estimate of the minimum position, based
  228.         on the quadratic approximation of equation (1) above, is given by
  229.  
  230.         x(min):
  231.  
  232.              xmin = - B-1.a
  233.  
  234.  
  235.              The function minimum at x(min) is given by y(min) :
  236.  
  237.              ymin = a(0) - a'.B-1.a
  238.  
  239.  
  240.              The position of the minimum in the original coordinate
  241.         system is obtained by back transformation of x(min) with the
  242.         matrix Q, giving the point p(min):
  243.  
  244.              pmin = p(0) - Q'.B-1.a
  245.  
  246.              As a test of the quality of the quadratic approximation (1),
  247.         there should be close agreement between the function values
  248.         y(min), y(p(min)), and y(centroid), and between the sets of
  249.         parameter values defining the points p(min) and the centroid.
  250.  
  251.              The variance-covariance matrix C is given by
  252.  
  253.              C = Q'.B-1.Q
  254.  
  255.              A diagonal element of C, multiplied by the mean square error
  256.         (MSE), gives the variance of the corresponding parameter:
  257.  
  258.              var(i) = cii * MSE
  259.  
  260.              MSE = ymin/(n - m)
  261.  
  262.         where the divisor (n - m) = # degrees freedom: n = # of data
  263.         points used in calcn of ymin, and m = # of free parameters; and
  264.         ymin = the sum of residuals squared
  265.  
  266.  
  267.  
  268.  
  269.  
  270.                                         4
  271.  
  272.  
  273.  
  274.  
  275.  
  276.  
  277.  
  278.  
  279.  
  280.              The standard deviation is the square root of the variance:
  281.  
  282.              std-dev(i) = sqrt(var(i))
  283.  
  284.              The quadratic fit can fail, giving negative values of the
  285.  
  286.         diagonal elements of the variance-covariance matrix or a value of
  287.         y(p(min)) > y(centroid).  Failure can come from too tight a
  288.         simplex, resulting in too small values of the B matrix elements,
  289.         or from too large a simplex, resulting in non-quadratic behavior
  290.         of the least squares function within the region defined by the
  291.         simplex.  The problem of too tight a simplex does not arise if
  292.         the precision of the floating point routines of the program is
  293.         sufficiently high (e.g., 14 digits).  The problem of non-quadratic
  294.         behavior commonly arises in the early stages of the fitting, when
  295.         the minimum has not yet been approached closely.  This can be
  296.         taken into account by modifying the fitting algorithm, or one can
  297.         be patient.  If a parameter moves so close to a bound that
  298.         expansion of the simplex toward the bound is not possible, then
  299.         that parameter should be fixed.
  300.  
  301.              It may be useful to cycle between short sessions of simplex
  302.         fitting and quadratic approximation, to more rapidly approach the
  303.         lsq function minimum.  Toward this end, the presumed more
  304.         accurate approximation of the function minimum given by p(min) is
  305.         used to replace the previously obtained centroid value of the
  306.         simplex.  That is, the simplex returned to the simplex fitting
  307.         routine would have M vertices at the "unit" positions on the
  308.         parameter axes and the (M+1)th vertex at p(min), not at the
  309.         "origin", ie the old centroid value.
  310.  
  311.  
  312.  
  313.  
  314.  
  315.  
  316.  
  317.  
  318.  
  319.  
  320.  
  321.  
  322.  
  323.  
  324.  
  325.  
  326.  
  327.  
  328.  
  329.  
  330.  
  331.  
  332.  
  333.  
  334.  
  335.  
  336.  
  337.                                         5
  338.  
  339.  
  340.  
  341.  
  342.