127 lines
5.5 KiB
Matlab
127 lines
5.5 KiB
Matlab
function fnnM = FalseNearestNeighborsSR(xV,tauV,mV,escape,theiler)
|
|
% fnnM = FalseNearestNeighbors(xV,tauV,mV,escape,theiler)
|
|
% FALSENEARESTNEIGHBORS computes the percentage of false nearest neighbors
|
|
% for a range of delays in 'tauV' and embedding dimensions in 'mV'.
|
|
% INPUT
|
|
% xV : Vector of the scalar time series
|
|
% tauV : A vector of the delay times.
|
|
% mV : A vector of the embedding dimension.
|
|
% escape : A factor of escaping from the neighborhood. Default=10.
|
|
% theiler : the Theiler window to exclude time correlated points in the
|
|
% search for neighboring points. Default=0.
|
|
% OUTPUT:
|
|
% fnnM : A matrix of size 'ntau' x 'nm', where 'ntau' is the number of
|
|
% given delays and 'nm' is the number of given embedding
|
|
% dimensions, containing the percentage of false nearest
|
|
% neighbors.
|
|
%========================================================================
|
|
% <FalseNearestNeighbors.m>, v 1.0 2010/02/11 22:09:14 Kugiumtzis & Tsimpiris
|
|
% This is part of the MATS-Toolkit http://eeganalysis.web.auth.gr/
|
|
|
|
%========================================================================
|
|
% Copyright (C) 2010 by Dimitris Kugiumtzis and Alkiviadis Tsimpiris
|
|
% <dkugiu@gen.auth.gr>
|
|
|
|
%========================================================================
|
|
% Version: 1.0
|
|
|
|
% LICENSE:
|
|
% This program is free software; you can redistribute it and/or modify
|
|
% it under the terms of the GNU General Public License as published by
|
|
% the Free Software Foundation; either version 3 of the License, or
|
|
% any later version.
|
|
%
|
|
% This program is distributed in the hope that it will be useful,
|
|
% but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
% GNU General Public License for more details.
|
|
%
|
|
% You should have received a copy of the GNU General Public License
|
|
% along with this program. If not, see http://www.gnu.org/licenses/>.
|
|
|
|
%=========================================================================
|
|
% Reference : D. Kugiumtzis and A. Tsimpiris, "Measures of Analysis of Time Series (MATS):
|
|
% A Matlab Toolkit for Computation of Multiple Measures on Time Series Data Bases",
|
|
% Journal of Statistical Software, in press, 2010
|
|
|
|
% Link : http://eeganalysis.web.auth.gr/
|
|
%=========================================================================
|
|
% Updates:
|
|
% Sietse Rispens, January 2012: use a different kdtree algorithm,
|
|
% searching k nearest neighbours, to improve performance
|
|
%
|
|
%=========================================================================
|
|
fthres = 0.1; % A factor of the data SD to be used as the maximum radius
|
|
% for searching for valid nearest neighbor.
|
|
propthres = fthres; % Limit for the proportion of valid points, i.e. points
|
|
% for which a nearest neighbor was found. If the proporion
|
|
% of valid points is beyond this limit, do not compute
|
|
% FNN.
|
|
|
|
n = length(xV);
|
|
if isempty(escape), escape=10; end
|
|
if isempty(theiler), theiler=0; end
|
|
% Rescale to [0,1] and add infinitesimal noise to have distinct samples
|
|
xmin = min(xV);
|
|
xmax = max(xV);
|
|
xV = (xV - xmin) / (xmax-xmin);
|
|
xV = AddNoise(xV,10^(-10));
|
|
ntau = length(tauV);
|
|
nm = length(mV);
|
|
fnnM = NaN*ones(ntau,nm);
|
|
for itau = 1:ntau
|
|
tau = tauV(itau);
|
|
for im=1:nm
|
|
m = mV(im);
|
|
nvec = n-m*tau; % to be able to add the component x(nvec+tau) for m+1
|
|
if nvec-theiler < 2
|
|
break;
|
|
end
|
|
xM = NaN*ones(nvec,m);
|
|
for i=1:m
|
|
xM(:,m-i+1) = xV(1+(i-1)*tau:nvec+(i-1)*tau);
|
|
end
|
|
% k-d-tree data structure of the training set for the given m
|
|
TreeRoot=kdtree_build(xM);
|
|
% For each target point, find the nearest neighbor, and check whether
|
|
% the distance increase over the escape distance by adding the next
|
|
% component for m+1.
|
|
idxV = NaN*ones(nvec,1);
|
|
distV = NaN*ones(nvec,1);
|
|
k0 = 2; % The initial number of nearest neighbors to look for
|
|
kmax = min(2*theiler + 2,nvec); % The maximum number of nearest neighbors to look for
|
|
for i=1:nvec
|
|
tarV = xM(i,:);
|
|
iV = [];
|
|
k=k0;
|
|
kmaxreached = 0;
|
|
while isempty(iV) && ~kmaxreached
|
|
[neiindV, neidisV] = kdtree_k_nearest_neighbors(TreeRoot,tarV,k);
|
|
% [neiM,neidisV,neiindV]=kdrangequery(TreeRoot,tarV,rthres*sqrt(m));
|
|
[oneidisV,oneiindV]=sort(neidisV);
|
|
neiindV = neiindV(oneiindV);
|
|
neidisV = neidisV(oneiindV);
|
|
iV = find(abs(neiindV(1)-neiindV(2:end))>theiler);
|
|
if k >= kmax
|
|
kmaxreached = 1;
|
|
elseif isempty(iV)
|
|
k = min(kmax,k*2);
|
|
end
|
|
end
|
|
if ~isempty(iV)
|
|
idxV(i) = neiindV(iV(1)+1);
|
|
distV(i) = neidisV(iV(1)+1);
|
|
end
|
|
end % for i
|
|
iV = find(~isnan(idxV));
|
|
nproper = length(iV);
|
|
% Compute fnn only if there is a sufficient number of target points
|
|
% having nearest neighbor (in R^m) within the threshold distance
|
|
if nproper>propthres*nvec
|
|
nnfactorV = 1+(xV(iV+m*tau)-xV(idxV(iV)+m*tau)).^2./distV(iV).^2;
|
|
fnnM(itau,im) = length(find(nnfactorV > escape^2))/nproper;
|
|
end
|
|
kdtree_delete(TreeRoot); % Free the pointer to k-d-tree
|
|
end
|
|
end
|