Estimation of diffusion coefficient and directional bias

error_ellipse.m

function [theta ax] =error_ellipse(varargin)
% OUTPUT: theta - rotation angle of the major longer axis of the ellipse
%                 with respect to the horizontal plane in rads

% error_ellipse - plot an error ellipse, or ellipsoid, defining confidence
% region
% error_ellipse(C22) - Given a 2x2 covariance matrix, plot the
% associated error ellipse, at the origin. It returns a graphics handle
% of the ellipse that was drawn.
%
% error_ellipse(C33) - Given a 3x3 covariance matrix, plot the
% associated error ellipsoid, at the origin, as well as its projections
% onto the three axes. Returns a vector of 4 graphics handles, for the
% three ellipses (in the X-Y, Y-Z, and Z-X planes, respectively) and for
% the ellipsoid.
%
% error_ellipse(C,MU) - Plot the ellipse, or ellipsoid, centered at MU,
% a vector whose length should match that of C (which is 2x2 or 3x3).
%
% error_ellipse(...,'Property1',Value1,'Name2',Value2,...) sets the
% values of specified properties, including:
% 'C' - Alternate method of specifying the covariance matrix
% 'mu' - Alternate method of specifying the ellipse (-oid) center
% 'conf' - A value betwen 0 and 1 specifying the confidence interval.
% the default is 0.5 which is the 50% error ellipse.
% 'scale' - Allow the plot the be scaled to difference units.
% 'style' - A plotting style used to format ellipses.
% 'clip' - specifies a clipping radius. Portions of the ellipse, -oid,
% outside the radius will not be shown.
%
% NOTES: C must be positive definite for this function to work properly.
ax = [];
FIGURE = 0;
default_properties = struct(...
    'C', [], ... % The covaraince matrix (required)
    'mu', [], ... % Center of ellipse (optional)
    'conf', 0.5, ... % Percent confidence/100
    'scale', 1, ... % Scale factor, e.g. 1e-3 to plot m as km
    'style', '', ... % Plot style
    'clip', inf, ... % Clipping radius
    'figure', FIGURE); %figure on
if length(varargin) >= 1 && isnumeric(varargin{1})
    default_properties.C = varargin{1};
    varargin(1) = [];
end

if length(varargin) >= 1 && isnumeric(varargin{1})
    default_properties.mu = varargin{1};
    varargin(1) = [];
end

if length(varargin) >= 1 && isnumeric(varargin{1})
    default_properties.conf = varargin{1};
    varargin(1) = [];
end

if length(varargin) >= 1 && isnumeric(varargin{1})
    default_properties.scale = varargin{1};
    varargin(1) = [];
end

if length(varargin) >= 1 && isnumeric(varargin{1})
    FIGUREON = varargin{1};
    varargin(1) = [];
end

if length(varargin) >= 1 && ~ischar(varargin{1})
    error('Invalid parameter/value pair arguments.')
end

prop = getopt(default_properties, varargin{:});
C = prop.C;

if isempty(prop.mu)
    mu = zeros(length(C),1);
else
    mu = prop.mu;
end

conf = prop.conf;
scale = prop.scale;

if conf <= 0 || conf >= 1
    error('conf parameter must be in range 0 to 1, exclusive')
end

[r,c] = size(C);
if r ~= c || (r ~= 2 && r ~= 3)
    error(['Don''t know what to do with ',num2str(r),'x',num2str(c),' matrix'])
end

