L = t*iL;
  c = exp(2*X(D+1,i))/sqrt(det(B));

  M(i) = sum(lb)*c;                                            <span class="comment">% predicted mean</span>
  V(:,i) = tL'*lb*c;                   <span class="comment">% inv(s) times input-output covariance</span>
  k(:,i) = 2*X(D+1,i)-sum(in.*in,2)/2;
<span class="keyword">end</span>

<span class="comment">% 3) Compute predictive covariance, non-central moments</span>
<span class="keyword">for</span> i=1:E
  ii = bsxfun(@rdivide,inp,exp(2*X(1:D,i)'));

  <sp