R Graphical Manual

Browse All

Last data update: 2014.03.03

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

DirichSampSat

R Documentation

Simulate samples from a Dirichlet prior or posterior under the saturated model

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

References

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

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="")

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/DirichSampSat.Rd_%03d_medium.png", width=480, height=480)
> ### Name: DirichSampSat
> ### Title: Simulate samples from a Dirichlet prior or posterior under the
> ###   saturated model
> ### Aliases: DirichSampSat
> ### Keywords: distribution
> 
> ### ** 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="")
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>