x0=mu(1);
y0=mu(2);
% Compute quantile for the desired percentile
k = sqrt(qchisq(conf,r)); % r is the number of dimensions (degrees of freedom)
if FIGURE
hold_state = get(gca,'nextplot');
end
if r==3 && c==3
    z0=mu(3);
    
    % Make the matrix has positive eigenvalues - else it's not a valid
    % covariance matrix!
    if any(eig(C) <=0)
        %error('The covariance matrix must be positive definite (it has non-positive eigenvalues)')
    end
    
    % C is 3x3; extract the 2x2 matricies, and plot the associated error
    % ellipses. They are drawn in space, around the ellipsoid; it may be
    % preferable to draw them on the axes.
    Cxy = C(1:2,1:2);
    Cyz = C(2:3,2:3);
    Czx = C([3 1],[3 1]);
    
    [x,y,z] = getpoints(Cxy,prop.clip);
    if FIGURE
    h1=plot3(x0+k*x,y0+k*y,z0+k*z,prop.style);
    hold on
    end
    [y,z,x] = getpoints(Cyz,prop.clip);
    if FIGURE
    h2=plot3(x0+k*x,y0+k*y,z0+k*z,prop.style);hold on
    end
    [z,x,y] = getpoints(Czx,prop.clip);
    if FIGURE
    h3=plot3(x0+k*x,y0+k*y,z0+k*z,prop.style);hold on
    end
    [eigvec,eigval] = eig(C);
    
    [X,Y,Z] = ellipsoid(0,0,0,1,1,1);
    XYZ = [X(:),Y(:),Z(:)]*sqrt(eigval)*eigvec';
    % eigen values corresponds to the two axis of the ellipse
    ax(1,1) = eigval(1,1);
    ax(1,2) = eigval(2,2);
    
    
    X(:) = scale*(k*XYZ(:,1)+x0);
    Y(:) = scale*(k*XYZ(:,2)+y0);
    Z(:) = scale*(k*XYZ(:,3)+z0);
    if FIGURE
        h4=surf(X,Y,Z);
        colormap gray
        alpha(0.3)
        camlight
        if nargout
            h=[h1 h2 h3 h4];
        end
    end
elseif r==2 && c==2
    % Make the matrix has positive eigenvalues - else it's not a valid
    % covariance matrix!
    if any(eig(C) <=0)
        %error('The covariance matrix must be positive definite (it has non-positive eigenvalues)')
    end
    % eigen values corresponds to the two axis of the ellipse
    [eigvec,eigval] = eig(C);
    ax(1,1) = eigval(1,1);
    ax(1,2) = eigval(2,2);
    [x,y,z, theta] = getpoints(C,prop.clip);
    if FIGURE
    h1=plot(scale*(x0+k*x),scale*(y0+k*y),prop.style);
    set(h1,'zdata',z+1)
    if nargout
       h=h1;
    end
    end
else
    error('C (covaraince matrix) must be specified as a 2x2 or 3x3 matrix)')
end
if FIGURE
axis equal
set(gca,'nextplot',hold_state);
end
%---------------------------------------------------------------
% getpoints - Generate x and y points that define an ellipse, given a 2x2
% covariance matrix, C. z, if requested, is all zeros with same shape as
% x and y.

function [x,y,z, theta] = getpoints(C,clipping_radius)

n=100; % Number of points around ellipse
p=0:pi/n:2*pi; % angles around a circle
[eigvec,eigval] = eig(C); % Compute eigen-stuff
theta =(cart2pol(eigvec(2, 1),eigvec(2,2)));
xy = [cos(p'),sin(p')] * sqrt(eigval) * eigvec'; % Transformation
x = xy(:,1);
y = xy(:,2);
z = zeros(size(x));

% Clip data to a bounding radius
if nargin >= 2
    r = sqrt(sum(xy.^2,2)); % Euclidian distance (distance from center)
    x(r > clipping_radius) = nan;
    y(r > clipping_radius) = nan;
    z(r > clipping_radius) = nan;
end

%---------------------------------------------------------------

function x=qchisq(P,n)
% QCHISQ(P,N) - quantile of the chi-square distribution.
if nargin<2
    n=1;
end

s0 = P==0;
s1 = P==1;
s = P>0 & P<1;
x = 0.5*ones(size(P));
x(s0) = -inf;
x(s1) = inf;
x(~(s0|s1|s))=nan;

for ii=1:14
    dx = -(pchisq(x(s),n)-P(s))./dchisq(x(s),n);
    x(s) = x(s)+dx;
    if all(abs(dx) < 1e-6)
        break;
    end
