在R中给定条件的情况下，如何向后求解未知问题

Question

对于一组输入(见下文(，当我运行以下R代码时，我会得到存储在ncp中的两个答案。但我想知道df2应该是什么，使得ncp(即abs(ncp[2] - ncp[1])(中这两个答案之间的区别是.15？

那么，其他的都固定了，df2应该是什么，abs(ncp[2] - ncp[1]) = .15应该是什么？这可以在R中完成吗？

alpha = c(.025, .975); df1 = 3; peta = .3   # The input
f <- function(alpha, q, df1, df2, ncp){     # Notice `ncp` is the unknown
alpha - suppressWarnings(pf(q = (peta / df1) / ((1 - peta)/df2), df1, df2, ncp, lower = FALSE))
}
ncp <- function(df2){      # Root finding: finds 2 `ncp` for a given `df2`
b <- sapply(c(alpha[1], alpha[2]),
function(x) uniroot(f, c(0, 1e7), alpha = x, q = peta, df1 = df1, df2 = df2)[[1]])
b / (b + (df2 + 4))
}
# Example of use:
ncp(df2 = 108) # Answers: 0.1498627 0.4100823
# What should `df2` be so that the difference between 2 answers is `.15`

Answer 1

stats包中的optim和optimize函数是处理此类问题的好方法。optimize更简单，处理一个单变量的情况。以下是一个示例解决方案：

# Define a function that we want to minimise the output of
ncp_diff <- function(df2, target = 0.15){
the_ncp <- ncp(df2)
return(abs(abs(the_ncp[2] - the_ncp[1]) - target))
}
# try it out - how far out is df2 = 108 from our target value?
ncp_diff(108) # 0.1102
# find the value of ncp_diff's argument that minimises its return:
optimize(ncp_diff, interval = c(0, 1000)) # minimum = 336.3956
ncp(336.3956) # 0.218, 0.368

请注意，虽然336.3956是的解决方案，但它不一定是的解决方案。你需要注意局部极小值。然而，处理这个问题超出了简单答案的范围。