Установка обратной функции

У меня есть функция, которая выглядит так: g(x) = f(x) - a^b/f(x)^b

g(x) - известная функция, предоставлен вектор данных.
f(x) - скрытый процесс.
a,b - параметры этой функции.

Из вышеизложенного получаем соотношение:
f(x) = inverse(g(x))

Моя цель — оптимизировать параметры a и b таким образом, чтобы f(x) была как можно ближе
к нормальному распределению. Если мы посмотрим на нормальный график f(x) Q-Q (прилагается), моя цель состоит в том, чтобы минимизировать расстояние между f(x) и прямой линией, представляющей нормальное распределение, путем оптимизации параметров a и < сильный>б.

Я написал следующий код:

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])  
}

У меня два вопроса:

  1. Как передать параметры a и b в функции?
  2. Как выполнить эту задачу оптимизации, то есть найти такие a и b, чтобы f(x) аппроксимировало нормальное распределение.

График нормы Q-Q f(x)


r function optimization inverse
person Paul    schedule 28.11.2013    source источник
comment
Хорошо, напишите, что такое G(x) (обратная функция). Это не будет f(x), так что ваша проблема неправильно сформулирована.   -  person Carl Witthoft    schedule 28.11.2013
comment
Карл – Например, g(x) – это дневная температура измерения в Нью-Йорке. Предположим, что мы можем объяснить это измерение по формуле g(x) = f(x) - a^b/f(x)^b. f(x) — неизвестный/скрытый процесс в измерении. f(x) = inverse(g(x)), которая может быть вычислена численно с помощью функции f_inverse для каждого значения g(x) (a и b известны).   -  person Paul    schedule 28.11.2013
comment
Не знаю, будет это иметь значение или нет, но поскольку a и f(x) оба возводятся в квадрат на b, не лучше ли было бы написать формулу как g(x) = f(x) - (a/f(x)) ^б? Таким образом, будет один вызов деления, затем один вызов экспоненты и, наконец, один вызов вычитания вместо 2 экспонент, 1 деления и 1 вычитания? Просто выбрасываю свои мысли туда.   -  person John Odom    schedule 03.12.2013

Ответы (1)


arrow_upward
2
arrow_downward

Не уверен, как вы смогли создать график QQ, поскольку предоставленные вами примеры не работают. Вы не указываете значения a и b и определяете f_func, но вызываете f_fun. В любом случае, вот мой ответ на ваши вопросы:

  1. Как передать параметры a и b в функции? - Просто передайте их в качестве аргументов функциям.
  2. Как выполнить эту оптимизационную задачу, то есть найти такие 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 и продолжать поиск, пока не найдете p-значение ниже 0,05 (учтите, что решения может и не быть)

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)

Эта последняя строка не работает, но без надлежащих данных и контекста это все, что я могу сделать. Надеюсь это поможет.

person Rohit Das    schedule 05.12.2013
comment
Рохит, возьмем f[1] = 0,9918427 и применим его к g_fun, когда a=1 и b=1 (ваш пример). g_fun = x - 1 ^ 1 / x ^ 1 = 0,9918427 - 1/0,9918427 = -0,01638169, но я должен получить -0,016339. Почему у нас такая разница? - person Paul; 05.12.2013
comment
Что вы имеете в виду, говоря, что вы должны получить -0,016339? - person Rohit Das; 06.01.2014