DNAprofiles (version 0.3.1)

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

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

Value

numeric (vector) of probabilities

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*\theta+(1-\theta)*p_i)/(1+(n-1)*\theta)$

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.

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