util_FitPowerLaw

util_FitPowerLaw fits a power law (to recession segments).
   dQ/dt = -a*Q^b
   Different options are available (e.g. linear fitting in loglog-space,
   non-linear, or fitting a fixed slope).

   INPUT
   Q_rec: streamflow (recession periods only)
   dQdt_rec: corresponding flow rate gradients
   fitting_type: specifies fitting procedure: linear, non-linear, linear
       with slope 1, linear with slope 2
   weights: weights to fit curve

   OUTPUT
   a: scaling parameter
   b: parameter of non-linearity
   error_flag: 0 (no error), 1 (warning), 2 (errorin data check), 3
       (error in signature calculation)
   error_str: string contraining error description

   EXAMPLE
   % load example data
   data = load('example/example_data/33029_daily.mat');
   Q = data.Q;
   t = data.t;
   flow_section = util_RecessionSegments(Q,t,'plot_results',true); % get recession segments
   [dQdt, Qm, flow_section, R2] = util_dQdt(Q, t, flow_section); % get flow rate gradient
   rec = ~isnan(Qm);
   [a, b] = util_FitPowerLaw(Qm(rec), dQdt(rec));
   [a, b] = util_FitPowerLaw(Qm(rec), dQdt(rec), 'fitting_type', 'linear', 'weights', R2(rec))

   Copyright (C) 2020
   This software is distributed under the GNU Public License Version 3.
   See <https://www.gnu.org/licenses/gpl-3.0.en.html> for details.

0001 function [a, b, error_flag, error_str] = ...
0002     util_FitPowerLaw(Q_rec, dQdt_rec, fitting_type, weights)
0003 %util_FitPowerLaw fits a power law (to recession segments).
0004 %   dQ/dt = -a*Q^b
0005 %   Different options are available (e.g. linear fitting in loglog-space,
0006 %   non-linear, or fitting a fixed slope).
0007 %
0008 %   INPUT
0009 %   Q_rec: streamflow (recession periods only)
0010 %   dQdt_rec: corresponding flow rate gradients
0011 %   fitting_type: specifies fitting procedure: linear, non-linear, linear
0012 %       with slope 1, linear with slope 2
0013 %   weights: weights to fit curve
0014 %
0015 %   OUTPUT
0016 %   a: scaling parameter
0017 %   b: parameter of non-linearity
0018 %   error_flag: 0 (no error), 1 (warning), 2 (errorin data check), 3
0019 %       (error in signature calculation)
0020 %   error_str: string contraining error description
0021 %
0022 %   EXAMPLE
0023 %   % load example data
0024 %   data = load('example/example_data/33029_daily.mat');
0025 %   Q = data.Q;
0026 %   t = data.t;
0027 %   flow_section = util_RecessionSegments(Q,t,'plot_results',true); % get recession segments
0028 %   [dQdt, Qm, flow_section, R2] = util_dQdt(Q, t, flow_section); % get flow rate gradient
0029 %   rec = ~isnan(Qm);
0030 %   [a, b] = util_FitPowerLaw(Qm(rec), dQdt(rec));
0031 %   [a, b] = util_FitPowerLaw(Qm(rec), dQdt(rec), 'fitting_type', 'linear', 'weights', R2(rec))
0032 %
0033 %   Copyright (C) 2020
0034 %   This software is distributed under the GNU Public License Version 3.
0035 %   See <https://www.gnu.org/licenses/gpl-3.0.en.html> for details.
0036 
0037 if nargin < 2
0038     error('Not enough input arguments.')
0039 elseif nargin < 3
0040     fitting_type = 'linear';
0041 elseif nargin < 4
0042     weights = ones(size(Q_rec));
0043 end
0044 
0045 % default setting reads as good data
0046 error_flag = 0;
0047 error_str = '';
0048 
0049 % data checks
0050 if all(isnan(Q_rec)) || all(isnan(dQdt_rec))
0051     a = NaN;
0052     b = NaN;
0053     error_flag = 1;
0054     error_str = ['Warning: Some recessions consist only of NaN values (possibly because eps > 0). ', error_str];
0055     return
0056 end
0057 
0058 if strcmp(fitting_type, 'linear') % weighted linear regression in log log space
0059     
0060     A = [ones(size(Q_rec)), log(Q_rec)];
0061     P = (weights.*A)\(weights.*log(-dQdt_rec));
0062     a = exp(P(1));
0063     b = P(2);
0064     
0065     %{
0066     linFcn = @(p,x) p(2).*x + p(1);
0067     p0 = [0.1, 1.0];
0068     fit_nonlin = fitnlm(log(Q_rec),log(-dQdt_rec),linFcn,p0,'Weight',weights);
0069     a = exp(fit_nonlin.Coefficients.Estimate(1));
0070     b = fit_nonlin.Coefficients.Estimate(2);
0071     %}
0072     
0073 elseif strcmp(fitting_type, 'nonlinear') % weighted nonlinear regression in lin space
0074     
0075     powerFcn = @(p,x) p(1).*x.^p(2);
0076     p0 = [0.1, 1.0];
0077     fit_nonlin = fitnlm(Q_rec,-dQdt_rec,powerFcn,p0,'Weight',weights);
0078     a = fit_nonlin.Coefficients.Estimate(1);
0079     b = fit_nonlin.Coefficients.Estimate(2);
0080     
0081 elseif strcmp(fitting_type, 'slope1') % linear regression in log log space with fixed slope 1
0082     
0083     A = [ones(size(Q_rec))];
0084     b = 1;
0085     B = log(-dQdt_rec) - b*log(Q_rec);
0086     a = exp((weights.*A)\(weights.*B));
0087     
0088 elseif strcmp(fitting_type, 'slope2') % linear regression in log log space with fixed slope 2
0089     
0090     A = [ones(size(Q_rec))];
0091     b = 2;
0092     B = log(-dQdt_rec) - b*log(Q_rec);
0093     a = exp((weights.*A)\(weights.*B));
0094     
0095 else
0096     error('Invalid fitting type.')
0097 end
0098 
0099 end
0100 
0101