{"action":"created","issue":{"url":"https://api.github.com/repos/pymc-devs/pymc/issues/178","repository_url":"https://api.github.com/repos/pymc-devs/pymc","labels_url":"https://api.github.com/repos/pymc-devs/pymc/issues/178/labels{/name}","comments_url":"https://api.github.com/repos/pymc-devs/pymc/issues/178/comments","events_url":"https://api.github.com/repos/pymc-devs/pymc/issues/178/events","html_url":"https://github.com/pymc-devs/pymc/issues/178","id":325155926,"number":178,"title":"Unable to define Dirichlet prior","user":{"login":"helmutsimon","id":14212177,"avatar_url":"https://avatars2.githubusercontent.com/u/14212177?v=4","gravatar_id":"","url":"https://api.github.com/users/helmutsimon","html_url":"https://github.com/helmutsimon","followers_url":"https://api.github.com/users/helmutsimon/followers","following_url":"https://api.github.com/users/helmutsimon/following{/other_user}","gists_url":"https://api.github.com/users/helmutsimon/gists{/gist_id}","starred_url":"https://api.github.com/users/helmutsimon/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/helmutsimon/subscriptions","organizations_url":"https://api.github.com/users/helmutsimon/orgs","repos_url":"https://api.github.com/users/helmutsimon/repos","events_url":"https://api.github.com/users/helmutsimon/events{/privacy}","received_events_url":"https://api.github.com/users/helmutsimon/received_events","type":"User","site_admin":false},"labels":[],"state":"closed","locked":false,"assignee":null,"assignees":[],"milestone":null,"comments":3,"created_at":"2018-05-22T06:25:39Z","updated_at":"2018-05-24T09:49:08Z","closed_at":"2018-05-24T09:49:08Z","author_association":"NONE","body":"Hello\r\n\r\nI am trying to build a simple model to estimate parameters of a multinomial distribution using a Dirichlet prior. Using the inbuilt Dirichlet distribution as in the following code seems to work OK as follows (JUpyter notebook also attached).\r\n\r\n```python\r\nimport numpy as np\r\nimport pymc\r\nimport sys\r\nfrom pymc import *\r\n\r\n\r\ndef simple_model(n, data, numseg, prior):\r\n    probs = Dirichlet('probs', prior)\r\n    process = Multinomial('process', numseg, probs, value=data, observed=True)\r\n    return locals()\r\n\r\n\r\nprint('numpy version ', np.__version__)\r\nprint('pymc version ', pymc.__version__)\r\nprint(sys.version)\r\nn = 4\r\ndata = np.array([5, 7, 1, 0])\r\nnumseg = np.sum(data)\r\nprior = np.ones(n) / n\r\nprint('MLE = ', data / np.sum(data))\r\nmodel = MCMC(simple_model(n, data, numseg, prior))\r\nmodel.sample(iter=100000, burn=1000, thin=10)\r\nmeans = model.stats()['probs']['mean']\r\nprint('\\nmeans = ', means)\r\nprint('sds = ', model.stats()['probs']['standard deviation'])\r\n```\r\n\r\nHowever, I am trying to replace the inbuilt Dirichlet distribution by a user-defined one (which I plan to modify later) as per below, and the samples never leave thee starting value. Am I doing something wrong? Environment is:\r\nnumpy version  1.11.3\r\npymc version  2.3.6\r\n3.5.4 |Continuum Analytics, Inc.| (default, Aug 14 2017, 12:43:10) \r\n[GCC 4.2.1 Compatible Apple LLVM 6.0 (clang-600.0.57)]\r\n\r\n```python\r\nimport numpy as np\r\nfrom pymc import *\r\n\r\n\r\ndef simple_model(n, data, numseg, prior):\r\n    \r\n    @stochastic(dtype='float64')\r\n    def probs(value=prior, pms=prior):\r\n        return stats.dirichlet.logpdf(value, pms)\r\n    \r\n    process = Multinomial('process', numseg, probs, value=data, observed=True)\r\n    return locals()\r\n\r\n\r\nn = 4\r\ndata = np.array( [5, 7, 1, 0])\r\nnumseg = np.sum(data)\r\nprior = np.ones(n) / n\r\nprint('MLE = ', data / np.sum(data))\r\nprint('numpy version ', np.__version__)\r\nmodel = MCMC(simple_model(n, data, numseg, prior))\r\nmodel.sample(iter=1000, burn=100, thin=10)\r\nmeans = model.stats()['probs']['mean']\r\nprint('\\nmeans = ', means)\r\nprint('sds = ', model.stats()['probs']['standard deviation'])\r\n[Pymc2_model_issues.txt](https://github.com/pymc-devs/pymc/files/2025215/Pymc2_model_issues.txt)\r\n```"},"comment":{"url":"https://api.github.com/repos/pymc-devs/pymc/issues/comments/391656026","html_url":"https://github.com/pymc-devs/pymc/issues/178#issuecomment-391656026","issue_url":"https://api.github.com/repos/pymc-devs/pymc/issues/178","id":391656026,"user":{"login":"helmutsimon","id":14212177,"avatar_url":"https://avatars2.githubusercontent.com/u/14212177?v=4","gravatar_id":"","url":"https://api.github.com/users/helmutsimon","html_url":"https://github.com/helmutsimon","followers_url":"https://api.github.com/users/helmutsimon/followers","following_url":"https://api.github.com/users/helmutsimon/following{/other_user}","gists_url":"https://api.github.com/users/helmutsimon/gists{/gist_id}","starred_url":"https://api.github.com/users/helmutsimon/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/helmutsimon/subscriptions","organizations_url":"https://api.github.com/users/helmutsimon/orgs","repos_url":"https://api.github.com/users/helmutsimon/repos","events_url":"https://api.github.com/users/helmutsimon/events{/privacy}","received_events_url":"https://api.github.com/users/helmutsimon/received_events","type":"User","site_admin":false},"created_at":"2018-05-24T09:49:08Z","updated_at":"2018-05-24T09:49:08Z","author_association":"NONE","body":"Thanks. The stick-breaking transormation did the job.  If anyone else comes across this, the corrected code is as follows:\r\n\r\n```\r\nwith Model() as model1:\r\n    prior = np.ones(n) / n\r\n    \r\n    def dirich_logpdf(value=prior):\r\n        return -n * gammaln(1/n) + (-1 + 1/n) * tt.log(value).sum()\r\n    \r\n    stick = distributions.transforms.StickBreaking()\r\n    probs = DensityDist('probs', dirich_logpdf, shape=n, testval=np.array(prior), transform=stick)\r\n    data = np.array([5, 7, 1, 0])\r\n    sfs_obs = Multinomial('sfs_obs', n=np.sum(data), p=probs, observed=data)\r\n```"}}
{"actor":{"display_login":"helmutsimon"}}