% compute inverse covariance matrix recursively

function S=Sigmainv(S,y,N);

S = N/(N-1)*(S-S*y*y'*S/(N-1+y'*S*y));

