r - 拟合反函数

Question

我有一个函数，它看起来像:
g(x) = f(x) - a^b/f(x)^b

g(x) - 已知函数，提供数据向量。
f(x) - 隐藏进程。
a,b - 此函数的参数。

从上面我们得到关系:
f(x) = 逆(g(x))

我的目标是优化参数 和 b 使得 f(x) 尽可能接近
为正态分布。如果我们查看 f(x) Q-Q 正态图(附后)，我的目的是通过优化参数 来最小化 f(x) 到代表正态分布的直线之间的距离和 b .

我写了下面的代码:

g_fun <- function(x) {x - a^b/x^b}

inverse = function (f, lower = 0, upper = 2000) {
      function (y) uniroot((function (x) f(x) - y), lower = lower, upper = upper)[1]
}


f_func = inverse(function(x) g_fun(x))
enter code here

# let's made up an example 
# g(x) values are known 
g <- c(-0.016339, 0.029646, -0.0255258, 0.003352, -0.053258, -0.018971, 0.005172,  
       0.067114, 0.026415, 0.051062)  

# Calculate f(x) by using the inverse of g(x), when a=a0 and b=b0
for (i in 1:10) {  
  f[i] <- f_fun(g[i])  
}

我有两个问题:

如何将参数 a 和 b 传递给函数？

如何执行这个优化任务，意味着找到 a 和 b 使得 f(x) 接近正态分布。

Answer 1

不确定您是如何生成 Q-Q 图的，因为您提供的示例不起作用。您没有指定 a 和 b 的值，而是定义 f_func 而是调用 f_fun。无论如何，这是我对您问题的回答:

如何将参数 a 和 b 传递给函数？ - 只需将它们作为
函数的参数。

如何执行这个优化任务，意味着找到 a 和 b 使得 f(x) 接近正态分布？ - 以相同的方式完成任何优化任务。定义一个成本函数，然后将其最小化。

这是修改后的代码:我添加了 a 和 b 作为参数，删除了反函数并将其合并到 f_func 中，它现在可以接受向量输入，因此不需要 for 循环。

g_fun <- function(x,a,b) {x - a^b/x^b}

f_func = function(y,a,b,lower = 0, upper = 2000){
  sapply(y,function(z) { uniroot(function(x) g_fun(x,a,b) - z, lower = lower, upper = upper)$root})
} 

# g(x) values are known 
g <- c(-0.016339, 0.029646, -0.0255258, 0.003352, -0.053258, -0.018971, 0.005172,  
       0.067114, 0.026415, 0.051062)  
f <- f_func(g,1,1) # using a = 1 and b = 1
#[1] 0.9918427 1.0149329 0.9873386 1.0016774 0.9737270 0.9905320 1.0025893
#[8] 1.0341199 1.0132947 1.0258569

f_func(g,2,10)
 [1] 1.876408 1.880554 1.875578 1.878138 1.873094 1.876170 1.878304 1.884049
 [9] 1.880256 1.882544

现在对于优化部分，它取决于您所说的 f(x) 近似正态分布的含义。如果需要，您可以比较 qq 行的均方误差。另外，既然你说近似，那么接近多少就足够了？您可以使用 shapiro.test 并继续搜索，直到找到低于 0.05 的 p 值(请注意可能没有解决方案)

shapiro.test(f_func(g,1,2))$p
[1] 0.9484821

cost <- function(x,y) shapiro.test(f_func(g,x,y))$p

现在我们有了一个成本函数，我们如何去最小化它。有许多不同的方法可以进行数值优化。看看优化函数 http://stat.ethz.ch/R-manual/R-patched/library/stats/html/optim.html .

optim(c(1,1),cost)

最后一行不起作用，但没有适当的数据和上下文，这是我所能做的。希望这可以帮助。