> library(pacman)
> p_load(dplyr, kohonen, GGally)

SOM（Self Organizing Maps，自組織映射）本質(zhì)上是一種只有輸入層--隱藏層的神經(jīng)網(wǎng)絡(luò)祭芦。輸入層神經(jīng)元的數(shù)量是由輸入數(shù)據(jù)的維度決定的，一個(gè)神經(jīng)元對(duì)應(yīng)一個(gè)特征，隱藏層中的一個(gè)節(jié)點(diǎn)代表一個(gè)需要聚成的類(lèi)蒂窒。訓(xùn)練時(shí)采用“競(jìng)爭(zhēng)學(xué)習(xí)”的方式，每個(gè)輸入的樣例在隱藏層中找到一個(gè)和它最匹配的節(jié)點(diǎn)荞怒，稱(chēng)為它的激活節(jié)點(diǎn)洒琢，也叫“winning neuron”。 緊接著用隨機(jī)梯度下降法更新激活節(jié)點(diǎn)的參數(shù)褐桌。同時(shí)衰抑，和激活節(jié)點(diǎn)臨近的點(diǎn)也根據(jù)它們距離激活節(jié)點(diǎn)的遠(yuǎn)近而適當(dāng)?shù)馗聟?shù)。兩種常見(jiàn)鄰域函數(shù)：bubble function和Gaussian function荧嵌。

1呛踊、SOM步驟總結(jié)

1.創(chuàng)建節(jié)點(diǎn)網(wǎng)格。
2.隨機(jī)給每個(gè)節(jié)點(diǎn)分配權(quán)重（數(shù)據(jù)集中每個(gè)變量一個(gè)權(quán)重（很小的隨機(jī)數(shù)））啦撮。
3.隨機(jī)選擇一行谭网，并計(jì)算其與網(wǎng)格中每個(gè)節(jié)點(diǎn)權(quán)重的距離（相似度，通常為歐式距離）赃春。
4.把此行放到權(quán)重與該行距離最小的節(jié)點(diǎn)中（BMU愉择，best matching unit）。
5.更新BMU（基本思想是：越靠近優(yōu)勝節(jié)點(diǎn)织中，更新幅度越大锥涕；越遠(yuǎn)離優(yōu)勝節(jié)點(diǎn)，更新幅度越邢梁稹）及其鄰域內(nèi)節(jié)點(diǎn)的權(quán)重（取決于鄰域函數(shù)）层坠。
6.重復(fù)步驟3-5，迭代指定次數(shù)刁笙。

2破花、kohonen包最重要的四個(gè)函數(shù)

som()谦趣、xyf()、supersom()旧乞、somgrid()

簡(jiǎn)單說(shuō)蔚润，som()和xyf()是supersom()的封裝版本，分別對(duì)應(yīng)單層SOM和雙層SOM尺栖，如果是兩層以上的多層SOM嫡纠，必須使用supersom()。somgrid()函數(shù)用于建立SOM網(wǎng)絡(luò)延赌。

3除盏、實(shí)例

> data(flea, package = "GGally")
> str(flea)

## 'data.frame':    74 obs. of  7 variables:
##  $ species: Factor w/ 3 levels "Concinna","Heikert.",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ tars1  : int  191 185 200 173 171 160 188 186 174 163 ...
##  $ tars2  : int  131 134 137 127 118 118 134 129 131 115 ...
##  $ head   : int  53 50 52 50 49 47 54 51 52 47 ...
##  $ aede1  : int  150 147 144 144 153 140 151 143 144 142 ...
##  $ aede2  : int  15 13 14 16 13 15 14 14 14 15 ...
##  $ aede3  : int  104 105 102 97 106 99 98 110 116 95 ...

74行7列，1列為因子型挫以，其他全為整數(shù)型者蠕。

> DataExplorer::profile_missing(flea)

##   feature num_missing pct_missing
## 1 species           0           0
## 2   tars1           0           0
## 3   tars2           0           0
## 4    head           0           0
## 5   aede1           0           0
## 6   aede2           0           0
## 7   aede3           0           0

