Simulated Annealing (Monte Carlo) sloužící k dekompozici cen balíčků televizních stanic

dekompozice_cen_TV_balicku.R

####################################################################
# Simulated Annealing je Monte Carlo metoda, která vychází z 
# fyzikálního pozorování - při kalení oceli je vhodné její
# teplotu snižovat postupně, protože molekuly při pomalém tuhnutí
# zaujímají prostorovou konfiguraci s minimální energetickou
# náročností vazeb - to se projevuje v pevnosti (a křehkosti kovu)
# annealing si půjčuje myšlenku, že postupným zkracováním 
# vzdálenosti pátrače od dosud nejlepšího lokálního minima se 
# velmi často dosáhne globálního minima, nebo alespoň minima,
# které je mu blízké.
###################################################################
# skript slouží k odhadu cen jednotlivých stanic schovaných v 
# balíčcích poskytovatelů televizního vysílání
# poskytovatelé neuvádějí cenu za stanici, ale pouze za celý balíček
# hlavní účel skriptu je stanovit cenové hladiny stanic pro 
# ACBC conjoint analýzu
#################################################################### 
setwd("./dekompozice cen")
for (c in 1:10) { # cyklus pro vícenásobné spuštění annealingu
#################################################################
# ceny.csv jsou ceny jednotlivých balíčků, ke stažení jsou zde: 
# https://dl.dropboxusercontent.com/u/5817117/ceny.csv
# balicky.csv jsou zastoupeni stanic v baliccich, ke stažení zde:
# https://dl.dropboxusercontent.com/u/5817117/balicky.csv
#################################################################
ceny<-read.csv2("ceny.csv",as.is=T)
ceny<-ceny[order(ceny$balicek),]
balicky<-read.csv2("balicky.csv",as.is=T)
balicky[is.na(balicky)]<-0
#set.seed(1234) # pro zopakovani simulace, pokud chcete pustit více simulací za sebou, tak je to třeba vypnout, jinak se spočítá pokaždé stejně
 
maxcena<-200 # maximální cena jedné stanice
# pro jednoduchost je minimální cena vždy 0
step<-2 # maximální krok, kterým se annealing odchyluje od současného lokálního optima
pokusu<-100000 # kolik kroků simulace se provede 100 000 už dává velmi dobré výsledky, pro hrubou orientaci stačí 1000

####################################################
# varianta "logaritmické schody"
# an<-sort(rep(2:round(step,0),pokusu),decreasing=T)
# an<-round(log(an),0)
####################################################

###########################
# varianta "pevného kroku"
#an<-rep(1,pokusu)
###########################

##############################################
# varianta logaritmického poklesu kroku
# sample si s tím poradí tak, že část veličin 
#generuje zaokrouhlených nahoru a část dolu 
#(asi v proporci, ale nemám to ověřené)
##############################################
an<-(log(pokusu:1)/log(pokusu)*step)
an[an<1]<-1
plot(an) #graf jak se zmenšuje krok = tuhne ocel
balicku<-dim(balicky)[1]
asum<-1e10 # něco velkého, co určitě bude větší, než součet odchylek annealingu od reálných balíčků
asum_e<-1e10 # dtto
balicky$best<-round(maxcena/2) #pro balicky, kde je vice nulovych cen je lepsi zadat vychozi hodnotu best<-0

pb <- txtProgressBar(style=3) # inicializace progressbaru
j<-0 # v kolikátém jsem kroku cyklu
delka_cyklu<-length(an) # kolik kroků cyklus má
for (i in an) {
	j<-j+1
	p<-j/delka_cyklu # procento pro porgress bar
	setTxtProgressBar(pb, p) # update progressbaru
	print(paste("asum",asum,"asum_e",asum_e)) # jak daleko je lokální optimum od předlohy
	
	################################################################
	# tohle je vlastní annealing - k optimu přičtu náhodný vektor
	# inkrementů, které se ((obvykle) logaritmicky) zmenšují v čase
	# podstata annealingu je, že na začátku je dobrodruh a na
	# konci je konzerva, která se drží u posledního optima
	################################################################
	balicky$try<-balicky$best+((sample(2,balicku,T)-1)*2-1)*(sample(i,balicku,replace=T))
	############################################################
	# tohle je z lenosti - když annealingem vyjedu z hranic, 
	# tak to oříznu zpátky
	############################################################
	balicky$try[balicky$try<0]<-0
	balicky$try[balicky$try>maxcena]<-maxcena
	############################################################
	# průběžné statistiky o tom, jak vypadá rozložení cen stanic 
	# a taky, jaký je aktuální krok
	############################################################
	print(paste(min(balicky$try),median(balicky$try),max(balicky$try),"@",i))
	############################################################
	# tohle by asi chtělo uklidit do funkce, je to výpočet
	# vzdálenosti simulovaného balíčku od předlohy a ukládání
	# dosud nejlepšího výsledku
	############################################################
	csum<-as.data.frame(colSums(balicky[,2:(length(balicky)-1)]*balicky$try))
	csum<-data.frame("balicek"=rownames(csum),"try"=csum[,1])
	csum<-csum[order(csum$balicek),]
	mceny<-merge(ceny,csum,by="balicek")
	##############################################################
	# vypadá to, že s eukleidovskou vzdáleností (diff_e) metoda
	# lépe konverguje, nicméně počítám obě, absolutní hodnoty jsou
	# poměrně názorné pro laiky
	##############################################################
	diff<-abs(mceny$cena-mceny$try)
	diff_e<-diff*diff
	if (sqrt(sum(diff_e))<asum_e) {
		ceny$abest<-mceny$try
		asum<-sum(diff)
		asum_e<-sqrt(sum(diff_e))
		ceny$diff<-diff
		ceny$diff_e<-diff_e
		balicky$best<-balicky$try
		}
	}
	close(pb) # konec progressbaru annealing doběhl

##################################
# zobrazení a uložení výsledků
##################################
ceny
show<-data.frame("balicek"=balicky[,1],"cena"=balicky$best)
show<-show[order(show[,2]),]
show
asum
asum_e
write.csv2(show,paste("dekompozice_cen_",asum,"_",c,".csv",sep=""))
write.csv2(ceny,paste("dekompozice_balicky_",asum,"_",c,".csv",sep=""))
alarm() # akustická signál, že je hotovo
} # ukončení cyklu, pokud annealing spouštíme opakovaně, číslo běhu se ukládá do názvu souboru