Re: Help with BUGS model
- From: vontressms@xxxxxx
- Date: Wed, 22 Aug 2007 06:16:43 -0700
On Aug 22, 3:27 am, "sci.stat.math" <fra...@xxxxxxxx> wrote:
This is a quote from an article by W. DuMouchel et al. 1999 "BAYESIAN
DATA MINING IN LARGE FREQUENCY TABLES":
[Quote]:
... to achieve the interpretability of the relative risk measures but
also to adjust properly for sampling variation...Assume that each
observed count Nij is a draw from a Poisson distribution with unknown
mean mu[i], and that interest centers on the ratios lambda[i]= mu[i]/
E[i]. But rather than treat the M values of lambda[i] as unrelated
constants, assume that each lambda is drawn from a common prior
distribution. This distribution is assumed be a mixture of two gamma
distributions.
...the marginal distribution of each N is a mixture of negative
binomial distributions; and second, the
posterior distribution of each lambda is a mixture of two gamma
distributions with modified parameters:
Prob(N = n) = P f(n; a1, b1, E) + (1 - P) f(n; a2, b2, E)
with f(n; a, b, E) = (1 + b/E)^-n (1 + E/b)^-a Gamma(a + n) /
Gamma(a)n!
[End Quote]
I have is values of N and E. Now I have tried to perform it in winBUGS
using the model below that assumes lambda is drawn from one gamma
distribution (but not a mixture):
model
{
for (i in 1:K)
{
N[i]~ dpois( mu[i] )
mu[i]<- lambda[i]*E[i]
lambda[i] ~ dgamma( a, b )}
a ~dunif(0.000001, 100)
b ~ dunif(0.000001, 100 )
}
It works OK I can have an estimate of lambdas and the parameters a and
b but how can I modify my BUGS model to have lambdas drawn from the
mixture of two gamma distributions like the extract above states?
Thanks
This seems to work as an extension of your code:
model
{
for (i in 1:K)
{
Signal[i] <- mu[i]/E[i]
N[i]~ dpois( mu[i] )
mu[i]<- lambda[i]*E[i]
lambda[i] <-P*lambda1[i]+(1-P)*lambda2[i] #result of equation 4.
lambda1[i] ~ dgamma( a1, b1 )
lambda2[i] ~ dgamma( a2, b2 )
}
a1 ~dunif(0.000001, 100)
b1 ~ dunif(0.000001, 100 )
a2 ~dunif(0.001, 100)
b2 ~ dunif(0.001, 100 )
P~dbeta(10,100)}
}
#data
list(K=5, N=c(1,2,3,4,5), E=c(1,1,1,1,1) )
Note that to get Equation 7 of the DuMouchel paper (1999 American
Statistician), you have to put N in the data list. That is a trick in
WinBugs to get a posterior distribution since N is defined as a
stochastic node and as data. This makes the distribution of lambda
conditional on N. The data list above is just made up.
I added a signal node which is the MCMC distribution estimate of EBGM.
You can look at the posterior distribution of EBGM using the MCMC
sample. The mean of that distribution is the EBGM and the 5'th and
95'th percentiles are the credible region for EBGM.
DuMouchel advocates using MLE's for the estimates of a1, b1, a2, b2
and P, rather than as random hyperparameters, so you have departed
from the paper a little by letting them be random. Still, your
departure may be a neat extension. You could also add the MLE's of a1,
b1, a2, b2 and P to the data list, and leave them as stochastic nodes
to judge how uncertianty around the MLE's effect the posterior
distribution of the signal.
As an aside, I think that you will find your approach to be difficult
to implement in a production environment since you have to build the
data list with actual data for each data set. I guess you could use R
and BRUGS to manage it, but that would be a maintainance nightmare.
An advantage of your approach though is that you can extend the
definitions of lambda to loglinear models, which is a limitation in
the 1999 paper.
Mark
.
- Follow-Ups:
- Re: Help with BUGS model
- From: sci.stat.math
- Re: Help with BUGS model
- From: sci.stat.math
- Re: Help with BUGS model
- References:
- Help with BUGS model
- From: sci.stat.math
- Help with BUGS model
- Prev by Date: Re: K Nearest Neighbors in time series forecasting
- Next by Date: Re: Afonso's "fiducial tail probability"
- Previous by thread: Help with BUGS model
- Next by thread: Re: Help with BUGS model
- Index(es):
Relevant Pages
|
