# R/log_stirling2.R In rhierbaps: Clustering Genetic Sequence Data Using the HierBAPS Algorithm

#### Documented in log_stirling2

```#' log_stirling2
#'
#' @param n number of objects
#' @param k number of partitions
#'
#' @return log of the Stirling number of the second kind
#'
#'
log_stirling2 <- function(n, k){
if(!is.numeric(n)) stop("n is not numeric!")
if(!is.numeric(k)) stop("k is not numeric!")
if(k>n) stop("k must be less than n!")

v <- n/k
G <- lambertW(-v*exp(-v))

lS2 <- log(sqrt((v-1)/(v*(1-G)))) +
(n-k)*(log(v-1)-log(v-G)) +
n*log(k)-k*log(n) +
k*(1-G) +
lchoose(n, k)

return(lS2)
}

# This function was written by Ben Bolker and taken from https://stat.ethz.ch/pipermail/r-help/2003-November/042793.html
lambertW = function(z,b=0,maxiter=10,eps=.Machine\$double.eps,
min.imag=1e-9) {
if (any(round(Re(b)) != b))
stop("branch number for W must be an integer")
if (!is.complex(z) && any(z<0)) z=as.complex(z)
##
## p = (1 - 2*abs(b)).*sqrt(2*e*z + 2);
## w = (11/72)*p;
## w = (w - 1/3).*p;
## w = (w + 1).*p - 1
##
## first-order version suffices:
##
w = (1 - 2*abs(b))*sqrt(2*exp(1)*z + 2) - 1
## asymptotic expansion at 0 and Inf
##
v = log(z + as.numeric(z==0 & b==0)) + 2*pi*b*1i;
v = v - log(v + as.numeric(v==0))
## choose strategy for initial guess
##
c = abs(z + exp(-1));
c = (c > 1.45 - 1.1*abs(b));
c = c | (b*Im(z) > 0) | (!Im(z) & (b == 1))
w = (1 - c)*w + c*v
## Halley iteration
##
for (n in 1:maxiter) {
p = exp(w)
t = w*p - z
f = (w != -1)
t = f*t/(p*(w + f) - 0.5*(w + 2.0)*t/(w + f))
w = w - t
if (abs(Re(t)) < (2.48*eps)*(1.0 + abs(Re(w)))
&& abs(Im(t)) < (2.48*eps)*(1.0 + abs(Im(w))))
break
}
if (n==maxiter) warning(paste("iteration limit (",maxiter,
") reached, result of W may be inaccurate",sep=""))
if (all(Im(w)<min.imag)) w = as.numeric(w)
return(w)
}
```

## Try the rhierbaps package in your browser

Any scripts or data that you put into this service are public.

rhierbaps documentation built on Jan. 11, 2020, 9:33 a.m.