{"action":"created","issue":{"url":"https://api.github.com/repos/pymc-devs/pymc3/issues/593","labels_url":"https://api.github.com/repos/pymc-devs/pymc3/issues/593/labels{/name}","comments_url":"https://api.github.com/repos/pymc-devs/pymc3/issues/593/comments","events_url":"https://api.github.com/repos/pymc-devs/pymc3/issues/593/events","html_url":"https://github.com/pymc-devs/pymc3/issues/593","id":41182879,"number":593,"title":"helper class for Gaussians","user":{"login":"karlnapf","id":699071,"avatar_url":"https://avatars.githubusercontent.com/u/699071?v=3","gravatar_id":"","url":"https://api.github.com/users/karlnapf","html_url":"https://github.com/karlnapf","followers_url":"https://api.github.com/users/karlnapf/followers","following_url":"https://api.github.com/users/karlnapf/following{/other_user}","gists_url":"https://api.github.com/users/karlnapf/gists{/gist_id}","starred_url":"https://api.github.com/users/karlnapf/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/karlnapf/subscriptions","organizations_url":"https://api.github.com/users/karlnapf/orgs","repos_url":"https://api.github.com/users/karlnapf/repos","events_url":"https://api.github.com/users/karlnapf/events{/privacy}","received_events_url":"https://api.github.com/users/karlnapf/received_events","type":"User","site_admin":false},"labels":[],"state":"open","locked":false,"assignee":null,"milestone":null,"comments":4,"created_at":"2014-08-26T15:09:33Z","updated_at":"2015-04-30T16:29:24Z","closed_at":null,"body":"Hi there,\r\nmany samplers are based on Gaussian distributions, that need to be sampled from and whose log-pdf needs to be evaluated. If the Gaussian's covariance is not diagonal, this requires some computation. I think in the current pymc3 branch, things are sometimes re-computed (see #592), so what about adding a helper class for Gaussians that just decompose the covariance once, and then cheap to sample/log-pdf from?\r\n\r\nI use:\r\n```python\r\nclass Gaussian():\r\n    \"\"\"\r\n    Helper class to sample from and evaluate log-pdf of a multivariate Gaussian,\r\n    using efficient Cholesky based representation (Cholesky only computed once)\r\n    \"\"\"\r\n    def __init__(self, mu, Sigma, is_cholesky=False):\r\n        self.mu = mu\r\n        self.is_cholesky = is_cholesky\r\n        \r\n        if self.is_cholesky:\r\n            self.L = Sigma\r\n        else:\r\n            self.L = np.linalg.cholesky(Sigma)\r\n            \r\n        self.dimension = len(mu)\r\n    \r\n    def log_pdf(self, X):\r\n        # duck typing for shape\r\n        if len(np.shape(X)) == 1:\r\n            X = X.reshape(1, len(X))\r\n        \r\n        log_determinant_part = -sum(np.log(np.diag(self.L)))\r\n        \r\n        quadratic_parts = np.zeros(len(X))\r\n        for i in range(len(X)):\r\n            x = X[i] - self.mu\r\n            \r\n            # solve y=K^(-1)x = L^(-T)L^(-1)x\r\n            y = sp.linalg.solve_triangular(self.L, x.T, lower=True)\r\n            y = sp.linalg.solve_triangular(self.L.T, y, lower=False)\r\n            quadratic_parts[i] = -0.5 * x.dot(y)\r\n            \r\n        const_part = -0.5 * len(self.L) * np.log(2 * np.pi)\r\n        \r\n        return const_part + log_determinant_part + quadratic_parts\r\n    \r\n    def sample(self, n=1):\r\n        V = np.random.randn(self.dimension, n)\r\n\r\n        # map to our desired Gaussian and transpose to have row-wise vectors\r\n        return self.L.dot(V).T + self.mu\r\n```"},"comment":{"url":"https://api.github.com/repos/pymc-devs/pymc3/issues/comments/97869495","html_url":"https://github.com/pymc-devs/pymc3/issues/593#issuecomment-97869495","issue_url":"https://api.github.com/repos/pymc-devs/pymc3/issues/593","id":97869495,"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":"2015-04-30T16:29:23Z","updated_at":"2015-04-30T16:29:23Z","body":"Ah, I see. If I understand you correctly, your logpdf can't be used with `pymc3` yet because `theano` does not support triangular solve yet."}}