{"action":"created","issue":{"url":"https://api.github.com/repos/pymc-devs/pymc3/issues/1163","repository_url":"https://api.github.com/repos/pymc-devs/pymc3","labels_url":"https://api.github.com/repos/pymc-devs/pymc3/issues/1163/labels{/name}","comments_url":"https://api.github.com/repos/pymc-devs/pymc3/issues/1163/comments","events_url":"https://api.github.com/repos/pymc-devs/pymc3/issues/1163/events","html_url":"https://github.com/pymc-devs/pymc3/issues/1163","id":158838408,"number":1163,"title":"ADVI, NUTS and Metropolis produce significantly different results","user":{"login":"kayhan-batmanghelich","id":12301533,"avatar_url":"https://avatars.githubusercontent.com/u/12301533?v=3","gravatar_id":"","url":"https://api.github.com/users/kayhan-batmanghelich","html_url":"https://github.com/kayhan-batmanghelich","followers_url":"https://api.github.com/users/kayhan-batmanghelich/followers","following_url":"https://api.github.com/users/kayhan-batmanghelich/following{/other_user}","gists_url":"https://api.github.com/users/kayhan-batmanghelich/gists{/gist_id}","starred_url":"https://api.github.com/users/kayhan-batmanghelich/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/kayhan-batmanghelich/subscriptions","organizations_url":"https://api.github.com/users/kayhan-batmanghelich/orgs","repos_url":"https://api.github.com/users/kayhan-batmanghelich/repos","events_url":"https://api.github.com/users/kayhan-batmanghelich/events{/privacy}","received_events_url":"https://api.github.com/users/kayhan-batmanghelich/received_events","type":"User","site_admin":false},"labels":[],"state":"open","locked":false,"assignee":null,"assignees":[],"milestone":null,"comments":16,"created_at":"2016-06-07T04:55:46Z","updated_at":"2016-07-07T21:08:08Z","closed_at":null,"body":"Dear PyMC developers,\r\n\r\nI am trying to develop simple mixed effect model using pymc (see the code below). I have tried NUTS, ADVI, Metropolis for inference. Aside from variations in speed for this  model (Time : AVDI ~= Metrolopis << NUTS ), the results are significantly different. To be more specific, I am interested to estimate $h^2$ which is basically proportion of variance explained by the random effect. \r\n\r\nThis is the code:\r\n\r\n```\r\nimport numpy as np, pymc3 as pm, theano.tensor as T\r\n\r\n\r\nmixedEffect_model = pm.Model()\r\n\r\nwith pm.Model() as mixedEffect_model:\r\n    ### hyperpriors\r\n    tau = pm.Gamma('tau', 1.e-3, 1.e-3)\r\n    sigma2 = 1.0/tau\r\n    h2 = pm.Uniform('h2')\r\n    beta_0 = pm.Uniform('beta_0', lower=-1000, upper=1000)   # a replacement for improper prior\r\n    w = pm.Uniform('w', lower=-1000, upper=1000, shape=(M,1))   # a replacement for improper prior\r\n    \r\n    z = pm.Normal('z', mu=0, sd= (h2*sigma2)**0.5 , shape=(N,1))\r\n    g = T.dot(L,z)\r\n    \r\n    y = pm.Normal('y', mu = beta_0+g + T.dot(X,w), sd= ((1-h2)*sigma2)**0.5 , shape=(N,1) , observed=y )\r\n    \r\n```\r\n\r\nADVI:\r\n```\r\nwith mixedEffect_model:\r\n    v_params = pm.variational.advi(n=5000)\r\n\r\nwith mixedEffect_model:\r\n    trace = pm.variational.sample_vp(v_params, draws=5000)\r\n\r\n_ = hist(trace['h2'],100)\r\n```\r\n\r\n![image](https://cloud.githubusercontent.com/assets/12301533/15840000/d61231f4-2c14-11e6-9166-093fd8c860be.png)\r\n\r\n\r\nMetropolis:\r\n```\r\nwith mixedEffect_model:\r\n    start = pm.find_MAP()\r\n    step = pm.Metropolis()\r\n    trace = pm.sample(5000, step, start=start)\r\n\r\n_ = hist(trace['h2'],100)\r\n```\r\n\r\n![image](https://cloud.githubusercontent.com/assets/12301533/15840009/ec4a6176-2c14-11e6-9349-a1cd61e3a97c.png)\r\n\r\n\r\nNUTS (slow):\r\n```\r\nwith mixedEffect_model:\r\n    start = pm.find_MAP()\r\n    step = pm.NUTS()\r\n    trace = pm.sample(5000, step, start=start)\r\n    \r\n_ = hist(trace['h2'],100)    \r\n```\r\n\r\n![image](https://cloud.githubusercontent.com/assets/12301533/15842995/2e1c2b36-2c2c-11e6-8367-3480819839a0.png)\r\n\r\nGiven that `w` and `beta_0` are treated as parameter with (approximately) improper prior, I expect the results to be very close to maximum likelihood estimation.  Using maximum likelihood estimation (more specifically Restricted maximum Likelihood estimation), I expect values around 0.5 for $h^2$ with standard error of around 0.07. The results is very different. It is very different from ML and also very different across different inference engines. NUTS is closer to ML but given that the results are drastically different, I don't know how to explain it.\r\n\r\nAny idea?"},"comment":{"url":"https://api.github.com/repos/pymc-devs/pymc3/issues/comments/231207816","html_url":"https://github.com/pymc-devs/pymc3/issues/1163#issuecomment-231207816","issue_url":"https://api.github.com/repos/pymc-devs/pymc3/issues/1163","id":231207816,"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},"created_at":"2016-07-07T21:08:08Z","updated_at":"2016-07-07T21:08:08Z","body":"@hvasbath I see. But I \"...would need to add a likelihood variable to your model that contains the model likelihood and pass it to the input.\" Is the likelihood variable being pass to the following need to be the one with observed?\r\n```\r\nwith ATMIP_test:\r\n    step = atmcmc.ATMCMC(n_chains=n_chains, tune_interval=tune_interval,\r\n                         likelihood_name=ATMIP_test.deterministics[0].name)\r\n```\r\nor I can also do \r\n`llk = pm.Potential('like', like)`"}}
{"actor":{"display_login":"junpenglao"}}