home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / e / e032 / 3.ddi / FILES / STATISTI.PAK / NONLINEA.M < prev    next >
Encoding:
Text File  |  1992-07-29  |  14.5 KB  |  441 lines
  1. (* :Title: Nonlinear Curve Fitting *)
  2.  
  3. (* :Context: Statistics`NonlinearFit` *)
  4.  
  5. (* :Author: John M. Novak *)
  6.  
  7. (* :Summary: A package to perform nonlinear curve fitting,
  8.     using a variety of algorithms.
  9. *)
  10.  
  11. (* :Package Version: 1.1 *)
  12.  
  13. (* :Mathematica Version: 2.0 *)
  14.  
  15. (* :History:
  16.     Version 1.0, Oct. 91 by John M. Novak
  17.     Version 1.1, Feb. 92 by John M. Novak - major revisions
  18. *)
  19.  
  20. (* :Keywords:
  21.     curve fitting, nonlinear regression
  22. *)
  23.  
  24. (* :Sources:
  25.     Janhunen, Pekka, NonlinearFit`, (a Mathematica package), 1990.
  26.     Press, William, et. al., "Numerical Recipies in Pascal",
  27.         pp. 572-580, Cambridge University Press (Cambridge, 1989).
  28.     Withoff, Dave, DataAnalysis`NonlinearRegression`,
  29.         (a Mathematica package), 1989.
  30. *)
  31.  
  32. (* :Discussion:
  33.     
  34. *)
  35.  
  36. BeginPackage["Statistics`NonlinearFit`",
  37.     "Utilities`FilterOptions`"]
  38.  
  39. NonlinearFit::usage =
  40. "NonlinearFit[data,model,vars,params,(opts)] uses iterative
  41. methods to fit the data to the model containing the given
  42. variables and parameters.  Paramaters may be expressed as
  43. a list of symbols or {symbol, startingvalue} pairs.  The
  44. data are given as {{x1, y1, ..., f1},{x2, y2, ... f2},...};
  45. some variations on this are also valid (see documentation.)";
  46.  
  47. Weights::usage =
  48. "Weights is an option for NonlinearFit.  One can pass in an array
  49. of weights, one for each point in the data set.";
  50.  
  51. ShowProgress::usage =
  52. "ShowProgress is an option for NonlinearFit.  If True, it will
  53. print a summary of intermediate results for each iteration
  54. of the algorithm.";
  55.  
  56. LevenbergMarquardt::usage =
  57. "LevenbergMarquardt is a value for the Method option of NonlinearFit.
  58. This causes NonlinearFit to use the Levenberg-Marquardt method of
  59. minimizing chi^2."
  60.  
  61. Begin["`Private`"]
  62.  
  63. (* first, a quick private utility... *)
  64.  
  65. numberQ[n_] := NumberQ[N[n]]
  66.  
  67. (* NonlinearFit options and messages *)
  68.  
  69. Options[NonlinearFit] =
  70.     {MaxIterations -> 30,
  71.     Method -> LevenbergMarquardt,
  72.     ShowProgress -> False,
  73.     Weights -> Automatic,
  74.     WorkingPrecision -> $MachinePrecision,
  75.     AccuracyGoal -> 1,
  76.     PrecisionGoal -> 3};
  77.  
  78. NonlinearFit::badits =
  79. "MaxIterations -> `1` is not set to a positive integer;
  80. setting to 30.";
  81.  
  82. NonlinearFit::badmeth =
  83. "Method `1` is not available for NonlinearFit;  using
  84. SteepestDescent method.";
  85.  
  86. NonlinearFit::toomany =
  87. "Too many variables have been provided for the size of the data;
  88. given `1`-tuples for data, given `2` variables.";
  89.  
  90. NonlinearFit::badweights =
  91. "There is not a valid weight for each data point!
  92. Setting all weights to 1.";
  93.  
  94. NonlinearFit::badparams =
  95. "The form of the parameters `1` given to NonlinearFit is not correct.
  96. Parameters should be of form {parameter}, {parameter, startpoint}, or
  97. {parameter, startpoint, minrange, maxrange}.";
  98.  
  99. NonlinearFit::badwork =
  100. "The value `1` given for the WorkingPrecision is not a positive
  101. integer.  Setting to $MachinePrecision...";
  102.  
  103. NonlinearFit::badprec =
  104. "The value `1` given for the PrecisionGoal is not a positive integer
  105. or Automatic.  Setting to $MachinePrecision - 10...";
  106.  
  107. NonlinearFit::badacc =
  108. "The value `1` given for the AccuracyGoal is not a positive integer
  109. or Automatic.  Setting to $MachinePrecision - 10...";
  110.  
  111. NonlinearFit::degrees =
  112. "Warning: The data set is smaller than the number of parameters; this is likely
  113. to generate odd answers...";
  114.  
  115. NonlinearFit::lmnoconv =
  116. "Warning: NonlinearFit did not converge in `1` iterations. The
  117. returned value may not be at the minimum.";
  118.  
  119. NonlinearFit::lmoutofrng =
  120. "Warning: the LevenbergMarquardt search exited the range given to the
  121. parameters.  The returned value may not be at the minimum.";
  122.  
  123. NonlinearFit::lmpnocon =
  124. "Warning: the values of the parameters given to NonlinearFit do not
  125. appear to have converged.  The returned value may not be at the
  126. minimum.";
  127.  
  128. NonlinearFit::lmovfl =
  129. "Warning: numeric overflow occurred during this search.  This may mean
  130. that the search is being started from an inappropriate point, that your
  131. model is unstable, or that insufficient precision is being used for these
  132. calculations.  The returned value may not be at the minimum.";
  133.  
  134. NonlinearFit::fitfail =
  135. "The fitting algorithm failed utterly, returning `1`.";
  136.  
  137. (* NonlinearFit arguments:
  138.     data - vector or matrix of numbers;
  139.     model - the model, a symbolic expression;
  140.     vars - a list of symbols, expressing the variables in model;
  141.     inputparams - a list of symbol and {symbol, startingpoint}
  142.         pairs, expressing the parameters and their starting
  143.         points in the model;
  144.     opts - the options
  145. *)
  146.  
  147. NonlinearFit[d_,m_,var_Symbol,params_,opts___] :=
  148.     NonlinearFit[d,m,{var},params,opts]
  149.  
  150. NonlinearFit[d_,m_,v_,
  151.         param:(_Symbol |
  152.             {_Symbol,(_?numberQ | Infinity | -Infinity)..}
  153.             ),
  154.         opts___] :=
  155.     NonlinearFit[d,m,v,{param},opts]
  156.  
  157. NonlinearFit[data:{{__?numberQ, {__?numberQ}..}..},
  158.         m_,vars_,params_,opts___] :=
  159.     NonlinearFit[
  160.         Map[If[Head[#] === List, First[#], #]&,data,{2}],
  161.         m,vars,params,opts,Weights -> 
  162.             Map[If[Head[#] === List, Last[#], 1]&,
  163.                 Map[Drop[#,Length[vars]]&,data]
  164.             ]
  165.     ]
  166.  
  167. NonlinearFit[data_?(VectorQ[N[#],NumberQ]&),
  168.         m_,vars_,params_,opts___] :=
  169.     NonlinearFit[Transpose[Append[(Evaluate[Table[#,
  170.                 {Length[vars]}]]&)[Range[Length[data]]],data]],
  171.             m,vars,params,opts]
  172.  
  173. (* The following routine will catch whether the fit could be
  174. performed;  if not, it will allow the function to return
  175. unevaluated, rather than returning a meaningless error value.
  176. (Instead, it returns a meaningless unevaluated expression...:-)
  177. *)
  178.  
  179. NonlinearFit[data:_?(MatrixQ[N[#],NumberQ]&),
  180.         model_,vars:{__Symbol},
  181.         inputparams:{(_Symbol |
  182.             {_Symbol,(_?numberQ | Infinity | -Infinity)..})..},
  183.         opts:((_Rule | RuleDelayed)...)] :=
  184.     With[{fitted =
  185.             nonlinearfit[data,model,vars,inputparams,opts]},
  186.         fitted/;fitted =!= $Failed
  187.     ]
  188.  
  189. nonlinearfit[data_, model_, vars_, inputparams_, opts___] :=
  190.     Module[{maxits,method,progress,weights,params,workprec,accgoal,
  191.             pts,datalen,varlen = Length[vars],minrng,maxrng,min,max,x,
  192.             fmopts,out,start,precgoal,response,s},
  193.     (* Get options *)
  194.         {maxits,method,progress,weights,workprec,accgoal,precgoal} =
  195.             {MaxIterations, Method, ShowProgress, Weights,
  196.                 WorkingPrecision, AccuracyGoal, PrecisionGoal}/.
  197.             {opts}/.Options[NonlinearFit];
  198.         fmopts = FilterOptions[FindMinimum,opts];
  199.     (* Option and Argument checking *)
  200.     (* - check input data, separate out response and coords *)
  201.         If[(datalen = Length[First[data]]) - varlen < 1,
  202.                 Message[NonlinearFit::toomany,datalen,varlen];
  203.                 Return[$Failed]];
  204.         response = Transpose[Map[Drop[#,varlen]&, data]];
  205.         pts = Map[Drop[#,varlen - datalen]&, data];
  206.     (* - set and check weights *)
  207.         If[weights === Automatic,
  208.             weights = Table[1,{#1},{#2}]& @@ Dimensions[response]];
  209.         If[VectorQ[N[weights],NumberQ],
  210.             weights = Table[weights,{Length[response]}]];
  211.         If[Head[weights =!= List],
  212.             weights = Map[weights, response, {2}]];
  213.         If[!MatrixQ[N[weights],NumberQ] ||
  214.                 Dimensions[weights] =!= Dimensions[response],
  215.             Message[NonlinearFit::badweights];
  216.             weights = Table[1,{#1},{#2}]& @@ Dimensions[response]
  217.         ];
  218.     (* - check parameters, start points, ranges *)
  219.         {params,start,minrng,maxrng} = Transpose[
  220.             Map[Replace[#,
  221.                 {(x_Symbol | {x_Symbol}) -> {x,{1},-Infinity, Infinity},
  222.                 {x_Symbol, s_?numberQ} -> {x,{s},-Infinity,Infinity},
  223.                 {x_Symbol, min_?numberQ, max_?numberQ} ->
  224.                     {x,startposns[min,max],min,max},
  225.                 {x_Symbol, min_?numberQ, Infinity} ->
  226.                     {x,{min + 1, min}, Infinity},
  227.                 {x_Symbol, -Infinity, max_?numberQ} ->
  228.                     {x, {max - 1}, -Infinity, max},
  229.                 {x_Symbol, -Infinity, Infinity} ->
  230.                     {x,{1}, -Infinity, Infinity},
  231.                 {x_Symbol, s_?numberQ, min:(-Infinity | _?numberQ),
  232.                         max:(_?numberQ | Infinity)} ->
  233.                     {x,s,min,max},
  234.                 _ -> {$Failed,$Failed,$Failed,$Failed}}]&,
  235.             inputparams]
  236.         ];
  237.         If[(Select[start,(# === $Failed)&] =!= {}) ||
  238.                 Not[And @@ Thread[N[minrng <= start <= maxrng]]],
  239.             Message[NonlinearFit::badparams, inputparams];
  240.             Return[$Failed]
  241.         ];
  242.         start = Flatten[Outer[List,##]& @@ start,Length[start] - 1];
  243.         If[Length[start] === 1,
  244.             start = First[start],
  245.             start = findbeststart[start, pts, First[response],
  246.                         First[weights], model, vars, params]
  247.         ];
  248.     (* - check degrees of freedom *)
  249.         If[Length[params] > Length[pts],
  250.             Message[NonlinearFit::degrees]
  251.         ];
  252.     (* - check assorted options *)
  253.         If[!(IntegerQ[maxits] && Positive[maxits]),
  254.             Message[NonlinearFit::badits,maxits];
  255.                 maxits = 30];
  256.         If[!(IntegerQ[workprec] && Positive[workprec]),
  257.             Message[NonlinearFit::badwork, workprec];
  258.                 workprec = $MachinePrecision];
  259.         If[precgoal === Automatic, precgoal = workprec - 10];
  260.         If[accgoal === Automatic, accgoal = workprec - 10];
  261.         If[!(IntegerQ[precgoal] && Positive[precgoal]),
  262.             Message[NonlinearFit::badprec, precgoal];
  263.                 precgoal = $MachinePrecision - 10];
  264.         If[!(IntegerQ[accgoal] && Positive[accgoal]),
  265.             Message[NonlinearFit::badacc, accgoal];
  266.                 accgoal = $MachinePrecision - 10];
  267.         If[!MemberQ[{LevenbergMarquardt,FindMinimum},method],
  268.             Message[NonlinearFit::badmeth,method];
  269.                 method = LevenbergMarquardt];
  270.         progress = TrueQ[progress];
  271.     (* Call proper fitting routine *)
  272.     (* Implementation note: for now, we thread across all the
  273.             responses;  none of the algorithms in use can be made
  274.             more efficient when all are done at once.  This should
  275.             change if an algorithm that can take advantage of this
  276.             is used... *)
  277.         out = MapThread[
  278.             nlfit[method, pts, #1, #2, model, vars, params, start, minrng,
  279.                 maxrng, progress, maxits, workprec, precgoal, accgoal,
  280.                 fmopts]&,
  281.             {response,weights}];
  282.     (* check output *)
  283.         If[!MatchQ[out, {{__Rule}..}],
  284.             Message[NonlinearFit::fitfail,out];
  285.                 Return[$Failed]
  286.         ];
  287.         If[Length[out] === 1, First[out], out]
  288.     ]
  289.  
  290. (* SteepestDescent method;  creates symbolic form of chi-squared,
  291.     then uses built-in function FindMinimum to minimize this, given
  292.     the starting point.
  293. *)
  294.  
  295. nlfit[FindMinimum, data_, response_, weights_, model_, vars_,
  296.         parameters_, start_, minr_, maxr_, prog_, mi_, wp_, pg_, ag_,
  297.         opts___] :=
  298.     Module[{modelF, chisq, tmp, inc, ranges},
  299.         ranges = Map[If[FreeQ[#,DirectedInfinity[_]],#,Take[#,2]]&,
  300.                 Transpose[{parameters, start, minr, maxr}]
  301.         ];
  302.         modelF = Function[vars, model];
  303.            chisq = Plus @@ MapThread[
  304.                     (#3 (#2 - modelF @@ #1)^2) &,
  305.                 {data,response,weights}];
  306.         inc = 1;
  307.         If[prog,
  308.             Last[FindMinimum[
  309.                 (If[Head[First[parameters]] =!= Symbol,
  310.                     printprogress[inc++,tmp = chisq, parameters]
  311.                 ]; tmp),
  312.                 Evaluate[Sequence @@ ranges],
  313.                 MaxIterations -> mi,
  314.                 WorkingPrecision -> wp,
  315.                 AccuracyGoal -> ag + 1,
  316.                 PrecisionGoal -> pg + 1,
  317.                 opts,
  318.                 Evaluate[Gradient -> Map[D[chisq,#]&,parameters]]
  319.             ]],
  320.             Last[FindMinimum[chisq,
  321.                 Evaluate[Sequence @@ ranges],
  322.                 MaxIterations -> mi,
  323.                 WorkingPrecision -> wp,
  324.                 AccuracyGoal -> ag + 1,
  325.                 PrecisionGoal -> pg + 1,
  326.                 opts]]
  327.         ]
  328.     ]
  329.  
  330. nlfit[LevenbergMarquardt, data_, response_, weights_, model_, vars_,
  331.         params_, start_, minr_, maxr_, progress_, maxits_, 
  332.         wp_, pg_, ag_, ___] :=
  333.     Module[{ndata = N[data, wp], nresp = N[response, wp],
  334.             nweights = N[weights, wp], guess = N[start, wp], 
  335.             nminr = N[minr, wp], nmaxr = N[maxr, wp], lambda = 1/1000,
  336.             maplocs = Transpose[{#,#}]&[Range[Length[params]]],
  337.             alpha, beta, chisq, flag, lastflag, its, deltaguess,
  338.             tmpchi, newchi, newvals, vals, derivs, rngflag, oldguess,
  339.             accgoal = 10^-ag, precgoal = 10^-pg, ovflag = False},
  340.     (* initialize derivatives and function *)
  341.         vals = Apply[Function[vars,model],ndata,{1}];
  342.         derivs = Function[Evaluate[params,
  343.             Map[D[vals,#]&,params]]];
  344.         vals = Function[Evaluate[params, vals]];
  345.     (* initialize other variables *)
  346.         {alpha, beta} =
  347.             calcstep[response, nweights, vals @@ guess, derivs @@ guess,wp];
  348.         chisq = calcchisq[response, nweights, vals @@ guess,wp];
  349.         tmpchi = chisq; newchi = -1; its = 1; oldguess = guess;
  350.     (* main loop *)
  351.         ovflag = Check[While[!((flag = ((0 <= chisq - newchi < accgoal) || 
  352.                     (0 <= chisq - newchi < precgoal newchi))) &&
  353.                     lastflag) &&
  354.                 its <= maxits &&
  355.                 (rngflag = And @@ Thread[N[nminr < guess < nmaxr,wp]]),
  356.             lastflag = flag;
  357.             chisq = tmpchi;
  358.             oldguess = guess;
  359.             ++its;
  360.             deltaguess = LinearSolve[
  361.                     MapAt[# (1 + lambda)&,alpha,maplocs],
  362.                     beta];
  363.             newvals = vals @@ (guess + deltaguess);
  364.             newchi = calcchisq[response, nweights, newvals,wp];
  365.             If[progress,
  366.                 printprogress[its - 1, chisq, guess]
  367.             ];
  368.             If[newchi >= chisq,
  369.                 lambda = 10 lambda,
  370.                 lambda = lambda/10;
  371.                     guess = guess + deltaguess;
  372.                     tmpchi = newchi;
  373.                     {alpha, beta} = calcstep[response, nweights,
  374.                         newvals, derivs @@ guess,wp]
  375.             ]
  376.         ],True,General::ovfl];
  377.     (* check convergence of parameters *)
  378.     (* - Warning: this is not a truly appropriate test, but I
  379.             couldn't think of one that worked the way I wanted.  May
  380.             revise this after talking to Jerry Keiper... --John N. *)
  381.         If[Not[And @@ Thread[Abs[oldguess - guess] < 1]] && rngflag,
  382.             Message[NonlinearFit::lmpnocon]
  383.         ];
  384.     (* check other error possibilities *)
  385.         If[!rngflag,
  386.             Message[NonlinearFit::lmoutofrng]
  387.         ];
  388.         If[its > maxits,
  389.             Message[NonlinearFit::lmnoconv,maxits]
  390.         ];
  391.         If[TrueQ[ovflag],
  392.             Message[NonlinearFit::lmovfl]
  393.         ];
  394.     (* turn parameters into rules *)
  395.         Apply[Rule, Transpose[{params, guess}],{1}]
  396.     ]
  397.  
  398. calcchisq[resp_,weight_, vals_,wp_] :=
  399.     N[Plus @@ (weight (resp - vals)^2),wp]
  400.  
  401. calcstep[resp_, weight_, vals_, derivs_,wp_] :=
  402.     Module[{betatmp,hold},
  403.         betatmp = weight (resp - vals);
  404.         N[{Outer[Dot[weight List @@ #1, List @@ #2]&,
  405.                 #, #]&[hold @@ derivs]/.hold->List,
  406.             Map[betatmp . # &, derivs]},wp]
  407.     ]
  408.  
  409. (* utility to choose starting positions based on a minimum and a maximum
  410.     for the parameter *)
  411.  
  412. startposns[min_,max_] := {min + (1/3)(max - min), min + (2/3)(max - min)}
  413.  
  414. (* utility to choose a starting position from those derived with the
  415.     above utility by picking the point with the smallest chi squared. *)
  416.  
  417. findbeststart[starts_, data_, response_, weights_, model_, vars_, params_] :=
  418.     Module[{chis, chisq, pairs},
  419.         chisq = Function[params, 
  420.                     Evaluate[Plus @@ MapThread[
  421.                             (#3 (#2 - Function[vars, model] @@ #1)^2) &,
  422.                         {data,response,weights}]
  423.                     ]
  424.                 ];
  425.         pairs = Transpose[{starts,chis = N[Apply[chisq, starts,{1}]]}];
  426.         First[First[Select[pairs, #[[2]] == Min[chis] &]]]
  427.     ]
  428.  
  429.  
  430. (* utility to allow easy formatting of the print progress step... *)
  431.  
  432. printprogress[iteration_, chisq_, params_] :=
  433.     Print[StringForm[
  434.         "Iteration:`1`  ChiSquared:`2`  Parameters:`3`",
  435.             iteration, chisq, params]]
  436.  
  437.  
  438. End[]
  439.  
  440. EndPackage[]
  441.