From: Christopher on
Peter, that helps a lot. But now I have a new problem. The '.01' is actually an element in an array, so the real equation reads

sigma = [.01 .05 .1 .2 .25];
for i = 1:5
berder(2,i) = 1-erfc((-.5/sqrt(sigma(i)/3))/sqrt(2))/2
end

where there are a few values of sigma and this equation is being used in a loop. only the lowest value of sigma, .01, is generating this error. But now that I am using the vpa function over the sigma array, it returns a symbolic expression that can't be stored an type double array.



Peter Perkins <Peter.Perkins(a)MathRemoveThisWorks.com> wrote in message <hs12fr$nk$1(a)fred.mathworks.com>...
> On 5/6/2010 11:19 PM, Christopher wrote:
> > I'm doing a statistics project and when using the erfc, I get a value my
> > professor tells me is wrong.
> > "1-erfc((-.5/sqrt(.01/3))/sqrt(2))/2 = 0"
> >
> > I know that the real value is close to 0, but isn't actually 0. It's
> > something on the order of 2*10^-15 or so.
> > How can I get MatLab to give the the correct value?
>
> Christopher, this is kind of a funny question. Most people start with
> erf, or with the Gaussian cumulative distribution function, and have a
> similar computational problem (where the correct result is very small),
> and then discover erfc. But OK, let's assume you really started with erfc.
>
> First, you need to think about double precision arithmetic. Consider
> this simpler situation:
>
> >> x = 1e-18
> x =
> 1.000000000000000e-018
> >> y = 1 - (1-x)
> y =
> 0
>
> Is this just a display thing? No, y really _is_ exactly zero:
>
> >> isequal(y,0)
> ans =
> 1
>
> The question you first need to figure out is, why is y zero? Hint: the
> outer subtraction is not at fault.
>
> Next: as a student, you likely have the Symbolic Toolbox, so do this in
> variable precision arithmetic (which, for your purposes, is cheating):
>
> >> vpa('erfc((-.5/sqrt(.01/3))/sqrt(2))/2',32)
> ans =
> 0.99999999999999999764642970492981
>
> So the difference between your erfc value and 1, the thing corresponding
> to (1-x) above, is on the order of 10-18. Now you need to figure out
> how to compute 1 - erfc(...)/2 without ending up subtracting 1 from 1.
> It will help to draw a picture of what erfc is (the upper tail of an
> integral) and what 1-erfc/2 is.
>
> You likely also have the Statistics Toolbox. Once you figure out the
> above, you might draw a picture of what the Gaussian CDF represents, and
> then try to figure out these lines in the NORMCDF function:
>
> z = (x-mu) ./ sigma;
> p = 0.5 * erfc(-z ./ sqrt(2));
>
> _That's_ the direction most people end up with erfc from.
>
> Hope this helps.
From: Christopher on
Never mind the previous question. I figured it out. Instead of

>>1-erfc(-z/2^.5)

like I was doing, I used

>>erfc(z/2^.5)

Thanks for the help!
From: Peter Perkins on
On 5/7/2010 2:20 PM, Christopher wrote:
> Peter, that helps a lot. But now I have a new problem. The '.01' is
> actually an element in an array, so the real equation reads
> sigma = [.01 .05 .1 .2 .25];
> for i = 1:5
> berder(2,i) = 1-erfc((-.5/sqrt(sigma(i)/3))/sqrt(2))/2
> end
>
> where there are a few values of sigma and this equation is being used in
> a loop. only the lowest value of sigma, .01, is generating this error.
> But now that I am using the vpa function over the sigma array, it
> returns a symbolic expression that can't be stored an type double array.

Christopher, if you're using vpa, and still subtracting erfc from 1,
then you've missed completely the point of my and of Bruno's post.
You're on your own.