1 # Created by Octave 3.6.2, Tue Jun 12 20:02:41 2012 UTC <root@t61>
13 # name: <cell-element>
17 general linear regression
19 [p,y_var,r,p_var]=LinearRegression(F,y)
20 [p,y_var,r,p_var]=LinearRegression(F,y,weight)
22 determine the parameters p_j (j=1,2,...,m) such that the function
23 f(x) = sum_(i=1,...,m) p_j*f_j(x) fits as good as possible to the
24 given values y_i = f(x_i)
27 F n*m matrix with the values of the basis functions at the support points
28 in column j give the values of f_j at the points x_i (i=1,2,...,n)
29 y n column vector of given values
30 weight n column vector of given weights
33 p m vector with the estimated values of the parameters
34 y_var estimated variance of the error
35 r weighted norm of residual
36 p_var estimated variance of the parameters p_j
40 # name: <cell-element>
44 general linear regression
49 # name: <cell-element>
56 # name: <cell-element>
60 ADSMAX Alternating directions method for direct search optimization.
61 [x, fmax, nf] = ADSMAX(FUN, x0, STOPIT, SAVIT, P) attempts to
62 maximize the function FUN, using the starting vector x0.
63 The alternating directions direct search method is used.
65 x = vector yielding largest function value found,
66 fmax = function value at x,
67 nf = number of function evaluations.
68 The iteration is terminated when either
69 - the relative increase in function value between successive
70 iterations is <= STOPIT(1) (default 1e-3),
71 - STOPIT(2) function evaluations have been performed
72 (default inf, i.e., no limit), or
73 - a function value equals or exceeds STOPIT(3)
74 (default inf, i.e., no test on function values).
75 Progress of the iteration is not shown if STOPIT(5) = 0 (default 1).
76 If a non-empty fourth parameter string SAVIT is present, then
77 `SAVE SAVIT x fmax nf' is executed after each inner iteration.
78 By default, the search directions are the co-ordinate directions.
79 The columns of a fifth parameter matrix P specify alternative search
80 directions (P = EYE is the default).
81 NB: x0 can be a matrix. In the output argument, in SAVIT saves,
82 and in function calls, x has the same shape as x0.
83 ADSMAX(fun, x0, STOPIT, SAVIT, P, P1, P2,...) allows additional
84 arguments to be passed to fun, via feval(fun,x,P1,P2,...).
86 N. J. Higham, Optimization by direct search in matrix computations,
87 SIAM J. Matrix Anal. Appl, 14(2): 317-333, 1993.
88 N. J. Higham, Accuracy and Stability of Numerical Algorithms,
89 Second edition, Society for Industrial and Applied Mathematics,
90 Philadelphia, PA, 2002; sec. 20.5.
94 # name: <cell-element>
98 ADSMAX Alternating directions method for direct search optimization.
102 # name: <cell-element>
109 # name: <cell-element>
113 battery.m: repeatedly call bfgs using a battery of
114 start values, to attempt to find global min
115 of a nonconvex function
118 func: function to mimimize
119 args: args of function
120 minarg: argument to minimize w.r.t. (usually = 1)
121 startvals: kxp matrix of values to try for sure (don't include all zeros, that's automatic)
122 max iters per start value
123 number of additional random start values to try
125 OUTPUT: theta - the best value found - NOT iterated to convergence
129 # name: <cell-element>
137 # name: <cell-element>
144 # name: <cell-element>
148 bfgsmin: bfgs or limited memory bfgs minimization of function
150 Usage: [x, obj_value, convergence, iters] = bfgsmin(f, args, control)
152 The function must be of the form
153 [value, return_2,..., return_m] = f(arg_1, arg_2,..., arg_n)
154 By default, minimization is w.r.t. arg_1, but it can be done
155 w.r.t. any argument that is a vector. Numeric derivatives are
156 used unless analytic derivatives are supplied. See bfgsmin_example.m
160 * f: name of function to minimize (string)
161 * args: a cell array that holds all arguments of the function
162 The argument with respect to which minimization is done
164 * control: an optional cell array of 1-8 elements. If a cell
165 array shorter than 8 elements is provided, the trailing elements
166 are provided with default values.
167 * elem 1: maximum iterations (positive integer, or -1 or Inf for unlimited (default))
169 0 = no screen output (default)
170 1 = only final results
171 2 = summary every iteration
172 3 = detailed information
173 * elem 3: convergence criterion
174 1 = strict (function, gradient and param change) (default)
175 0 = weak - only function convergence required
176 * elem 4: arg in f_args with respect to which minimization is done (default is first)
177 * elem 5: (optional) Memory limit for lbfgs. If it's a positive integer
178 then lbfgs will be use. Otherwise ordinary bfgs is used
179 * elem 6: function change tolerance, default 1e-12
180 * elem 7: parameter change tolerance, default 1e-6
181 * elem 8: gradient tolerance, default 1e-5
185 * obj_value: the value of f() at x
186 * convergence: 1 if normal conv, other values if not
187 * iters: number of iterations performed
189 Example: see bfgsmin_example.m
193 # name: <cell-element>
197 bfgsmin: bfgs or limited memory bfgs minimization of function
202 # name: <cell-element>
209 # name: <cell-element>
217 # name: <cell-element>
226 # name: <cell-element>
233 # name: <cell-element>
237 -- Function File: [S,V,N] brent_line_min ( F,DF,ARGS,CTL )
238 Line minimization of f along df
240 Finds minimum of f on line x0 + dx*w | a < w < b by bracketing.
241 a and b are passed through argument ctl.
246 * F : string : Name of function. Must return a real value
248 * ARGS : cell : Arguments passed to f or RxC : f's only
249 argument. x0 must be at ARGS{ CTL(2) }
251 * CTL : 5 : (optional) Control variables, described
257 * S : 1 : Minimum is at x0 + s*dx
259 * V : 1 : Value of f at x0 + s*dx
261 * NEV : 1 : Number of function evaluations
266 * CTL(1) : Upper bound for error on s
269 * CTL(2) : Position of minimized argument in args
272 * CTL(3) : Maximum number of function evaluations
281 Default values will be used if ctl is not passed or if nan values
287 # name: <cell-element>
291 Line minimization of f along df
296 # name: <cell-element>
303 # name: <cell-element>
307 -- Function File: cauchy (N, R, X, F )
308 Return the Taylor coefficients and numerical differentiation of a
309 function F for the first N-1 coefficients or derivatives using the
310 fft. N is the number of points to evaluate, R is the radius of
311 convergence, needs to be chosen less then the smallest singularity,
312 X is point to evaluate the Taylor expansion or differentiation.
315 If X is a scalar, the function F is evaluated in a row vector of
316 length N. If X is a column vector, F is evaluated in a matrix of
317 length(x)-by-N elements and must return a matrix of the same size.
319 d = cauchy(16, 1.5, 0, @(x) exp(x));
320 => d(2) = 1.0000 # first (2-1) derivative of function f (index starts from zero)
325 # name: <cell-element>
329 Return the Taylor coefficients and numerical differentiation of a
334 # name: <cell-element>
341 # name: <cell-element>
345 c = cdiff (func,wrt,N,dfunc,stack,dx) - Code for num. differentiation
346 = "function df = dfunc (var1,..,dvar,..,varN) .. endfunction
348 Returns a string of octave code that defines a function 'dfunc' that
349 returns the derivative of 'func' with respect to it's 'wrt'th
352 The derivatives are obtained by symmetric finite difference.
354 dfunc()'s return value is in the same format as that of ndiff()
356 func : string : name of the function to differentiate
358 wrt : int : position, in argument list, of the differentiation
361 N : int : total number of arguments taken by 'func'.
362 If N=inf, dfunc will take variable argument list.
365 dfunc : string : Name of the octave function that returns the
366 derivatives. Default:['d',func]
368 stack : string : Indicates whether 'func' accepts vertically
369 (stack="rstack") or horizontally (stack="cstack")
370 arguments. Any other string indicates that 'func'
371 does not allow stacking. Default:''
373 dx : real : Step used in the symmetric difference scheme.
376 See also : ndiff, eval, todisk
380 # name: <cell-element>
384 c = cdiff (func,wrt,N,dfunc,stack,dx) - Code for num.
388 # name: <cell-element>
395 # name: <cell-element>
399 -- Function File: [X0,V,NEV] cg_min ( F,DF,ARGS,CTL )
400 NonLinear Conjugate Gradient method to minimize function F.
405 * F : string : Name of function. Return a real value
407 * DF : string : Name of f's derivative. Returns a (R*C) x 1
410 * ARGS: cell : Arguments passed to f.
411 * CTL : 5-vec : (Optional) Control variables, described
417 * X0 : matrix : Local minimum of f
419 * V : real : Value of f in x0
421 * NEV : 1 x 2 : Number of evaluations of f and of df
426 * CTL(1) : 1 or 2 : Select stopping criterion amongst :
428 * CTL(1)==0 : Default value
430 * CTL(1)==1 : Stopping criterion : Stop search when value
431 doesn't improve, as tested by ctl(2) > Deltaf/max(|f(x)|,1)
432 where Deltaf is the decrease in f observed in the last
433 iteration (each iteration consists R*C line searches).
435 * CTL(1)==2 : Stopping criterion : Stop search when updates
436 are small, as tested by ctl(2) > max { dx(i)/max(|x(i)|,1) |
437 i in 1..N } where dx is the change in the x that occured in
440 * CTL(2) : Threshold used in stopping tests.
443 * CTL(2)==0 : Default value
445 * CTL(3) : Position of the minimized argument in args
448 * CTL(3)==0 : Default value
450 * CTL(4) : Maximum number of function evaluations
453 * CTL(4)==0 : Default value
455 * CTL(5) : Type of optimization:
457 * CTL(5)==1 : "Fletcher-Reves" method
459 * CTL(5)==2 : "Polak-Ribiere" (Default)
461 * CTL(5)==3 : "Hestenes-Stiefel" method
463 CTL may have length smaller than 4. Default values will be used if
464 ctl is not passed or if nan values are given.
469 function r=df( l ) b=[1;0;-1]; r = -( 2*l{1} - 2*b +
470 rand(size(l{1}))); endfunction
471 function r=ff( l ) b=[1;0;-1]; r = (l{1}-b)' * (l{1}-b);
475 [x0,v,nev]=cg_min( "ff", "df", ll, ctl )
476 Comment: In general, BFGS method seems to be better performin in
477 many cases but requires more computation per iteration
480 http://en.wikipedia.org/wiki/Nonlinear_conjugate_gradient
486 # name: <cell-element>
490 NonLinear Conjugate Gradient method to minimize function F.
494 # name: <cell-element>
501 # name: <cell-element>
505 [lb, idx, ridx, mv] = cpiv_bard (v, m[, incl])
507 v: column vector; m: matrix; incl (optional): index. length (v)
508 must equal rows (m). Finds column vectors w and l with w == v + m *
509 l, w >= 0, l >= 0, l.' * w == 0. Chooses idx, w, and l so that
510 l(~idx) == 0, l(idx) == -inv (m(idx, idx)) * v(idx), w(idx) roughly
511 == 0, and w(~idx) == v(~idx) + m(idx, ~idx).' * l(idx). idx indexes
512 at least everything indexed by incl, but l(incl) may be < 0. lb:
513 l(idx) (column vector); idx: logical index, defined above; ridx:
514 ~idx & w roughly == 0; mv: [m, v] after performing a Gauss-Jordan
515 'sweep' (with gjp.m) on each diagonal element indexed by idx.
516 Except the handling of incl (which enables handling of equality
517 constraints in the calling code), this is called solving the
518 'complementary pivot problem' (Cottle, R. W. and Dantzig, G. B.,
519 'Complementary pivot theory of mathematical programming', Linear
520 Algebra and Appl. 1, 102--125. References for the current
521 algorithm: Bard, Y.: Nonlinear Parameter Estimation, p. 147--149,
522 Academic Press, New York and London 1974; Bard, Y., 'An eclectic
523 approach to nonlinear programming', Proc. ANU Sem. Optimization,
524 Canberra, Austral. Nat. Univ.).
528 # name: <cell-element>
532 [lb, idx, ridx, mv] = cpiv_bard (v, m[, incl])
537 # name: <cell-element>
544 # name: <cell-element>
548 -- Function File: INFO = residmin_stat (F, P, X, Y, SETTINGS)
549 Frontend for computation of statistics for fitting of values,
550 computed by a model function, to observed values.
552 Please refer to the description of `residmin_stat'. The only
553 differences to `residmin_stat' are the additional arguments X
554 (independent values) and Y (observations), that the model function
555 F, if provided, has a second obligatory argument which will be set
556 to X and is supposed to return guesses for the observations (with
557 the same dimensions), and that the possibly user-supplied function
558 for the jacobian of the model function has also a second
559 obligatory argument which will be set to X.
561 See also: residmin_stat
567 # name: <cell-element>
571 Frontend for computation of statistics for fitting of values, computed
576 # name: <cell-element>
583 # name: <cell-element>
587 [x,v,nev,h,args] = d2_min(f,d2f,args,ctl,code) - Newton-like minimization
589 Minimize f(x) using 1st and 2nd derivatives. Any function w/ second
590 derivatives can be minimized, as in Newton. f(x) decreases at each
591 iteration, as in Levenberg-Marquardt. This function is inspired from the
592 Levenberg-Marquardt algorithm found in the book "Numerical Recipes".
595 f : string : Cost function's name
597 d2f : string : Name of function returning the cost (1x1), its
598 differential (1xN) and its second differential or it's
599 pseudo-inverse (NxN) (see ctl(5) below) :
601 [v,dv,d2v] = d2f (x).
603 args : list : f and d2f's arguments. By default, minimize the 1st
604 or matrix : argument.
606 ctl : vector : Control arguments (see below)
609 code : string : code will be evaluated after each outer loop that
610 produced some (any) improvement. Variables visible from
611 "code" include "x", the best parameter found, "v" the
612 best value and "args", the list of all arguments. All can
613 be modified. This option can be used to re-parameterize
614 the argument space during optimization
616 CONTROL VARIABLE ctl : (optional). May be a struct or a vector of length
617 ---------------------- 5 or less where NaNs are ignored. Default values
622 ftol, f N/A : Stop search when value doesn't improve, as tested by
624 f > Deltaf/max(|f(x)|,1)
626 where Deltaf is the decrease in f observed in the last
627 iteration. <10*sqrt(eps)>
629 utol, u N/A : Stop search when updates are small, as tested by
631 u > max { dx(i)/max(|x(i)|,1) | i in 1..N }
633 where dx is the change in the x that occured in the last
636 dtol, d N/A : Stop search when derivative is small, as tested by
640 crit, c ctl(1) : Set one stopping criterion, 'ftol' (c=1), 'utol' (c=2)
641 or 'dtol' (c=3) to the value of by the 'tol' option. <1>
643 tol, t ctl(2) : Threshold in termination test chosen by 'crit' <10*eps>
645 narg, n ctl(3) : Position of the minimized argument in args <1>
646 maxev,m ctl(4) : Maximum number of function evaluations <inf>
647 maxout,m : Maximum number of outer loops <inf>
648 id2f, i ctl(5) : 0 if d2f returns the 2nd derivatives, 1 if <0>
649 it returns its pseudo-inverse.
651 verbose, v N/A : Be more or less verbose (quiet=0) <0>
655 # name: <cell-element>
659 [x,v,nev,h,args] = d2_min(f,d2f,args,ctl,code) - Newton-like minimization
664 # name: <cell-element>
671 # name: <cell-element>
675 function prt = dcdp (f, p, dp, func[, bounds])
677 This is an interface to __dfdp__.m, similar to dfdp.m, but for
678 functions only of parameters 'p', not of independents 'x'. See
681 dfpdp is more general and is meant to be used instead of dcdp in
686 # name: <cell-element>
690 function prt = dcdp (f, p, dp, func[, bounds])
695 # name: <cell-element>
702 # name: <cell-element>
706 de_min: global optimisation using differential evolution
708 Usage: [x, obj_value, nfeval, convergence] = de_min(fcn, control)
710 minimization of a user-supplied function with respect to x(1:D),
711 using the differential evolution (DE) method based on an algorithm
712 by Rainer Storn (http://www.icsi.berkeley.edu/~storn/code.html)
713 See: http://www.softcomputing.net/tevc2009_1.pdf
718 fcn string : Name of function. Must return a real value
719 control vector : (Optional) Control variables, described below
724 x vector : parameter vector of best solution
725 obj_value scalar : objective function value of best solution
726 nfeval scalar : number of function evaluations
727 convergence : 1 = best below value to reach (VTR)
728 0 = population has reached defined quality (tol)
729 -1 = some values are close to constraints/boundaries
730 -2 = max number of iterations reached (maxiter)
731 -3 = max number of functions evaluations reached (maxnfe)
733 Control variable: (optional) may be named arguments (i.e. "name",value
734 ---------------- pairs), a struct, or a vector, where
737 XVmin : vector of lower bounds of initial population
738 *** note: by default these are no constraints ***
739 XVmax : vector of upper bounds of initial population
740 constr : 1 -> enforce the bounds not just for the initial population
741 const : data vector (remains fixed during the minimization)
742 NP : number of population members
743 F : difference factor from interval [0, 2]
744 CR : crossover probability constant from interval [0, 1]
745 strategy : 1 --> DE/best/1/exp 7 --> DE/best/1/bin
746 2 --> DE/rand/1/exp 8 --> DE/rand/1/bin
747 3 --> DE/target-to-best/1/exp 9 --> DE/target-to-best/1/bin
748 4 --> DE/best/2/exp 10--> DE/best/2/bin
749 5 --> DE/rand/2/exp 11--> DE/rand/2/bin
750 6 --> DEGL/SAW/exp else DEGL/SAW/bin
751 refresh : intermediate output will be produced after "refresh"
752 iterations. No intermediate output will be produced
754 VTR : Stopping criterion: "Value To Reach"
755 de_min will stop when obj_value <= VTR.
756 Use this if you know which value you expect.
757 tol : Stopping criterion: "tolerance"
758 stops if (best-worst)/max(1,worst) < tol
759 This stops basically if the whole population is "good".
760 maxnfe : maximum number of function evaluations
761 maxiter : maximum number of iterations (generations)
763 The algorithm seems to work well only if [XVmin,XVmax] covers the
764 region where the global minimum is expected.
765 DE is also somewhat sensitive to the choice of the
766 difference factor F. A good initial guess is to choose F from
767 interval [0.5, 1], e.g. 0.8.
768 CR, the crossover probability constant from interval [0, 1]
769 helps to maintain the diversity of the population and is
770 rather uncritical but affects strongly the convergence speed.
771 If the parameters are correlated, high values of CR work better.
772 The reverse is true for no correlation.
773 Experiments suggest that /bin likes to have a slightly
775 The number of population members NP is also not very critical. A
776 good initial guess is 10*D. Depending on the difficulty of the
777 problem NP can be lower than 10*D or must be higher than 10*D
778 to achieve convergence.
797 Example to find the minimum of the Rosenbrock saddle:
798 ----------------------------------------------------
800 function result = f(x);
801 result = 100 * (x(2) - x(1)^2)^2 + (1 - x(1))^2;
807 [x, obj_value, nfeval, convergence] = de_min (@f, ctl);
809 Keywords: global-optimisation optimisation minimisation
813 # name: <cell-element>
817 de_min: global optimisation using differential evolution
822 # name: <cell-element>
829 # name: <cell-element>
833 -- Function File: DX = deriv (F, X0)
834 -- Function File: DX = deriv (F, X0, H)
835 -- Function File: DX = deriv (F, X0, H, O)
836 -- Function File: DX = deriv (F, X0, H, O, N)
837 Calculate derivate of function F.
839 F must be a function handle or the name of a function that takes X0
840 and returns a variable of equal length and orientation. X0 must be
841 a numeric vector or scalar.
843 H defines the step taken for the derivative calculation. Defaults
846 O defines the order of the calculation. Supported values are 2
847 (h^2 order) or 4 (h^4 order). Defaults to 2.
849 N defines the derivative order. Defaults to the 1st derivative of
850 the function. Can be up to the 4th derivative.
852 Reference: Numerical Methods for Mathematics, Science, and
853 Engineering by John H. Mathews.
858 # name: <cell-element>
862 Calculate derivate of function F.
866 # name: <cell-element>
873 # name: <cell-element>
877 function prt = dfdp (x, f, p, dp, func[, bounds])
878 numerical partial derivatives (Jacobian) df/dp for use with leasqr
879 --------INPUT VARIABLES---------
880 x=vec or matrix of indep var(used as arg to func) x=[x0 x1 ....]
881 f=func(x,p) vector initialsed by user before each call to dfdp
882 p= vec of current parameter values
883 dp= fractional increment of p for numerical derivatives
884 dp(j)>0 central differences calculated
885 dp(j)<0 one sided differences calculated
886 dp(j)=0 sets corresponding partials to zero; i.e. holds p(j) fixed
887 func=function (string or handle) to calculate the Jacobian for,
888 e.g. to calc Jacobian for function expsum prt=dfdp(x,f,p,dp,'expsum')
889 bounds=two-column-matrix of lower and upper bounds for parameters
890 If no 'bounds' options is specified to leasqr, it will call
891 dfdp without the 'bounds' argument.
892 ----------OUTPUT VARIABLES-------
893 prt= Jacobian Matrix prt(i,j)=df(i)/dp(j)
894 ================================
896 dfxpdp is more general and is meant to be used instead of dfdp in
901 # name: <cell-element>
905 function prt = dfdp (x, f, p, dp, func[, bounds])
906 numerical partial derivative
910 # name: <cell-element>
917 # name: <cell-element>
921 function jac = dfpdp (p, func[, hook])
923 Returns Jacobian of func (p) with respect to p with finite
924 differencing. The optional argument hook is a structure which can
925 contain the following fields at the moment:
927 hook.f: value of func(p) for p as given in the arguments
929 hook.diffp: positive vector of fractional steps from given p in
930 finite differencing (actual steps may be smaller if bounds are
931 given). The default is .001 * ones (size (p)).
933 hook.diff_onesided: logical vector, indexing elements of p for
934 which only one-sided differences should be computed (faster); even
935 if not one-sided, differences might not be exactly central if
936 bounds are given. The default is false (size (p)).
938 hook.fixed: logical vector, indexing elements of p for which zero
939 should be returned instead of the guessed partial derivatives
940 (useful in optimization if some parameters are not optimized, but
943 hook.lbound, hook.ubound: vectors of lower and upper parameter
944 bounds (or -Inf or +Inf, respectively) to be respected in finite
945 differencing. The consistency of bounds is not checked.
949 # name: <cell-element>
953 function jac = dfpdp (p, func[, hook])
958 # name: <cell-element>
965 # name: <cell-element>
969 function jac = dfxpdp (x, p, func[, hook])
971 Returns Jacobian of func (p, x) with respect to p with finite
972 differencing. The optional argument hook is a structure which can
973 contain the following fields at the moment:
975 hook.f: value of func(p, x) for p and x as given in the arguments
977 hook.diffp: positive vector of fractional steps from given p in
978 finite differencing (actual steps may be smaller if bounds are
979 given). The default is .001 * ones (size (p));
981 hook.diff_onesided: logical vector, indexing elements of p for
982 which only one-sided differences should be computed (faster); even
983 if not one-sided, differences might not be exactly central if
984 bounds are given. The default is false (size (p)).
986 hook.fixed: logical vector, indexing elements of p for which zero
987 should be returned instead of the guessed partial derivatives
988 (useful in optimization if some parameters are not optimized, but
991 hook.lbound, hook.ubound: vectors of lower and upper parameter
992 bounds (or -Inf or +Inf, respectively) to be respected in finite
993 differencing. The consistency of bounds is not checked.
997 # name: <cell-element>
1001 function jac = dfxpdp (x, p, func[, hook])
1006 # name: <cell-element>
1013 # name: <cell-element>
1017 USAGE [alpha,c,rms] = expfit( deg, x1, h, y )
1019 Prony's method for non-linear exponential fitting
1021 Fit function: \sum_1^{deg} c(i)*exp(alpha(i)*x)
1023 Elements of data vector y must correspond to
1024 equidistant x-values starting at x1 with stepsize h
1026 The method is fully compatible with complex linear
1027 coefficients c, complex nonlinear coefficients alpha
1028 and complex input arguments y, x1, non-zero h .
1029 Fit-order deg must be a real positive integer.
1031 Returns linear coefficients c, nonlinear coefficients
1032 alpha and root mean square error rms. This method is
1033 known to be more stable than 'brute-force' non-linear
1034 least squares fitting.
1037 x0 = 0; step = 0.05; xend = 5; x = x0:step:xend;
1038 y = 2*exp(1.3*x)-0.5*exp(2*x);
1039 error = (rand(1,length(y))-0.5)*1e-4;
1040 [alpha,c,rms] = expfit(2,x0,step,y+error)
1050 The fit is very sensitive to the number of data points.
1051 It doesn't perform very well for small data sets.
1052 Theoretically, you need at least 2*deg data points, but
1053 if there are errors on the data, you certainly need more.
1055 Be aware that this is a very (very,very) ill-posed problem.
1056 By the way, this algorithm relies heavily on computing the
1057 roots of a polynomial. I used 'roots.m', if there is
1058 something better please use that code.
1060 Demo for a complex fit-function:
1061 deg= 2; N= 20; x1= -(1+i), x= linspace(x1,1+i/2,N).';
1063 y= (2+i)*exp( (-1-2i)*x ) + (-1+3i)*exp( (2+3i)*x );
1064 A= 5e-2; y+= A*(randn(N,1)+randn(N,1)*i); % add complex noise
1065 [alpha,c,rms]= expfit( deg, x1, h, y )
1069 # name: <cell-element>
1073 USAGE [alpha,c,rms] = expfit( deg, x1, h, y )
1078 # name: <cell-element>
1085 # name: <cell-element>
1093 # name: <cell-element>
1102 # name: <cell-element>
1109 # name: <cell-element>
1113 -- Function File: [X] = fmins(F,X0,OPTIONS,GRAD,P1,P2, ...)
1114 Find the minimum of a funtion of several variables. By default
1115 the method used is the Nelder&Mead Simplex algorithm
1117 Example usage: fmins(inline('(x(1)-5).^2+(x(2)-8).^4'),[0;0])
1121 A string containing the name of the function to minimize
1124 A vector of initial parameters fo the function F.
1127 Vector with control parameters (not all parameters are used) options(1) - Show progress (if 1, default is 0, no progress)
1128 options(2) - Relative size of simplex (default 1e-3)
1129 options(6) - Optimization algorithm
1130 if options(6)==0 - Nelder & Mead simplex (default)
1131 if options(6)==1 - Multidirectional search Method
1132 if options(6)==2 - Alternating Directions search
1134 if options(6)==0 && options(5)==0 - regular simplex
1135 if options(6)==0 && options(5)==1 - right-angled simplex
1136 Comment: the default is set to "right-angled simplex".
1137 this works better for me on a broad range of problems,
1138 although the default in nmsmax is "regular simplex"
1139 options(10) - Maximum number of function evaluations
1142 Unused (For compatibility with Matlab)
1145 Optional parameters for function F
1151 # name: <cell-element>
1155 Find the minimum of a funtion of several variables.
1159 # name: <cell-element>
1166 # name: <cell-element>
1170 -- Function File: X = fminsearch (FUN, X0)
1171 -- Function File: [X, FVAL] = fminsearch (FUN, X0, OPTIONS, GRAD, P1,
1173 Find the minimum of a funtion of several variables. By default
1174 the method used is the Nelder&Mead Simplex algorithm.
1176 See also: fmin, fmins, nmsmax
1182 # name: <cell-element>
1186 Find the minimum of a funtion of several variables.
1190 # name: <cell-element>
1197 # name: <cell-element>
1201 [x,v,flag,out,df,d2f] = fminunc_compat (f,x,opt,...) - M*tlab-like optimization
1203 Imitation of m*tlab's fminunc(). The optional 'opt' argument is a struct,
1204 e.g. produced by 'optimset()'. 'fminunc_compat' has been deprecated in
1205 favor of 'fminunc', which is now part of core Octave. This function
1206 will possibly be removed from future versions of the 'optim' package.
1210 Diagnostics, [off|on] : Be verbose
1211 Display , [off|iter|notify|final]
1212 : Be verbose unless value is "off"
1213 GradObj , [off|on] : Function's 2nd return value is derivatives
1214 Hessian , [off|on] : Function's 2nd and 3rd return value are
1215 derivatives and Hessian.
1216 TolFun , scalar : Termination criterion (see 'ftol' in minimize())
1217 TolX , scalar : Termination criterion (see 'utol' in minimize())
1218 MaxFunEvals, int : Max. number of function evaluations
1219 MaxIter , int : Max. number of algorithm iterations
1221 These non-m*tlab are provided to facilitate porting code to octave:
1222 -----------------------
1223 "MinEquiv" , [off|on] : Don't minimize 'fun', but instead return the
1224 option passed to minimize().
1226 "Backend" , [off|on] : Don't minimize 'fun', but instead return
1227 [backend, opt], the name of the backend
1228 optimization function that is used and the
1229 optional arguments that will be passed to it. See
1230 the 'backend' option of minimize().
1232 This function is a front-end to minimize().
1236 # name: <cell-element>
1240 [x,v,flag,out,df,d2f] = fminunc_compat (f,x,opt,.
1244 # name: <cell-element>
1251 # name: <cell-element>
1257 m: matrix; k, l: row- and column-index of pivot, l defaults to k.
1259 Gauss-Jordon pivot as defined in Bard, Y.: Nonlinear Parameter
1260 Estimation, p. 296, Academic Press, New York and London 1974. In
1261 the pivot column, this seems not quite the same as the usual
1262 Gauss-Jordan(-Clasen) pivot. Bard gives Beaton, A. E., 'The use of
1263 special matrix operators in statistical calculus' Research Bulletin
1264 RB-64-51 (1964), Educational Testing Service, Princeton, New Jersey
1265 as a reference, but this article is not easily accessible. Another
1266 reference, whose definition of gjp differs from Bards by some
1267 signs, is Clarke, R. B., 'Algorithm AS 178: The Gauss-Jordan sweep
1268 operator with detection of collinearity', Journal of the Royal
1269 Statistical Society, Series C (Applied Statistics) (1982), 31(2),
1274 # name: <cell-element>
1283 # name: <cell-element>
1290 # name: <cell-element>
1294 -- Function File: Df = jacobs (X, F)
1295 -- Function File: Df = jacobs (X, F, HOOK)
1296 Calculate the jacobian of a function using the complex step method.
1298 Let F be a user-supplied function. Given a point X at which we
1299 seek for the Jacobian, the function `jacobs' returns the Jacobian
1300 matrix `d(f(1), ..., df(end))/d(x(1), ..., x(n))'. The function
1301 uses the complex step method and thus can be applied to real
1304 The optional argument HOOK is a structure with additional options.
1305 HOOK can have the following fields:
1306 * `h' - can be used to define the magnitude of the complex step
1307 and defaults to 1e-20; steps larger than 1e-3 are not allowed.
1309 * `fixed' - is a logical vector internally usable by some
1310 optimization functions; it indicates for which elements of X
1311 no gradient should be computed, but zero should be returned.
1315 f = @(x) [x(1)^2 + x(2); x(2)*exp(x(1))];
1316 Df = jacobs ([1, 2], f)
1321 # name: <cell-element>
1325 Calculate the jacobian of a function using the complex step method.
1329 # name: <cell-element>
1336 # name: <cell-element>
1340 function [f,p,cvg,iter,corp,covp,covr,stdresid,Z,r2]=
1341 leasqr(x,y,pin,F,{stol,niter,wt,dp,dFdp,options})
1343 Levenberg-Marquardt nonlinear regression of f(x,p) to y(x).
1346 Optional parameters are in braces {}.
1347 x = vector or matrix of independent variables.
1348 y = vector or matrix of observed values.
1349 wt = statistical weights (same dimensions as y). These should be
1350 set to be proportional to (sqrt of var(y))^-1; (That is, the
1351 covariance matrix of the data is assumed to be proportional to
1352 diagonal with diagonal equal to (wt.^2)^-1. The constant of
1353 proportionality will be estimated.); default = ones( size (y)).
1354 pin = vec of initial parameters to be adjusted by leasqr.
1355 dp = fractional increment of p for numerical partial derivatives;
1356 default = .001*ones(size(pin))
1357 dp(j) > 0 means central differences on j-th parameter p(j).
1358 dp(j) < 0 means one-sided differences on j-th parameter p(j).
1359 dp(j) = 0 holds p(j) fixed i.e. leasqr wont change initial guess: pin(j)
1360 F = name of function in quotes or function handle; the function
1361 shall be of the form y=f(x,p), with y, x, p of the form y, x, pin
1363 dFdp = name of partial derivative function in quotes or function
1364 handle; default is 'dfdp', a slow but general partial derivatives
1365 function; the function shall be of the form
1366 prt=dfdp(x,f,p,dp,F[,bounds]). For backwards compatibility, the
1367 function will only be called with an extra 'bounds' argument if the
1368 'bounds' option is explicitely specified to leasqr (see dfdp.m).
1369 stol = scalar tolerance on fractional improvement in scalar sum of
1370 squares = sum((wt.*(y-f))^2); default stol = .0001;
1371 niter = scalar maximum number of iterations; default = 20;
1372 options = structure, currently recognized fields are 'fract_prec',
1373 'max_fract_change', 'inequc', 'bounds', and 'equc'. For backwards
1374 compatibility, 'options' can also be a matrix whose first and
1375 second column contains the values of 'fract_prec' and
1376 'max_fract_change', respectively.
1377 Field 'options.fract_prec': column vector (same length as 'pin')
1378 of desired fractional precisions in parameter estimates.
1379 Iterations are terminated if change in parameter vector (chg)
1380 relative to current parameter estimate is less than their
1381 corresponding elements in 'options.fract_prec' [ie. all (abs
1382 (chg) < abs (options.fract_prec .* current_parm_est))] on two
1383 consecutive iterations, default = zeros().
1384 Field 'options.max_fract_change': column vector (same length as
1385 'pin) of maximum fractional step changes in parameter vector.
1386 Fractional change in elements of parameter vector is constrained to
1387 be at most 'options.max_fract_change' between sucessive iterations.
1388 [ie. abs(chg(i))=abs(min([chg(i)
1389 options.max_fract_change(i)*current param estimate])).], default =
1391 Field 'options.inequc': cell-array containing up to four entries,
1392 two entries for linear inequality constraints and/or one or two
1393 entries for general inequality constraints. Initial parameters
1394 must satisfy these constraints. Either linear or general
1395 constraints may be the first entries, but the two entries for
1396 linear constraints must be adjacent and, if two entries are given
1397 for general constraints, they also must be adjacent. The two
1398 entries for linear constraints are a matrix (say m) and a vector
1399 (say v), specifying linear inequality constraints of the form
1400 `m.' * parameters + v >= 0'. If the constraints are just bounds,
1401 it is suggested to specify them in 'options.bounds' instead,
1402 since then some sanity tests are performed, and since the
1403 function 'dfdp.m' is guarantied not to violate constraints during
1404 determination of the numeric gradient only for those constraints
1405 specified as 'bounds' (possibly with violations due to a certain
1406 inaccuracy, however, except if no constraints except bounds are
1407 specified). The first entry for general constraints must be a
1408 differentiable vector valued function (say h), specifying general
1409 inequality constraints of the form `h (p[, idx]) >= 0'; p is the
1410 column vector of optimized paraters and the optional argument idx
1411 is a logical index. h has to return the values of all constraints
1412 if idx is not given, and has to return only the indexed
1413 constraints if idx is given (so computation of the other
1414 constraints can be spared). If a second entry for general
1415 constraints is given, it must be a function (say dh) which
1416 returnes a matrix whos rows contain the gradients of the
1417 constraint function h with respect to the optimized parameters.
1418 It has the form jac_h = dh (vh, p, dp, h, idx[, bounds]); p is
1419 the column vector of optimized parameters, and idx is a logical
1420 index --- only the rows indexed by idx must be returned (so
1421 computation of the others can be spared). The other arguments of
1422 dh are for the case that dh computes numerical gradients: vh is
1423 the column vector of the current values of the constraint
1424 function h, with idx already applied. h is a function h (p) to
1425 compute the values of the constraints for parameters p, it will
1426 return only the values indexed by idx. dp is a suggestion for
1427 relative step width, having the same value as the argument 'dp'
1428 of leasqr above. If bounds were specified to leasqr, they are
1429 provided in the argument bounds of dh, to enable their
1430 consideration in determination of numerical gradients. If dh is
1431 not specified to leasqr, numerical gradients are computed in the
1432 same way as with 'dfdp.m' (see above). If some constraints are
1433 linear, they should be specified as linear constraints (or
1434 bounds, if applicable) for reasons of performance, even if
1435 general constraints are also specified.
1436 Field 'options.bounds': two-column-matrix, one row for each
1437 parameter in 'pin'. Each row contains a minimal and maximal value
1438 for each parameter. Default: [-Inf, Inf] in each row. If this
1439 field is used with an existing user-side function for 'dFdp'
1440 (see above) the functions interface might have to be changed.
1441 Field 'options.equc': equality constraints, specified the same
1442 way as inequality constraints (see field 'options.inequc').
1443 Initial parameters must satisfy these constraints.
1444 Note that there is possibly a certain inaccuracy in honoring
1445 constraints, except if only bounds are specified.
1446 _Warning_: If constraints (or bounds) are set, returned guesses
1447 of corp, covp, and Z are generally invalid, even if no constraints
1448 are active for the final parameters. If equality constraints are
1449 specified, corp, covp, and Z are not guessed at all.
1450 Field 'options.cpiv': Function for complementary pivot algorithm
1451 for inequality constraints, default: cpiv_bard. No different
1452 function is supplied.
1455 f = column vector of values computed: f = F(x,p).
1456 p = column vector trial or final parameters. i.e, the solution.
1457 cvg = scalar: = 1 if convergence, = 0 otherwise.
1458 iter = scalar number of iterations used.
1459 corp = correlation matrix for parameters.
1460 covp = covariance matrix of the parameters.
1461 covr = diag(covariance matrix of the residuals).
1462 stdresid = standardized residuals.
1463 Z = matrix that defines confidence region (see comments in the source).
1464 r2 = coefficient of multiple determination, intercept form.
1466 Not suitable for non-real residuals.
1469 Bard, Nonlinear Parameter Estimation, Academic Press, 1974.
1470 Draper and Smith, Applied Regression Analysis, John Wiley and Sons, 1981.
1474 # name: <cell-element>
1478 function [f,p,cvg,iter,corp,covp,covr,stdresid,Z,r2]=
1483 # name: <cell-element>
1490 # name: <cell-element>
1494 [a,fx,nev] = line_min (f, dx, args, narg, h, nev_max) - Minimize f() along dx
1497 f : string : Name of minimized function
1498 dx : matrix : Direction along which f() is minimized
1499 args : cell : Arguments of f
1500 narg : integer : Position of minimized variable in args. Default=1
1501 h : scalar : Step size to use for centered finite difference
1502 approximation of first and second derivatives. Default=1E-3.
1503 nev_max : integer : Maximum number of function evaluations. Default=30
1506 a : scalar : Value for which f(x+a*dx) is a minimum (*)
1507 fx : scalar : Value of f(x+a*dx) at minimum (*)
1508 nev : integer : Number of function evaluations
1510 (*) The notation f(x+a*dx) assumes that args == {x}.
1512 Reference: David G Luenberger's Linear and Nonlinear Programming
1516 # name: <cell-element>
1520 [a,fx,nev] = line_min (f, dx, args, narg, h, nev_max) - Minimize f() along dx
1525 # name: <cell-element>
1532 # name: <cell-element>
1536 -- Function File: X = linprog (F, A, B)
1537 -- Function File: X = linprog (F, A, B, AEQ, BEQ)
1538 -- Function File: X = linprog (F, A, B, AEQ, BEQ, LB, UB)
1539 -- Function File: [X, FVAL] = linprog (...)
1540 Solve a linear problem.
1546 (both f and x are column vectors) subject to
1552 If not specified, AEQ and BEQ default to empty matrices.
1554 If not specified, the lower bound LB defaults to minus infinite
1555 and the upper bound UB defaults to infinite.
1563 # name: <cell-element>
1567 Solve a linear problem.
1571 # name: <cell-element>
1578 # name: <cell-element>
1582 MDSMAX Multidirectional search method for direct search optimization.
1583 [x, fmax, nf] = MDSMAX(FUN, x0, STOPIT, SAVIT) attempts to
1584 maximize the function FUN, using the starting vector x0.
1585 The method of multidirectional search is used.
1587 x = vector yielding largest function value found,
1588 fmax = function value at x,
1589 nf = number of function evaluations.
1590 The iteration is terminated when either
1591 - the relative size of the simplex is <= STOPIT(1)
1593 - STOPIT(2) function evaluations have been performed
1594 (default inf, i.e., no limit), or
1595 - a function value equals or exceeds STOPIT(3)
1596 (default inf, i.e., no test on function values).
1597 The form of the initial simplex is determined by STOPIT(4):
1598 STOPIT(4) = 0: regular simplex (sides of equal length, the default),
1599 STOPIT(4) = 1: right-angled simplex.
1600 Progress of the iteration is not shown if STOPIT(5) = 0 (default 1).
1601 If a non-empty fourth parameter string SAVIT is present, then
1602 `SAVE SAVIT x fmax nf' is executed after each inner iteration.
1603 NB: x0 can be a matrix. In the output argument, in SAVIT saves,
1604 and in function calls, x has the same shape as x0.
1605 MDSMAX(fun, x0, STOPIT, SAVIT, P1, P2,...) allows additional
1606 arguments to be passed to fun, via feval(fun,x,P1,P2,...).
1608 This implementation uses 2n^2 elements of storage (two simplices), where x0
1609 is an n-vector. It is based on the algorithm statement in [2, sec.3],
1610 modified so as to halve the storage (with a slight loss in readability).
1613 [1] V. J. Torczon, Multi-directional search: A direct search algorithm for
1614 parallel machines, Ph.D. Thesis, Rice University, Houston, Texas, 1989.
1615 [2] V. J. Torczon, On the convergence of the multidirectional search
1616 algorithm, SIAM J. Optimization, 1 (1991), pp. 123-145.
1617 [3] N. J. Higham, Optimization by direct search in matrix computations,
1618 SIAM J. Matrix Anal. Appl, 14(2): 317-333, 1993.
1619 [4] N. J. Higham, Accuracy and Stability of Numerical Algorithms,
1620 Second edition, Society for Industrial and Applied Mathematics,
1621 Philadelphia, PA, 2002; sec. 20.5.
1625 # name: <cell-element>
1629 MDSMAX Multidirectional search method for direct search optimization.
1633 # name: <cell-element>
1640 # name: <cell-element>
1644 [x,v,nev,...] = minimize (f,args,...) - Minimize f
1647 f : string : Name of function. Must return a real value
1648 args : list or : List of arguments to f (by default, minimize the first)
1649 matrix : f's only argument
1652 x : matrix : Local minimum of f. Let's suppose x is M-by-N.
1653 v : real : Value of f in x0
1654 nev : integer : Number of function evaluations
1655 or 1 x 2 : Number of function and derivative evaluations (if
1656 derivatives are used)
1659 Extra arguments are either a succession of option-value pairs or a single
1660 list or struct of option-value pairs (for unary options, the value in the
1663 OPTIONS : DERIVATIVES Derivatives may be used if one of these options
1664 --------------------- uesd. Otherwise, the Nelder-Mean (see
1665 nelder_mead_min) method is used.
1667 'd2f', d2f : Name of a function that returns the value of f, of its
1668 1st and 2nd derivatives : [fx,dfx,d2fx] = feval (d2f, x)
1669 where fx is a real number, dfx is 1x(M*N) and d2fx is
1670 (M*N)x(M*N). A Newton-like method (d2_min) will be used.
1672 'hess' : Use [fx,dfx,d2fx] = leval (f, args) to compute 1st and
1673 2nd derivatives, and use a Newton-like method (d2_min).
1675 'd2i', d2i : Name of a function that returns the value of f, of its
1676 1st and pseudo-inverse of second derivatives :
1677 [fx,dfx,id2fx] = feval (d2i, x) where fx is a real
1678 number, dfx is 1x(M*N) and d2ix is (M*N)x(M*N).
1679 A Newton-like method will be used (see d2_min).
1681 'ihess' : Use [fx,dfx,id2fx] = leval (f, args) to compute 1st
1682 derivative and the pseudo-inverse of 2nd derivatives,
1683 and use a Newton-like method (d2_min).
1685 NOTE : df, d2f or d2i take the same arguments as f.
1687 'order', n : Use derivatives of order n. If the n'th order derivative
1688 is not specified by 'df', 'd2f' or 'd2i', it will be
1689 computed numerically. Currently, only order 1 works.
1691 'ndiff' : Use a variable metric method (bfgs) using numerical
1694 OPTIONS : STOPPING CRITERIA Default is to use 'tol'
1695 ---------------------------
1696 'ftol', ftol : Stop search when value doesn't improve, as tested by
1698 ftol > Deltaf/max(|f(x)|,1)
1700 where Deltaf is the decrease in f observed in the last
1701 iteration. Default=10*eps
1703 'utol', utol : Stop search when updates are small, as tested by
1705 tol > max { dx(i)/max(|x(i)|,1) | i in 1..N }
1707 where dx is the change in the x that occured in the last
1710 'dtol',dtol : Stop search when derivatives are small, as tested by
1712 dtol > max { df(i)*max(|x(i)|,1)/max(v,1) | i in 1..N }
1714 where x is the current minimum, v is func(x) and df is
1715 the derivative of f in x. This option is ignored if
1716 derivatives are not used in optimization.
1720 'maxev', m : Maximum number of function evaluations <inf>
1722 'narg' , narg : Position of the minimized argument in args <1>
1723 'isz' , step : Initial step size (only for 0 and 1st order method) <1>
1724 Should correspond to expected distance to minimum
1725 'verbose' : Display messages during execution
1727 'backend' : Instead of performing the minimization itself, return
1728 [backend, control], the name and control argument of the
1729 backend used by minimize(). Minimimzation can then be
1730 obtained without the overhead of minimize by calling, if
1731 a 0 or 1st order method is used :
1733 [x,v,nev] = feval (backend, args, control)
1735 or, if a 2nd order method is used :
1737 [x,v,nev] = feval (backend, control.d2f, args, control)
1741 # name: <cell-element>
1749 # name: <cell-element>
1756 # name: <cell-element>
1760 [x0,v,nev] = nelder_mead_min (f,args,ctl) - Nelder-Mead minimization
1762 Minimize 'f' using the Nelder-Mead algorithm. This function is inspired
1763 from the that found in the book "Numerical Recipes".
1767 f : string : Name of function. Must return a real value
1768 args : list : Arguments passed to f.
1769 or matrix : f's only argument
1770 ctl : vector : (Optional) Control variables, described below
1775 x0 : matrix : Local minimum of f
1776 v : real : Value of f in x0
1777 nev : number : Number of function evaluations
1779 CONTROL VARIABLE : (optional) may be named arguments (i.e. "name",value
1780 ------------------ pairs), a struct, or a vector of length <= 6, where
1781 NaN's are ignored. Default values are written <value>.
1784 ftol,f N/A : Stopping criterion : stop search when values at simplex
1785 vertices are all alike, as tested by
1787 f > (max_i (f_i) - min_i (f_i)) /max(max(|f_i|),1)
1789 where f_i are the values of f at the vertices. <10*eps>
1791 rtol,r N/A : Stop search when biggest radius of simplex, using
1792 infinity-norm, is small, as tested by :
1794 ctl(2) > Radius <10*eps>
1796 vtol,v N/A : Stop search when volume of simplex is small, tested by
1800 crit,c ctl(1) : Set one stopping criterion, 'ftol' (c=1), 'rtol' (c=2)
1801 or 'vtol' (c=3) to the value of the 'tol' option. <1>
1803 tol, t ctl(2) : Threshold in termination test chosen by 'crit' <10*eps>
1805 narg ctl(3) : Position of the minimized argument in args <1>
1806 maxev ctl(4) : Maximum number of function evaluations. This number <inf>
1807 may be slightly exceeded.
1808 isz ctl(5) : Size of initial simplex, which is : <1>
1810 { x + e_i | i in 0..N }
1812 Where x == args{narg} is the initial value
1813 e_0 == zeros (size (x)),
1814 e_i(j) == 0 if j != i and e_i(i) == ctl(5)
1815 e_i has same size as x
1817 Set ctl(5) to the distance you expect between the starting
1818 point and the minimum.
1820 rst ctl(6) : When a minimum is found the algorithm restarts next to
1821 it until the minimum does not improve anymore. ctl(6) is
1822 the maximum number of restarts. Set ctl(6) to zero if
1823 you know the function is well-behaved or if you don't
1824 mind not getting a true minimum. <0>
1826 verbose, v Be more or less verbose (quiet=0) <0>
1830 # name: <cell-element>
1834 [x0,v,nev] = nelder_mead_min (f,args,ctl) - Nelder-Mead minimization
1839 # name: <cell-element>
1846 # name: <cell-element>
1850 NMSMAX Nelder-Mead simplex method for direct search optimization.
1851 [x, fmax, nf] = NMSMAX(FUN, x0, STOPIT, SAVIT) attempts to
1852 maximize the function FUN, using the starting vector x0.
1853 The Nelder-Mead direct search method is used.
1855 x = vector yielding largest function value found,
1856 fmax = function value at x,
1857 nf = number of function evaluations.
1858 The iteration is terminated when either
1859 - the relative size of the simplex is <= STOPIT(1)
1861 - STOPIT(2) function evaluations have been performed
1862 (default inf, i.e., no limit), or
1863 - a function value equals or exceeds STOPIT(3)
1864 (default inf, i.e., no test on function values).
1865 The form of the initial simplex is determined by STOPIT(4):
1866 STOPIT(4) = 0: regular simplex (sides of equal length, the default)
1867 STOPIT(4) = 1: right-angled simplex.
1868 Progress of the iteration is not shown if STOPIT(5) = 0 (default 1).
1869 STOPIT(6) indicates the direction (ie. minimization or
1870 maximization.) Default is 1, maximization.
1871 set STOPIT(6)=-1 for minimization
1872 If a non-empty fourth parameter string SAVIT is present, then
1873 `SAVE SAVIT x fmax nf' is executed after each inner iteration.
1874 NB: x0 can be a matrix. In the output argument, in SAVIT saves,
1875 and in function calls, x has the same shape as x0.
1876 NMSMAX(fun, x0, STOPIT, SAVIT, P1, P2,...) allows additional
1877 arguments to be passed to fun, via feval(fun,x,P1,P2,...).
1879 N. J. Higham, Optimization by direct search in matrix computations,
1880 SIAM J. Matrix Anal. Appl, 14(2): 317-333, 1993.
1881 C. T. Kelley, Iterative Methods for Optimization, Society for Industrial
1882 and Applied Mathematics, Philadelphia, PA, 1999.
1886 # name: <cell-element>
1890 NMSMAX Nelder-Mead simplex method for direct search optimization.
1894 # name: <cell-element>
1901 # name: <cell-element>
1905 -- Function File: [P, FY, CVG, OUTP] = nonlin_curvefit (F, PIN, X, Y)
1906 -- Function File: [P, FY, CVG, OUTP] = nonlin_curvefit (F, PIN, X, Y,
1908 Frontend for nonlinear fitting of values, computed by a model
1909 function, to observed values.
1911 Please refer to the description of `nonlin_residmin'. The only
1912 differences to `nonlin_residmin' are the additional arguments X
1913 (independent values, mostly, but not necessarily, an array of the
1914 same dimensions or the same number of rows as Y) and Y (array of
1915 observations), the returned value FY (final guess for observed
1916 values) instead of RESID, that the model function has a second
1917 obligatory argument which will be set to X and is supposed to
1918 return guesses for the observations (with the same dimensions),
1919 and that the possibly user-supplied function for the jacobian of
1920 the model function has also a second obligatory argument which
1923 See also: nonlin_residmin
1929 # name: <cell-element>
1933 Frontend for nonlinear fitting of values, computed by a model function,
1938 # name: <cell-element>
1945 # name: <cell-element>
1949 -- Function File: [P, OBJF, CVG, OUTP] = nonlin_min (F, PIN)
1950 -- Function File: [P, OBJF, CVG, OUTP] = nonlin_min (F, PIN, SETTINGS)
1951 Frontend for constrained nonlinear minimization of a scalar
1952 objective function. The functions supplied by the user have a
1953 minimal interface; any additionally needed constants can be
1954 supplied by wrapping the user functions into anonymous functions.
1956 The following description applies to usage with vector-based
1957 parameter handling. Differences in usage for structure-based
1958 parameter handling will be explained in a separate section below.
1960 F: objective function. It gets a column vector of real parameters
1961 as argument. In gradient determination, this function may be
1962 called with an informational second argument, whose content
1963 depends on the function for gradient determination.
1965 PIN: real column vector of initial parameters.
1967 SETTINGS: structure whose fields stand for optional settings
1968 referred to below. The fields can be set by `optimset()' with
1969 Octave versions 3.3.55 or greater; with older Octave versions, the
1970 fields must be set directly as structure-fields in the correct
1973 The returned values are the column vector of final parameters P,
1974 the final value of the objective function OBJF, an integer CVG
1975 indicating if and how optimization succeeded or failed, and a
1976 structure OUTP with additional information, curently with only one
1977 field: NITER, the number of iterations. CVG is greater than zero
1978 for success and less than or equal to zero for failure; its
1979 possible values depend on the used backend and currently can be
1980 `0' (maximum number of iterations exceeded), `1' (fixed number of
1981 iterations completed, e.g. in stochastic optimizers), `2'
1982 (parameter change less than specified precision in two consecutive
1983 iterations), `3' (improvement in objective function less than
1984 specified), or `-4' (algorithm got stuck).
1988 `Algorithm': String specifying the backend. Currently available
1989 are `"lm_feasible"' (default) and `"siman"'. They are described in
1990 separate sections below.
1992 `objf_grad': Function computing the gradient of the objective
1993 function with respect to the parameters, assuming residuals are
1994 reshaped to a vector. Default: finite differences. Will be called
1995 with the column vector of parameters and an informational structure
1996 as arguments. The structure has the fields `f': value of objective
1997 function for current parameters, `fixed': logical vector
1998 indicating which parameters are not optimized, so these partial
1999 derivatives need not be computed and can be set to zero, `diffp',
2000 `diff_onesided', `lbound', `ubound': identical to the user
2001 settings of this name, `plabels': 1-dimensional cell-array of
2002 column-cell-arrays, each column with labels for all parameters,
2003 the first column contains the numerical indices of the parameters.
2004 The default gradient function will call the objective function
2005 with the second argument set with fields `f': as the `f' passed to
2006 the gradient function, `plabels': cell-array of 1x1 cell-arrays
2007 with the entries of the column-cell-arrays of `plabels' as passed
2008 to the jacobian function corresponding to current parameter,
2009 `side': `0' for one-sided interval, `1' or `2', respectively, for
2010 the sides of a two-sided interval, and `parallel': logical scalar
2011 indicating parallel computation of partial derivatives.
2013 `objf_hessian': Function computing the Hessian of the objective
2014 function with respect to the parameters. The default is backend
2015 specific. Will be called with the column vector of parameters as
2018 `diffp': column vector of fractional intervals (doubled for
2019 central intervals) supposed to be used by gradient functions
2020 performing finite differencing. Default: `.001 * ones (size
2021 (parameters))'. The default gradient function will use these as
2022 absolute intervals for parameters with value zero.
2024 `diff_onesided': logical column vector indicating that one-sided
2025 intervals should be used by gradient functions performing finite
2026 differencing. Default: `false (size (parameters))'.
2028 `complex_step_derivative_objf', `complex_step_derivative_inequc',
2029 `complex_step_derivative_equc': logical scalars, default: false.
2030 Estimate gradient of objective function, general inequality
2031 constraints, and general equality constraints, respectively, with
2032 complex step derivative approximation. Use only if you know that
2033 your objective function, function of general inequality
2034 constraints, or function of general equality constraints,
2035 respectively, is suitable for this. No user function for the
2036 respective gradient must be specified.
2038 `cstep': scalar step size for complex step derivative
2039 approximation. Default: 1e-20.
2041 `fixed': logical column vector indicating which parameters should
2042 not be optimized, but kept to their inital value. Fixing is done
2043 independently of the backend, but the backend may choose to fix
2044 additional parameters under certain conditions.
2046 `lbound', `ubound': column vectors of lower and upper bounds for
2047 parameters. Default: `-Inf' and `+Inf', respectively. The bounds
2048 are non-strict, i.e. parameters are allowed to be exactly equal to
2049 a bound. The default gradient function will respect bounds (but no
2050 further inequality constraints) in finite differencing.
2052 `inequc': Further inequality constraints. Cell-array containing up
2053 to four entries, two entries for linear inequality constraints
2054 and/or one or two entries for general inequality constraints.
2055 Either linear or general constraints may be the first entries, but
2056 the two entries for linear constraints must be adjacent and, if
2057 two entries are given for general constraints, they also must be
2058 adjacent. The two entries for linear constraints are a matrix (say
2059 `m') and a vector (say `v'), specifying linear inequality
2060 constraints of the form `m.' * parameters + v >= 0'. The first
2061 entry for general constraints must be a differentiable vector
2062 valued function (say `h'), specifying general inequality
2063 constraints of the form `h (p[, idx]) >= 0'; `p' is the column
2064 vector of optimized paraters and the optional argument `idx' is a
2065 logical index. `h' has to return the values of all constraints if
2066 `idx' is not given. It may choose to return only the indexed
2067 constraints if `idx' is given (so computation of the other
2068 constraints can be spared); in this case, the additional setting
2069 `inequc_f_idx' has to be set to `true'. In gradient determination,
2070 this function may be called with an informational third argument,
2071 whose content depends on the function for gradient determination.
2072 If a second entry for general inequality constraints is given, it
2073 must be a function computing the jacobian of the constraints with
2074 respect to the parameters. For this function, the description of
2075 `dfdp' above applies, with 2 exceptions: 1) it is called with 3
2076 arguments since it has an additional argument `idx' -- a logical
2077 index -- at second position, indicating which rows of the jacobian
2078 must be returned (if the function chooses to return only indexed
2079 rows, the additional setting `inequc_df_idx' has to be set to
2080 `true'). 2) the default jacobian function calls `h' with 3
2081 arguments, since the argument `idx' is also supplied. Note that
2082 specifying linear constraints as general constraints will generally
2083 waste performance, even if further, non-linear, general constraints
2086 `equc': Equality constraints. Specified the same way as inequality
2087 constraints (see `inequc'). The respective additional settings are
2088 named `equc_f_idx' and `equc_df_idx'.
2090 `cpiv': Function for complementary pivoting, usable in algorithms
2091 for constraints. Default: cpiv_bard. Only the default function is
2092 supplied with the package.
2094 `TolFun': Minimum fractional improvement in objective function in
2095 an iteration (abortion criterium). Default: .0001.
2097 `MaxIter': Maximum number of iterations (abortion criterium).
2098 Default: backend-specific.
2100 `fract_prec': Column Vector, minimum fractional change of
2101 parameters in an iteration (abortion criterium if violated in two
2102 consecutive iterations). Default: backend-specific.
2104 `max_fract_change': Column Vector, enforced maximum fractional
2105 change in parameters in an iteration. Default: backend-specific.
2107 `Display': String indicating the degree of verbosity. Default:
2108 `"off"'. Possible values are currently `"off"' (no messages) and
2109 `"iter"' (some messages after each iteration). Support of this
2110 setting and its exact interpretation are backend-specific.
2112 `debug': Logical scalar, default: `false'. Will be passed to the
2113 backend, which might print debugging information if true.
2115 Structure-based parameter handling
2117 The setting `param_order' is a cell-array with names of the
2118 optimized parameters. If not given, and initial parameters are a
2119 structure, all parameters in the structure are optimized. If
2120 initial parameters are a structure, it is an error if
2121 `param_order' is not given and there are any non-structure-based
2122 configuration items or functions.
2124 The initial parameters PIN can be given as a structure containing
2125 at least all fields named in `param_order'. In this case the
2126 returned parameters P will also be a structure.
2128 Each user-supplied function can be called with the argument
2129 containing the current parameters being a structure instead of a
2130 column vector. For this, a corresponding setting must be set to
2131 `true': `objf_pstruct' (objective function), `objf_grad_pstruct'
2132 (gradient of objective function), `objf_hessian_pstruct' (hessian
2133 of objective function), `f_inequc_pstruct' (general inequality
2134 constraints), `df_inequc_pstruct' (jacobian of general inequality
2135 constraints), `f_equc_pstruct' (general equality constraints), and
2136 `df_equc_pstruct' (jacobian of general equality constraints). If a
2137 gradient (jacobian) function is configured in such a way, it must
2138 return the entries (columns) of the gradient (jacobian) as fields
2139 of a structure under the respective parameter names. If the
2140 hessian function is configured in such a way, it must return a
2141 structure (say `h') with fields e.g. as `h.a.b = value' for
2142 `value' being the 2nd partial derivative with respect to `a' and
2143 `b'. There is no need to also specify the field `h.b.a' in this
2146 Similarly, for specifying linear constraints, instead of the matrix
2147 (called `m' above), a structure containing the rows of the matrix
2148 in fields under the respective parameter names can be given. In
2149 this case, rows containing only zeros need not be given.
2151 The vector-based settings `lbound', `ubound', `fixed', `diffp',
2152 `diff_onesided', `fract_prec', and `max_fract_change' can be
2153 replaced by the setting `param_config'. It is a structure that can
2154 contain fields named in `param_order'. For each such field, there
2155 may be subfields with the same names as the above vector-based
2156 settings, but containing a scalar value for the respective
2157 parameter. If `param_config' is specified, none of the above
2158 vector/matrix-based settings may be used.
2160 Additionally, named parameters are allowed to be non-scalar real
2161 arrays. In this case, their dimensions are given by the setting
2162 `param_dims', a cell-array of dimension vectors, each containing
2163 at least two dimensions; if not given, dimensions are taken from
2164 the initial parameters, if these are given in a structure. Any
2165 vector-based settings or not structure-based linear constraints
2166 then must correspond to an order of parameters with all parameters
2167 reshaped to vectors and concatenated in the user-given order of
2168 parameter names. Structure-based settings or structure-based
2169 initial parameters must contain arrays with dimensions reshapable
2170 to those of the respective parameters.
2172 Description of backends
2176 A Levenberg/Marquardt-like optimizer, attempting to honour
2177 constraints throughout the course of optimization. This means that
2178 the initial parameters must not violate constraints (to find an
2179 initial feasible set of parameters, e.g. Octaves `sqp' can be
2180 used, by specifying an objective function which is constant or
2181 which returns the quadratic distance to the initial values). If the
2182 constraints need only be honoured in the result of the
2183 optimization, Octaves `sqp' may be preferable. The Hessian is
2184 either supplied by the user or is approximated by the BFGS
2187 Returned value CVG will be `2' or `3' for success and `0' or `-4'
2188 for failure (see above for meaning).
2190 Backend-specific defaults are: `MaxIter': 20, `fract_prec': `zeros
2191 (size (parameters))', `max_fract_change': `Inf' for all parameters.
2193 Interpretation of `Display': if set to `"iter"', currently only
2194 information on applying `max_fract_change' is printed.
2198 A simulated annealing (stochastic) optimizer, changing all
2199 parameters at once in a single step, so being suitable for
2200 non-bound constraints.
2202 No gradient or hessian of the objective function is used. The
2203 settings `MaxIter', `fract_prec', `TolFun', and `max_fract_change'
2206 Accepts the additional settings `T_init' (initial temperature,
2207 default 0.01), `T_min' (final temperature, default 1.0e-5), `mu_T'
2208 (factor of temperature decrease, default 1.005), `iters_fixed_T'
2209 (iterations within one temperature step, default 10),
2210 `max_rand_step' (column vector or structure-based configuration of
2211 maximum random steps for each parameter, default 0.005 * PIN),
2212 `stoch_regain_constr' (if `true', regain constraints after a
2213 random step, otherwise take new random value until constraints are
2214 met, default false), `trace_steps' (set field `trace' of OUTP with
2215 a matrix with a row for each step, first column iteration number,
2216 second column repeat number within iteration, third column value
2217 of objective function, rest columns parameter values, default
2218 false), and `siman_log' (set field `log' of OUTP with a matrix
2219 with a row for each iteration, first column temperature, second
2220 column value of objective function, rest columns numbers of tries
2221 with decrease, no decrease but accepted, and no decrease and
2224 Steps with increase `diff' of objective function are accepted if
2225 `rand (1) < exp (- diff / T)', where `T' is the temperature of the
2228 If regaining of constraints failed, optimization will be aborted
2229 and returned value of CVG will be `0'. Otherwise, CVG will be `1'.
2231 Interpretation of `Display': if set to `"iter"', an informational
2232 line is printed after each iteration.
2238 # name: <cell-element>
2242 Frontend for constrained nonlinear minimization of a scalar objective
2247 # name: <cell-element>
2254 # name: <cell-element>
2258 -- Function File: [P, RESID, CVG, OUTP] = nonlin_residmin (F, PIN)
2259 -- Function File: [P, RESID, CVG, OUTP] = nonlin_residmin (F, PIN,
2261 Frontend for nonlinear minimization of residuals returned by a
2264 The functions supplied by the user have a minimal interface; any
2265 additionally needed constants (e.g. observed values) can be
2266 supplied by wrapping the user functions into anonymous functions.
2268 The following description applies to usage with vector-based
2269 parameter handling. Differences in usage for structure-based
2270 parameter handling will be explained in a separate section below.
2272 F: function returning the array of residuals. It gets a column
2273 vector of real parameters as argument. In gradient determination,
2274 this function may be called with an informational second argument,
2275 whose content depends on the function for gradient determination.
2277 PIN: real column vector of initial parameters.
2279 SETTINGS: structure whose fields stand for optional settings
2280 referred to below. The fields can be set by `optimset()' with
2281 Octave versions 3.3.55 or greater; with older Octave versions, the
2282 fields must be set directly as structure-fields in the correct
2285 The returned values are the column vector of final parameters P,
2286 the final array of residuals RESID, an integer CVG indicating if
2287 and how optimization succeeded or failed, and a structure OUTP
2288 with additional information, curently with only one field: NITER,
2289 the number of iterations. CVG is greater than zero for success and
2290 less than or equal to zero for failure; its possible values depend
2291 on the used backend and currently can be `0' (maximum number of
2292 iterations exceeded), `2' (parameter change less than specified
2293 precision in two consecutive iterations), or `3' (improvement in
2294 objective function -- e.g. sum of squares -- less than specified).
2298 `Algorithm': String specifying the backend. Default:
2299 `"lm_svd_feasible"'. The latter is currently the only backend
2300 distributed with this package. It is described in a separate
2303 `dfdp': Function computing the jacobian of the residuals with
2304 respect to the parameters, assuming residuals are reshaped to a
2305 vector. Default: finite differences. Will be called with the column
2306 vector of parameters and an informational structure as arguments.
2307 The structure has the fields `f': value of residuals for current
2308 parameters, reshaped to a column vector, `fixed': logical vector
2309 indicating which parameters are not optimized, so these partial
2310 derivatives need not be computed and can be set to zero, `diffp',
2311 `diff_onesided', `lbound', `ubound': identical to the user
2312 settings of this name, `plabels': 1-dimensional cell-array of
2313 column-cell-arrays, each column with labels for all parameters,
2314 the first column contains the numerical indices of the parameters.
2315 The default jacobian function will call the model function with
2316 the second argument set with fields `f': as the `f' passed to the
2317 jacobian function, `plabels': cell-array of 1x1 cell-arrays with
2318 the entries of the column-cell-arrays of `plabels' as passed to
2319 the jacobian function corresponding to current parameter, `side':
2320 `0' for one-sided interval, `1' or `2', respectively, for the
2321 sides of a two-sided interval, and `parallel': logical scalar
2322 indicating parallel computation of partial derivatives.
2324 `diffp': column vector of fractional intervals (doubled for
2325 central intervals) supposed to be used by jacobian functions
2326 performing finite differencing. Default: `.001 * ones (size
2327 (parameters))'. The default jacobian function will use these as
2328 absolute intervals for parameters with value zero.
2330 `diff_onesided': logical column vector indicating that one-sided
2331 intervals should be used by jacobian functions performing finite
2332 differencing. Default: `false (size (parameters))'.
2334 `complex_step_derivative_f', `complex_step_derivative_inequc',
2335 `complex_step_derivative_equc': logical scalars, default: false.
2336 Estimate Jacobian of model function, general inequality
2337 constraints, and general equality constraints, respectively, with
2338 complex step derivative approximation. Use only if you know that
2339 your model function, function of general inequality constraints,
2340 or function of general equality constraints, respectively, is
2341 suitable for this. No user function for the respective Jacobian
2344 `cstep': scalar step size for complex step derivative
2345 approximation. Default: 1e-20.
2347 `fixed': logical column vector indicating which parameters should
2348 not be optimized, but kept to their inital value. Fixing is done
2349 independently of the backend, but the backend may choose to fix
2350 additional parameters under certain conditions.
2352 `lbound', `ubound': column vectors of lower and upper bounds for
2353 parameters. Default: `-Inf' and `+Inf', respectively. The bounds
2354 are non-strict, i.e. parameters are allowed to be exactly equal to
2355 a bound. The default jacobian function will respect bounds (but no
2356 further inequality constraints) in finite differencing.
2358 `inequc': Further inequality constraints. Cell-array containing up
2359 to four entries, two entries for linear inequality constraints
2360 and/or one or two entries for general inequality constraints.
2361 Either linear or general constraints may be the first entries, but
2362 the two entries for linear constraints must be adjacent and, if
2363 two entries are given for general constraints, they also must be
2364 adjacent. The two entries for linear constraints are a matrix (say
2365 `m') and a vector (say `v'), specifying linear inequality
2366 constraints of the form `m.' * parameters + v >= 0'. The first
2367 entry for general constraints must be a differentiable vector
2368 valued function (say `h'), specifying general inequality
2369 constraints of the form `h (p[, idx]) >= 0'; `p' is the column
2370 vector of optimized paraters and the optional argument `idx' is a
2371 logical index. `h' has to return the values of all constraints if
2372 `idx' is not given. It may choose to return only the indexed
2373 constraints if `idx' is given (so computation of the other
2374 constraints can be spared); in this case, the additional setting
2375 `inequc_f_idx' has to be set to `true'. In gradient determination,
2376 this function may be called with an informational third argument,
2377 whose content depends on the function for gradient determination.
2378 If a second entry for general inequality constraints is given, it
2379 must be a function computing the jacobian of the constraints with
2380 respect to the parameters. For this function, the description of
2381 `dfdp' above applies, with 2 exceptions: 1) it is called with 3
2382 arguments since it has an additional argument `idx' -- a logical
2383 index -- at second position, indicating which rows of the jacobian
2384 must be returned (if the function chooses to return only indexed
2385 rows, the additional setting `inequc_df_idx' has to be set to
2386 `true'). 2) the default jacobian function calls `h' with 3
2387 arguments, since the argument `idx' is also supplied. Note that
2388 specifying linear constraints as general constraints will generally
2389 waste performance, even if further, non-linear, general constraints
2392 `equc': Equality constraints. Specified the same way as inequality
2393 constraints (see `inequc').
2395 `cpiv': Function for complementary pivoting, usable in algorithms
2396 for constraints. Default: cpiv_bard. Only the default function is
2397 supplied with the package.
2399 `weights': Array of weights for the residuals. Dimensions must
2402 `TolFun': Minimum fractional improvement in objective function
2403 (e.g. sum of squares) in an iteration (abortion criterium).
2406 `MaxIter': Maximum number of iterations (abortion criterium).
2407 Default: backend-specific.
2409 `fract_prec': Column Vector, minimum fractional change of
2410 parameters in an iteration (abortion criterium if violated in two
2411 consecutive iterations). Default: backend-specific.
2413 `max_fract_change': Column Vector, enforced maximum fractional
2414 change in parameters in an iteration. Default: backend-specific.
2416 `Display': String indicating the degree of verbosity. Default:
2417 `"off"'. Possible values are currently `"off"' (no messages) and
2418 `"iter"' (some messages after each iteration). Support of this
2419 setting and its exact interpretation are backend-specific.
2421 `plot_cmd': Function enabling backend to plot results or
2422 intermediate results. Will be called with current computed
2423 residuals. Default: plot nothing.
2425 `debug': Logical scalar, default: `false'. Will be passed to the
2426 backend, which might print debugging information if true.
2428 Structure-based parameter handling
2430 The setting `param_order' is a cell-array with names of the
2431 optimized parameters. If not given, and initial parameters are a
2432 structure, all parameters in the structure are optimized. If
2433 initial parameters are a structure, it is an error if
2434 `param_order' is not given and there are any non-structure-based
2435 configuration items or functions.
2437 The initial parameters PIN can be given as a structure containing
2438 at least all fields named in `param_order'. In this case the
2439 returned parameters P will also be a structure.
2441 Each user-supplied function can be called with the argument
2442 containing the current parameters being a structure instead of a
2443 column vector. For this, a corresponding setting must be set to
2444 `true': `f_pstruct' (model function), `dfdp_pstruct' (jacobian of
2445 model function), `f_inequc_pstruct' (general inequality
2446 constraints), `df_inequc_pstruct' (jacobian of general inequality
2447 constraints), `f_equc_pstruct' (general equality constraints), and
2448 `df_equc_pstruct' (jacobian of general equality constraints). If a
2449 jacobian-function is configured in such a way, it must return the
2450 columns of the jacobian as fields of a structure under the
2451 respective parameter names.
2453 Similarly, for specifying linear constraints, instead of the matrix
2454 (called `m' above), a structure containing the rows of the matrix
2455 in fields under the respective parameter names can be given. In
2456 this case, rows containing only zeros need not be given.
2458 The vector-based settings `lbound', `ubound', `fixed', `diffp',
2459 `diff_onesided', `fract_prec', and `max_fract_change' can be
2460 replaced by the setting `param_config'. It is a structure that can
2461 contain fields named in `param_order'. For each such field, there
2462 may be subfields with the same names as the above vector-based
2463 settings, but containing a scalar value for the respective
2464 parameter. If `param_config' is specified, none of the above
2465 vector/matrix-based settings may be used.
2467 Additionally, named parameters are allowed to be non-scalar real
2468 arrays. In this case, their dimensions are given by the setting
2469 `param_dims', a cell-array of dimension vectors, each containing
2470 at least two dimensions; if not given, dimensions are taken from
2471 the initial parameters, if these are given in a structure. Any
2472 vector-based settings or not structure-based linear constraints
2473 then must correspond to an order of parameters with all parameters
2474 reshaped to vectors and concatenated in the user-given order of
2475 parameter names. Structure-based settings or structure-based
2476 initial parameters must contain arrays with dimensions reshapable
2477 to those of the respective parameters.
2479 Description of backends (currently only one)
2483 A Levenberg/Marquardt algorithm using singular value decomposition
2484 and featuring constraints which must be met by the initial
2485 parameters and are attempted to be kept met throughout the
2488 Parameters with identical lower and upper bounds will be fixed.
2490 Returned value CVG will be `0', `2', or `3'.
2492 Backend-specific defaults are: `MaxIter': 20, `fract_prec': `zeros
2493 (size (parameters))', `max_fract_change': `Inf' for all parameters.
2495 Interpretation of `Display': if set to `"iter"', currently
2496 `plot_cmd' is evaluated for each iteration, and some further
2497 diagnostics may be printed.
2499 Specific option: `lm_svd_feasible_alt_s': if falling back to
2500 nearly gradient descent, do it more like original
2501 Levenberg/Marquardt method, with descent in each gradient
2502 component; for testing only.
2504 See also: nonlin_curvefit
2510 # name: <cell-element>
2514 Frontend for nonlinear minimization of residuals returned by a model
2519 # name: <cell-element>
2526 # name: <cell-element>
2530 -- Function File: XMIN = nrm(F,X0)
2531 Using X0 as a starting point find a minimum of the scalar function
2532 F. The Newton-Raphson method is used.
2537 # name: <cell-element>
2541 Using X0 as a starting point find a minimum of the scalar function F.
2545 # name: <cell-element>
2552 # name: <cell-element>
2556 Problems for testing optimizers. Documentation is in the code.
2560 # name: <cell-element>
2564 Problems for testing optimizers.
2568 # name: <cell-element>
2575 # name: <cell-element>
2579 opt = optimset_compat (...) - manipulate m*tlab-style options structure
2581 This function returns a m*tlab-style options structure that can be used
2582 with the fminunc() function. 'optimset_compat' has been deprecated in
2583 favor of 'optimset', which is now part of core Octave. This function
2584 will possibly be removed from future versions of the 'optim' package.
2586 INPUT : Input consist in one or more structs followed by option-value
2587 pairs. The option that can be passed are those of m*tlab's 'optimset'.
2588 Whether fminunc() accepts them is another question (see fminunc()).
2590 Two extra options are supported which indicate how to use directly octave
2591 optimization tools (such as minimize() and other backends):
2593 "MinEquiv", [on|off] : Tell 'fminunc()' not to minimize 'fun', but
2594 instead return the option passed to minimize().
2596 "Backend", [on|off] : Tell 'fminunc()' not to minimize 'fun', but
2597 instead return the [backend, opt], the name of the
2598 backend optimization function that is used and the
2599 optional arguments that will be passed to it. See
2600 the 'backend' option of minimize().
2604 # name: <cell-element>
2608 opt = optimset_compat (.
2612 # name: <cell-element>
2619 # name: <cell-element>
2623 ex = poly_2_ex (l, f) - Extremum of a 1-var deg-2 polynomial
2625 l : 3 : Values of variable at which polynomial is known.
2626 f : 3 : f(i) = Value of the degree-2 polynomial at l(i).
2628 ex : 1 : Value for which f reaches its extremum
2630 Assuming that f(i) = a*l(i)^2 + b*l(i) + c = P(l(i)) for some a, b, c,
2631 ex is the extremum of the polynome P.
2636 # name: <cell-element>
2640 ex = poly_2_ex (l, f) - Extremum of a 1-var deg-2 polynomial
2645 # name: <cell-element>
2652 # name: <cell-element>
2656 [y,dy] = polyconf(p,x,s)
2658 Produce prediction intervals for the fitted y. The vector p
2659 and structure s are returned from polyfit or wpolyfit. The
2660 x values are where you want to compute the prediction interval.
2662 polyconf(...,['ci'|'pi'])
2664 Produce a confidence interval (range of likely values for the
2665 mean at x) or a prediction interval (range of likely values
2666 seen when measuring at x). The prediction interval tells
2667 you the width of the distribution at x. This should be the same
2668 regardless of the number of measurements you have for the value
2669 at x. The confidence interval tells you how well you know the
2670 mean at x. It should get smaller as you increase the number of
2671 measurements. Error bars in the physical sciences usually show
2672 a 1-alpha confidence value of erfc(1/sqrt(2)), representing
2673 one standandard deviation of uncertainty in the mean.
2675 polyconf(...,1-alpha)
2677 Control the width of the interval. If asking for the prediction
2678 interval 'pi', the default is .05 for the 95% prediction interval.
2679 If asking for the confidence interval 'ci', the default is
2680 erfc(1/sqrt(2)) for a one standard deviation confidence interval.
2683 [p,s] = polyfit(x,y,1);
2684 xf = linspace(x(1),x(end),150);
2685 [yf,dyf] = polyconf(p,xf,s,'ci');
2686 plot(xf,yf,'g-;fit;',xf,yf+dyf,'g.;;',xf,yf-dyf,'g.;;',x,y,'xr;data;');
2687 plot(x,y-polyval(p,x),';residuals;',xf,dyf,'g-;;',xf,-dyf,'g-;;');
2691 # name: <cell-element>
2695 [y,dy] = polyconf(p,x,s)
2700 # name: <cell-element>
2707 # name: <cell-element>
2711 function [A,REF,HMAX,H,R,EQUAL] = polyfitinf(M,N,K,X,Y,EPSH,MAXIT,REF0)
2713 Best polynomial approximation in discrete uniform norm
2717 M : degree of the fitting polynomial
2718 N : number of data points
2719 X(N) : x-coordinates of data points
2720 Y(N) : y-coordinates of data points
2721 K : character of the polynomial:
2722 K = 0 : mixed parity polynomial
2723 K = 1 : odd polynomial ( X(1) must be > 0 )
2724 K = 2 : even polynomial ( X(1) must be >= 0 )
2725 EPSH : tolerance for leveling. A useful value for 24-bit
2726 mantissa is EPSH = 2.0E-7
2727 MAXIT : upper limit for number of exchange steps
2728 REF0(M2): initial alternating set ( N-vector ). This is an
2729 OPTIONAL argument. The length M2 is given by:
2730 M2 = M + 2 , if K = 0
2731 M2 = integer part of (M+3)/2 , if K = 1
2732 M2 = 2 + M/2 (M must be even) , if K = 2
2736 A : polynomial coefficients of the best approximation
2737 in order of increasing powers:
2738 p*(x) = A(1) + A(2)*x + A(3)*x^2 + ...
2739 REF : selected alternating set of points
2740 HMAX : maximum deviation ( uniform norm of p* - f )
2741 H : pointwise approximation errors
2742 R : total number of iterations
2743 EQUAL : success of failure of algorithm
2745 EQUAL=0 : convergence not acheived
2746 EQUAL=-1: input error
2747 EQUAL=-2: algorithm failure
2749 Relies on function EXCH, provided below.
2752 M = 5; N = 10000; K = 0; EPSH = 10^-12; MAXIT = 10;
2753 X = linspace(-1,1,N); % uniformly spaced nodes on [-1,1]
2754 k=1; Y = abs(X).^k; % the function Y to approximate
2755 [A,REF,HMAX,H,R,EQUAL] = polyfitinf(M,N,K,X,Y,EPSH,MAXIT);
2756 p = polyval(A,X); plot(X,Y,X,p) % p is the best approximation
2758 Note: using an even value of M, e.g., M=2, in the example above, makes
2759 the algorithm to fail with EQUAL=-2, because of collocation, which
2760 appears because both the appriximating function and the polynomial are
2761 even functions. The way aroung it is to approximate only the right half
2762 of the function, setting K = 2 : even polynomial. For example:
2764 N = 10000; K = 2; EPSH = 10^-12; MAXIT = 10; X = linspace(0,1,N);
2766 k = 2*i-1; Y = abs(X).^k;
2769 [~,~,HMAX] = polyfitinf(M,N,K,X,Y,EPSH,MAXIT);
2770 approxerror(i,j) = HMAX;
2773 disp('Table 3.1 from Approximation theory and methods, M.J.D.POWELL, p. 27');
2776 disp(' '); format short g;
2777 disp([(2.^(1:4))' approxerror']);
2781 Computation of the polynomial that best approximates the data (X,Y)
2782 in the discrete uniform norm, i.e. the polynomial with the minimum
2783 value of max{ | p(x_i) - y_i | , x_i in X } . That polynomial, also
2784 known as minimax polynomial, is obtained by the exchange algorithm,
2785 a finite iterative process requiring, at most,
2787 ( ) iterations ( usually p = M + 2. See also function EXCH ).
2789 since this number can be very large , the routine may not converge
2790 within MAXIT iterations . The other possibility of failure occurs
2791 when there is insufficient floating point precision for the input
2794 CREDITS: This routine was developed and modified as
2795 computer assignments in Approximation Theory courses by
2796 Prof. Andrew Knyazev, University of Colorado Denver, USA.
2798 Team Fall 98 (Revision 1.0):
2799 Chanchai Aniwathananon
2804 Team Spring 11 (Revision 1.1): Manuchehr Aminian
2806 The algorithm and the comments are based on a FORTRAN code written
2807 by Joseph C. Simpson. The code is available on Netlib repository:
2808 http://www.netlib.org/toms/501
2809 See also: Communications of the ACM, V14, pp.355-356(1971)
2813 1) A may contain the collocation polynomial
2814 2) If MAXIT is exceeded, REF contains a new reference set
2815 3) M, EPSH and REF can be altered during the execution
2816 4) To keep consistency to the original code , EPSH can be
2817 negative. However, the use of REF0 is *NOT* determined by
2818 EPSH< 0, but only by its inclusion as an input parameter.
2820 Some parts of the code can still take advantage of vectorization.
2822 Revision 1.0 from 1998 is a direct human translation of
2823 the FORTRAN code http://www.netlib.org/toms/501
2824 Revision 1.1 is a clean-up and technical update.
2825 Tested on MATLAB Version 7.11.0.584 (R2010b) and
2826 GNU Octave Version 3.2.4
2830 # name: <cell-element>
2834 function [A,REF,HMAX,H,R,EQUAL] = polyfitinf(M,N,K,X,Y,EPSH,MAXIT,REF0)
2839 # name: <cell-element>
2846 # name: <cell-element>
2850 -- Function File: [ P, OBJ_VALUE, CONVERGENCE, ITERS, NEVS] = powell
2852 Multidimensional minimization (direction-set method). Implements a
2853 direction-set (Powell's) method for multidimensional minimization
2854 of a function without calculation of the gradient [1, 2]
2859 * F: name of function to minimize (string or handle), which
2860 should accept one input variable (see example for how to pass
2861 on additional input arguments)
2863 * P0: An initial value of the function argument to minimize
2865 * OPTIONS: an optional structure, which can be generated by
2866 optimset, with some or all of the following fields:
2867 - MaxIter: maximum iterations (positive integer, or -1
2868 or Inf for unlimited (default))
2870 - TolFun: minimum amount by which function value must
2871 decrease in each iteration to continue (default is 1E-8)
2873 - MaxFunEvals: maximum function evaluations (positive
2874 integer, or -1 or Inf for unlimited (default))
2876 - SearchDirections: an n*n matrix whose columns contain
2877 the initial set of (presumably orthogonal) directions to
2878 minimize along, where n is the number of elements in the
2879 argument to be minimized for; or an n*1 vector of
2880 magnitudes for the initial directions (defaults to the
2881 set of unit direction vectors)
2886 y = @(x, s) x(1) ^ 2 + x(2) ^ 2 + s;
2887 o = optimset('MaxIter', 100, 'TolFun', 1E-10);
2889 [x_optim, y_min, conv, iters, nevs] = powell(@(x) y(x, s), [1 0.5], o); %pass y wrapped in an anonymous function so that all other arguments to y, which are held constant, are set
2890 %should return something like x_optim = [4E-14 3E-14], y_min = 1, conv = 1, iters = 2, nevs = 24
2896 * P: the minimizing value of the function argument
2898 * OBJ_VALUE: the value of F() at P
2900 * CONVERGENCE: 1 if normal convergence, 0 if not
2902 * ITERS: number of iterations performed
2904 * NEVS: number of function evaluations
2909 1. Powell MJD (1964), An efficient method for finding the
2910 minimum of a function of several variables without
2911 calculating derivatives, `Computer Journal', 7 :155-162
2913 2. Press, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP
2914 (1992). `Numerical Recipes in Fortran: The Art of Scientific
2915 Computing' (2nd Ed.). New York: Cambridge University Press
2921 # name: <cell-element>
2925 Multidimensional minimization (direction-set method).
2929 # name: <cell-element>
2936 # name: <cell-element>
2940 -- Function File: INFO = residmin_stat (F, P, SETTINGS)
2941 Frontend for computation of statistics for a residual-based
2944 SETTINGS is a structure whose fields can be set by `optimset' with
2945 Octave versions 3.3.55 or greater; with older Octave versions, the
2946 fields must be set directly and in the correct case. With SETTINGS
2947 the computation of certain statistics is requested by setting the
2948 fields `ret_<name_of_statistic>' to `true'. The respective
2949 statistics will be returned in a structure as fields with name
2950 `<name_of_statistic>'. Depending on the requested statistic and on
2951 the additional information provided in SETTINGS, F and P may be
2952 empty. Otherwise, F is the model function of an optimization (the
2953 interface of F is described e.g. in `nonlin_residmin', please see
2954 there), and P is a real column vector with parameters resulting
2955 from the same optimization.
2957 Currently, the following statistics (or general information) can be
2960 `dfdp': Jacobian of model function with respect to parameters.
2962 `covd': Covariance matrix of data (typically guessed by applying a
2963 factor to the covariance matrix of the residuals).
2965 `covp': Covariance matrix of final parameters.
2967 `corp': Correlation matrix of final parameters.
2971 The functionality of the interface is similar to
2972 `nonlin_residmin'. In particular, structure-based, possibly
2973 non-scalar, parameters and flagging parameters as fixed are
2974 possible. The following settings have the same meaning as in
2975 `nonlin_residmin' (please refer to there): `param_order',
2976 `param_dims', `f_pstruct', `dfdp_pstruct', `diffp',
2977 `diff_onesided', `complex_step_derivative', `cstep', `fixed', and
2978 `weights'. Similarly, `param_config' can be used, but only with
2979 fields corresponding to the settings `fixed', `diffp', and
2982 `dfdp' can be set in the same way as in `nonlin_residmin', but
2983 alternatively may already contain the computed Jacobian of the
2984 model function at the final parameters in matrix- or
2985 structure-form. Users may pass information on the result of the
2986 optimization in `residuals' (self-explaining) and `covd'
2987 (covariance matrix of data). In many cases the type of objective
2988 function of the optimization must be specified in `objf';
2989 currently, there is only a backend for the type "wls" (weighted
2992 Backend-specific information
2994 The backend for `objf == "wls"' (currently the only backend)
2995 computes `cord' (due to user request or as a prerequisite for
2996 `covp' and `corp') as a diagonal matrix by assuming that the
2997 variances of data points are proportional to the reciprocal of the
2998 squared `weights' and guessing the factor of proportionality from
2999 the residuals. If `covp' is not defined (e.g. because the Jacobian
3000 has no full rank), it makes an attempt to still compute its
3001 uniquely defined elements, if any, and to find the additional
3002 defined elements (being `1' or `-1'), if any, in `corp'.
3004 See also: curvefit_stat
3010 # name: <cell-element>
3014 Frontend for computation of statistics for a residual-based
3019 # name: <cell-element>
3026 # name: <cell-element>
3030 Rosenbrock function - used to create example obj. fns.
3032 Function value and gradient vector of the rosenbrock function
3033 The minimizer is at the vector (1,1,..,1),
3034 and the minimized value is 0.
3038 # name: <cell-element>
3042 Rosenbrock function - used to create example obj.
3046 # name: <cell-element>
3053 # name: <cell-element>
3061 # name: <cell-element>
3070 # name: <cell-element>
3077 # name: <cell-element>
3081 Plain run, just to make sure ######################################
3082 Minimum wrt 'x' is y0
3086 # name: <cell-element>
3090 Plain run, just to make sure ######################################
3095 # name: <cell-element>
3102 # name: <cell-element>
3106 [x,v,niter] = feval (optim_func, "testfunc","dtestf", xinit);
3110 # name: <cell-element>
3114 [x,v,niter] = feval (optim_func, "testfunc","dtestf", xinit);
3119 # name: <cell-element>
3126 # name: <cell-element>
3130 [xlev,vlev,nlev] = feval(optim_func, "ff", "dff", xinit) ;
3134 # name: <cell-element>
3138 [xlev,vlev,nlev] = feval(optim_func, "ff", "dff", xinit) ;
3143 # name: <cell-element>
3150 # name: <cell-element>
3154 [xlev,vlev,nlev] = feval (optim_func, "ff", "dff", xinit, "extra", extra) ;
3155 [xlev,vlev,nlev] = feval \
3156 (optim_func, "ff", "dff", list (xinit, obsmat, obses));
3160 # name: <cell-element>
3164 [xlev,vlev,nlev] = feval (optim_func, "ff", "dff", xinit, "extra", extra) ;
3169 # name: <cell-element>
3176 # name: <cell-element>
3180 Plain run, just to make sure ######################################
3181 Minimum wrt 'x' is y0
3182 [xlev,vlev,nlev] = feval (optim_func, "ff", "dff", {x0,y0,1});
3187 # name: <cell-element>
3191 Plain run, just to make sure ######################################
3196 # name: <cell-element>
3203 # name: <cell-element>
3207 Plain run, just to make sure ######################################
3208 Minimum wrt 'x' is y0
3212 # name: <cell-element>
3216 Plain run, just to make sure ######################################
3221 # name: <cell-element>
3225 test_nelder_mead_min_1
3228 # name: <cell-element>
3232 Use vanilla nelder_mead_min
3236 # name: <cell-element>
3240 Use vanilla nelder_mead_min
3245 # name: <cell-element>
3249 test_nelder_mead_min_2
3252 # name: <cell-element>
3257 Test using volume #################################################
3261 # name: <cell-element>
3266 Test using volume #################################################
3271 # name: <cell-element>
3278 # name: <cell-element>
3286 # name: <cell-element>
3295 # name: <cell-element>
3302 # name: <cell-element>
3306 -- Function File: vfzero (FUN, X0)
3307 -- Function File: vfzero (FUN, X0, OPTIONS)
3308 -- Function File: [X, FVAL, INFO, OUTPUT] = vfzero (...)
3309 A variant of `fzero'. Finds a zero of a vector-valued multivariate
3310 function where each output element only depends on the input
3311 element with the same index (so the Jacobian is diagonal).
3313 FUN should be a handle or name of a function returning a column
3314 vector. X0 should be a two-column matrix, each row specifying two
3315 points which bracket a zero of the respective output element of
3318 If X0 is a single-column matrix then several nearby and distant
3319 values are probed in an attempt to obtain a valid bracketing. If
3320 this is not successful, the function fails. OPTIONS is a structure
3321 specifying additional options. Currently, `vfzero' recognizes
3322 these options: `"FunValCheck"', `"OutputFcn"', `"TolX"',
3323 `"MaxIter"', `"MaxFunEvals"'. For a description of these options,
3324 see *note optimset: doc-optimset.
3326 On exit, the function returns X, the approximate zero and FVAL,
3327 the function value thereof. INFO is a column vector of exit flags
3328 that can have these values:
3330 * 1 The algorithm converged to a solution.
3332 * 0 Maximum number of iterations or function evaluations has
3335 * -1 The algorithm has been terminated from user output
3338 * -5 The algorithm may have converged to a singular point.
3340 OUTPUT is a structure containing runtime information about the
3341 `fzero' algorithm. Fields in the structure are:
3343 * iterations Number of iterations through loop.
3345 * nfev Number of function evaluations.
3347 * bracketx A two-column matrix with the final bracketing of the
3348 zero along the x-axis.
3350 * brackety A two-column matrix with the final bracketing of the
3351 zero along the y-axis.
3353 See also: optimset, fsolve
3359 # name: <cell-element>
3363 A variant of `fzero'.
3367 # name: <cell-element>
3374 # name: <cell-element>
3378 -- Function File: [P, S] = wpolyfit (X, Y, DY, N)
3379 Return the coefficients of a polynomial P(X) of degree N that
3380 minimizes `sumsq (p(x(i)) - y(i))', to best fit the data in the
3381 least squares sense. The standard error on the observations Y if
3382 present are given in DY.
3384 The returned value P contains the polynomial coefficients suitable
3385 for use in the function polyval. The structure S returns
3386 information necessary to compute uncertainty in the model.
3388 To compute the predicted values of y with uncertainty use
3389 [y,dy] = polyconf(p,x,s,'ci');
3390 You can see the effects of different confidence intervals and
3391 prediction intervals by calling the wpolyfit internal plot
3392 function with your fit:
3393 feval('wpolyfit:plt',x,y,dy,p,s,0.05,'pi')
3394 Use DY=[] if uncertainty is unknown.
3396 You can use a chi^2 test to reject the polynomial fit:
3397 p = 1-chi2cdf(s.normr^2,s.df);
3398 p is the probability of seeing a chi^2 value higher than that which
3399 was observed assuming the data are normally distributed around the
3400 fit. If p < 0.01, you can reject the fit at the 1% level.
3402 You can use an F test to determine if a higher order polynomial
3404 [poly1,S1] = wpolyfit(x,y,dy,n);
3405 [poly2,S2] = wpolyfit(x,y,dy,n+1);
3406 F = (S1.normr^2 - S2.normr^2)/(S1.df-S2.df)/(S2.normr^2/S2.df);
3407 p = 1-f_cdf(F,S1.df-S2.df,S2.df);
3408 p is the probability of observing the improvement in chi^2 obtained
3409 by adding the extra parameter to the fit. If p < 0.01, you can
3410 reject the lower order polynomial at the 1% level.
3412 You can estimate the uncertainty in the polynomial coefficients
3414 dp = sqrt(sumsq(inv(s.R'))'/s.df)*s.normr;
3415 but the high degree of covariance amongst them makes this a
3416 questionable operation.
3418 -- Function File: [P, S, MU] = wpolyfit (...)
3419 If an additional output `mu = [mean(x),std(x)]' is requested then
3420 the X values are centered and normalized prior to computing the
3421 fit. This will give more stable numerical results. To compute a
3422 predicted Y from the returned model use `y = polyval(p,
3425 -- Function File: wpolyfit (...)
3426 If no output arguments are requested, then wpolyfit plots the data,
3427 the fitted line and polynomials defining the standard error range.
3430 x = linspace(0,4,20);
3431 dy = (1+rand(size(x)))/2;
3432 y = polyval([2,3,1],x) + dy.*randn(size(x));
3435 -- Function File: wpolyfit (..., 'origin')
3436 If 'origin' is specified, then the fitted polynomial will go
3437 through the origin. This is generally ill-advised. Use with
3440 Hocking, RR (2003). Methods and Applications of Linear Models.
3441 New Jersey: John Wiley and Sons, Inc.
3444 See also: polyfit, polyconf
3449 # name: <cell-element>
3453 Return the coefficients of a polynomial P(X) of degree N that minimizes
3458 # name: <cell-element>
3465 # name: <cell-element>
3469 [ret1, ret2] = wrap_f_dfdp (f, dfdp, varargin)
3471 f and dftp should be the objective function (or "model function" in
3472 curve fitting) and its jacobian, respectively, of an optimization
3473 problem. ret1: f (varagin{:}), ret2: dfdp (varargin{:}). ret2 is
3474 only computed if more than one output argument is given. This
3475 manner of calling f and dfdp is needed by some optimization
3480 # name: <cell-element>
3484 [ret1, ret2] = wrap_f_dfdp (f, dfdp, varargin)
3489 # name: <cell-element>
3496 # name: <cell-element>
3500 [x,s] = wsolve(A,y,dy)
3502 Solve a potentially over-determined system with uncertainty in
3507 Use QR decomposition for increased accuracy. Estimate the
3508 uncertainty for the solution from the scatter in the data.
3510 The returned structure s contains
3512 normr = sqrt( A x - y ), weighted by dy
3513 R such that R'R = A'A
3514 df = n-p, n = rows of A, p = columns of A
3516 See polyconf for details on how to use s to compute dy.
3517 The covariance matrix is inv(R'*R). If you know that the
3518 parameters are independent, then uncertainty is given by
3519 the diagonal of the covariance matrix, or
3521 dx = sqrt(N*sumsq(inv(s.R'))')
3523 where N = normr^2/df, or N = 1 if df = 0.
3525 Example 1: weighted system
3527 A=[1,2,3;2,1,3;1,1,1]; xin=[1;2;3];
3528 dy=[0.2;0.01;0.1]; y=A*xin+randn(size(dy)).*dy;
3529 [x,s] = wsolve(A,y,dy);
3530 dx = sqrt(sumsq(inv(s.R'))');
3533 Example 2: weighted overdetermined system y = x1 + 2*x2 + 3*x3 + e
3535 A = fullfact([3,3,3]); xin=[1;2;3];
3536 y = A*xin; dy = rand(size(y))/50; y+=dy.*randn(size(y));
3537 [x,s] = wsolve(A,y,dy);
3538 dx = s.normr*sqrt(sumsq(inv(s.R'))'/s.df);
3541 Note there is a counter-intuitive result that scaling the
3542 uncertainty in the data does not affect the uncertainty in
3543 the fit. Indeed, if you perform a monte carlo simulation
3544 with x,y datasets selected from a normal distribution centered
3545 on y with width 10*dy instead of dy you will see that the
3546 variance in the parameters indeed increases by a factor of 100.
3547 However, if the error bars really do increase by a factor of 10
3548 you should expect a corresponding increase in the scatter of
3549 the data, which will increase the variance computed by the fit.
3553 # name: <cell-element>
3557 [x,s] = wsolve(A,y,dy)