inst/extdata/nptfc.R

In xiangjiexue/npfixedcomp: Non-parametric methods for estimating mixing distribution

# Beta is the degrees of freedom. Assume here df > 2

# npt family
npt = function(v, w=1) {
  if(class(v) == "npt") { v$w = v$w * w; v }
  else structure(list(v=v, w=w), class="npt")
}

length.npt = function(x) length(x$v)

weights.npt = function(x, beta) x$w

gridpoints.npt = function(x, beta, grid=100) {
  rx = range(x$v)
  breaks = pmax(ceiling(diff(rx) / (5 * sqrt(beta / (beta - 2)))), 5)   # number of breaks
  r = whist(x$v, x$w, breaks=breaks, freq=FALSE, plot=FALSE)
  i = r$density != 0
  i = i | c(i[-1],FALSE) | c(FALSE,i[-length(i)])  # include neighbours
  m = sum(i)
  k = pmax(ceiling(grid / m), 10)           # at least 10 in each interval
  d = r$breaks[2] - r$breaks[1]
  s = r$breaks[-length(r$breaks)][i]
  c(rx[1], rep(s, rep(k,m)) + d * (1:k-0.5)/k, rx[2])
}

initial.npt = function(x, beta=NULL, mix=NULL, kmax=NULL) {
  if(is.null(beta)) stop("Error. need to specify degrees of freedom.")
  if(!is.finite(beta)) stop("degrees of freedom is Inf, consider using npnorm")
  if(is.null(mix) || is.null(mix$pt)) {
    if(is.null(kmax)) breaks = pmax(ceiling(diff(range(x$v)) / (5 * sqrt(beta / (beta - 2)))), 5)
    else {
      if(kmax == 1)
        return(list(beta=beta,
                    mix=disc(sum(x$v*x$w) / sum(rep(x$w,len=length(x$v))))))
      breaks = seq(min(x$v), max(x$v), len=pmin(20, kmax))
    }
    r = whist(x$v, x$w, breaks=breaks, freq=FALSE, plot=FALSE)
    i = r$density != 0 
    mix = disc(r$mids[i], r$density[i])
  }
  list(beta=beta, mix=mix)
}

valid.npt = function(x, beta, theta) beta > 2

suppspace.npt = function(x, beta) c(-Inf,Inf)

logd.npt = function(x, beta, pt, which=c(1,0,0)) {
  # It does not make sense in this case to have derivative wrt beta.
  # df needs to be specified or use default as a constant  
  dl = vector("list", 3)
  names(dl) = c("ld","db","dt")
  x.pt = npfixedcomp:::dt1(x$v, beta, pt, FALSE)
  if (which[1] == 1){
    dl$ld = log(x.pt)
    dim(dl$ld) = c(length(x$v), length(pt))
  }
  if (which[3] == 1){
    dl$dt = npfixedcomp:::ddtdmu1(x$v, beta, pt) / x.pt - rep(pt, each = length(x$v))
    dim(dl$dt) = c(length(x$v), length(pt))
  }
  dl
}

plot.npt = function(x, mix, beta, breaks=NULL, col="red", len=100,
                       add=FALSE, border.col=NULL, border.lwd=1,
                       fill="lightgrey", main, lwd=2, lty=1, xlab="Data",
                       ylab="Density", 
                       components=c("proportions","curves","null"),
                       lty.components=2, lwd.components=2, ...) {
  components = match.arg(components)
  m = length(mix$pt)
  z = sort(unique(round(mix$pt +
                        rep(sqrt(beta / (beta - 2))*c(-10*exp(10*10:0), seq(-4,4,len=len),
                                   10*exp(10*10:0)),
                            each=m),
                        ceiling(-log10(sqrt(beta / (beta - 2))/20)))))
  nz = length(z)
  dj = npfixedcomp:::dt1(x = z, ncp = mix$pt, df=beta, lg = FALSE) * rep(mix$pr, rep(nz,length(mix$pr)))
  d = rowSums(dj)
  if(missing(main))
    main = substitute("npt (" * df ~ "=" ~ xxx * ")",
                      list(xxx=beta))
  if(add || missing(x) || length(x) == 0) {
    if(add) lines(z, d, col=col, lwd=lwd, lty=lty)
    else {
      plot(0, 0, type="n", xlim=range(mix$pt) + sqrt(beta / (beta - 2)) * c(-3,3), ylim=range(d),
           xlab=xlab, ylab=ylab, frame.plot=FALSE,
           main=main, ...)
      lines(range(z), rep(0,2), col="darkgrey")
      lines(z, d, col=col, lwd=lwd, lty=lty) 
    }
  }
  else {
    if(is.null(breaks)) breaks = 10 + round(sqrt(length(x$v)))
    whist(x$v, x$w, breaks=breaks, freq=FALSE, 
          xlab=xlab, ylab=ylab, main=main, col=fill, border=border.col, 
          lwd=border.lwd, ylim=range(d), ...)
    lines(z, d, col=col, lwd=lwd, lty=lty)
  }
  if(components != "null") {
    points(mix$pt, rep(0,length(mix$pt)), col=col)
    switch(components,
           proportions = {
             segments(mix$pt, rep(0,m), y1=mix$pr*max(d), col=col, lwd=3)
           },
           curves = {
             if(ncol(dj) > 1) {
               j2 = 1:floor(ncol(dj)/2) * 2
               for(j in 1:ncol(dj))     # slow if too many components
                 lines(z, dj[,j], col=col, lty=lty.components,
                       lwd=lwd.components)
             }
           }
           )
  }
}




