{"action":"created","issue":{"url":"https://api.github.com/repos/stan-dev/math/issues/191","labels_url":"https://api.github.com/repos/stan-dev/math/issues/191/labels{/name}","comments_url":"https://api.github.com/repos/stan-dev/math/issues/191/comments","events_url":"https://api.github.com/repos/stan-dev/math/issues/191/events","html_url":"https://github.com/stan-dev/math/issues/191","id":112684897,"number":191,"title":"Add squared exponential covariance function for gaussian processes","user":{"login":"syclik","id":425751,"avatar_url":"https://avatars.githubusercontent.com/u/425751?v=3","gravatar_id":"","url":"https://api.github.com/users/syclik","html_url":"https://github.com/syclik","followers_url":"https://api.github.com/users/syclik/followers","following_url":"https://api.github.com/users/syclik/following{/other_user}","gists_url":"https://api.github.com/users/syclik/gists{/gist_id}","starred_url":"https://api.github.com/users/syclik/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/syclik/subscriptions","organizations_url":"https://api.github.com/users/syclik/orgs","repos_url":"https://api.github.com/users/syclik/repos","events_url":"https://api.github.com/users/syclik/events{/privacy}","received_events_url":"https://api.github.com/users/syclik/received_events","type":"User","site_admin":false},"labels":[],"state":"open","locked":false,"assignee":null,"milestone":{"url":"https://api.github.com/repos/stan-dev/math/milestones/6","html_url":"https://github.com/stan-dev/math/milestones/v2.8.0++","labels_url":"https://api.github.com/repos/stan-dev/math/milestones/6/labels","id":1251012,"number":6,"title":"v2.8.0++","description":"Math library version after v2.8.0.","creator":{"login":"syclik","id":425751,"avatar_url":"https://avatars.githubusercontent.com/u/425751?v=3","gravatar_id":"","url":"https://api.github.com/users/syclik","html_url":"https://github.com/syclik","followers_url":"https://api.github.com/users/syclik/followers","following_url":"https://api.github.com/users/syclik/following{/other_user}","gists_url":"https://api.github.com/users/syclik/gists{/gist_id}","starred_url":"https://api.github.com/users/syclik/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/syclik/subscriptions","organizations_url":"https://api.github.com/users/syclik/orgs","repos_url":"https://api.github.com/users/syclik/repos","events_url":"https://api.github.com/users/syclik/events{/privacy}","received_events_url":"https://api.github.com/users/syclik/received_events","type":"User","site_admin":false},"open_issues":12,"closed_issues":5,"state":"open","created_at":"2015-08-12T02:36:48Z","updated_at":"2015-10-27T13:50:07Z","due_on":null,"closed_at":null},"comments":25,"created_at":"2015-10-21T21:38:49Z","updated_at":"2015-10-27T17:33:57Z","closed_at":null,"body":"## Feature Request\r\n\r\n__Version__\r\n\r\n2.8.0++\r\n\r\n__Description__\r\n\r\nAdd a function for the squared exponential covariance function to make GPs easier in Stan.\r\nThe covariance function:\r\n$$ k(x_i, x_j) = \\sigma^2 \\exp\\left(-\\frac{(x_i - x_j)^2}{2l^2}\\right) $$\r\n\r\ngpstuff calls this function `gpcf_sexp`.\r\ngpy calls this function `GPy.kern.src.rbf`\r\n\r\nWe should call this something reasonable.\r\n\r\n@avehtari + @cunni: are these usually parameterized by the variance (sigma^2) and the length-scale (l)?\r\n\r\n__Example__\r\n\r\nWe still need a reasonable name for the function, but the usage would look something like:\r\n```\r\ndata {\r\n  int N;\r\n  vector[N] x;\r\n}\r\nparameters {\r\n  real<lower = 0> variance;  // can't really use variance as a name here\r\n  real<lower = 0> lengthscale; \r\n}\r\ntransformed parameters {\r\n  matrix[N, N] cov;\r\n  cov <- squared_exponential(X, variance, lengthscale);\r\n}\r\n...\r\n```\r\n\r\nThat look about right?\r\n\r\n__Additional Information__\r\n\r\nWhat should we name this function?"},"comment":{"url":"https://api.github.com/repos/stan-dev/math/issues/comments/151580472","html_url":"https://github.com/stan-dev/math/issues/191#issuecomment-151580472","issue_url":"https://api.github.com/repos/stan-dev/math/issues/191","id":151580472,"user":{"login":"bob-carpenter","id":3383807,"avatar_url":"https://avatars.githubusercontent.com/u/3383807?v=3","gravatar_id":"","url":"https://api.github.com/users/bob-carpenter","html_url":"https://github.com/bob-carpenter","followers_url":"https://api.github.com/users/bob-carpenter/followers","following_url":"https://api.github.com/users/bob-carpenter/following{/other_user}","gists_url":"https://api.github.com/users/bob-carpenter/gists{/gist_id}","starred_url":"https://api.github.com/users/bob-carpenter/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/bob-carpenter/subscriptions","organizations_url":"https://api.github.com/users/bob-carpenter/orgs","repos_url":"https://api.github.com/users/bob-carpenter/repos","events_url":"https://api.github.com/users/bob-carpenter/events{/privacy}","received_events_url":"https://api.github.com/users/bob-carpenter/received_events","type":"User","site_admin":false},"created_at":"2015-10-27T17:33:57Z","updated_at":"2015-10-27T17:33:57Z","body":"Thanks.  The birthday problem's a bit simpler than many we'll get\nbecause there's only a single length scale.   Below is what I came\nup with using the full matrix approach, though I'm sure I botched the\ngenerated quantities --- I still don't really understand GPs very well,\nbut am too lazy to figure out what this really needs to be.  Don't we\nneed a covariance matrix involving both x and x_tilde?  And then we\nonly want to generate y_tilde?  So what I wrote can't be right.  But\nthen neither can Michael's because he's redeclaring the data variable y\nin a way I couldn't figure out.\n\n- Bob\n\nP.S.  Using four binary sums leads to four expressions per entry,\nwhereas it'd be better to just have a sum().  Michael's approach has\nthe same issue with sum not being efficient as a chain compared to\nthe vectorized sum because of the additional subexpressions generated.\nOn the other hand, given everything else, it may not even be noticeable.\n\n\nfunctions {\n  matrix cov_fun(vector[] x, real sigma, real l) {\n    return cov_sq_exp(x, l[1], sigma[1])\n      + cov_sq_exp(x, l[2], sigma[2]) \n      + cov_sin_exp(x, pi() / 7, l[3]) .* cov_sq_exp(x1, x2, l[4], sigma[3])\n      + cov_sin_exp(x, pi() / 365, l[5]) .* cov_sq_exp(x, l[6], sigma[4])\n      + cov_delta(x);\n  }\n}\ndata {\n  int<lower=1> n_train;\n  int<lower=1> n_dim; // for the birthday problem n_dim = 1\n  matrix[n_data, n_dim] x; // this corresponds to time for the birthday problem\n  vector[n_train] y; // this corresponds to births\n  int<lower=1> n_test;\n  matrix[n_test, n_dim] x_tilde;\n}\nparameters {\n  vector sigma[4];\n  vector l[6];\n}\nmodel {\n  theta ~ cauchy(0, 2.5);\n  y ~ multi_normal(rep_vector(0, n_train), cov_fun(x, sigma, l));\n}\ngenerated quantities {\n  real y_tilde[n_test];\n  y_tilde <- multi_normal_rng(rep_vector(0, n_train + n_test),\n                              cov_fun(append_rows(x, x_tilde), sigma, l));\n}\n\n\n> On Oct 27, 2015, at 1:05 PM, Michael Betancourt <notifications@github.com> wrote:\n> \n> Here’s my proposal as for the birthday problem as discussed at the meeting.\n> Note that I left out the 5th, 6th, and 7th kernels as they do not seem to induce\n> a positive-definite covariance so I must not be understanding correctly.\n> \n> I didn’t bother to be careful about length scales and noises and their squares.\n> \n> function {\n> gp_cov(real x1[], real x2[], real theta[]) {\n> real sigma1;\n> real l1;\n> real sigma2;\n> real l2;\n> real sigma3;\n> real l31;\n> real l32;\n> real sigma4;\n> real l41;\n> real l42;\n> real sigma5;\n> real l5;\n> real sigma6;\n> real l6;\n> real epsilon;\n> \n> // My kingdom for a declare and initialize outside of the parameters block!\n> sigma1 <- theta[1];\n> l1 <- theta[2];\n> sigma2 <- theta[3];\n> l2 <- theta[4];\n> sigma3 <- theta[5];\n> l31 <- theta[6];\n> l32 <- theta[7];\n> sigma4 <- theta[8];\n> l41 <- theta[9];\n> l42 <- theta[10];\n> sigma5 <- theta[11];\n> l5 <- theta[12];\n> sigma6 <- theta[13];\n> l6 <- theta[14];\n> epsilon <- theta[15];\n> \n> return sigma1 * cov_sq_exp(x1, x2, l1)\n> + sigma2 * cov_sq_exp(x1, x2, l2) \n> + sigma3 * cov_sin_exp(x1, x2, pi() / 7, l31) * cov_sq_exp(x1, x2, l32)\n> + sigma4 * cov_sin_exp(x1, x2, pi() / 365, l41) * cov_sq_exp(x1, x2, l42)\n> + epsilon * cov_delta(x1, x2);\n> }\n> }\n> \n> data {\n> int<lower=1> n_train;\n> int<lower=1> n_dim; // for the birthday problem n_dim = 1\n> \n> matrix[n_data, n_dim] x; // this corresponds to time for the birthday problem\n> vector[n_train] y; // this corresponds to births\n> \n> int<lower=1> n_test;\n> matrix[n_test, n_dim] x_tilde;\n> }\n> \n> transformed data {\n> int n_hyper;\n> n_hyper=15;\n> }\n> \n> parameters {\n> real<lower=0> theta[n_hyper];\n> }\n> \n> model {\n> theta ~ cauchy(0, 2.5);\n> y ~ gp(gp_cov, x, theta);\n> }\n> \n> generated quantities {\n> real y[n_test];\n> y ~ gp_rng(gp_cov, x_test, theta);\n> }\n> \n> On Oct 27, 2015, at 9:15 AM, Michael Betancourt <betanalpha@gmail.com> wrote:\n> \n> > \n> >> 1. all cov functions should start cov_ . See\n> >> http://www.gaussianprocess.org/gpml/code/matlab/cov/ for a list of\n> >> covariance functions. One can not really hope to cover the space, so as\n> >> extensible as possible is a sensible goal.\n> > \n> > Covariance functions are going to be a rabbit hole because we\n> > also have to consider their algebraic structure — i.e., we can\n> > add and multiply covariance functions to generate new covariance\n> > functions. Given the utility of adding a bunch of covariance functions\n> > together, this isn’t a pedantic point. The problem is in implementation,\n> > and I’m not sure how we’d do it outside of expression templates or\n> > something of that sort.\n> > \n> >> 2. kernels will likely be desired much more often than gp (kernel methods\n> >> are a big class), so I favor accessible kernel functions, not just inside\n> >> the GP likelihood.\n> > \n> > But the Stan modeling language is probabilistic in nature, and the only\n> > decent measures floating around are GPs. So what kernel method is\n> > not going to end with a multivariate gaussian or a related marginal or \n> > conditional?\n> > \n> \n> —\n> Reply to this email directly or view it on GitHub.\n> \n\n"}}