]> Creatis software - CreaPhase.git/blob - octave_packages/optim-1.2.0/doc-cache
Add a useful package (from Source forge) for octave
[CreaPhase.git] / octave_packages / optim-1.2.0 / doc-cache
1 # Created by Octave 3.6.2, Tue Jun 12 20:02:41 2012 UTC <root@t61>
2 # name: cache
3 # type: cell
4 # rows: 3
5 # columns: 58
6 # name: <cell-element>
7 # type: sq_string
8 # elements: 1
9 # length: 16
10 LinearRegression
11
12
13 # name: <cell-element>
14 # type: sq_string
15 # elements: 1
16 # length: 717
17  general linear regression
18
19  [p,y_var,r,p_var]=LinearRegression(F,y)
20  [p,y_var,r,p_var]=LinearRegression(F,y,weight)
21
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)
25
26  parameters
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
31
32  return values
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
37
38
39
40 # name: <cell-element>
41 # type: sq_string
42 # elements: 1
43 # length: 27
44  general linear regression
45
46
47
48
49 # name: <cell-element>
50 # type: sq_string
51 # elements: 1
52 # length: 6
53 adsmax
54
55
56 # name: <cell-element>
57 # type: sq_string
58 # elements: 1
59 # length: 1909
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.
64         Output arguments:
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,...).
85      Reference:
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.
91
92
93
94 # name: <cell-element>
95 # type: sq_string
96 # elements: 1
97 # length: 69
98 ADSMAX  Alternating directions method for direct search optimization.
99
100
101
102 # name: <cell-element>
103 # type: sq_string
104 # elements: 1
105 # length: 7
106 battery
107
108
109 # name: <cell-element>
110 # type: sq_string
111 # elements: 1
112 # length: 474
113  battery.m: repeatedly call bfgs using a battery of 
114  start values, to attempt to find global min
115  of a nonconvex function
116
117  INPUTS:
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
124
125  OUTPUT: theta - the best value found - NOT iterated to convergence
126
127
128
129 # name: <cell-element>
130 # type: sq_string
131 # elements: 1
132 # length: 9
133  battery.
134
135
136
137 # name: <cell-element>
138 # type: sq_string
139 # elements: 1
140 # length: 7
141 bfgsmin
142
143
144 # name: <cell-element>
145 # type: sq_string
146 # elements: 1
147 # length: 1747
148  bfgsmin: bfgs or limited memory bfgs minimization of function
149
150  Usage: [x, obj_value, convergence, iters] = bfgsmin(f, args, control)
151
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
157  for methods.
158
159  Arguments:
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
163         MUST be a vector
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))
168         * elem 2: verbosity
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
182
183  Returns:
184  * x: the minimizer
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
188
189  Example: see bfgsmin_example.m
190
191
192
193 # name: <cell-element>
194 # type: sq_string
195 # elements: 1
196 # length: 63
197  bfgsmin: bfgs or limited memory bfgs minimization of function
198
199
200
201
202 # name: <cell-element>
203 # type: sq_string
204 # elements: 1
205 # length: 15
206 bfgsmin_example
207
208
209 # name: <cell-element>
210 # type: sq_string
211 # elements: 1
212 # length: 16
213  initial values
214
215
216
217 # name: <cell-element>
218 # type: sq_string
219 # elements: 1
220 # length: 16
221  initial values
222
223
224
225
226 # name: <cell-element>
227 # type: sq_string
228 # elements: 1
229 # length: 14
230 brent_line_min
231
232
233 # name: <cell-element>
234 # type: sq_string
235 # elements: 1
236 # length: 1186
237  -- Function File: [S,V,N] brent_line_min ( F,DF,ARGS,CTL )
238      Line minimization of f along df
239
240      Finds minimum of f on line  x0 + dx*w | a < w < b  by bracketing.
241      a and b are passed through argument ctl.
242
243 Arguments
244 ---------
245
246         * F     : string : Name of function. Must return a real value
247
248         * ARGS  : cell   : Arguments passed to f or RxC    : f's only
249           argument. x0 must be at ARGS{ CTL(2) }
250
251         * CTL   : 5      : (optional) Control variables, described
252           below.
253
254 Returned values
255 ---------------
256
257         * S   : 1        : Minimum is at x0 + s*dx
258
259         * V   : 1        : Value of f at x0 + s*dx
260
261         * NEV : 1        : Number of function evaluations
262
263 Control Variables
264 -----------------
265
266         * CTL(1)       : Upper bound for error on s
267           Default=sqrt(eps)
268
269         * CTL(2)       : Position of minimized argument in args
270           Default= 1
271
272         * CTL(3)       : Maximum number of function evaluations
273           Default= inf
274
275         * CTL(4)       : a
276           Default=-inf
277
278         * CTL(5)       : b
279           Default= inf
280
281      Default values will be used if ctl is not passed or if nan values
282 are given.
283
284
285
286
287 # name: <cell-element>
288 # type: sq_string
289 # elements: 1
290 # length: 32
291 Line minimization of f along df
292
293
294
295
296 # name: <cell-element>
297 # type: sq_string
298 # elements: 1
299 # length: 6
300 cauchy
301
302
303 # name: <cell-element>
304 # type: sq_string
305 # elements: 1
306 # length: 766
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.
313      For example,
314
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.
318
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)
321
322
323
324
325 # name: <cell-element>
326 # type: sq_string
327 # elements: 1
328 # length: 80
329 Return the Taylor coefficients and numerical differentiation of a
330 function F for
331
332
333
334 # name: <cell-element>
335 # type: sq_string
336 # elements: 1
337 # length: 5
338 cdiff
339
340
341 # name: <cell-element>
342 # type: sq_string
343 # elements: 1
344 # length: 1377
345  c = cdiff (func,wrt,N,dfunc,stack,dx) - Code for num. differentiation
346    = "function df = dfunc (var1,..,dvar,..,varN) .. endfunction
347  
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
350  argument.
351
352  The derivatives are obtained by symmetric finite difference.
353
354  dfunc()'s return value is in the same format as that of  ndiff()
355
356  func  : string : name of the function to differentiate
357
358  wrt   : int    : position, in argument list, of the differentiation
359                   variable.                                Default:1
360
361  N     : int    : total number of arguments taken by 'func'. 
362                   If N=inf, dfunc will take variable argument list.
363                                                          Default:wrt
364
365  dfunc : string : Name of the octave function that returns the
366                    derivatives.                   Default:['d',func]
367
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:''
372
373  dx    : real   : Step used in the symmetric difference scheme.
374                                                   Default:10*sqrt(eps)
375
376  See also : ndiff, eval, todisk
377
378
379
380 # name: <cell-element>
381 # type: sq_string
382 # elements: 1
383 # length: 54
384  c = cdiff (func,wrt,N,dfunc,stack,dx) - Code for num.
385
386
387
388 # name: <cell-element>
389 # type: sq_string
390 # elements: 1
391 # length: 6
392 cg_min
393
394
395 # name: <cell-element>
396 # type: sq_string
397 # elements: 1
398 # length: 2479
399  -- Function File: [X0,V,NEV] cg_min ( F,DF,ARGS,CTL )
400      NonLinear Conjugate Gradient method to minimize function F.
401
402 Arguments
403 ---------
404
405         * F   : string   : Name of function. Return a real value
406
407         * DF  : string   : Name of f's derivative. Returns a (R*C) x 1
408           vector
409
410         * ARGS: cell     : Arguments passed to f.
411         * CTL   : 5-vec    : (Optional) Control variables, described
412           below
413
414 Returned values
415 ---------------
416
417         * X0    : matrix   : Local minimum of f
418
419         * V     : real     : Value of f in x0
420
421         * NEV   : 1 x 2    : Number of evaluations of f and of df
422
423 Control Variables
424 -----------------
425
426         * CTL(1)       : 1 or 2 : Select stopping criterion amongst :
427
428         * CTL(1)==0    : Default value
429
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).
434
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
438           the last iteration.
439
440         * CTL(2)       : Threshold used in stopping tests.
441           Default=10*eps
442
443         * CTL(2)==0    : Default value
444
445         * CTL(3)       : Position of the minimized argument in args
446           Default=1
447
448         * CTL(3)==0    : Default value
449
450         * CTL(4)       : Maximum number of function evaluations
451           Default=inf
452
453         * CTL(4)==0    : Default value
454
455         * CTL(5)       : Type of optimization:
456
457         * CTL(5)==1    : "Fletcher-Reves" method
458
459         * CTL(5)==2    : "Polak-Ribiere" (Default)
460
461         * CTL(5)==3    : "Hestenes-Stiefel" method
462
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.
465
466 Example:
467 --------
468
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);
472 endfunction
473      ll = { [10; 2; 3] };
474      ctl(5) = 3;
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
478
479      See also: bfgsmin,
480 http://en.wikipedia.org/wiki/Nonlinear_conjugate_gradient
481
482
483
484
485
486 # name: <cell-element>
487 # type: sq_string
488 # elements: 1
489 # length: 59
490 NonLinear Conjugate Gradient method to minimize function F.
491
492
493
494 # name: <cell-element>
495 # type: sq_string
496 # elements: 1
497 # length: 9
498 cpiv_bard
499
500
501 # name: <cell-element>
502 # type: sq_string
503 # elements: 1
504 # length: 1197
505  [lb, idx, ridx, mv] = cpiv_bard (v, m[, incl])
506
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.).
525
526
527
528 # name: <cell-element>
529 # type: sq_string
530 # elements: 1
531 # length: 48
532  [lb, idx, ridx, mv] = cpiv_bard (v, m[, incl])
533
534
535
536
537 # name: <cell-element>
538 # type: sq_string
539 # elements: 1
540 # length: 13
541 curvefit_stat
542
543
544 # name: <cell-element>
545 # type: sq_string
546 # elements: 1
547 # length: 749
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.
551
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.
560
561      See also: residmin_stat
562
563
564
565
566
567 # name: <cell-element>
568 # type: sq_string
569 # elements: 1
570 # length: 80
571 Frontend for computation of statistics for fitting of values, computed
572 by a mode
573
574
575
576 # name: <cell-element>
577 # type: sq_string
578 # elements: 1
579 # length: 6
580 d2_min
581
582
583 # name: <cell-element>
584 # type: sq_string
585 # elements: 1
586 # length: 2782
587  [x,v,nev,h,args] = d2_min(f,d2f,args,ctl,code) - Newton-like minimization
588
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".
593
594  ARGUMENTS :
595  f    : string : Cost function's name
596
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) :
600
601                  [v,dv,d2v] = d2f (x).
602
603  args : list   : f and d2f's arguments. By default, minimize the 1st
604      or matrix : argument.
605
606  ctl  : vector : Control arguments (see below)
607       or struct
608
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
615
616  CONTROL VARIABLE ctl : (optional). May be a struct or a vector of length
617  ---------------------- 5 or less where NaNs are ignored. Default values
618                         are written <value>.
619  FIELD  VECTOR
620  NAME    POS
621
622  ftol, f N/A    : Stop search when value doesn't improve, as tested by
623
624                    f > Deltaf/max(|f(x)|,1)
625
626              where Deltaf is the decrease in f observed in the last
627              iteration.                                     <10*sqrt(eps)>
628
629  utol, u N/A    : Stop search when updates are small, as tested by
630
631                    u > max { dx(i)/max(|x(i)|,1) | i in 1..N }
632
633              where  dx is the change in the x that occured in the last
634              iteration.                                              <NaN>
635
636  dtol, d N/A    : Stop search when derivative is small, as tested by
637  
638                    d > norm (dv)                                     <eps>
639
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>
642
643  tol, t  ctl(2) : Threshold in termination test chosen by 'crit'  <10*eps>
644
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.
650
651  verbose, v N/A : Be more or less verbose (quiet=0)                    <0>
652
653
654
655 # name: <cell-element>
656 # type: sq_string
657 # elements: 1
658 # length: 75
659  [x,v,nev,h,args] = d2_min(f,d2f,args,ctl,code) - Newton-like minimization
660
661
662
663
664 # name: <cell-element>
665 # type: sq_string
666 # elements: 1
667 # length: 4
668 dcdp
669
670
671 # name: <cell-element>
672 # type: sq_string
673 # elements: 1
674 # length: 268
675  function prt = dcdp (f, p, dp, func[, bounds])
676
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
679  dfdp.m.
680
681  dfpdp is more general and is meant to be used instead of dcdp in
682  optimization.
683
684
685
686 # name: <cell-element>
687 # type: sq_string
688 # elements: 1
689 # length: 48
690  function prt = dcdp (f, p, dp, func[, bounds])
691
692
693
694
695 # name: <cell-element>
696 # type: sq_string
697 # elements: 1
698 # length: 6
699 de_min
700
701
702 # name: <cell-element>
703 # type: sq_string
704 # elements: 1
705 # length: 4432
706  de_min: global optimisation using differential evolution
707
708  Usage: [x, obj_value, nfeval, convergence] = de_min(fcn, control)
709
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
714
715
716  Arguments:  
717  ---------------
718  fcn        string : Name of function. Must return a real value
719  control    vector : (Optional) Control variables, described below
720          or struct
721
722  Returned values:
723  ----------------
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)
732
733  Control variable:   (optional) may be named arguments (i.e. "name",value
734  ----------------    pairs), a struct, or a vector, where
735                      NaN's are ignored.
736
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
753                 if refresh is < 1
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)
762
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
774        larger CR than /exp.
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.
779
780  Default Values:
781  ---------------
782  XVmin = [-2];
783  XVmax = [ 2];
784  constr= 0;
785  const = [];
786  NP    = 10 *D
787  F     = 0.8;
788  CR    = 0.9;
789  strategy = 12;
790  refresh  = 0;
791  VTR   = -Inf;
792  tol   = 1.e-3;
793  maxnfe  = 1e6;
794  maxiter = 1000;
795
796
797  Example to find the minimum of the Rosenbrock saddle:
798  ----------------------------------------------------
799  Define f as:
800                     function result = f(x);
801                       result = 100 * (x(2) - x(1)^2)^2 + (1 - x(1))^2;
802                     end
803  Then type:
804
805         ctl.XVmin = [-2 -2];
806         ctl.XVmax = [ 2  2];
807         [x, obj_value, nfeval, convergence] = de_min (@f, ctl);
808
809  Keywords: global-optimisation optimisation minimisation
810
811
812
813 # name: <cell-element>
814 # type: sq_string
815 # elements: 1
816 # length: 58
817  de_min: global optimisation using differential evolution
818
819
820
821
822 # name: <cell-element>
823 # type: sq_string
824 # elements: 1
825 # length: 5
826 deriv
827
828
829 # name: <cell-element>
830 # type: sq_string
831 # elements: 1
832 # length: 820
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.
838
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.
842
843      H defines the step taken for the derivative calculation. Defaults
844      to 1e-7.
845
846      O defines the order of the calculation. Supported values are 2
847      (h^2 order) or 4 (h^4 order). Defaults to 2.
848
849      N defines the derivative order. Defaults to the 1st derivative of
850      the function. Can be up to the 4th derivative.
851
852      Reference: Numerical Methods for Mathematics, Science, and
853      Engineering by John H. Mathews.
854
855
856
857
858 # name: <cell-element>
859 # type: sq_string
860 # elements: 1
861 # length: 33
862 Calculate derivate of function F.
863
864
865
866 # name: <cell-element>
867 # type: sq_string
868 # elements: 1
869 # length: 4
870 dfdp
871
872
873 # name: <cell-element>
874 # type: sq_string
875 # elements: 1
876 # length: 1050
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 ================================
895
896  dfxpdp is more general and is meant to be used instead of dfdp in
897  optimization.
898
899
900
901 # name: <cell-element>
902 # type: sq_string
903 # elements: 1
904 # length: 80
905  function prt = dfdp (x, f, p, dp, func[, bounds])
906  numerical partial derivative
907
908
909
910 # name: <cell-element>
911 # type: sq_string
912 # elements: 1
913 # length: 5
914 dfpdp
915
916
917 # name: <cell-element>
918 # type: sq_string
919 # elements: 1
920 # length: 1099
921  function jac = dfpdp (p, func[, hook])
922
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:
926
927  hook.f: value of func(p) for p as given in the arguments
928
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)).
932
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)).
937
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
941  are 'fixed').
942
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.
946
947
948
949 # name: <cell-element>
950 # type: sq_string
951 # elements: 1
952 # length: 40
953  function jac = dfpdp (p, func[, hook])
954
955
956
957
958 # name: <cell-element>
959 # type: sq_string
960 # elements: 1
961 # length: 6
962 dfxpdp
963
964
965 # name: <cell-element>
966 # type: sq_string
967 # elements: 1
968 # length: 1115
969  function jac = dfxpdp (x, p, func[, hook])
970
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:
974
975  hook.f: value of func(p, x) for p and x as given in the arguments
976
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));
980
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)).
985
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
989  are 'fixed').
990
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.
994
995
996
997 # name: <cell-element>
998 # type: sq_string
999 # elements: 1
1000 # length: 44
1001  function jac = dfxpdp (x, p, func[, hook])
1002
1003
1004
1005
1006 # name: <cell-element>
1007 # type: sq_string
1008 # elements: 1
1009 # length: 6
1010 expfit
1011
1012
1013 # name: <cell-element>
1014 # type: sq_string
1015 # elements: 1
1016 # length: 1621
1017  USAGE  [alpha,c,rms] = expfit( deg, x1, h, y )
1018
1019  Prony's method for non-linear exponential fitting
1020
1021  Fit function:   \sum_1^{deg} c(i)*exp(alpha(i)*x)
1022
1023  Elements of data vector y must correspond to
1024  equidistant x-values starting at x1 with stepsize h
1025
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.
1030
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.
1035
1036  Example
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)
1041
1042   alpha =
1043     2.0000
1044     1.3000
1045   c =
1046     -0.50000
1047      2.00000
1048   rms = 0.00028461
1049
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.
1054
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.
1059
1060  Demo for a complex fit-function:
1061  deg= 2; N= 20; x1= -(1+i), x= linspace(x1,1+i/2,N).';
1062  h = x(2) - x(1)
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 )
1066
1067
1068
1069 # name: <cell-element>
1070 # type: sq_string
1071 # elements: 1
1072 # length: 48
1073  USAGE  [alpha,c,rms] = expfit( deg, x1, h, y )
1074
1075
1076
1077
1078 # name: <cell-element>
1079 # type: sq_string
1080 # elements: 1
1081 # length: 4
1082 fmin
1083
1084
1085 # name: <cell-element>
1086 # type: sq_string
1087 # elements: 1
1088 # length: 19
1089  alias for fminbnd
1090
1091
1092
1093 # name: <cell-element>
1094 # type: sq_string
1095 # elements: 1
1096 # length: 19
1097  alias for fminbnd
1098
1099
1100
1101
1102 # name: <cell-element>
1103 # type: sq_string
1104 # elements: 1
1105 # length: 5
1106 fmins
1107
1108
1109 # name: <cell-element>
1110 # type: sq_string
1111 # elements: 1
1112 # length: 1412
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
1116
1117      Example usage:   fmins(inline('(x(1)-5).^2+(x(2)-8).^4'),[0;0])
1118
1119      *Inputs*
1120     F
1121           A string containing the name of the function to minimize
1122
1123     X0
1124           A vector of initial parameters fo the function F.
1125
1126     OPTIONS
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
1133           options(5)
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
1140
1141     GRAD
1142           Unused (For compatibility with Matlab)
1143
1144     P1,P2, ...
1145           Optional parameters for function F
1146
1147
1148
1149
1150
1151 # name: <cell-element>
1152 # type: sq_string
1153 # elements: 1
1154 # length: 51
1155 Find the minimum of a funtion of several variables.
1156
1157
1158
1159 # name: <cell-element>
1160 # type: sq_string
1161 # elements: 1
1162 # length: 10
1163 fminsearch
1164
1165
1166 # name: <cell-element>
1167 # type: sq_string
1168 # elements: 1
1169 # length: 300
1170  -- Function File: X = fminsearch (FUN, X0)
1171  -- Function File: [X, FVAL] = fminsearch (FUN, X0, OPTIONS, GRAD, P1,
1172           P2, ...)
1173      Find the minimum of a funtion of several variables.  By default
1174      the method used is the Nelder&Mead Simplex algorithm.
1175
1176      See also: fmin, fmins, nmsmax
1177
1178
1179
1180
1181
1182 # name: <cell-element>
1183 # type: sq_string
1184 # elements: 1
1185 # length: 51
1186 Find the minimum of a funtion of several variables.
1187
1188
1189
1190 # name: <cell-element>
1191 # type: sq_string
1192 # elements: 1
1193 # length: 14
1194 fminunc_compat
1195
1196
1197 # name: <cell-element>
1198 # type: sq_string
1199 # elements: 1
1200 # length: 1594
1201  [x,v,flag,out,df,d2f] = fminunc_compat (f,x,opt,...) - M*tlab-like optimization
1202
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.
1207
1208  Supported options
1209  -----------------
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
1220
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().
1225
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().
1231
1232  This function is a front-end to minimize().
1233
1234
1235
1236 # name: <cell-element>
1237 # type: sq_string
1238 # elements: 1
1239 # length: 50
1240  [x,v,flag,out,df,d2f] = fminunc_compat (f,x,opt,.
1241
1242
1243
1244 # name: <cell-element>
1245 # type: sq_string
1246 # elements: 1
1247 # length: 3
1248 gjp
1249
1250
1251 # name: <cell-element>
1252 # type: sq_string
1253 # elements: 1
1254 # length: 829
1255  m = gjp (m, k[, l])
1256
1257  m: matrix; k, l: row- and column-index of pivot, l defaults to k.
1258
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),
1270  166--168.
1271
1272
1273
1274 # name: <cell-element>
1275 # type: sq_string
1276 # elements: 1
1277 # length: 21
1278  m = gjp (m, k[, l])
1279
1280
1281
1282
1283 # name: <cell-element>
1284 # type: sq_string
1285 # elements: 1
1286 # length: 6
1287 jacobs
1288
1289
1290 # name: <cell-element>
1291 # type: sq_string
1292 # elements: 1
1293 # length: 1030
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.
1297
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
1302      analytic functions.
1303
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.
1308
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.
1312
1313      For example:
1314
1315           f = @(x) [x(1)^2 + x(2); x(2)*exp(x(1))];
1316           Df = jacobs ([1, 2], f)
1317
1318
1319
1320
1321 # name: <cell-element>
1322 # type: sq_string
1323 # elements: 1
1324 # length: 67
1325 Calculate the jacobian of a function using the complex step method.
1326
1327
1328
1329 # name: <cell-element>
1330 # type: sq_string
1331 # elements: 1
1332 # length: 6
1333 leasqr
1334
1335
1336 # name: <cell-element>
1337 # type: sq_string
1338 # elements: 1
1339 # length: 7546
1340 function [f,p,cvg,iter,corp,covp,covr,stdresid,Z,r2]=
1341                    leasqr(x,y,pin,F,{stol,niter,wt,dp,dFdp,options})
1342
1343  Levenberg-Marquardt nonlinear regression of f(x,p) to y(x).
1344
1345  Version 3.beta
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
1362    as described above.
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 =
1390    Inf*ones().
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.
1453
1454           OUTPUT VARIABLES
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.
1465
1466  Not suitable for non-real residuals.
1467
1468  References:
1469  Bard, Nonlinear Parameter Estimation, Academic Press, 1974.
1470  Draper and Smith, Applied Regression Analysis, John Wiley and Sons, 1981.
1471
1472
1473
1474 # name: <cell-element>
1475 # type: sq_string
1476 # elements: 1
1477 # length: 80
1478 function [f,p,cvg,iter,corp,covp,covr,stdresid,Z,r2]=
1479                    leasqr(
1480
1481
1482
1483 # name: <cell-element>
1484 # type: sq_string
1485 # elements: 1
1486 # length: 8
1487 line_min
1488
1489
1490 # name: <cell-element>
1491 # type: sq_string
1492 # elements: 1
1493 # length: 806
1494  [a,fx,nev] = line_min (f, dx, args, narg, h, nev_max) - Minimize f() along dx
1495
1496  INPUT ----------
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
1504
1505  OUTPUT ---------
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
1509
1510  (*) The notation f(x+a*dx) assumes that args == {x}.
1511
1512  Reference: David G Luenberger's Linear and Nonlinear Programming
1513
1514
1515
1516 # name: <cell-element>
1517 # type: sq_string
1518 # elements: 1
1519 # length: 79
1520  [a,fx,nev] = line_min (f, dx, args, narg, h, nev_max) - Minimize f() along dx
1521
1522
1523
1524
1525 # name: <cell-element>
1526 # type: sq_string
1527 # elements: 1
1528 # length: 7
1529 linprog
1530
1531
1532 # name: <cell-element>
1533 # type: sq_string
1534 # elements: 1
1535 # length: 591
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.
1541
1542      Finds
1543
1544           min (f' * x)
1545
1546      (both f and x are column vectors) subject to
1547
1548           A   * x <= b
1549           Aeq * x  = beq
1550           lb <= x <= ub
1551
1552      If not specified, AEQ and BEQ default to empty matrices.
1553
1554      If not specified, the lower bound LB defaults to minus infinite
1555      and the upper bound UB defaults to infinite.
1556
1557      See also: glpk
1558
1559
1560
1561
1562
1563 # name: <cell-element>
1564 # type: sq_string
1565 # elements: 1
1566 # length: 23
1567 Solve a linear problem.
1568
1569
1570
1571 # name: <cell-element>
1572 # type: sq_string
1573 # elements: 1
1574 # length: 6
1575 mdsmax
1576
1577
1578 # name: <cell-element>
1579 # type: sq_string
1580 # elements: 1
1581 # length: 2358
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.
1586         Output arguments:
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)
1592                  (default 1e-3),
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,...).
1607
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).
1611
1612  References:
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.
1622
1623
1624
1625 # name: <cell-element>
1626 # type: sq_string
1627 # elements: 1
1628 # length: 70
1629 MDSMAX  Multidirectional search method for direct search optimization.
1630
1631
1632
1633 # name: <cell-element>
1634 # type: sq_string
1635 # elements: 1
1636 # length: 8
1637 minimize
1638
1639
1640 # name: <cell-element>
1641 # type: sq_string
1642 # elements: 1
1643 # length: 4103
1644  [x,v,nev,...] = minimize (f,args,...) - Minimize f
1645
1646  ARGUMENTS
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
1650
1651  RETURNED VALUES
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)
1657  
1658
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
1661  struct is ignored).
1662  
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.
1666  
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.
1671
1672  'hess'         : Use [fx,dfx,d2fx] = leval (f, args) to compute 1st and
1673                   2nd derivatives, and use a Newton-like method (d2_min).
1674
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).
1680
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).
1684
1685             NOTE : df, d2f or d2i take the same arguments as f.
1686  
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.
1690  
1691  'ndiff'        : Use a variable metric method (bfgs) using numerical
1692                   differentiation.
1693
1694  OPTIONS : STOPPING CRITERIA  Default is to use 'tol'
1695  ---------------------------
1696  'ftol', ftol   : Stop search when value doesn't improve, as tested by
1697
1698               ftol > Deltaf/max(|f(x)|,1)
1699
1700                  where Deltaf is the decrease in f observed in the last
1701                  iteration.                                 Default=10*eps
1702
1703  'utol', utol   : Stop search when updates are small, as tested by
1704
1705               tol > max { dx(i)/max(|x(i)|,1) | i in 1..N }
1706
1707                  where  dx is the change in the x that occured in the last
1708                  iteration.
1709
1710  'dtol',dtol    : Stop search when derivatives are small, as tested by
1711
1712               dtol > max { df(i)*max(|x(i)|,1)/max(v,1) | i in 1..N }
1713
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.
1717
1718  MISC. OPTIONS
1719  -------------
1720  'maxev', m     : Maximum number of function evaluations             <inf>
1721
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
1726
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 :
1732
1733               [x,v,nev] = feval (backend, args, control)
1734                    
1735                   or, if a 2nd order method is used :
1736
1737               [x,v,nev] = feval (backend, control.d2f, args, control)
1738
1739
1740
1741 # name: <cell-element>
1742 # type: sq_string
1743 # elements: 1
1744 # length: 11
1745  [x,v,nev,.
1746
1747
1748
1749 # name: <cell-element>
1750 # type: sq_string
1751 # elements: 1
1752 # length: 15
1753 nelder_mead_min
1754
1755
1756 # name: <cell-element>
1757 # type: sq_string
1758 # elements: 1
1759 # length: 2747
1760  [x0,v,nev] = nelder_mead_min (f,args,ctl) - Nelder-Mead minimization
1761
1762  Minimize 'f' using the Nelder-Mead algorithm. This function is inspired
1763  from the that found in the book "Numerical Recipes".
1764
1765  ARGUMENTS
1766  ---------
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
1771       or struct
1772
1773  RETURNED VALUES
1774  ---------------
1775  x0  : matrix   : Local minimum of f
1776  v   : real     : Value of f in x0
1777  nev : number   : Number of function evaluations
1778  
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>.
1782   OPT.   VECTOR
1783   NAME    POS
1784  ftol,f  N/A    : Stopping criterion : stop search when values at simplex
1785                   vertices are all alike, as tested by 
1786
1787                    f > (max_i (f_i) - min_i (f_i)) /max(max(|f_i|),1)
1788
1789                   where f_i are the values of f at the vertices.  <10*eps>
1790
1791  rtol,r  N/A    : Stop search when biggest radius of simplex, using
1792                   infinity-norm, is small, as tested by :
1793
1794               ctl(2) > Radius                                     <10*eps>
1795
1796  vtol,v  N/A    : Stop search when volume of simplex is small, tested by
1797             
1798               ctl(2) > Vol
1799
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>
1802
1803  tol, t ctl(2)  : Threshold in termination test chosen by 'crit'  <10*eps>
1804
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>
1809
1810                 { x + e_i | i in 0..N } 
1811  
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
1816
1817                 Set ctl(5) to the distance you expect between the starting
1818                 point and the minimum.
1819
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>
1825
1826  verbose, v     Be more or less verbose (quiet=0)                      <0>
1827
1828
1829
1830 # name: <cell-element>
1831 # type: sq_string
1832 # elements: 1
1833 # length: 70
1834  [x0,v,nev] = nelder_mead_min (f,args,ctl) - Nelder-Mead minimization
1835
1836
1837
1838
1839 # name: <cell-element>
1840 # type: sq_string
1841 # elements: 1
1842 # length: 6
1843 nmsmax
1844
1845
1846 # name: <cell-element>
1847 # type: sq_string
1848 # elements: 1
1849 # length: 1967
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.
1854         Output arguments:
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)
1860                  (default 1e-3),
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,...).
1878  References:
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.
1883
1884
1885
1886 # name: <cell-element>
1887 # type: sq_string
1888 # elements: 1
1889 # length: 66
1890 NMSMAX  Nelder-Mead simplex method for direct search optimization.
1891
1892
1893
1894 # name: <cell-element>
1895 # type: sq_string
1896 # elements: 1
1897 # length: 15
1898 nonlin_curvefit
1899
1900
1901 # name: <cell-element>
1902 # type: sq_string
1903 # elements: 1
1904 # length: 1010
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,
1907           SETTINGS)
1908      Frontend for nonlinear fitting of values, computed by a model
1909      function, to observed values.
1910
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
1921      will be set to X.
1922
1923      See also: nonlin_residmin
1924
1925
1926
1927
1928
1929 # name: <cell-element>
1930 # type: sq_string
1931 # elements: 1
1932 # length: 80
1933 Frontend for nonlinear fitting of values, computed by a model function,
1934 to obser
1935
1936
1937
1938 # name: <cell-element>
1939 # type: sq_string
1940 # elements: 1
1941 # length: 10
1942 nonlin_min
1943
1944
1945 # name: <cell-element>
1946 # type: sq_string
1947 # elements: 1
1948 # length: 15376
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.
1955
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.
1959
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.
1964
1965      PIN: real column vector of initial parameters.
1966
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
1971      case.
1972
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).
1985
1986      SETTINGS:
1987
1988      `Algorithm': String specifying the backend. Currently available
1989      are `"lm_feasible"' (default) and `"siman"'. They are described in
1990      separate sections below.
1991
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.
2012
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
2016      argument.
2017
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.
2023
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))'.
2027
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.
2037
2038      `cstep': scalar step size for complex step derivative
2039      approximation. Default: 1e-20.
2040
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.
2045
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.
2051
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
2084      are also specified.
2085
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'.
2089
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.
2093
2094      `TolFun': Minimum fractional improvement in objective function in
2095      an iteration (abortion criterium). Default: .0001.
2096
2097      `MaxIter': Maximum number of iterations (abortion criterium).
2098      Default: backend-specific.
2099
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.
2103
2104      `max_fract_change': Column Vector, enforced maximum fractional
2105      change in parameters in an iteration. Default: backend-specific.
2106
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.
2111
2112      `debug': Logical scalar, default: `false'. Will be passed to the
2113      backend, which might print debugging information if true.
2114
2115      Structure-based parameter handling
2116
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.
2123
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.
2127
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
2144      example.
2145
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.
2150
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.
2159
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.
2171
2172      Description of backends
2173
2174      "lm_feasible"
2175
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
2185      algorithm.
2186
2187      Returned value CVG will be `2' or `3' for success and `0' or `-4'
2188      for failure (see above for meaning).
2189
2190      Backend-specific defaults are: `MaxIter': 20, `fract_prec': `zeros
2191      (size (parameters))', `max_fract_change': `Inf' for all parameters.
2192
2193      Interpretation of `Display': if set to `"iter"', currently only
2194      information on applying `max_fract_change' is printed.
2195
2196      "siman"
2197
2198      A simulated annealing (stochastic) optimizer, changing all
2199      parameters at once in a single step, so being suitable for
2200      non-bound constraints.
2201
2202      No gradient or hessian of the objective function is used. The
2203      settings `MaxIter', `fract_prec', `TolFun', and `max_fract_change'
2204      are not honoured.
2205
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
2222      rejected.
2223
2224      Steps with increase `diff' of objective function are accepted if
2225      `rand (1) < exp (- diff / T)', where `T' is the temperature of the
2226      current iteration.
2227
2228      If regaining of constraints failed, optimization will be aborted
2229      and returned value of CVG will be `0'. Otherwise, CVG will be `1'.
2230
2231      Interpretation of `Display': if set to `"iter"', an informational
2232      line is printed after each iteration.
2233
2234
2235
2236
2237
2238 # name: <cell-element>
2239 # type: sq_string
2240 # elements: 1
2241 # length: 79
2242 Frontend for constrained nonlinear minimization of a scalar objective
2243 function.
2244
2245
2246
2247 # name: <cell-element>
2248 # type: sq_string
2249 # elements: 1
2250 # length: 15
2251 nonlin_residmin
2252
2253
2254 # name: <cell-element>
2255 # type: sq_string
2256 # elements: 1
2257 # length: 12987
2258  -- Function File: [P, RESID, CVG, OUTP] = nonlin_residmin (F, PIN)
2259  -- Function File: [P, RESID, CVG, OUTP] = nonlin_residmin (F, PIN,
2260           SETTINGS)
2261      Frontend for nonlinear minimization of residuals returned by a
2262      model function.
2263
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.
2267
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.
2271
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.
2276
2277      PIN: real column vector of initial parameters.
2278
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
2283      case.
2284
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).
2295
2296      SETTINGS:
2297
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
2301      section below.
2302
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.
2323
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.
2329
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))'.
2333
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
2342      must be specified.
2343
2344      `cstep': scalar step size for complex step derivative
2345      approximation. Default: 1e-20.
2346
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.
2351
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.
2357
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
2390      are also specified.
2391
2392      `equc': Equality constraints. Specified the same way as inequality
2393      constraints (see `inequc').
2394
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.
2398
2399      `weights': Array of weights for the residuals. Dimensions must
2400      match.
2401
2402      `TolFun': Minimum fractional improvement in objective function
2403      (e.g. sum of squares) in an iteration (abortion criterium).
2404      Default: .0001.
2405
2406      `MaxIter': Maximum number of iterations (abortion criterium).
2407      Default: backend-specific.
2408
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.
2412
2413      `max_fract_change': Column Vector, enforced maximum fractional
2414      change in parameters in an iteration. Default: backend-specific.
2415
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.
2420
2421      `plot_cmd': Function enabling backend to plot results or
2422      intermediate results. Will be called with current computed
2423      residuals. Default: plot nothing.
2424
2425      `debug': Logical scalar, default: `false'. Will be passed to the
2426      backend, which might print debugging information if true.
2427
2428      Structure-based parameter handling
2429
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.
2436
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.
2440
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.
2452
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.
2457
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.
2466
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.
2478
2479      Description of backends (currently only one)
2480
2481      "lm_svd_feasible"
2482
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
2486      optimization.
2487
2488      Parameters with identical lower and upper bounds will be fixed.
2489
2490      Returned value CVG will be `0', `2', or `3'.
2491
2492      Backend-specific defaults are: `MaxIter': 20, `fract_prec': `zeros
2493      (size (parameters))', `max_fract_change': `Inf' for all parameters.
2494
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.
2498
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.
2503
2504      See also: nonlin_curvefit
2505
2506
2507
2508
2509
2510 # name: <cell-element>
2511 # type: sq_string
2512 # elements: 1
2513 # length: 78
2514 Frontend for nonlinear minimization of residuals returned by a model
2515 function.
2516
2517
2518
2519 # name: <cell-element>
2520 # type: sq_string
2521 # elements: 1
2522 # length: 3
2523 nrm
2524
2525
2526 # name: <cell-element>
2527 # type: sq_string
2528 # elements: 1
2529 # length: 153
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.
2533
2534
2535
2536
2537 # name: <cell-element>
2538 # type: sq_string
2539 # elements: 1
2540 # length: 69
2541 Using X0 as a starting point find a minimum of the scalar function F.
2542
2543
2544
2545 # name: <cell-element>
2546 # type: sq_string
2547 # elements: 1
2548 # length: 14
2549 optim_problems
2550
2551
2552 # name: <cell-element>
2553 # type: sq_string
2554 # elements: 1
2555 # length: 64
2556  Problems for testing optimizers. Documentation is in the code.
2557
2558
2559
2560 # name: <cell-element>
2561 # type: sq_string
2562 # elements: 1
2563 # length: 33
2564  Problems for testing optimizers.
2565
2566
2567
2568 # name: <cell-element>
2569 # type: sq_string
2570 # elements: 1
2571 # length: 15
2572 optimset_compat
2573
2574
2575 # name: <cell-element>
2576 # type: sq_string
2577 # elements: 1
2578 # length: 1208
2579  opt = optimset_compat (...)         - manipulate m*tlab-style options structure
2580  
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.
2585
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()).
2589  
2590  Two extra options are supported which indicate how to use directly octave
2591  optimization tools (such as minimize() and other backends):
2592
2593  "MinEquiv", [on|off] : Tell 'fminunc()' not to minimize 'fun', but
2594                         instead return the option passed to minimize().
2595
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().
2601
2602
2603
2604 # name: <cell-element>
2605 # type: sq_string
2606 # elements: 1
2607 # length: 25
2608  opt = optimset_compat (.
2609
2610
2611
2612 # name: <cell-element>
2613 # type: sq_string
2614 # elements: 1
2615 # length: 9
2616 poly_2_ex
2617
2618
2619 # name: <cell-element>
2620 # type: sq_string
2621 # elements: 1
2622 # length: 353
2623   ex = poly_2_ex (l, f)       - Extremum of a 1-var deg-2 polynomial
2624
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).
2627  
2628  ex : 1 : Value for which f reaches its extremum
2629  
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.
2632
2633
2634
2635
2636 # name: <cell-element>
2637 # type: sq_string
2638 # elements: 1
2639 # length: 69
2640   ex = poly_2_ex (l, f)       - Extremum of a 1-var deg-2 polynomial
2641
2642
2643
2644
2645 # name: <cell-element>
2646 # type: sq_string
2647 # elements: 1
2648 # length: 8
2649 polyconf
2650
2651
2652 # name: <cell-element>
2653 # type: sq_string
2654 # elements: 1
2655 # length: 1439
2656  [y,dy] = polyconf(p,x,s)
2657
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.
2661
2662  polyconf(...,['ci'|'pi'])
2663
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.
2674
2675  polyconf(...,1-alpha)
2676
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.
2681
2682  Example:
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-;;');
2688
2689
2690
2691 # name: <cell-element>
2692 # type: sq_string
2693 # elements: 1
2694 # length: 26
2695  [y,dy] = polyconf(p,x,s)
2696
2697
2698
2699
2700 # name: <cell-element>
2701 # type: sq_string
2702 # elements: 1
2703 # length: 10
2704 polyfitinf
2705
2706
2707 # name: <cell-element>
2708 # type: sq_string
2709 # elements: 1
2710 # length: 4746
2711  function [A,REF,HMAX,H,R,EQUAL] = polyfitinf(M,N,K,X,Y,EPSH,MAXIT,REF0)
2712
2713    Best polynomial approximation in discrete uniform norm
2714
2715    INPUT VARIABLES:
2716
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
2733
2734    OUTPUT VARIABLES:
2735
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
2744                    EQUAL=1 :  succesful
2745                    EQUAL=0 :  convergence not acheived
2746                    EQUAL=-1:  input error
2747                    EQUAL=-2:  algorithm failure
2748
2749    Relies on function EXCH, provided below.
2750
2751    Example: 
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
2757
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: 
2763
2764  N = 10000; K = 2; EPSH = 10^-12; MAXIT = 10;  X = linspace(0,1,N);
2765  for i = 1:2
2766      k = 2*i-1; Y = abs(X).^k;
2767      for j = 1:4
2768          M = 2^j;
2769          [~,~,HMAX] = polyfitinf(M,N,K,X,Y,EPSH,MAXIT);
2770          approxerror(i,j) = HMAX;
2771      end
2772  end
2773  disp('Table 3.1 from Approximation theory and methods, M.J.D.POWELL, p. 27');
2774  disp(' ');
2775  disp('            n          K=1          K=3'); 
2776  disp(' '); format short g;
2777  disp([(2.^(1:4))' approxerror']);
2778
2779    ALGORITHM:
2780
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,
2786       n
2787     (   ) iterations ( usually p = M + 2. See also function EXCH ).
2788       p
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
2792    data chosen.
2793
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.
2797
2798    Team Fall 98 (Revision 1.0):
2799            Chanchai Aniwathananon
2800            Crhistopher Mehl
2801            David A. Duran
2802            Saulo P. Oliveira
2803
2804    Team Spring 11 (Revision 1.1): Manuchehr Aminian
2805
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)
2810
2811    NOTES:
2812
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.
2819
2820    Some parts of the code can still take advantage of vectorization.  
2821
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
2827
2828
2829
2830 # name: <cell-element>
2831 # type: sq_string
2832 # elements: 1
2833 # length: 73
2834  function [A,REF,HMAX,H,R,EQUAL] = polyfitinf(M,N,K,X,Y,EPSH,MAXIT,REF0)
2835
2836
2837
2838
2839 # name: <cell-element>
2840 # type: sq_string
2841 # elements: 1
2842 # length: 6
2843 powell
2844
2845
2846 # name: <cell-element>
2847 # type: sq_string
2848 # elements: 1
2849 # length: 2656
2850  -- Function File: [ P, OBJ_VALUE, CONVERGENCE, ITERS, NEVS] = powell
2851           (F, P0, CONTROL)
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]
2855
2856 Arguments
2857 ---------
2858
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)
2862
2863         * P0: An initial value of the function argument to minimize
2864
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))
2869
2870              -  TolFun: minimum amount by which function value must
2871                decrease in each iteration to continue (default is 1E-8)
2872
2873              -  MaxFunEvals: maximum function evaluations  (positive
2874                integer, or -1 or Inf for unlimited (default))
2875
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)
2882
2883 Examples
2884 --------
2885
2886           y = @(x, s) x(1) ^ 2 + x(2) ^ 2 + s;
2887           o = optimset('MaxIter', 100, 'TolFun', 1E-10);
2888           s = 1;
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
2891
2892
2893 Returns:
2894 --------
2895
2896         * P: the minimizing value of the function argument
2897
2898         * OBJ_VALUE: the value of F() at P
2899
2900         * CONVERGENCE: 1 if normal convergence, 0 if not
2901
2902         * ITERS: number of iterations performed
2903
2904         * NEVS: number of function evaluations
2905
2906 References
2907 ----------
2908
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
2912
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
2916           (Section 10.5)
2917
2918
2919
2920
2921 # name: <cell-element>
2922 # type: sq_string
2923 # elements: 1
2924 # length: 53
2925 Multidimensional minimization (direction-set method).
2926
2927
2928
2929 # name: <cell-element>
2930 # type: sq_string
2931 # elements: 1
2932 # length: 13
2933 residmin_stat
2934
2935
2936 # name: <cell-element>
2937 # type: sq_string
2938 # elements: 1
2939 # length: 3188
2940  -- Function File: INFO = residmin_stat (F, P, SETTINGS)
2941      Frontend for computation of statistics for a residual-based
2942      minimization.
2943
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.
2956
2957      Currently, the following statistics (or general information) can be
2958      requested:
2959
2960      `dfdp': Jacobian of model function with respect to parameters.
2961
2962      `covd': Covariance matrix of data (typically guessed by applying a
2963      factor to the covariance matrix of the residuals).
2964
2965      `covp': Covariance matrix of final parameters.
2966
2967      `corp': Correlation matrix of final parameters.
2968
2969      Further SETTINGS
2970
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
2980      `diff_onesided'.
2981
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
2990      least squares).
2991
2992      Backend-specific information
2993
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'.
3003
3004      See also: curvefit_stat
3005
3006
3007
3008
3009
3010 # name: <cell-element>
3011 # type: sq_string
3012 # elements: 1
3013 # length: 73
3014 Frontend for computation of statistics for a residual-based
3015 minimization.
3016
3017
3018
3019 # name: <cell-element>
3020 # type: sq_string
3021 # elements: 1
3022 # length: 10
3023 rosenbrock
3024
3025
3026 # name: <cell-element>
3027 # type: sq_string
3028 # elements: 1
3029 # length: 195
3030  Rosenbrock function - used to create example obj. fns.
3031
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.
3035
3036
3037
3038 # name: <cell-element>
3039 # type: sq_string
3040 # elements: 1
3041 # length: 50
3042  Rosenbrock function - used to create example obj.
3043
3044
3045
3046 # name: <cell-element>
3047 # type: sq_string
3048 # elements: 1
3049 # length: 13
3050 samin_example
3051
3052
3053 # name: <cell-element>
3054 # type: sq_string
3055 # elements: 1
3056 # length: 16
3057  dimensionality
3058
3059
3060
3061 # name: <cell-element>
3062 # type: sq_string
3063 # elements: 1
3064 # length: 16
3065  dimensionality
3066
3067
3068
3069
3070 # name: <cell-element>
3071 # type: sq_string
3072 # elements: 1
3073 # length: 14
3074 test_fminunc_1
3075
3076
3077 # name: <cell-element>
3078 # type: sq_string
3079 # elements: 1
3080 # length: 92
3081  Plain run, just to make sure ######################################
3082  Minimum wrt 'x' is y0
3083
3084
3085
3086 # name: <cell-element>
3087 # type: sq_string
3088 # elements: 1
3089 # length: 80
3090  Plain run, just to make sure ######################################
3091  Minimum wr
3092
3093
3094
3095 # name: <cell-element>
3096 # type: sq_string
3097 # elements: 1
3098 # length: 10
3099 test_min_1
3100
3101
3102 # name: <cell-element>
3103 # type: sq_string
3104 # elements: 1
3105 # length: 63
3106  [x,v,niter] = feval (optim_func, "testfunc","dtestf", xinit);
3107
3108
3109
3110 # name: <cell-element>
3111 # type: sq_string
3112 # elements: 1
3113 # length: 63
3114  [x,v,niter] = feval (optim_func, "testfunc","dtestf", xinit);
3115
3116
3117
3118
3119 # name: <cell-element>
3120 # type: sq_string
3121 # elements: 1
3122 # length: 10
3123 test_min_2
3124
3125
3126 # name: <cell-element>
3127 # type: sq_string
3128 # elements: 1
3129 # length: 60
3130  [xlev,vlev,nlev] = feval(optim_func, "ff", "dff", xinit) ;
3131
3132
3133
3134 # name: <cell-element>
3135 # type: sq_string
3136 # elements: 1
3137 # length: 60
3138  [xlev,vlev,nlev] = feval(optim_func, "ff", "dff", xinit) ;
3139
3140
3141
3142
3143 # name: <cell-element>
3144 # type: sq_string
3145 # elements: 1
3146 # length: 10
3147 test_min_3
3148
3149
3150 # name: <cell-element>
3151 # type: sq_string
3152 # elements: 1
3153 # length: 166
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));
3157
3158
3159
3160 # name: <cell-element>
3161 # type: sq_string
3162 # elements: 1
3163 # length: 80
3164  [xlev,vlev,nlev] = feval (optim_func, "ff", "dff", xinit, "extra", extra) ;
3165  [x
3166
3167
3168
3169 # name: <cell-element>
3170 # type: sq_string
3171 # elements: 1
3172 # length: 10
3173 test_min_4
3174
3175
3176 # name: <cell-element>
3177 # type: sq_string
3178 # elements: 1
3179 # length: 173
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});
3183  ctl.df = "dff";
3184
3185
3186
3187 # name: <cell-element>
3188 # type: sq_string
3189 # elements: 1
3190 # length: 80
3191  Plain run, just to make sure ######################################
3192  Minimum wr
3193
3194
3195
3196 # name: <cell-element>
3197 # type: sq_string
3198 # elements: 1
3199 # length: 15
3200 test_minimize_1
3201
3202
3203 # name: <cell-element>
3204 # type: sq_string
3205 # elements: 1
3206 # length: 92
3207  Plain run, just to make sure ######################################
3208  Minimum wrt 'x' is y0
3209
3210
3211
3212 # name: <cell-element>
3213 # type: sq_string
3214 # elements: 1
3215 # length: 80
3216  Plain run, just to make sure ######################################
3217  Minimum wr
3218
3219
3220
3221 # name: <cell-element>
3222 # type: sq_string
3223 # elements: 1
3224 # length: 22
3225 test_nelder_mead_min_1
3226
3227
3228 # name: <cell-element>
3229 # type: sq_string
3230 # elements: 1
3231 # length: 29
3232  Use vanilla nelder_mead_min
3233
3234
3235
3236 # name: <cell-element>
3237 # type: sq_string
3238 # elements: 1
3239 # length: 29
3240  Use vanilla nelder_mead_min
3241
3242
3243
3244
3245 # name: <cell-element>
3246 # type: sq_string
3247 # elements: 1
3248 # length: 22
3249 test_nelder_mead_min_2
3250
3251
3252 # name: <cell-element>
3253 # type: sq_string
3254 # elements: 1
3255 # length: 70
3256
3257  Test using volume #################################################
3258
3259
3260
3261 # name: <cell-element>
3262 # type: sq_string
3263 # elements: 1
3264 # length: 70
3265
3266  Test using volume #################################################
3267
3268
3269
3270
3271 # name: <cell-element>
3272 # type: sq_string
3273 # elements: 1
3274 # length: 13
3275 test_wpolyfit
3276
3277
3278 # name: <cell-element>
3279 # type: sq_string
3280 # elements: 1
3281 # length: 34
3282           x         y          dy
3283
3284
3285
3286 # name: <cell-element>
3287 # type: sq_string
3288 # elements: 1
3289 # length: 34
3290           x         y          dy
3291
3292
3293
3294
3295 # name: <cell-element>
3296 # type: sq_string
3297 # elements: 1
3298 # length: 6
3299 vfzero
3300
3301
3302 # name: <cell-element>
3303 # type: sq_string
3304 # elements: 1
3305 # length: 1964
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).
3312
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
3316      FUN.
3317
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.
3325
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:
3329
3330         * 1 The algorithm converged to a solution.
3331
3332         * 0 Maximum number of iterations or function evaluations has
3333           been reached.
3334
3335         * -1 The algorithm has been terminated from user output
3336           function.
3337
3338         * -5 The algorithm may have converged to a singular point.
3339
3340      OUTPUT is a structure containing runtime information about the
3341      `fzero' algorithm.  Fields in the structure are:
3342
3343         * iterations Number of iterations through loop.
3344
3345         * nfev Number of function evaluations.
3346
3347         * bracketx A two-column matrix with the final bracketing of the
3348           zero along the x-axis.
3349
3350         * brackety A two-column matrix with the final bracketing of the
3351           zero along the y-axis.
3352
3353      See also: optimset, fsolve
3354
3355
3356
3357
3358
3359 # name: <cell-element>
3360 # type: sq_string
3361 # elements: 1
3362 # length: 21
3363 A variant of `fzero'.
3364
3365
3366
3367 # name: <cell-element>
3368 # type: sq_string
3369 # elements: 1
3370 # length: 8
3371 wpolyfit
3372
3373
3374 # name: <cell-element>
3375 # type: sq_string
3376 # elements: 1
3377 # length: 2922
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.
3383
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.
3387
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.
3395
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.
3401
3402      You can use an F test to determine if a higher order polynomial
3403      improves the fit:
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.
3411
3412      You can estimate the uncertainty in the polynomial coefficients
3413      themselves using
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.
3417
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,
3423      (x-mu(1))/mu(2)'
3424
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.
3428
3429      Example
3430           x = linspace(0,4,20);
3431           dy = (1+rand(size(x)))/2;
3432           y = polyval([2,3,1],x) + dy.*randn(size(x));
3433           wpolyfit(x,y,dy,2);
3434
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
3438      caution.
3439
3440      Hocking, RR (2003). Methods and Applications of Linear Models.
3441      New Jersey: John Wiley and Sons, Inc.
3442
3443
3444    See also: polyfit, polyconf
3445
3446
3447
3448
3449 # name: <cell-element>
3450 # type: sq_string
3451 # elements: 1
3452 # length: 80
3453 Return the coefficients of a polynomial P(X) of degree N that minimizes
3454 `sumsq (
3455
3456
3457
3458 # name: <cell-element>
3459 # type: sq_string
3460 # elements: 1
3461 # length: 11
3462 wrap_f_dfdp
3463
3464
3465 # name: <cell-element>
3466 # type: sq_string
3467 # elements: 1
3468 # length: 387
3469  [ret1, ret2] = wrap_f_dfdp (f, dfdp, varargin)
3470
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
3476  functions.
3477
3478
3479
3480 # name: <cell-element>
3481 # type: sq_string
3482 # elements: 1
3483 # length: 48
3484  [ret1, ret2] = wrap_f_dfdp (f, dfdp, varargin)
3485
3486
3487
3488
3489 # name: <cell-element>
3490 # type: sq_string
3491 # elements: 1
3492 # length: 6
3493 wsolve
3494
3495
3496 # name: <cell-element>
3497 # type: sq_string
3498 # elements: 1
3499 # length: 1736
3500  [x,s] = wsolve(A,y,dy)
3501
3502  Solve a potentially over-determined system with uncertainty in
3503  the values. 
3504
3505      A x = y +/- dy
3506
3507  Use QR decomposition for increased accuracy.  Estimate the 
3508  uncertainty for the solution from the scatter in the data.
3509
3510  The returned structure s contains
3511
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
3515
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 
3520
3521     dx = sqrt(N*sumsq(inv(s.R'))')
3522
3523  where N = normr^2/df, or N = 1 if df = 0.
3524
3525  Example 1: weighted system
3526
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'))');
3531     res = [xin, x, dx]
3532
3533  Example 2: weighted overdetermined system  y = x1 + 2*x2 + 3*x3 + e
3534
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);
3539     res = [xin, x, dx]
3540
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.
3550
3551
3552
3553 # name: <cell-element>
3554 # type: sq_string
3555 # elements: 1
3556 # length: 24
3557  [x,s] = wsolve(A,y,dy)
3558
3559
3560
3561
3562
3563