r - 自动求解一个 `pt`的方程为 `ncp`

Question

我想知道是否可以有效地改变下面代码中的ncp，使x变成.025 和 .975(在舍入误差内)。

x <- pt(q = 5, df = 19, ncp = ?)

------------

澄清

q = 5 和 df = 19(以上)只是两个假设的数字，所以 q 和 df 可以是任何其他两个数字。我期望的是一个函数/例程，它将 q 和 df 作为输入。

Answer 1

uniroot 有什么问题？

f <- function (ncp, alpha) pt(q = 5, df = 19, ncp = ncp) - alpha

par(mfrow = c(1,2))
curve(f(ncp, 0.025), from = 5, to = 10, xname = "ncp", main = "0.025")
abline(h = 0)
curve(f(ncp, 0.975), from = 0, to = 5, xname = "ncp", main = "0.975")
abline(h = 0)

所以对于0.025的情况，根在(7, 8)；对于 0.975 情况，根位于 (2, 3)。

uniroot(f, c(7, 8), alpha = 0.025)$root
#[1] 7.476482

uniroot(f, c(2, 3), alpha = 0.975)$root
#[1] 2.443316

------------

(经过一些讨论...)

好的，现在我明白你的最终目标了。您希望将此方程求解器实现为一个函数，输入为 q 和 df。所以他们是未知的，但固定的。他们可能来自实验。

理想情况下，如果存在解析解，即 ncp 可以写成根据 q、df 和 的公式>alpha，那太好了。然而，这对于 t 分布是不可能的。

数值解是方法，但 uniroot 不是用于此目的的好选择，因为它依赖于“plot - view - guess - specification”。 loki 的回答也很粗糙，但有一些改进。这是一个网格搜索，具有固定的步长。从接近 0 的值开始，比如 0.001，然后增加该值并检查近似误差。当此错误无法减少时，我们停止。

这真正引发了使用牛顿法或拟牛顿法进行数值优化的想法。在一维情况下，我们可以使用函数 optimize。它在搜索中采用变步长搜索，因此比固定步长搜索收敛速度更快。

让我们将函数定义为:

ncp_solver <- function (alpha, q, df) {
  ## objective function: we minimize squared approximation error
  obj_fun <- function (ncp, alpha = alpha, q = q, df = df) {
    (pt(q = q, df = df, ncp = ncp) - alpha) ^ 2
    }
  ## now we call `optimize`
  oo <- optimize(obj_fun, interval = c(-37.62, 37.62), alpha = alpha, q = q, df = df)
  ## post processing
  oo <- unlist(oo, use.names = FALSE)  ## list to numerical vector
  oo[2] <- sqrt(oo[2])  ## squared error to absolute error
  ## return
  setNames(oo, c("ncp", "abs.error"))
  }

请注意，-37.62/37.62 被选为 ncp 的下限/上限，因为它是 t 分布支持的最大值在 R 中(阅读 ?dt)。

例如，让我们试试这个功能。如果您在问题中给出了 q = 5 和 df = 19:

ncp_solver(alpha = 0.025, q = 5, df = 19)
#         ncp    abs.error
#7.476472e+00 1.251142e-07 

结果是一个命名向量，带有ncp 和绝对逼近误差。

同样我们可以这样做:

ncp_solver(alpha = 0.975, q = 5, df = 19)
#         ncp    abs.error
#2.443347e+00 7.221928e-07 

------------

跟进

Is it possible that in the function ncp_solver(), alpha takes a c(.025, .975) together?

为什么不把它包装起来进行“矢量化”:

sapply(c(0.025, 0.975), ncp_solver, q = 5, df = 19)

#                  [,1]         [,2]
#ncp       7.476472e+00 2.443347e+00
#abs.error 1.251142e-07 7.221928e-07

How come 0.025 gives upper bound of confidence interval, while 0.975 gives lower bound of confidence interval? Should this relationship reversed?

