return statistical scalar to user

[jbraw]

With these two lines in paresto, you are doing in the main stat calculation. Fstat = finv(alpha, n_est, n_data-n_est);

theta_conf = sqrt(2*n_est/(n_data-n_est)*Fstat*r.f*diag(inv_hess));

I think it would be nice to pass back this scalar as well.

b = 2n_est/(n_data-n_est)Fstat*r.f

Then the user doesn't have to redo this calculation outside of paresto when the user wants to do something else with a confidence interval. For example, when doing a 2-d problem, I often want to plot the ellipse:

deltheta' H del\theta \leq b

you provide H already (est.d2_dtheta2). It would be nice to have b as well. Not sure what to call it.

est.b, est.level, est.prob_level

What do you think?

Jim Rawlings

[jbraw]

Also, it might be nice to provide the user with ellipse.m. That's a small function that's handy for plotting a 2-d ellipse.

octave:1> help ellipse [x, y, major, minor, bbox] = ellipse (amat, level, n, shift)

Given a 2x2 matrix, generate ellipse data for plotting. The arguments N and SHIFT are optional. If N is an empty matrix, a default value of 100 is used.

[jaeandersson]

Good idea. I'm not sure what to call it, or what it's physical interpretation is. Plotting the ellipse makes a lot of sense. There could be a helper function for plotting the uncertainty ellipse for a given set of parameters.

[jbraw]

Given the ellipse calling list, I would lean toward calling the new parameter est.level.

If we make ellipse.m available, I don't think you need a helper function. Just call

ellipse(est.d2f_dtheta2, est.level, n, est.theta(ind_est))

The user would have to know that he's estimated two and only two parameters for this call to work. ellipse.m has its own error checking and reporting so we don't have to do much in paresto.m.

Jim