R软件拟合统计分布并绘制曲线

在统计分析中， 经常会用到统计分布拟合， 这里介绍MASS程序包的 fitdistr(), 该函数可以用来拟合  "beta", "cauchy", "chi-squared", "exponential", "f", "gamma", "geometric", "log-normal", "lognormal", "logistic", "negative binomial", "normal", "Poisson", "t" 以及 "weibull" 的概率密度分布， 现举例如下：

library(MASS)

#### Gammadistribution

x <-rgamma(100, shape =5, rate =0.1)

res <- fitdistr(x, "gamma")

h <-hist(x)

xfit <-seq(min(x), max(x), by=(max(x)-min(x))/100)

yfit <-dgamma(xfit,res[[1]][1],res[[1]][2])

yfit <- yfit*diff(h$mids[1:2])*length(xfit)

lines(xfit, yfit, col="blue", lwd=2)

####"beta"                

x <-rbeta(100, shape1 =5, shape2 =2)

res <- fitdistr(x, "beta", start=list(shape1 =4, shape2 =1))

h <-hist(x)

xfit <-seq(min(x), max(x), by=(max(x)-min(x))/100)

yfit <-dbeta(xfit,res[[1]][1],res[[1]][2])

yfit <- yfit*diff(h$mids[1:2])*length(xfit)

lines(xfit, yfit, col="blue", lwd=2)

####"cauchy",            

"cauchy",

x <-rcauchy(100, location =200, scale=300)

res <- fitdistr(x, "cauchy")

h <-hist(x)

xfit <-seq(min(x), max(x), by=(max(x)-min(x))/100)

yfit <-dcauchy(xfit,res[[1]][1],res[[1]][2])

yfit <- yfit*diff(h$mids[1:2])*length(xfit)

lines(xfit, yfit, col="blue", lwd=2)

####"chi-squared",        

"chi-squared",

x <-rchisq(100, df=5, ncp =2)

res <- fitdistr(x, "chi-squared", start=list(df=2, ncp =1))

h <-hist(x)

xfit <-seq(min(x), max(x), by=(max(x)-min(x))/100)

yfit <-dchisq(xfit,res[[1]][1],res[[1]][2])

yfit <- yfit*diff(h$mids[1:2])*length(xfit)

lines(xfit, yfit, col="blue", lwd=2)

####"exponential",        

"exponential"

x <-rexp(100, rate =10)

res <- fitdistr(x, "exponential")

h <-hist(x)

xfit <-seq(min(x), max(x), by=(max(x)-min(x))/100)

yfit <-dexp(xfit,res[[1]][1])

yfit <- yfit*diff(h$mids[1:2])*length(xfit)

lines(xfit, yfit, col="blue", lwd=2)

####"f",                  

"f"

x <-rf(100, df1 =60, df2 =50)

res <- fitdistr(x, "f", start=list(df1 =10, df2 =2))

h <-hist(x)

xfit <-seq(min(x), max(x), by=(max(x)-min(x))/100)

yfit <-df(xfit,res[[1]][1], res[[1]][2])

yfit <- yfit*diff(h$mids[1:2])*length(xfit)

lines(xfit, yfit, col="blue", lwd=2)

####"lognormal",          

"lognormal"

x <-rlnorm(100, meanlog =5, sdlog =1.5)

res <- fitdistr(x, "lognormal")

h <-hist(x)

xfit <-seq(min(x), max(x), by=(max(x)-min(x))/100)

yfit <-dlnorm(xfit,res[[1]][1], res[[1]][2])

yfit <- yfit*diff(h$mids[1:2])*length(xfit)

lines(xfit, yfit, col="blue", lwd=2)

####"logistic",          

"logistic",

x <-rlogis(100, location =5, scale=2)

res <- fitdistr(x, "logistic", start=list(location=5, scale=2))

h <-hist(x)

xfit <-seq(min(x), max(x), by=(max(x)-min(x))/100)

yfit <-dlogis(xfit,res[[1]][1], res[[1]][2])

yfit <- yfit*diff(h$mids[1:2])*length(xfit)

lines(xfit, yfit, col="blue", lwd=2)

#### "negativebinomial",  

"negative binomial",

x <-rnbinom(100, size =300, prob =.3)

res <- fitdistr(x, "negative binomial")

h <-hist(x)

xfit <-floor(seq(min(x), max(x), by=(max(x)-min(x))/100))

yfit <-dnbinom(xfit,size = res[[1]][1], mu = res[[1]][2])

yfit <- yfit*diff(h$mids[1:2])*length(xfit)

lines(xfit, yfit, col="blue", lwd=2)

####"normal",            

"normal"

x <-rnorm(100, mean=15, sd=2)

res <- fitdistr(x, "normal")

h <-hist(x)

xfit <-seq(min(x), max(x), by=(max(x)-min(x))/100)

yfit <-dnorm(xfit, mean= res[[1]][1], sd= res[[1]][2])

yfit <- yfit*diff(h$mids[1:2])*length(xfit)

lines(xfit, yfit, col="blue", lwd=2)

####"Poisson",            

"Poisson",

x <-rpois(100, lambda =400)

res <- fitdistr(x, "Poisson")

h <-hist(x)

xfit <-floor(seq(min(x), max(x), by=(max(x)-min(x))/100))

yfit <-dpois(xfit,res[[1]][1])

yfit <- yfit*diff(h$mids[1:2])*length(xfit)

lines(xfit, yfit, col="blue", lwd=2)

####"t"                  

"t"

x <-rt(100, df=5)

res <- fitdistr(x, "t")

h <-hist(x)

xfit <-seq(min(x), max(x), by=(max(x)-min(x))/100)

yfit <-dt(xfit,res[[1]][3])  

yfit <- yfit*diff(h$mids[1:2])*length(xfit)

lines(xfit, yfit, col="blue", lwd=2)

####"weibull"            

"weibull"

x <-rweibull(100, shape =5, scale=9)

res <- fitdistr(x, "weibull", start=list(shape =1, scale=1))

h <-hist(x)

xfit <-seq(min(x), max(x), by=(max(x)-min(x))/100)

yfit <-dweibull(xfit,res[[1]][1], res[[1]][2])

yfit <- yfit*diff(h$mids[1:2])*length(xfit)

lines(xfit, yfit, col="blue", lwd=2)