improve numerical accuracy when solving linear systems, correct dates…

… modified and update matlab toolbox too

DESCRIPTION

@@ -1,6 +1,6 @@
 name: statistics-resampling
-version: 5.5.12
-date: 2024-05-15
+version: 5.5.13
+date: 2024-05-17
 author: Andrew Penn <[email protected]>
 maintainer: Andrew Penn <[email protected]>
 title: A statistics package with a variety of resampling tools

inst/boot.m

@@ -65,7 +65,7 @@
 %        vs. Smoothing; Proceedings of the Section on Statistics & the 
 %        Environment. Alexandria, VA: American Statistical Association.
 %
-%  boot (version 2024.04.16)
+%  boot (version 2024.04.24)
 %  Author: Andrew Charles Penn
 %  https://www.researchgate.net/profile/Andrew_Penn/
 %

inst/boot1way.m

@@ -121,7 +121,7 @@
 %  [1] Efron, and Tibshirani (1993) An Introduction to the Bootstrap. 
 %        New York, NY: Chapman & Hall
 %
-%  boot1way (version 2023.10.04)
+%  boot1way (version 2024.04.24)
 %  Bootstrap tests for comparing independent groups in a one-way layout
 %  Author: Andrew Charles Penn
 %  https://www.researchgate.net/profile/Andrew_Penn/

inst/bootbayes.m

@@ -124,7 +124,7 @@
 %  [3] Liu, Gelman & Zheng (2015). Simulation-efficient shortest probability
 %        intervals. Statistics and Computing, 25(4), 809–819. 
 %
-%  bootbayes (version 2023.07.06)
+%  bootbayes (version 2024.05.17)
 %  Author: Andrew Charles Penn
 %  https://www.researchgate.net/profile/Andrew_Penn/
 %
@@ -374,9 +374,9 @@
     bootstat = cell2mat (cellfun (bootfun, num2cell (Y, 1)', ...
                                  'UniformOutput', false));
   else
-    bootfun = @(w) lmfit (X, Y, diag (w), L);
+    bootfun = @(w) lmfit (X, Y, w, L);
     original = bootfun (ones (n, 1) / n);
-    bootstat = cell2mat (cellfun (bootfun, num2cell (W, 1), ...
+    bootstat = cell2mat (cellfun (bootfun, num2cell (sqrt (W), 1), ...
                                   'UniformOutput', false));
   end
 
@@ -419,25 +419,30 @@
 
 % FUNCTION TO FIT THE LINEAR MODEL
 
-function param = lmfit (X, y, W, L)
+function param = lmfit (X, y, w, L)
 
-  % Get model coefficients by solving the linear equation by matrix arithmetic.
+  % Get model coefficients by solving the linear equation.
   % Solve the linear least squares problem using the Moore-Penrose pseudo
   % inverse (pinv) to make it more robust to the situation where X is singular.
-  % If optional arument W is provided, it should be a diagonal matrix of
-  % weights or a positive definite covariance matrix
+  % If optional arument w is provided, it should be equal to square root of the
+  % weights we want to apply to the regression.
   n = numel (y);
-  if ( (nargin < 3) || isempty (W) )
-    % If no weights are provided, create an identity matrix
-    W = eye (n);
+  if ( (nargin < 3) || isempty (w) )
+    % If no weights are provided, create a vector of ones
+    w = ones (n, 1);
   end
   if ( (nargin < 4) || isempty (L) )
     % If no hypothesis matrix (L) is provided, set L to 1
     L = 1;
   end
 
-  % Solve linear equation to minimize weighted least squares
-  b = pinv (X' * W * X) * (X' * W * y);
+  % Solve the linear equation to minimize weighted least squares, where the
+  % weights are equal to w.^2. The whitening transformation actually implemented
+  % avoids using an n * n matrix and has improved accuracy over the solution
+  % using the equivalent normal equation:
+  %   b = pinv (X' * W * X) * (X' * W * y);
+  % Where W is the diagonal matrix of weights (i.e. W = diag (w.^2))
+  b = pinv (bsxfun (@times, w, X)) * (w .* y);
   param = L' * b;
 
 end
@@ -454,7 +459,7 @@ function print_output (stats, nboot, prob, prior, p, L, method, intercept_only)
                      '******************************\n\n'));
     fprintf ('Bootstrap settings: \n');
     if (intercept_only)