無(wú)缺失值。

> ggpairs(flea, axisLabels = "none", aes(col = species), 
+         upper = list(continuous = ggally_density, 
+                      combo = ggally_box_no_facet)) +
+   theme_bw()

變量相關(guān)性

> # xdim/ydim：網(wǎng)格尺寸
> # topo：六邊形或矩形掐松，Hexagonal或Rectangular
> # neighbourhood.fct：鄰近函數(shù)踱侣，bubble或gaussian
> # toroidal：是否為環(huán)形
> som.grid <- somgrid(xdim = 5, ydim = 5, topo = "hexagonal",
+                     neighbourhood.fct = "bubble", toroidal = F)
> som.grid

## $pts
##         x         y
##  [1,] 1.5 0.8660254
##  [2,] 2.5 0.8660254
##  [3,] 3.5 0.8660254
##  [4,] 4.5 0.8660254
##  [5,] 5.5 0.8660254
##  [6,] 1.0 1.7320508
##  [7,] 2.0 1.7320508
##  [8,] 3.0 1.7320508
##  [9,] 4.0 1.7320508
## [10,] 5.0 1.7320508
## [11,] 1.5 2.5980762
## [12,] 2.5 2.5980762
## [13,] 3.5 2.5980762
## [14,] 4.5 2.5980762
## [15,] 5.5 2.5980762
## [16,] 1.0 3.4641016
## [17,] 2.0 3.4641016
## [18,] 3.0 3.4641016
## [19,] 4.0 3.4641016
## [20,] 5.0 3.4641016
## [21,] 1.5 4.3301270
## [22,] 2.5 4.3301270
## [23,] 3.5 4.3301270
## [24,] 4.5 4.3301270
## [25,] 5.5 4.3301270
## 
## $xdim
## [1] 5
## 
## $ydim
## [1] 5
## 
## $topo
## [1] "hexagonal"
## 
## $neighbourhood.fct
## [1] bubble
## Levels: bubble gaussian
## 
## $toroidal
## [1] FALSE
## 
## attr(,"class")
## [1] "somgrid"

標(biāo)準(zhǔn)化，有助于使每個(gè)特征對(duì)于計(jì)算相似度（距離）的貢獻(xiàn)相同大磺。

> flea.scale <- flea %>% 
+   as_tibble() %>%
+   # 去掉species變量
+   select(-species) %>% 
+   # 標(biāo)準(zhǔn)化
+   scale()
> str(flea.scale)

##  num [1:74, 1:6] 0.467 0.263 0.773 -0.145 -0.213 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : NULL
##   ..$ : chr [1:6] "tars1" "tars2" "head" "aede1" ...
##  - attr(*, "scaled:center")= Named num [1:6] 177.3 124 50.4 134.8 13 ...
##   ..- attr(*, "names")= chr [1:6] "tars1" "tars2" "head" "aede1" ...
##  - attr(*, "scaled:scale")= Named num [1:6] 29.41 8.48 2.75 10.35 2.14 ...
##   ..- attr(*, "names")= chr [1:6] "tars1" "tars2" "head" "aede1" ...

> # 標(biāo)準(zhǔn)化的兩個(gè)屬性
> # 新數(shù)據(jù)也需要使用這個(gè)屬性抡句，否則結(jié)果將不一致
> attr(flea.scale, "scaled:center")

##     tars1     tars2      head     aede1     aede2     aede3 
## 177.25676 123.95946  50.35135 134.81081  12.98649  95.37838

> attr(flea.scale, "scaled:scale")

##     tars1     tars2      head     aede1     aede2     aede3 
## 29.412541  8.481146  2.751998 10.350932  2.142162 14.304614

> # rlen：迭代次數(shù)
> # alpha：學(xué)習(xí)速率，默認(rèn)從0.05下降到0.01
> flea.som <- som(flea.scale, grid = som.grid, rlen = 5000, alpha = c(0.05, 0.01))
> flea.som

## SOM of size 5x5 with a hexagonal topology.
## Training data included.

畫(huà)圖：

