an class="keyword">for</span> 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);
  <span class="keyword">end</span>

  S(i,i) = S(i,i) + 1e-6;          <span class="comment">% add small jitter for numerical reasons</span>

<span class="keyword">end</span>

<span class="comment">% 4) Centralize moments</span>