end

%---------------------------------------------------------------

function F=pchisq(x,n)
% PCHISQ(X,N) - Probability function of the chi-square distribution.
if nargin<2
    n=1;
end
F=zeros(size(x));

if rem(n,2) == 0
    s = x>0;
    k = 0;
    for jj = 0:n/2-1;
        k = k + (x(s)/2).^jj/factorial(jj);
    end
    F(s) = 1-exp(-x(s)/2).*k;
else
    for ii=1:numel(x)
        if x(ii) > 0
            F(ii) = quadl(@dchisq,0,x(ii),1e-6,0,n);
        else
            F(ii) = 0;
        end
    end
end

%---------------------------------------------------------------
function f=dchisq(x,n)
% DCHISQ(X,N) - Density function of the chi-square distribution.
if nargin<2
    n=1;
end
f=zeros(size(x));
s = x>=0;
f(s) = x(s).^(n/2-1).*exp(-x(s)/2)./(2^(n/2)*gamma(n/2));

%---------------------------------------------------------------
function properties = getopt(properties,varargin)
%GETOPT - Process paired optional arguments as
%'prop1',val1,'prop2',val2,...
%
% getopt(properties,varargin) returns a modified properties structure,
% given an initial properties structure, and a list of paired arguments.
% Each argumnet pair should be of the form property_name,val where
% property_name is the name of one of the field in properties, and val is
% the value to be assigned to that structure field.
%
% No validation of the values is performed.
%
% EXAMPLE:
% properties =
% struct('zoom',1.0,'aspect',1.0,'gamma',1.0,'file',[],'bg',[]);
% properties = getopt(properties,'aspect',0.76,'file','mydata.dat')
% would return:
% properties =
% zoom: 1
% aspect: 0.7600
% gamma: 1
% file: 'mydata.dat'
% bg: []
%
% Typical usage in a function:
% properties = getopt(properties,varargin{)
% Process the properties (optional input arguments)
prop_names = fieldnames(properties);
TargetField = [];
for ii=1:length(varargin)
    arg = varargin{ii};
    if isempty(TargetField)
        if ~ischar(arg)
            error('Propery names must be character strings');
        end
        f = find(strcmp(prop_names, arg));
        if isempty(f)
            error('%s ',['invalid property ''',arg,'''; must be one of:'],prop_names{:});
        end
        TargetField = arg;
    else
        % properties.(TargetField) = arg; % Ver 6.5 and later only
        properties = setfield(properties, TargetField, arg); % Ver 6.1 friendly
        TargetField = '';
    end
end
if ~isempty(TargetField)
    error('Property names and values must be specified in pairs.');
end

estimateDiffCoefficients.m

function [D1, D2, D3, bias] = estimateDiffCoefficients (x, y, samplingFreq)

% x, y: horizontal and vertical traces of drift samples
%   D1: diffusion coefficient obtained from the slope of the regression
%       line between <d^2> and deltaT.
%   D2: diffusion coefficient obtained from the slope of the regression
%       line between Gaussian area and deltaT.
%   D3: diffusion coefficient obtained from the slope of the regression
%       line after removing bias term using PCA.
% BiasMeanX, BiasMeanY: Directional bias
%
% example:
% [x, y] = generate2DRandomWalk (41, 1000, 100, 1000, 60, 10, .02);
% [D1, D2, D3, bias] = estimateDiffCoefficients (x, y, 1000);
% % plot estimated bias direction
% plot(atan(bias.biasMeanY./bias.biasMeanX).*180/pi) %


dimensions = 2;

[~, nT] = size(x);
tm = (0:nT-1)./samplingFreq;   % create a time vector for plotting
Dx = cell(nT,1);
Dy = cell(nT,1);

for dt = 1:nT-1
    dx = x(:,1+dt:end) - x(:,1:end-dt);
    dy = y(:,1+dt:end) - y(:,1:end-dt);          % pool displacements
    
    Dx{dt+1} = [Dx{dt+1} dx(:)'];
    Dy{dt+1} = [Dy{dt+1} dy(:)'];
end

% displacements at zero time
dispSquared = zeros(1, nT);
area = zeros(1, nT);
area_unbiased = zeros(1, nT);
area_unbiased(1) = 0;

for dt = 2:nT
    % method-1
    dispSquared(dt) = mean(Dx{dt}.^2 + Dy{dt}.^2);
    
    % method-2 finding std along the principle axes
    n = length(Dx{dt});
    mat = [Dx{dt}; Dy{dt}];
    sigmas = eig(mat*mat')./n;
    area(dt) = sum(sigmas);
    
    % method-3 (removing the mean (bias))
    bias.biasMeanX(dt-1) = mean(Dx{dt});
    bias.biasMeanY(dt-1) = mean(Dy{dt});
    bias.biasStdX(dt-1) = std(Dx{dt});
    bias.biasStdY(dt-1) = std(Dy{dt});
    
    mat = [Dx{dt} - bias.biasMeanX(dt-1); ...
        Dy{dt} - bias.biasMeanY(dt-1)];
    
    sigmas = eig(mat*mat')./n;
    area_unbiased(dt) = sum(sigmas);
end

% compute diffusion coefficent
ind = 1:round(nT / 2);

slp1 = regress(dispSquared(ind)', tm(ind)');
D1 = slp1/(2*dimensions);

% DEBUGGING: plot regression
%figure;
%plot(tm, dispSquared, 'b-', tm(ind), tm(ind)*slp1, 'r:');
%legend('Actual', sprintf('Slope %f', slp1));
%title('Code by Murat; uses intervals up to T / 2');

slp2 = regress(area(ind)', tm(ind)');
D2 = slp2/(2*dimensions);

slp3 = regress(area_unbiased(ind)', tm(ind)');
D3 = slp3/(2*dimensions);

run_comparison.m

set(0, 'DefaultAxesLineWidth', 2);
set(0, 'DefaultLineLineWidth', 2);
set(0, 'DefaultAxesFontSize', 16);

numb_iter = 500;
len_iter = 350;

%% Scenario 1: no bias
diff_co = 40;
bias_spd = 0;
values = zeros(numb_iter, 7);
for j = 1:numb_iter
    [x, y] = simulate_random_walk(diff_co, len_iter);
    
    % save trace
    if 1 == j
        plot(x, y);
        title(sprintf('Random walk; diffusion coefficient = %d; bias = %d', diff_co, bias_spd));
        print(gcf, 'trace1.png', '-dpng', '-r300');
    end
    
    % Murat's code
    [d1, d2, d3, bias] = estimateDiffCoefficients(x', y', 1000);
    
    % Martina's code
    a = struct('x', x, 'y', y);
    [Area, Dsq, Time, AreaCoef, DCoef, angle, axisRatio, xyCorrelation, Mean, std_th] = SimpleCalculate2DArea(x, y);
    
    mt = numel(bias.biasMeanX) + 1;
    bias_norm = mean(sqrt(((bias.biasMeanX ./ (2:mt) * 1000) .^ 2) + ((bias.biasMeanY ./ (2:mt) * 1000) .^ 2)));
    
    values(j, :) = [d1 d2 d3 bias_norm AreaCoef DCoef std_th];
    
    %fprintf('Diffusion coefficients are D1=%4.2f, D2=%4.2f, D3=%4.2f, D4=%4.2f, D5=%4.2f\n', d1, d2, d3, DCoef, AreaCoef)
end

% display mean values
disp(mean(values));

% plot
hist(values(:, 1), 20);
xlabel('Diffusion coefficient');
title(sprintf('Standard estimation; actual = %d; bias = %d', diff_co, bias_spd));
print(gcf, 'est1img1.png', '-dpng', '-r300');

hist(values(:, 2), 20);
xlabel('Diffusion coefficient');
title(sprintf('Elliptical estimation; actual = %d; bias = %d', diff_co, bias_spd));
print(gcf, 'est1img2.png', '-dpng', '-r300');

hist(values(:, 3), 20);
xlabel('Diffusion coefficient');
title(sprintf('Unbiased elliptical estimation; actual = %d; bias = %d', diff_co, bias_spd));
print(gcf, 'est1img3.png', '-dpng', '-r300');

hist(values(:, 4), 20);
xlabel('Bias');
title(sprintf('Removed bias; actual = %d; bias = %d', diff_co, bias_spd));
print(gcf, 'est1img4.png', '-dpng', '-r300');

hist(values(:, 7), 20);
xlabel('Diffusion coefficient');
title(sprintf('Standard estimation (Martina); actual = %d; bias = %d', diff_co, bias_spd));
print(gcf, 'est1img5.png', '-dpng', '-r300');

%% Scenario 2: with bias
diff_co = 40;
bias_spd = 150;
values = zeros(numb_iter, 7);
for j = 1:numb_iter
    [x, y] = simulate_random_walk(diff_co, len_iter, 45, bias_spd, bias_spd / 4);
    
    % save trace
    if 1 == j
        plot(x, y);
        title(sprintf('Random walk; diffusion coefficient = %d; bias = %d', diff_co, bias_spd));
        print(gcf, 'trace2.png', '-dpng', '-r300');
    end
    
    % Murat's code
    [d1, d2, d3, bias] = estimateDiffCoefficients(x', y', 1000);
    
    % Martina's code
    a = struct('x', x, 'y', y);
    [Area, Dsq, Time, AreaCoef, DCoef, angle, axisRatio, xyCorrelation, Mean, std_th] = SimpleCalculate2DArea(x, y);
    
    mt = numel(bias.biasMeanX) + 1;
    bias_norm = mean(sqrt(((bias.biasMeanX ./ (2:mt) * 1000) .^ 2) + ((bias.biasMeanY ./ (2:mt) * 1000) .^ 2)));
    
    values(j, :) = [d1 d2 d3 bias_norm AreaCoef DCoef std_th];
    
    %fprintf('Diffusion coefficients are D1=%4.2f, D2=%4.2f, D3=%4.2f, D4=%4.2f, D5=%4.2f\n', d1, d2, d3, DCoef, AreaCoef)
end

% display mean values
disp(mean(values));

% plot
hist(values(:, 1), 20);
xlabel('Diffusion coefficient');
title(sprintf('Standard estimation; actual = %d; bias = %d', diff_co, bias_spd));
print(gcf, 'est2img1.png', '-dpng', '-r300');

hist(values(:, 2), 20);
xlabel('Diffusion coefficient');
title(sprintf('Elliptical estimation; actual = %d; bias = %d', diff_co, bias_spd));
print(gcf, 'est2img2.png', '-dpng', '-r300');

hist(values(:, 3), 20);
xlabel('Diffusion coefficient');
title(sprintf('Unbiased elliptical estimation; actual = %d; bias = %d', diff_co, bias_spd));
print(gcf, 'est2img3.png', '-dpng', '-r300');

hist(values(:, 4), 20);
xlabel('Bias');
title(sprintf('Removed bias; actual = %d; bias = %d', diff_co, bias_spd));
print(gcf, 'est2img4.png', '-dpng', '-r300');

hist(values(:, 7), 20);
xlabel('Diffusion coefficient');
title(sprintf('Standard estimation (Martina); actual = %d; bias = %d', diff_co, bias_spd));
print(gcf, 'est2img5.png', '-dpng', '-r300');

%% Scenario 3: no bias
diff_co = 80;
bias_spd = 0;
values = zeros(numb_iter, 7);
for j = 1:numb_iter
    [x, y] = simulate_random_walk(diff_co, len_iter);
    
    % save trace
    if 1 == j
        plot(x, y);
        title(sprintf('Random walk; diffusion coefficient = %d; bias = %d', diff_co, bias_spd));
        print(gcf, 'trace3.png', '-dpng', '-r300');
    end
    
    % Murat's code
    [d1, d2, d3, bias] = estimateDiffCoefficients(x', y', 1000);
    
    % Martina's code
    a = struct('x', x, 'y', y);
    [Area, Dsq, Time, AreaCoef, DCoef, angle, axisRatio, xyCorrelation, Mean, std_th] = SimpleCalculate2DArea(x, y);
    
    mt = numel(bias.biasMeanX) + 1;
    bias_norm = mean(sqrt(((bias.biasMeanX ./ (2:mt) * 1000) .^ 2) + ((bias.biasMeanY ./ (2:mt) * 1000) .^ 2)));
    
    values(j, :) = [d1 d2 d3 bias_norm AreaCoef DCoef std_th];
    
    %fprintf('Diffusion coefficients are D1=%4.2f, D2=%4.2f, D3=%4.2f, D4=%4.2f, D5=%4.2f\n', d1, d2, d3, DCoef, AreaCoef)
end

% display mean values
disp(mean(values));

% plot
hist(values(:, 1), 20);
xlabel('Diffusion coefficient');
title(sprintf('Standard estimation; actual = %d; bias = %d', diff_co, bias_spd));
print(gcf, 'est3img1.png', '-dpng', '-r300');

hist(values(:, 2), 20);
xlabel('Diffusion coefficient');
title(sprintf('Elliptical estimation; actual = %d; bias = %d', diff_co, bias_spd));
print(gcf, 'est3img2.png', '-dpng', '-r300');

hist(values(:, 3), 20);
xlabel('Diffusion coefficient');
title(sprintf('Unbiased elliptical estimation; actual = %d; bias = %d', diff_co, bias_spd));
print(gcf, 'est3img3.png', '-dpng', '-r300');

hist(values(:, 4), 20);
xlabel('Bias');
title(sprintf('Removed bias; actual = %d; bias = %d', diff_co, bias_spd));
print(gcf, 'est3img4.png', '-dpng', '-r300');

hist(values(:, 7), 20);
xlabel('Diffusion coefficient');
title(sprintf('Standard estimation (Martina); actual = %d; bias = %d', diff_co, bias_spd));
print(gcf, 'est3img5.png', '-dpng', '-r300');

%% Scenario 4: with bias
diff_co = 80;
bias_spd = 150;
values = zeros(numb_iter, 7);
for j = 1:numb_iter
    [x, y] = simulate_random_walk(diff_co, len_iter, 45, bias_spd, bias_spd / 4);
    
    % save trace
    if 1 == j
        plot(x, y);
        title(sprintf('Random walk; diffusion coefficient = %d; bias = %d', diff_co, bias_spd));
        print(gcf, 'trace4.png', '-dpng', '-r300');
    end
    
    % Murat's code
    [d1, d2, d3, bias] = estimateDiffCoefficients(x', y', 1000);
    
    % Martina's code
    a = struct('x', x, 'y', y);
    [Area, Dsq, Time, AreaCoef, DCoef, angle, axisRatio, xyCorrelation, Mean, std_th] = SimpleCalculate2DArea(x, y);
    
    mt = numel(bias.biasMeanX) + 1;
    bias_norm = mean(sqrt(((bias.biasMeanX ./ (2:mt) * 1000) .^ 2) + ((bias.biasMeanY ./ (2:mt) * 1000) .^ 2)));
    
    values(j, :) = [d1 d2 d3 bias_norm AreaCoef DCoef std_th];
    
    %fprintf('Diffusion coefficients are D1=%4.2f, D2=%4.2f, D3=%4.2f, D4=%4.2f, D5=%4.2f\n', d1, d2, d3, DCoef, AreaCoef)
end

% display mean values
disp(mean(values));

% plot
hist(values(:, 1), 20);
xlabel('Diffusion coefficient');
title(sprintf('Standard estimation; actual = %d; bias = %d', diff_co, bias_spd));
print(gcf, 'est4img1.png', '-dpng', '-r300');

hist(values(:, 2), 20);
xlabel('Diffusion coefficient');
title(sprintf('Elliptical estimation; actual = %d; bias = %d', diff_co, bias_spd));
print(gcf, 'est4img2.png', '-dpng', '-r300');

hist(values(:, 3), 20);
xlabel('Diffusion coefficient');
title(sprintf('Unbiased elliptical estimation; actual = %d; bias = %d', diff_co, bias_spd));
print(gcf, 'est4img3.png', '-dpng', '-r300');

hist(values(:, 4), 20);
xlabel('Bias');
title(sprintf('Removed bias; actual = %d; bias = %d', diff_co, bias_spd));
print(gcf, 'est4img4.png', '-dpng', '-r300');

hist(values(:, 7), 20);
xlabel('Diffusion coefficient');
title(sprintf('Standard estimation (Martina); actual = %d; bias = %d', diff_co, bias_spd));
print(gcf, 'est4img5.png', '-dpng', '-r300');

SimpleCalculate2DArea.m

function [Area, Dsq, Time, AreaCoef, DCoef, angle, axisRatio, xyCorrelation, Mean, std_th] = SimpleCalculate2DArea(fix_x, fix_y, varargin)
%SIMPLECALCULATE2DAREA Summary of this function goes here
%   This function calculates the are covered by the 2D displacement distribution
%   at each delta t.
%
% function [Area] = DiffusionCoefficient2D (Fix,varargin)
%
% INPUT:
% Fix   : matrix with the eye movements data. This matrix should be
%         organized in the the following way:
%         Fix{sbj} is an array containing each fixation as a structure.
%         A given drift's horizontal and vertical angles (in arcmin)
%         can be accessed via Fix{sbj}(fx).x and Fix{sbj}(fx).x. where fx
%         is the fixation index
%
% Properties: Optional properties
%             'MaxTime'       : Maximum temporal interval
%                               (default: 256 ms)
%             'SpanLimit'     : Maximum span of an event to be considered
%                               events with a higher span will be discarded
%                               (default: 60 arcmin)
%             'FigKey'        : Set FigKey at one to generate a figure
%                               (default: 0)
%             'TimeInterval'  : Temporal segment of the event considered
%                               if it is set at 1 all the event segment is selected
%                               otherwise a 2x1 vector will indicate the temporal
%                               start and end of the segment
%                              (default: [1])
% OUTPUT:
% Area           : Area (armin^2) covered by the displacements distribution at each delta t.
% Dsq            : Displacement squared (arcmin) at each delta t
% Time           : The temporal bins over which the area has been calculated.
% NFix           : Number of events considered for each delta t.
% RegCoef        : Regression coefficient for the fit of the area vs. time.
% theta          : degrees of rotations of the major axis of the ellips for the
%                 final value of delta t
% ax             : major and minor axis at the final value of delta t (the ratio
%                 between the two axis gives a measure of the distribution symmetry.
% xyCorrelation  :
%
% HISTORY:
% This function is based on the DiffusionCoefficients_DPI.m function by
% Murat, which is based on Xutao diffusion analysis function
%
% 2013 @APLAB

% set default parameters
NT = 256;
TIMEINTERVAL = 1; %ms
FSAMPLING = 1000;
DIMENSIONS = 2;

% Interpret the user parameters
k = 1;
Properties = varargin;
while k <= length(Properties) && ischar(Properties{k})
    switch (Properties{k})
        case 'MaxTime'
            NT = Properties{k + 1};
            Properties(k:k+1) = [];
        case 'TimeInterval'
            TIMEINTERVAL = Properties{k + 1};
            Properties(k:k+1) = [];
        otherwise
            k = k + 1;
    end
end

if length(TIMEINTERVAL)==2
    if TIMEINTERVAL(2)<NT
        warning(sprintf('Maximum Time %.0f ms is larger than the Time Interval considered [%.0f-%.0f]',...
            NT, TIMEINTERVAL(1),TIMEINTERVAL(2)));
        NT = TIMEINTERVAL(2)-TIMEINTERVAL(1)-1;
        warning(sprintf('Maximum time has been reset at %.0f' , NT));
    end
end

% time space
Time = linspace(1,NT-1,NT-1)./FSAMPLING;

% process eye movements
[Area, Dsq, angle, axisRatio, xyCorrelation, Mean, std_th] = probEyeMov(fix_x, fix_y, Time, TIMEINTERVAL);

% regress the area and find the regression coefficients
RegCoef = regress(Area', Time');
RegCoef1 = regress(Dsq', Time');

% DEBUGGING: plot regression
%figure;
%plot(Time, Dsq, 'b-', Time, Time*RegCoef1(1), 'r:');
%legend('Actual', sprintf('Slope %f', RegCoef1(1)));
%title('Code by Martina; intervals up to MaxTime (default 255)');


% rate of increase per sec
AreaCoef = (RegCoef(1)/(2*DIMENSIONS));
DCoef = (RegCoef1(1)/(2*DIMENSIONS));

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
function  [Area, Dsq, angle, axisRatio, xyCorrelation, Mean, std_th] = probEyeMov(fix_x, fix_y, time, time_interval)

Mean = [];

% loop through all delta t
for dt = 1:length(time);
    d_x= [];
    d_y = [];
    dd2 = [];
    
    % get movement
    if length(time_interval)==2
        hor = fix_x(time_interval(1):time_interval(2));
        ver = fix_y(time_interval(1):time_interval(2));
    else
        hor = fix_x(time_interval(1):end);
        ver = fix_y(time_interval(1):end);
    end

    % pick all the possible couples of values based on nT and the
    % maximum time interval
    enD = length(hor)-dt;
    d_x = [ d_x (hor(1+dt:enD+dt)-hor(1:enD)) ];
    d_y = [ d_y (ver(1+dt:enD+dt)-ver(1:enD)) ];
    dd2 = [ dd2 (hor(1+dt:enD+dt)-hor(1:enD)).^2 + (ver(1+dt:enD+dt)-ver(1:enD)).^2 ];
    
    % calculate the area of the 2D distribution
    c = corrcoef(d_x',d_y');
    % area at .68 without multiplication by pi (to make it comparable with
    % the d)
    Area(dt) = 2*std(d_x)*std(d_y)*sqrt(1-c(2,1)^2);
    Dsq(dt) = mean(dd2);%mean(d_x.^2+d_y.^2);
    Mean(dt,1:2) = [mean(d_x) mean(d_y)];
    C = cov(d_x,d_y);
    [theta, ax] = error_ellipse(C, [mean(d_x); mean(d_y)],'conf', .95, 'style', 'w');
    axisRatio(dt) = max(ax)/min(ax);
    angle(dt) = rad2deg(theta);
    % corrleation coefficient of the distributions on the two axis (if the
    % distribution is not symmeterical the correlation value is higher)
    xyCorrelation(dt) = c(2,1);
end
th = cart2pol(d_x, d_y);
std_th = std(th);

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

end

end

simulate_random_walk.m

function [x, y] = simulate_random_walk(diffusion_coefficient, samples, bias_dir, bias_spd, bias_spd_std, freq)
%SIMULATE_RANDOM_WALK Generate random walk motion
%   Diffusion coefficient: determines the standard deviation per second
%   Samples: the number of samples to simulate
%   Bias_dir: optional, the direction of bias in degrees
%   Bias_spd: optional, the speed of bias movement (unit/s)
%   Bias_spd_std: optional, the standard deviation of the bias movement
%   (unit/s)
%   Freq: the sampling frequency, default to 1 kHz

% set frequency
if nargin < 6
    freq = 1000;
end

tau = 1 / freq;
dimensions = 2;

% generate deltas
d = normrnd(0, sqrt(diffusion_coefficient * dimensions * tau), samples, 2);

% generate bias
if nargin >= 5
    bias_mov = normrnd(bias_spd * tau, sqrt(bias_spd_std) * tau, samples, 2) * diag([cosd(bias_dir) sind(bias_dir)]);
    d = d + bias_mov;
end

x = cumsum(d(:, 1));
y = cumsum(d(:, 2));

end