Behavioral states affect eye position tuning in the parietal cortex of macaques.

1. function [rho, beta, theta, p] = vectorCorrelation(x, y, u, v, varargin)
2. % VECCORRELATION 2D vector correlation between two vector sets
3. %
4. % [rho, beta, theta, p] = VECCORRELATION(x, y, u, v)
5. % Computes the correlation between two paired 2D vector sets [x,y] and
6. % [u,v] using the approach from Hanson et al., 1992.
7. %
8. % Inputs
9. % x, y - Numeric vectors (Nx1 or 1xN) containing the first set of
10. % paired coordinates. x(i) and y(i) form the i-th 2D vector.
11. % u, v - Numeric vectors (same length as x and y) containing the
12. % second set of paired coordinates. u(i) and v(i) form the
13. % i-th 2D vector to compare to [x(i),y(i)].
14. % Name/value optional inputs (use quoted names):
15. % 'pvalue' - 'default' (compute p with t-transform) or 'shuffle'
16. % (compute empirical p by random shuffling). Default:
17. % 'default'.
18. % 'nShuffle' - Positive integer number of permutations used when
19. % 'pvalue' is 'shuffle'. Default: 1000.
20. %
21. % Outputs
22. % rho - Scalar correlation coefficient in [-1, +1]. Positive values
23. % indicate rotation-like alignment; negative values indicate
24. % reflection-like alignment.
25. % beta - Scale factor relating the magnitudes of the two vector sets.
26. % Example: Beta = 2 indicates that [x,y] vector fields must
27. % be scaled by a factor of 2 to match [u,v].
28. % theta - Rotation/reflection angle in degrees (computed via atand).
29. % p - Two-tailed p-value for rho. If 'pvalue' is 'default', p is
30. % computed from a t-transform. If 'pvalue' is 'shuffle', p
31. % is an empirical two-tailed p-value from the permutation
32. % null distribution (add-one correction applied).
33. %
34. % Notes
35. % - All input vectors must have equal length N >= 3 for the t-transform.
36. % - When using 'shuffle', the implementation applies the same random
37. % permutation to x and y (preserving their pairing) and a separate
38. % same random permutation to u and v. Each shuffle yields a rho; the
39. % empirical p is the fraction of |rho_null| >= |rho_observed|.
40. % - For reproducible permutations, set the RNG state before calling:
41. % rng(seed)
42. %
43. % Examples
44. % % Default p via t-transform
45. % [rho,beta,theta,p] = veccorrelation(x,y,u,v);
46. %
47. % % Empirical p with 5000 permutations
48. % [rho,beta,theta,p] = veccorrelation(x,y,u,v,'pvalue','shuffle','nShuffle',5000);
49. %
50. % References
51. % ----------
52. % Hanson, B., Klink, K., Matsuura, K., Robeson, S. M., & Willmott, C. J. (1992).
53. % Vector correlation: Review, exposition, and geographic application.
54. % Annals of the Association of American Geographers, 82(1), 103–116.
55. % https://doi.org/10.1111/j.1467-8306.1992.tb01900.x
56. %
57. % Please cite
58. % Hadjidimitrakis, K., Vaccari, F. E., De Vitis, M., Filippini, M., Diomedi, S., & Fattori, P. (2026).
59. % Spontaneous oculomotor behavior sharpens eye position signals in parietal cortex.
60. % [Manuscript under review].
62. narginchk(4, inf);
64. % Default options
65. opts.pvalue = 'default';
66. opts.nShuffle = 1000;
68. % Parse optional name/value pairs
69. if ~isempty(varargin)
70. if mod(length(varargin),2) ~= 0
71. error('Optional parameters must be provided as name/value pairs.');
72. end
73. for k = 1:2:length(varargin)
74. name = varargin{k};
75. val = varargin{k+1};
76. validateattributes(name, {'char','string'}, {'scalartext'});
77. switch lower(char(name))
78. case 'pvalue'
79. validateattributes(val, {'char','string'}, {'scalartext'});
80. valstr = lower(char(val));
81. if ~ismember(valstr, {'default','shuffle'})
82. error('''pvalue'' must be ''default'' or ''shuffle''.');
83. end
84. opts.pvalue = valstr;
85. case 'nshuffle'
86. validateattributes(val, {'numeric'}, {'scalar','integer','positive'});
87. opts.nShuffle = double(val);
88. otherwise
89. error('Unknown option name ''%s''.', char(name));
90. end
91. end
92. end
94. n = length(x);
95. if length(y) ~= n || length(u) ~= n || length(v) ~= n
96. error('All input vectors must have the same length.');
97. end
99. % Compute observed rho, beta, theta
100. [rho, beta, theta] = veccorrelation_core_compute(x, y, u, v);
102. switch opts.pvalue
103. case 'default'
104. % Default p-value by t-transform (two-tailed)
105. tstat = rho * sqrt(length(x) - 2) / sqrt(1 - rho^2);
106. p = 2 * (1 - tcdf(abs(tstat), length(x) - 2)); % two-tailed
107. case 'shuffle'
108. % Shuffle-based p-value
111. nShuffle = opts.nShuffle;
112. rho_null = zeros(nShuffle, 1);
114. % For reproducibility a user can set rng before calling this function.
115. for i = 1:nShuffle
116. perm_uv = randperm(n); % same permutation for u and v
117. u_sh = u(perm_uv);
118. v_sh = v(perm_uv);
119. [rho_i, ~, ~] = veccorrelation_core_compute(x, y, u_sh, v_sh);
120. rho_null(i) = rho_i;
121. end
123. % Two-tailed empirical p-value: proportion of |rho_null| >= |rho_observed|
124. count_extreme = sum(abs(rho_null) >= abs(rho));
125. p = (count_extreme + 1) / (nShuffle + 1); % add-one correction
127. end
129. end
131. function [rho, beta, theta] = veccorrelation_core_compute(x, y, u, v)
132. % Core computations (returns rho, beta, theta)
133. sigmax = std(x, 1);
134. sigmay = std(y, 1);
135. sigmau = std(u, 1);
136. sigmav = std(v, 1);
138. tmp = cov(x, u, 1);
139. sigmaxu = tmp(1,2);
140. tmp = cov(x, v, 1);
141. sigmaxv = tmp(1,2);
142. tmp = cov(y, u, 1);
143. sigmayu = tmp(1,2);
144. tmp = cov(y, v, 1);
145. sigmayv = tmp(1,2);
147. ksi = (sigmaxu * sigmayv) - (sigmaxv * sigmayu);
148. % handle degenerate ksi == 0 (avoid division by zero)
149. if ksi == 0
150. s = 0;
151. disp("WARNING: (sigmaxu * sigmayv) - (sigmaxv * sigmayu) == 0, computations might be altered")
152. else
153. s = ksi / abs(ksi);
154. end
156. a = sigmaxu^2 + sigmayv^2 + sigmaxv^2 + sigmayu^2 + (2 * s * ksi);
157. b = (sigmax^2 + sigmay^2) * (sigmau^2 + sigmav^2);
159. if b <= 0
160. rho = NaN;
161. else
162. rho = s * sqrt(a / b);
163. end
165. if (sigmax^2 + sigmay^2) <= 0
166. beta = NaN;
167. else
168. beta = s * rho * sqrt(sigmau^2 + sigmav^2) / sqrt(sigmax^2 + sigmay^2);
169. end
171. denom = (sigmaxu - s * sigmayv);
172. if denom == 0
173. theta = NaN;
174. else
175. theta = atand((sigmaxv - s * sigmayu) / denom);
176. end
177. end

vectorCorrelation.m at commit e8bd6f9, under BSD-3-Clause · at the source