#==============================================================================
# Obtendo a estimativa de maxima verossimilhanca de (X|alpha)~G( alpha; 1)
#==============================================================================
set.seed(5432)               # para gerar sempre a mesma sequencia
n       <- 20                # tamanho da amostra
alpha   <- 4                 # valor verdadeiro (hipotese)
x       <- rgamma(n,alpha,1) # amostra aleatoria
sumlogx <- sum(log(x))
sumx    <- sum(x)
xbar    <- mean(x)
sdx     <- sd(x)
sumlogx
sumx
xbar    
sdx

27.0388233398705

81.0468607244784

4.05234303622392

1.30934186435061

#==============================================================================
# Utilizando o metodo da bissecao
#==============================================================================
tol     <- 1e-10
kmax    <- 1000
k       <- 1
a       <- max(1e-10,xbar-5*sdx)
b       <- xbar+5*sdx
fa      <- -n*digamma(a) + sumlogx
fb      <- -n*digamma(b) + sumlogx
while(k <= kmax){
    c   <- (a+b)/2
    fc  <- -n*digamma(c) + sumlogx
    if(!fc || ((b-a)/2 < tol) )
        break;
    k <- k+1
    if(sign(fc)==sign(fa)){
        a  <- c
        fa <- fc
    }else{
        b  <- c
        fb <- fc
    }
}
k  # numero de iteracoes ate a convergencia
c  # estimativa de alpha

37

4.35421834117983

#==============================================================================
# Figura ilustrativa 1
#==============================================================================
xseq <- seq(1,15,length=1000) 
yseq <- -n*digamma(xseq)+sumlogx

par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(10,5,5),
    mar=c(5,5,2,2.5),xpd=T,cex.main=2.0)
plot(xseq,yseq,type="l",lwd=2,xlab=expression(alpha),xlim=c(0,16),
     ylab=expression(u(alpha)),col="red")
lines(c(0,16),c(0,0),lwd=3,col="black",lty=3)
lines(rep(c,2),c(min(yseq),max(yseq)),lwd=3,col="darkblue")

#==============================================================================
# Figura ilustrativa 2
#==============================================================================
xseq2 <- seq(0.001,16,length=10000) 
yseq2 <- -n*lgamma(xseq2)+(xseq2-1)*sumlogx - sumx

par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(10,5,5),
    mar=c(5,5,2,2.5),xpd=T,cex.main=2.0)
plot(xseq2,yseq2,type="l",lwd=2,xlab=expression(alpha),xlim=c(0,16),
     ylab=expression(paste("ln ", L(alpha),sep="")),col="red")
lines(rep(c,2),c(min(yseq2),max(yseq2)),lwd=3,col="darkblue")