In [1]:
# rm(list=ls())
options(OutDec = ",")
#==============================================================================
# Exemplo da aplicacao do algoritmo Esperanca-Maximizacao
# Modelo Poisson com observação faltante
#==============================================================================
# Gerando dados artificiais
#==============================================================================
n <- 100 # Tamanho da amostra
set.seed(2345) # Semente
beta <- 5
tau <- runif(n,0.5,4)
x <- rpois(n,tau)
y <- rpois(n,beta*tau)
In [2]:
#==============================================================================
# As estimativas se nao tivessemos dados faltantes
#==============================================================================
beta.hat <- mean(y)/mean(x)
tau.hat <- (x+y)/(beta.hat+1)
print(c(beta,beta.hat))
print(cbind(tau,tau.hat))
[1] 5,000000 4,491736
tau tau.hat
[1,] 0,9086021 0,9104590
[2,] 1,1825879 1,2746426
[3,] 2,9801014 2,5492852
[4,] 0,6205911 0,5462754
[5,] 2,1629203 2,7313770
[6,] 1,5295645 1,8209180
[7,] 2,6163286 2,3671934
[8,] 3,2596929 3,2776524
[9,] 1,9254120 2,9134688
[10,] 2,9917959 5,0985704
[11,] 1,0576072 1,6388262
[12,] 1,6990550 2,5492852
[13,] 0,7837161 0,3641836
[14,] 1,0247577 0,7283672
[15,] 1,9367184 1,8209180
[16,] 1,6563481 2,1851016
[17,] 2,6490742 2,7313770
[18,] 2,0197395 2,5492852
[19,] 2,8632761 3,4597442
[20,] 1,1607834 1,8209180
[21,] 3,9302707 5,0985704
[22,] 2,2313993 2,7313770
[23,] 3,3652788 2,9134688
[24,] 2,2621914 3,2776524
[25,] 1,7267421 2,5492852
[26,] 0,9907801 1,4567344
[27,] 3,6244448 3,2776524
[28,] 2,6322905 3,4597442
[29,] 2,5238203 2,3671934
[30,] 3,5101145 3,2776524
[31,] 1,0887114 0,5462754
[32,] 2,3059721 3,2776524
[33,] 3,8800979 4,3702032
[34,] 2,5782177 3,0955606
[35,] 0,7237824 0,7283672
[36,] 3,4220357 3,6418360
[37,] 2,4823936 2,7313770
[38,] 3,0178928 2,5492852
[39,] 0,8311088 0,7283672
[40,] 2,4949370 1,8209180
[41,] 3,1948007 3,2776524
[42,] 2,3792959 1,8209180
[43,] 2,3223436 3,0955606
[44,] 3,2925017 2,9134688
[45,] 1,5809458 0,9104590
[46,] 3,9443465 3,0955606
[47,] 3,2782473 3,8239278
[48,] 2,3599059 3,0955606
[49,] 1,5272572 1,4567344
[50,] 1,1610151 1,6388262
[51,] 2,2482110 2,1851016
[52,] 3,3489085 3,4597442
[53,] 2,6387077 4,0060196
[54,] 1,2127710 0,9104590
[55,] 2,5998391 3,6418360
[56,] 2,0714851 2,1851016
[57,] 1,7103704 2,7313770
[58,] 3,1354192 3,0955606
[59,] 2,6269831 2,5492852
[60,] 0,5226903 1,2746426
[61,] 1,7032974 1,2746426
[62,] 2,8188487 2,7313770
[63,] 2,6219105 2,1851016
[64,] 1,4009743 2,1851016
[65,] 0,7963026 1,0925508
[66,] 3,3604122 2,9134688
[67,] 1,0657819 1,2746426
[68,] 1,0296527 0,9104590
[69,] 1,7779494 2,3671934
[70,] 3,7249433 4,1881114
[71,] 3,2659476 2,9134688
[72,] 0,5563507 0,9104590
[73,] 2,6264871 2,5492852
[74,] 0,6416980 0,5462754
[75,] 3,6084877 3,2776524
[76,] 2,2958531 2,7313770
[77,] 2,4894541 3,0955606
[78,] 1,6843026 2,7313770
[79,] 3,6305934 4,1881114
[80,] 0,8824017 1,4567344
[81,] 3,5453178 5,2806622
[82,] 1,2327081 1,2746426
[83,] 2,5594922 2,9134688
[84,] 1,7871659 1,6388262
[85,] 1,3167711 0,5462754
[86,] 2,1498119 2,7313770
[87,] 1,8682297 2,0030098
[88,] 1,4086294 1,4567344
[89,] 1,7585881 2,1851016
[90,] 3,5880890 5,0985704
[91,] 2,5021175 3,6418360
[92,] 2,9723305 2,3671934
[93,] 0,9036160 1,8209180
[94,] 3,9443753 5,8269375
[95,] 2,4859370 1,8209180
[96,] 0,6793552 0,3641836
[97,] 1,8643601 1,8209180
[98,] 0,5715798 0,3641836
[99,] 1,0222745 1,4567344
[100,] 1,7606231 3,4597442
In [3]:
#==============================================================================
# Com x1 faltante; algoritmo EM
#==============================================================================
# Passo E ja foi calculado
# Passo M
M <- 10 # Numero de iteracoes do algoritmo
betaM <- rep(NA,M)
tauM <- matrix(NA,M,n)
betaM[1] <- 1 # Valor inicial
tauM[1,] <- rep(1,n)
sumy <- sum(y)
sumx <- sum(x[2:n])
xplusy <- x[2:n]+y[2:n]
for(i in 2:M){
betaM[i] <- sumy/(sumx+tauM[i-1,1])
tauM[i,1] <- (tauM[i-1,1]+y[1])/(betaM[i]+1)
tauM[i,2:n] <- xplusy/(betaM[i]+1)
}
print(c(beta,betaM))
print(rbind(tau,tauM))
[1] 5,000000 1,000000 4,491736 4,493398 4,493706 4,493763 4,493774 4,493775
[9] 4,493776 4,493776 4,493776
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
tau 0,9086021 1,182588 2,980101 0,6205911 2,162920 1,529565 2,616329 3,259693
1,0000000 1,000000 1,000000 1,0000000 1,000000 1,000000 1,000000 1,000000
0,9104590 1,274643 2,549285 0,5462754 2,731377 1,820918 2,367193 3,277652
0,8938837 1,274257 2,548514 0,5461101 2,730550 1,820367 2,366477 3,276660
0,8908165 1,274185 2,548371 0,5460795 2,730397 1,820265 2,366344 3,276477
0,8902489 1,274172 2,548344 0,5460738 2,730369 1,820246 2,366320 3,276443
0,8901439 1,274170 2,548339 0,5460727 2,730364 1,820242 2,366315 3,276436
0,8901244 1,274169 2,548339 0,5460726 2,730363 1,820242 2,366314 3,276435
0,8901209 1,274169 2,548338 0,5460725 2,730363 1,820242 2,366314 3,276435
0,8901202 1,274169 2,548338 0,5460725 2,730363 1,820242 2,366314 3,276435
0,8901201 1,274169 2,548338 0,5460725 2,730363 1,820242 2,366314 3,276435
[,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16]
tau 1,925412 2,991796 1,057607 1,699055 0,7837161 1,0247577 1,936718 1,656348
1,000000 1,000000 1,000000 1,000000 1,0000000 1,0000000 1,000000 1,000000
2,913469 5,098570 1,638826 2,549285 0,3641836 0,7283672 1,820918 2,185102
2,912587 5,097027 1,638330 2,548514 0,3640734 0,7281468 1,820367 2,184440
2,912424 5,096742 1,638238 2,548371 0,3640530 0,7281059 1,820265 2,184318
2,912394 5,096689 1,638221 2,548344 0,3640492 0,7280984 1,820246 2,184295
2,912388 5,096679 1,638218 2,548339 0,3640485 0,7280970 1,820242 2,184291
2,912387 5,096677 1,638218 2,548339 0,3640484 0,7280967 1,820242 2,184290
2,912387 5,096677 1,638218 2,548338 0,3640483 0,7280967 1,820242 2,184290
2,912387 5,096677 1,638218 2,548338 0,3640483 0,7280967 1,820242 2,184290
2,912387 5,096677 1,638218 2,548338 0,3640483 0,7280967 1,820242 2,184290
[,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
tau 2,649074 2,019740 2,863276 1,160783 3,930271 2,231399 3,365279 2,262191
1,000000 1,000000 1,000000 1,000000 1,000000 1,000000 1,000000 1,000000
2,731377 2,549285 3,459744 1,820918 5,098570 2,731377 2,913469 3,277652
2,730550 2,548514 3,458697 1,820367 5,097027 2,730550 2,912587 3,276660
2,730397 2,548371 3,458503 1,820265 5,096742 2,730397 2,912424 3,276477
2,730369 2,548344 3,458467 1,820246 5,096689 2,730369 2,912394 3,276443
2,730364 2,548339 3,458461 1,820242 5,096679 2,730364 2,912388 3,276436
2,730363 2,548339 3,458459 1,820242 5,096677 2,730363 2,912387 3,276435
2,730363 2,548338 3,458459 1,820242 5,096677 2,730363 2,912387 3,276435
2,730363 2,548338 3,458459 1,820242 5,096677 2,730363 2,912387 3,276435
2,730363 2,548338 3,458459 1,820242 5,096677 2,730363 2,912387 3,276435
[,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32]
tau 1,726742 0,9907801 3,624445 2,632290 2,523820 3,510115 1,0887114 2,305972
1,000000 1,0000000 1,000000 1,000000 1,000000 1,000000 1,0000000 1,000000
2,549285 1,4567344 3,277652 3,459744 2,367193 3,277652 0,5462754 3,277652
2,548514 1,4562935 3,276660 3,458697 2,366477 3,276660 0,5461101 3,276660
2,548371 1,4562119 3,276477 3,458503 2,366344 3,276477 0,5460795 3,276477
2,548344 1,4561968 3,276443 3,458467 2,366320 3,276443 0,5460738 3,276443
2,548339 1,4561940 3,276436 3,458461 2,366315 3,276436 0,5460727 3,276436
2,548339 1,4561935 3,276435 3,458459 2,366314 3,276435 0,5460726 3,276435
2,548338 1,4561934 3,276435 3,458459 2,366314 3,276435 0,5460725 3,276435
2,548338 1,4561934 3,276435 3,458459 2,366314 3,276435 0,5460725 3,276435
2,548338 1,4561934 3,276435 3,458459 2,366314 3,276435 0,5460725 3,276435
[,33] [,34] [,35] [,36] [,37] [,38] [,39] [,40]
tau 3,880098 2,578218 0,7237824 3,422036 2,482394 3,017893 0,8311088 2,494937
1,000000 1,000000 1,0000000 1,000000 1,000000 1,000000 1,0000000 1,000000
4,370203 3,095561 0,7283672 3,641836 2,731377 2,549285 0,7283672 1,820918
4,368881 3,094624 0,7281468 3,640734 2,730550 2,548514 0,7281468 1,820367
4,368636 3,094450 0,7281059 3,640530 2,730397 2,548371 0,7281059 1,820265
4,368590 3,094418 0,7280984 3,640492 2,730369 2,548344 0,7280984 1,820246
4,368582 3,094412 0,7280970 3,640485 2,730364 2,548339 0,7280970 1,820242
4,368580 3,094411 0,7280967 3,640484 2,730363 2,548339 0,7280967 1,820242
4,368580 3,094411 0,7280967 3,640483 2,730363 2,548338 0,7280967 1,820242
4,368580 3,094411 0,7280967 3,640483 2,730363 2,548338 0,7280967 1,820242
4,368580 3,094411 0,7280967 3,640483 2,730363 2,548338 0,7280967 1,820242
[,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48]
tau 3,194801 2,379296 2,322344 3,292502 1,5809458 3,944346 3,278247 2,359906
1,000000 1,000000 1,000000 1,000000 1,0000000 1,000000 1,000000 1,000000
3,277652 1,820918 3,095561 2,913469 0,9104590 3,095561 3,823928 3,095561
3,276660 1,820367 3,094624 2,912587 0,9101834 3,094624 3,822770 3,094624
3,276477 1,820265 3,094450 2,912424 0,9101324 3,094450 3,822556 3,094450
3,276443 1,820246 3,094418 2,912394 0,9101230 3,094418 3,822517 3,094418
3,276436 1,820242 3,094412 2,912388 0,9101212 3,094412 3,822509 3,094412
3,276435 1,820242 3,094411 2,912387 0,9101209 3,094411 3,822508 3,094411
3,276435 1,820242 3,094411 2,912387 0,9101209 3,094411 3,822508 3,094411
3,276435 1,820242 3,094411 2,912387 0,9101208 3,094411 3,822508 3,094411
3,276435 1,820242 3,094411 2,912387 0,9101208 3,094411 3,822508 3,094411
[,49] [,50] [,51] [,52] [,53] [,54] [,55] [,56]
tau 1,527257 1,161015 2,248211 3,348909 2,638708 1,2127710 2,599839 2,071485
1,000000 1,000000 1,000000 1,000000 1,000000 1,0000000 1,000000 1,000000
1,456734 1,638826 2,185102 3,459744 4,006020 0,9104590 3,641836 2,185102
1,456294 1,638330 2,184440 3,458697 4,004807 0,9101834 3,640734 2,184440
1,456212 1,638238 2,184318 3,458503 4,004583 0,9101324 3,640530 2,184318
1,456197 1,638221 2,184295 3,458467 4,004541 0,9101230 3,640492 2,184295
1,456194 1,638218 2,184291 3,458461 4,004533 0,9101212 3,640485 2,184291
1,456193 1,638218 2,184290 3,458459 4,004532 0,9101209 3,640484 2,184290
1,456193 1,638218 2,184290 3,458459 4,004532 0,9101209 3,640483 2,184290
1,456193 1,638218 2,184290 3,458459 4,004532 0,9101208 3,640483 2,184290
1,456193 1,638218 2,184290 3,458459 4,004532 0,9101208 3,640483 2,184290
[,57] [,58] [,59] [,60] [,61] [,62] [,63] [,64]
tau 1,710370 3,135419 2,626983 0,5226903 1,703297 2,818849 2,621911 1,400974
1,000000 1,000000 1,000000 1,0000000 1,000000 1,000000 1,000000 1,000000
2,731377 3,095561 2,549285 1,2746426 1,274643 2,731377 2,185102 2,185102
2,730550 3,094624 2,548514 1,2742568 1,274257 2,730550 2,184440 2,184440
2,730397 3,094450 2,548371 1,2741854 1,274185 2,730397 2,184318 2,184318
2,730369 3,094418 2,548344 1,2741722 1,274172 2,730369 2,184295 2,184295
2,730364 3,094412 2,548339 1,2741697 1,274170 2,730364 2,184291 2,184291
2,730363 3,094411 2,548339 1,2741693 1,274169 2,730363 2,184290 2,184290
2,730363 3,094411 2,548338 1,2741692 1,274169 2,730363 2,184290 2,184290
2,730363 3,094411 2,548338 1,2741692 1,274169 2,730363 2,184290 2,184290
2,730363 3,094411 2,548338 1,2741692 1,274169 2,730363 2,184290 2,184290
[,65] [,66] [,67] [,68] [,69] [,70] [,71] [,72]
tau 0,7963026 3,360412 1,065782 1,0296527 1,777949 3,724943 3,265948 0,5563507
1,0000000 1,000000 1,000000 1,0000000 1,000000 1,000000 1,000000 1,0000000
1,0925508 2,913469 1,274643 0,9104590 2,367193 4,188111 2,913469 0,9104590
1,0922201 2,912587 1,274257 0,9101834 2,366477 4,186844 2,912587 0,9101834
1,0921589 2,912424 1,274185 0,9101324 2,366344 4,186609 2,912424 0,9101324
1,0921476 2,912394 1,274172 0,9101230 2,366320 4,186566 2,912394 0,9101230
1,0921455 2,912388 1,274170 0,9101212 2,366315 4,186558 2,912388 0,9101212
1,0921451 2,912387 1,274169 0,9101209 2,366314 4,186556 2,912387 0,9101209
1,0921450 2,912387 1,274169 0,9101209 2,366314 4,186556 2,912387 0,9101209
1,0921450 2,912387 1,274169 0,9101208 2,366314 4,186556 2,912387 0,9101208
1,0921450 2,912387 1,274169 0,9101208 2,366314 4,186556 2,912387 0,9101208
[,73] [,74] [,75] [,76] [,77] [,78] [,79] [,80]
tau 2,626487 0,6416980 3,608488 2,295853 2,489454 1,684303 3,630593 0,8824017
1,000000 1,0000000 1,000000 1,000000 1,000000 1,000000 1,000000 1,0000000
2,549285 0,5462754 3,277652 2,731377 3,095561 2,731377 4,188111 1,4567344
2,548514 0,5461101 3,276660 2,730550 3,094624 2,730550 4,186844 1,4562935
2,548371 0,5460795 3,276477 2,730397 3,094450 2,730397 4,186609 1,4562119
2,548344 0,5460738 3,276443 2,730369 3,094418 2,730369 4,186566 1,4561968
2,548339 0,5460727 3,276436 2,730364 3,094412 2,730364 4,186558 1,4561940
2,548339 0,5460726 3,276435 2,730363 3,094411 2,730363 4,186556 1,4561935
2,548338 0,5460725 3,276435 2,730363 3,094411 2,730363 4,186556 1,4561934
2,548338 0,5460725 3,276435 2,730363 3,094411 2,730363 4,186556 1,4561934
2,548338 0,5460725 3,276435 2,730363 3,094411 2,730363 4,186556 1,4561934
[,81] [,82] [,83] [,84] [,85] [,86] [,87] [,88]
tau 3,545318 1,232708 2,559492 1,787166 1,3167711 2,149812 1,868230 1,408629
1,000000 1,000000 1,000000 1,000000 1,0000000 1,000000 1,000000 1,000000
5,280662 1,274643 2,913469 1,638826 0,5462754 2,731377 2,003010 1,456734
5,279064 1,274257 2,912587 1,638330 0,5461101 2,730550 2,002404 1,456294
5,278768 1,274185 2,912424 1,638238 0,5460795 2,730397 2,002291 1,456212
5,278713 1,274172 2,912394 1,638221 0,5460738 2,730369 2,002271 1,456197
5,278703 1,274170 2,912388 1,638218 0,5460727 2,730364 2,002267 1,456194
5,278701 1,274169 2,912387 1,638218 0,5460726 2,730363 2,002266 1,456193
5,278701 1,274169 2,912387 1,638218 0,5460725 2,730363 2,002266 1,456193
5,278701 1,274169 2,912387 1,638218 0,5460725 2,730363 2,002266 1,456193
5,278701 1,274169 2,912387 1,638218 0,5460725 2,730363 2,002266 1,456193
[,89] [,90] [,91] [,92] [,93] [,94] [,95] [,96]
tau 1,758588 3,588089 2,502118 2,972330 0,903616 3,944375 2,485937 0,6793552
1,000000 1,000000 1,000000 1,000000 1,000000 1,000000 1,000000 1,0000000
2,185102 5,098570 3,641836 2,367193 1,820918 5,826938 1,820918 0,3641836
2,184440 5,097027 3,640734 2,366477 1,820367 5,825174 1,820367 0,3640734
2,184318 5,096742 3,640530 2,366344 1,820265 5,824848 1,820265 0,3640530
2,184295 5,096689 3,640492 2,366320 1,820246 5,824787 1,820246 0,3640492
2,184291 5,096679 3,640485 2,366315 1,820242 5,824776 1,820242 0,3640485
2,184290 5,096677 3,640484 2,366314 1,820242 5,824774 1,820242 0,3640484
2,184290 5,096677 3,640483 2,366314 1,820242 5,824774 1,820242 0,3640483
2,184290 5,096677 3,640483 2,366314 1,820242 5,824773 1,820242 0,3640483
2,184290 5,096677 3,640483 2,366314 1,820242 5,824773 1,820242 0,3640483
[,97] [,98] [,99] [,100]
tau 1,864360 0,5715798 1,022274 1,760623
1,000000 1,0000000 1,000000 1,000000
1,820918 0,3641836 1,456734 3,459744
1,820367 0,3640734 1,456294 3,458697
1,820265 0,3640530 1,456212 3,458503
1,820246 0,3640492 1,456197 3,458467
1,820242 0,3640485 1,456194 3,458461
1,820242 0,3640484 1,456193 3,458459
1,820242 0,3640483 1,456193 3,458459
1,820242 0,3640483 1,456193 3,458459
1,820242 0,3640483 1,456193 3,458459
In [4]:
# rm(list=ls())
#==============================================================================
# graphics.off()
#==============================================================================
# Fim
#==============================================================================