# GFI in STAN

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!

You could calculate it in the generated quantities block with

generated quantities {
int I = mu <= 0.1;
}


in which case the average of that over the MCMC draws is an estimate of the sought probability. Something that simple can also easily be calculated in R afterward after extracting the draws of mu.

Hi all,

I want to separate my R code and STAN code but unfortunately am getting a “PARSER EXPECTED: whitespace to end of file. FOUND AT line 1:” error.

My R code looks like this

library(rstan)
library(rstanarm)
library("car")

set.seed(1)

vn<-c(5,25,125)
vmu<-c(0.1,1,5)

coverage=matrix(rep(NA,9),3,3)
expectedlength=matrix(rep(NA,9),3,3)

vpval=numeric(0)

for(i in 1:3){
n=vn[i]
for(j in 1:3){
mu0=vmu[j]
myoutput = foreach(i=1:1000,.combine=cbind)%dopar%{
x<-rnorm(n,mu0,mu0^(0.75))
meanx=mean(x)
varx=var(x)

m.stan <- stan(file = "gfi.stan",iter = 1000, chains=1)
#data from stan processed
fitted_model = stan_model(model_code = gfi_model)

f_g = sampling(fitted_model, data = ds)

plot(f_g,show_density=TRUE)

summary_stats<-summary(f_g,pars="mu",probs=c(.05,.95))$summary summary_stats list_of_draws <- extract(f_g) list_of_draws$mu

if (list_of_draws$mu>summary_stats[4] && list_of_draws$mu<summary_stats[5]){
indicator<-1
}else{
indicator <-0
}

CI_length<-summary_stats[5]-summary_stats[4]

Exp_.1<-abs(summary_stats[1]-j)
plot(summary_stats[4],pars=c("mu"))
qqPlot(list_of_draws\$mu)

}
vpval=rbind(vpval,output[1,])
coverage[i,j]=mean(output[3,]) #coverage
expectedleangth[i,j]=mean(output[4,]) #expected length
}
}

x_8<-rnorm(vn[3],vmu[3],vmu[3]^(0.75))
new_x<-cbind(x,x_1,x_2)
new_x1<-cbind(x_3,x_4,x_5)
new_x2<-cbind(x_6,x_7,x_8)
output<-rbind(new_x,new_x1,new_x2)

ds <- list(Y=x, N = length(x))



and my STAN code looks like this

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);
}

"


Am not sure what is going wrong. I have tried adding a newline character to my STAN code but that is not resolving the issue.

That Stan code is valid, although it is easier to keep it in an external file. I think you need to upgrade your Stan.

I have upgraded rstan and rstanarm to 2.18.2 but that still doesn’t seem to fix the issue. Additionally, What is the cleanest function to call a STAN model if I wish to run a STAN model on data I have generated in R?