[NumPy-Tickets] [NumPy] #1801: Failure: test_noncentral_f in random

[NumPy Trac]

#1801: Failure: test_noncentral_f in random
--------------------------+-------------------------------------------------
 Reporter:  rgommers      |       Owner:  somebody
     Type:  defect        |      Status:  new     
 Priority:  normal        |   Milestone:  1.6.0   
Component:  numpy.random  |     Version:  1.5.1   
 Keywords:                |  
--------------------------+-------------------------------------------------

Comment(by cgohlke):

 The following code returns unexpected results with any of my 64 bit msvc9
 builds of numpy 1.5.1 or 1.6.x, and also with EPD 7.0-1 64 bit:

 {{{
 >>> import numpy as np
 >>> np.random.seed(1234567890)
 >>> np.random.noncentral_f(dfnum=5, dfden=2, nonc=1.0, size=(3, 2))
 array([[ 1.62003345,  1.7253997 ],
        [ 0.96735921,  0.42933718],
        [ 0.71714872,  6.24979552]])
 }}}

 Expected result:

 {{{
 array([[ 1.405981  ,  0.34207973],
        [ 3.57715069,  7.92632663],
        [ 0.43741599,  1.17742088]])
 }}}

 The code fails whether MKL is used/linked or not.

 The 32 bit msvc9 and 64 bit "Windows Server 2003 R2 Platform SDK" builds
 are OK. So this seems to be a 64 bit msvc9 specific issue.

 Other numpy.random tests pass with the same seed. Specifically, chisquare
 and noncentral_chisquare return expected results.

 The noncentral_f function is implemented in
 numpy\random\mtrand\distributions.c as

 {{{
 double rk_noncentral_f(rk_state *state, double dfnum, double dfden, double
 nonc)
 {
     return ((rk_noncentral_chisquare(state, dfnum, nonc)*dfden) /
             (rk_chisquare(state, dfden)*dfnum));
 }
 }}}

 I changed the function to the following, which does return the correct
 results:

 {{{
 double rk_noncentral_f(rk_state *state, double dfnum, double dfden, double
 nonc)
 {
     double t1, t2;
     t1 = rk_noncentral_chisquare(state, dfnum, nonc) * dfden;
     t2 = rk_chisquare(state, dfden) * dfnum;
     t1 /= t2;
     return t1;
 }
 }}}

 So far I only tested this with numpy-1.5.1.win-amd64-py2.6.

-- 
Ticket URL: <http://projects.scipy.org/numpy/ticket/1801#comment:1>
NumPy <http://projects.scipy.org/numpy>
My example project