Sample program

stan code

functions {
  vector fun(vector x, real a, real b, real c) {
    return a * square(x) + b * x + c;
  }
}
data {
  int<lower=0> N;
  vector[N] y;
  vector[N] x;
}
parameters {
  real b;
  real<lower=0> c;
  real<lower=0> sigma;
}
transformed parameters {
  vector[N] mu;
  real a = 0;
  mu = fun(x, a, b, c);
}
model {
  a ~ uniform(-5, 5);
  b ~ uniform(-5, 5);
  c ~ uniform(-30, 30);
  //sigma ~ exponential(100);
  sigma ~ normal(0, 1); // half normal
  y ~ normal(mu, sigma);
}
generated quantities {
  real xpred = uniform_rng(-5, 5);
  real ypred = normal_rng(fun([xpred]', a, b, c)[1], sigma);
}

data:

{
  "N": 4,
  "x": [
    9,
    11,
    16,
    18
 ],
  "y": [
    20,
    22,
    50,
    52
  ]
}