R: Simulate samples from a Dirichlet prior or posterior under...
DirichSampHWE
R Documentation
Simulate samples from a Dirichlet prior or posterior under HWE
Description
Function to simulate samples from the HWE 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 allele frequencies in the order p_{1},p_{2},...,p_{k}.
Usage
DirichSampHWE(nvec, bvec0, nsim)
Arguments
nvec
vector of genotype frequencies in the order n_{11}, n_{12},..., n_{1k},n_{22} ..., n_{2k},..., n_{kk}.
bvec0
vector of length k 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 containing samples for the genotype frequencies, in the order p_{1}, p_{12},..., p_{k}.
Author(s)
Jon Wakefield (jonno@u.washington).
References
Wakefield, J. (2010). Bayesian methods for examining Hardy-Weinberg
equilibrium. Biometrics; Vol 66:257-65
See Also
DirichSampSat, DirichNormSat, DirichNormHWE
Examples
# First sample from the prior
PriorSampHWE <- DirichSampHWE(nvec=rep(0,10),bvec0=rep(1,4),nsim=1000)
par(mfrow=c(1,1))
hist(PriorSampHWE$pvec[,1],xlab="p1",main="")
# Now sample from the posterior
data(DiabRecess)
PostSampHWE <- DirichSampHWE(nvec=DiabRecess,bvec0=rep(1,4),nsim=1000)
par(mfrow=c(1,1))
hist(PostSampHWE$pvec[,1],xlab="p1",main="")
Results
R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(HWEBayes)
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/HWEBayes/DirichSampHWE.Rd_%03d_medium.png", width=480, height=480)
> ### Name: DirichSampHWE
> ### Title: Simulate samples from a Dirichlet prior or posterior under HWE
> ### Aliases: DirichSampHWE
> ### Keywords: distribution
>
> ### ** Examples
>
> # First sample from the prior
> PriorSampHWE <- DirichSampHWE(nvec=rep(0,10),bvec0=rep(1,4),nsim=1000)
> par(mfrow=c(1,1))
> hist(PriorSampHWE$pvec[,1],xlab="p1",main="")
> # Now sample from the posterior
> data(DiabRecess)
> PostSampHWE <- DirichSampHWE(nvec=DiabRecess,bvec0=rep(1,4),nsim=1000)
> par(mfrow=c(1,1))
> hist(PostSampHWE$pvec[,1],xlab="p1",main="")
>
>
>
>
>
> dev.off()
null device
1
>