7
$\begingroup$

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.)

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Oleksandr R. Feb 25 '13 at 17:45
  • $\begingroup$ @OleksandrR. and what is the reason (must be a bit slow today:) $\endgroup$ – Ajasja Feb 25 '13 at 17:48
  • 1
    $\begingroup$ 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. $\endgroup$ – Oleksandr R. Feb 25 '13 at 17:51
  • 2
    $\begingroup$ 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. $\endgroup$ – Oleksandr R. Feb 25 '13 at 17:59
  • 1
    $\begingroup$ Given your code, isn't this what you are trying to do? absCorr[lists_, n_] := Length@lists AbsoluteCorrelationFunction[#, n] & /@ Transpose@lists // Total $\endgroup$ – Rojo Feb 25 '13 at 19:56
2
$\begingroup$

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

$\endgroup$
1
$\begingroup$

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

$\endgroup$
1
$\begingroup$

TL/DR:

temporaltestData = TemporalData[{testdata},Automatic];
r = CorrelationFunction[temporaltestData, {0, 1000}]["Values"];

Recap of Question

In your question, you have provided the code

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}]

and you found that this code finds the autocorrelations in the first and second components of your time series just fine. However, you also wanted to find the lagged cross-correlations between the first and second components.

Attempted solution

You attempted to do this with

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

but you found that this doesn't work. I got an error that said

CorrelationFunction::tmpln: Uniformly spaced sampling expected in {{-3.38558,-24.5187},{-5.09168,-0.21692},{7.14122,14.2209},{-9.21681,19.2432},{8.15215,14.9583},{0.911419,-8.53165},{-19.1743,-22.1},{10.885,14.1329},{5.9896,3.49815},{9.51289,32.5084},<<31>>,{-1.15955,25.1969},{11.6961,13.6815},{-13.5907,-14.7231},{15.2179,17.0268},{4.61082,11.2426},{4.04496,-2.09691},{-0.278494,17.8067},{-4.02344,10.879},{19.5057,8.53984},<<9950>>}.

What I believe is happening is that mathematica is interpreting the first component as a time. You can try to get around this by prepending the times:

timedTestData = Prepend[testdata[[#]], #] & /@ Range@Length@testdata
r = CorrelationFunction[timedTestData, {0, 100}]

but now I get the error

CorrelationFunction::bdlag: The lag specification {0,100} should be a symbol, an integer with magnitude less than the length of the data, or a range specification indicating such integers.

which doesn't make sense.

Real solution

However, there is a way to find the correlation, as you can see by going to the help for CorrelationFunction, and then going to Examples->Scope->Empirical Estimates, and looking at the last example in this section. The example is

proc = ARProcess[{{{.3, .1}, {.6, .3}}}, {{1, .3}, {.3, .6}}];
data = RandomFunction[proc, {100}];
cov = CorrelationFunction[data, {-6, 6}];

They have no problem taking the correlation of the multidimensional data set. This is because their data is a TemporalData object. So we must make one of our own. This can be done with

temporaltestData = TemporalData[{testdata},Automatic];

The brackets around testdata are necessary because otherwise TemporalData will treat testdata as 10000 time series of length two instead of one 2D time series of length 10000. The second argument, Automatic, just assigns consecutive integer times to each data point. Once the data is packaged in the TemporalData object, CorrelationFunction has a much easier time knowing what to do with it:

r = CorrelationFunction[temporaltestdata, {0, 1000}]["Values"]

Here, CorrelationFunction returns a TimeSeries object, which is a special kind of TemporalData object. If we just want the list of lagged autocovariance matrices, we can access them by applying our time series object to "Values". Now we have our lagged autocovariances. Notice for example:

Norm@Chop[r[[;; , 1, 1]] - x]

gives 0

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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