Sampling from Gaussian distribution by using stan

Modeling
1

Here, I would like to sample from a univariate normal distribution with mean mu and std sigma:
I know that it is simple to sample by simply using y ~ normal(mu,sigma).

But, as an exercise, I only want to use the proportional part of the Guassian distribution (kernel), and check whether stan works or not.

I made code as bellow, but it does not work…

Can I anyone help me to fix this???

functions {
  real gaussian(real x, real mu, real sigma) {
    return (1 / sigma) * exp( - (((x - mu)/sigma)^2) / 2 );
  }
}

data {
  real mu;      
  real<lower=0> sigma;
  real x;
}

parameters {
  real y;      
}

model {
  // Option 1
  // y ~ normal(mu, sigma);
  // Option 2
  target += gaussian(x, mu, sigma);
}

EDIT: @maxbiostat edited this post for syntax highlighting.

1 Like
2

I think you need to delete real x and use target += gaussian(y, mu, sigma). Finally, Stan works on the log scale of the pdf. So you need something like

return - log(sigma) - ((x - mu)/sigma)^2) / 2;

(Don’t trust my code!)

2 Likes
4

The example code ?rstan would help.

stanmodelcode <- "
data {
  int<lower=0> N;
  real y[N];
} 

parameters {
  real mu;
} 

model {
  target += normal_lpdf(mu | 0, 10);
  target += normal_lpdf(y  | mu, 1);
} 
"

y <- rnorm(20) 
dat <- list(N = 20, y = y); 
fit <- stan(model_code = stanmodelcode, model_name = "example", 
            data = dat, iter = 2012, chains = 3, verbose = TRUE,
            sample_file = file.path(tempdir(), 'norm.csv')) 
print(fit)

# extract samples 
e <- extract(fit, permuted = FALSE) # return a list of arrays 
str(e)
2 Likes