使用 R 评估重叠分布的积分误差

Integration error using R to evaluate overlapping distributions

此问题是对原始线程 (https://stats.stackexchange.com/questions/12209/percentage-of-overlapping-regions-of-two-normal-distributions)

的跟进

我修改了上面线程中的原始代码以标记图,但仅此而已。我 运行 遇到了使用一组输入而不是另一组输入的代码问题。作为 R 初学者,我只是在寻求帮助以了解原因。

mu1 <- 5
sd1 <- 1

mu2 <- 3
sd2 <- .75

这个数据集有效。

当我更改数字时,我遇到了一个问题:

mu1 <- .3439
sd1 <- .0005

mu2 <- .3420
sd2 <- .00075

明显有重叠,但积分不是计算。

当数字改变时,积分如何失败?我该如何补救?我在 xs 中没有足够的 x 点来计算积分吗?

完整代码如下:

min.f1f2 <- function(x, mu1, mu2, sd1, sd2) {
  f1 <- dnorm(x, mean=mu1, sd=sd1)
  f2 <- dnorm(x, mean=mu2, sd=sd2)
  pmin(f1, f2)
}


mu1 <- .3439
sd1 <- .0005

mu2 <- .3420
sd2 <- .00075

xs <- seq(min(mu1 - 3*sd1, mu2 - 3*sd2), max(mu1 + 3*sd1, mu2 + 3*sd2), .00001)
f1 <- dnorm(xs, mean=mu1, sd=sd1)
f2 <- dnorm(xs, mean=mu2, sd=sd2)

plot(xs, f1, type="l", ylim=c(0, max(f1,f2)), ylab="density")
lines(xs, f2, lty="dotted")
ys <- min.f1f2(xs, mu1=mu1, mu2=mu2, sd1=sd1, sd2=sd2)
xs <- c(xs, xs[1])
ys <- c(ys, ys[1])
polygon(xs, ys, col="red")


### integrate to find % overlap
iP <- integrate(min.f1f2, -Inf, Inf, mu1=mu1, mu2=mu2, sd1=sd1, sd2=sd2)
VaLue <- iP$value
VaLue <- sprintf("%.1f %%", 100*VaLue)

#percentage on plot (in middle of overlap.. half of ys height )
text(((mu2 + 3*sd2)+(mu1 - 3*sd1))/2, max(ys)/2, VaLue)
#label f1
text(mu1+sd1,max(f1),"HOLE")
#label f2
text(mu2+sd2,max(f2),"PEG")

我发现无论出于何种原因,当接近无穷大(或什至远高于或低于 3 sigma 尾部...)时,小数字的积分都会失败。

我只更新了积分的限制,从 -inf, inf 到 min(mu1 - 3*sd1, mu2 - 3*sd2), max(mu1 + 3*sd1, mu2 + 3*sd2) 和我再有一个近似值。当然,它并不完美,但我真的只关心我的应用程序的几个数字。

真的很想更好地理解为什么只有比例发生变化时这首先失败了。有什么想法吗?

更新了下面的代码。

min.f1f2 <- function(x, mu1, mu2, sd1, sd2) {
  f1 <- dnorm(x, mean=mu1, sd=sd1)
  f2 <- dnorm(x, mean=mu2, sd=sd2)
  pmin(f1, f2)
}

#HOLE 
mu1 <- .3439
sd1 <- .0005
#Peg
mu2 <- .3420
sd2 <- .00075

xs <- seq(min(mu1 - 3*sd1, mu2 - 3*sd2), max(mu1 + 3*sd1, mu2 + 3*sd2), .00001)
f1 <- dnorm(xs, mean=mu1, sd=sd1)
f2 <- dnorm(xs, mean=mu2, sd=sd2)

plot(xs, f1, type="l", ylim=c(0, max(f1,f2)), ylab="density")
lines(xs, f2, lty="dotted")
ys <- min.f1f2(xs, mu1=mu1, mu2=mu2, sd1=sd1, sd2=sd2)
xs <- c(xs, xs[1])
ys <- c(ys, ys[1])
polygon(xs, ys, col="red")


### integrate to find % overlap
iP <- integrate(min.f1f2, min(mu1 - 3*sd1, mu2 - 3*sd2), max(mu1 + 3*sd1, mu2 + 3*sd2), mu1=mu1, mu2=mu2, sd1=sd1, sd2=sd2)
VaLue <- iP$value
VaLue <- sprintf("%.1f %%", 100*VaLue)

#percentage on plot (in middle of overlap.. )
text(((mu2 + 3*sd2)+(mu1 - 3*sd1))/2, max(ys)/3.5, VaLue)
#label f1
text(mu1+sd1,max(f1),"HOLE")
#label f2
text(mu2+sd2,max(f2),"PEG")