R语言如何用gamlss从最佳拟合分布中生成随机数？

6ojccjat 于 2023-03-15 发布在其他

关注(0)|答案(1)|浏览(351)

目标

查找数据的最佳拟合分布，然后从此分布生成随机数。

示例

使用R中的gamlss包，我发现最佳拟合分布是"Skew exponential power (Azzalini type 1)"：

library(gamlss)
library(gamlss.dist)
library(gamlss.add)

m1 <- fitDist(mtcars$wt, k = 2, type = "realAll", trace = FALSE, try.gamlss = TRUE)

  
summary(m1)
# *******************************************************************
#   Family:  c("SEP1", "Skew exponential power (Azzalini type 1)") 
# 
# Call:  gamlssML(formula = y, family = DIST[i]) 
# 
# Fitting method: "nlminb" 
# 
# 
# Coefficient(s):
#   Estimate   Std. Error     t value Pr(>|t|)
# eta.mu     3.440000000  0.000149516 23007.56651  < 2e-16
# eta.sigma -1.856665040  0.861826733    -2.15434 0.031214
# eta.nu     0.150728244           NA          NA       NA
# eta.tau   -3.524272086           NA          NA       NA
# 
# eta.mu    ***
#   eta.sigma *  
#   eta.nu       
# eta.tau      
# ---
#   Signif. codes:  
#   0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Degrees of Freedom for the fit: 4 Residual Deg. of Freedom   28 
# Global Deviance:     -27.4747 
# AIC:     -19.4747 
# SBC:     -13.6117 
# Warning message:
#   In sqrt(diag(object$vcov)) : NaNs produced

错误

但是sigma和tau的值是负数，当我把这些值提供给rSEP1()生成一个随机数时，它抛出了下面的错误：

rSEP1(1, mu = 3.44, sigma = -1.8, nu = 0.15, tau = -3.5)
# Error in rSEP1(1, mu = 3.44, sigma = -1.8, nu = 0.15, tau = -3.5) : 
#   sigma must be positive

这些值是否已转换？如何为rSEP1()提供正确的输入？

r

来源：https://stackoverflow.com/questions/75738376/how-to-generate-random-numbers-from-the-best-fitted-distribution-with-gamlss

1条答案

按热度按时间

kgsdhlau1#

如果查看参数的链接函数，您将看到：

SEP1()

#> GAMLSS Family: SEP1 Skew exponential power (Azzalini type 1) 
#> Link function for mu   : identity 
#> Link function for sigma: log 
#> Link function for nu   : identity 
#> Link function for tau  : log

因此，fitDist返回的sigma和tau值是您要插入到rSEP1中的数字的 log。换句话说，从模型运行rSEP1的正确方法是：

rSEP1(1, mu = 3.44, sigma = exp(-1.8), nu = 0.15, tau = exp(-3.5))
#> [1] 3.460991

为了证明这一点，我们从定义的rSEP1创建一组随机数，其中sigma = 0.5，tau = 0.5。如果我们随后找到最佳拟合分布，我们应该得到sigma和tau都接近log（0.5）的结果。由于log(0.5)约为-0.69，我们预计这两个参数的值约为-0.69：

set.seed(1)
testvec <- rSEP1(10000, mu = 3.5, sigma = 0.5, nu = 2.5, tau = 0.5)
m2 <- fitDist(testvec, k = 2, type = "realAll", trace = FALSE, try.gamlss = TRUE)

m2
#> Family:  c("SEP1", "Skew exponential power (Azzalini type 1)") 
#> Fitting method: "nlminb" 
#> 
#> Call:  gamlssML(formula = y, family = DIST[i]) 
#> 
#> Mu Coefficients:
#> [1]  3.5
#> Sigma Coefficients:
#> [1]  -0.6828
#> Nu Coefficients:
#> [1]  2.421
#> Tau Coefficients:
#> [1]  -0.6952

将sigma和tau的指数代入dSEP1，我们可以得到测试向量密度的近乎完美的拟合：

d <- density(testvec)
plot(d$x, d$y)
lines(d$x, dSEP1(d$x, mu = 3.5, sigma = exp(-0.68), nu = 2.42, tau = exp(-0.69)),
     col = "red", lwd = 2, lty = 2)

值得指出的是，示例中mtcars$wt的实际拟合非常糟糕，sigma jyst的值很小，这意味着从分布中提取的大多数随机值将接近mtcars$wt的均值，数据集中只有32个点，这使得自动精确拟合参数分布非常困难。

赞(0）回复(0）举报 2023-03-15

我来回答

R语言如何用gamlss从最佳拟合分布中生成随机数？

目标

示例

错误

1条答案

相关问题

热门标签

最新问答

R语言 如何用gamlss从最佳拟合分布中生成随机数？

目标

示例

错误

1条答案

相关问题

热门标签

最新问答

R语言如何用gamlss从最佳拟合分布中生成随机数？