2014-02-18 3 views
1

В следующем коде я пытаюсь сравнить гамма-распределение имитируемых данных с 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)); 

ответ

0
library(sfsmisc); 

#Find Error 
finderror<-function(densx,estDensity){ 
    newy<-abs(densx$y-estDensity(densx$x)); 
    integrate.xy(densx$x,newy); 
} 

estDensity<-function(x){ 
    dgamma(x,mleEst[1],mleEst[2]); 
} 

finderror(densx,estDensity); 

Это решило его.