/usr/share/octave/packages/optim-1.3.0/expfit.m is in octave-optim 1.3.0-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 | ## Copyright (C) 2000 Gert Van den Eynde <na.gvandeneynde@na-net.ornl.gov>
## Copyright (C) 2002 Rolf Fabian <fabian@tu-cottbus.de>
##
## This program is free software; you can redistribute it and/or modify it under
## the terms of the GNU General Public License as published by the Free Software
## Foundation; either version 3 of the License, or (at your option) any later
## version.
##
## This program is distributed in the hope that it will be useful, but WITHOUT
## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
## details.
##
## You should have received a copy of the GNU General Public License along with
## this program; if not, see <http://www.gnu.org/licenses/>.
## USAGE [alpha,c,rms] = expfit( deg, x1, h, y )
##
## Prony's method for non-linear exponential fitting
##
## Fit function: \sum_1^{deg} c(i)*exp(alpha(i)*x)
##
## Elements of data vector y must correspond to
## equidistant x-values starting at x1 with stepsize h
##
## The method is fully compatible with complex linear
## coefficients c, complex nonlinear coefficients alpha
## and complex input arguments y, x1, non-zero h .
## Fit-order deg must be a real positive integer.
##
## Returns linear coefficients c, nonlinear coefficients
## alpha and root mean square error rms. This method is
## known to be more stable than 'brute-force' non-linear
## least squares fitting.
##
## Example
## x0 = 0; step = 0.05; xend = 5; x = x0:step:xend;
## y = 2*exp(1.3*x)-0.5*exp(2*x);
## error = (rand(1,length(y))-0.5)*1e-4;
## [alpha,c,rms] = expfit(2,x0,step,y+error)
##
## alpha =
## 2.0000
## 1.3000
## c =
## -0.50000
## 2.00000
## rms = 0.00028461
##
## The fit is very sensitive to the number of data points.
## It doesn't perform very well for small data sets.
## Theoretically, you need at least 2*deg data points, but
## if there are errors on the data, you certainly need more.
##
## Be aware that this is a very (very,very) ill-posed problem.
## By the way, this algorithm relies heavily on computing the
## roots of a polynomial. I used 'roots.m', if there is
## something better please use that code.
##
## Demo for a complex fit-function:
## deg= 2; N= 20; x1= -(1+i), x= linspace(x1,1+i/2,N).';
## h = x(2) - x(1)
## y= (2+i)*exp( (-1-2i)*x ) + (-1+3i)*exp( (2+3i)*x );
## A= 5e-2; y+= A*(randn(N,1)+randn(N,1)*i); % add complex noise
## [alpha,c,rms]= expfit( deg, x1, h, y )
function [a,c,rms] = expfit(deg,x1,h,y)
% Check input
if deg<1, error('expfit: deg must be >= 1'); end
if ~h, error('expfit: vanishing stepsize h'); end
if ( N=length( y=y(:) ) ) < 2*deg
error('expfit: less than %d samples',2*deg);
end
% Solve for polynomial coefficients
A = hankel( y(1:N-deg),y(N-deg:N) );
s = - A(:,1:deg) \ A(:,deg+1);
% Compose polynomial
p = flipud([s;1]);
% Compute roots
u = roots(p);
% nonlinear coefficients
a = log(u)/h;
% Compose second matrix A(i,j) = u(j)^(i-1)
A = ( ones(N,1) * u(:).' ) .^ ( [0:N-1]' * ones(1,deg) );
% Solve linear system
f = A\y;
% linear coefficients
c = f./exp( a*x1 );
% Compute rms of y - approx
% where approx(i) = sum_k ( c(k) * exp(x(i)*a(k)) )
if nargout > 2
x = x1+h*[0:N-1];
approx = exp( x(:) * a(:).' ) * c(:);
rms = sqrt( sumsq(approx - y) );
end
endfunction
%!demo % same as in help - part
%! deg= 2; N= 20; x1= -(1+i), x= linspace(x1,1+i/2,N).';
%! h = x(2) - x(1)
%! y= (2+i)*exp( (-1-2i)*x ) + (-1+3i)*exp( (2+3i)*x );
%! A= 5e-2; y+= A*(randn(N,1)+randn(N,1)*i); % add complex noise
%! [alpha,c,rms]= expfit( deg, x1, h, y )
%!demo % demo for stepsize with negative real part
%! deg= 2; N= 20; x1= +3*(1+i), x= linspace(x1,1+i/3,N).';
%! h = x(2) - x(1)
%! y= (2+i)*exp( (-1-2i)*x ) + (-1+3i)*exp( (2+3i)*x );
%! A= 5e-2; y+= A*(randn(N,1)+randn(N,1)*i); % add complex noise
%! [alpha,c,rms]= expfit( deg, x1, h, y )
%!demo
%! x0 = 1.5; step = 0.05; xend = 5;
%! a = [1.3, 2]';
%! c = [2, -0.5]';
%! v = 1e-4;
%!
%! x = x0:step:xend;
%! y = exp (x(:) * a(:).') * c(:);
%! err = randn (size (y)) * v;
%! plot (x, y + err);
%!
%! [a_out, c_out, rms] = expfit (2, x0, step, y+err)
|