system reliability modeling using Negative log-gamma priors

[maxwell]

Dear STAn users, 

The following STAN code models a series system reliability where pass/fail data are collected at both components and system level. We specify a Negative log gamma (nlg) prior for each of the component prior. Clearly, that also implies a prior for the system reliability since R_syst = prod_i=1 to n R_compi. Note that if X ~ gamma(alpha, beta), then R = exp(-X) ~ NLG(alpha, beta). I hope that's enough as a background.

Bcca_mod <- '

data {

real<lower=0> t0[2]; // time which the priors are specified

real<lower=0> al[2]; // vectors of al for prior

real<lower=0> bet[2]; // vectors of bet values 

int<lower=1> N; // length of data time points

int<lower=1> nc; // number of component

int<lower=0> nt; // number of test at each time points

real<lower=0> pi_m; // pi mixing proportion

int yc[N,nc]; // component data -- component data by columns

int ys[N]; // system data

real<lower=0> t[N]; // time of data collection

}

parameters {

real<lower=0> x0[nc];   // R(0) = exp(-xo) ## defines the negative log-gamma prior (NLG) prior

real<lower=0> xdel[nc]; // Delta = exp(-del)

}

transformed parameters {

real y0[nc];

real yd[nc];

real ytc[N,nc]; 

vector[nc] Rc[N];

real Rs[N];

real b[nc];

for(i in 1:nc)

{

y0[i] <- -x0[i] - log(-expm1(-x0[i]));  

yd[i] <- -xdel[i] - log(-expm1(-xdel[i])); 

b[i] <- -log1p(exp(-(yd[i]-y0[i]))+exp(-yd[i]));

for(j in 1:N)

{

ytc[j,i]<- (t[j]/(t0[2]-t0[1]))*b[i]+ y0[i];

Rc[j,i] <- inv_logit(ytc[j,i]);

}

}

// system reliability

for(i in 1:N)

Rs[i] <- prod(Rc[i]);

}

// done defining the transformed varaibles

model {

x0 ~ gamma(al[1],bet[1]);

xdel ~ gamma(al[2],bet[2]);

for(j in 1:N)

{

for(i in 1:nc)

yc[j,i] ~ binomial(nt,Rc[j,i]); // component part

ys[j] ~ binomial(nt,Rs[j]); // system part   

}

}

'

nc=3

t0=c(0,1)

nt=20

t=seq(0,1.5,by=.5)

n_ac=1

pam=c(2/nc,.001,5.421,.05)

al=c(2,2.269)/nc    ## shape parameters -- component priors

bet=c(.01,.0431)    ## scale parameters -- component priors

yc=matrix(rep(20,4*nc),ncol=nc)  ### component test data

yc[,1]=c(19,19,18,17)

pi_m=.4

ys=c(19,18,17,16) ## system test data

20 tests are collected at each of 4 time points for each components (data in yc) and system (data in ys)

Note that only component 1 has failures and the remaining components have no failures. 

I get the follwing error message over and over again. Anyone has an idea about what is going on? My guess is I may running into loss of precision with some of my parameters. Any suggestions about ways to ameliorate the code? I am actually able to simulate from the posterior using MH.

SAMPLING FOR MODEL 'Bcca_mod' NOW (CHAIN 1).

Iteration:     1 / 50000 [  0%]  (Adapting)

Informational Message: The parameter state is about to be Metropolis rejected due to the following underlying, non-fatal (really) issue (and please ignore that what comes next might say 'error'): Error in function stan::prob::binomial_log(N4stan5agrad3varE): Probability parameter is 1.#QNAN:0, but must be finite!

If the problem persists across multiple draws, you might have a problem with an initial state or a gradient somewhere.

 If the problem does not persist, the resulting samples will still be drawn from the posterior.

Informational Message: The parameter state is about to be Metropolis rejected due to the following underlying, non-fatal (really) issue (and please ignore that what comes next might say 'error'): Error in function stan::prob::binomial_log(N4stan5agrad3varE): Probability parameter is 1.#QNAN:0, but must be finite!

If the problem persists across multiple draws, you might have a problem with an initial state or a gradient somewhere.

 If the problem does not persist, the resulting samples will still be drawn from the posterior.

Thank you

[Bob Carpenter]

The message indicates that NaN (not-a-number)
values into one of the binomials:

> yc[j,i] ~ binomial(nt,Rc[j,i]); // component part

> ys[j] ~ binomial(nt,Rs[j]); // system part

You have a lot of arithmetic in the model, some of which is
numerically unstable near boundaries (in the sense of
having possible underflows or overflows -- for instance, with
values less than -700 or so, exp underflows to 0).

Given their definitions:

> for(i in 1:nc)
> ...

> for(j in 1:N) {
> ytc[j,i]<- (t[j]/(t0[2]-t0[1]))*b[i]+ y0[i];
> Rc[j,i] <- inv_logit(ytc[j,i]);
> }

and

> for(i in 1:N)
> Rs[i] <- prod(Rc[i]);

what's probably happening is that one or more of the ytc values are
infinite or NaN. They can be NaN if the t0[2] - t0[1] values
are NaN, or if any of the other numbers are NaN.

