Fwd: qfunc/qfuncinv implementation in communications package

[Richardson, Anthony]

 

------ Original Message ------
From: Alec Teal
Date: 1/24/2013 1:48 PM
To: Richardson, Anthony;
Cc: address@hidden;
Subject: Re: qfunc/qfuncinv implementation in communications package

 

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

 

 

_____________________________________

I did not mean to imply that one was theoretically better than the other.  I meant that you practically run into machine epsilon problems with the current implementation.  qfunc is a wrapper around normcdf.  you will get better results wrapping it around erfc.  Try the line of test code in the original posting.  By the way, I believe octave uses the erfc implementation from C library.  Not sure how normcdf is implemented.

Tony