In [1]:
# rm(list=ls())
options(OutDec = ",")
#==============================================================================
# Exemplo - Metodo de Newton-Raphson
# y_i ~ Poisson( exp(beta_1 + beta_2*x_i) )
# Obtendo as estimativas de maxima verossimilhanca de beta_1 e beta_2
#==============================================================================
u.poi <- function(beta,y,x){
if(length(beta)!=2 || !is.vector(beta)) return(NA)
if(length(y)!=length(x)) return(NA)
z <- rep(NA,2)
z[1] <- sum( y - exp(beta[1]+beta[2]*x) )
z[2] <- sum((y - exp(beta[1]+beta[2]*x))*x )
return(z)
}
In [2]:
#==============================================================================
# Matriz hessiana
#==============================================================================
hessiana.poi <- function(beta,y,x){
if(length(beta)!=2 || !is.vector(beta)) return(NA)
if(length(y)!=length(x)) return(NA)
H <- matrix(0,2,2)
H[1,1] <- -sum(exp(beta[1]+beta[2]*x) )
H[1,2] <- -sum(exp(beta[1]+beta[2]*x)*x )
H[2,1] <- H[1,2]
H[2,2] <- -sum(exp(beta[1]+beta[2]*x)*(x^2))
return(H)
}
In [3]:
#==============================================================================
# Newton-Raphson
#==============================================================================
newton.raphson.poi <- function(beta0=c(0,0),y,x,
nmax.iter=1000,tolerancia=1.0e-9){
if(length(beta0)!=2 || !is.vector(beta0)) return(NA)
beta <- matrix(NA,(nmax.iter+1),2)
beta[1,] <- beta0
for(i in 2:(nmax.iter+1)){
print(paste("Iteracao:",i-1,sep=""))
H <- hessiana.poi(beta[i-1,],y,x)
beta[i,] <- beta[i-1,]-solve(H)%*%u.poi(beta[i-1,],y,x)
if(sqrt(sum( (beta[i,]-beta[i-1,])^2 ) ) <tolerancia)
return(beta[i,])
}
print("Numero maximo de iteracoes atingido!")
return(beta[(nmax.iter+1),])
}
In [4]:
#==============================================================================
# Exemplo
#==============================================================================
belg.aids <- data.frame(cases=c(12,14,33,50,67,74,123,141,
165,204,253,246,240),year=1:13)
y <- belg.aids$cases
x <- belg.aids$year
In [5]:
#==============================================================================
# Valores iniciais diferentes
#==============================================================================
print("===============================")
z1 <- newton.raphson.poi(beta0=c(0,0),y=y,x=x)
print("===============================")
z2 <- newton.raphson.poi(beta0=c(10,0),y=y,x=x)
print("===============================")
print(z1)
print("===============================")
print(z2)
print("===============================")
print(u.poi(z1,y,x)) # Verificando o zero!
print("===============================")
print(u.poi(z2,y,x)) # Verificando o zero!
[1] "===============================" [1] "Iteracao:1" [1] "Iteracao:2" [1] "Iteracao:3" [1] "Iteracao:4" [1] "Iteracao:5" [1] "Iteracao:6" [1] "Iteracao:7" [1] "Iteracao:8" [1] "Iteracao:9" [1] "Iteracao:10" [1] "Iteracao:11" [1] "Iteracao:12" [1] "Iteracao:13" [1] "Iteracao:14" [1] "Iteracao:15" [1] "Iteracao:16" [1] "Iteracao:17" [1] "Iteracao:18" [1] "Iteracao:19" [1] "Iteracao:20" [1] "Iteracao:21" [1] "Iteracao:22" [1] "Iteracao:23" [1] "Iteracao:24" [1] "Iteracao:25" [1] "Iteracao:26" [1] "Iteracao:27" [1] "Iteracao:28" [1] "Iteracao:29" [1] "Iteracao:30" [1] "Iteracao:31" [1] "Iteracao:32" [1] "Iteracao:33" [1] "Iteracao:34" [1] "Iteracao:35" [1] "Iteracao:36" [1] "Iteracao:37" [1] "Iteracao:38" [1] "Iteracao:39" [1] "Iteracao:40" [1] "Iteracao:41" [1] "Iteracao:42" [1] "Iteracao:43" [1] "Iteracao:44" [1] "Iteracao:45" [1] "Iteracao:46" [1] "Iteracao:47" [1] "Iteracao:48" [1] "Iteracao:49" [1] "Iteracao:50" [1] "Iteracao:51" [1] "Iteracao:52" [1] "Iteracao:53" [1] "Iteracao:54" [1] "Iteracao:55" [1] "Iteracao:56" [1] "Iteracao:57" [1] "Iteracao:58" [1] "Iteracao:59" [1] "Iteracao:60" [1] "Iteracao:61" [1] "Iteracao:62" [1] "Iteracao:63" [1] "Iteracao:64" [1] "Iteracao:65" [1] "Iteracao:66" [1] "Iteracao:67" [1] "Iteracao:68" [1] "Iteracao:69" [1] "Iteracao:70" [1] "Iteracao:71" [1] "Iteracao:72" [1] "Iteracao:73" [1] "Iteracao:74" [1] "Iteracao:75" [1] "Iteracao:76" [1] "Iteracao:77" [1] "Iteracao:78" [1] "Iteracao:79" [1] "Iteracao:80" [1] "Iteracao:81" [1] "Iteracao:82" [1] "Iteracao:83" [1] "Iteracao:84" [1] "Iteracao:85" [1] "Iteracao:86" [1] "Iteracao:87" [1] "Iteracao:88" [1] "Iteracao:89" [1] "Iteracao:90" [1] "Iteracao:91" [1] "Iteracao:92" [1] "Iteracao:93" [1] "Iteracao:94" [1] "Iteracao:95" [1] "Iteracao:96" [1] "Iteracao:97" [1] "Iteracao:98" [1] "Iteracao:99" [1] "Iteracao:100" [1] "Iteracao:101" [1] "Iteracao:102" [1] "Iteracao:103" [1] "Iteracao:104" [1] "Iteracao:105" [1] "Iteracao:106" [1] "Iteracao:107" [1] "Iteracao:108" [1] "Iteracao:109" [1] "Iteracao:110" [1] "Iteracao:111" [1] "Iteracao:112" [1] "Iteracao:113" [1] "Iteracao:114" [1] "Iteracao:115" [1] "Iteracao:116" [1] "Iteracao:117" [1] "Iteracao:118" [1] "Iteracao:119" [1] "Iteracao:120" [1] "Iteracao:121" [1] "Iteracao:122" [1] "Iteracao:123" [1] "Iteracao:124" [1] "Iteracao:125" [1] "Iteracao:126" [1] "Iteracao:127" [1] "Iteracao:128" [1] "Iteracao:129" [1] "Iteracao:130" [1] "Iteracao:131" [1] "Iteracao:132" [1] "Iteracao:133" [1] "Iteracao:134" [1] "Iteracao:135" [1] "Iteracao:136" [1] "Iteracao:137" [1] "Iteracao:138" [1] "Iteracao:139" [1] "Iteracao:140" [1] "Iteracao:141" [1] "Iteracao:142" [1] "Iteracao:143" [1] "Iteracao:144" [1] "Iteracao:145" [1] "Iteracao:146" [1] "Iteracao:147" [1] "Iteracao:148" [1] "Iteracao:149" [1] "Iteracao:150" [1] "Iteracao:151" [1] "Iteracao:152" [1] "Iteracao:153" [1] "Iteracao:154" [1] "Iteracao:155" [1] "Iteracao:156" [1] "Iteracao:157" [1] "Iteracao:158" [1] "Iteracao:159" [1] "Iteracao:160" [1] "Iteracao:161" [1] "Iteracao:162" [1] "Iteracao:163" [1] "Iteracao:164" [1] "Iteracao:165" [1] "Iteracao:166" [1] "Iteracao:167" [1] "Iteracao:168" [1] "Iteracao:169" [1] "Iteracao:170" [1] "Iteracao:171" [1] "Iteracao:172" [1] "Iteracao:173" [1] "Iteracao:174" [1] "Iteracao:175" [1] "Iteracao:176" [1] "Iteracao:177" [1] "Iteracao:178" [1] "Iteracao:179" [1] "Iteracao:180" [1] "Iteracao:181" [1] "Iteracao:182" [1] "Iteracao:183" [1] "Iteracao:184" [1] "Iteracao:185" [1] "Iteracao:186" [1] "Iteracao:187" [1] "Iteracao:188" [1] "Iteracao:189" [1] "Iteracao:190" [1] "Iteracao:191" [1] "Iteracao:192" [1] "Iteracao:193" [1] "Iteracao:194" [1] "Iteracao:195" [1] "Iteracao:196" [1] "Iteracao:197" [1] "Iteracao:198" [1] "Iteracao:199" [1] "Iteracao:200" [1] "Iteracao:201" [1] "Iteracao:202" [1] "Iteracao:203" [1] "Iteracao:204" [1] "Iteracao:205" [1] "Iteracao:206" [1] "Iteracao:207" [1] "Iteracao:208" [1] "Iteracao:209" [1] "Iteracao:210" [1] "Iteracao:211" [1] "Iteracao:212" [1] "Iteracao:213" [1] "Iteracao:214" [1] "Iteracao:215" [1] "Iteracao:216" [1] "Iteracao:217" [1] "Iteracao:218" [1] "Iteracao:219" [1] "Iteracao:220" [1] "Iteracao:221" [1] "Iteracao:222" [1] "Iteracao:223" [1] "Iteracao:224" [1] "Iteracao:225" [1] "Iteracao:226" [1] "Iteracao:227" [1] "Iteracao:228" [1] "Iteracao:229" [1] "Iteracao:230" [1] "Iteracao:231" [1] "Iteracao:232" [1] "Iteracao:233" [1] "Iteracao:234" [1] "Iteracao:235" [1] "Iteracao:236" [1] "Iteracao:237" [1] "Iteracao:238" [1] "Iteracao:239" [1] "Iteracao:240" [1] "Iteracao:241" [1] "Iteracao:242" [1] "Iteracao:243" [1] "Iteracao:244" [1] "Iteracao:245" [1] "Iteracao:246" [1] "Iteracao:247" [1] "Iteracao:248" [1] "Iteracao:249" [1] "Iteracao:250" [1] "Iteracao:251" [1] "Iteracao:252" [1] "Iteracao:253" [1] "Iteracao:254" [1] "Iteracao:255" [1] "Iteracao:256" [1] "Iteracao:257" [1] "Iteracao:258" [1] "Iteracao:259" [1] "Iteracao:260" [1] "===============================" [1] "Iteracao:1" [1] "Iteracao:2" [1] "Iteracao:3" [1] "Iteracao:4" [1] "Iteracao:5" [1] "Iteracao:6" [1] "Iteracao:7" [1] "Iteracao:8" [1] "Iteracao:9" [1] "Iteracao:10" [1] "Iteracao:11" [1] "===============================" [1] 3,1405895 0,2021212 [1] "===============================" [1] 3,1405895 0,2021212 [1] "===============================" [1] -4,085621e-13 -5,325518e-12 [1] "===============================" [1] -5,115908e-13 -5,712764e-12
In [6]:
#==============================================================================
# Aproximando o erro padrao das estimativas
#==============================================================================
ep <- sqrt(diag(solve(-hessiana.poi(z1,y,x))))
print(ep)
[1] 0,078247000 0,007771492
In [7]:
#==============================================================================
# Estatistica de teste t e p-valor
#==============================================================================
tstat <- z1/ep
print("===============================")
print(tstat)
print("===============================")
options(digits=15)
print(2*(1-pnorm(tstat)))
[1] "===============================" [1] 40,13687 26,00803 [1] "===============================" [1] 0 0
In [8]:
#==============================================================================
# Utilizando a funcao do R: glm
#==============================================================================
ajuste <- glm(y ~ x,family=poisson("log"),start=c(10,0))
print("===============================")
print(ajuste$coef)
print("===============================")
summary(ajuste)
[1] "==============================="
(Intercept) x
3,140589536182758 0,202121203364528
[1] "==============================="
Call:
glm(formula = y ~ x, family = poisson("log"), start = c(10, 0))
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3,14058953618276 0,07824638053811 40,13719 < 2,22e-16 ***
x 0,20212120336453 0,00777144617881 26,00818 < 2,22e-16 ***
---
Signif. codes: 0 ‘***’ 0,001 ‘**’ 0,01 ‘*’ 0,05 ‘.’ 0,1 ‘ ’ 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 872,20578809760 on 12 degrees of freedom
Residual deviance: 80,68648552677 on 11 degrees of freedom
AIC: 166,3698159985
Number of Fisher Scoring iterations: 10
In [9]:
# rm(list=ls())
#==============================================================================
# Fim
#==============================================================================