{"action":"created","issue":{"url":"https://api.github.com/repos/pymc-devs/pymc3/issues/1501","repository_url":"https://api.github.com/repos/pymc-devs/pymc3","labels_url":"https://api.github.com/repos/pymc-devs/pymc3/issues/1501/labels{/name}","comments_url":"https://api.github.com/repos/pymc-devs/pymc3/issues/1501/comments","events_url":"https://api.github.com/repos/pymc-devs/pymc3/issues/1501/events","html_url":"https://github.com/pymc-devs/pymc3/issues/1501","id":187579952,"number":1501,"title":"Out-of-range value for uniform prior","user":{"login":"sebp","id":74460,"avatar_url":"https://avatars.githubusercontent.com/u/74460?v=3","gravatar_id":"","url":"https://api.github.com/users/sebp","html_url":"https://github.com/sebp","followers_url":"https://api.github.com/users/sebp/followers","following_url":"https://api.github.com/users/sebp/following{/other_user}","gists_url":"https://api.github.com/users/sebp/gists{/gist_id}","starred_url":"https://api.github.com/users/sebp/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/sebp/subscriptions","organizations_url":"https://api.github.com/users/sebp/orgs","repos_url":"https://api.github.com/users/sebp/repos","events_url":"https://api.github.com/users/sebp/events{/privacy}","received_events_url":"https://api.github.com/users/sebp/received_events","type":"User","site_admin":false},"labels":[{"id":76646,"url":"https://api.github.com/repos/pymc-devs/pymc3/labels/defects","name":"defects","color":"e10c02","default":false}],"state":"open","locked":false,"assignee":null,"assignees":[],"milestone":null,"comments":11,"created_at":"2016-11-06T18:43:10Z","updated_at":"2016-12-01T14:51:57Z","closed_at":null,"body":"# Problem Description\r\n\r\nWhile attempting to port code related to [Bayesian estimation of the correlation coefficient](http://www.philippsinger.info/?p=581) from PyMC2 and PyMC3, I noticed that `rho` (the correlation coefficient) can be out of range, which in turn causes NaN in the multivariate normal distribution when taking the log of a negative value ([line 117](https://github.com/pymc-devs/pymc3/blob/master/pymc3/distributions/multivariate.py#L117) of pymc3/distributions/multivariate.py).\r\n\r\nThe model is defined as below:\r\n```python\r\nfrom theano.printing import Print\r\nimport theano.tensor as T\r\nimport pymc3 as pm\r\nimport numpy as np\r\n\r\n\r\ndef covariance(sigma, rho):\r\n    C = T.fill_diagonal(T.alloc(rho, 2, 2), 1.)\r\n    S = T.diag(sigma)\r\n    M = S.dot(C).dot(S)\r\n    return M\r\n\r\n\r\ndef analyze_robust(data):\r\n    with pm.Model() as model:\r\n        # priors\r\n        sigma = pm.Uniform('sigma', lower=0, upper=1000, shape=2,\r\n                           testval=[0.133434, 5.9304],  # init with mad\r\n                           transform=None)\r\n        rho = pm.Uniform('r', lower=-1, upper=1,\r\n                         testval=-0.2144021,  # init with Spearman's correlation\r\n                         transform=None)\r\n\r\n        # print values for debugging\r\n        rho_p = Print('rho')(rho)\r\n        sigma_p = Print('sigma')(sigma)\r\n\r\n        cov = pm.Deterministic('cov', covariance(sigma_p, rho_p))\r\n        cov_p = Print('cov')(cov)\r\n\r\n        mult_norm = pm.MvNormal('mult_norm', mu=[10.1, 79.],  # set mu to median\r\n                                cov=cov_p, observed=data.T)\r\n\r\n    return model\r\n\r\n\r\ns = [9.92, 94, 9.94, 79, 9.97, 78, 9.93, 83, 9.90, 77, 9.93, 76, 10.00, 74, 9.97, 87, 10.00, 86, 10.01,\r\n     83, 10.08, 75, 10.09, 74, 10.15, 92, 10.15, 69, 10.17, 79, 10.17, 71, 10.19, 80, 10.30, 80, 10.31,\r\n     77, 10.34, 87]\r\nx = [float(a) for i, a in enumerate(s) if i % 2 == 0]\r\ny = [float(a) for i, a in enumerate(s) if i % 2 != 0]\r\ndata = np.array([x, y])\r\n\r\nx.append(9.5)\r\ny.append(115.)\r\n\r\nrobust_model = analyze_robust(data)\r\nwith robust_model:\r\n    robust_trace = pm.sample(500, tune=250, step=pm.Metropolis(),\r\n                             random_seed=21412, progressbar=False)\r\n```\r\n\r\nNotice that the prior on `rho` is uniform between -1 and 1.\r\n\r\n## Result\r\n\r\nRunning the code with Theano's [NanGuardMode](http://deeplearning.net/software/theano/library/compile/nanguardmode.html) enabled results in NaN after a couple of iterations:\r\n\r\n```\r\nrho __str__ = -0.2144021\r\nsigma __str__ = [ 0.133434  5.9304  ]\r\ncov __str__ = [[  1.78046324e-02  -1.69660025e-01]\r\n [ -1.69660025e-01   3.51696442e+01]]\r\nsigma __str__ = [ 0.133434  5.9304  ]\r\nrho __str__ = -0.2144021\r\nsigma __str__ = [ 0.133434  5.9304  ]\r\nrho __str__ = -0.2144021\r\nrho __str__ = -0.6753367200324893\r\ncov __str__ = [[  1.78046324e-02  -1.69660025e-01]\r\n [ -1.69660025e-01   3.51696442e+01]]\r\ncov __str__ = [[  1.78046324e-02  -5.34405423e-01]\r\n [ -5.34405423e-01   3.51696442e+01]]\r\nrho __str__ = -0.2144021\r\nsigma __str__ = [ 0.133434  5.9304  ]\r\nsigma __str__ = [-0.08992942  6.58322907]\r\ncov __str__ = [[  1.78046324e-02  -1.69660025e-01]\r\n [ -1.69660025e-01   3.51696442e+01]]\r\ncov __str__ = [[  8.08730093e-03   1.26931614e-01]\r\n [  1.26931614e-01   4.33389049e+01]]\r\nsigma __str__ = [ 0.133434  5.9304  ]\r\nrho __str__ = -0.2144021\r\nsigma __str__ = [ 0.133434  5.9304  ]\r\nrho __str__ = -0.2144021\r\nrho __str__ = -1.3866936871442477\r\ncov __str__ = [[  1.78046324e-02  -1.69660025e-01]\r\n [ -1.69660025e-01   3.51696442e+01]]\r\ncov __str__ = [[  1.78046324e-02  -1.09731428e+00]\r\n [ -1.09731428e+00   3.51696442e+01]]\r\nTraceback (most recent call last):\r\n  File \"/home/sebp/miniconda3/envs/pymc3/lib/python3.5/site-packages/theano/gof/link.py\", line 1012, in f\r\n    wrapper(i, node, *thunks)\r\n  File \"/home/sebp/miniconda3/envs/pymc3/lib/python3.5/site-packages/theano/compile/nanguardmode.py\", line 307, in nan_check\r\n    do_check_on(x[0], node, fn, False)\r\n  File \"/home/sebp/miniconda3/envs/pymc3/lib/python3.5/site-packages/theano/compile/nanguardmode.py\", line 272, in do_check_on\r\n    raise AssertionError(msg)\r\nAssertionError: NaN detected\r\nNanGuardMode found an error in the output of a node in this variable:\r\nElemwise{log,no_inplace} [id A] ''   \r\n |Elemwise{true_div,no_inplace} [id B] ''   \r\n   |TensorConstant{1.0} [id C]\r\n   |Det [id D] ''   \r\n     |MatrixInverse [id E] ''   \r\n       |Print{message='cov', attrs=('__str__',), global_fn=<function _print_fn at 0x7fa10a171ea0>} [id F] ''   \r\n         |Dot22 [id G] 'cov'   \r\n           |Dot22 [id H] ''   \r\n           | |Diag [id I] ''   \r\n           | | |Print{message='sigma', attrs=('__str__',), global_fn=<function _print_fn at 0x7fa10a171ea0>} [id J] ''   \r\n           | |   |sigma_shared [id K]\r\n           | |FillDiagonal [id L] ''   \r\n           |   |Alloc [id M] ''   \r\n           |   | |Print{message='rho', attrs=('__str__',), global_fn=<function _print_fn at 0x7fa10a171ea0>} [id N] ''   \r\n           |   | | |Subtensor{int64} [id O] ''   \r\n           |   | |   |Subtensor{int64:int64:} [id P] ''   \r\n           |   | |   | |inarray1 [id Q]\r\n           |   | |   | |Constant{0} [id R]\r\n           |   | |   | |Constant{1} [id S]\r\n           |   | |   |Constant{0} [id R]\r\n           |   | |TensorConstant{2} [id T]\r\n           |   | |TensorConstant{2} [id T]\r\n           |   |TensorConstant{1.0} [id U]\r\n           |Diag [id I] ''   \r\n\r\n\r\n\r\nDuring handling of the above exception, another exception occurred:\r\n\r\nTraceback (most recent call last):\r\n  File \"/home/sebp/Projekte/bayesian/robust_cor.py\", line 49, in <module>\r\n    robust_trace = pm.sample(500, tune=250, step=pm.Metropolis(), random_seed=21412, progressbar=False)\r\n  File \"/home/sebp/miniconda3/envs/pymc3/lib/python3.5/site-packages/pymc3/sampling.py\", line 157, in sample\r\n    return sample_func(**sample_args)\r\n  File \"/home/sebp/miniconda3/envs/pymc3/lib/python3.5/site-packages/pymc3/sampling.py\", line 167, in _sample\r\n    for strace in sampling:\r\n  File \"/home/sebp/miniconda3/envs/pymc3/lib/python3.5/site-packages/pymc3/sampling.py\", line 249, in _iter_sample\r\n    point = step.step(point)\r\n  File \"/home/sebp/miniconda3/envs/pymc3/lib/python3.5/site-packages/pymc3/step_methods/compound.py\", line 16, in step\r\n    point = method.step(point)\r\n  File \"/home/sebp/miniconda3/envs/pymc3/lib/python3.5/site-packages/pymc3/step_methods/arraystep.py\", line 142, in step\r\n    apoint = self.astep(bij.map(point))\r\n  File \"/home/sebp/miniconda3/envs/pymc3/lib/python3.5/site-packages/pymc3/step_methods/metropolis.py\", line 137, in astep\r\n    q_new = metrop_select(self.delta_logp(q, q0), q, q0)\r\n  File \"/home/sebp/miniconda3/envs/pymc3/lib/python3.5/site-packages/theano/compile/function_module.py\", line 859, in __call__\r\n    outputs = self.fn()\r\n  File \"/home/sebp/miniconda3/envs/pymc3/lib/python3.5/site-packages/theano/gof/link.py\", line 1014, in f\r\n    raise_with_op(node, *thunks)\r\n  File \"/home/sebp/miniconda3/envs/pymc3/lib/python3.5/site-packages/theano/gof/link.py\", line 314, in raise_with_op\r\n    reraise(exc_type, exc_value, exc_trace)\r\n  File \"/home/sebp/miniconda3/envs/pymc3/lib/python3.5/site-packages/six.py\", line 685, in reraise\r\n    raise value.with_traceback(tb)\r\n  File \"/home/sebp/miniconda3/envs/pymc3/lib/python3.5/site-packages/theano/gof/link.py\", line 1012, in f\r\n    wrapper(i, node, *thunks)\r\n  File \"/home/sebp/miniconda3/envs/pymc3/lib/python3.5/site-packages/theano/compile/nanguardmode.py\", line 307, in nan_check\r\n    do_check_on(x[0], node, fn, False)\r\n  File \"/home/sebp/miniconda3/envs/pymc3/lib/python3.5/site-packages/theano/compile/nanguardmode.py\", line 272, in do_check_on\r\n    raise AssertionError(msg)\r\nAssertionError: NaN detected\r\nNanGuardMode found an error in the output of a node in this variable:\r\nElemwise{log,no_inplace} [id A] ''   \r\n |Elemwise{true_div,no_inplace} [id B] ''   \r\n   |TensorConstant{1.0} [id C]\r\n   |Det [id D] ''   \r\n     |MatrixInverse [id E] ''   \r\n       |Print{message='cov', attrs=('__str__',), global_fn=<function _print_fn at 0x7fa10a171ea0>} [id F] ''   \r\n         |Dot22 [id G] 'cov'   \r\n           |Dot22 [id H] ''   \r\n           | |Diag [id I] ''   \r\n           | | |Print{message='sigma', attrs=('__str__',), global_fn=<function _print_fn at 0x7fa10a171ea0>} [id J] ''   \r\n           | |   |sigma_shared [id K]\r\n           | |FillDiagonal [id L] ''   \r\n           |   |Alloc [id M] ''   \r\n           |   | |Print{message='rho', attrs=('__str__',), global_fn=<function _print_fn at 0x7fa10a171ea0>} [id N] ''   \r\n           |   | | |Subtensor{int64} [id O] ''   \r\n           |   | |   |Subtensor{int64:int64:} [id P] ''   \r\n           |   | |   | |inarray1 [id Q]\r\n           |   | |   | |Constant{0} [id R]\r\n           |   | |   | |Constant{1} [id S]\r\n           |   | |   |Constant{0} [id R]\r\n           |   | |TensorConstant{2} [id T]\r\n           |   | |TensorConstant{2} [id T]\r\n           |   |TensorConstant{1.0} [id U]\r\n           |Diag [id I] ''   \r\n\r\n\r\nApply node that caused the error: Elemwise{log,no_inplace}(Elemwise{true_div,no_inplace}.0)\r\nToposort index: 58\r\nInputs types: [TensorType(float64, scalar)]\r\nInputs shapes: [()]\r\nInputs strides: [()]\r\nInputs values: [array(-0.5779160437553532)]\r\nOutputs clients: [[Elemwise{add,no_inplace}(Elemwise{mul,no_inplace}.0, Elemwise{log,no_inplace}.0)]]\r\n\r\nBacktrace when the node is created(use Theano flag traceback.limit=N to make it longer):\r\n  File \"/home/sebp/Projekte/bayesian/robust_cor.py\", line 47, in <module>\r\n    robust_model = analyze_robust(data)\r\n  File \"/home/sebp/Projekte/bayesian/robust_cor.py\", line 32, in analyze_robust\r\n    cov=cov_p, observed=data.T)\r\n  File \"/home/sebp/miniconda3/envs/pymc3/lib/python3.5/site-packages/pymc3/distributions/distribution.py\", line 27, in __new__\r\n    return model.Var(name, dist, data)\r\n  File \"/home/sebp/miniconda3/envs/pymc3/lib/python3.5/site-packages/pymc3/model.py\", line 309, in Var\r\n    distribution=dist, model=self)\r\n  File \"/home/sebp/miniconda3/envs/pymc3/lib/python3.5/site-packages/pymc3/model.py\", line 592, in __init__\r\n    self.logp_elemwiset = distribution.logp(data)\r\n  File \"/home/sebp/miniconda3/envs/pymc3/lib/python3.5/site-packages/pymc3/distributions/multivariate.py\", line 117, in logp\r\n    result = k * tt.log(2 * np.pi) + tt.log(1. / det(tau))\r\n\r\nDebugprint of the apply node: \r\nElemwise{log,no_inplace} [id A] <TensorType(float64, scalar)> ''   \r\n |Elemwise{true_div,no_inplace} [id B] <TensorType(float64, scalar)> ''   \r\n   |TensorConstant{1.0} [id C] <TensorType(float32, scalar)>\r\n   |Det [id D] <TensorType(float64, scalar)> ''   \r\n     |MatrixInverse [id E] <TensorType(float64, matrix)> ''   \r\n       |Print{message='cov', attrs=('__str__',), global_fn=<function _print_fn at 0x7fa10a171ea0>} [id F] <TensorType(float64, matrix)> ''   \r\n         |Dot22 [id G] <TensorType(float64, matrix)> 'cov'   \r\n```\r\n\r\nEventually, `rho` becomes less than -1 (**rho __str__ = -1.3866936871442477**) and the computed covariance matrix is not positive definite anymore, thus its determinate is negative, which causes log to produce NaN.\r\n\r\nGiven a uniform prior between -1 and 1, I would expect that a value smaller -1 would be impossible.\r\n\r\n## Enabling NanGuardMode\r\n\r\nYou need to modify the `theano.function` call in [pymc3.step_methods.metropolis.delta_logp](pymc3/step_methods/metropolis.py) to use Theano's `NanGuardMode` class:\r\n```python\r\nfrom theano.compile.nanguardmode import NanGuardMode\r\nf = theano.function([inarray1, inarray0], logp1 - logp0,\r\n                    mode=NanGuardMode(True, False, False))  # only raise error on NaN, not Inf or big value\r\n```\r\n\r\n\r\n## Setup\r\nAnaconda on Fedora 24 with\r\n\r\n- Python 3.5.2\r\n- pymc3 3.0.rc2\r\n- numpy: 1.11.2\r\n- scipy: 0.18.1\r\n- theano 0.8.2\r\n"},"comment":{"url":"https://api.github.com/repos/pymc-devs/pymc3/issues/comments/264192597","html_url":"https://github.com/pymc-devs/pymc3/issues/1501#issuecomment-264192597","issue_url":"https://api.github.com/repos/pymc-devs/pymc3/issues/1501","id":264192597,"user":{"login":"sebp","id":74460,"avatar_url":"https://avatars.githubusercontent.com/u/74460?v=3","gravatar_id":"","url":"https://api.github.com/users/sebp","html_url":"https://github.com/sebp","followers_url":"https://api.github.com/users/sebp/followers","following_url":"https://api.github.com/users/sebp/following{/other_user}","gists_url":"https://api.github.com/users/sebp/gists{/gist_id}","starred_url":"https://api.github.com/users/sebp/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/sebp/subscriptions","organizations_url":"https://api.github.com/users/sebp/orgs","repos_url":"https://api.github.com/users/sebp/repos","events_url":"https://api.github.com/users/sebp/events{/privacy}","received_events_url":"https://api.github.com/users/sebp/received_events","type":"User","site_admin":false},"created_at":"2016-12-01T14:51:57Z","updated_at":"2016-12-01T14:51:57Z","body":"Following @fonnesbeck's [suggestion](https://github.com/pymc-devs/pymc3/issues/1501#issuecomment-261326343) is probably the way to go, add a class similar to NanGuardMode that warns about NaN values, but does not exit the application.\r\n\r\nRegarding MultivariateT, I tried what you suggested by fixing `nu` to  100 and running the sampler for 50000 iterations. I don't get an error (or NaN values), but the posteriors for `sigma` and `r` cover almost the complete range. From what I can tell, the implementation of the log-likelihood of the MultivariateT distribution looks okay."}}
{"actor":{"display_login":"sebp"}}