> par(mfrow = c(2, 3))
> plot.type <- c("codes", "changes", "counts", "quality", "dist.neighbours", "mapping")
> # 根據(jù)plot.type中的類(lèi)型杠愧，依次畫(huà)圖待榔，并按2×3排列
> purrr::walk(plot.type, ~ plot(flea.som, type = ., shape = "straight"))

依次畫(huà)圖

"changes" - Training progress：展示訓(xùn)練過(guò)程，距離隨著迭代減少的趨勢(shì)流济，判斷迭代是否足夠锐锣，最后趨于平穩(wěn)比較好。
"codes" - Codes plot：查看SOM中心點(diǎn)的變化趨勢(shì)绳瘟。
"counts" - Counts plot：展示每個(gè)SOM中心點(diǎn)包含的樣本數(shù)目雕憔。可以跟“mapping”一起看糖声，“counts”顏色越淺橘茉，對(duì)應(yīng)的“mapping”數(shù)量越多。
"dist.neighbours" - Neighbours distance plot：鄰近距離姨丈，查看潛在邊界點(diǎn)，顏色越深表示與周邊點(diǎn)差別越大擅腰，越可能是邊界點(diǎn)蟋恬。
"mapping" - Mapping plot：展示每個(gè)樣本的映射。
"quality" - Quality plot：計(jì)量SOM中心點(diǎn)的內(nèi)斂性和質(zhì)量趁冈，距離越小展示得越好歼争。

"property"：每個(gè)單元的屬性可以計(jì)算并顯示在顏色代碼中拜马。用來(lái)可視化一個(gè)特定對(duì)象與映射中所有單元的相似性，以顯示所有單元和映射到它們的對(duì)象的平均相似性沐绒。

> getCodes(flea.som) %>% 
+   as_tibble() %>% 
+   # property：屬性值
+   # main：標(biāo)題
+   purrr::iwalk(~ plot(flea.som, type = "property", property = .,
+                main = .y, shape = "straight"))

tars1

針對(duì)每個(gè)特征會(huì)有一幅圖俩莽，就不一一展示了。

4乔遮、SOM中心點(diǎn)相關(guān)的樣本數(shù)量

> table(flea.som$unit.classif)

## 
##  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 23 24 25 
##  3  3  4  2  2  3  3  7  1  3  3  2  2  2  4  3  3  2  4  5  4  2  4  3

因?yàn)槎x的網(wǎng)格為5×5扮超，所以一共25個(gè)中心點(diǎn)。

> # 不能用flea.scale蹋肮，因?yàn)閠ibble沒(méi)有行名
> code.class <- tibble(name = rownames(flea),
+                      class = flea.som$unit.classif)
> head(code.class)

## # A tibble: 6 x 2
##   name  class
##   <chr> <dbl>
## 1 1         3
## 2 2         8
## 3 3         2
## 4 4        13
## 5 5        12
## 6 6        18

5出刷、SOM結(jié)果進(jìn)一步聚類(lèi)

> # 轉(zhuǎn)換為數(shù)據(jù)框再轉(zhuǎn)換為矩陣
> mydata <- as.matrix(as.data.frame(flea.som$codes))
> wss <- (nrow(mydata) - 1) * sum(apply(mydata, 2, var)) 
> for (i in 2:15) wss[i] <- sum(kmeans(mydata, centers = i)$withinss)
> 
> plot(1:15, wss, type = "b", xlab = "聚類(lèi)數(shù)量",
+      ylab = "類(lèi)內(nèi)平方和", main = "類(lèi)內(nèi)平方和 (WCSS)")

選擇聚類(lèi)數(shù)量

選擇曲線(xiàn)逐漸開(kāi)始平緩的第一個(gè)拐點(diǎn)，本例選擇3類(lèi)坯辩。

