)'));
  
  for j=1:i
    R = s*diag(exp(-2*X(1:D,i))+exp(-2*X(1:D,j)))+eye(D);
    t = 1/sqrt(det(R));
    ij = bsxfun(@rdivide,inp,exp(2*X(1:D,j)'));
    L = exp(bsxfun(@plus,k(:,i),k(:,j)')+maha(ii,-ij,R\s/2));
    S(i,j) = t*beta(:,i)'*L*beta(:,j); S(j,i) = S(i,j);
  end
  
  S(i,i) = S(i,i) + 1e-6;          % add small jitter for numerical reasons
  
end

% 4) Centralize moments
S = S - M*M';

##### SOURCE END #####
--></body></html>