[SciPy-user] Possible bug in scipy.stats.rv_discrete

Michael Nandris mnandris@btinternet....
Tue Aug 14 07:08:45 CDT 2007


Hi,

I think there may be an issue with the way rv_discrete orders its output when it encounters zeros in the input, causing non-zero probabilities to creep up towards the end of the output. If you are attempting to track the accumulation of states in an n-state Markov chain, this is a problem! 

I have had a look at the API but can't figure it. Any help at fixing this would be much appreciated.

regards

M.Nandris 

###########
# test.py

from __future__ import division
from scipy.stats import rv_discrete

# looks at output of scipy.stats.rv_discrete 

STATES = [0,1,2,3]
SIZE = 10000

def count(inpt):
    opt = dict(zip( STATES, (0,0,0,0) )) 
    for i in inpt:
    opt[i]+=1
    while opt:
    k,v = opt.popitem()
    print k, ' : ', v/SIZE # probability should reflect that of the input distribution 

def bugdemo():
    test0 = rv_discrete( name='sample', values=( STATES, [ 0.3, 0.4, 0.2, 0.1 ] ) ).rvs( size=SIZE )
    count(test0)
    print 'Output .3,.4,.2,.1 approximately matches input in the correct order: the problem only occurs when zeros are included in the initial distribution (see below)';print
    test1 = rv_discrete( name='sample', values=( STATES, [ 0.5, 0.4, 0.0, 0.1 ] ) ).rvs( size=SIZE )
    count(test1)
    print 'State 1 and State 2 have been mixed up.'; print 
    test2 = rv_discrete( name='sample', values=( STATES, [ 0.6, 0.0, 0.0, 0.4 ] ) ).rvs( size=SIZE )
    count(test2)
    print 'State 2 is sampled with the probability that State 0 should be sampled.'; print 
    test3 = rv_discrete( name='sample', values=( STATES, [ 0.0, 1.0, 0.0, 0.0 ] ) ).rvs( size=SIZE )
    count(test3)
    print 'State 3 appears to have the probability State 1 should have.'
    

if __name__=='__main__':
    bugdemo()
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://projects.scipy.org/pipermail/scipy-user/attachments/20070814/7678a8e2/attachment.html 


More information about the SciPy-user mailing list