Entropia
Jak wiadomo, oprócz swojego termodynamicznego znaczenia, entropia jest równie¿ wykorzystywana w teorii informacji, gdzie przyjmuje nastêpuj±c± postaæ:
\(S = -\sum\limits_i p(x_i)\log p(x_i)\),
przy czym \(p(x_i)\) to po prostu prawdopodobieñstwo zaj¶cia zdarzenia \(x_i\).
W ramach biblioteki entropy (po uprzednim zainstalowaniu komend± install.packages("entropy")) warto¶æ entropii dla liczby zliczeñ otrzymujemy za pomoc± funkcji entropy().
## [1] 0## [1] 0.6931472## [1] 0.6931472## [1] 0## [1] 1.098612## [1] 0.6931472## [1] 1.098612Charakterystyczn± cech± tej funkcji jest przyjmowanie warto¶ci maksymalnej dla rozk³adu jednorodnego (tzn. gdy prawdopodobieñstwo ka¿dego zadrzenia jest takie samo), natomiast minimalnej w przypadku, gdy wszystkie obserwacje s± jednego typu.
# PRZYK£AD 8.2
library(ggplot2)
N <- 100
n1 <- 0:N
df.ent <- data.frame(n1 = n1, S = apply(cbind(n1, N - n1), 1, entropy))
df.ent## n1 S
## 1 0 0.00000000
## 2 1 0.05600153
## 3 2 0.09803911
## 4 3 0.13474217
## 5 4 0.16794415
## 6 5 0.19851524
## 7 6 0.22696752
## 8 7 0.25363895
## 9 8 0.27876937
## 10 9 0.30253782
## 11 10 0.32508297
## 12 11 0.34651534
## 13 12 0.36692499
## 14 13 0.38638671
## 15 14 0.40496349
## 16 15 0.42270909
## 17 16 0.43966988
## 18 17 0.45588621
## 19 18 0.47139349
## 20 19 0.48622296
## 21 20 0.50040242
## 22 21 0.51395667
## 23 22 0.52690796
## 24 23 0.53927634
## 25 24 0.55107993
## 26 25 0.56233514
## 27 26 0.57305692
## 28 27 0.58325884
## 29 28 0.59295332
## 30 29 0.60215168
## 31 30 0.61086430
## 32 31 0.61910066
## 33 32 0.62686946
## 34 33 0.63417864
## 35 34 0.64103548
## 36 35 0.64744664
## 37 36 0.65341819
## 38 37 0.65895568
## 39 38 0.66406413
## 40 39 0.66874809
## 41 40 0.67301167
## 42 41 0.67685855
## 43 42 0.68029200
## 44 43 0.68331491
## 45 44 0.68592980
## 46 45 0.68813881
## 47 46 0.68994376
## 48 47 0.69134610
## 49 48 0.69234697
## 50 49 0.69294717
## 51 50 0.69314718
## 52 51 0.69294717
## 53 52 0.69234697
## 54 53 0.69134610
## 55 54 0.68994376
## 56 55 0.68813881
## 57 56 0.68592980
## 58 57 0.68331491
## 59 58 0.68029200
## 60 59 0.67685855
## 61 60 0.67301167
## 62 61 0.66874809
## 63 62 0.66406413
## 64 63 0.65895568
## 65 64 0.65341819
## 66 65 0.64744664
## 67 66 0.64103548
## 68 67 0.63417864
## 69 68 0.62686946
## 70 69 0.61910066
## 71 70 0.61086430
## 72 71 0.60215168
## 73 72 0.59295332
## 74 73 0.58325884
## 75 74 0.57305692
## 76 75 0.56233514
## 77 76 0.55107993
## 78 77 0.53927634
## 79 78 0.52690796
## 80 79 0.51395667
## 81 80 0.50040242
## 82 81 0.48622296
## 83 82 0.47139349
## 84 83 0.45588621
## 85 84 0.43966988
## 86 85 0.42270909
## 87 86 0.40496349
## 88 87 0.38638671
## 89 88 0.36692499
## 90 89 0.34651534
## 91 90 0.32508297
## 92 91 0.30253782
## 93 92 0.27876937
## 94 93 0.25363895
## 95 94 0.22696752
## 96 95 0.19851524
## 97 96 0.16794415
## 98 97 0.13474217
## 99 98 0.09803911
## 100 99 0.05600153
## 101 100 0.00000000
Istotne jest równie¿ to, ¿e entropia pokazuje kierunek zachodzenia zmian w uk³adzie - mówi o tym tzw. zasada maksymalizacji entropii. Zgodnie z ni±, wszystkie procesy (niewymuszane z zewn±trz) prowadz± do maksymalizacji warto¶ci entropii w uk³adzie, co de facto oznacza, ¿e uk³ad d±¿y do stanu równowagi.
W poni¿szym przyk³adzie mamy do czynienia z modelem Isinga na ³añcuchu - bardzo prostym modelem magnetyka, zaproponowanym w latach 20-tych XX w. W uk³adzie tym w ka¿dym wê¼le ³añcuha znajduje siê spin, który mo¿e przyj±æ jedn± z dwóch orientacji - w górê lub w dó³. W rozpatrywanym przez nas przypadku spin oddzia³uje jedynie z swoimi najbli¿szymi s±siadem a si³a oddzia³ywania jest sta³a i wynosi \(J\). W efekcie hamiltonian uk³adu to
\(H = -J \sum\limits_i s_i s_{i+1}\)
Symulacje startujemy z ca³kowicie uporz±dkowanego ³añcucha i dokonujemy termalizacji za pomoc± algorytmu Metropolisa.
# PRZYK£AD 8.3
# Hamiltonian modelu Isinga
H <- function(S, i, J) {
-J * S[i]*(S[i - 1] + S[i + 1])
}
# Algorytm Metropolisa
metro <- function(S, i, J, T) {
testS <- S
testS[i] <- -testS[i]
# Zmiana energii
dE <- H(testS, i, J) - H(S, i, J)
# Testowanie zmiany energii
ifelse(dE < 0, S <- testS, ifelse(exp(-dE/(kB * T)) > runif(1), S <- testS, S <- S))
S
}
# stale
J <- 1
kB <- 1
# parametry
N <- 50
T <- 5
h <- 0
# liczba krokow
M <- 1000
# macierz kolejnych zmian
mS <- matrix(0, M, N)
# poczatkowe ustawienie spinow
S0 <- S <- rep(1, N)
# glowna petla termalizacji
for(j in 1:M) {
S <- metro(S, sample(2:(N-1), 1), J, T)
mS[j,] <- S
}
plus <- sapply(1:M, function(x) sum(mS[x,] == 1))
S.df <- data.frame(n = 1:M, ent = sapply(1:M, function(x) entropy(c(plus[x], N - plus[x]))))Nastêpnie wizualizujemy stan pocz±tkowy i koñcowy ³añcucha.
# PRZYK£AD 8.4
library(grid)
S.plot <- data.frame(x = 1:length(S), y = -S + 1, yend = S + 1, spin = S)
S0.plot <- data.frame(x = 1:length(S0), y = -S0 + 1, yend = S0 + 1, spin = S0)
g <- ggplot()
s0 <- geom_segment(data = S0.plot, aes(x = x, y = y, xend = x, yend = yend, colour=spin), arrow = arrow(length=unit(0.5,"cm")), size=1.5)
s <- geom_segment(data = S.plot, aes(x = x, y = y-4, xend = x, yend = yend-4, colour=spin), arrow = arrow(length=unit(0.5,"cm")), size=1.5)
sc <- scale_y_continuous(limits = c(-7,5.5))
g + s0 + s + sc
Na koniec sprawdzamy w jaki sposób zmienia³a siê entropia uk³adu w czasie.
Informacja wzajemna
Na bazie pojêcia entropii powsta³o wiele miar teorioinformacyjnych, które maj± za zadanie wskazywaæ mo¿liwo¶æ uzyskania informacji o okre¶lonym rozk³adzie (lub zmiennej) pod warunkiem posiadania wiadomo¶ci o innym. Jedn± z takich miar jest informacja wzajemna (mutual information), zdefiniowana dla dyskretnych zmiennych losowych \(X\) i \(Y\) jako
\(I(X, Y) = \sum\limits_{y \in Y}\sum\limits_{x \in X} p(x,y) \log\frac{p(x,y)}{p(x)p(y)}\),
przy czym \(p(x,y)\) to ³±czny rozk³ad prawdopodobieñstwa \(X\) i \(Y\), natomiast \(p(x)\) i \(p(y)\) to rozk³ady brzegowe (czyli rozk³ady zmiennych \(X\) i \(Y\)).
Informacja wzajemna okre¶la, ile informacji o \(X\) mo¿na poznaæ, znaj±c \(Y\), czyli o ile poznanie jednej z tych zmiennych zmniejsza niepewno¶æ o drugiej. Je¶li zmienne \(X\) i Y s± niezale¿ne, to ich wzajemna informacja jest zerowa (znajomo¶æ jednej nie mówi niczego o drugiej). Je¶li \(X\) i \(Y\) s± identyczne, to ka¿da zawiera pe³n± wiedzê o drugiej. Wtedy informacja wzajemna jest równa entropii \(X\) lub \(Y\) (tekst z https://pl.wikipedia.org/wiki/Informacja_wzajemna).
W R informacja wzajemna poprzez funkcjê mi.empirical() z biblioteki entropy. Jako argument wymagana jest tabela kontygencji. W pierwszym przypadku liczymy informacjê wzajemn± dla powtórzonych danych (tutaj: rozk³ad jednorodny).
# PRZYK£AD 8.6
x <- sample(20, 10000, replace=TRUE)
y <- x
T1 <- table(x, y)
T1.mi <- mi.empirical(T1)
heatmap(T1, Rowv = NA, Colv = NA)
Nastêpnie modyfikujemy trochê uk³ad poprzez dodanie/odjêcie warto¶ci 1.
# PRZYK£AD 8.7
x <- sample(2:19, 10000, replace=TRUE)
y <- x + sample(c(-1,0,1), 10000, replace=TRUE)
T2 <- table(x, y)
T2.mi <- mi.empirical(T2)
heatmap(T2, Rowv = NA, Colv = NA)
Wreszcie liczymy informacjê wzajemn± dwóch nieskorelowanych ze sob± rozk³adów.
# PRZYK£AD 8.8
x <- sample(20, 10000, replace=TRUE)
y <- sample(20, 10000, replace=TRUE)
T3 <- table(x, y)
T3.mi <- mi.empirical(T3)
heatmap(T3, Rowv = NA, Colv = NA)
Warto na koñcu dokonaæ porównania warto¶ci dla poszczególnych rozk³adów
| Przypadek 1 | Przypadek 2 | Przypadek 3 |
|---|---|---|
| 2.9945337 | 1.8639747 | 0.0176282 |
oraz unormowania ich przez warto¶æ maksymaln±, czyli w tym przypadku \(\log 20\)
| Przypadek 1 | Przypadek 2 | Przypadek 3 |
|---|---|---|
| 0.9995999 | 0.62221 | 0.0058844 |
Dywergencja Kullbacka-Leibera
Czêsto konieczne jest porównanie "blisko¶ci" rozk³adów prawdopodobieñstwa dwóch zmiennych - jedn± z wykorzystywanych miar jest dywergencja Kullbacka-Leibera, która dla dwóch dyskretnych rozk³adów prawdopodobieñstwa jest zdefiniowana jako
\(D_{KL} (P \parallel Q) = -\sum\limits_i P(i) \log \frac{Q(i)}{P(i)}\)
Dywergencja Kullbacka-Leibera jest realizowana w R poprzez funkcjê KLdiv() z biblioteki flexmix (do zainstalowania). W poni¿szym przyk³adzie liczmy dywergencjê KL dla dwóch rozk³adów Gaussa, przy czym istotne jest, ¿e ogólnie miara ta nie jest symetryczna.
## Loading required package: latticex <- seq(-10, 10, 0.01)
y1 <- dnorm(x, 0, 1)
y2 <- dnorm(x, 0, 2)
df <- data.frame(x = x, y1 = y1, y2 = y2, KL1 = log(y1/y2)*y1, KL2 = log(y2/y1)*y2)
g <- ggplot(df, aes(x = x))
l1 <- geom_line(aes(y = y1), colour="red")
l2 <- geom_line(aes(y = y2), colour="blue")
a1 <- geom_area(aes(y = KL1), fill="red", alpha=0.2)
a2 <- geom_area(aes(y = KL2), fill="blue", alpha=0.2)
theme_set(theme_bw())
g + l1 + l2 + a1 + a2
## y1 y2
## y1 0.0000000 0.3181466
## y2 0.8058743 0.0000000Dywergencja KL jest ca³k± z zacieniowanych obszarów. Poni¿ej podobny przyk³ad, ilustruj±cy jak zwiêksza siê warto¶æ dywergencji wraz ze wzajemnym przesuwaniem siê rozk³adów.
# PRZYK£AD 8.10
x <- seq(-6, 6, 0.01)
y1 <- dnorm(x, 0, 1)
y2 <- dnorm(x, 1, 1)
y3 <- dnorm(x, 2, 1)
y4 <- dnorm(x, 3, 1)
df <- data.frame(x = x, y1 = y1, y2 = y2, y3 = y3, y4 = y4, KL1 = log(y1/y2)*y1, KL2 = log(y1/y3)*y1, KL3 = log(y1/y4)*y1)
g <- ggplot(df, aes(x = x))
l1 <- geom_line(aes(y = y1), colour="black")
l2 <- geom_line(aes(y = y2), colour="red")
l3 <- geom_line(aes(y = y3), colour="blue")
l4 <- geom_line(aes(y = y4), colour="darkgreen")
KL1 <- geom_area(aes(y = KL1), fill="red", alpha=0.2)
KL2 <- geom_area(aes(y = KL2), fill="blue", alpha=0.2)
KL3 <- geom_area(aes(y = KL3), fill="darkgreen", alpha=0.2)
g + l1 + l2 + l3 + l4 + KL3 + KL2 + KL1
## y1 y2 y3 y4
## y1 0.0000000 0.4999997 1.9999667 4.4984514
## y2 0.4999988 0.0000000 0.4999707 1.9986722
## y3 1.9997687 0.4998996 0.0000000 0.4988334
## y4 4.4882143 1.9925855 0.4969263 0.0000000Dywergencjê mo¿na wykorzystaæ równie¿ do bardzo "egzotycznych" zadañ. W poni¿szym przyk³adzie mamy do czyniania z trzema tekstami, pochodz±cymi odpowiednio z dwóch opowiadañ kryminalnych The Murders in the Rue Morgue Edgara Allana Poe'a (1841), The Adventure of the "Gloria Scott" Sir Arthura Conan Doyle'a (1893) oraz pracy naukowej On the Origin of Species Karola Darwina (1859). Ka¿da pozycja zawiera oko³o 1000 s³ów, które najpierw zostaj± oczyszczone ze znaków interpunkcyjnych, potem sprowadzone do ma³ych liter, a na koniec wyznaczany jest histogram ich wystêpowania. Efektywne rozk³ady prawdopodobieñstwa s³u¿a do wyznaczenia dywergencji KL pomiêdzy poszczególnymi pozycjami.
# PRZYK£AD 8.11
rue <- read.table("http://www.fizyka.pw.edu.pl/~julas/CSAR/rue.txt")
rue.t <- t(rue)
rue.t <- sapply(1:length(rue.t), function(x) gsub("([.,:;'()?!])", "", rue.t[x]))
rue.t <- sapply(1:length(rue.t), function(x) tolower(rue.t[x]))
rue.t.tab <- as.data.frame(table(rue.t))
colnames(rue.t.tab) <- c("word", "freq")
gloria <- read.table("http://www.fizyka.pw.edu.pl/~julas/CSAR/gloria.txt")
gloria.t <- t(gloria)
gloria.t <- sapply(1:length(gloria.t), function(x) gsub("([.,:;'()?!])", "", gloria.t[x]))
gloria.t <- sapply(1:length(gloria.t), function(x) tolower(gloria.t[x]))
gloria.t.tab <- as.data.frame(table(gloria.t))
colnames(gloria.t.tab) <- c("word", "freq")
origin <- read.table("http://www.fizyka.pw.edu.pl/~julas/CSAR/origin.txt")
origin.t <- t(origin)
origin.t <- sapply(1:length(origin.t), function(x) gsub("([.,:;'()?!])", "", origin.t[x]))
origin.t <- sapply(1:length(origin.t), function(x) tolower(origin.t[x]))
origin.t.tab <- as.data.frame(table(origin.t))
colnames(origin.t.tab) <- c("word", "freq")
hist.max <- max(gloria.t.tab$freq, rue.t.tab$freq, origin.t.tab$freq)
rue.h <- hist(rue.t.tab$freq, breaks = 1:(hist.max + 1))
density.mat <- cbind(rue.h$density, gloria.h$density, origin.h$density)
colnames(density.mat) <- c("Rue", "Gloria", "Origin")
KLmat <- KLdiv(density.mat)
print(KLmat)## Rue Gloria Origin
## Rue 0.00000000 0.06300841 0.1110916
## Gloria 0.07227568 0.00000000 0.1101920
## Origin 0.17588277 0.13931440 0.0000000g <- ggplot()
h3 <- geom_histogram(data = origin.t.tab, aes(x = freq,..density.., fill="Origin of Species"), alpha=0.2, binwidth=1)
h2 <- geom_histogram(data = gloria.t.tab, aes(x = freq,..density.., fill="Gloria Scott"), alpha=0.2, binwidth=1)
h1 <- geom_histogram(data = rue.t.tab, aes(x = freq,..density.., fill="Rue Morgue"), alpha=0.2, binwidth=1)
s1 <- scale_fill_manual("Text", breaks = c("Rue Morgue", "Gloria Scott", "Origin of Species"), values = c("red", "blue", "green"))
print(g + s1 + h1 + h2 + h3)