{"action":"created","issue":{"url":"https://api.github.com/repos/pymc-devs/pymc3/issues/981","repository_url":"https://api.github.com/repos/pymc-devs/pymc3","labels_url":"https://api.github.com/repos/pymc-devs/pymc3/issues/981/labels{/name}","comments_url":"https://api.github.com/repos/pymc-devs/pymc3/issues/981/comments","events_url":"https://api.github.com/repos/pymc-devs/pymc3/issues/981/events","html_url":"https://github.com/pymc-devs/pymc3/issues/981","id":133663727,"number":981,"title":"Stochastic variable is not being sampled properly in pymc3","user":{"login":"junpenglao","id":12952641,"avatar_url":"https://avatars.githubusercontent.com/u/12952641?v=3","gravatar_id":"","url":"https://api.github.com/users/junpenglao","html_url":"https://github.com/junpenglao","followers_url":"https://api.github.com/users/junpenglao/followers","following_url":"https://api.github.com/users/junpenglao/following{/other_user}","gists_url":"https://api.github.com/users/junpenglao/gists{/gist_id}","starred_url":"https://api.github.com/users/junpenglao/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/junpenglao/subscriptions","organizations_url":"https://api.github.com/users/junpenglao/orgs","repos_url":"https://api.github.com/users/junpenglao/repos","events_url":"https://api.github.com/users/junpenglao/events{/privacy}","received_events_url":"https://api.github.com/users/junpenglao/received_events","type":"User","site_admin":false},"labels":[],"state":"closed","locked":false,"assignee":null,"milestone":null,"comments":5,"created_at":"2016-02-15T09:28:21Z","updated_at":"2016-02-16T10:38:34Z","closed_at":"2016-02-15T23:30:56Z","body":"Hi there,\r\nI am fairly new to Bayesian computation and pymc3, so this is probably a very naive question.\r\nI am trying to adapt a pymc2 code to pymc3. It is a latent-mixture model from Lee and Wagenmakers' [Bayesian Cognitive Modeling - A Pracitcal Course](https://webfiles.uci.edu/mdlee/LeeWagenmakers2013_Free.pdf) on page 93 (Chapter 6 example 7). While the pymc2 code runs without any problem, there is a problem in the pymc3 model, one of the stochastic variable is not being sampled properly...\r\n\r\npymc2 code:\r\n```python\r\nk = np.array([26,23,39,34,28,23,36,32,29,31,40,23,40,39,25,40,34,29,38,37,33,30,37,31,40\r\n,34,38,23,39,19,35,26,38,18,31,37,40,34,33,40,36,35,21,37,23,40,39,37,36,34\r\n,31,32,34,30,40,24,33,40,40,23,39,22,23,32,37,40,25,34,40,27,35,32,40,36,40\r\n,29,30,35,39,28,27,40,32,40,25,30,26,37,17,38,30,33,36,33,27,38,34,40,20,40\r\n,25,38,20,37,32,40,21,40,34,40,31,37,27,34,36,36,21,26])\r\nt1 = np.array([1,0,1,0,0,0,1,1,0,0,1,0,1,1,0,1,0,0,1,1,0,0,\r\n               1,0,1,0,1,0,1,0,1,0,1,0,0,1,1,0,0,1,1,1,0,0,\r\n               0,1,0,1,1,0,0,1,0,1,1,0,0,1,1,0,1,0,1,0,0,1,\r\n               0,0,1,0,1,0,1,0,1,0,1,0,1,0,0,1,0,1,0,1,0,1,\r\n               0,1,1,0,1,0,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,\r\n               0,1,0,1,0,1,0,1])\r\n\r\np = len(k) # number of people\r\nn = 40     # number of questions\r\n\r\nimport pymc as pymc\r\n\r\nmub = pymc.Beta(\"mub\",alpha=1,beta=1)\r\nmud = pymc.HalfNormal(\"mud\",tau=1/.5**2)\r\nlambdab = pymc.Uniform(\"lambdab\",lower=5,upper=50)\r\nlambdac = pymc.Uniform(\"lambdac\",lower=5,upper=50)\r\npsi = pymc.Beta(\"psi\",alpha=5,beta=5)\r\n\r\nzi = pymc.Bernoulli(\"zi\",p=psi,size=p)\r\n\r\n@pymc.deterministic\r\ndef muc(mub=mub,mud=mud):\r\n    return pymc.invlogit(pymc.logit(mub)+mud)\r\n\r\ntheta1 = pymc.Beta(\"theta1\",alpha=mub*lambdab,beta=(1-mub)*lambdab)\r\ntheta2 = pymc.Beta(\"theta2\",alpha=muc*lambdac,beta=(1-muc)*lambdac)\r\n# theta = theta1*(zi==0)+theta2*(zi==1)\r\ndef thetafunc(theta1=theta1,theta2=theta2,zi=zi):\r\n    return theta1*(zi==0)+theta2*(zi==1)\r\ntheta = pymc.Deterministic(eval = thetafunc,\r\n                  name = 'theta',\r\n                  parents = {'theta1': theta1,\r\n                          'theta2': theta2,\r\n                          'zi': zi},\r\n                  doc = 'Theta',\r\n                  trace = False,\r\n                  verbose = 0,\r\n                  dtype=float,\r\n                  plot=False)\r\n\r\nkij = pymc.Binomial(\"kij\",p=theta, n=n, value=k, observed = True)\r\n\r\ntestmodel = pymc.Model([theta1,theta2,psi,zi,kij])\r\nmcmc = pymc.MCMC(testmodel)\r\nmcmc.sample(iter = 20000, burn = 50, thin = 2)\r\npymc.Matplot.plot(mcmc)\r\n```\r\n\r\npymc3 code:\r\n```python\r\nimport pymc3 as pm\r\n\r\nwith pm.Model() as model7:\r\n    # prior\r\n    mub = pm.Beta(\"mub\",alpha=1,beta=1)\r\n    mud = pm.HalfNormal(\"mud\",tau=1/.5**2)\r\n    lambdab = pm.Uniform(\"lambdab\",lower=5,upper=50)\r\n    lambdac = pm.Uniform(\"lambdac\",lower=5,upper=50)\r\n    psi = pm.Beta(\"psi\",alpha=5,beta=5)\r\n    # psi = pm.Uniform(\"psi\",lower=0,upper=1)\r\n    zi = pm.Bernoulli(\"zi\",p=psi,shape=p)\r\n    \r\n    muc = pm.Deterministic(\"muc\",1/(1+T.exp(T.log(1/mub-1)-mud)))\r\n    theta1 = pm.Beta(\"theta1\",alpha=mub*lambdab,beta=(1-mub)*lambdab)\r\n    theta2 = pm.Beta(\"theta2\",alpha=muc*lambdac,beta=(1-muc)*lambdac)\r\n    \r\n    theta = pm.Deterministic(\"theta\",theta1*(T.eq(zi,0))+theta2*(T.eq(zi,1)))\r\n    # theta = theta1*(T.eq(zi,0))+theta2*(T.eq(zi,1))\r\n    # observed\r\n    kij = pm.Binomial(\"kij\",p=theta, n=n, observed = k)\r\n\r\n    start = pm.find_MAP()\r\n    start['zi']=np.random.binomial(1,.5,p)\r\n    \r\n#     step1 = pm.Metropolis(zi)\r\n#     step2 = pm.NUTS()\r\n#     trace7=pm.sample(1e4, start = start, step = [step1,step2], model=model7)\r\n    \r\n    step = pm.Metropolis()\r\n    trace7=pm.sample(1e4, start = start, step = step, model=model7)\r\n    \r\n    # step1=pm.NUTS(psi)\r\n    # step2=pm.Metropolis([theta1,theta2])\r\n    # trace7=pm.sample(1e4, step=[step1,step2], model=model7)\r\n    \r\nburnin=0\r\npm.traceplot(trace7[burnin:],varnames=['psi','theta1','theta2']);\r\nplt.show()\r\n```\r\n\r\n(I tried different sampling method but the result is the same.)\r\n\r\nThe problem seems to be that: the grouping variable vector \"zi\" was not updated in the sampling. It stuck in the initial value and never being update in each sampling.\r\n\r\n[You can have a look at the output from the notebook here](http://nbviewer.jupyter.org/gist/junpenglao/d92a1db8ee7ba4353a21)"},"comment":{"url":"https://api.github.com/repos/pymc-devs/pymc3/issues/comments/184618053","html_url":"https://github.com/pymc-devs/pymc3/issues/981#issuecomment-184618053","issue_url":"https://api.github.com/repos/pymc-devs/pymc3/issues/981","id":184618053,"user":{"login":"twiecki","id":674200,"avatar_url":"https://avatars.githubusercontent.com/u/674200?v=3","gravatar_id":"","url":"https://api.github.com/users/twiecki","html_url":"https://github.com/twiecki","followers_url":"https://api.github.com/users/twiecki/followers","following_url":"https://api.github.com/users/twiecki/following{/other_user}","gists_url":"https://api.github.com/users/twiecki/gists{/gist_id}","starred_url":"https://api.github.com/users/twiecki/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/twiecki/subscriptions","organizations_url":"https://api.github.com/users/twiecki/orgs","repos_url":"https://api.github.com/users/twiecki/repos","events_url":"https://api.github.com/users/twiecki/events{/privacy}","received_events_url":"https://api.github.com/users/twiecki/received_events","type":"User","site_admin":false},"created_at":"2016-02-16T10:38:34Z","updated_at":"2016-02-16T10:38:34Z","body":"@junpenglao I just pushed an update to the `gibbs` branch that adds Gibbs sampling to the `BinaryMetropolis` sampler. Can you update and move to that branch and try again using automatic step method assignment?\r\n\r\n![image](https://cloud.githubusercontent.com/assets/674200/13073944/ce4ea10c-d4a1-11e5-9670-f137cc0d4509.png)\r\n"}}