Ey.given.x {sonicLength}

E-Step for Abundance Estimation

Description

Take the E-step of the EM algorithm for estimating the number of sonicants in a sample

Usage

Ey.given.x(x, theta, phi)
pr.y.given.x(x, theta, phi, kmax=20 )

Arguments

x

a zero-one indicator matrix whose rows correspond to unique lengths with rownames indicating those lengths

theta

a vector of the abundance estimates. length(theta)==ncol(x)

phi

a vector of the probabilities of sonicant lengths. length(phi)==nrow(x)

kmax

highest count to bother with (all higher values are globbed together in the result)

Details

Supposing Poisson sampling of sonicants, Ey.given.x(x,theta,phi)[i,j] gives the expected value of the number of sonicants given that x[i,j] distinct sonicant lengths were observed. pr.y.given.x(x,theta,phi.kmax)[k,j] gives the probability of k sonicants (censored at kmax+1) of insertion site j

Value

Ey.given.x() is a double matrix of the same dimension as x. pr.y.given.x(...,kmax) is a double matrix of dimension c( kmax+1, ncol(x) )

Author(s)

Charles C. Berry [e-mail]

Examples

mat <- diag(10)
mat[ ,10 ] <- 1.0
phi1 <- prop.table( rep(1,10))
theta1 <- 1:10
Ey.given.x( mat, theta1, phi1 )
pr.mat <- pr.y.given.x( mat, theta1, phi1 )
## Estimate Seen plus Unseen sites
popsize.chao <- function(tab) sum(tab)+tab[1]*(tab[1]-1)/(2*(tab[2]+1))
## evaluate at expected counts
popsize.chao( rowSums(pr.mat) ) 
## average randomly sampled counts
cnt.samp <- function(x) sample( seq_along(x) , 1 ,pr=x )
mean(replicate(100,popsize.chao( table(apply(pr.mat,2, cnt.samp) ))))