If you use regstats to estimate the coefficient standard errors, here's what you get using the Hald data:
load hald
s = regstats(heat,ingredients,'linear');
sqrt(diag(s.covb * (8/13)))' % 8/13 converts unbiased variance to mle
If you have another estimator that is the mle for some probability model, you could compute the second derivative of the log likelihood function and use that to estimate the standard errors. Here's an example using the normal distribution so it gives roughly the same answer as above:
loglik = @(b) sum(log(normpdf(heat,b(1)+ingredients*b(2:end),sqrt(8/13)*sqrt(s.mse))));
H = zeros(5);
b = s.beta;
db = 0.001;
for j=1:5,
ej = j==(1:5)';
H(j,j) = (loglik(b+ej*db)-2*loglik(b)+loglik(b-ej*db))/db^2;
for k=j+1:5
ek = k==(1:5)';
H(j,k) = ( loglik(b+ej*db+ek*db) + loglik(b) ...
- loglik(b+ej*db) - loglik(b+ek*db))/db^2;
H(k,j) = H(j,k);
end
end
sqrt(diag(inv(-H)))'
