5
$\begingroup$

I have a Bayesian model here in R/WinBUGS. This data "pumps" has two columns, time and fail. Each observation is a pump and fail[1] tells how many failures there were in time[1]. I want to know how I can write this same model in SAS using Proc MCMC. Thank you!

pump <- function() {
    for(i in 1:10){
        fail[i] ~ dpois(lam[i])
        theta[i] ~ dgamma(a,b)
        lam[i] <- theta[i]*time[i]
    }
    a ~ dgamma(1.5,.25) 
    b ~ dgamma(1.5,.25)
}

filename <- file.path(tempdir(),'pump.bug') 
write.model(pump,filename)
time <- pumps[,1]
fail <- pumps[,2]

data <- c('time','fail')
parameters <- c('fail','lam','theta','a','b')
pump.sim <- bugs(data, inits=NULL, parameters, model.file='pump.bug',
    n.iter=20000, n.burnin=1000, n.chains=1, n.thin=1, debug=T) 
$\endgroup$
3
  • $\begingroup$ Could you do a little bit of explaining what the code does in order to make it clearer to those of us who might be good at SAS but not so much familiar with what the WinBUGS code is doing? On first pass, I'd suggest passing Proc MCMC and taking a look at Proc Lifereg with the MCMC option instead. $\endgroup$
    – Fomite
    Commented Mar 28, 2012 at 4:05
  • 1
    $\begingroup$ Epi, it's all specified in the pump function: $a$ and $b$ are hyperparameters for $\theta$ (themselves with Gamma priors) and $\theta$ is a parameter for a Poisson model for the probability of failure (as a function of time). $\endgroup$
    – whuber
    Commented Mar 28, 2012 at 15:33
  • $\begingroup$ THis code does a Bayesian analysis of the survival time for a given pump to break. The data has 10 obs, so the first for loop is setting the model and priors for each obs. Then a and b are given hyper-priors. Time and Fail are the two variable in the data set. Hope this helps clarify. $\endgroup$
    – Lynn
    Commented Apr 4, 2012 at 18:54

1 Answer 1

1
$\begingroup$

See if something like this works:

proc mcmc
 data = /* your data here */
 nmc = 20000, nbi = 1000; /* iterations and burn-in */

 /* Specify parameters -- separate statements means they're proposed separately */
 parms a;
 parms b;
 parms theta;

 /* Specify priors */
 hyper a ~ gamma(1.5, scale = .25);
 hyper b ~ gamma(1.5, scale = .25); 
 prior theta ~ gamma(a, scale = b);

 /* Sampling model */
 lam = theta*time;
 model fail ~ poisson(lam);
run;
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.