# nptfc family
nptfc = function(v, pi0 = 0, mu0 = 0, df, w = 1){
  if (length(pi0) != length(mu0))
    stop("Lengths of pi0 and mu0 do not match")
  structure(list(v=v, w = w, df = df, pi0 = pi0, mu0 = mu0), 
            class="nptfc")
}

length.nptfc = function(x) length(x$v)

weights.nptfc = function(x, beta) x$w

gridpoints.nptfc = function(x, beta, grid=100) {
  rx = range(x$v)
  breaks = pmax(ceiling(diff(rx) / (5 * sqrt(beta / (beta - 2)))), 5)   # number of breaks
  r = whist(x$v, x$w, breaks=breaks, freq=FALSE, plot=FALSE)
  i = r$density != 0
  i = i | c(i[-1],FALSE) | c(FALSE,i[-length(i)])  # include neighbours
  m = sum(i)
  k = pmax(ceiling(grid / m), 10)           # at least 10 in each interval
  d = r$breaks[2] - r$breaks[1]
  s = r$breaks[-length(r$breaks)][i]
  c(rx[1], rep(s, rep(k,m)) + d * (1:k-0.5)/k, rx[2])
}

initial.nptfc = function(x, beta, mix=NULL, kmax=NULL) {
  if (is.null(beta)) beta = x$df
  if (beta != x$df) stop("Inconsistent structural parameter beta.")
  if(is.null(mix) || is.null(mix$pt)) {
    if(is.null(kmax)) breaks = pmax(ceiling(diff(range(x$v)) / (5 * sqrt(beta / (beta - 2)))), 5)
    else {
      if(kmax == 1)
        return(list(beta=beta, mix=disc(sum(x$v*x$w) / sum(rep(x$w,len=length(x$v))))))
      breaks = seq(min(x$v), max(x$v), len=pmin(20, kmax))
    }
    r = whist(x$v, x$w, breaks=breaks, freq=FALSE, plot=FALSE)
    # r$density = pmax(0, r$density - npfixedcomp:::precomputenptfc(r$mids, mu0 = x$mu0, pi0 = x$pi0, beta = beta))
    i = r$density != 0 
    mix = disc(r$mids[i], r$density[i])
  }
  list(beta=beta, mix=mix)
}

valid.nptfc = function(x, beta, theta) beta > 2

suppspace.nptfc = function(x, beta) c(-Inf,Inf)

logd.nptfc = function(x, beta, pt, which=c(1,0,0)) {
  dl = vector("list", 3)
  names(dl) = c("ld","db","dt")
  precompute = .rowSums(npfixedcomp:::dt1(x$v, beta, x$mu0, FALSE) * rep(x$pi0, each = length(x$v)), m = length(x$v), n = length(x$mu0))
  x.pt = npfixedcomp:::dt1(x$v, beta, pt, FALSE)
  if (which[1] == 1){
    dl$ld = log(x.pt * (1 - sum(x$pi0)) + precompute)
    dim(dl$ld) = c(length(x$v), length(pt))
  }
  if (which[3] == 1){
    dl$dt = (npfixedcomp:::ddtdmu1(x$v, beta, pt) - rep(pt, each = length(x$v)) * x.pt) * (1 - sum(x$pi0)) / (x.pt * (1 - sum(x$pi0)) + precompute)
    dim(dl$dt) = c(length(x$v), length(pt))
  }
  dl
}

