[Please include Stan program and accompanying data if possible]
Hello, everyone,
I am a beginner for Rstan, right now I encounter a problem regarding the Copula Bayesian estimation using Stan.
I try to estimate the joint probability of rainfall and streamflow and the parameter of Gumbel Copula. I assume the rainfall data comes from a Gamma distribution and the streamflow data comes from a GEV distribution. I found the Gumbel copula and GEV distribution are not available by default in Rstan, so I define the functions first.
functions{
real gev_pdf_log(real xx, real mu, real sigma, real xi){
real ystar;
real xy_p1;
ystar=(xx-mu)/sigma;
xy_p1 = xi* ystar +1;
return -log(sigma)-(1+1/xi)*log(xy_p1)-(xy_p1)^(-1/xi);
}
real gev_cdf(real xx, real mu, real sigma, real xi){
real ystar;
real xy_p1;
ystar=(xx-mu)/sigma;
xy_p1 = xi* ystar +1;
return exp(-(xy_p1)^(-1/xi));
}
real gumbel_copula(real u, real v, real theta) {
real neg_log_u = -log(u);
real neg_log_v = -log(v);
real theta_m1 = theta - 1;
real temp = neg_log_u ^ theta + neg_log_v ^ theta;
return theta_m1 * log(neg_log_u) + theta_m1 * log(neg_log_v)
+ neg_log_u + neg_log_v - pow(temp, inv(theta))
+ (-2+2/theta)*log(temp)
+ log1p(theta_m1*pow(temp,-1/theta));
}
}
data {
int<lower=0> N; // number of years
vector[N] x; // rainfall data series
vector[N] y; // streamflow data series
}
parameters {
real mu;
real <lower=0>sigma;
real xi;
real<lower=0> alpha;
real<lower=0> beta;
real<lower=0,upper=1> tau; // kendall’s tau
}
transformed parameters {
real<lower=1> theta = 1 / (1 - tau);
}
model {
vector[N] summands;
// priors
mu ~ normal(1000,100);
sigma ~ normal(750,100);
xi ~ uniform(-1,1);
alpha ~ normal(5,5);
beta ~ uniform(0,0.1);
//tau ~ gamma(0.01,0.01);
//tau ~ cauchy(0, 2.5);
// model
x ~ gamma(alpha,beta);
for (n in 1:N)
y[n] ~ gev_pdf(mu,sigma,xi);
for (n in 1:N)
summands[n] = gumbel_copula(gamma_cdf(x[n], alpha, beta),
gev_cdf(y[n], mu, sigma,xi), theta);
target += summands;
}
then I run it:
library(rstan)
dat <- read.csv("C:/Users/Administrator/Desktop/Rstan/downloaded/data.csv")
dat$N <- NULL
post <- stan("C:/Users/Administrator/Desktop/Rstan/downloaded/gumbel.stan", data = list(N = nrow(dat), x= dat$x, y = dat$y),iter=12000,warmup=2000,thin=10,chains=3)
print(post)
I got the following error:
Rejecting initial value:
Log probability evaluates to log(0), i.e. negative infinity.
Stan can't start sampling from this initial value.
Initialization between (-2, 2) failed after 100 attempts.
Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.
The priors I set for parameters first are non-informative distributions, but not work. Then I tried to use the MLE method to estimate the parameters of each distribution and set the prior distributions based on the MLE results which are shown in the Code. I don’t know the reason and have been confused by days`.
the Code and data are attached.Could you help me have a look? Thanks a lot!
data.csv (583 Bytes)