function tstatx = tstatapp(fun, beta, lambda);
% calculated tstats based on hessian and gradients
% inputs:  fun, beta,lambda
hh = myhessian(fun, beta, lambda);
dd = myjacobian(fun, beta, lambda);
robust = inv(hh) * (dd * dd') * inv(hh);
stdev = sqrt(diag(robust))';
tstatx = beta ./ stdev;