DirichSampSat: Simulate samples from a Dirichlet prior or posterior under...

In HWEBayes: Bayesian investigation of Hardy-Weinberg Equilibrium via estimation and testing.

Description

Function to simulate samples from the satuated Dirichlet model. Can be used for samples from the prior or the (conjugate) Dirichlet posterior, both in the k allele case. Samples are generated for the genotype frequencies in the order p_{11}, p_{12},..., p_{1k},p_{22} ..., p_{2k},..., p_{kk}, the allele frequencies, and the fixation indices.

Usage

DirichSampSat(nvec, bvec, nsim)

Arguments

nvec	vector of genotype frequencies in the order n_{11}, n_{12},..., n_{1k},n_{22} ..., n_{2k},..., n_{kk}.
bvec	vector of length k(k+1)/2 Dirichlet prior parameters, where k is the number of alleles.
nsim	number of samples to simulate from the prior/posterior.

Details

Uses the rdirichlet function from the MCMCpack library.

Value

pvec	matrix of size nsim \times k(k+1)/2 containing samples for the genotype frequencies, in the order p_{11}, p_{12},..., p_{1k},p_{22} ..., p_{2k},..., p_{kk}.
pmat	matrix of size nsim \times k(k+1)/2 \times k(k+1)/2 containing samples for the genotype probabilities.
pmarg	matrix of size nsim \times k containing samples for the allele frequencies, in the order p_{1},..., p_k.
fixind	matrix of size nsim \times k(k+1)/2 \times k(k+1)/2 containing samples for the fixation indices, contained in the lower diagonal, i.e. fixind[,i,j] for [i>j].

Author(s)

Jon Wakefield ([e-mail])

References

Wakefield, J. (2010). Bayesian methods for examining Hardy-Weinberg equilibrium. Biometrics; Vol 66:257-65

See Also

DirichSampHWE, DirichNormSat, DirichNormHWE

Examples

# First sample from the prior
PriorSampSat <- DirichSampSat(nvec=rep(0,10),bvec=rep(1,10),nsim=1000)
par(mfrow=c(1,2))
hist(PriorSampSat$pvec[,1],xlab="p1",main="")
hist(PriorSampSat$fixind[,2,1],xlab="f21",main="")
# Now sample from the posterior
data(DiabRecess)
PostSampSat <- DirichSampSat(nvec=DiabRecess,bvec=rep(1,10),nsim=1000)
par(mfrow=c(1,2))
hist(PostSampSat$pvec[,1],xlab="p1",main="")
hist(PostSampSat$fixind[,2,1],xlab="f21",main="")

HWEBayes documentation built on May 2, 2019, 11:07 a.m.