I have following code which give me error message:
Chain 1: Log probability evaluates to log(0), i.e. negative infinity.
Chain 1: Stan can’t start sampling from this initial value.
Chain 1:
Chain 1: Initialization between (-2, 2) failed after 100 attempts.
Chain 1: Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.
[1] “Error in sampler$call_sampler(args_list[[i]]) : Initialization failed.”
[1] “error occurred during calling the sampler; sampling not done”
The R-code is:
dat=read.table(“Loss Reserve Data.txt”); st=1;
N=dim(dat)[1] #number of rows
colnames(dat)=c(“Year”,“Y”,“poli”,“develop”)
y0=0
y=dat[,2]
burn=1000
burn1=burn+1
iters=10000
dat$poli<-factor(dat$poli)
dat$develop<-factor(dat$develop)
poli=as.integer(dat$poli) # rows i effect (policy year)
develop=as.integer(dat$develop) # columns j effect (development year)
I=nlevels(dat$poli);I1=I-1; #number of levels of policy (i)
J=nlevels(dat$develop);J1=J-1 #number of levels of development (j)
P=I1+6 #no of par in the model
#begin Rstan program
trialmodelA<- "
//Generalized Beta 2 Distribution
functions{
real beta2_log(real x, real a, real mu, real p, real q) {
return log(fabs(a))+(a*p-1)log(x)+lgamma(p+q)+(ap)lgamma(p+(1/a))+(ap)lgamma(q-(1/a))-(ap)log(mu)-(ap)lgamma§-(ap)*lgamma(q)-lgamma§-lgamma(q)-(p+q)log(1+((xtgamma(p+(1/a))tgamma(q-(1/a)))/(mutgamma§*tgamma(q)))^a);
}
}
//data block
data {
int<lower=1> N; // number of observations
real y[N]; // amounts of claims
real y0;
int<lower=1> I; //number of policy (i)
int<lower=1> J; //number of develop (j)
int<lower=1> I1; //dimension of alpha
int<lower=1> J1; //dimension of beta
int<lower=1, upper=I> poli[N]; //policy id
int<lower=1, upper=J> develop[N]; //develop id
}
parameters {
real mu0; //mean effect
real alpha[I1];
real gamma1;
real gamma2;
real a;
real p;
real q;
}
transformed parameters {
real mu[N];
real alpha18;
alpha18=-sum(alpha[1:I1]);
for (i in 1:N) {
if (develop[i]==1 && poli[i]!=18) {mu[i]=exp(mu0+alpha[poli[i]]+gamma1y0);} //case A with y0
if (develop[i]==18) {mu[i]=exp(mu0+alpha[poli[i]]+gamma1y[i-1]+gamma2mu[i-1]);} //case B with beta18
if (poli[i]==18) {mu[i]=exp(mu0+alpha18+gamma1y0);} //case C with alpha18
if (develop[i]!=1 && develop[i]!=18) {mu[i]=exp(mu0+alpha[poli[i]]+gamma1y[i-1]+gamma2mu[i-1]);} //else
}
}
model {
mu0~ gamma(1,1);
alpha[1:I1] ~ normal(0,1);
gamma1 ~ gamma(1,1);
gamma2 ~ gamma(1,1);
a ~ gamma(1,1);
p ~ gamma(1,1);
q ~ gamma(1,1);
for (i in 1:N){
y[i]~beta2(a, mu[i], p, q);
}
}
generated quantities {
real loglik[N];
real dev;
for (i in 1:N){
loglik[i]=log(fabs(a))+(a*p-1)log(y[i])+lgamma(p+q)+(ap)lgamma(p+(1/a))+(ap)lgamma(q-(1/a))-(ap)log(mu[i])-(ap)lgamma§-(ap)*lgamma(q)-lgamma§-lgamma(q)-(p+q)*log(1+((y[i]*tgamma(p+(1/a))*tgamma(q-(1/a)))/(mu[i]*tgamma§*tgamma(q)))^a);
}
dev = -2*sum(loglik[]);
}
"
datt = list(N=N,y=y,y0=y0,I=I,J=J,I1=I1,J1=J1,poli=poli,develop=develop)
fitA = stan(model_code = trialmodelA,
model_name = “Trial modelA”,
data = datt, iter = iters,
warmup=burn,
chains = 1)
The data is attached.Loss Reserve Data.txt (2.6 KB)