Ey.given.x {sonicLength} | R Documentation |
Take the E-step of the EM algorithm for estimating the number of sonicants in a sample
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. |
phi |
a vector of the probabilities of sonicant
lengths. |
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) ))))