Re: [S] Idle curiosity -- Summary of responses

Rob Creecy (rcreecy@Census.GOV)
Tue, 17 Feb 1998 10:46:42 -0500


Dave Krantz wrote:
> ..........
> Fourth, contrary to Jens Oehlschlaegel, I don't expect that Bayesians
> will do much laughing about this problem. Formulating a sensible
> hierarchical prior leading to a good Bayesian analysis for this
> problem strikes me as an interesting challenge.
>
>

Actually, this seemed very straightforward using a Bayesian approach,
so much so that I sat down and wrote the simple function below.
The idea is to parameterize the problem in terms of the 10 probabilities
that a person prefers number j, theta_j, j=1,10. The conjugate prior
to the multinomial is Dirichlet, with "prior counts" added to each cell.
Choosing priorcount = 1 gives a uniform density, choosing priorcount = 0
gives a prior density uniform in log(theta_j)'s

e.g see Gelman, Carlin, Stern & Rubin "Bayesian Data Analysis",p 76

The posterior is then also Dirichlet with priorcount added to each cell.
It is strightforward to then generate some large number of samples from
a Dirichlet, which can be done by generating gamma's (Gelman et al p
482),
and then testing to see if theta_7 is bigger than any other theta_j.

numpref<-
function (counts, nloop = 10000, biggest = 7, priorcount = 1)
{
# A simulation to test the hypothesis that the proportion
# of people preferring the number 7 (or biggest) is greater
# than the proportion for any other number.
# This simulation generates nloop samples from a Dirichlet
# distribution of dimension 10 (or length(counts))with
parameters counts+priorcount
# and counts the number of times 7 is the largest proportion
(numbiggest)
# The function returns the proportion of samples that 7 is NOT
the biggest
numbiggest <- 0
for (i in 1:nloop) {
# add 1 to each cell count for uniform prior
x <- rgamma(length(counts), counts + priorcount)
# x<-x/sum(x) yields sample from Dirichlet - not needed
since just looking for max
if (max(x) == x[biggest])
numbiggest <- numbiggest + 1
}
1 - numbiggest/nloop
}

> counts<-c( 9, 18, 24 , 17 , 38 , 30 , 70 ,12 ,11, 8 )
> counts
[1] 9 18 24 17 38 30 70 12 11 8
> sum(counts)
[1] 237
> numpref(counts)
[1] 0.0018
# Run again, there will be some variation because of the simulation
> numpref(counts)
[1] 0.0013
# If you don't like a priorcount of 1, try 0
>numpref(counts,priorcount=0)
> numpref(counts,priorcount=0)
[1]6e-004
> numpref(counts,priorcount=0)
[1]9e-004
-----------------------------------------------------------------------
This message was distributed by s-news@wubios.wustl.edu. To unsubscribe
send e-mail to s-news-request@wubios.wustl.edu with the BODY of the
message: unsubscribe s-news