In [1]:
# rm(list=ls())
options(OutDec = ",")
#==============================================================================
# Exemplo para encontrar o maximo de uma funcao envolvendo senos e cosenos
# pelo metodo de arrefecimento simulado (simulated anneling)
#==============================================================================
# Funcao do metodo de arrefecimento simulado
#==============================================================================
maxf <- function(fx,x0,temperatura,Delta,maxit,percent=0.95){
xmax <- rep(NA,maxit)
x <- x0
for(k in 1:maxit){
xprop <- runif(1,x-Delta,x+Delta)
prob <- min(1,exp((fx(xprop)-fx(x))/temperatura))
if(runif(1) < prob){
xmax[k] <- xprop
x <- xprop
}else{
xmax[k] <- x
}
temperatura <- percent*temperatura
}
return(xmax)
}
In [2]:
#==============================================================================
# Funcao de uma mistura discreta de duas normais
#==============================================================================
fx <- function(x){
0.3*dnorm(x,10,1)+0.7*dnorm(x,15,1.5)
}
In [3]:
#==============================================================================
# Encontrando o maximo em uma grade de valores
#==============================================================================
omega <- 0.3
mu1 <- 10
sigma1 <- 1
mu2 <- 15
sigma2 <- 1.5
xseq <- seq(5,22,length=100000)
yseq <- ( omega*dnorm(xseq,mu1,sigma1)
+ (1-omega)*dnorm(xseq,mu2,sigma2) )
xmax <- xseq[yseq==max(yseq)]
print(xmax)
print(max(yseq))
[1] 15,00001 [1] 0,1861735
In [4]:
#==============================================================================
# Grafico
#==============================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(14,6,0),
mar=c(4.5,5,2,1),cex.main=2.0,bty="n")
curve(fx(x),from=5,to=22,lwd=4,ylab=expression(f(x)),xlim=c(5,23),
ylim=c(0,0.20))
lines(c(xmax,xmax),c(0,max(yseq)),lty=2,lwd=4,col="red")
In [5]:
#==============================================================================
# Encontrando o maximo pelo arrefecimento simulado
#==============================================================================
set.seed(12345)
xmax.1 <- maxf(fx, 7,100,0.5,500)
print(xmax.1)
[1] 7,220904 7,481886 7,438367 7,263463 7,491168 7,025703 7,261388 [8] 7,152592 7,040736 6,719699 6,673427 7,138843 7,283385 7,481929 [15] 7,208396 7,501403 7,189116 7,059220 7,428015 7,545439 7,827632 [22] 8,254906 8,014588 7,574783 7,129837 7,594307 7,409335 7,641831 [29] 7,871603 7,807134 8,098702 8,584685 9,064464 9,513171 9,613528 [36] 9,801881 9,674817 9,223069 9,684516 9,694808 10,065256 9,573904 [43] 9,218416 9,544073 9,847645 10,275600 9,854414 10,069194 10,289310 [50] 9,884951 9,679416 10,153690 10,175059 10,312514 10,063631 10,173108 [57] 10,428379 10,723351 11,207377 10,716841 10,776298 10,366009 10,575982 [64] 11,020060 10,798702 10,332150 10,644917 11,089392 11,286763 10,884307 [71] 10,906085 10,622414 10,269422 10,759533 11,179180 11,506238 11,801909 [78] 11,510437 11,693488 12,039922 12,280582 12,632141 12,512565 13,005942 [85] 13,340223 13,809279 14,058501 13,937560 14,141761 14,437039 14,104586 [92] 13,800839 14,271255 14,118851 13,782238 14,036301 13,843066 14,327591 [99] 14,101911 13,845581 13,934173 13,557968 13,721609 13,314235 13,370333 [106] 13,166296 13,043966 13,484125 13,311189 13,437131 13,046447 12,639189 [113] 13,035578 13,442032 13,409851 13,216599 13,195802 13,431082 13,916453 [120] 13,888941 14,327752 14,734576 14,860907 15,298017 15,134807 14,758072 [127] 14,290564 14,666994 14,183848 13,799690 14,249415 14,560974 14,761139 [134] 14,379215 14,875259 14,522598 14,522598 14,770825 14,822217 14,658519 [141] 14,285607 13,883489 13,659511 13,659511 13,665836 13,871296 13,432193 [148] 13,734434 14,049311 14,104464 14,104464 14,312754 14,390168 14,100146 [155] 14,056897 13,957824 14,413197 14,413197 14,440214 14,440214 14,253117 [162] 14,372229 14,372229 14,187227 13,847874 13,847874 14,084513 13,692034 [169] 14,153194 14,143748 14,222977 14,634018 14,648406 14,943601 14,943601 [176] 15,128325 15,128325 15,284250 15,369435 14,904142 14,585095 14,565537 [183] 14,565537 14,936379 15,228976 15,259433 15,259433 15,405276 15,337425 [190] 15,145184 15,145184 14,951898 14,951898 15,070172 15,475112 15,475112 [197] 15,475112 15,303789 15,131308 15,136616 15,136616 14,901205 14,997405 [204] 14,997405 15,294516 15,294516 14,835296 15,075912 15,114590 15,084008 [211] 15,084008 14,779109 14,779109 15,259707 15,030123 15,153975 15,153975 [218] 15,153975 15,153975 15,181459 15,181459 15,181459 15,121766 15,015815 [225] 15,015815 15,015815 14,932921 14,937781 14,937781 14,968586 14,999846 [232] 14,999846 15,065987 15,065987 15,065987 14,885920 14,885920 14,885920 [239] 15,120810 15,120810 14,862750 14,862750 14,862750 14,862750 14,862750 [246] 15,142427 15,142427 15,108982 15,108982 15,108982 15,108982 15,108982 [253] 15,108982 15,108982 15,004369 15,004369 15,004369 15,004369 15,004369 [260] 15,004369 15,004369 15,004369 14,947454 14,947454 14,947454 14,947454 [267] 14,964089 14,964089 14,964089 14,964089 14,964089 14,964089 14,964089 [274] 14,964089 14,964089 14,964089 14,969991 14,969991 14,969991 14,969991 [281] 14,969991 15,035236 15,035236 15,035236 15,035236 15,035236 15,035236 [288] 15,035236 15,035236 15,035236 15,035236 15,035236 15,035236 15,035236 [295] 15,035236 15,035236 15,035236 15,035236 15,035236 15,035236 15,035236 [302] 15,035236 15,035236 15,035236 14,988480 15,003720 15,003720 15,003720 [309] 15,003720 15,003720 15,003720 15,003720 15,003720 15,003720 15,003720 [316] 15,003720 15,003720 15,003720 15,003720 15,003720 15,003720 15,003720 [323] 15,003720 15,003720 15,003720 15,003720 15,003720 15,003720 15,003720 [330] 15,003720 15,003720 15,003720 15,003720 15,003720 15,003720 15,003720 [337] 15,003720 15,003720 15,003720 15,003720 15,003720 15,003720 15,003720 [344] 15,003720 15,003720 15,003720 15,003720 15,003720 15,003720 15,003720 [351] 15,003720 15,003720 15,003720 15,003720 15,003720 15,003720 15,003720 [358] 15,003720 15,003720 15,003720 15,003720 15,003720 15,003720 15,003720 [365] 15,003720 15,003720 15,003720 15,003720 15,003720 15,003720 15,003720 [372] 15,003720 15,003720 15,003720 15,003720 15,003720 15,003720 15,003720 [379] 15,003720 15,003720 15,003720 15,003720 15,003720 15,003720 15,003720 [386] 15,003720 15,003720 15,003720 15,003720 15,003720 15,003720 15,003720 [393] 15,003720 15,003720 15,003720 15,003720 15,003720 15,003720 15,003720 [400] 15,003720 15,003720 15,004575 15,004575 15,004575 15,004575 15,004575 [407] 15,004575 15,004575 15,004575 15,004575 15,004575 15,004575 15,004575 [414] 15,004575 15,004575 14,996629 14,996629 14,996629 14,996629 14,996629 [421] 14,996629 14,996629 14,996629 14,996629 14,996629 14,996629 14,996629 [428] 14,996629 14,996629 14,996629 14,996629 14,996629 14,996629 14,996629 [435] 14,996629 14,996629 14,996629 14,996629 14,996629 14,996629 14,996629 [442] 14,996629 14,996629 14,996629 14,996629 14,996629 14,996629 14,996629 [449] 14,996629 14,996629 14,996629 14,996629 14,996629 14,996629 14,996629 [456] 14,996629 14,996629 14,996629 14,996629 14,996629 14,996629 14,996629 [463] 14,996629 14,996629 14,996629 14,996629 14,996629 14,996629 14,996629 [470] 14,996629 14,996629 14,996629 14,996629 14,996629 14,996629 14,996629 [477] 14,996629 14,996629 14,996629 14,996629 14,996629 14,996629 14,996629 [484] 14,996629 14,996629 14,996629 14,996629 14,996629 14,996629 14,996629 [491] 14,996629 14,996629 14,996629 14,996629 14,996629 14,996629 14,996629 [498] 14,996629 14,996629 14,996629
In [6]:
#==============================================================================
# Encontrando o maximo pelo arrefecimento simulado
#==============================================================================
xmax.2 <- maxf(fx,18,100,0.5,500)
print(xmax.2)
[1] 17,57776 17,08457 17,27063 16,94476 17,44433 17,37961 17,46685 17,77639 [9] 17,29221 17,19980 16,72263 16,33109 16,63240 16,61385 16,77062 17,04921 [17] 17,35627 16,96108 16,89839 16,70261 16,86281 17,07652 17,34262 17,18165 [25] 17,19825 17,28699 17,57117 17,40719 17,26380 16,80241 16,72960 17,02770 [33] 16,80624 16,66606 17,14881 17,38948 17,26688 17,64896 17,85910 17,83766 [41] 18,32097 18,24653 18,28586 18,75782 19,20536 19,68719 19,43158 18,94268 [49] 18,49158 18,70707 18,52581 18,51257 18,65962 19,03613 19,45517 19,93810 [57] 19,57652 19,96453 20,21540 19,94920 19,99413 20,02012 19,78427 19,98234 [65] 19,69334 19,44108 19,42626 19,31321 18,96565 19,38464 19,68326 19,55400 [73] 19,75761 19,37127 19,60652 19,82427 20,21608 20,39888 19,95980 19,98191 [81] 20,23368 19,82216 19,91897 19,79826 20,06598 20,14578 19,70335 19,71274 [89] 20,00136 19,50666 19,30579 18,90041 19,18880 19,45091 19,79985 20,25355 [97] 20,16403 20,08558 20,43727 20,26661 19,81746 19,90124 19,72226 19,47613 [105] 19,42325 19,01638 19,04695 18,89707 19,22533 19,23398 19,05735 19,14786 [113] 18,72916 19,05149 19,07789 19,36687 19,34274 19,37833 19,53423 19,22135 [121] 19,64771 19,72912 19,70674 19,46788 19,73143 19,70572 19,75839 19,77224 [129] 19,88314 20,20416 19,85270 19,76812 20,23041 20,65840 20,51612 20,28601 [137] 20,47603 20,57304 20,52216 20,12773 20,57663 20,91196 21,27280 21,13245 [145] 20,78738 21,19653 21,13471 21,45765 21,77071 21,82841 21,53490 21,14298 [153] 21,63273 21,65296 21,76737 21,30242 20,92849 21,32837 20,91537 21,38565 [161] 20,97122 21,17797 20,78583 21,25558 20,82969 20,76875 20,44090 20,56353 [169] 20,67177 20,54798 20,39834 19,97485 20,30718 20,53750 20,39373 19,92911 [177] 20,34147 20,51355 20,57900 20,59096 20,34044 20,71247 21,19745 20,83552 [185] 20,47059 20,87424 20,37697 20,51054 20,42209 20,34331 20,79211 20,32288 [193] 20,22414 19,80936 20,23601 19,75227 19,53415 19,62761 19,99191 20,20461 [201] 20,69399 21,17234 21,14927 21,32006 20,91629 20,90224 21,22456 21,08368 [209] 21,26349 21,18374 21,18161 20,97863 20,78176 20,74792 20,91648 21,25733 [217] 20,87971 21,33332 21,72154 21,71059 21,59591 21,93054 21,83582 21,51942 [225] 21,28339 20,83536 20,35116 20,83384 20,93991 20,74163 20,73535 20,67524 [233] 20,71039 20,99237 21,14289 21,27895 20,80770 20,78470 20,60341 20,84061 [241] 20,50895 20,21941 20,02937 19,86717 19,62341 19,62341 19,62341 19,62341 [249] 19,46972 19,24842 19,24842 19,28714 19,28714 19,01414 19,01414 18,75716 [257] 18,35794 18,32112 18,32112 17,99740 17,64872 17,24222 17,20867 17,20867 [265] 17,14513 17,14513 17,14513 17,14513 17,14513 17,14513 17,14513 17,05124 [273] 16,87659 16,61336 16,61336 16,61336 16,61336 16,61336 16,40107 16,40107 [281] 16,40107 16,05846 16,05846 16,05846 16,05846 16,05846 15,91273 15,91273 [289] 15,91273 15,91273 15,91273 15,91273 15,91273 15,67792 15,67792 15,67792 [297] 15,28246 15,28246 15,28246 15,28246 14,98928 14,98928 14,98928 14,98928 [305] 14,98928 14,98928 14,98928 14,98928 14,98928 14,98928 14,98928 14,98928 [313] 14,98928 14,98928 14,98928 14,98928 14,98928 14,98928 14,98928 14,98928 [321] 14,98928 14,98928 14,98928 14,98928 14,98928 14,98928 14,98928 14,98928 [329] 14,98928 14,98928 14,98928 14,98928 14,98928 14,98928 14,98928 14,98928 [337] 14,98928 14,98928 14,98928 14,98928 14,98928 14,98928 14,98928 14,98928 [345] 14,98928 14,98928 14,98928 14,99972 14,99972 14,99972 14,99972 14,99972 [353] 14,99972 14,99972 14,99972 14,99972 14,99972 14,99972 15,00219 15,00219 [361] 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 [369] 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 [377] 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 [385] 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 [393] 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 [401] 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 [409] 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 [417] 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 [425] 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 [433] 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 [441] 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 [449] 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 [457] 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 [465] 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 [473] 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 [481] 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 [489] 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 15,00219 [497] 15,00219 15,00219 15,00219 15,00219
In [7]:
#==============================================================================
# Grafico: iteracoes
#==============================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(5,14,5),
mar=c(4.5,5,1,1),cex.main=2.0,bty="n")
ts.plot(xmax.1,ylim=c(5,23),lwd=4,col="black",xlab="Iteração",
ylab=expression(g(x[i])))
lines(xmax.2,col="red",lwd=4)
In [8]:
# rm(list=ls())
#==============================================================================
# graphics.off()
#==============================================================================
# Fim
#==============================================================================