使用阶跃函数进行最大似然估计

Question

我想对一些数据拟合一个阶跃函数（两个参数）。下面的代码无法完成这项工作。我想知道是否 round() 参数是问题所在。然而，我也试着对参数进行划分，使参数的变化很小（如0.001），以引起显著的变化。但这并没有改变拟合的结果。有什么办法可以将这个函数正确地拟合到数据上吗？

dat <- c(rbinom(100, 100, 0.95), rbinom(50, 100, 0.01), rbinom(100, 100, 0.95))

plot(dat/100)

stepFnc <- function(parms, t) {
  par <- as.list(parms)
  (c(rep(1-(1e-5), par$t1), rep(1e-5, par$t2), rep(1-(1e-5), t)))[1:t]
}

lines(stepFnc(c(t1 = 50, t2 = 50), length(dat)))

loglik <- function(t1 = 50, t2 = 50) {
    fit <- snowStepCurve(parms = list(t1=round(t1,0), t2=round(t2,0)), t = length(dat))
    lines(fit)
    -sum(dbinom(x = dat, size = 100, prob = fit, log = T), na.rm = T)
}

mle <- bbmle::mle2(loglik)
mle@coef

lines(snowStepCurve(mle@coef, length(dat)), lwd = 2, lty = 2, col = "orange")

Answer 1

对于离散的x数据，我会用蛮力的方法。

x <- seq_along(dat)

foo <- function(x, lwr, upr) {
  y <- x
  y[x <= lwr | x > upr] <- mean(dat[x <= lwr | x > upr])
  y[x > lwr & x <= upr] <- mean(dat[x > lwr & x <= upr])
  y
}

SSE <- function(lwr, upr) {
  sum((dat - foo(x, lwr, upr))^ 2) 
}

limits <- expand.grid(lwr = x, upr = x)
limits <- limits[limits$lwr <= limits$upr,]
nrow(limits)

SSEvals <- mapply(SSE, limits$lwr, limits$upr)

id <- which(SSEvals == min(SSEvals))
optlims <- limits[id,]
meanouter <-  mean(dat[x <= optlims$lwr | x > optlims$upr])
meaninner <- mean(dat[x > optlims$lwr & x <= optlims$upr])

bar <- function(x) {
  y <- x
  y[x <= optlims$lwr | x > optlims$upr] <- meanouter
  y[x > optlims$lwr & x <= optlims$upr] <- meaninner
  y
}

plot(dat/100)
curve(bar(x) / 100, add = TRUE)