2

我想对我的数据是否符合特定的分布函数进行视觉评估。为此，我使用 R 生成分位数-分位数 (QQ) 图。分布函数非常具体，不在概率分布的标准列表中，所以我编写了自己的 R 函数来描述它。它在下面的代码中称为“DistFunc”，由两个伽马函数的比率组成。

简而言之，我在代码中所做的是从包含两列的文件“DistributionEstimate.txt”中读取数据。第 1 列是 x 值，第 2 列是 y 值。变量“a”和“b”是我之前在另一个程序中使用此分布函数与数据的最小二乘拟合确定的最佳拟合参数。然后我定义 DistFunc 并尝试使用 qqmath 函数绘制 QQ 图。

问题出现在这一点上。R 继续给我很多警告，说 DistFunc 返回的值超出了 'gammafn' 的范围，并且无法绘制任何内容。这很公平，因为我知道该函数包含一个靠近原点的极点。正如您在代码中看到的那样，我尝试规范化 DistFunc 以尝试将其转换为概率分布（我认为这是使用 qqmath 所必需的吗？），但是，这没有帮助。

你们中是否有人知道如何克服这个问题 - 例如，通过使用不需要归一化的不同绘图函数，或者将其转换为伪概率分布，而不会对结果造成太大影响？

我将非常感谢任何有用的输入！

install.packages('lattice')
library(lattice)
x<-read.table("C:/DistributionEstimate.txt", colClasses = c(rep("NULL",1),rep("numeric",1)), header = FALSE)
y<-read.table("C:/DistributionEstimate.txt", colClasses = c(rep("numeric",1),rep("NULL",1)), header = FALSE)
x<-sapply(x, as.numeric)
y<-sapply(y, as.numeric)
a<-16359727025.407821410;
b<-198838619.13262583836;
DistFunc <- function(k,ampl=a,stretch=b) {
    fdist<-ampl*gamma(k*stretch-1/2)/gamma(k*stretch+1)
    fnorm<-fdist/sum(fdist)
}
qqmath(DistFunc(x), y, col="blue", envelope=.95, xlab="Quantiles of the best-fit model", ylab="Quantiles of the data")
abline(0,1, col="red", lwd=2)
grid()
4

1 回答 1

3

QQ 图背后的想法是将被认为来自某个分布的观察值与您期望从相同大小的样本中从该分布中看到的值进行比较。

所以第一个问题是你同时拥有x和y价值观。QQ 图是单变量图。您正在将一组值与分布进​​行匹配。绘制(x,y)对的第二个维度由分布函数计算。

分布函数qqmath期望不是密度函数。它需要一个函数，将分位数转换为分布中的值。q*这与在 R 中工作的分布函数族相同，例如qnrom或qexp。(-Inf,Inf)该函数必须接受 0-1 范围内的数字，并将其转换为 forqnorm或(0, Inf)for分布域中的值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")

qqmath vc 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)
}

第一个参数是单参数密度函数。这将在集成期间使用。然后fromandto参数指定您的值在哪里具有非零概率。然后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

您会看到我们得到了非常相似的 0.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。

于 2014-05-26T01:33:26.557 回答