function [f,g,h] = crit(b,w)

k = length(b);
m = moment(b);
mbar = mean(m)';

f = mbar'/w*mbar;

if nargout > 1
    Q = dmoment(b);
    Qn = mean(Q,3);
    g = 2*Qn'/w*mbar;
end

if nargout > 2
    R = d2moment(b);
    Rn = mean(R,3);
    h = 2*(Qn'/w*Qn + kron(mbar'/w,eye(k))*Rn);
end