不足为奇。默认情况下，pt 计算较低的尾部概率。如果您想要“正确”的关系，请在 pt 中设置 lower.tail = FALSE:

ncp_solver <- function (alpha, q, df) {
  ## objective function: we minimize squared approximation error
  obj_fun <- function (ncp, alpha = alpha, q = q, df = df) {
    (pt(q = q, df = df, ncp = ncp, lower.tail = FALSE) - alpha) ^ 2
    }
  ## now we call `optimize`
  oo <- optimize(obj_fun, interval = c(-37.62, 37.62), alpha = alpha, q = q, df = df)
  ## post processing
  oo <- unlist(oo, use.names = FALSE)  ## list to numerical vector
  oo[2] <- sqrt(oo[2])  ## squared error to absolute error
  ## return
  setNames(oo, c("ncp", "abs.error"))
  }

现在你看到了:

ncp_solver(0.025, 5, 19)[[1]]  ## use "[[" not "[" to drop name
#[1] 2.443316

ncp_solver(0.975, 5, 19)[[1]]
#[1] 7.476492

--------

错误报告和修复

有人报告说上面的 ncp_solver 不稳定。例如:

ncp_solver(alpha = 0.025, q = 0, df = 98)
#      ncp abs.error 
#-8.880922  0.025000 

但另一方面，如果我们在这里用 uniroot 仔细检查:

f <- function (ncp, alpha) pt(q = 0, df = 98, ncp = ncp, lower.tail = FALSE) - alpha
curve(f(ncp, 0.025), from = -3, to = 0, xname = "ncp"); abline(h = 0)

uniroot(f, c(-2, -1.5), 0.025)$root
#[1] -1.959961

所以 ncp_solver 显然有问题。

事实证明我们不能使用太大的界限，c(-37.62, 37.62)。如果我们将它缩小到 c(-35, 35)，就没问题了。

此外，为了避免容差问题，我们可以将目标函数从平方误差更改为绝对误差:

ncp_solver <- function (alpha, q, df) {
  ## objective function: we minimize absolute approximation error
  obj_fun <- function (ncp, alpha = alpha, q = q, df = df) {
    abs(pt(q = q, df = df, ncp = ncp, lower.tail = FALSE) - alpha)
    }
  ## now we call `optimize`
  oo <- optimize(obj_fun, interval = c(-35, 35), alpha = alpha, q = q, df = df)
  ## post processing and return
  oo <- unlist(oo, use.names = FALSE)  ## list to numerical vector
  setNames(oo, c("ncp", "abs.error"))
  }

ncp_solver(alpha = 0.025, q = 0, df = 98)
#          ncp     abs.error 
#-1.959980e+00  9.190327e-07 

该死的，这是一个非常烦人的错误。但现在放松一下。

从pt获取警告信息的报告

我还从 pt 收到一些关于恼人警告消息的报告:

ncp_solver(0.025, -5, 19)
#          ncp     abs.error 
#-7.476488e+00  5.760562e-07
#Warning message:
#In pt(q = q, df = df, ncp = ncp, lower.tail = FALSE) :
#  full precision may not have been achieved in 'pnt{final}'

我不太清楚这里发生了什么，但同时我没有观察到误导性的结果。因此，我决定使用 suppressWarnings 抑制来自 pt 的警告:

ncp_solver <- function (alpha, q, df) {
  ## objective function: we minimize absolute approximation error
  obj_fun <- function (ncp, alpha = alpha, q = q, df = df) {
    abs(suppressWarnings(pt(q = q, df = df, ncp = ncp, lower.tail = FALSE)) - alpha)
    }
  ## now we call `optimize`
  oo <- optimize(obj_fun, interval = c(-35, 35), alpha = alpha, q = q, df = df)
  ## post processing and return
  oo <- unlist(oo, use.names = FALSE)  ## list to numerical vector
  setNames(oo, c("ncp", "abs.error"))
  }

ncp_solver(0.025, -5, 19)
#          ncp     abs.error 
#-7.476488e+00  5.760562e-07

好的，现在安静。