fixed_jacobian

functions {
    real inv_transform(real x_new, real center, real sigma) {
        return erf(x_new - center) / sigma + x_new;
    }
    real jacobian(real x_new, real center, real sigma) {
        return log_sum_exp((
            log(2 / sqrt(pi()))
            - (x-c)**2
            - log(sigma)
            ), 1
        );
    }
}
data {
    int N;
    vector[N] x;

    real center;
    real sigma;
}
parameters {
    real transformed_mu;
    real<lower=0> sd:
}
transformed parameters {
    real mu = inv_transform(transformed_mu, center, sigma);
}
model {
    mu ~ normal(0, 1);
    target += jacobian(transformed_mu, center, sigma);
    sd ~ normal(0, 1);
    x ~ normal(mu, sd);
}