GFI in STAN

#1

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!

0 Likes

#2

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.

0 Likes

#3

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)

vmad=numeric(0)
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)
      
      c(pvalue,mad,covered,length,meanx,varx)
    }
    vpval=rbind(vpval,output[1,])
    mad=rbinnd(mad,output[2,])
    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.

0 Likes

#4

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

0 Likes

#5

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?

0 Likes