inst/extdata/npnormfc.R

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

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

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

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

gridpoints.npnormfc = function(x, beta, grid=100) {
  rx = range(x$v)
  breaks = pmax(ceiling(diff(rx) / (5 * beta)), 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.npnormfc = function(x, beta=NULL, mix=NULL, kmax=NULL) {
  if (is.null(beta)) beta = x$stdev
  if (beta != x$stdev) 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 * beta)), 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:::dnpnorm1(r$mids, mu0 = x$mu0, pi0 = x$pi0, stdev = beta))
    i = r$density != 0 
    mix = disc(r$mids[i], r$density[i])
  }
  list(beta=beta, mix=mix)
}

valid.npnormfc = function(x, beta, theta) TRUE

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

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

makenpnormfc = function(data, mu0 = 0, pi0 = 0, stdev = 1, 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 = "npnormfc", num.iterations = 0, max.gradient = 0, convergence = 0, 
                                       ll = sum(dnorm(temp$v, mean = mu0, sd = stdev, log = TRUE) * temp$w), mix = list(pt = 0, pr = 1), beta = stdev, 
                                       pi0 = 1, mu0 = mu0)
                        class(result1) = "nspmix"
                      },
                      "0" = {
                        result1 = nspmix::cnm(npnorm(v = temp$v, w = temp$w), init = list(beta = stdev, mix = mix), ...)
                        result1$pi0 = 0
                        result1$mu0 = mu0
                        result1$family = "npnormfc"
                      },
                      {
                        result1 = nspmix::cnm(npnormfc(v = temp$v, mu0 = mu0, pi0 = pi0, stdev = stdev, w = temp$w), 
                                              init = list(beta = stdev, mix = mix), ...)
                        result1$family = "npnormfc"
                        result1$pi0 = pi0
                        result1$mu0 = mu0
                      }),
         {
           result1 = nspmix::cnm(npnormfc(temp$v, mu0 = mu0, pi0 = pi0, stdev = stdev, w = temp$w), 
                                 init = list(beta = stdev, mix = mix), ...)
           result1$family = "npnormfc"
           result1$pi0 = pi0
           result1$mu0 = mu0
         })
  
  npnormfc2npnorm(result1)
}

npnormfc2npnorm = function(result){
  if (result$family != "npnormfc"){
    stop("Input is not the family of npnormfc")
  }
  result1 = result
  result1$pi0 = NULL
  result1$mu0 = NULL
  result1$family = "npnorm"
  result1$mix = disc(pt = c(result$mix$pt, result$mu0), pr = c(result$mix$pr * (1 - sum(result$pi0)), result$pi0))
  result1
}

# the function value used in findnpnormfc
findnpnormfc.fn = function(ddd, data, pval, stdev, order){
  -2 * makenpnormfc(data, pi0 = ddd, stdev = stdev, order = order)$ll + pval
}

# the derivative function used in findnpnormfc
findnpnormfc.gr = function(ddd, data, pval, stdev, order){
  result1 = makenpnormfc(data, pi0 = ddd, stdev = stdev, order = order)
  -2 * (sum(dnorm(data) / dnpnorm1(data, pi0 = result1$mix$pr, mu0 = result1$mix$pt, stdev)) - length(data)) /
    (1 - ddd)
}

findnpnorm1fc = function(data, val = -log(length(data)), stdev = 1, 
                         order = -4, ...){
  result1 = makenpnormfc(data, pi0 = 0, stdev = stdev, order = order)
  result2 = makenpnormfc(data, pi0 = 1, stdev = stdev, order = order)
  
  if (-2 * (result2$ll - result1$ll) + val < 0){
    ans1 = result2
  }else{
    # give a relative good estimate
    pix = c(dnpnorm1(0, pi0 = result1$mix$pr, mu0 = result1$mix$pt, stdev) / dnorm(0))
    ans2 = NR(pix, fn = findnpnormfc.fn, gr = findnpnormfc.gr, 
              lower = 0, upper = 1, data = data, pval = val + 2 * result1$ll, 
              stdev = stdev, order = order)
    ans1 =  makenpnormfc(data, pi0 = ans2$par, stdev = stdev, order = order)
  }
  ans1
}