Print

---

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