Re: qfunc/qfuncinv implementation in communications package

[Alec Teal]

On 24/01/13 19:26, Richardson, Anthony wrote:

The function qfunc() in the communications package is currently implemented as:

 

y = 1-normcdf(x);

 

This returns 0 for values of x greater than approx. 8 (where normcdf() is close to 1).  I believe a better implementation is:

 

y = erfc(x/sqrt(2))/2;

 

It depends how the integral is done really, I suspect this one is far more accurate as for stats most stuff is simply "good enough" it need not be that accurate. I'll look at the integral it should be using Sampson's (Simpsons?) rule of integral estimating at least and an extrapolation when the (2nd or 4th one of the two!) ratio of differences in the approximations comes to 0.5 or 0.25 (respectively) IIRC (long time since I've done Numerical Methods at A-level, and I hated it!)

Should improve accuracy at least.

A simple test to compare the two implementations is:

 

x = [0:3:18]; [qfunc(x); erfc(x/sqrt(2))/2]

 

Can the implementation of qfunc() be changed please?

 

The function qfuncinv() suffers from similar problems (it returns Inf for small arguments).  I had hoped just rewriting it in terms of erfcinv() would fix this problem too, but I found that erfcinv() is implemented in the specfun package as:

 

y= erfinv(1 – x)

 

Maclaurin Expansion (special case of taylor's series, even though taylor's came first....), and do it in order of magnitude (add up all the small terms first) I'll look into doing this at some point

which exhibits problems for small values of x (due to the subtraction from 1).  A better implementation of erfcinv is needed, but I don’t know what that would be.

 

Tony Richardson

 

_______________________________________________
Help-octave mailing list
address@hidden
https://mailman.cae.wisc.edu/listinfo/help-octave