Hi,
I’m currently implementing Generalized fiducial inference in STAN and have the following code.
library(rstan)
library(rstanarm)
library(MASS)
set.seed(1)
x <- rnorm(100, 0.1,0.1^(0.75))
x
#x <- c(1,2,3,2,7,2,4,8,24,78,6,24,21,2,3,7,3,58,246)
ds <- list(Y=x, N = length(x))
gfi_model = "
functions{
real genfid_log(vector t,real mu){
vector[num_elements(t)] jac;
vector[num_elements(t)] log_likl;
for(i in 1:num_elements(t)){
jac[i]=fabs(1+(1.5)*((t[i]-mu)/mu));
log_likl[i]=(t[i]-mu)^2/(2*mu^(1.5))+(0.75)*log(mu);
}
return log(sum(jac))-sum(log_likl);
}
}
data {
int<lower=1> N; // number of observations
vector [N] Y;
}
parameters {
real<lower=0> mu;
}
model {
Y~genfid(mu);
}
generated quantities{
}
"
fitted_model = stan_model(model_code = gfi_model)
f_g = sampling(fitted_model, data = ds)
plot(f_g,show_density=TRUE)
plot(f_g,pars=c("mu"))
I want to calculate the probability that
p(\mu <= 0.1)
by using an indicator random variable. My hunch is that I would put this in the generated quantities section, but I am not quite sure how to notate it in STAN.
Any help would be appreciated!