Hi all,
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. And I assume that the location parameter of GEV distribution and the shape parameter of Gamma distribution have a linear relationship with time. The code and data are attached.
After running the Rstan model, I get 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 (-5, 5) failed after 100 attempts.
Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.”
I think the problem maybe the priors, but the initial values are only rarely necessary. Is there something wrong with what I am trying to do? Could you guys help me have a look? Thanks a lot!
> functions{
> real gev_lpdf(vector xx, vector mu, real sigma, real xi) {
> vector[rows(xx)] z;
> vector[rows(xx)] lp;
> int N;
> N = rows(xx);
> for(n in 1:N){
> z[n] = 1 + (xx[n] - mu[n]) * xi / sigma;
> lp[n] = (1 + (1 / xi)) * log(z[n]) + pow(z[n], -1/xi);
> }
> return -sum(lp) - N * log(sigma);
> }
>
> real gev_cdf(vector xx, vector mu, real sigma, real xi){ //fl is streamflow data series
> vector[rows(xx)] z;
> vector[rows(xx)] lp;
> int N;
> N = rows(xx);
> for(n in 1:N){
> z[n] = 1 + (xx[n] - mu[n]) * xi / sigma;
> lp[n]= exp(-pow(z[n],-1/xi));
> }
> return sum(lp);
> }
> 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));
> }
> real gumbel_copula_cdf(real u, real v, real theta){
> real neg_log_u = -log(u);
> real neg_log_v = -log(v);
> real temp = neg_log_u ^ theta + neg_log_v ^ theta;
> return -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 <lower=0>sigma;
> real xi;
> // real<lower=0> alpha;
> real<lower=0> beta;
> real mu0;
> real beta_mu;
> real alpha0;
> real beta_alpha;
> real<lower=0,upper=1> tau;
> }
> transformed parameters {
> real<lower=1> theta = 1 / (1 - tau);
> vector<lower=0>[N] alpha;
> vector[N] mu;
> for (n in 1:N)
> mu[n] = mu0 + beta_mu * n;
>
> for (n in 1:N)
> alpha[n] = alpha0 + beta_alpha * n;
> }
> model {
>
> vector[N] z;
> // model
>
> x~gamma(alpha,beta);
> y ~ gev(mu,sigma,xi);
>
> for (n in 1:N)
> z[n] = gumbel_copula(gamma_cdf(x[n], alpha, beta),
> gev_cdf(y, mu, sigma,xi), theta);
> target += z;
> // priors
> // mu ~ normal(1000,100);
> mu0~gamma(0.001,0.001);
> beta~gamma(0.001,0.001);
> alpha0~gamma(0.001,0.001);
> beta_alpha~gamma(0.001,0.001);
> sigma ~ normal(750,100);
> xi~ uniform(-1,1);
> alpha ~ normal(5,5);
> beta ~uniform(0,0.1);
> }
then I run it:
> library(rstan)
dat ← read.csv(“C:/Users/Administrator/Desktop/Rstan/data.csv”)
dat$N ← NULL
initf<-function()list(mu=1000,sigma=750,xi=0.5,alpha=5,beta=0.05)
post ← stan(“F:/paper_Copula/R_code/stan/gumbel_nonstationary.stan”, data = list(N = nrow(dat), x= dat$x, y = dat$y),init = initf,init_r=5,iter=12000,warmup=2000,thin=10,chains=3)
print(post)
data.csv (1.0 KB)