Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I want to compute the Autocorrelation Function (ACF) of a data table with $10^6$ entries. I know that there is a built-in function in Mathematica for that, but because I do not know how exactly it is defined, I tend to avoid it. Furthermore, I need the autocorrelation for Fluorescent Correlation Spectroscopy (FCS), which is specifically normalized:

$$G(\Delta t) = \frac{\displaystyle\sum_{i=0}^{M-m}I(i\tau)I(i\tau + m\tau)}{\langle I\rangle^2(M-m)}$$

where $\Delta t=m\tau$, $\tau$ is the step size, $M$ is the total numbers of steps.

Can you help me to define it, so that the computation finishes in less than a day?

This is how I generate my data with Mathematica

    n = 50.;
    radius = 6.*10.^-8.;
    k = 1.38*10.^-23.;
    T = 293.;
    Eta = 1.*10.^-3.;
    d = (k*T)/(6.*\[Pi]*Eta*radius);
    Deltat = 500000.;
    time = 10.*10.^-6.;
    Taud = (Omegar^2./(4.*d));
    Omegar = 200.*10.^-9.;
    Omegaz = 5.*Omegar;
    Io = 1;
    boundary = 5.*10.^-6.;
    initial = boundary/2.;
    step = Sqrt[2.*d*time];
    RandomWalk[x_] := 
      Accumulate[
       Join[{RandomReal[{-initial, initial}, 3]}, 
         RandomVariate[NormalDistribution[0, step], {x, 3}]]];
    p = Table[Mod[RandomWalk[Deltat], boundary, -initial], {i, n}];
    particleintensity = 
       Io Exp[(p^2).{-2/Omegar^2, -2/Omegar^2, -2/Omegaz^2}];
    int = Total[particleintensity];
    ListPlot[int, AxesOrigin -> {0., 0.}, AxesLabel -> {"t", "I"}, 
        PlotRange -> Full, PlotRangeClipping -> False, PlotStyle -> Red]

I tried with the definition from Wikipedia:

$$G(\tau)=\frac{\langle\delta I(t)\delta I(t+\tau)\rangle}{\langle I(t) \rangle^2}=\frac{\langle I(t)I(t+\tau) \rangle}{\langle I(t)\rangle^2}-1$$

which I implemented like this

    acff = Table[Mean[int.Take[PadLeft[int, t], Deltat]], {t, 1, Deltat}]
    acf = acff/Mean[int] - 1

but this still takes ridiculously long time to compute.

share|improve this question
    
Here's an autocorrelation function that you might be able to work with: demonstrations.wolfram.com/… –  bill s Apr 23 '13 at 20:41
    
How many time steps are in your data? I'm asking so that I can check if the code I've written is fast enough for you. –  RunnyKine Apr 23 '13 at 21:12
    
@bills the demonstration code was quite complex. I took the definition of acf, but this has no chance of completing, before my server runs out of free memory –  phidelio Apr 23 '13 at 21:18
3  
What's wrong with using ListCorrelate? –  Daniel Lichtblau Apr 23 '13 at 22:28
1  
@RunnyKine I'd need a better spec of the problem. Like what is I()? And is tau an integer? –  Daniel Lichtblau Apr 23 '13 at 22:56

2 Answers 2

up vote 2 down vote accepted

OK, I am stupid. ListCorrelate is what I want, indeed. I just have to use the correct parameters.

acf = norm2 ListCorrelate[int, int, {1, 1}, 0]*norm1 -1;

where norm1 and norm2 give me the normation I need and are defined as

norm1 = Table[1/(Deltat + 1 - m), {m, 0, Deltat}];
norm2 = 1/Mean[int]^2;

Thanks to Bill S, RunnyKine and Daniel Lichtblau for the crucial help.

share|improve this answer

As Daniel suggests, use ListCorrelate, It works very quickly:

r = RandomReal[{-1, 1}, 10^6];
ListCorrelate[r, r] // Timing

takes less than 0.02 seconds. You will want to read the documentation carefully to understand how to scale it as you desire.

share|improve this answer
    
To answer all questions about ListCorrelate. Yes it is very fast and represents the ACF curve very well. The problem is that the amplitude at $G(0)$ should be equal to $\frac{1}{n}$, where n is the number of particles. I see exactly the opposite behaviour. Furthermore, ListCorrelate[{a, b, c}, {a, b, c}, {1, 1}, 0] gives me {a^2 + b^2 + c^2, a b + b c, a c}. This is not what I would expect. I couldn't find in the documentation anything useful or comprehensable. More info at link –  phidelio Apr 24 '13 at 7:14
    
I forgot something. Yes $tau$ is an integer, but only in my code. Generally, it is measured in microseconds. –  phidelio Apr 24 '13 at 7:29

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.