 |
 |
| Hello SPM experts I am trying to implement the SPM approach for covariance component estimation. This is what I am doing: I have X:Design matrix Bhat: Estimate for parameter vector. y: Data V0: Starting value for the covariance matrix. My covariance matrix is Σ = o^2 * V , but I can solve for the ReML estimate of o^2 non-iterative: ohat^2 = 1/(n-p) * r' * V^-1 * r where n is the total number of scans, r is the vector of residuals and p is the length o Beta(number of columns in X) until convergence: invV = V^-1 Bhat = (X'*invV*X)^-1*X'*invV*y; r = y-X*Bhat; ohat = (r'*invV*r)/(n-p); P = invV-invV*X*(X'*invV*X)^-1*X'*invV; gi = -0.5*trace(P*Qi)+0.5*y'*P'*Qi*P*y;(gradient) Hi = -0.5*trace(P*Qi*P*Qi); (expectation of the Hessian) l = l + H^-1*g; V = l(1)*Q1+l(2)*Q2+...+l(k)*Qk; The Qis are covariance components that need to be especificy. Am I doing any thing wrong? In advance, thank you very much for your advice. -Jorge |