makenptfc = function(data, mu0 = 0, pi0 = 0, df, mix = NULL, order = FALSE, ...){
  if (length(mu0) != length(pi0) | sum(pi0) > 1 | any(pi0 < 0)){
    stop("Invalid inputs for mu0 or pi0!")
  }
  if (!order){
    temp = list(v = data, w = 1)
  }else{
    temp = bin(data, order = order)
  }
  switch(as.character(length(mu0)),
         "1" = switch(as.character(pi0),
                      "1" = {
                        result1 = list(family = "nptfc", num.iterations = 0, max.gradient = 0, convergence = 0, 
                                       ll = sum(dt1(x = temp$v, df = df, ncp = mu0, lg = TRUE) * temp$w), mix = list(pt = mu0, pr = 1), beta = df, 
                                       pi0 = 1, mu0 = mu0)
                        class(result1) = "nspmix"
                      },
                      "0" = {
                        result1 = nspmix::cnm(npt(v = temp$v, w = temp$w), init = list(beta = df, mix = mix), ...)
                        result1$pi0 = 0
                        result1$mu0 = mu0
                        result1$family = "nptfc"
                      },
                      {
                        result1 = nspmix::cnm(nptfc(v = temp$v, mu0 = mu0, pi0 = pi0, df = df, w = temp$w), 
                                              init = list(beta = df, mix = mix), ...)
                        result1$family = "nptfc"
                        result1$pi0 = pi0
                        result1$mu0 = mu0
                      }),
         {
           result1 = nspmix::cnm(nptfc(v = temp$v, mu0 = mu0, pi0 = pi0, df = df, w = temp$w), 
                                 init = list(beta = df, mix = mix), ...)
           result1$family = "nptfc"
           result1$pi0 = pi0
           result1$mu0 = mu0
         })
  
  nptfc2npt(result1)
}

#' convert nptfc object to npt object
#' 
#' convert nptfc object to npt object
#' 
#' @title convert nptfc object to npt object
#' @param result an object of family nptfc
#' @return an object of family npt
#' @export
nptfc2npt = function(result){
  if (result$family != "nptfc"){
    stop("Input is not the family of nptfc")
  }
  result1 = result
  result1$pi0 = NULL
  result1$mu0 = NULL
  result1$family = "npt"
  result1$mix = disc(pt = c(result$mix$pt, result$mu0), pr = c(result$mix$pr * (1 - sum(result$pi0)), result$pi0))
  result1
}

#' plot nptfc object
#' 
#' plot nptfc object
#' 
#' @title plot nptfc object
#' @param result an object of family nptfc
#' @param ... additional arguments passed to plot.npt or its underlying functions
#' @export
plotnptfc = function(result, ...){
  plot.npt(mix = result$mix, beta = result$beta, ...)
}


# the function value used in findnpnormfc
findnptfc.fn = function(d, data, pval, df, order){
  -2 * makenptfc(data, pi0 = d, df = df, order = order)$ll + pval
}

# the derivative function used in findnpnormfc
findnptfc.gr = function(d, data, pval, df, order){
  result1 = makenptfc(data, pi0 = d, df = df, order = order)
  -2 * (sum(dt1(x = data, df = df, ncp = 0, lg = FALSE) / 
              dnpt(data, pi0 = result1$mix$pr, mu0 = result1$mix$pt, df)) - length(data)) /
    (1 - d)
}

findnpt1fc = function(data, val = -log(length(data)), df, 
                      order = -4, ...){
  result1 = makenptfc(data, pi0 = 0, df = df, order = order)
  result2 = makenptfc(data, pi0 = 1, df = df, order = order)
  if (-2 * (result2$ll - result1$ll) + val < 0){
    ans1 = result2
  }else{
    # give a relative good estimate
    pix = c(dnpt(0, pi0 = result1$mix$pr,
                 mu0 = result1$mix$pt, df) / dt1(x = 0, df = df, ncp = 0, lg = FALSE))
    ans2 = NR(pix, fn = findnptfc.fn, gr = findnptfc.gr,
              lower = 0, upper = 1, data = data, pval = val + 2 * result1$ll, 
              df = df, order = order)
    ans1 = makenptfc(data, pi0 = ans2$par, df = df, order = order)
  }
  ans1
}