将分布拟合到 R 中的给定频率值

Posted

技术标签:

【中文标题】将分布拟合到 R 中的给定频率值【英文标题】:Fit distribution to given frequency values in R 【发布时间】:2015-05-17 05:39:45 【问题描述】:

我的频率值随时间变化(x 轴单位),如下图所示。经过一些归一化后,这些值可能被视为某些分布的密度函数的数据点。

问:假设这些频点来自 Weibull 分布T,我如何将最佳 Weibull 密度函数拟合到这些点从而推断分布@ 987654325@参数从哪里来的?

sample <- c(7787,3056,2359,1759,1819,1189,1077,1080,985,622,648,518,
            611,1037,727,489,432,371,1125,69,595,624)

plot(1:length(sample), sample, type = "l")
points(1:length(sample), sample)

更新。 为了防止被误解,我想补充一点解释。 我的频率值随时间变化(x 轴单位)我的意思是我有数据表明我有:

7787价值1的实现 3056 次价值 2 的实现 2359 个实现价值 3 ... 等。

实现我的目标的某种方式(我认为不正确)是创建一组这些实现:

# Loop to simulate values 
set.values <- c()
for(i in 1:length(sample))
  set.values <<- c(set.values, rep(i, times = sample[i]))


hist(set.values)
lines(1:length(sample), sample)
points(1:length(sample), sample)

并在set.values 上使用fitdistr

f2 <- fitdistr(set.values, 'weibull')
f2

为什么我认为这是不正确的方法以及为什么我要在R 中寻找更好的解决方案?

在上面介绍的分布拟合方法中,假设set.values 是我从分布T 实现的完整

