1
在下面的代碼中,我嘗試比較模擬數據的gamma分佈與使用該數據計算出的MLE。我希望能夠完成將產生兩種密度之間的空間並將其報告爲我的錯誤的積分。查找密度與分佈之間的區域
圖表工作得很好,所以我不確定錯誤計算的原因。
非常感謝您的幫助!
library("maxLik");
#GAMMA MLE Process
GammaLogLikelihood<- function(t){
# log likelihood for Gamma(alpha,scale=beta)
alpha <- t[1]
beta <- t[2]
loglik <- sum(dgamma(x,t[1],scale=1/t[2],log=TRUE))
return(loglik)
}
GetGammaParameters<-function(x){
start<- mean(x)
GammaEst<-maxLik(GammaLogLikelihood, start=c(start,start))
return(GammaEst$estimate)
}
#Simulation
x<-rgamma(100,3,2);
mleEst<-GetGammaParameters(x);
#Graph
color2 <- rgb(1,0,0,0.2)
color1 <- rgb(0,0,1,0.2)
xax<-seq(0,max(x),.01);
plot(density(x),type="l",xlim=c(0,max(x)),ylim=c(0,1.1));
lines(xax,dgamma(xax,mleEst[1],mleEst[2]),type="l",lty=2);
polygon(density(x),density=-1,col=color1);
polygon(c(xax,max(x)),c(dgamma(xax,mleEst[1],mleEst[2]),0),density=-1,col=color2);
#Find Error
finderror<-function(data,est,l,u){
integrand<-function(x){
abs(data(x)-est(x));
}
integrate(integrand, lower = l, upper = u)
}
dataDensity<-function(x){
density(x)
}
estDensity<-function(x){
dgamma(x,mleEst[1],mleEst[2]);
}
finderror(dataDensity,estDensity,min(x),max(x));