estep: E-Step for Abundance Estimation
In sonicLength: Estimating Abundance of Clones from DNA Fragmentation Data

Description Usage Arguments Details Value Author(s) Examples

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

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

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

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

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

Charles C. Berry ccberry@users.r-forge.r-project.org

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