在我原来的问题中,我知道密度曲线 第一部分的点 - 我不知道它的尾巴,我想 估计尾部(以及整个密度函数

【问题讨论】:

我已经用直方图更新了我的答案。 你知道密度曲线第一部分结束和尾部开始的确切值吗?您的样本以值 22 结束:我可以假设尾部从 23 开始吗? 恐怕我不明白(我不知道我可以在这里使用“分布尾部”的正式定义)。我的最终目标是计算分布T 的变量的期望值。也许有理由假设第一部分(上面直方图中 1. 和 2. 点之间的部分)是线性的,而后一部分 - Weibull(Weibull 是我从向我提供数据的人那里得到的假设。我不会我不会为此赌上我的性命,但我倾向于假设相同。) 你说:“在我原来的问题中,我知道密度曲线第一部分的点”。 “第一部分”到底是什么意思? “第一部分”在什么值处停止?您还说:“我不知道它的尾巴,我想估计尾巴(以及整个密度函数)”。为此,您需要(一个标准)选择尾部开始的位置。 我想我已经回答过了。我的解决方案在哪些方面不是您想要的? 【参考方案1】:

假设数据来自 Weibull 分布,您可以像这样得到形状和尺度参数的估计值:

sample <- c(7787,3056,2359,1759,1819,1189,1077,1080,985,622,648,518,
        611,1037,727,489,432,371,1125,69,595,624)
 f<-fitdistr(sample, 'weibull')
 f

如果你不确定它是否是分布式的 Weibull,我建议使用 ks.test。这将测试您的数据是否来自假设分布。鉴于您对数据性质的了解,您可以测试几个选定的分布,看看哪一个效果最好。

对于您的示例,这将如下所示:

 ks = ks.test(sample, "pweibull", shape=f$estimate[1], scale=f$estimate[2])
 ks

p 值不显着,因此您不会拒绝数据来自 Weibull 分布的假设。

更新:Weibull 或指数直方图看起来与您的数据非常匹配。我认为指数分布给你一个更好的拟合。帕累托分布是另一种选择。

f<-fitdistr(sample, 'weibull')
z<-rweibull(10000, shape= f$estimate[1],scale= f$estimate[2])
hist(z)

f<-fitdistr(sample, 'exponential')
z = rexp(10000, f$estimate[1]) 
hist(z)

【讨论】:

嗯,我承认这个答案是正确的,恐怕我犯了一个错误。 fitdistr 函数将值(此处:来自sample 向量的值)视为分布T 中的实现(换句话说:点绘制drom 分布T),而不是:某些分布的密度函数曲线的数据点。看到当我使用估计的 shapescale 参数从估计的 Tthen 绘制这些点的密度点时(不是的情况我的问题),我最终得到像 this 这样的密度,其中 x 轴值不正确。 您所说的“某种分布的密度函数曲线的数据点”是什么意思?在您的问题中,您说您认为是威布尔。 pdf 适用于具有估计参数的 Weibull。如果要将其与图表进行比较,则需要将其与 hist(sample) 进行比较。您上面的图表看起来不像 pdf。 嗨@TinaW,请参考我刚刚添加到我的问题中的更新。 是什么让你认为这是 Weibull 分布式? 我认为只有尾巴是。【参考方案2】:

这是一个更好的尝试,就像之前它使用optim 来找到限制在框中的一组值的最佳值(由optim 调用中的lowerupper 向量定义)。请注意,除了 Weibull 分布形状参数之外,它还缩放 x 和 y 作为优化的一部分,因此我们有 3 个参数需要优化。

不幸的是,当使用所有点时,它几乎总是在约束框的边缘找到一些东西,这向我表明 Weibull 可能并不适合所有数据。问题在于这两点——它们太大了。您会在第一个绘图中看到尝试拟合所有数据。

如果我放弃前两点并仅拟合其余部分,我们会得到更好的拟合。您可以在 第二个情节 中看到这一点。我认为这是一个很好的拟合,无论如何它是约束框内部的局部最小值。

library(optimx)
sample <- c(60953,7787,3056,2359,1759,1819,1189,1077,1080,985,622,648,518,
            611,1037,727,489,432,371,1125,69,595,624)
t.sample <- 0:22

s.fit <- sample[3:23]
t.fit <- t.sample[3:23]

wx <- function(param)  
  res <- param[2]*dweibull(t.fit*param[3],shape=param[1])
  return(res)
 
minwx <- function(param)
  v <- s.fit-wx(param)
  sqrt(sum(v*v))


p0 <- c(1,200,1/20)
paramopt <- optim(p0,minwx,gr=NULL,lower=c(0.1,100,0.01),upper=c(1.1,5000,1))

popt <- paramopt$par
popt
rms <- paramopt$value
tit <- sprintf("Weibull - Shape:%.3f xscale:%.1f  yscale:%.5f rms:%.1f",popt[1],popt[2],popt[3],rms)

plot(t.sample[2:23], sample[2:23], type = "p",col="darkred")
lines(t.fit, wx(popt),col="blue")
title(main=tit)

【讨论】:

嗨@Mike Wise,感谢您的关注和这个完整的例子!正如您所看到的,通过这种方式很难拟合曲线 - 在我看来,拟合的曲线不能很好地拟合,因为它不够“弯曲”。我相信它应该更像来自here 的蓝色环,不是吗? 哇,我刚刚意识到我认为只有尾巴是威布尔可能是一个很好的观点!谢谢你。我会在几天内进一步调查它和你的解决方案。 我还有一些想法,明天或今晚可能会尝试一下。 试图一次拟合两个 Weibull 来处理前两个点,但无法收敛。 您可以通过稍微改变 x 和 y 比例来获得其他合适的效果。了解更多关于时间尺度(起源是什么等)会很有帮助。如果这是我的项目,我可能会在这些拟合上进行引导以获得参数范围和分布。【参考方案3】:

您可以直接计算最大似然参数,如here所述。

# Defining the error of the implicit function
k.diff <- function(k, vec)
  x2 <- seq(length(vec))
  abs(k^-1+weighted.mean(log(x2), w = sample)-weighted.mean(log(x2), 
                                                            w = x2^k*sample))


# Setting the error to "quite zero", fulfilling the equation
k <- optimize(k.diff, vec=sample, interval=c(0.1,5), tol=10^-7)$min

# Calculate lambda, given k
l <- weighted.mean(seq(length(sample))^k, w = sample)

# Plot
plot(density(rep(seq(length(sample)),sample)))
x <- 1:25
lines(x, dweibull(x, shape=k, scale= l))

【讨论】:

在我运行我的代码之前它不起作用。不知道为什么。错误消息是:k 我收到错误消息:as.double(w) 中的错误:无法将“闭包”类型强制转换为“双”类型的向量 您好@user1965813,谢谢您的回答!我能够重现您的代码。我还复制了删除第一个元素的示例代码(因为在讨论中,有人认为第一个点不“适合”其余部分,我倾向于这种想法)see here。然后我比较了these dendisty plots 的形状,似乎 Mike 的解决方案在这种情况下给出了更合适的结果。不过,非常感谢您分享这种方法!

以上是关于将分布拟合到 R 中的给定频率值的主要内容,如果未能解决你的问题,请参考以下文章

将 Weibull 累积分布拟合到 R 中的质量传递数据

如何将频率分布转换为R中的概率分布

R中的对数对数概率图

R语言分布的卡方拟合优度检验

拟合 3 参数 Weibull 分布

使用 R 中的矩法拟合 Weibull