MathPark
2023-11-18
Génération de graphes avec le modèle d’Erdös-Rényi
Choisir un nombre de noeuds :
n <- 5
n## [1] 5Choisir une probabilité de connexion :
p <- 0.5
p## [1] 0.5Créer une matrice d’adjacence “vide” :
A <- matrix(0, n, n)
A## [,1] [,2] [,3] [,4] [,5]
## [1,] 0 0 0 0 0
## [2,] 0 0 0 0 0
## [3,] 0 0 0 0 0
## [4,] 0 0 0 0 0
## [5,] 0 0 0 0 0Nombre de toutes les paires de noeuds :
N <- n*(n-1)/2
N## [1] 10Tier N réalisations d’une loi de Bernoulli de paramètre
p (= jouer N fois à pile/face avec probabilité
p) :
connexions <- rbinom(N, 1, p)
connexions## [1] 0 0 0 0 1 0 0 1 0 0Remplir la matrice d’adjacence :
A[lower.tri(A)] <- connexions
A## [,1] [,2] [,3] [,4] [,5]
## [1,] 0 0 0 0 0
## [2,] 0 0 0 0 0
## [3,] 0 1 0 0 0
## [4,] 0 0 1 0 0
## [5,] 0 0 0 0 0… et la symétriser :
A <- A + t(A)
A## [,1] [,2] [,3] [,4] [,5]
## [1,] 0 0 0 0 0
## [2,] 0 0 1 0 0
## [3,] 0 1 0 1 0
## [4,] 0 0 1 0 0
## [5,] 0 0 0 0 0rgraphER <- function(n, p){
A <- matrix(0, n, n) # matrice d'adjacence à remplir
N <- n*(n-1)/2 # nb de paires de noeuds
connexions <- rbinom(N, 1, p) # jouer N fois à pile/face avec proba p
A[lower.tri(A)] <- connexions # remplir la matrice d'adjacence
A <- A + t(A) # symétriser la matrice
return(A)
}rgraphER(10, .9)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] 0 1 1 0 1 1 0 1 1 1
## [2,] 1 0 0 1 1 1 1 1 1 1
## [3,] 1 0 0 1 0 1 1 0 1 1
## [4,] 0 1 1 0 1 1 1 1 1 1
## [5,] 1 1 0 1 0 1 1 1 1 1
## [6,] 1 1 1 1 1 0 1 1 1 0
## [7,] 0 1 1 1 1 1 0 1 1 1
## [8,] 1 1 0 1 1 1 1 0 1 1
## [9,] 1 1 1 1 1 1 1 1 0 0
## [10,] 1 1 1 1 1 0 1 1 0 0Réseau social
Charger un jeu de données contenant les connexions d’élèves d’une école primaire durant 2 jours :
load('PrimarySchoolNetwork.RData')La matrice d’adjacence s’appelle primary`
Taille de la matrice d’adjacence :
dim(primary)## [1] 242 242Visualiser le réseau :
library(igraph)##
## Attaching package: 'igraph'## The following objects are masked from 'package:stats':
##
## decompose, spectrum## The following object is masked from 'package:base':
##
## unionprimaryIgraph <- graph_from_adjacency_matrix(primary, mode="undirected")
plot(primaryIgraph)
Degrés de noeuds
Le degré du noeud \(i\) est le nombre de ses connexions :
\[ d_i = \sum_{j=1}^n A_{i,j} \]
degresPrimary <- rowSums(primary)
degresPrimary## [1] 83 47 82 64 37 71 94 39 95 76 32 96 53 77 57 111 98 63
## [19] 48 61 84 53 46 98 44 95 81 49 65 76 83 59 106 82 76 73
## [37] 93 33 39 56 41 72 40 71 23 47 49 92 96 26 46 55 59 79
## [55] 70 134 123 68 101 100 58 97 78 86 51 106 63 71 122 44 79 70
## [73] 97 39 32 56 54 20 111 38 72 34 104 111 63 47 35 58 53 36
## [91] 70 54 54 46 82 33 55 85 123 110 49 57 124 75 104 50 43 40
## [109] 103 63 43 86 103 79 116 99 121 58 95 68 33 68 118 59 46 86
## [127] 56 32 82 22 58 90 84 95 68 80 108 58 73 38 47 52 35 60
## [145] 72 49 56 51 77 92 21 72 78 38 112 79 30 128 91 77 118 68
## [163] 85 96 31 65 53 70 62 44 129 70 72 80 69 90 70 89 66 64
## [181] 75 33 36 74 72 24 29 29 65 54 34 49 54 83 113 36 39 117
## [199] 78 40 96 90 92 89 111 61 47 89 88 51 88 80 30 29 97 35
## [217] 35 45 35 81 89 33 87 120 30 81 30 71 98 40 87 62 97 21
## [235] 97 111 53 109 34 44 71 64plot(table(degresPrimary), xlab='degré', ylab='fréquence', main='Primary school dataset')
Densité d’arêtes
Densité d’arêtes d’un graphe :
\(\alpha=\frac2{n(n-1)}\sum_{i<j}A_{i,j}\)
Densité d’un graphe d’Erdös-Rényi
Simuler un graphe :
A <- rgraphER(500, 0.5)Calculer sa densité :
mean(A[lower.tri(A)])## [1] 0.5000641Dans un graphe Erdös-Rényi \(\alpha\approx\) p.
Plus n est grand, plus \(\alpha\) est proche de p.
D’après la loi des grands nombres, on a même :
\[ \alpha_n \stackrel{p.s.}{\longrightarrow} p \quad \text{ lorsque } n\to\infty. \]
Densité du réseau social
Densité du graphe social :
densPrim <- mean(primary[lower.tri(primary)])
densPrim## [1] 0.2852097Simuler un graphe d’Erdös-Rényi de même densité :
primER <- rgraphER(242, densPrim)Calculons les degrés \(d_i\) de noeuds pour ce graphe :
degresPrimER <- rowSums(primER)
plot(table(degresPrimER), xlab='degré', ylab='fréquence', main='Erdös-Rényi graphe')
Comparons-les au réseau social :
v1 <- tabulate(degresPrimary)
v2 <- tabulate(degresPrimER)
m <- max(c(length(v1), length(v2)))
deg <- matrix(0, 2, m)
deg[1, 1:length(v1)] <- v1
deg[2, 1:length(v2)] <- v2
barplot(deg, beside = TRUE, col = c('blue', 'orange'), border=c('blue', 'orange'), main = 'Distribution des degrés des noeuds')
legend(300, 14, legend=c('Primary School', 'Erdös-Rényi'), col = c('blue', 'orange'), pch=16)
Conclusion : il est très peu probable que ce réseau social a été généré par un modèle Erdös-Rényi.
Nombre de triangles
Pour chaque noeuds, nombre de triangles auxquel appartient le noeud :
count_triangles(primaryIgraph)## [1] 1804 620 1794 1046 456 1267 1953 540 2125 1534 421 2154 688 1517 857
## [16] 2763 2019 992 541 1010 1625 830 573 2037 512 2018 1769 680 1140 1413
## [31] 1785 777 2441 1550 1498 1409 1865 299 541 713 432 1307 366 1225 234
## [46] 549 600 2041 2188 260 605 885 858 1409 1220 3528 3132 1156 2347 2071
## [61] 898 2058 1719 1734 818 2424 952 1287 3091 614 1567 1349 2145 372 313
## [76] 766 593 145 2479 377 1251 287 2380 2618 763 495 386 905 823 364
## [91] 1305 642 662 469 1550 398 724 1727 3073 2597 596 709 3168 1328 2350
## [106] 629 501 496 2366 996 519 1611 2294 1477 2917 2255 3120 879 2280 1240
## [121] 369 1169 3065 934 652 1547 857 337 1704 229 911 1933 1969 2232 1317
## [136] 1495 2661 992 1393 464 793 823 381 930 1235 671 754 747 1547 1933
## [151] 210 1305 1425 360 2857 1302 267 3359 1897 1426 2866 1184 1803 2154 313
## [166] 1018 668 1140 1018 530 3633 1404 1072 1401 1111 1975 1196 2099 935 1061
## [181] 1101 352 451 1472 1152 215 260 254 1074 795 320 517 820 1432 3028
## [196] 432 389 3020 1319 412 2101 1978 2022 1650 2713 752 545 2054 1679 755
## [211] 2008 1518 306 259 1911 410 394 461 410 1386 1866 301 1784 3108 248
## [226] 1334 348 1227 2162 421 1728 909 2327 153 2242 2774 563 2703 319 504
## [241] 1262 997Nombre moyen de triangles auxquels appartient un noeud du graphe :
mean(count_triangles(primaryIgraph))## [1] 1286.281Dans un graphe d’Erdös-Rényi :
primER <- rgraphER(242, densPrim)
primERIgraph <- graph_from_adjacency_matrix(primER, mode = "undirected")
mean(count_triangles(primERIgraph))## [1] 673.7231Spectral clustering
Exemple de graphe avec deux composantes connexes :
A <- matrix(0, 10,10)
A[matrix(c(1,3, 1,7, 2,6, 1,9, 2,10, 8,10, 2,8, 3,5, 3,7, 4,6, 4,8, 5,9, 6,8, 7,9), ncol=2, byrow=T)] <- 1
A <- A + t(A)
G <- graph_from_adjacency_matrix(A, "undirected")
plot(G)
D <- diag(rowSums(A))
L <- D - A
L## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] 3 0 -1 0 0 0 -1 0 -1 0
## [2,] 0 3 0 0 0 -1 0 -1 0 -1
## [3,] -1 0 3 0 -1 0 -1 0 0 0
## [4,] 0 0 0 2 0 -1 0 -1 0 0
## [5,] 0 0 -1 0 2 0 0 0 -1 0
## [6,] 0 -1 0 -1 0 3 0 -1 0 0
## [7,] -1 0 -1 0 0 0 3 0 -1 0
## [8,] 0 -1 0 -1 0 -1 0 4 0 -1
## [9,] -1 0 0 0 -1 0 -1 0 3 0
## [10,] 0 -1 0 0 0 0 0 -1 0 2spectre <- eigen(L)
spectre$values## [1] 5.000000e+00 5.000000e+00 4.414214e+00 4.000000e+00 3.000000e+00
## [6] 3.000000e+00 2.000000e+00 1.585786e+00 7.993606e-15 6.217249e-15plot(spectre$values)
spectre$vectors[ , 9:10]## [,1] [,2]
## [1,] 0.000000e+00 4.472136e-01
## [2,] -4.472136e-01 0.000000e+00
## [3,] 1.732422e-16 4.472136e-01
## [4,] -4.472136e-01 0.000000e+00
## [5,] -4.397343e-17 4.472136e-01
## [6,] -4.472136e-01 2.681443e-17
## [7,] 6.221988e-17 4.472136e-01
## [8,] -4.472136e-01 2.681443e-17
## [9,] -1.043136e-16 4.472136e-01
## [10,] -4.472136e-01 2.681443e-17round(spectre$vectors[ , 9:10], digits=2)## [,1] [,2]
## [1,] 0.00 0.45
## [2,] -0.45 0.00
## [3,] 0.00 0.45
## [4,] -0.45 0.00
## [5,] 0.00 0.45
## [6,] -0.45 0.00
## [7,] 0.00 0.45
## [8,] -0.45 0.00
## [9,] 0.00 0.45
## [10,] -0.45 0.00Exemple de graphe avec deux communautés, mais une composante connexe seulement :
A <- matrix(0, 10, 10)
A[matrix(c(3,6, 1,3, 1,7, 2,6, 1,9, 2,10, 8,10, 2,8, 3,5, 3,7, 4,6, 4,8, 5,9, 6,8, 7,9), ncol=2, byrow=TRUE)] <- 1
A <- A + t(A)
A## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] 0 0 1 0 0 0 1 0 1 0
## [2,] 0 0 0 0 0 1 0 1 0 1
## [3,] 1 0 0 0 1 1 1 0 0 0
## [4,] 0 0 0 0 0 1 0 1 0 0
## [5,] 0 0 1 0 0 0 0 0 1 0
## [6,] 0 1 1 1 0 0 0 1 0 0
## [7,] 1 0 1 0 0 0 0 0 1 0
## [8,] 0 1 0 1 0 1 0 0 0 1
## [9,] 1 0 0 0 1 0 1 0 0 0
## [10,] 0 1 0 0 0 0 0 1 0 0G <- graph_from_adjacency_matrix(A, "undirected")
plot(G)
D <- diag(rowSums(A))
L <- D - A
spectre <- eigen(L)
spectre$values## [1] 5.927380e+00 5.000000e+00 4.686982e+00 4.000000e+00 3.487430e+00
## [6] 3.000000e+00 2.000000e+00 1.632697e+00 2.655107e-01 2.664535e-15plot(spectre$values)
round(spectre$vectors[ , 9:10], digits=2)## [,1] [,2]
## [1,] -0.33 0.32
## [2,] 0.33 0.32
## [3,] -0.21 0.32
## [4,] 0.31 0.32
## [5,] -0.33 0.32
## [6,] 0.20 0.32
## [7,] -0.33 0.32
## [8,] 0.33 0.32
## [9,] -0.36 0.32
## [10,] 0.38 0.32clust <- kmeans(spectre$vectors[ , 9:10], centers = 2)$cluster
clust## [1] 1 2 1 2 1 2 1 2 1 2Dprim <- diag(rowSums(primary))
Lprim <- Dprim - primary
spectrePrim <- eigen(Lprim)
plot(spectrePrim$values)
k <- 10
clust <- kmeans(spectrePrim$vectors[ , (242-k+1):242], centers = k, nstart=500)$cluster
clust## [1] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
## [26] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 6 3 4 3 3 3 3 7
## [51] 3 3 3 3 5 5 5 5 5 5 5 3 5 5 5 5 5 3 5 5 3 3 5 6 5
## [76] 5 5 2 5 5 5 3 5 5 5 5 3 5 5 5 5 5 5 3 5 5 5 5 5 5
## [101] 3 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 3 10 5 5 5 5
## [126] 5 5 10 5 8 5 5 5 5 5 5 5 5 5 5 5 5 10 5 5 5 5 5 5 5
## [151] 1 3 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
## [176] 5 5 5 5 3 5 5 5 5 5 6 6 6 5 3 6 5 5 5 5 3 6 5 5 6
## [201] 3 5 5 5 5 5 5 5 5 5 5 5 10 6 5 10 10 6 10 5 5 6 5 5 6
## [226] 5 10 5 5 5 5 5 5 9 5 5 5 5 5 5 5 5table(clust)## clust
## 1 2 3 4 5 6 7 8 9 10
## 1 1 65 1 151 12 1 1 1 8Laplacien normalisé
La matrice laplacienne normalisée est définie par
\(L = I - D^{-1/2} A D^{-1/2}\)
où \(I\) est la matrice d’identité et \(D^{-1/2}\) la matrice diagonale dont les entrées diagonales valent \(1/\sqrt{d_i}\).
Dm12 <- diag(1 / sqrt(diag(Dprim)))
LN <- diag(242) - Dm12 %*% primary %*% Dm12
spectreLN <- eigen(LN)
plot(spectreLN$values)
k <- 10
clust <- kmeans(spectreLN$vectors[ , (242-k+1):242], centers = k, nstart=100)$cluster
table(clust)## clust
## 1 2 3 4 5 6 7 8 9 10
## 26 24 22 26 22 24 27 24 21 26info <- read.table("metadata_primaryschool.txt")
table(info$V2, clust)## clust
## 1 2 3 4 5 6 7 8 9 10
## 1A 0 23 0 0 0 0 0 0 0 0
## 1B 25 0 0 0 0 0 0 0 0 0
## 2A 0 0 0 0 0 0 0 23 0 0
## 2B 0 0 0 0 0 0 26 0 0 0
## 3A 0 0 0 16 0 0 0 0 7 0
## 3B 0 0 0 10 0 0 0 0 12 0
## 4A 0 0 21 0 0 0 0 0 0 0
## 4B 0 0 0 0 0 23 0 0 0 0
## 5A 0 0 0 0 21 0 0 0 0 1
## 5B 0 0 0 0 0 0 0 0 0 24
## Teachers 1 1 1 0 1 1 1 1 2 1