Something wrong with correlation dimension(Grassberger Proccacia) code

Hi all, I downloaded this code for calculating correlation dimension from pudn.com, but it isnt giving desired output for any particular input. Kindly run it and help fix the errors . Im even having a lot of trouble understanding this code as it didnt come with any proper comments, kindly help in that perspective too.

It would be much better if someone could write a basic Grassberger Procaccia correlation code for N number of points distributed in m dimensional space. I dont want it for any time series, rather simply for N number of points distributed in M dimensional space. Please see these -- http://en.wikipedia.org/wiki/Correlation_dimension

http://www.physics.unlv.edu/~thanki/thesis/node14.html

F=130*rand(1,8192);
F=F';
m=6;%Dimension
c=300;%The number of each dimension vector
y=zeros(c,m);
t=18;%Time delay
p=10;%Vector spacing

for ii=1:c
    for qq=1:m
        y(ii,qq)=F((ii-1)*p+(qq-1)*t+1);
    end
end
%re-building the above  phase space

n=0;
%Dot pitch
for i=1:c-1
    for j=i+1:c
        n=n+1;
        ee(n)=dist(y(i,:),y(j,:)');
    end
end

mm=0;
pl=0;%Correlation dimension

for eps=min(ee):0.9:min(ee)+(max(ee)-min(ee))/9 %eps:Given the very small number of points

      mm=mm+1;%How many points
      nn=0;
      for kk=1:n
          if ee(kk)<=eps
              nn=nn+1;
          end
      end
     pl(mm)=log(nn/(c*c))/log(eps);

end
 plot(pl);

[Walter Roberson]

Please show the error message and which line it is occurring at, and please show the size() of all variables referred to on that line.

[Vinita]

Walter actually the code isnt printing out correlation dimension(I mean it should at least print a value for corelation dimension).

And if I play with values of m (the Dimension) ,then also theres no visible effect in the plot(pl).

Also I want to modify this code.I dont want it for any time series, rather simply for N number of points distributed in M dimensional space. I have been looking for a basic Grassberger Procaccia correlation code for N number of points distributed in m dimensional space.

[Vinita]

Please help..i still havent been able to fix it

[Walter Roberson]

There is nothing in that code that attempts to display the correlation dimension numerically. That is the design, rather than an error. In order to "Kindly run it and help fix the errors" you need to tell us which MATLAB error message you are encountering, or you must show us where the calculation does not agree with what other evidence shows the correct calculation to be.

Changing to code to do something other than what it was designed to do is _not_ fixing an error.

If you want something displayed, display p1 instead of plotting it.

[Vinita]

Walter Im really not that good in matlab. So being a biologist,i have been struggling a lot since a week to write something as complicated as correlation dimension.

My problem is that :

I have N(50,000) points, and I have all the pairwise distance information between these points calculated and stored in a separate file. Now I want to find correlation dimension for these points. And Im totally unaware as to how to proceed for my case.

Kindly help

Answer 1

[Vinita]

Guys if theres some info from my side that I have been unable to provide,atleast let me know of that. !!!!!!!!

[Walter Roberson]

A lot of people are asleep now. I should be too (it is 4 AM my time.)

Answer 2

[Vishnu Sreekumar]

Hi Vinita,

If you have the distances stored in a separate file, all you need to do is plot the cumulative counts on a log-log plot and calculate the slope. The slope gives you the estimate of the correlation dimension plot. Here's how you do it: Assume G is your distance matrix (store it as an upper or lower triangle matrix since it is a symmetric matrix). What I'm doing below isn't the most efficient way to do it, this is just a quick pointer:

[a1,a2] = size(G);
N = a1*a2;
DistVect = reshape(G',N,1); 
%just making it a single vector, this isn't necessary if you know how to code well, which I don't. ;)

DistVect = DistVect(find(DistVect>10^-8)); 
%You don't want to count 0 distances

[g1,xout1] = hist(log(DistVect),1000); %choose how many ever bins you want

%Calculate the number of elements in each bin
n_elements = histc(log(DistVect),xout1);

%Calculate the cumulative sum of these elements using cumsum
c_elements = cumsum(n_elements);

A1 = xout1;
B1 = log(c_elements);
plot(A1,B1);
xlabel('log(r)')
ylabel('log(C(r))');
title('Correlation Dimension')

%The slope of this plot will give you the correlation dimension estimate. The data need not be single scaled of course (depending on the structure of your data, which is what you want to find here). So you might need to estimate the dimensionalities at different length scales in that case (just calculate the slopes in the regions of interest).

You say that you have 50,000 points. I'm not sure how exactly you stored them. You need a (50,000*49,999/2)*8 bytes (or close to 9GB) of memory to do that. Matlab will not handle 9GB of contiguous memory as far as I understand if you are on a PC with the usual 4 or 8GB of RAM (I'm not an expert either, others can help verify this). So unless you get much more memory, I doubt you will be able to work with this large matrix. You should try to break this computation down in such a way that you don't have to read in and work with the whole 9GB at the same time.

Thanks, Vishnu

[Walter Roberson]

The 64 bit versions of MATLAB can handle 9 GB of contiguous memory. If the system involved did not have that much physical memory then it would have to have been set up to use swap space, which is usually like running on a computer 100 times slower.

Answer 3

[Vishnu Sreekumar]

Btw the code you attached in your original post seems to be doing the Takens' delay embedding procedure to reconstruct the dynamical properties of the phase space (of which the correlation dimension is one). Did you want to do Takens' delay embedding? If not, you don't need that code.

[Vinita]

Hi Vishnu.

Im using this command to load all files in present directory

files = dir('*.out');

for i=1:length(files)

eval(['load ' files(i).name ]);

end

Each of these 50 files(suppose named file1,file2....) are 1x10000.

Now I want to modify each of these 10000 values (and do so for each of the 50 files)as : -

for j=1:50

for i=1:10000

file(j)(i)= file(j)(i)/(1-file(j)(i));

end

end.

But this code is failing.I get this error and dont know how to fix it :-

??? Error: ()-indexing must appear last in an index

expression.

[Walter Roberson]

Do not use eval() if it can be avoided.

load(files(i).name)

For the other difficulties, see http://matlab.wikia.com/wiki/FAQ#How_can_I_create_variables_A1.2C_A2.2C....2CA10_in_a_loop.3F

[Vinita]

Thanks walter. But it isnt helping. If i simply run this code for a single file

for i=1:10000

file1(i)= file1(i)/(1-file1(i));

end

Its working. But when i introduce a second loop so as to cover all the 50 files,it fails.

Please help

[Vinita]

Have searched the forums but I really cant understand it. Please help.I have been trying a lot but it isnt fixing.

[Vinita]

WORKING >> for i=1:10000

X10k1(i)= X10k1(i)/(1-X10k1(i));

end

NOT WORKING >>for j=1:50

for i=1:10000

p=X10k(j);

p(i)= p(i)/(1-p(i));

end

end

??? Undefined function or method 'X10k' for input arguments

of type 'double'.

[Vinita]

Somebody ??