将边缘分布拟合到直方图中的示例R

Question

有人可以告诉我如何将多项式边际分布拟合到我的数据中吗？我做了二项和二项二项式，但我想看看如何拟合一个多项式。如果这是你知道该怎么做的，我也会对尝试伽玛感兴趣。

这是我迄今为止所做的。

nodes <- read.table("https://web.stanford.edu/~hastie/CASI_files/DATA/nodes.txt", 
      header = T) 

nodes %>% 
ggplot(aes(x=x/n))+ 
    geom_histogram(bins = 30)+ 
    theme_bw()+ 
    labs(x = "nodes", 
     n = "p=x/n") 

# log-likelihood function 
ll <- function(alpha, beta) { 
x <- nodes$x 
total <- nodes$n 
-sum(VGAM::dbetabinom.ab(x, total, alpha, beta, log = TRUE)) 
} 

# maximum likelihood estimation 
m <- mle(ll, start = list(alpha = 1, beta = 10), method = "L-BFGS-B", 
lower = c(0.0001, .1)) 
ab <- coef(m) 
alpha0 <- ab[1] 
beta0 <- ab[2] 

nodes %>% 
    ggplot() + 
    geom_histogram(aes(x/n, y = ..density..), bins= 30) + 
    stat_function(fun = function(x) dbeta(x, alpha0, beta0), color = "red", 
       size = 1) + 
    xlab("p=x/n") 

这里是另一个适合

ll <- function(a){ 
    x <- nodes$x 
    total <- nodes$n 
    -sum(stats::dbinom(x, total, a, log = TRUE)) 
} 

#stats::dbinom() 
m <- mle(ll, start = list(a=.5), method = "L-BFGS-B", 
lower = c(0.0001, .1)) 

a = coef(m) 

nodes %>% 
    ggplot() + 
    geom_histogram(aes(x/n, y = ..density..), bins=40) + 
    stat_function(fun = function(x) dbeta(x, a, 1), color = "red", 
       size = 1) + 
    xlab("proportion x/n") 

Answer 1

拟合gamma分布：

data(iris) 
library(MASS) ##for the fitdistr function 

fit.params <- fitdistr(iris$Sepal.Length, "gamma", lower = c(0, 0)) 

ggplot(data = iris) + 
geom_histogram(data = as.data.frame(x), aes(x=iris$Sepal.Length, y=..density..)) + 
geom_line(aes(x=iris$Sepal.Length, 
y=dgamma(iris$Sepal.Length,fit.params$estimate["shape"], 
fit.params$estimate["rate"])), color="red", size = 1) + 
theme_classic() 

你可能阿尔斯o喜欢使用汽车包装中的qqp函数来查看分位数的分布。以下是几个示例：

library(car) 
qqp(iris$Sepal.Length, "norm") ##normal distribution 

qqp(iris$Sepal.Length, "lnorm") ##log-normal distribution 

gamma <- fitdistr(iris$Sepal.Length, "gamma") 
qqp(iris$Sepal.Length, "gamma", shape = gamma$estimate[[1]], 
rate = gamma$estimate[[2]]) ##gamma distribution 

nbinom <- fitdistr(iris$Sepal.Length, "Negative Binomial") 
qqp(iris$Sepal.Length, "nbinom", size = nbinom$estimate[[1]], 
mu = nbinom$estimate[[2]]) ##negative binomial distribution 

您可以对ggplots或qqPlot使用fitdistr函数。它支持很多不同的发行版。看看？fitdistr