首页 > 解决方案 > 用负比例拟合 log pearson III 的问题

问题描述

我想对我在 R 中的一些数据点执行 log Pearson III 拟合。我遵循了这个指南 链接。但是当我的偏度 (g) 为负时我遇到了一个问题(当然参数“scale”也是负的,因为在 scale 的计算中使用了“sign(g)”)。“FAdist”的分布不适用于负比例参数,我需要它作为“fitdist”(来自 fitdistrplus)的起始值。有些页面说参数形状和比例仅在皮尔逊三世(或广义伽马分布)中是正的,而其他的则不是,我已经没有想法了。代码是:

library(fitdistrplus)
library(FAdist)
library(e1071) 
#data
df <-(92.8, 53.2, 112.0, 164.0, 132.0,  69.9, 140.0,  48.3, 123.0 ,24.6, 179.0,  55.1,  31.3,  17.0, 111.0,  35.4, 133.0, 505.0, 303.0, 121.5, 203.0, 198.0, 250.0, 232.0, 185.0, 222.0, 191.0, 238.0,  53.0, 121.0, 106.4, 347.3, 186.4,  89.1, 131.4 ,53.2, 252.6)
# log of df
df<-log(df)
#Pearson 3 Sample moments
 m <- mean(df)
v <- var(df)
s <- sd(df)
g <- e1071::skewness(df, type=1)
n <- length(df)
#Pearson 3 Parameter estimation
my.shape <- (2/g)^2
my.scale <- sqrt(v)/sqrt(my.shape)*sign(g) # modified as recommended by Carl Schwarz, this is negative =(
my.thres <- m-(my.shape*my.scale)
# All parameter together
my.param <- list(shape=my.shape, scale=my.scale, thres=my.thres )
# fit dist from the "fitdistrplus" and "lgamma3" from "FAdist"
q.p3 <- fitdist(data = caudales,distr = "lgamma3", start = my.param)

给我以下

fitdist(data = df, distr = "lgamma3", start = my.param, : 函数 mle 估计参数失败,错误代码 100

标签: rdistributiongamma-distribution

解决方案


我已经解决了问题。我使用包“PearsonDS”而不是“FAdist”来计算 log pearson III,为什么?因为后者接受偏度 (g) 的负值和连续的“比例”参数的负值(如果您查看 Pearson III 中 PearsonDS 的文档,它们会在参数中添加一个绝对值)。我必须对 Pearson III 的方程进行一些更改才能正确使用包“fitdistrplus”(以适应样本数据的分布):

最终的代码是正确的(重要提示:数据必须是原始值而不是日志(数据)),但平均值、VAR 和 G 它们需要用日志(数据)计算:

library(FAdist); library(fitdistrplus)
library("e1071")
library(PearsonDS)
load("data")
## PP for precipitation
pp <- as.numeric(data$pp)
## IMPORTANT : DATA MUST BE THE ORIGINAL VALUES NOT THE LOG(DATA)
data <- (pp)
m <- mean(log(data))
v <- var(log(data))
s <- sd(log(data))
g <- e1071::skewness(log(data), type=1)
n <- length(data)
#g <- g*(sqrt(n*(n-1))/(n-2))*(1+8.5/n)
my.shape <- (2/g)^2
my.scale <- (sqrt(v/my.shape))*sign(g) # modified as recommended by Carl Schwarz
my.thres <- m-my.shape*my.scale
#load functions
dlPIII<-function(x, shape, location, scale) {PearsonDS::dpearsonIII(log(x), shape, location, scale, log=FALSE)/x}
plPIII<-function(q, shape, location, scale) {PearsonDS::ppearsonIII(log(q), shape, location, scale, lower.tail = TRUE, log.p = FALSE)}
qlPIII<-function(p, shape, location, scale) {exp(PearsonDS::qpearsonIII(p, shape, location, scale, lower.tail = TRUE, log.p = FALSE))}
# Fit the distribution
f.pl3 <- fitdist(data, distr="lPIII", method="mle", start=list(shape=my.shape,location=my.thres ,scale=my.scale )

推荐阅读