10

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

问：假设这些频率点来自 Weibull 分布T，如何将最佳 Weibull 密度函数拟合到这些点，从而从中推断出分布T参数？ 

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)

在此处输入图像描述

并fitdistr用于set.values：

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

为什么我认为这是不正确的方式，为什么我要寻找更好的解决方案R？

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

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

4

3 回答 3

3

首先尝试所有点

第二次尝试，第一点下降 这是一个更好的尝试，就像之前它用于optim查找限制在框中的一组值的最佳值（由调用中的lowerandupper向量定义optim）。请注意，除了 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)
于 2015-05-03T10:33:36.337 回答
3

您可以直接计算最大似然参数，如此处所述。

# 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))
于 2015-05-06T07:43:02.180 回答
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)
于 2015-05-03T09:37:35.257 回答