[Scipy-svn] r5074 - trunk/scipy/stats

scipy-svn at scipy.org scipy-svn at scipy.org
Wed Nov 12 21:30:47 EST 2008


Author: josef
Date: 2008-11-12 20:30:45 -0600 (Wed, 12 Nov 2008)
New Revision: 5074

Modified:
   trunk/scipy/stats/distributions.py
Log:
correct function _drv2_moment which didn't do anything before 

Modified: trunk/scipy/stats/distributions.py
===================================================================
--- trunk/scipy/stats/distributions.py	2008-11-13 01:45:30 UTC (rev 5073)
+++ trunk/scipy/stats/distributions.py	2008-11-13 02:30:45 UTC (rev 5074)
@@ -3267,15 +3267,36 @@
     return sum(exp(self.xk * t[newaxis,...]) * self.pk, axis=0)
 
 def _drv2_moment(self, n, *args):
+    '''non-central moment of discrete distribution'''
+    #many changes, originally not even a return
     tot = 0.0
     diff = 1e100
-    pos = self.a
+    #pos = self.a
+    pos = max(0, self.a)
     count = 0
-    while (pos <= self.b) and ((pos >= (self.b + self.a)/2.0) and \
+    #handle cases with infinite support 
+    ulimit = max(1000, (min(self.b,1000) + max(self.a,-1000))/2.0 )
+    llimit = min(-1000, (min(self.b,1000) + max(self.a,-1000))/2.0 )
+    
+    while (pos <= self.b) and ((pos <= ulimit) or \
                                (diff > self.moment_tol)):
-        diff = pos**n * self._pdf(pos,*args)
+        diff = pos**n * self.pmf(pos,*args) 
+        # use pmf because _pmf does not check support in randint
+        #     and there might be problems ? with correct self.a, self.b at this stage
         tot += diff
         pos += self.inc
+        count += 1
+        
+    if self.a < 0: #handle case when self.a = -inf
+        diff = 1e100
+        pos = -self.inc
+        while (pos >= self.a) and ((pos >= llimit) or \
+                                   (diff > self.moment_tol)):
+            diff = pos**n * self.pmf(pos,*args)  #using pmf instead of _pmf
+            tot += diff
+            pos -= self.inc
+            count += 1
+    return tot
 
 def _drv2_ppfsingle(self, q, *args):  # Use basic bisection algorithm
     b = self.invcdf_b




More information about the Scipy-svn mailing list