pr.next.alleles: Probability of seeing next alleles (Dirichlet sampling)

In DNAprofiles: DNA Profiling Evidence Analysis

Description

Probability of seeing next alleles (Dirichlet sampling)

Usage

pr.next.alleles(ij, seen, fr, theta = 0)

Arguments

ij	integer matrix with allele numbers
seen	integer matrix with alleles already seen
fr	numeric vector with allelic proportions
theta	numeric background relatedness

Details

When a population is subdivided into subpopulations, consecutively sampled alleles are not independent draws. This function implements the Dirichlet formula which states that after sampling n alleles, of which m are of type A_i, the probability that the next allele is of type A_i equals:

(m*θ+(1-θ)*p_i)/(1+(n-1)*θ)

The alleles already sampled are passed as the rows of the matrix seen, while the corresponding row of ij specifies for which alleles the probability of sampling is computed. The numer of rows of ij has to be equal to the number of rows of seen.

Value

numeric (vector) of probabilities

See Also

pr.next.allele, rmp

Examples

# compute the pr. of seeing 1,1 after 1,1,1,1
# when theta=0 this is simply p_1^2
pr.next.alleles(t(c(1,1)),seen=t(c(1,1,1,1)),fr=c(1/4,3/4),theta=0)
# but when theta>0, the pr. of seeing more 1's increases slightly
pr.next.alleles(t(c(1,1)),seen=t(c(1,1,1,1)),fr=c(1/4,3/4),theta=0.05)

# pr. distribution of (ordered!) genotypes after seeing 1,1,1
ij=matrix(c(1,1,1,2,2,2),ncol=2,byrow=TRUE)
seen=matrix(1,nrow=3,ncol=3,byrow=TRUE)
pr.next.alleles(ij,seen,fr=c(1/4,3/4),theta=0) # theta=0
pr.next.alleles(ij,seen,fr=c(1/4,3/4),theta=0.1) # theta=0.1

p0 <- pr.next.alleles(ij,seen,fr=c(1/4,3/4),theta=0)
stopifnot(all.equal(p0[1]+2*p0[2]+p0[3],1))

p1 <- pr.next.alleles(ij,seen,fr=c(1/4,3/4),theta=0.05)
stopifnot(all.equal(p1[1]+2*p1[2]+p1[3],1))

DNAprofiles documentation built on Jan. 15, 2017, 9:27 p.m.