Newbie question: how to do Bayesian parameter estimation with RStan?

I have a function

likelihood <- function(theta) { # theta is a list of N real numbers 
    # ... 
    # some black box here that depends on theta 
    # ... 

I also priors (uniform and normal) on the different numbers in theta.

I want to estimate the posterior distribution of theta according to my likelihood function and the priors that I have set. I am under the impression that I can do this with RStan, although from reading the manual and trying out several tutorials, it is very unclear to me how.

I’m pretty new though and I’m sure this is a pretty naive question, but could someone please explain or point me to a relevant tutorial?

(For what it’s worth, I consider this to be pretty easily done with python’s emcee package)

This (usually) isn’t too difficult, but you need to get familiar with the Stan language and the approach to estimation that is facilitated by Stan. Chapters 22 and 23 of the Stan manual discuss custom probability functions. Basically, you need to translate this R function to a Stan function. Then use rstan::expose_stan_function so that you can call it and make sure it matches your R implementation. You don’t have to combine everything into a single parameter vector theta like you would if you were calling optim; you can have multiple arguments to your function.

For the parameters that have priors, you just have to do

target += normal_lpdf(that_parameters | mean, sd);

For the parameters that have uniform priors (which is usually a bad idea), you don’t have to do anything to alter the posterior kernel (in log units).

1 Like

You will have to initialise your thetas, give them a prior and then express your likelihood in terms of distribution around equations/functional forms

real theta[N];
real<lower=0> sigma;
theta ~ normal(0, 1);
y ~ normal(theta[1] … theta[N], sigma); // likelihood

If you want the actual LL check this