Re: [Numpy-discussion] multinomial question
> Alan G Isaac wrote: >> I would think that >> multinomial(1,prob,size=ntrials).sum(axis=0) >> would be equivalent to >> multinomial(ntrials,prob) >> but the first gives a surprising result. (See below.) >> Explanation? On Wed, 05 Dec 2007, Robert Kern apparently wrote: > Pretty much anyone who derives their binomial distribution algorithm from > http://www.unc.edu/~gfish/fcmc.html is also wrong. > SVN now has a bound such that CDF(bound) is within 1e-16 (or so) of 1.0. Thanks! Alan ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] multinomial question
Alan G Isaac wrote: > I would think that > multinomial(1,prob,size=ntrials).sum(axis=0) > would be equivalent to > multinomial(ntrials,prob) > but the first gives a surprising result. (See below.) > Explanation? Pretty much anyone who derives their binomial distribution algorithm from http://www.unc.edu/~gfish/fcmc.html is also wrong. SVN now has a bound such that CDF(bound) is within 1e-16 (or so) of 1.0. -- Robert Kern "I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] multinomial question
Alan G Isaac wrote: > I would think that > multinomial(1,prob,size=ntrials).sum(axis=0) > would be equivalent to > multinomial(ntrials,prob) > but the first gives a surprising result. (See below.) > Explanation? A bug in rk_binomial_inversion(). Unfortunately, this looks like a logical bug in the sources I was deriving this code from. The "safety bound" on the search inversion search cuts out too early. Now that I re-examine it, it looks the bound (whichever of the multiple choice of bounds one could use) could always be legitimately exceeded, so there shouldn't be a bound at all. I'll have to dive deeper to figure out what is going on. This makes me grumpy. -- Robert Kern "I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] multinomial question
I would think that multinomial(1,prob,size=ntrials).sum(axis=0) would be equivalent to multinomial(ntrials,prob) but the first gives a surprising result. (See below.) Explanation? Thank you, Alan Isaac >>> ntrials = 10 >>> prob = N.arange(100,dtype=N.float32)/4950 >>> multinomial(1,prob,size=ntrials).sum(axis=0) array([ 0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0, 990, 1058, 996, 1102, 1088, 1137, 1144, 1150, 1196, 1198, 1272, 1273, 1268, 1265, 1380, 1336, 1371, 1405, 1348, 1420, 1515, 1506, 1571, 1499, 1556, 1517, 1603, 1691, 1696, 1763, 1622, 1716, 1722, 1785, 1866, 1799, 1918, 1818, 1868, 1938, 2010, 1916, 1950, 1983, 2062, 2079, 2224, 2165, 2136, 2156, 2215, 2118, 2309, 2389, 2377, 2423, 2328, 2325, 2469]) >>> multinomial(ntrials,prob) array([ 0, 27, 33, 55, 104, 104, 116, 153, 166, 181, 189, 199, 244, 262, 259, 303, 330, 343, 373, 360, 371, 437, 423, 470, 460, 550, 551, 497, 517, 593, 623, 623, 648, 660, 638, 718, 713, 754, 784, 831, 804, 868, 902, 851, 918, 932, 945, 972, 966, 1025, 1005, 1038, 1075, 1046, 1121, 1069, 1121, 1152, 1209, 1148, 1196, 1261, 1288, 1304, 1250, 1324, 1348, 1430, 1370, 1419, 1388, 1364, 1473, 1414, 1511, 1523, 1583, 1574, 1575, 1575, 1613, 1559, 1665, 1666, 1712, 1728, 1715, 1709, 1846, 1774, 1819, 1869, 1886, 1963, 1837, 1906, 1983, 1867, 1968, 1916]) >>> ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
