[Scipy-svn] r5280 - in trunk/scipy/stats: . tests

scipy-svn@scip... scipy-svn@scip...
Sat Dec 20 00:53:25 CST 2008


Author: josef
Date: 2008-12-20 00:53:22 -0600 (Sat, 20 Dec 2008)
New Revision: 5280

Modified:
   trunk/scipy/stats/stats.py
   trunk/scipy/stats/tests/test_stats.py
Log:
ttest_rel: some cleanup, rewrite main formula

Modified: trunk/scipy/stats/stats.py
===================================================================
--- trunk/scipy/stats/stats.py	2008-12-19 20:34:02 UTC (rev 5279)
+++ trunk/scipy/stats/stats.py	2008-12-20 06:53:22 UTC (rev 5280)
@@ -1950,7 +1950,6 @@
     zerodivproblem = svar == 0
     t = (x1-x2)/np.sqrt(svar*(1.0/n1 + 1.0/n2))  # N-D COMPUTATION HERE!!!!!!
     t = np.where(zerodivproblem, 1.0, t)           # replace NaN t-values with 1.0
-    #probs = betai(0.5*df,0.5,float(df)/(df+t*t))
     probs = distributions.t.sf(np.abs(t),df)*2
 
     if not np.isscalar(t):
@@ -2017,13 +2016,11 @@
     df = float(n-1)
     d = (a-b).astype('d')
 
-    #denom is just var(d)/df
-    denom = np.sqrt((n*np.add.reduce(d*d,axis) - np.add.reduce(d,axis)**2) /df)
+    #denom is just sqrt(var(d)/df)
+    denom = np.sqrt(np.var(d,axis,ddof=1)/float(n))
     zerodivproblem = denom == 0
-    t = np.add.reduce(d, axis) / denom      # N-D COMPUTATION HERE!!!!!!
+    t = np.mean(d, axis) / denom
     t = np.where(zerodivproblem, 1.0, t)    # replace NaN t-values with 1.0
-    t = np.where(zerodivproblem, 1.0, t)    # replace NaN t-values with 1.0
-    #probs = betai(0.5*df,0.5,float(df)/(df+t*t))
     probs = distributions.t.sf(np.abs(t),df)*2
     if not np.isscalar(t):
         probs = np.reshape(probs, t.shape)

Modified: trunk/scipy/stats/tests/test_stats.py
===================================================================
--- trunk/scipy/stats/tests/test_stats.py	2008-12-19 20:34:02 UTC (rev 5279)
+++ trunk/scipy/stats/tests/test_stats.py	2008-12-20 06:53:22 UTC (rev 5280)
@@ -1054,5 +1054,24 @@
                               np.linspace(1,100,110)+20-0.1)),
         np.array((0.20818181818181825, 0.017981441789762638)))
 
+def test_ttest_rel():
+    #regression test
+    tr,pr = 0.81248591389165692, 0.41846234511362157
+    tpr = ([tr,-tr],[pr,pr])
+
+    rvs1 = np.linspace(1,100,100)
+    rvs2 = np.linspace(1.01,99.989,100)
+    rvs1_2D = np.array([np.linspace(1,100,100), np.linspace(1.01,99.989,100)])
+    rvs2_2D = np.array([np.linspace(1.01,99.989,100), np.linspace(1,100,100)])
+
+    t,p = stats.ttest_rel(rvs1, rvs2, axis=0)
+    assert_array_almost_equal([t,p],(tr,pr))
+    t,p = stats.ttest_rel(rvs1_2D.T, rvs2_2D.T, axis=0)
+    assert_array_almost_equal([t,p],tpr)
+    t,p = stats.ttest_rel(rvs1_2D, rvs2_2D, axis=1)
+    assert_array_almost_equal([t,p],tpr)
+
+
+
 if __name__ == "__main__":
     run_module_suite()



More information about the Scipy-svn mailing list