+ l.*bdX(:,i,d);
    dtdX = inp*(-bsxfun(@times,iR(:,d),2.*exp(2*X(d,i))*iR(d,:)));
    dlbt = lb'*dtdX + dlb'*t;
    dVdX(:,i,d,i) = (dlbt'*c + cdX(d)*(lb'*t)');
  end % d
  
  dMdi(i,:,:) = c*(tlbdi - tlb);
  dMdt(i,:,i) = c*liK';
  dMdX(i,1:D,i) = cdX.*sum(beta(:,i).*l) + c.*dslb;
  v = bsxfun(@rdivide,inp,exp(X(1:D,i)'));
  k(:,i) = 2*X(D+1,i)-sum(v.*v,2)/2;
  V(:,i) = t'*lb*c;                                   % input-output covariance
  
  for d = 1:D
    dVds(d,i,:,:) = c*bsxfun(@times,t,t(:,d))'*tlb