CircularRegression - Linear-circular regression.
There are three possible methods to estimate the regression parameters:
1) Use the method of Kempter et al. (2012). The slope maximizes the resultant length
of the residual distribution, the phase minimizes the circular variance of the
residual distribution.
2) Fit 'barber-pole' lines to the data. The model consists of three parallel lines
y = a.x + b + k*2pi (k in {-1,0,1}). The distance of a point to the model is
defined as the minimal distance to either of the three lines. The best-fit model
is computed using a least squared error approach.
3) Compute Thiel-Sen estimators after reorganizing angles around the regression line
obtained above (move by multiples of 2.pi).
Regression can be performed either on a list of (x,phi) pairs, or using a probability
distribution of phi for each value of x.
USAGE
[beta,R2,p] = CircularRegression(x,angles,<options>)
x values for linear independent variable
angles values for circular dependent variable (in radians)
<options> optional list of property-value pairs (see table below)
[beta,R2,p] = CircularRegression(x,angles,p,<options>)
x N linear variable bins
angles M circular variable bins
p MxN matrix where the nth column is the angular probability
distribution for the nth value of x (angular bins appear
in ascending order from top to bottom)
<options> optional list of property-value pairs (see table below)
=========================================================================
Properties Values
-------------------------------------------------------------------------
'slope' search start value for slope (default = 0)
'method' regression method: 'k' for Kempter (default), 'pole' for
barber pole, or 'ts' for Theil-Sen (but 'ts' cannot be
computed for probability distribution matrices)
=========================================================================
OUTPUT
beta regression slope and intercept
R2 coefficient of determination
p p-value associated with H0: slope = 0 (Kempter only)
SEE
See also CircularVariance, CircularConfidenceIntervals, Concentration,
ConcentrationTest.0001 function [beta,R2,p] = CircularRegression(x,angles,p,varargin) 0002 0003 %CircularRegression - Linear-circular regression. 0004 % 0005 % There are three possible methods to estimate the regression parameters: 0006 % 0007 % 1) Use the method of Kempter et al. (2012). The slope maximizes the resultant length 0008 % of the residual distribution, the phase minimizes the circular variance of the 0009 % residual distribution. 0010 % 0011 % 2) Fit 'barber-pole' lines to the data. The model consists of three parallel lines 0012 % y = a.x + b + k*2pi (k in {-1,0,1}). The distance of a point to the model is 0013 % defined as the minimal distance to either of the three lines. The best-fit model 0014 % is computed using a least squared error approach. 0015 % 0016 % 3) Compute Thiel-Sen estimators after reorganizing angles around the regression line 0017 % obtained above (move by multiples of 2.pi). 0018 % 0019 % Regression can be performed either on a list of (x,phi) pairs, or using a probability 0020 % distribution of phi for each value of x. 0021 % 0022 % USAGE 0023 % 0024 % [beta,R2,p] = CircularRegression(x,angles,<options>) 0025 % 0026 % x values for linear independent variable 0027 % angles values for circular dependent variable (in radians) 0028 % <options> optional list of property-value pairs (see table below) 0029 % 0030 % [beta,R2,p] = CircularRegression(x,angles,p,<options>) 0031 % 0032 % x N linear variable bins 0033 % angles M circular variable bins 0034 % p MxN matrix where the nth column is the angular probability 0035 % distribution for the nth value of x (angular bins appear 0036 % in ascending order from top to bottom) 0037 % <options> optional list of property-value pairs (see table below) 0038 % 0039 % ========================================================================= 0040 % Properties Values 0041 % ------------------------------------------------------------------------- 0042 % 'slope' search start value for slope (default = 0) 0043 % 'method' regression method: 'k' for Kempter (default), 'pole' for 0044 % barber pole, or 'ts' for Theil-Sen (but 'ts' cannot be 0045 % computed for probability distribution matrices) 0046 % ========================================================================= 0047 % 0048 % OUTPUT 0049 % 0050 % beta regression slope and intercept 0051 % R2 coefficient of determination 0052 % p p-value associated with H0: slope = 0 (Kempter only) 0053 % 0054 % SEE 0055 % 0056 % See also CircularVariance, CircularConfidenceIntervals, Concentration, 0057 % ConcentrationTest. 0058 0059 % Copyright (C) 2012 by Michaël Zugaro 0060 % 0061 % This program is free software; you can redistribute it and/or modify 0062 % it under the terms of the GNU General Public License as published by 0063 % the Free Software Foundation; either version 3 of the License, or 0064 % (at your option) any later version. 0065 0066 % Default values 0067 method = 'k'; 0068 slope = 0; 0069 0070 % Check parameters 0071 if nargin < 2, 0072 error('Incorrect number of parameters (type ''help <a href="matlab:help CircularRegression">CircularRegression</a>'' for details).'); 0073 end 0074 0075 if ~isdvector(x), 0076 error('Incorrect x values (type ''help <a href="matlab:help CircularRegression">CircularRegression</a>'' for details).'); 0077 end 0078 if ~isdvector(angles), 0079 error('Incorrect angles (type ''help <a href="matlab:help CircularRegression">CircularRegression</a>'' for details).'); 0080 end 0081 0082 if nargin >= 3 0083 if isdmatrix(p), 0084 % Two lists of bins and a probability matrix 0085 if size(p,1) ~= length(angles) || size(p,2) ~= length(x), 0086 error('Incoherent parameter sizes (type ''help <a href="matlab:help CircularRegression">CircularRegression</a>'' for details).'); 0087 end 0088 else 0089 % Two lists of samples (x and angles) 0090 varargin = {p,varargin{:}}; 0091 p = []; 0092 end 0093 else 0094 p = []; 0095 end 0096 0097 % Parse parameter list 0098 for i = 1:2:length(varargin), 0099 if ~ischar(varargin{i}), 0100 error(['Parameter ' num2str(i+2) ' is not a property (type ''help <a href="matlab:help CircularRegression">CircularRegression</a>'' for details).']); 0101 end 0102 switch(lower(varargin{i})), 0103 case 'method', 0104 method = varargin{i+1}; 0105 if ~isastring(method,'pole','k','ts'), 0106 error('Incorrect value for property ''method'' (type ''help <a href="matlab:help CircularRegression">CircularRegression</a>'' for details).'); 0107 end 0108 case 'slope', 0109 slope = varargin{i+1}; 0110 if ~isdscalar(slope), 0111 error('Incorrect value for property ''slope'' (type ''help <a href="matlab:help CircularRegression">CircularRegression</a>'' for details).'); 0112 end 0113 otherwise, 0114 parameters = varargin{i:end}; 0115 if ~isa(parameters,'cell'), parameters = {parameters}; end 0116 break; 0117 end 0118 end 0119 0120 if strcmp(method,'ts') && ~isempty(p), 0121 error('Thiel-Sen regression cannot be computed for probability distribution matrices (type ''help <a href="matlab:help CircularRegression">CircularRegression</a>'' for details).'); 0122 end 0123 0124 isradians(angles); 0125 if ~isempty(p), 0126 % Replicate x values (one copy for each angular bin) 0127 [m,n] = size(p); 0128 x = reshape(x,1,length(x)); 0129 x = repmat(x,m,1); 0130 x = x(:); 0131 % Angular bins (one copy for each x bin) 0132 phi = angles(:); 0133 phi = repmat(phi,n,1); 0134 prob = p(:); 0135 else 0136 x = x(:); 0137 phi = angles(:); 0138 prob = ones(size(phi)); 0139 prob = prob/length(prob); 0140 end 0141 0142 beta = [nan nan]; 0143 R2 = nan; 0144 p = nan; 0145 0146 if isempty(angles) || isempty(x), return; end 0147 0148 % Store data in global variables (necessary for minimization) 0149 global CircularRegression_x CircularRegression_phi CircularRegression_prob; 0150 CircularRegression_x = x; 0151 CircularRegression_phi = phi; 0152 CircularRegression_prob = prob; 0153 0154 % ------------------------------------------------------------------------------------------------ 0155 % Compute barber-pole regression 0156 % (this is a first step for Thiel-Sen) 0157 0158 if strcmp(method,'pole') || strcmp(method,'ts'), 0159 % Minimize residual error (least squared error) 0160 [beta,RSE] = fminsearch(@ResidualSquaredError,[slope 0]); 0161 0162 % Compute coefficient of determination 0163 TSE = norm(prob.*(phi-CircularMean(phi)))^2; 0164 R2 = 1-RSE/TSE; 0165 % p-value not implemented 0166 end 0167 0168 % ------------------------------------------------------------------------------------------------ 0169 % Compute 'Kempter' regression 0170 0171 if strcmp(method,'k'), 0172 % Maximize resultant length 0173 [a,R] = fminsearch(@ResultantLength,slope); 0174 0175 % Estimate phase shift 0176 C = sum(prob.*cos(phi-a*x)); 0177 S = sum(prob.*sin(phi-a*x)); 0178 phi0 = atan2(S,C); 0179 0180 beta = [a phi0]; 0181 0182 % Compute coefficient of determination 0183 theta = mod(abs(a)*x,2*pi); 0184 phi_bar = atan2(sum(prob.*sin(phi)),sum(prob.*cos(phi))); 0185 theta_bar = atan2(sum(prob.*sin(theta)),sum(prob.*cos(theta))); 0186 N = sum(prob.*sin(phi-phi_bar).*sin(theta-theta_bar)); 0187 D = sqrt(sum(prob.*sin(phi-phi_bar).^2)*sum(prob.*sin(theta-theta_bar).^2)); 0188 rho = N/D; 0189 R2 = rho^2; 0190 0191 % Compute p-value 0192 n = length(phi); 0193 lambda_02 = 1/n*sum(prob.*sin(theta-theta_bar).^2); 0194 lambda_20 = 1/n*sum(prob.*sin(phi-phi_bar).^2); 0195 lambda_22 = 1/n*sum(prob.*sin(phi-phi_bar).^2.*prob.*sin(theta-theta_bar).^2); 0196 z = rho*sqrt(n*lambda_20*lambda_02/lambda_22); 0197 p = 1-erf(abs(z)/sqrt(2)); 0198 0199 end 0200 0201 % ------------------------------------------------------------------------------------------------ 0202 % Compute Theil-Sen estimator 0203 0204 if strcmp(method,'ts'), 0205 % Move phi values by multiples of 2.pi to get them as close as possible to the regression line 0206 d = phi-beta(1)*x+beta(2); 0207 rd = sign(d).*floor(abs(d/(2*pi))); 0208 phi = phi+rd*2*pi; 0209 0210 % Keep at most 500 random samples 0211 n = length(x); 0212 if n == 1, 0213 beta = [nan nan]; 0214 R2 = nan; 0215 return; 0216 end 0217 pp = randperm(n); 0218 n = min([n 500]); 0219 x0 = x(pp(1:n)); 0220 phi0 = phi(pp(1:n)); 0221 0222 % Compute pairwise slopes 0223 slopes = nan(nchoosek(n,2),1); 0224 k = 1; 0225 for i = 1:n-1, 0226 for j = (i+1):n, 0227 slopes(k) = (phi0(i)-phi0(j))/(x0(i)-x0(j)); 0228 k = k + 1; 0229 end 0230 end 0231 0232 % Take median 0233 beta(1) = median(slopes); 0234 beta(2) = median(phi-beta(1)*x); 0235 % Compute coefficient of determination 0236 RSE = sum((phi-(beta(1)*x+beta(2))).^2); 0237 R2 = 1-RSE/TSE; 0238 % p-value not implemented 0239 end 0240 0241 % -------------------------------------------------------------------------------------------------------------- 0242 0243 function RSE = ResidualSquaredError(beta) 0244 0245 % Retrieve data 0246 global CircularRegression_x CircularRegression_phi; 0247 x = CircularRegression_x; 0248 phi = CircularRegression_phi; 0249 0250 % Model parameters: slope and intercept 0251 a = beta(1); 0252 b = beta(2); 0253 0254 % Three lines 0255 y1 = a*x+b; 0256 y2 = y1+2*pi; 0257 y3 = y1-2*pi; 0258 0259 % Squared error 0260 d = min(prob.*[(phi-y1).^2 (phi-y2).^2 (phi-y3).^2],[],2); 0261 RSE = sum(d); 0262 0263 % -------------------------------------------------------------------------------------------------------------- 0264 0265 function R = ResultantLength(beta) 0266 0267 % Retrieve data 0268 global CircularRegression_x CircularRegression_phi CircularRegression_prob; 0269 x = CircularRegression_x; 0270 phi = CircularRegression_phi; 0271 prob = CircularRegression_prob; 0272 0273 % Model parameters: slope 0274 a = beta(1); 0275 0276 n = length(x); 0277 C = 1/n*sum(prob.*cos(phi-a*x)); 0278 S = 1/n*sum(prob.*sin(phi-a*x)); 0279 R = sqrt(C.^2+S.^2); 0280 0281 % We must actually return -R because Matlab can only minimize (but not maximize...) 0282 R = -R;