From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: (qmail 28063 invoked by alias); 17 Nov 2008 09:14:36 -0000 Received: (qmail 4834 invoked by uid 22791); 16 Nov 2008 23:36:49 -0000 X-Spam-Check-By: sourceware.org X-IronPort-AV: E=Sophos;i="4.33,615,1220191200"; d="h'?c'?scan'208";a="3276229" From: To: Date: Mon, 17 Nov 2008 09:14:00 -0000 Subject: Multinomial CDF Message-ID: Accept-Language: en-US Content-Language: en-US acceptlanguage: en-US Content-Type: multipart/mixed; boundary="_004_E5DDD62B446B4E4B9399046EAD0144A6242A014AEXNSWMBX02nexus_" MIME-Version: 1.0 Mailing-List: contact gsl-discuss-help@sourceware.org; run by ezmlm Precedence: bulk List-Id: List-Subscribe: List-Archive: List-Post: List-Help: , Sender: gsl-discuss-owner@sourceware.org X-SW-Source: 2008-q4/txt/msg00021.txt.bz2 --_004_E5DDD62B446B4E4B9399046EAD0144A6242A014AEXNSWMBX02nexus_ Content-Type: text/plain; charset="us-ascii" Content-Transfer-Encoding: quoted-printable Content-length: 808 Hi, I have written a piece of code to implement the multinomial CDF. GSL curre= ntly provides random variates from a multinomial and the multinomial PDF. T= he CDF is implemented as the approximation from, "A Representation for Multinomial Cumulative Distribution Functions", Bruc= e Levin, The Annals of Statistics, v.9, n.5, pp.1123-1126, 1981 The approximation is meant to be pretty good. Because the multinomial is a vector random variate the "P" version of the C= DF is implemented. The P version corresponds to P(X1<=3Dn1,X2<=3Dn2,...,XK<=3DnK) In keeping with the pdf function already in GSL no checking of the consiste= ncy/validity of the inputs is done. Code attached along with a small test function. Glenn Stone CSIRO Mathematical and Information Sciences North Ryde +61 2 9325 3259 --_004_E5DDD62B446B4E4B9399046EAD0144A6242A014AEXNSWMBX02nexus_ Content-Type: text/plain; name="test.c" Content-Description: test.c Content-Disposition: attachment; filename="test.c"; size=830; creation-date="Mon, 17 Nov 2008 09:52:53 GMT"; modification-date="Mon, 17 Nov 2008 09:52:53 GMT" Content-Transfer-Encoding: base64 Content-length: 1127 I2luY2x1ZGUgPHN0ZGlvLmg+CiNpbmNsdWRlICJtdWx0aW5vbWlhbC5oIgoK LyogdGhpcyBpcyBhIHRlc3QgZXhhbXBsZSB0YWtlbiBmcm9tIHRoZSBwYXBl ciAsc2VlIG11bHRpbm9taWFsLmMKICogZ2NjIC1XYWxsIC1zdGF0aWMgLWcg dGVzdC5jIG11bHRpbm9taWFsLmMgLWxnc2wgLWxtCiAqIG9yIAogKiBnY2Mg LVdhbGwgLVhsaW5rZXIgLS1hbGxvdy1zaGxpYi11bmRlZmluZWQgLWcgdGVz dC5jIG11bHRpbm9taWFsLmMgLWxnc2wgLWxtCiAqIChOb3Qgc3VyZSB3aGF5 IGl0IG5lZWRzIHRoZSBhbGxvdy1zaGxpYi11bmRlZmluZWQsIGJ1dCBvdGhl cndpc2UgbG90cyBvZiB1bnJlc29sdmVkIGNibGFzIGVycm9ycykKICovCgpp bnQgbWFpbigpIHsKICBkb3VibGUgcFsxMl0sIHJlcywgbGIsIHViOwogIHVu c2lnbmVkIGludCBuWzEyXTsKICBpbnQgaTsKCiAgZm9yKGk9MDsgaTwxMjsg aSsrICkgewogICAgcFtpXSA9IDEvMTIuMDsKICAgIG5baV0gPSAzOwogIH0K CiAgcmVzID0gY2RmX211bHRpbm9taWFsX1AoMTIsMTIscCxuKTsKCiAgLyog UmVzdWx0IHNob3VsZCBiZSAwLjgzNjcgKi8KICBwcmludGYoIiVnXG4iLCBy ZXMpOwoKICAvKiBSZXN1bHQgc2hvdWxkIGJlICgwLjgzNDAsMC44NDYxKSAq LwogIGNkZl9tdWx0aW5vbWlhbF9CTWJvdW5kc1AoMTIsMTIscCxuLCZsYiwm dWIpOwogIHByaW50ZigiKCVnLCVnKVxuIiwgbGIsIHViKTsKCgogIHJlcyA9 IGNkZl9tdWx0aW5vbWlhbF9QbWF4KDEyLDEyLDMpOwogIC8qIFJlc3VsdCBz aG91bGQgYmUgMC44MzY3ICovCiAgcHJpbnRmKCIlZ1xuIiwgcmVzKTsKCiAg cmV0dXJuKDApOwp9CiAgCgogIAo= --_004_E5DDD62B446B4E4B9399046EAD0144A6242A014AEXNSWMBX02nexus_ Content-Type: text/plain; name="multinomial.c" Content-Description: multinomial.c Content-Disposition: attachment; filename="multinomial.c"; size=3786; creation-date="Mon, 17 Nov 2008 09:54:23 GMT"; modification-date="Mon, 17 Nov 2008 09:54:23 GMT" Content-Transfer-Encoding: base64 Content-length: 5133 I2luY2x1ZGUgPG1hdGguaD4KI2luY2x1ZGUgPGdzbC9nc2xfbWF0aC5oPgoj aW5jbHVkZSA8Z3NsL2dzbF9lcnJuby5oPgojaW5jbHVkZSA8Z3NsL2dzbF9j ZGYuaD4KI2luY2x1ZGUgPGdzbC9nc2xfcmFuZGlzdC5oPgoKLyogQ29tcHV0 ZXMgdGhlIChsb2cpIENERiBvZiBhIG11bHRpbm9taWFsIGRpc3RyaWJ1dGlv biB3aXRoIHBhcmFtZXRlcnMgcCBhbmQgTgogKiBpLmUuIFAoWDE8PW4xLFgy PD1uMiwuLi4sWEs8PW5LKQogKiBVc2VzIGEgYXBwcm94aW1hdGlvbiBmcm9t IAogKiAiQSBSZXByZXNlbnRhdGlvbiBmb3IgTXVsdGlub21pYWwgQ3VtdWxh dGl2ZSBEaXN0cmlidXRpb24gRnVuY3Rpb25zIiwgCiAqIEJydWNlIExldmlu LCBUaGUgQW5uYWxzIG9mIFN0YXRpc3RpY3MsIHYuOSwgbi41LCBwcC4xMTIz LTExMjYsIDE5ODEKICovCgpkb3VibGUKY2RmX211bHRpbm9taWFsX2xuUCAo Y29uc3Qgc2l6ZV90IEssIGNvbnN0IHVuc2lnbmVkIGludCBOLAogICAgICAg ICAgICAgICAgICAgICAgICAgICBjb25zdCBkb3VibGUgcFtdLCBjb25zdCB1 bnNpZ25lZCBpbnQgbltdKQp7CiAgc2l6ZV90IGs7CiAgZG91YmxlIGxvZ19j ZGYgPSAwLjA7CiAgZG91YmxlIHM7CiAgZG91YmxlIGdhbW1hMT0wLjAsIGdh bW1hMj0wLjAsIHN1bV9zMj0wLjAsIHN1bV9tdT0wLjA7CiAgZG91YmxlIG11 LCBzMiwgc3AsIHBjZGYsIG1yLCBtZjIsIG1mMywgbWY0LCBtdTIsIG11Mywg bXU0LCB4LCB4MiwgUFdOOwogIAogIHMgPSAoZG91YmxlKSBOOwogIGxvZ19j ZGYgPSAtbG9nKGdzbF9yYW5fcG9pc3Nvbl9wZGYoTixzKSk7CgogIC8qIFRo aXMgaXMgdGhlIFAoVz1OKSBiaXQgKi8KICBmb3Ioaz0wOyBrPEs7IGsrKykg ewogICAgc3AgPSBzKnBba107CgogICAgcGNkZiA9IGdzbF9jZGZfcG9pc3Nv bl9QKG5ba10sc3ApOwogICAgbG9nX2NkZiArPSBsb2cocGNkZik7CgogICAg bXUgPSBzcCooMS1nc2xfcmFuX3BvaXNzb25fcGRmKG5ba10sc3ApL3BjZGYp OwogICAgczIgPSBtdS0obltrXS1tdSkqKHNwLW11KTsKCiAgICAvKiBmYWN0 b3JpYWwgbW9tZW50cz8gKi8KICAgIG1yID0gbltrXTsKICAgIG1mMj0gc3Aq bXUtbXIqKHNwLW11KTsKCiAgICBtciAqPSBuW2tdLTE7CiAgICBtZjMgPSBz cCptZjItbXIqKHNwLW11KTsKICAgIAogICAgbXIgKj0gbltrXS0yOwogICAg bWY0ID0gc3AqbWYzLW1yKihzcC1tdSk7CgogICAgLyogQ2VudHJhbCBNb21l bnRzICovCiAgICBtdTIgPSBtZjIrbXUqKDEtbXUpOwogICAgbXUzID0gbWYz K21mMiooMy0zKm11KSttdSooMSttdSooLTMrMiptdSkpOwogICAgbXU0ID0g bWY0K21mMyooNi00Km11KSttZjIqKDcrbXUqKC0xMis2Km11KSkrbXUqKDEr bXUqKC00K211Kig2LTMqbXUpKSk7CgogICAgLyogYWNjdW11bGF0ZSBjb2Vm IHNrZXduZXNzIGFuZCBleGNlc3MgKi8KICAgIGdhbW1hMSArPSBtdTM7CiAg ICBnYW1tYTIgKz0gbXU0LTMqczIqczI7CiAgICBzdW1fbXUgKz0gbXU7CiAg ICBzdW1fczIgKz0gczI7CiAgfQogIAogIHNwID0gc3FydChzdW1fczIpOwog IGdhbW1hMSAvPSBzdW1fczIqc3A7CiAgZ2FtbWEyIC89IHN1bV9zMipzdW1f czI7CgogIHggPSAoTi1zdW1fbXUpL3NwOwogIHgyID0geCp4OwogIFBXTiA9 IC14Mi8yCiAgICArbG9nKDErZ2FtbWExLzYqeCooeDItMykrZ2FtbWEyLzI0 Kih4Mip4Mi02KngyKzMpKwoJIGdhbW1hMSpnYW1tYTEvNzIqKCgoeDItMTUp KngyKzQ1KSp4Mi0xNSkpCiAgICAtbG9nKDIqTV9QSSkvMiAtbG9nKHNwKTsK CiAgbG9nX2NkZiArPSBQV047CgogIHJldHVybiBsb2dfY2RmOwp9Cgpkb3Vi bGUKY2RmX211bHRpbm9taWFsX1AgKGNvbnN0IHNpemVfdCBLLCBjb25zdCB1 bnNpZ25lZCBpbnQgTiwKICAgICAgICAgICAgICAgICAgICAgICAgIGNvbnN0 IGRvdWJsZSBwW10sIGNvbnN0IHVuc2lnbmVkIGludCBuW10pCnsKICByZXR1 cm4gZXhwIChjZGZfbXVsdGlub21pYWxfbG5QIChLLCBOLCBwLCBuKSk7Cn0K Ci8qIENvbXB1dGVzIHRoZSBCb25mZXJyb25pLU1hbGxvd3MgYm91bmRzIGZv ciB0aGUgQ0RGICovCgp2b2lkCmNkZl9tdWx0aW5vbWlhbF9CTWJvdW5kc1Ag KGNvbnN0IHNpemVfdCBLLCBjb25zdCB1bnNpZ25lZCBpbnQgTiwKCQkJICAg Y29uc3QgZG91YmxlIHBbXSwgY29uc3QgdW5zaWduZWQgaW50IG5bXSwgZG91 YmxlICpsYiwgZG91YmxlICp1YikKewogIGRvdWJsZSBMQiA9IDEuMCwgVUIg PSAxLjA7CiAgc2l6ZV90IGs7CiAgCiAgZm9yKGs9MDsgazxLOyBrKyspIHsK ICAgIExCIC09IGdzbF9jZGZfYmlub21pYWxfUShuW2tdLHBba10sTik7CiAg ICBVQiAqPSBnc2xfY2RmX2Jpbm9taWFsX1AobltrXSxwW2tdLE4pOyAKICB9 CgogICpsYiA9IExCOwogICp1YiA9IFVCOwp9CiAgCgogIAovKiBTcGVjaWFs IGNhc2Ugd2l0aCBhbGwgcCBlcXVhbCBhbmQgYWxsIG4gZXF1YWwKICogZWku IFAobWF4X2sgWF9rIDw9IG4pIHdoZW4gcF9rID0gcCBmb3IgYWxsIGsKICov Cgpkb3VibGUKY2RmX211bHRpbm9taWFsX2xuUG1heCAoY29uc3Qgc2l6ZV90 IEssIGNvbnN0IHVuc2lnbmVkIGludCBOLAoJCQljb25zdCB1bnNpZ25lZCBp bnQgbikKewogIGRvdWJsZSBsb2dfY2RmID0gMC4wOwogIGRvdWJsZSBwLCBz OwogIGRvdWJsZSBnYW1tYTEsIGdhbW1hMjsKICBkb3VibGUgbXUsIHMyLCBz cCwgcGNkZiwgbXIsIG1mMiwgbWYzLCBtZjQsIG11MiwgbXUzLCBtdTQsIHgs IHgyLCBQV047CiAgCiAgcCA9IDEvKChkb3VibGUpIEspOwogIHMgPSAoZG91 YmxlKSBOOwoKICBsb2dfY2RmID0gLWxvZyhnc2xfcmFuX3BvaXNzb25fcGRm KE4scykpOwoKICBzcCA9IHMqcDsKICBwY2RmID0gZ3NsX2NkZl9wb2lzc29u X1AobixzcCk7CiAgbG9nX2NkZiArPSBLKmxvZyhwY2RmKTsKCiAgLyogVGhp cyBpcyB0aGUgUChXPU4pIGJpdCAqLwogIG11ID0gc3AqKDEtZ3NsX3Jhbl9w b2lzc29uX3BkZihuLHNwKS9wY2RmKTsKICBzMiA9IG11LShuLW11KSooc3At bXUpOwoKICAvKiBmYWN0b3JpYWwgbW9tZW50cz8gKi8KICBtciA9IG47CiAg bWYyPSBzcCptdS1tciooc3AtbXUpOwoKICBtciAqPSBuLTE7CiAgbWYzID0g c3AqbWYyLW1yKihzcC1tdSk7CiAgICAKICBtciAqPSBuLTI7CiAgbWY0ID0g c3AqbWYzLW1yKihzcC1tdSk7CgogICAgLyogQ2VudHJhbCBNb21lbnRzICov CiAgbXUyID0gbWYyK211KigxLW11KTsKICBtdTMgPSBtZjMrbWYyKigzLTMq bXUpK211KigxK211KigtMysyKm11KSk7CiAgbXU0ID0gbWY0K21mMyooNi00 Km11KSttZjIqKDcrbXUqKC0xMis2Km11KSkrbXUqKDErbXUqKC00K211Kig2 LTMqbXUpKSk7CgogICAgLyogYWNjdW11bGF0ZSBjb2VmIHNrZXduZXNzIGFu ZCBleGNlc3MgKi8KICBzcCA9IHNxcnQoSypzMik7CiAgZ2FtbWExID0gbXUz LyhzcCpzMik7CiAgZ2FtbWEyID0gKG11NC0zKnMyKnMyKS8oczIqczIqSyk7 CgogIHggPSAoTi1LKm11KS9zcDsKICB4MiA9IHgqeDsKICBQV04gPSAteDIv MgogICAgK2xvZygxK2dhbW1hMS82KngqKHgyLTMpK2dhbW1hMi8yNCooeDIq eDItNip4MiszKSsKCSBnYW1tYTEqZ2FtbWExLzcyKigoKHgyLTE1KSp4Mis0 NSkqeDItMTUpKQogICAgLWxvZygyKk1fUEkpLzIgLWxvZyhzcCk7CgogIGxv Z19jZGYgKz0gUFdOOwoKICByZXR1cm4gbG9nX2NkZjsKfQoKZG91YmxlCmNk Zl9tdWx0aW5vbWlhbF9QbWF4IChjb25zdCBzaXplX3QgSywgY29uc3QgdW5z aWduZWQgaW50IE4sCgkJICAgICAgY29uc3QgdW5zaWduZWQgaW50IG4pCnsK ICByZXR1cm4gZXhwIChjZGZfbXVsdGlub21pYWxfbG5QbWF4IChLLCBOLCBu KSk7Cn0K --_004_E5DDD62B446B4E4B9399046EAD0144A6242A014AEXNSWMBX02nexus_ Content-Type: text/plain; name="multinomial.h" Content-Description: multinomial.h Content-Disposition: attachment; filename="multinomial.h"; size=604; creation-date="Mon, 17 Nov 2008 09:53:42 GMT"; modification-date="Mon, 17 Nov 2008 09:53:42 GMT" Content-Transfer-Encoding: base64 Content-length: 822 ZXh0ZXJuIGRvdWJsZQpjZGZfbXVsdGlub21pYWxfbG5QIChjb25zdCBzaXpl X3QgSywgY29uc3QgdW5zaWduZWQgaW50IE4sCgkJICAgICBjb25zdCBkb3Vi bGUgcFtdLCBjb25zdCB1bnNpZ25lZCBpbnQgbltdKTsKCmV4dGVybiBkb3Vi bGUKY2RmX211bHRpbm9taWFsX1AgKGNvbnN0IHNpemVfdCBLLCBjb25zdCB1 bnNpZ25lZCBpbnQgTiwKCQkgICBjb25zdCBkb3VibGUgcFtdLCBjb25zdCB1 bnNpZ25lZCBpbnQgbltdKTsKCmV4dGVybiB2b2lkCmNkZl9tdWx0aW5vbWlh bF9CTWJvdW5kc1AgKGNvbnN0IHNpemVfdCBLLCBjb25zdCB1bnNpZ25lZCBp bnQgTiwKCQkJICAgY29uc3QgZG91YmxlIHBbXSwgY29uc3QgdW5zaWduZWQg aW50IG5bXSwgZG91YmxlICpsYiwgZG91YmxlICp1Yik7CgoKZXh0ZXJuIGRv dWJsZQpjZGZfbXVsdGlub21pYWxfbG5QbWF4IChjb25zdCBzaXplX3QgSywg Y29uc3QgdW5zaWduZWQgaW50IE4sCgkJCWNvbnN0IHVuc2lnbmVkIGludCBu KTsKCmV4dGVybiBkb3VibGUKY2RmX211bHRpbm9taWFsX1BtYXggKGNvbnN0 IHNpemVfdCBLLCBjb25zdCB1bnNpZ25lZCBpbnQgTiwKCQkgICAgICBjb25z dCB1bnNpZ25lZCBpbnQgbik7Cg== --_004_E5DDD62B446B4E4B9399046EAD0144A6242A014AEXNSWMBX02nexus_--