fix few errors

test_xicor.m

@@ -3,20 +3,63 @@
 x = rand(1,100);
 y = rand(1,100);
 xi = xicor(x,y);
-if abs(xi) < 0.3
-    disp('Verification passed.');
-else
-    warning('Verification failed');
-end
+my_assert(abs(xi) < 0.3);
+
 %% Basic test: totally dependent
 x = randn(1,100);
 y = x.^2;
 xi = xicor(x,y);
-assert(xi > 0.9);
+my_assert(xi > 0.9);
 
 %% Break tests
 testCase = matlab.unittest.TestCase.forInteractiveUse;
 verifyError(testCase,@() xicor('0'),"err1:MoreInputsRequired")
 verifyError(testCase,@() xicor('0',2),"err2:TypeError")
 verifyError(testCase,@() xicor([1 2],[2 4 5]),"err3:IncorrectLength")
-verifyError(testCase,@() xicor([1 2 3],[2 4 5],'symmetric',2),"err2:TypeError")
+verifyError(testCase,@() xicor([1 2 3],[2 4 5],'symmetric',2),"err2:TypeError")
+
+%% With - without ties comparison
+n = 4000;
+x = 1:n;
+y = cos(x./30) + 4*randn(1,n);
+[xi, p] = xicor(x,y);
+
+idx = randi(n);
+x2 = [x x(idx)];
+y2 = [y y(idx)];
+[xi2, p2] = xicor(x2,y2);
+
+my_assert(abs(xi - xi2) < 0.1);
+
+%% Theoretical p_value vs. null
+n = 4000;
+n_boot = 10000;
+
+x = 1:n;
+y = cos(x./30) + 5*randn(1,n);
+
+[xi, p] = xicor(x,y);
+
+x_boot = randn(n_boot, n);
+y_boot = randn(n_boot, n);
+xi_boot = nan(1,n_boot);
+for i_boot=1:n_boot
+   xi_boot(i_boot) = xicor(x_boot(i_boot,:),y_boot(i_boot,:));
+end
+
+p_boot = sum(xi_boot > xi)/n_boot;
+
+subplot(121);plot(x,y);
+subplot(122);
+h = histogram(sqrt(n)*xi_boot);hold on;
+plot(xi*sqrt(n)*[1 1],[0 max(h.Values)],'r','LineWidth',1.5);
+
+my_assert(abs(p - p_boot) < 0.1);
+%% Basic assert function
+function my_assert(cond)
+if cond
+    disp('Verification passed.');
+else
+    warning('Verification failed');
+end
+end

xicor.m

@@ -136,20 +136,11 @@
             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));
+            d = 1/n^3*sum(l .* (n - l));
 
-            t = (a - 2*b + c^2)/d^2;
+            tau = sqrt((a - 2*b + c^2)/d^2);
 
-            p = 1 - normcdf(sqrt(n)*xi,0,sqrt(t));
-%            qr = sort(r);            
-%            ind = 1:n;
-%            ind2 = 2*n - 2*ind + 1;
-%            a = mean(ind2*qfr*qfr)/n;
-%            c = mean(ind2*qfr)/n;
-%            cq = cumsum(qfr);
-%            m = (cq + (n - ind)*qfr)/n;
-%            b = mean(m^2);
-%            v = (ai - 2*b + ci^2)/(CU^2);
+            p = 1 - normcdf(sqrt(n)*xi,0,tau);
         end
     case 'permutation'
         xi_perm = nan(1,n_perm);
@@ -160,7 +151,7 @@
             end
         else
             for i_perm=1:n_perm
-                x_perm = x(rand_perm(n));
+                x_perm = x(randperm(n));
                 xi_perm(i_perm) = (compute_xi(x_perm,y) + ...
                                    compute_xi(y,x_perm))/2;
             end