// poisson_regression.stan

data {
    int<lower=0> n; // number of observations
    int<lower=0> p; // number of covariates
    int<lower=0> m; // number of grid points
    array[n] int<lower=0> y; // response variables
    matrix[n,p] X; // matrix of covariates
    matrix[m,p] grid; // matrix of grid points
    real<lower=0> a;
    real<lower=0> b;
}
transformed data {
    real t_c = (2*a+p-1)/(2*b);
}
parameters {
    vector[p] beta;
}
model {
    sqrt(dot_self(beta)*t_c) ~ student_t(2*a+p-1, 0, 1);
    target += poisson_log_glm_lpmf( y | X, 0, beta );
}
generated quantities {
    vector[m] fn_vals;
    for (i in 1:m)
      fn_vals[i] = exp( dot_product(beta,grid[i]) );
}