s w.r.t length-scales</span>
  dldX = bsxfun(@times,l,bsxfun(@times,t,2.*exp(2*X(1:D,i)')).*t./2);

  M(i) = sum(lb)*c;                                            <span class="comment">% predicted mean</span>

  iK2beta = K2\beta(:,i);
  dMds(i,:,:) = c*t'*tlb/2-iR*M(i)/2;
  dMdX(i,D+2,i) = -c*sum(l.*(2*exp(2*X(D+2,i))*iK2beta));           <span class="comment">% OK</span>
  dMdX(i,D+1,i) = -dMdX(i,(i-1)*(D+2)+D+2);

  dVdX(:,i,D+2,i) = -((l.*(2*exp(2*X(D+2,i))*iK2beta))'*t*c)';
  dVdX(:,i,D+1,i) = -dVdX(:,