ggplot2 - Example of fitting marginal distributions to histogram in R -
could show me how fit polynomial marginal distribution data? have done binomial , beta binomial, see how fit polynomial. interested in trying gamma if know how do.
this have done far.
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") 
here fit
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)) = 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") 
for fitting gamma distribution:
data(iris) library(mass) ##for 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() you might take @ distribution of quantiles using qqp function in car package. here few examples:
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 you can use fitdistr function ggplots or qqplots. supports lots of different distributions. take @ ?fitdistr
Comments
Post a Comment