> # hclust聚類(lèi)后剪枝為3類(lèi)
> som.cluster <- cutree(hclust(dist(mydata)), 3)
> 
> # 定義色塊顏色
> cluster.palette <- function(x, alpha = 0.6) {
+   n = length(unique(x)) * 2
+   rainbow(n, start = 1/3, end = 3/3, alpha = alpha)[seq(n, 0, -2)]
+ }
> cluster.palette.init <- cluster.palette(som.cluster)
> bgcol <- cluster.palette.init[som.cluster]
> 
> plot(flea.som, type="codes", bgcol = bgcol, main = "Clusters", codeRendering="lines")
> add.cluster.boundaries(flea.som, som.cluster)

重新聚類(lèi)

查看數(shù)據(jù)所在新聚的類(lèi)：

> code.class <- code.class %>% 
+   bind_cols(new_class = som.cluster[code.class$class])
> head(code.class)

## # A tibble: 6 x 3
##   name  class new_class
##   <chr> <dbl>     <int>
## 1 1         3         1
## 2 2         8         1
## 3 3         2         1
## 4 4        13         1
## 5 5        12         2
## 6 6        18         2

6馁龟、新數(shù)據(jù)上應(yīng)用SOM

> new.data <- tibble(tars1 = c(120, 200),
+                    tars2 = c(125, 120),
+                    head = c(52, 48),
+                    aede1 = c(140, 128),
+                    aede2 = c(12, 14),
+                    aede3 = c(100, 85))
> 
> new.flea <- new.data %>% 
+   # 使用之前的屬性值標(biāo)準(zhǔn)化
+   scale(center = attr(flea.scale, "scaled:center"),
+         scale = attr(flea.scale, "scaled:scale")) %>% 
+   # 預(yù)測(cè)
+   predict(flea.som, newdata = .)
> 
> plot(flea.som, type = "counts", classif = new.flea, shape = "round")

新數(shù)據(jù)映射

新數(shù)據(jù)所在新聚的類(lèi)：

> # 預(yù)測(cè)的聚類(lèi)
> new.flea$unit.classif

## [1] 11 14

> # 新聚類(lèi)
> som.cluster[new.flea$unit.classif]

## V11 V14 
##   2   3

7、與K-Means的比較

（1）K-Means需要事先定下類(lèi)的個(gè)數(shù)漆魔，也就是K的值坷檩。SOM則不用，隱藏層中的某些節(jié)點(diǎn)可以沒(méi)有任何輸入數(shù)據(jù)屬于它改抡。所以矢炼，K-Means受初始化的影響要比較大。
（2）K-means為每個(gè)輸入數(shù)據(jù)找到一個(gè)最相似的類(lèi)后雀摘，只更新這個(gè)類(lèi)的參數(shù)裸删。SOM則會(huì)更新臨近的節(jié)點(diǎn)。所以K-mean受noise data的影響比較大阵赠，SOM的準(zhǔn)確性可能會(huì)比k-means低（因?yàn)橐哺铝伺R近節(jié)點(diǎn)）涯塔。
（3）SOM的可視化比較好。

8清蚀、練習(xí)

1.設(shè)定 topo=rectangular,toroidal=T匕荸，重新運(yùn)行 SOM 比較。

> som.grid2 <- somgrid(xdim = 5, ydim = 5, topo = "rectangular",
+                      neighbourhood.fct = "bubble", toroidal = T)
> flea.som2 <- som(flea.scale, grid = som.grid2, 
+                  rlen = 5000, alpha = c(0.05, 0.01))
> 
> par(mfrow = c(1, 2))
> plot(flea.som, type = "mapping", shape = "straight", main = "Hexagonal")
> plot(flea.som2, type = "mapping", shape = "straight", main = "Rectangular")

比較

2.利用同一個(gè) somgrid枷邪，迭代次數(shù)改為 10000 次榛搔，alpha 設(shè)為 c(0.1, 0.001) 來(lái)做SOM。

> flea.som3 <- som(flea.scale, grid = som.grid, 
+                  rlen = 10000, alpha = c(0.1, 0.001))
> 
> par(mfrow = c(1, 2))
> plot(flea.som, type = "mapping", shape = "straight", main = "5000次")
> plot(flea.som3, type = "mapping", shape = "straight", main = "10000次")

比較