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

Description Usage Arguments Details Value Author(s) References See Also Examples

View source: R/DirichSampSat.R

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

1
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 (jonno@u.washington.edu)

References

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

See Also

DirichSampHWE, DirichNormSat, DirichNormHWE

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
# 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 30, 2017, 2:34 a.m.