log(0) 的调整函数

Question

我为泊松回归编写了一个函数。数据集发现有一些计数数据，其中 5 个条目是 y = 0。我想计算偏差残差，根据我的函数中的公式：

devianceResiduals <- sign(y - fittedValuesFullModell) * sqrt(2 * y *
                           log(y / fittedValuesFullModell) - 2 *
                            (y - fittedValuesFullModell))

我的问题是我得到 NaN，因为 log(y = 0) = -inf。因此尝试编写一个循环，使用 2 个不同的论坛来计算偏差残差。如果 y = 0 由于日志规则，公式简化为 - 2 * (y - compatibleValuesFullModell))。 如果 y > 0 我想使用我的常规公式。我如何调整我在计算偏差残差时没有得到 NaN 的函数？ （德国cmets对不起）

data(discoveries)
disc <- data.frame(count=as.numeric(discoveries),
                   year=seq(0,(length(discoveries)-1),1))

yearSqr <- disc$year^2


poisMod<-function(formula, data){

  # Definiere Regressionsformel
  form <- formula(formula)

  # DataFrame wird erzeugt
  model <- model.frame(formula, data = data)

  #  Designmatrix wird erzeugt
  x <- model.matrix(formula, data = data)

  # Response Variable erzeugt
  y <- model.response(model)

  # Designmatrix f?r "Nullmodell-Sch?tzung" (nur Intercept wird gesch?tzt)
  # wird erzeugt
  xNull <- as.matrix(x[, 1])

  # Parameteranzahl die es zu sch?tzen gilt im "Nullmodell"
  # (nur Intercept gesch?tzt)
  parNullModell <- rep(0, 1)

  # Parameteranzahl die es zu sch?tzen gilt im im "Fullmodell" (Fullmodell
  # hei?t, dass Koeffizienten f?r jede erkl?rende Variable aus der Formula
  # gesch?tzt werden)
  parFullModell <- rep(0, ncol(x))

  # Funktionaufruf wird gespeichert
  call <- match.call()

  # Maximierung der Loglikelihood durch Sch?tzung der Koeffzienten, Hessematrix
  # soll ebenfalls mit angegeben werden
  optimierung <- optim(par =  parFullModell, fn = logLike,
                       x = x, y = y, hessian = TRUE)
  # Hessematix aus Optimierung
  hesseMatrix <- optimierung$hessian

  # Koeffizienten aus der Sch?tzung des "Fullmodells"
  koefFullModell <- round(optimierung$par, 4)

  # Koeffizienten aus der Sch?tzung des "Nullmodells"
  koefNullModell <- round(optim(par = parNullModell, fn = logLike,x = xNull,
                                y = y, method = "Brent", lower = -1000 ,
                                upper = 1000)$par, 4)

  # Bennenung der Koeffizienten entsprechend ihrer erkl?renden Variablen
  names(koefFullModell) <- colnames(x)

  # Loglikelihood des "FUllmodells" bestimmt
  logLikeFullModell <- logLike(y, x, koefFullModell)

  # Loglikelihood des "NullModells"
  logLikeNullModell <- logLike(y, xNull, koefNullModell)

  # AIC Kriterium f?r das "FullModell"
  aic <- round(2 * (logLikeFullModell) + 2 * length(koefFullModell), 1)

  # Degrees of freedom vom "Fullmodell"
  dofFullModell <- nrow(x) - ncol(x)

  #Degrees of freedom vom "Nullmodell"
  dofNullModell<- nrow(xNull) - ncol(xNull)

  # Fitted values des "Nullmodells"
  fittedValuesNullModell <- exp(xNull %*% koefNullModell)

  # Fitted values des "Fullmodells"
  fittedValuesFullModell <- exp(x %*% koefFullModell)

  # Residual deviance bestimmt
  residualDeviance <- round(2 * sum(y * log(y / fittedValuesFullModell) -
                           (y - fittedValuesFullModell),na.rm = TRUE), 1)

  # Null deviance bestimmt
  nullDeviance <- round(2 * sum(y * log(y / fittedValuesNullModell) -
                       (y - fittedValuesNullModell),na.rm = TRUE), 1)

  # Deviance residuals bestimmt
  devianceResiduals <- sign(y - fittedValuesFullModell) * sqrt(2 * y *
                           log(y / fittedValuesFullModell) - 2 *
                            (y - fittedValuesFullModell))

  # y und fitted values des "Fullmodells" in dataframe gepackt, damit ggplo2
  # sp?ter darauf angwendet werden kann, da ggplot2 nicht auf diese Klasse
  # poisMod angewendet werden kann
  ggPlotData <- data.frame(obs = y, fitted = fittedValuesFullModell)

  # Liste mit allen Kennzahlen
  result <- list("coefficients" = koefFullModell,
                 "call" = call,
                 "aic" = aic,
                 "dofNullModell" = dofNullModell,
                 "dofFullModell" = dofFullModell,
                 "x" = x,
                 "y" = y,
                 "fittedValuesNullModell" = fittedValuesNullModell,
                 "fittedValuesFullModell" = fittedValuesFullModell,
                 "residualDeviance" = residualDeviance,
                 "nullDeviance" = nullDeviance,
                 "hesseMatrix" = hesseMatrix,
                 "optimierung" = optimierung,
                 "devianceResiduals" = devianceResiduals,
                 "ggPlotData" = ggPlotData,
                 "model" = model)

  # poisMod Klasse erzeugt
  class(result) <- "poisMod"

  # Result liste ausgegeben
  return(result)
}


p2 <- poisMod(count~year+yearSqr,data=disc)

Answer 1

这不是一个解决方案，只是一个评论。 基本上是log(0) = Inf。所以你可以定义你想要的无限大/小。例如，您可以将 0 替换为 0.000000001 或更小的值，该值接近于零。比如log(1E-99) = -227.9559