# rm(list=ls()) 
options(OutDec = ",")
#==============================================================================
# Exemplo:
#     Regressao linear intrinseca.
#==============================================================================
# Analise: Dados artificiais do modelo gama
#==============================================================================
library(MASS)
dataxy      <- read.table("../dados/d10_reg_lin_intrinseca.txt",header=T)
y           <- dataxy[,2]
x           <- dataxy[,3]
n           <- length(y)

#==============================================================================
# Ajuste da Regressaom Linear Intrinseca por Minimos Qadrados
#==============================================================================
fitymq  <- lm(y ~ x)
rhomq   <- as.vector(fitymq$coef[2])
betamq  <- as.vector(fitymq$coef[1]/fitymq$coef[2])  # b_1/ b_2
bcov    <- matrix(vcov(fitymq),2,2,T)
sdrho   <- sqrt(bcov[2,2]) # erro padrao do rho  (EMQ)
# Aplicacao do metodo delta para estimar a variancia de beta

# Erro padrao do beta (EMQ)
sdbeta  <- as.vector(sqrt(bcov[1,1]/(fitymq$coef[2]^2) 
         + (fitymq$coef[1]^2)*bcov[2,2]/(fitymq$coef[2]^4)
         - 2*fitymq$coef[1]*bcov[1,2]/(fitymq$coef[2]^3))) 
print(c(rhomq,sdrho))
print(c(betamq,sdbeta))

[1] 2,426104 1,591482
[1] -1,707724  8,689001

#==============================================================================
# Estimacao por Maxima Verossimilhanca
#==============================================================================
# Analisando para rho: (beta fixo). Podemos pular esta parte.
#==============================================================================
myfrho  <- function(rho,x,y,beta){
    n   <- length(y)
    out <- -(-rho*sum(log(beta+x))-n*lgamma(rho)
           +(rho-1)*sum(log(y))-sum(y/(beta+x)))
    return(out)
}
fitrho <- optimize(myfrho,lower=0,upper=100,x=x,y=y,beta=-4.7198)
print(fitrho)

$minimum
[1] 3,151288

$objective
[1] 82,91605

#==============================================================================
# Analisando para beta: (rho fixo). Podemos pular esta parte.
#==============================================================================
myfbeta <- function(beta,x,y,rho){
    n   <- length(y)
    out <- -(-rho*sum(log(beta+x))-n*lgamma(rho)
           +(rho-1)*sum(log(y))-sum(y/(beta+x)))
    return(out)
}
fitbeta <- optimize(myfbeta,lower=max(-x),upper=100,x=x,y=y,rho=3.151)
print(fitbeta)

$minimum
[1] -4,718759

$objective
[1] 82,91605

#==============================================================================
# Estimacao por Maxima Verossimilhanca: rho e beta 
# Estamos na verdade minimizando o negativo da verossimilhanca
#==============================================================================
myf <- function(w,x,y){
    rhox  <- w[1]
    betax <- w[2]
    n     <- length(y)
    out   <- -(-rhox*sum(log(betax+x))-n*lgamma(rhox)
             +(rhox-1)*sum(log(y))-sum(y/(betax+x)))
    return(out)
}
fityx    <- optim(c(rhomq,betamq),myf,x=x,y=y,lower=c(0,max(-x)), 
                     method="L-BFGS-B",hessian=T)
print(fityx)
H        <- matrix(NA,2,2)
H[1,1]   <- -n*trigamma(fityx$par[1])
H[2,2]   <- ( fityx$par[1]*sum(1/((fityx$par[2]+x)^2))
            -2*sum(y/((fityx$par[2]+x)^3)))
H[1,2]   <- -sum(1/(fityx$par[2]+x))
H[2,1]   <- H[1,2]
V        <- solve(-H)
G        <- solve(fityx$hessian)
sdrhomv  <- sqrt(V[1,1]) # erro padrao do rho  (EMV)
sdbetamv <- sqrt(V[2,2]) # erro padrao do beta (EMV)
rhomv    <- fityx$par[1]
betamv   <- fityx$par[2]
print(c(rhomv,sdrhomv,sqrt(G[1,1])))
print(c(betamv,sdbetamv,sqrt(G[2,2])))

$par
[1]  3,150896 -4,718504

$value
[1] 82,91605

$counts
function gradient 
      13       13 

$convergence
[1] 0

$message
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

$hessian
         [,1]     [,2]
[1,] 7,459183 2,241969
[2,] 2,241969 0,855704

[1] 3,1508965 0,7942620 0,7942616
[1] -4,718504  2,345026  2,345025