@StaffanBetner I was playing around with the continuous bernoulli and if p is a vector I needed to divide the normalizing constant by the sample size.
The stan code. Note that in Stan 2.29 you can declare the same function name with different inputs, one is a vector of p and the other is a real p. Does this seem right?
real continuous_bernoulli_lpdf(vector y, vector lambda) {
int N = num_elements(y);
real lp = N * log2();
real lp_c = 0;
int counter = 0;
for (n in 1:N) {
lp += y[n] * log(lambda[n]) + (1 - y[n]) * log1m(lambda[n]);
if (lambda[n] != 0.5) {
counter += 1;
lp_c += log( atanh(1 - 2 * lambda[n]) / (1 - 2 * lambda[n]) );
}
}
return lp + lp_c / counter;
}
real continuous_bernoulli_lpdf(vector y, real lambda) {
int N = num_elements(y);
real lp = N * log2() + sum(y * log(lambda) + (1 - y) * log1m(lambda));
if (lambda != 0.5)
lp += log( atanh(1 - 2 * lambda) / (1 - 2 * lambda) );
return lp;
}