Improve docstrings. Fix small bugs.

xicor.m

@@ -1,9 +1,9 @@
-function [xi,p] = xicor(x,y,varargin)
+function [xi, p] = xicor(x, y, varargin)
 %XICOR Computes Chaterjee's xi correlation between x and y variables
 %
-%   [xi,p] = xicor(x,y,method)
-%   Returns the xi-correlationa as well as the corresponding p-value
-%   for the pair of variables x and y.
+%   [xi, p] = xicor(x, y)
+%   Returns the xi-correlation with the corresponding p-value for the pair
+%   of variables x and y.
 %
 %   Input arguments:
 %  
@@ -14,11 +14,17 @@
 %
 %   Name-value arguments:
 %  
-%   'method'         Method to be used to compute the correlation. Options
-%                    Options are
+%   'symmetric'      If true xi is computed as (r(x,y)+r(y,x))/2.
+%                    Default: false.
+%
+%   'p_val_method'   Method to be used to compute the p-value.
+%                    Options: 'theoretical' or 'permutation'.
+%                    Default: 'theoretical'.
+%
+%   'n_perm'         Number of permutations when p_val_method is
+%                    'permutation'. 
+%                    Default: 1000.
 %
-%   'symmetric'      If true xi is computed as (r(x,y)+r(y,x))/2. Default is 
-%                    false.
 %  
 %   Output arguments:
 %  
@@ -27,10 +33,12 @@
 %   'p'              Estimated p-value.
 %
 %
-%   
 %   Notes
 %   -----
-%   The xi-correlation is not symmetric. 
+%   This is an independent implementation of the method largely based on 
+%   the R-package developed by the original authors [3].
+%   The xi-correlation is not symmetric by default. 
+%   Check [2] for a potential improvement over the current implementation.
 %
 %
 %   References
@@ -42,6 +50,9 @@
 %   [2] Zhexiao Lin* and Fang Han†, On boosting the power of Chatterjee’s
 %   rank correlation, arXiv, 2021. https://arxiv.org/abs/2108.06828
 %
+%   [3] XICOR R package. 
+%   https://cran.r-project.org/web/packages/XICOR/index.html
+%
 %
 %   Example
 %   ---------      
@@ -55,38 +66,37 @@
 %   David Romero-Bascones, [email protected]
 %   Biomedical Engineering Department, Mondragon Unibertsitatea, 2022
 
-% Initial checks
 if nargin == 1
-    error('err1:MoreInputsRequired','xicor requires at least two inputs');
+    error('err1:MoreInputsRequired', 'xicor requires at least 2 inputs.');
 end
 
 parser = inputParser;
-addRequired(parser,'x');
-addRequired(parser,'y');
-addOptional(parser,'symmetric',false)
-addOptional(parser,'p_value_method','theoretical')
-addOptional(parser,'n_perm',100)
+addRequired(parser, 'x');
+addRequired(parser, 'y');
+addOptional(parser, 'symmetric', false)
+addOptional(parser, 'p_val_method', 'theoretical')
+addOptional(parser, 'n_perm', 1000)
 
 parse(parser,x,y,varargin{:})
 
 x = parser.Results.x;
 y = parser.Results.y;
 symmetric = parser.Results.symmetric;
-p_value_method = parser.Results.p_value_method;
+p_val_method = parser.Results.p_val_method;
 n_perm = parser.Results.n_perm;
 
 if ~isnumeric(x) || ~isnumeric(y)
-    error('err2:TypeError','x and y are must be numeric');
+    error('err2:TypeError', 'x and y are must be numeric.');
 end
 
 n = length(x);
 
 if n ~= length(y)
-    error('err3:IncorrectLength','x and y must have the same length');
+    error('err3:IncorrectLength', 'x and y must have the same length.');
 end
 
 if ~islogical(symmetric)
-    error('err2:TypeError','symmetric must be true or false');
+    error('err2:TypeError', 'symmetric must be true or false.');
 end
 
 % Check for NaN values
@@ -105,61 +115,62 @@
 
 if n < 10
     warning(['Running xicor with only ', num2str(n),...
-        ' points. This might produce unstable results.']);
+             ' points. This might produce unstable results.']);
 end
 
 [xi, r, l] = compute_xi(x, y);
 
 if symmetric
-    xi = (xi + compute_xi(y,x))/2;
+    xi = (xi + compute_xi(y, x))/2;
 end
 
 % If only one output return xi
 if nargout <= 1
     return
 end
 
-if strcmp(p_value_method, 'numeric') && symmetric==true
-    p_value_method = 'permutation';
+if ~strcmp(p_val_method, 'permutation') && symmetric==true
+    error('err2:TypeError', ...
+          'p_val_method when symmetric==true must be permutation.');
 end
 
 % Compute p-values (only valid for large n)
-switch p_value_method
+switch p_val_method
     case 'theoretical'
         if length(unique(y)) == n
-            p = 1 - normcdf(sqrt(n)*xi,0,sqrt(2/5));                
+            p = 1 - normcdf(sqrt(n)*xi, 0, sqrt(2/5));                
         else
             u = sort(r);
             v = cumsum(u);
             i = 1:n;
 
-            a = 1/n^4*sum((2*n -2*i +1) .* u.^2);
-            b = 1/n^5*sum((v + (n - i) .* u).^2);
-            c = 1/n^3*sum((2*n -2*i +1) .* u);
-            d = 1/n^3*sum(l .* (n - l));
+            a = 1/n^4 * sum((2*n -2*i +1) .* u.^2);
+            b = 1/n^5 * sum((v + (n - i) .* u).^2);
+            c = 1/n^3 * sum((2*n -2*i +1) .* u);
+            d = 1/n^3 * sum(l .* (n - l));
 
             tau = sqrt((a - 2*b + c^2)/d^2);
 
-            p = 1 - normcdf(sqrt(n)*xi,0,tau);
+            p = 1 - normcdf(sqrt(n)*xi, 0, tau);
         end
     case 'permutation'
-        xi_perm = nan(1,n_perm);
+        xi_perm = nan(1, n_perm);
 
         if symmetric
             for i_perm=1:n_perm
-                xi_perm(i_perm) = compute_xi(x(rand_perm(n)),y);
+                x_perm = x(randperm(n));
+                xi_perm(i_perm) = (compute_xi(x_perm, y) + ...
+                                   compute_xi(y, x_perm))/2;
             end
         else
             for i_perm=1:n_perm
-                x_perm = x(randperm(n));
-                xi_perm(i_perm) = (compute_xi(x_perm,y) + ...
-                                   compute_xi(y,x_perm))/2;
+                xi_perm(i_perm) = compute_xi(x(randperm(n)), y);
             end            
         end
 
         p = sum(xi_perm > xi)/n_perm;
     otherwise
-        error("Wrong p_value_method. Use 'numeric' or 'permutation'");        
+        error("Wrong p_value_method. Use 'theoretical' or 'permutation'");        
 end
 
 function [xi, r, l] = compute_xi(x,y)