R/polysFast.R

In wCorr: Weighted Correlations

#' @importFrom minqa bobyqa
#' @importFrom stats dnorm
#' @importFrom stats optimize
#' @importFrom stats cor
#' @importFrom stats weighted.mean
polysFast <- function(x, M, w, ML=FALSE) {
  M <- as.numeric(as.factor(M))
  uM <- sort(unique(M))
  mapTheta <- function(v) {
    vv <- cumsum(c(v[1],exp(v[-1])))
    c(NA,-Inf,vv,Inf)
  }
  theta0 <- sapply(uM[-length(uM)],function(z) qnorm(weighted.mean(M<=z, w)) )
  if(ML) {
    temp = w/sum(w)
    temp2 = fixxFast(x, temp)
    temp3 = sum(temp*dnorm(temp2,log=TRUE))
    bob <- suppressWarnings(bobyqa(par=c(atanh(cor(x,M)),imapThetaFast2(theta0)),
                  fn=optFFast, x=temp2,  w=temp, M=M, temp3 = temp3))
    return(tanh(bob$par[1]))
  } else {
    temp = w/sum(w)
    temp2 = fixxFast(x, temp)
    temp3 = sum(temp*dnorm(temp2, log=TRUE))
    values = mainF(x, M, w, theta0)
    opt <- suppressWarnings(optimize(optFcFast, interval=unlist(values[1]),
                    x=temp2, w=temp, theta0=(imapThetaFast2(theta0)),
                    M=M, temp3= temp3))
    return( tanh(opt$minimum) )
  }
}