-        fprintf (' Function: sum (w .* Y)\n');
+        fprintf (' Function: sum (W * Y)\n');
     else
       if ( (numel(L) > 1) || (L ~= 1) )
         fprintf (' Function: L'' * pinv (X'' * W * X) * (X'' * W * y)\n');

inst/bootcdf.m

@@ -27,7 +27,7 @@
 %
 %     '[x, F, P] = bootcdf (...)' also returns the distribution of P values.
 %
-%  bootcdf (version 2023.07.05)
+%  bootcdf (version 2024.04.21)
 %  Author: Andrew Charles Penn
 %  https://www.researchgate.net/profile/Andrew_Penn/
 %

inst/bootci.m

@@ -131,7 +131,7 @@
 %        vs. Smoothing; Proceedings of the Section on Statistics & the 
 %        Environment. Alexandria, VA: American Statistical Association.
 %
-%  bootci (version 2023.07.04)
+%  bootci (version 2024.05.13)
 %  Author: Andrew Charles Penn
 %  https://www.researchgate.net/profile/Andrew_Penn/
 %

inst/bootclust.m

@@ -123,7 +123,7 @@
 %        vs. Smoothing; Proceedings of the Section on Statistics & the 
 %        Environment. Alexandria, VA: American Statistical Association.
 %
-%  bootclust (version 2024.04.02)
+%  bootclust (version 2024.05.16)
 %  Author: Andrew Charles Penn
 %  https://www.researchgate.net/profile/Andrew_Penn/
 %

inst/bootknife.m

@@ -159,7 +159,7 @@
 %  [9] Gleason, J.R. (1988) Algorithms for Balanced Bootstrap Simulations. 
 %        The American Statistician. Vol. 42, No. 4 pp. 263-266
 %
-%  bootknife (version 2023.07.04)
+%  bootknife (version 2024.05.15)
 %  Author: Andrew Charles Penn
 %  https://www.researchgate.net/profile/Andrew_Penn/
 %

inst/bootlm.m

@@ -474,7 +474,7 @@
 %       installed and loaded, then these computations will be automatically
 %       accelerated by parallel processing on platforms with multiple processors
 %
-%  bootlm (version 2023.09.01)
+%  bootlm (version 2024.05.17)
 %  Author: Andrew Charles Penn
 %  https://www.researchgate.net/profile/Andrew_Penn/
 %
@@ -1732,7 +1732,7 @@
   % Get model coefficients by solving the linear equation. The number of free
   % parameters (i.e. intercept + coefficients) is equal to n - dfe (i.e. the
   % number of columns in X).
-  b = X \ Y;                 % Equivalent to inv (X' * X) * (X' * y);
+  b = pinv (X) * Y;                 % Equivalent to inv (X' * X) * (X' * y);
 
   % Get fitted values
   fit = X * b;
