numerical noise question

[John H Palmieri]

How can I tell if numerical noise is actually noise, or if it is
indicative of a bug? For example, with one or two OS/processor
combinations, I get this (from chmm.pyx):

sage: m.viterbi([0,1,10,10,1])
Expected:
([0, 0, 1, 1, 0], -9.0604285688230899)
Got:
([0, 0, 1, 1, 0], -9.0604285688230917)

Can I tell how accurate this is actually supposed to be? I can
certainly just change the doctest to

([0, 0, 1, 1, 0], -9.0604285688230...)

but if the code is actually supposed to be accurate to a few more
decimal places, this is concealing a bug. Sometimes I can look at the
code and easily figure out its supposed accuracy, but more frequently
I can't. So what should be done in cases like this?

--
John

[Dr. David Kirkby]

There are as I see it two problems here with this test.

1) Do we know what a value for a high-precision result? I get rather frustrated
seeing doc tests where the only justification given for a result seems to be

"what I got with my computer".

2) Assuming the "expected" result is correct, then the value "got" on Solaris
has a relative error of

In[12]:= (-9.0604285688230899+9.0604285688230917)/-9.0604285688230899

-16
Out[12]= -1.96057 10

A relative error 1.96 x 10^-16 does not seem unreasonable to me for double
precision floating point maths. You can't generally expect to better.

However, assuming we permit a doctest of -9.0604285688230..., then a value of
-9.0604285688230 would be acceptable.

That's a relative error of

In[13]:= (-9.0604285688230899+9.0604285688230)/-9.0604285688230899

-15
Out[13]= 9.99889 10

which is 51 times worst!

In other words, in this case, to permit a small amount of numerical noise (~ 2 x
10-16), we would have to permit more than 50x what we actually got. That's why I
asked John about this test, as doing that does not seem reasonable to me.

It would be really good if we know what the exact (or at least a high precision
result was though). I really hate tests where the reason for the chosen number
seems only based on what someone got on their machine.

Dave

[William Stein]

I'm the upstream author of 100% of this code, and know precisely what
it does, which is a sequence of floating point ops and calls to the
math library, which *of course* can be machine dependent.

Put in dots.

-- William

[YannLC]

On Aug 12, 10:40 pm, "Dr. David Kirkby" <david.kir...@onetel.net>
wrote:

FWIW, here is the exact expected result in (obtained from my own
implementation, no proof but seems to be the same):
([0, 0, 1, 1, 0], -5/2*log(pi) - 15/2*log(2) - 1)

and in high precision:

([0, 0, 1, 1, 0],

-9.0604285688230902559878092893189710396844880399901933346892)

Yann

[John H Palmieri]

Great, that's very helpful. (My guess was that it was just a sequence
of floating point ops, as you say, but I wanted to check, to make sure
we were preserving the integrity of Sage's doctests.)

--
John

[Dr. David Kirkby]

Thank you Yann.

Dave