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.

Mathematica 9 has a new CorrelationFunction. Sadly the joy does not last long, as I can't get it to work with vectors. I would like to make a velocity correlation function. That is, given a list of n vectors I would like to calculate:

velocity correlation function

For example, this is almost instantaneous:

testdata = RandomVariate[BinormalDistribution[{10, 20}, 0.4], {10000}];
AbsoluteTiming[
x = CorrelationFunction[testdata[[All, 1]], {0, 1000}];
y = CorrelationFunction[testdata[[All, 2]], {0, 1000}];
]
ListPlot[{x, y}, PlotRange -> {0, 1}]

(the test data here is not the best as it is by definition uncorrelated)

Here I would like the autocorrelation of 2D vectors evenly spaced in time. But this does not work:

r = CorrelationFunction[testdata, {0, 100}];

How can I make an autocorrelation for vectors? Or do I have to go the manual way?

ClearAll[AbsoluteAutocorrelationT, AutocorrelationT]

AbsoluteAutocorrelationT[l_, 0] := Total[(#.#) & /@ l];
AbsoluteAutocorrelationT[l_, n_] := 
  Total@Table[l[[i]].l[[i + n]], {i, Length[l] - n}];
AutocorrelationT[l_, range_] := Block[{t0, ti},
   t0 = AbsoluteAutocorrelationT[l, 0];
   ti = AbsoluteAutocorrelationT[l, #] & /@ (Range @@ range);
   (ti/t0)];

AbsoluteTiming[r = AutocorrelationT[testdata, {1, 1000}];]

Which is 200 times slower. I know I can compile this (but than I have to worry about overflows in the sum, etc.)

share|improve this question
    
Easy to guess why this isn't supported directly, but it does look like something of an oversight not to at least mention it in the documentation. One obvious way around it would be to interpolate and resample the data. –  Oleksandr R. Feb 25 '13 at 17:45
    
@OleksandrR. and what is the reason (must be a bit slow today:) –  Ajasja Feb 25 '13 at 17:48
1  
Correlation is probably being done using an FFT, which requires evenly spaced data. The good news is an FFT is fast, so you can happily upsample your data significantly, especially if you're able to arrive at a power of two length without absurdly excessive padding. –  Oleksandr R. Feb 25 '13 at 17:51
2  
No, the 1-d case is assumed to be evenly spaced to start with. Your 2-d input is assumed to be of the form {{t1, x1}, ..., {tn, xn}}. If I misunderstood your question and you want the n-dimensional autocorrelation for values sampled on a regular grid, rather than that for 1-d but unevenly spaced data, you can get it yourself using ListCorrelate. –  Oleksandr R. Feb 25 '13 at 17:59
1  
Given your code, isn't this what you are trying to do? absCorr[lists_, n_] := Length@lists AbsoluteCorrelationFunction[#, n] & /@ Transpose@lists // Total –  Rojo Feb 25 '13 at 19:56

2 Answers 2

When each $v_i = (v_{i,1}, v_{i,2}, \ldots, v_{i,k})$ is a $k$-vector we may interchange the summations to obtain

$$\sum_{i=1}^{n-t} v_i\cdot v_{i+t} = \sum_{j=1}^k\sum_{i=1}^{n-t} v_{i,j} v_{i+t,j},$$

whence the desired correlation is the sum of the correlations of the components. Thus:

kernel = PDF[NormalDistribution[0, 1]] /@ Range[-5, 5, .005]; (* Induces autocorrelation *)
testdata = ListConvolve[kernel, #] & /@
  Transpose[RandomVariate[BinormalDistribution[{10, 20}, 0.4], {10^4}]];
c = Sum[ListCorrelate[v, v, {1, 1}], {v, testdata}]; // AbsoluteTiming

$\{0.0030002, \text{Null}\}$

ListLinePlot[c[[;; Length[c]/2]]/c[[1]]]

Line plot of ACF

The timings scale as expected ($k=2$ shown):

timing vs. size

share|improve this answer

The following is (only) six times faster:

l = RandomVariate[BinormalDistribution[{10, 20}, 0.4], {10000}];

c[l_, 0] := ListCorrelate[l, l, {1, -1}, {}, Dot, Plus, 1];
c[l_, s_List] := Module[{m = c[l, 0][[1]]}, 
   ListCorrelate[l[[;;Length@l - #]], l[[1+#;;]], {1, -1}, {}, Dot, Plus, 1] / m & /@ Range @@ s];

ListLinePlot@Table[ AbsoluteTiming[AutocorrelationT[l, {1, n}];][[1]]/
                    AbsoluteTiming[c[l, {1, n}];][[1]], {n, 30}]

Mathematica graphics

share|improve this answer

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.