@@ -1746,12 +1746,8 @@
   % Calculate the unscaled covariance matrix (i.e. inv (X'*X )) and the Hat
   % matrix (i.e. X*(X'*X)^−1*X') by QR decomposition
   if (nargout > 3)
-    [Q, R] = qr (X, 0);      % Economy-sized QR decomposition
-    if ISOCTAVE
-      ucov = chol2inv (R);
-    else
-      ucov = inv (R' * R);
-    end
+    [Q, R] = qr (X, 0);     % Economy-sized QR decomposition
+    ucov = pinv (R' * R);   % Instead of pinv (X' * X)
     hat = Q * Q';
   end
 

inst/bootstrp.m

@@ -106,7 +106,7 @@
 %        vs. Smoothing; Proceedings of the Section on Statistics & the 
 %        Environment. Alexandria, VA: American Statistical Association.
 %
-%  bootstrp (version 2024.04.23)
+%  bootstrp (version 2024.05.15)
 %  Author: Andrew Charles Penn
 %  https://www.researchgate.net/profile/Andrew_Penn/
 %

inst/bootwild.m

@@ -112,7 +112,7 @@
 %  [7] Sellke, Bayarri and Berger (2001) Calibration of p-values for Testing
 %        Precise Null Hypotheses. Am Stat. 55(1), 62-71
 %
-%  bootwild (version 2023.07.05)
+%  bootwild (version 2024.05.17)
 %  Author: Andrew Charles Penn
 %  https://www.researchgate.net/profile/Andrew_Penn/
 %
@@ -282,13 +282,12 @@
     rand ('seed', seed);
   end
 
-  % Compute unscaled covariance matrix
+  % Compute unscaled covariance matrix by QR decomposition (instead of using
+  % the less accurate method of getting to the solution directly with normal
+  % equations)
   [Q, R] = qr (X, 0);        % Economy-sized QR decomposition
-  if ISOCTAVE
-    ucov = chol2inv (R);     % Instead of pinv (X' * X). Faster
-  else
-    ucov = inv (R' * R);     % Instead of pinv (X' * X). Slower
-  end
+  ucov = pinv (R' * R);      % Instead of pinv (X' * X)
+
 
   % Create least squares anonymous function for bootstrap
   bootfun = @(y) lmfit (X, y, ucov, clusters, L, ISOCTAVE);
@@ -306,7 +305,7 @@
   if (~ isempty (IC))
     s = s(IC, :);  % Enforce clustering/blocking
   end
-  yf = X * (X \ y);
+  yf = X * pinv (X) * y;
   r = y - yf;
   rs = bsxfun (@times, r, s);
   Y = bsxfun (@plus, yf, rs);
@@ -383,7 +382,7 @@
 
   % Solve linear equation to minimize least squares and compute the
   % regression coefficients (b) 
-  b = X \ y;                 % Instead of inv (X' * X) * (X' * y);
+  b = pinv (X) * y;                 % Instead of inv (X' * X) * (X' * y);
 
   % Calculate heteroscedasticity-consistent (HC) or cluster robust (CR) standard 
   % errors for the regression coefficients. When the number of observations
@@ -395,10 +394,10 @@
   yf = X * b;
   u = y - yf;
   if ( (nargin < 3) || isempty (clusters) )
-    % Heteroscedasticity-Consistent (HC0) standard errors
+    % For Heteroscedasticity-Consistent (HC0) standard errors
     meat = X' * diag (u.^2) * X;
   else
-    % Cluster Robust (CR0) standard errors
+    % For Cluster Robust (CR0) standard errors
     Sigma = cellfun (@(g) X(g,:)' * u(g) * u(g)' * X(g,:), ...
                   num2cell (clusters, 1), 'UniformOutput', false);
     meat = sum (cat (3, Sigma{:}), 3);
@@ -458,9 +457,9 @@ function print_output (stats, nboot, alpha, p, L, method)
                      '************************************\n\n'));
     fprintf ('Bootstrap settings: \n');
     if ( (numel(L) > 1) || (L ~= 1) )
-      fprintf (' Function: L'' * X \\ y\n');
+      fprintf (' Function: L'' * pinv (X) * y\n');
     else
-      fprintf (' Function: X \\ y\n');
+      fprintf (' Function: pinv (X) * y\n');
     end
     fprintf (' Resampling method: Wild %sbootstrap-t\n', method)
     fprintf (' Number of resamples: %u \n', nboot)

matlab/statistics-resampling.prj

@@ -1,13 +1,13 @@
 <deployment-project plugin="plugin.toolbox" plugin-version="1.0">
-  <configuration build-checksum="3003727527" file="Y:\Documents\GitHub\statistics-resampling\matlab\statistics-resampling.prj" location="Y:\Documents\GitHub\statistics-resampling\matlab" name="statistics-resampling" target="target.toolbox" target-name="Package Toolbox">
+  <configuration build-checksum="3022689601" file="Y:\Documents\GitHub\statistics-resampling\matlab\statistics-resampling.prj" location="Y:\Documents\GitHub\statistics-resampling\matlab" name="statistics-resampling" target="target.toolbox" target-name="Package Toolbox">
     <param.appname>statistics-resampling</param.appname>
     <param.authnamewatermark>Andrew Penn</param.authnamewatermark>
     <param.email>[email protected]</param.email>
     <param.company>University of Sussex, UK</param.company>
     <param.summary>Statistical analysis using resampling methods</param.summary>
     <param.description>The statistics-resampling package is an Octave package and Matlab toolbox that can be used to perform a wide variety of statistics tasks using non-parametric resampling methods. In particular, the functions included can be used to estimate bias, uncertainty (standard errors and confidence intervals), prediction error, and calculate p-values for null hypothesis significance tests. Variations of the resampling methods are included that improve the accuracy of the statistics for small samples and samples with complex dependence structures.</param.description>
     <param.screenshot>Y:\Documents\GitHub\statistics-resampling\doc\icon.png</param.screenshot>
-    <param.version>5.5.12</param.version>
+    <param.version>5.5.13</param.version>
     <param.output>${PROJECT_ROOT}\statistics-resampling.mltbx</param.output>
     <param.products.name />
     <param.products.id />