So I'd try diagnosing by printing ytc and Rc values in
the code to see where they're going wrong, then maybe doing
the same for the inputs you think may not be right.

I'd also recomment using binomial_logit, where you don't
need to do the inverse logit transform --- it's more numerically
stable.

- Bob

> --
> You received this message because you are subscribed to the Google Groups "stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> For more options, visit https://groups.google.com/groups/opt_out.

[Roger Zoh]

Hi Bob,

Thanks for your suggestions. How do you print the value of the ytc and Rc in the code?

 

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

For more options, visit https://groups.google.com/groups/opt_out.

--
You received this message because you are subscribed to the Google Groups "stan users mailing list" group.

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

For more options, visit https://groups.google.com/groups/opt_out.

--
Roger Zoh 

Postdoc.

Department of Statistics

and Nutrition

Blocker 525B and Kleberg 318

Texas A&M University

[Bob Carpenter]

You can just use the print() statement in Stan.
It accepts any sequence of expressions and strings.
So you can do

print("ytc=",ytc);

to print out all of ytc, or for a finer-grained print
with indexes, something like:

for (i in N)
for (j in nc)
print("ytc",i,j,ytc[i,j]);

where I and J are replaced with whatever the sizes are.

- Bob

> b[i] <- -log1p(exp(-(yd[i]-y0[i]))+__exp(-yd[i]));

> for(j in 1:N)
> {
> ytc[j,i]<- (t[j]/(t0[2]-t0[1]))*b[i]+ y0[i];
> Rc[j,i] <- inv_logit(ytc[j,i]);
> }
> }
> // system reliability
> for(i in 1:N)
> Rs[i] <- prod(Rc[i]);
> }
> // done defining the transformed varaibles
>
> model {
>
> x0 ~ gamma(al[1],bet[1]);
> xdel ~ gamma(al[2],bet[2]);
> for(j in 1:N)
> {
> for(i in 1:nc)
> yc[j,i] ~ binomial(nt,Rc[j,i]); // component part
>
> ys[j] ~ binomial(nt,Rs[j]); // system part
> }
> }
> '
> nc=3
> t0=c(0,1)
> nt=20
> t=seq(0,1.5,by=.5)
> n_ac=1
> pam=c(2/nc,.001,5.421,.05)
> al=c(2,2.269)/nc ## shape parameters -- component priors
> bet=c(.01,.0431) ## scale parameters -- component priors

> yc=matrix(rep(20,4*nc),ncol=__nc) ### component test data

> yc[,1]=c(19,19,18,17)
> pi_m=.4
> ys=c(19,18,17,16) ## system test data
>
> 20 tests are collected at each of 4 time points for each components (data in yc) and system (data in ys)
> Note that only component 1 has failures and the remaining components have no failures.
>
> I get the follwing error message over and over again. Anyone has an idea about what is going on? My guess is I may
> running into loss of precision with some of my parameters. Any suggestions about ways to ameliorate the code? I am
> actually able to simulate from the posterior using MH.
>
> SAMPLING FOR MODEL 'Bcca_mod' NOW (CHAIN 1).
> Iteration: 1 / 50000 [ 0%] (Adapting)
> Informational Message: The parameter state is about to be Metropolis rejected due to the following underlying,
> non-fatal
> (really) issue (and please ignore that what comes next might say 'error'): Error in function

> stan::prob::binomial_log(__N4stan5agrad3varE): Probability parameter is 1.#QNAN:0, but must be finite!

> If the problem persists across multiple draws, you might have a problem with an initial state or a gradient
> somewhere.
> If the problem does not persist, the resulting samples will still be drawn from the posterior.
>
> Informational Message: The parameter state is about to be Metropolis rejected due to the following underlying,
> non-fatal
> (really) issue (and please ignore that what comes next might say 'error'): Error in function

> stan::prob::binomial_log(__N4stan5agrad3varE): Probability parameter is 1.#QNAN:0, but must be finite!

> If the problem persists across multiple draws, you might have a problem with an initial state or a gradient
> somewhere.
> If the problem does not persist, the resulting samples will still be drawn from the posterior.
>
> Thank you
>
> --
> You received this message because you are subscribed to the Google Groups "stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to

> stan-users+unsubscribe@__googlegroups.com <mailto:stan-users%2Bunsu...@googlegroups.com>.
> For more options, visit https://groups.google.com/__groups/opt_out <https://groups.google.com/groups/opt_out>.

>
>
> --
> You received this message because you are subscribed to the Google Groups "stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to

> stan-users+unsubscribe@__googlegroups.com <mailto:stan-users%2Bunsu...@googlegroups.com>.
> For more options, visit https://groups.google.com/__groups/opt_out <https://groups.google.com/groups/opt_out>.

>
>
>
>
> --
> Roger Zoh
> Postdoc.
> Department of Statistics
> and Nutrition
> Blocker 525B and Kleberg 318
> Texas A&M University
>

> --
> You received this message because you are subscribed to the Google Groups "stan users mailing list" group.

> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.