vb_exponential_estimate

 MSE estimation of exponential function parameters [a, b, c]
    [a, b, c, err] = vb_exponential_estimate(y,x)

    y = a + b * exp( c * x )
  
  err = mean((y - a - b*exp( c * x )).^2)/mean(y.^2)

  if mode=1 is given, a = 0 is assumed:
    [a, b, c] = vb_exponential_estimate(y,x, mode)
    y = b * exp( c * x ) 

 Masa-aki Sato 2009-12-22

 Copyright (C) 2011, ATR All Rights Reserved.
 License : New BSD License(see VBMEG_LICENSE.txt)

0001 function    [a, b, c, err, E, A, C] = vb_exponential_estimate(y,x,mode)
0002 % MSE estimation of exponential function parameters [a, b, c]
0003 %    [a, b, c, err] = vb_exponential_estimate(y,x)
0004 %
0005 %    y = a + b * exp( c * x )
0006 %
0007 %  err = mean((y - a - b*exp( c * x )).^2)/mean(y.^2)
0008 %
0009 %  if mode=1 is given, a = 0 is assumed:
0010 %    [a, b, c] = vb_exponential_estimate(y,x, mode)
0011 %    y = b * exp( c * x )
0012 %
0013 % Masa-aki Sato 2009-12-22
0014 %
0015 % Copyright (C) 2011, ATR All Rights Reserved.
0016 % License : New BSD License(see VBMEG_LICENSE.txt)
0017 
0018 %  Niter : # of iteration for MSE optimization = 5000 [Default]
0019 
0020 if nargin < 3, mode = 0; end;
0021 
0022 Niter  = 5000; 
0023 errmin = 1e-10;
0024 ddmin  = 1e-10;
0025 
0026 % --- MSE for ( log(y) = c*x + B )
0027 [a, b, c] = exponent_init(y,x);
0028 ex = exp(c*x);
0029 
0030 yy = mean(y.^2);
0031 
0032 ddmin  = ddmin * yy;
0033 %errmin = errmin* yy;
0034 
0035 A = zeros(Niter,1);
0036 C = zeros(Niter,1);
0037 E = zeros(Niter,1);
0038 
0039 % Quasi-Newton iteration
0040 for n=1:Niter
0041     % error
0042     dy = (a + b*ex - y);
0043     
0044     % E(c)''
0045     da  = b*mean( dy.*ex.*(x.^2) );
0046     db  = (b^2)*mean( (ex.*x).^2 );
0047     ddE = da + db;
0048     
0049     if ddE < ddmin, 
0050         ddE = db; 
0051         fprintf('-')
0052     end;
0053     
0054     % E(c)'
0055     dE = b * mean( dy.*x.*ex );
0056     % dc = - E'/E''
0057     c  = c - dE/ddE;
0058     
0059     ex = exp(c*x);
0060     
0061     % y = a + b*ex
0062     [b,a] = slope_estimate(y,ex,mode);
0063     
0064     A(n) = a;
0065     C(n) = c;
0066     E(n) = mean((y - a - b*ex).^2)/yy;
0067     
0068     if n > 1 && abs(E(n-1) - E(n)) < errmin, 
0069         break; 
0070     end;
0071     
0072 end
0073 
0074 err = sum(abs(y - a - b*ex))/sum(abs(y));
0075 
0076 if nargout < 5, return; end;
0077 
0078 A = A(1:n);
0079 C = C(1:n);
0080 E = E(1:n);
0081 
0082 %fprintf('\n MSE err for y = %g\n',E(end))
0083 
0084 return
0085 
0086 %
0087 % --- MSE for ( log(y) = c*x + B )
0088 %     y = a + b * exp( c * x )
0089 %
0090 function    [a,b,c] = exponent_init(y,x)
0091 
0092 ix = find(y > 0); 
0093 
0094 %  log(y) = c * x + B
0095 xx = x(ix);
0096 ly = log(y(ix));
0097 
0098 c = slope_estimate(ly, xx);
0099 
0100 ex = exp(c*x);
0101 
0102 %  y = a + b * exp( c * x ) = a + b * ex
0103 [b,a] = slope_estimate(y, ex, 1);
0104 
0105 return
0106 %
0107 % --- Slope estimation: MSE for ( y = c*x + b )
0108 %
0109 function    [c, b] = slope_estimate(y,x,mode)
0110 % slope estimation
0111 %   [c, b] = slope_estimate(y,x)
0112 %   [c, b] = slope_estimate(y,x,mode)
0113 % mode = 0: [default]
0114 %   y = c * x + b
0115 %   dy = (y - <y>) = c * (x -<x>) = c * dx
0116 %   y = c * x + (<y> - c * <x>)
0117 %   c = <dy*dx>/<dx.^2>
0118 %   b = <y> - c * <x>
0119 % mode = 1:
0120 %   y = c * x
0121 %   c = <y*x>/<x.^2>
0122 
0123 if nargin==2 || mode~=1
0124     % mean
0125     mx = mean(x);
0126     my = mean(y);
0127     
0128     % zero mean
0129     dx = x - mx;
0130     dy = y - my;
0131     c = sum(dy.*dx)/sum(dx.^2);
0132     b = my - c * mx;
0133 else
0134     c = sum(y.*x)/sum(x.^2);
0135     b = 0;
0136 end
0137 
0138 return