1## Copyright (C) 2021 Nicholas R. Jankowski <jankowskin@asme.org> 2## 3## This program is free software; you can redistribute it and/or modify it under 4## the terms of the GNU General Public License as published by the Free Software 5## Foundation; either version 3 of the License, or (at your option) any later 6## version. 7## 8## This program is distributed in the hope that it will be useful, but WITHOUT 9## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 10## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 11## details. 12## 13## You should have received a copy of the GNU General Public License along with 14## this program; if not, see <http://www.gnu.org/licenses/>. 15 16## -*- texinfo -*- 17## @deftypefn {Function File} {@var{mu} =} expfit (@var{s}) 18## @deftypefnx {Function File} {[@var{mu}, @var{ci}] =} expfit (@var{s}) 19## @deftypefnx {Function File} {[@var{mu}, @var{ci}] =} expfit (@var{s}, @var{alpha}) 20## @deftypefnx {Function File} {@dots{} =} expfit (@var{s}, @var{alpha}, @var{c}) 21## @deftypefnx {Function File} {@dots{} =} expfit (@var{s}, @var{alpha}, @var{c}, @var{f}) 22## 23## Estimate the mean of the exponential probability distribution function from 24## which sample data @var{s} has been taken. @var{s} is expected to be a 25## non-negative vector. If @var{s} is an array, the mean will be computed for 26## each column of @var{s}. If any elements of @var{s} are NaN, that vector's 27## mean will be returned as NaN. 28## 29## If the optional output variable @var{ci} is requested, @code{expfit} will 30## also return the confidence interval bounds for the estimate as a two element 31## column vector. If @var{s} is an array, each column of data will have a 32## confidence interval returned as a two row array. 33## 34## The optional scalar input @var{alpha} can be used to define the 35## (1-@var{alpha}) confidence interval to be applied to all estimates as a 36## value between 0 and 1. The default is 0.05, resulting in a 0.95 or 95% CI. 37## Any invalid values for alpha will return NaN for both CI bounds. 38## 39## The optional input @var{c} is a logical or numeric array of zeros and ones 40## the same size as @var{s}, used to right-censor individual elements of 41## @var{s}. A value of 1 indicates the data should be censored from 42## the mean estimation. Any nonzero values in @var{c} are treated as a 1. 43## 44## The optional input @var{f} is a numeric array the same size as @var{s}, used 45## to specify occurrence frequencies for the elements in @var{s}. Values of 46## @var{f} need not be integers. Any NaN elements in the frequency array will 47## produce a NaN output for @var{mu}. 48## 49## Options can be skipped by using [] to revert to the default. 50## 51## Matlab incompatibility: Matlab's @code{expfit} produces unpredictable results 52## for some cases with higher dimensions (specifically 1 x m x n x ... arrays). 53## Octave's implementation allows for n-D arrays, consistently performing 54## calculations on individual column vectors. Additionally, @var{c} and @var{f} 55## can be used with arrays of any size, whereas Matlab only allows their use 56## when @var{s} is a vector. 57## 58## @end deftypefn 59## 60## @seealso{expcdf, expinv, explike, exppdf, exprnd, expstat} 61## @seealso{expstat, exprnd, expcdf, expinv} 62 63function [m, v] = expfit (s, alpha = 0.05, c = [], f = []) 64 65 ## Check arguments 66 if (nargin ==0 || nargin > 4 || nargout > 2) 67 print_usage (); 68 endif 69 70 if ! (isnumeric (s) || islogical (s)) 71 s = double(s); 72 endif 73 74 ## guarantee working with column vectors 75 if isvector (s) 76 77 s = s(:); 78 endif 79 80 if any (s(:) < 0) 81 error("expfit: input data S cannot be negative"); 82 endif 83 84 sz_s = size (s); 85 86 if (isempty (alpha)) 87 alpha = 0.05; 88 elseif !(isscalar (alpha)) 89 error ("expfit: ALPHA must be a scalar quantity"); 90 endif 91 92 if (isempty (c) && isempty (f)) 93 ##simple case without f or c, shortcut other validations 94 m = mean (s, 1); 95 96 if (nargout == 2) 97 S = sum (s, 1); 98 v = [2*S ./ chi2inv(1 - alpha/2, 2*sz_s(1));... 99 2*S ./ chi2inv(alpha/2, 2*sz_s(1))]; 100 endif 101 else 102 103 ## input validation for c and f 104 105 if (isempty (c)) 106 ##expand to full c with values that don't affect results 107 c = zeros (sz_s); 108 elseif (! (isnumeric(c) || islogical (c))) 109 #check for incorrect f type 110 error ("expfit: C must be a numeric or logical array") 111 elseif (isvector (c)) 112 ## guarantee working with a column vector 113 c = c(:); 114 endif 115 116 if (isempty (f)) 117 ##expand to full c with values that don't affect results 118 f = ones (sz_s); 119 elseif (! (isnumeric(f) || islogical (f))) 120 #check for incorrect f type 121 error ("expfit: F must be a numeric or logical array") 122 elseif (isvector (f)) 123 ## guarantee working with a column vector 124 f = f(:); 125 endif 126 127 #check that size of c and f match s 128 if !(isequal(size (c), sz_s)) 129 error("expfit: C must be the same size as S"); 130 elseif (! isequal(size (f), sz_s)) 131 error("expfit: F must be the same size as S"); 132 endif 133 134 ## trivial case where c and f have no effect 135 if (all (c(:) == 0 & f(:) == 1)) 136 137 m = mean (s, 1); 138 139 if (nargout == 2) 140 S = sum (s, 1); 141 v = [2*S ./ chi2inv(1 - alpha/2, 2*sz_s(1));... 142 2*S ./ chi2inv(alpha/2, 2*sz_s(1))]; 143 endif 144 145 ## no censoring, just adjust sample counts for f 146 elseif (all (c(:) == 0)) 147 148 S = sum (s.*f, 1); 149 n = sum (f, 1); 150 m = S ./ n; 151 152 if (nargout == 2) 153 v = [2*S ./ chi2inv(1 - alpha/2, 2*n);... 154 2*S ./ chi2inv(alpha/2, 2*n)]; 155 endif 156 157 ## censoring, but no sample counts adjustment 158 elseif (all (f(:) == 1)) 159 160 c = logical(c); ##convert any numeric c's to 0s and 1s 161 S = sum (s, 1); 162 r = sz_s(1) - sum (c, 1); 163 m = S ./ r; 164 165 if (nargout == 2) 166 v = [2*S ./ chi2inv(1 - alpha/2, 2*r);... 167 2*S ./ chi2inv(alpha/2, 2*r)]; 168 endif 169 170 ## both censoring and sample count adjustment 171 else 172 173 c = logical(c); ##convert any numeric c's to 0s and 1s 174 S = sum (s.*f , 1); 175 r = sum (f.*(!c), 1); 176 m = S ./ r; 177 178 if (nargout == 2) 179 v = [2*S ./ chi2inv(1 - alpha/2, 2*r);... 180 2*S ./ chi2inv(alpha/2, 2*r)]; 181 endif 182 endif 183 184 ## compatibility check, NaN for columns where all c's or f's remove all samples 185 null_columns = all (c) | ! all (f); 186 m(null_columns) = NaN; 187 188 if (nargout == 2) 189 v(:,null_columns) = NaN; 190 endif 191 endif 192 193endfunction 194 195##tests for mean 196%!assert (expfit (1), 1) 197%!assert (expfit (1:3), 2) 198%!assert (expfit ([1:3]'), 2) 199%!assert (expfit (1:3, []), 2) 200%!assert (expfit (1:3, [], [], []), 2) 201%!assert (expfit (magic (3)), [5 5 5]) 202%!assert (expfit (cat (3, magic (3), 2*magic (3))), cat (3,[5 5 5], [10 10 10])) 203%!assert (expfit (1:3, 0.1, [0 0 0], [1 1 1]), 2) 204%!assert (expfit ([1:3]', 0.1, [0 0 0]', [1 1 1]'), 2) 205%!assert (expfit (1:3, 0.1, [0 0 0]', [1 1 1]'), 2) 206%!assert (expfit (1:3, 0.1, [1 0 0], [1 1 1]), 3) 207%!assert (expfit (1:3, 0.1, [0 0 0], [4 1 1]), 1.5) 208%!assert (expfit (1:3, 0.1, [1 0 0], [4 1 1]), 4.5) 209%!assert (expfit (1:3, 0.1, [1 0 1], [4 1 1]), 9) 210%!assert (expfit (1:3, 0.1, [], [-1 1 1]), 4) 211%!assert (expfit (1:3, 0.1, [], [0.5 1 1]), 2.2) 212%!assert (expfit (1:3, 0.1, [1 1 1]), NaN) 213%!assert (expfit (1:3, 0.1, [], [0 0 0]), NaN) 214%!assert (expfit (reshape (1:9, [3 3])), [2 5 8]) 215%!assert (expfit (reshape (1:9, [3 3]), [], eye(3)), [3 7.5 12]) 216%!assert (expfit (reshape (1:9, [3 3]), [], 2*eye(3)), [3 7.5 12]) 217%!assert (expfit (reshape (1:9, [3 3]), [], [], [2 2 2; 1 1 1; 1 1 1]), [1.75 4.75 7.75]) 218%!assert (expfit (reshape (1:9, [3 3]), [], [], [2 2 2; 1 1 1; 1 1 1]), [1.75 4.75 7.75]) 219%!assert (expfit (reshape (1:9, [3 3]), [], eye(3), [2 2 2; 1 1 1; 1 1 1]), [3.5 19/3 31/3]) 220 221##tests for confidence intervals 222%!assert ([~,v] = expfit (1:3, 0), [0; Inf]) 223%!assert ([~,v] = expfit (1:3, 2), [Inf; 0]) 224%!assert ([~,v] = expfit (1:3, 0.1, [1 1 1]), [NaN; NaN]) 225%!assert ([~,v] = expfit (1:3, 0.1, [], [0 0 0]), [NaN; NaN]) 226%!assert ([~,v] = expfit (1:3, -1), [NaN; NaN]) 227%!assert ([~,v] = expfit (1:3, 5), [NaN; NaN]) 228#!assert ([~,v] = expfit ([1:3;1:3], -1), NaN(2, 3)] 229#!assert ([~,v] = expfit ([1:3;1:3], 5), NaN(2, 3)] 230%!assert ([~,v] = expfit (1:3), [0.830485728373393; 9.698190330474096], 1000*eps) 231%!assert ([~,v] = expfit (1:3, 0.1), [0.953017262058213; 7.337731146400207], 1000*eps) 232%!assert ([~,v] = expfit ([1:3;2:4]), ... 233%! [0.538440777613095, 0.897401296021825, 1.256361814430554; ... 234%! 12.385982973214016, 20.643304955356694, 28.900626937499371], 1000*eps) 235%!assert ([~,v] = expfit ([1:3;2:4], [], [1 1 1; 0 0 0]), ... 236%! 100*[0.008132550920455, 0.013554251534091, 0.018975952147727; ... 237%! 1.184936706156216, 1.974894510260360, 2.764852314364504], 1000*eps) 238%!assert ([~,v] = expfit ([1:3;2:4], [], [], [3 3 3; 1 1 1]), ... 239%! [0.570302756652583, 1.026544961974649, 1.482787167296715; ... 240%! 4.587722594914109, 8.257900670845396, 11.928078746776684], 1000*eps) 241%!assert ([~,v] = expfit ([1:3;2:4], [], [0 0 0; 1 1 1], [3 3 3; 1 1 1]), ... 242%! [0.692071440311161, 1.245728592560089, 1.799385744809018; ... 243%! 8.081825275395081, 14.547285495711145, 21.012745716027212], 1000*eps) 244 245%!test 246%! s = reshape (1:8, [4 2]); 247%! s(4) = NaN; 248%! [m,v] = expfit (s); 249%! assert ({m, v}, {[NaN, 6.5], [NaN, 2.965574334593430;NaN, 23.856157493553368]}, 1000*eps); 250 251%!test 252%! s = magic (3); 253%! c = [0 1 0; 0 1 0; 0 1 0]; 254%! f = [1 1 0; 1 1 0; 1 1 0]; 255%! [m,v] = expfit (s, [], c, f); 256%! assert ({m, v}, {[5 NaN NaN], [[2.076214320933482; 24.245475826185242],NaN(2)]}, 1000*eps); 257 258## input validation 259%!error expfit () 260%!error expfit (1,2,3,4,5) 261%!error [a b c] = expfit (1) 262%!error <ALPHA must be a scalar quantity> expfit (1, [1 2]) 263%!error <input data S cannot be negative> expfit ([-1 2 3 4 5]) 264%!error <C must be a numeric or logical array> expfit ([1:5], [], "test") 265%!error <F must be a numeric or logical array> expfit ([1:5], [], [], "test") 266%!error <C must be the same size as S> expfit ([1:5], [], [0 0 0 0]) 267%!error <F must be the same size as S> expfit ([1:5], [], [], [1 1 1 1]) 268