R/psi_complex.R

In siamakz/BNPseq: Bayesian Nonparametric Differential Expression Analysis of Sequencing Count Data

psi_complex <- function(z)
{
    siz <-  dim(z)
    z <- c(z)
    zz <- z
    f <-  0*z 
    p <- which(Re(z)<0.5)
    if (length(p)>0)  z[p]<- 1-z[p]
    
    g<- 607/128 
    
    c <-  c(0.99999999999999709182,
           57.156235665862923517,
           -59.597960355475491248,
           14.136097974741747174,
           -0.49191381609762019978,
           0.33994649984811888699e-4,
           0.46523628927048575665e-4,
           -0.98374475304879564677e-4,
           0.15808870322491248884e-3,
           -0.21026444172410488319e-3,
           0.21743961811521264320e-3,
           -0.16431810653676389022e-3,
           0.84418223983852743293e-4,
           -0.26190838401581408670e-4,
           0.36899182659531622704e-5)
    
    
    n <- 0
    d <- 0
    for (k in seq(length(c), 2, -1))
    {
        dz <- 1/(z+k-2)
        dd <- c[k]*dz
        d <- d+dd
        n <- n-dd*dz
    }
    
    d <- d+c[1]
    gg <- z+g-0.5
    
    f <-  log(gg) + (n/d - g/gg) 
    
    if (length(p)>0)
        f[p] <-  f[p]-pi/tan(pi*zz[p])
    
    p <- which((round(zz) == zz) & (Re(zz) <= 0) & (Im(zz) == 0))
    if (length(p)>0)  f(p) <-  Inf
    
    dim(f) <- siz
    f
}