Dear Tor,
Here are some comments re. your AR query.
In the paper Friston et al. (2000) "To smooth or not to smooth?"
equation A.3 shows the Yule-Walker relations. This shows
how the AR coefficients, 'a', can be computed from the sample
autocorrelation matrix using, for example, the Yule-Walker method:
apoly=aryule(timeseries,p);
a=-apoly(2:p+1);
where p is the model order.
Once you have computed 'a' you can use equation A.2 to compute
k:
A=tril(toeplitz([0,a(1:p)]));
k=inv(eye(p+1)-A);
You can then compute V:
V=k*k';
This V is
an AR-model-based-estimate of the true autocorrelation
matrix (showing the correlations at lags zero to p). It
should closely resemble the sample autocorrelation matrix
(that you can compute using xcorr - use [xc,lag]=xcorr(timeseries,'coeff'),
so that the sample correlation is 1 at lag zero and remember that
lag zero occurs half-way along the vector xc).
The higher the AR model order, the closer the resemblance.
Check out the attached matlab script.
Hope this helps.
All the best,
Will.
Tor Dessart Wager wrote:
> Dear SPM list,
>
> I'm trying to estimate the autocorrelation functions for my data using
> MATLAB and the methods described in "To Smooth or Not to Smooth?" (2000).
>
> I'm running into trouble defining the Va matrix. I'd like to use an
> AR(16) or AR(32) model, but I'm unclear about the notation in Appendix A.
>
> I assumed that rho(t) is the autocorrelation function, which is what I
> thought would be given by xcorr(timeseries). If that were so, then Va =
> toeplitz(xcorr(ts)), and there's no need for the AR model. So rho(t) must
> be something more complex. I did:
> if xc is xcorr(timeseries)...
> xc = xc/max(xc);
> xc = xc(200:end);
> V = toeplitz(xc); % A.1
>
> ~ ~ ~ ~
> Another option is to estimate the 'autoregression coefficients', a, using
> eq. A.3. In this equation, I take the inverse of the autocorrelation
> matrix, i think (?), and 'whiten' or deconvolve the autoregressive
> function, it looks like. I'm assuming that rho(i) are the elements of
> xcorr(timeseries), which gives me the autoregressive function.
>
> I can compute this, if M is toeplitz(xc), with:
> a = inv(M)*xc;
> a = inv(M)*[xc(2:end); 0]; % A.3
>
> ...but I'm not sure which. is rho(1) the first element of xc? If so,
> this should be 1, as the autocorrelation at time 0 = rho(1) = 1. Since
> the diagonals of the V matrix are 1, and the first off-diagonal is rho(1),
> I'm a little confused about exactly what to use. maybe rho(1) is the
> autocorrelation function at timeshift = 1? If so, I need to add an
> element to the xc vector to make it long enough to multiply; thus the
> second of these two options. One clarifying question: what is matrix M
> called in this case, if it's not the autocorrelation matrix V?
>
> If I use the elements of a to make a triangular matrix, as in equation
> A.2, and then compute the K and V matrices from that,
>
> A = zeros(size(ts,1),size(ts,1));
> for i = 1:size(ts,1)-1
> A = A + diag(ones(200-i,1)*a(i),i);
> end
> K = inv(1 - A);
> V = K * K';
>
> In this case, the diagonals of V do not equal one, and the matrix doesn't
> look much like what might be expected from an MRI autocorrelation
> function, so I must be doing something wrong.
>
> ~ ~ ~ ~
> In neither of these potential methods have I used the Yule-Walker method,
> which I thought was necessary. I can estimate the AR model using:
>
> coeff = aryule(timeseries,16) for an AR(16) model
>
> but this gives me coefficients for the polynomial, and I'm not sure how to
> convert those into an autocorrelation matrix V.
>
> In summary, how can I compute Va?
>
> I know this message may be very silly and confusing, but any help on this
> would be gratefully appreciated!
>
> Thanks,
> Tor
>
> _____________________________
> Tor Wager
> Department of Psychology
> University of Michigan
> Cognition and Perception Area
> 525 East University
> Ann Arbor, MI 48109-1109
>
> Office: 734-936-1295
> Home: 734-995-8975
> Email: [log in to unmask]
> _____________________________
--
William D. Penny
Wellcome Department of Cognitive Neurology
University College London
12 Queen Square
London WC1N 3BG
Tel: 020 7833 7478
FAX: 020 7813 1420
Email: [log in to unmask]
URL: http://www.fil.ion.ucl.ac.uk/~wpenny/
|