Как построить график Q-Q данных по сравнению с пользовательской теоретической функцией

Идея графика QQ состоит в том, чтобы сравнить наблюдения, которые, как считается, формируют определенное распределение, со значениями, которые вы ожидаете увидеть из этого распределения в выборке того же размера.

Итак, первая проблема заключается в том, что у вас есть значения x и y. QQ-график является одномерным графиком. Вы сопоставляете один набор значений с распределением. Второе измерение для построения (x,y) пар вычисляется функцией распределения.

Функция распределения, которую ожидает qqmath, не является функцией плотности. Нужна функция, которая будет превращать квантили в значения из распределения. Это то же самое, что семейство функций распределения q* работает в R, например qnrom или qexp. Функция должна принимать число в диапазоне 0-1 и преобразовывать его в значение в домене распределения (-Inf,Inf) для qnorm или (0, Inf) для qexp. Во время построения qqmath передаст этой функции список квантилей и получит список ожидаемых значений. Затем он построит список ожидаемых значений по сравнению с (отсортированными) наблюдаемыми значениями.

В качестве примера я просто собираюсь использовать функцию qexp в качестве «пользовательской» квантильной функции. Обратите внимание, что

myDist<-function(x) {
    qexp(x, 5)
}

set.seed(15)
x <- rexp(100, 5)
qqmath(~x, distribution=myDist, main="qqmath")

И это точно так же, как

exp.x <- myDist(ppoints(length(x)))
xyplot(sort(x)~exp.x, main="xyplot")

Я думаю, что одна из ваших проблем заключается в том, что DistFunc больше похоже на плотность, чем на функцию квантиля. Чтобы перейти от функции плотности к вероятностям, вы должны интегрировать. Вот вспомогательная функция, чтобы попытаться создать функцию q-like для произвольной функции плотности.

getq <- function(density, from, to, steps=1000) {
    x <- seq(from=from, to=to, length.out=steps) 
    y <- mapply(function(a,b)integrate(density,a,b)$value, x[-steps], x[-1])
    approxfun(c(0,cumsum(y)),x)
}

Первый параметр представляет собой однопараметрическую функцию плотности. Это будет использоваться во время интеграции. Затем параметры from и to указывают, где ваши значения имеют ненулевые вероятности. Тогда steps — это количество точек, где мы будем выполнять интегрирование. Затем мы используем approxfun для интерполяции между фактически рассчитанным количеством точек и точкой, запрошенной последней функцией q. Давайте посмотрим, как это работает со стандартной плотностью. Снова мы будем использовать экспоненциальную, скорость 5, плотность

myq <- getq(function(x) dexp(x,5), 0, 4)

Обратите внимание, что мы создаем анонимную функцию, чтобы обернуть dexp параметром скорости, поэтому наша плотность принимает только один параметр. Здесь мы просто переходим от 0 к 4, потому что к этому моменту мы почти достигаем вероятности 1,0. Теперь мы можем использовать эту функцию как стандартную qexp

> qexp(.5,5)
[1] 0.1386294
> myq(.5)
[1] 0.1386388

Вы видите, что мы получаем очень похожие ответы для .5. Так что, похоже, это работает. Итак, это один из быстрых способов преобразовать функцию плотности в функцию квантилей, если ваша функция квантилей не имеет красивой закрытой формы.

И последняя проблема, которую я вижу, заключается в том, что ваши значения a и b огромны. Использование их внутри функции gamma быстро приведет к числам, с которыми R не справится. Теперь вы делите одно gamma на другое, поэтому есть надежда, что они несколько компенсируются, но обычно вы сталкиваетесь с переполнением, используя стандартные версии. Таким образом, хитрость заключается в том, чтобы вычислять большие значения в логарифмическом масштабе, а затем exp(), когда вы все сделали, чтобы вернуться к естественному масштабу. Таким образом, вы можете изменить свою функцию на

DistFunc <- function(k,ampl=a,stretch=b) {
    fdist <- exp(log(ampl) + lgamma(k*stretch-1/2) - lgamma(k*stretch+1))
    fnorm <- fdist/sum(fdist)
}

Обратите внимание, что lgamma — это гамма-функция в логарифмическом масштабе. Но с вашими значениями a и b даже этого в большинстве случаев недостаточно. Я не уверен, как вы можете использовать числа из этой функции с учетом ваших параметров. Я также не уверен, каков, по вашему мнению, диапазон вашего дистрибутива. Я не мог найти способ сделать это интегрированным с 1, как это должно быть в хорошей функции плотности.