鬼滅の刃キャラクタ顔形状の形態測定学的分析およびGMMでのクラスタリング

##序論
###背景
　日本における漫画は，テーマや読者層の多様さにおいても，また発行部数においても世界で群を抜いており，今や世界中の言語に翻訳され輸出されています．また，近年では漫画 を媒介とした国際交流が図られるなど，日本を代表するポップカルチャーとして広く浸透していて，その中でも大きな社会現象を起こした作品が『鬼滅の刃』です．このような人気漫画がどのような特徴を持つか明らかにすることで，漫画文化の発展に役立つのではないかと考えました．
###目的
　本研究では近年人気になった『鬼滅の刃』の特徴を明らかにし，日本の漫画の中でどこに位置付けられるか調べるために顔の作画に着目しました．なぜなら，漫画顔形状には漫画に作者の受けた影響や，出版社の傾向など，様々な多様性が存在するからです．
###概要
　日本の様々な漫画の形状を取得し，形態測定学的手法を用いて規格化を行い．多変量データに変換した顔形状をクラスタリングします．クラスタリングには様々な手法が存在しますが，本研究では混合ガウスモデル
(Gaussian Mixture Model：GMM)とBICでのモデル選択を行います．

##解析対象
###漫画材料
　本研究ではNHK『BSマンガ夜話』で題材とされた作品，文化庁「日本のメディア芸術100選マンガ部門」に選ばれた作品，文化庁「メディア芸術祭」，各出版社の漫画賞受賞作品を中心に119作品を材料とし，各マンガ作品に登場するキャラクタの正面向き真顔をスキャンにより取得しました．
　材料となるマンガ作品の中には数十年に渡って連載されているものもあるため，執筆されている時期によってその特徴なども変わっている可能性が考えられます．本研究ではそれらの影響を排除するため，各マンガ作品の1巻を対象としました.また短編集のような形式の作品は取得できる顔画像が少なくなってしまうため，表題作以外の掲載作品も顔画像の取得対象としました．本研究では鬼滅の刃に登場する主要キャラクターを題材としているため，鬼滅の刃のみ全巻を対象にしています．


##顔形状の計測と規格化
###材料
　解析対象としたキャラクタは各作品の主役とそれに準ずる役割のキャラクタです.正面真顔が見当たらないキャラクタについては別々のコマから顔の正面部位を集め，瞳孔や，耳元，鼻，アゴの位置などを基に真顔の正面画像を目視によりphoto shopで合成しました．また，取得できる顔画像が男女それぞれ2名に達しない場合はその他のキャラクタも顔画像の取得の対象とすることで，男性キャラクタ231人，女性キャラクタ186人，計417人のキャラクタの顔画像を取得しました.

###座標点の取得
　取得したマンガ顔画像を座標データとして扱うため，瞳孔や顔，目，眉，鼻，口の輪郭などが含まれる特徴点を自作アプリケーションにより取得しました.
本研究の特徴点の設定は，合成顔を作成するツールである顔画像合成システム「FUTON システム[1]」の特徴点の定義を一部改変して用いた.FUTONシステムでは顔の特徴点を以下のように定義しています.
　目の特徴点は，瞳の中心点，目の内側端点と外側端点，両端点を結ぶ直線を4等分した垂線と目との境界線の交点.対応点は左目9点，右目9点で計18点.
　眉の特徴点は，眉の内側端点と外側端点，目の中心点から両端点へ伸ばした直線と，その角を4等分するよう目の中心点から伸ばした直線と眉の上側/下側の境界との交点.対応点は左眉8点，右眉8点の計16点.
　鼻の特徴点は，鼻の右側上点，右側下点，右側中点，左側上点，左側下点，左側中点，下側端点，頂点.口の中心から，両目の中心を結んだ直線へ降ろした垂線の交点を鼻の上側端点.対応点は計9点となる.口の特徴点は，口の左側端，右側端.両端点を結ぶ直線を3等分しそれぞれ垂線を引く，この直線と口との境界線の交点.対応点は計14点.
　輪郭の特徴点は，あごの先端の点.あごの先端と鼻の上端を結んだ直線を1/4，1/2，5/8，3/4，13/16，7/8に分割する垂線と輪郭の交点12点，さらに，7/8 に分割した点と輪郭の交点の2点，この線分を1/8，1/4，3/8，5/8，3/4，7/8に分割する垂線と顔輪郭の交点6点.あごの先端の点と鼻の上端点を結ぶ線分を反転させ，鼻の上端から輪郭の上部へ伸ばしたものを，左右それぞれに30◦，60◦，90◦ 回転させた直線と輪郭の交点 7 点.対応点は計 26 点.
　マンガの顔表現は幅広く，本研究で用いた顔画像においても鼻や唇の描画されている部分が作品によってまちまちなので，本研究では鼻と唇の特徴点を減らすことにしました.鼻の特徴点は，自動的に計算された鼻の上端点と鼻の頂点の2点としました.また，口の特徴点は，口の両端点を結ぶ直線を3等分しそれぞれ垂線を引き，この直線と口との境界線上唇と下唇の境界3点ずつ，計5点.取得した特徴点は左右の目を9点ずつ，左右の眉を8点ずつ，口を5点，鼻を2点，輪郭を26点，計67点です.
　取得した顔座標は大きさや向きがさまざまであったため，Generalized Procrustes Analysis[2]を用いて標準化しました.これにより，顔座標データは134次元(67点×2($x,y$座標))の多次元正規分布するデータとして扱うことができます.

####GPA
　幾何学的形態測定学の特徴は特徴点を多変量正規分布させることができ，形状を直接統計的に扱うことが可能となります.幾何学的形態測定学を用いることで，形状そのものを変更することなく，複数の形状の規格化をすることが出来ます.以下に形状の規格化の手順を説明します.形状は形態的もしくは機能的な特徴点の座標値として表現され, n個の2次元平面上の特徴点であれば, $2\times n$次元のデータになります.形態上の特徴点を$Z$とし，各形状を重心でセンタリングを行います.重心サイズ (S) は重心 $Z_{1},Z_{2},...,Z_{n}$ から各特徴点間の平方距離の和の平方根として定義されます.
$$
S=\sqrt{n\times\sum_{i-1}^{n}||(Z_i-\bar{Z})||^2}
$$
重心サイズを同一とすることで，サイズの規格化を行います.サイズの違いによる特徴点の分布の違いを減らし，さらに，サイズと重心を規格化されたサンプル間で, 特徴点間の2乗距離の和を最小にするように最適化計算 Generalized Procrustes Analysisを行うことにより，形状を回転し整列を行い，形状の規格化を行う.形状の回転による整列化の手順は次のようになります.2 次元特徴点の場合，複素平面上の$n$個の特徴点$Z_{1},Z_{2},...,Z_{n}$がつくる配置を図形と呼び，ある図形に対して，その重心$Z$へセンタリングを行ったものを前形態と呼ばれます.
$$
\bar{Z}=\frac{1}{n}\times\sum_{i=1}^{n}Z_i
$$
前形態に対してさらに重心サイズSによるスケーリングを行ったものを前形状と呼ばれます．前形状への変換後の座標は$\tau_{i}(i = 1, 2, ..., n)$と表されます.
$$
\tau_i=\frac{Z_i-\bar{Z}}{S}
$$
この定義式から，前形状は，
$$
\sum_{i=1}^{n}\tau_i=0
$$
$$
\sum_{i=1}^{n}||\tau_i||^2=1
$$
を満たします.式は，前形状がある超平面の上にあることを，そして式はそれが半径1の単位超球に含まれることを意味しています.よって，超平面と単位超球が交わる超円の円周が前形状の存在する空間となります.前形状が含まれるこの超円を前形状空間といいます.しかし，前形状空間は回転移動を考慮していないので，前形状 $\tau$ を角 $\theta(0 ≦\theta< 2|\pi)$ だけ回転させて得られるすべての前形状 $\theta(\tau)$の集合 $O(\tau)$を考えます.この集合$O(\tau)$を軌道と呼ばれます.ある同じ軌道に含まれる2つの異なる前形状は，適当な回転移動により必ず一致するので，同値類としてのこの軌道はこの軌道の集合をKendall形状空間と呼ばれます.Kendall形状空間内の軌道，すなわち形状間の距離を定義することが出来れば，形状間の差を定量化できます.回転移動に関して不変な形状とみなせます.ここでの形状とは，複素射影空間の中の多次元座標です.前形状空間は球であるから，球上の2点τ~1~,τ~2~間の最短距離は，2点間の球面上の測地線距離によって表される.この 2 点が半径 1 の球の中心となす角$\phi$とすると$\phi$と測地線距離 $d(\tau_1 , \tau_2)$ との関係は，弧度法の定義により $\phi = d(\tau_{1},\tau_{2})/1 = d(\tau_{1},\tau_{2})$となる.よって，角$\phi(\tau_{1},\tau_{2})$が分かれば$d(\tau_{1},\tau_{2}) = cos−1(||\sum_{k=1}^{n} \tau_{1k} \times\tau_{2k}^{∗} ||)$ が求まります.2つの前形状$\tau_{1},\tau_{2}$ それぞれの軌道 $o_{1},o_{2}$ 内での回転を考えたとき，測地線距離$d(\tau_{1},\tau_{2})$が最小になるように，回転角 $\theta_{1},\theta_{2}$ を決める.このときの最小値 $d(\tau_{1},\tau_{2})$は，
$$
d(o_1,o_2)=cos^{-1}(||\sum_{k=1}^{n}\tau_{1k}\times\tau_{2k}^{*}||)
$$
と計算できます.上の式で与えられる形状空間距離$d(o_{1},o_{2})$を，プロクラステス距離と呼びます.プロクラステス距離はKendall形状空間の表面に分布する点の距離を求めるので非ユーグリッド距離となります.これを，接平面に射影すると近似的に線形として扱うことが可能となります.これらにより特徴点が線形空間で多変量正規分布します.

###主成分分析
　主成分分析(Principle Component Analysis)は多次元のデータを次元圧縮する方法です．
　本研究では整列された特徴点座標データは134次元あり，そのままでは解釈が困難なため，主成分分析を行い次元を圧縮しました．累積寄与率が72.88%となる第7主成分までを対象にしました．

###年代の規格化
####多次元尺度構成法
　多次元尺度構成（Multi-Dimensional Scaling: MDS）は、対象間の類似度を視覚化する多変量解析の一手法です．主成分分析の様に分類対象物の関係を低次元空間における点の布置で表現する手法です.
行列$X$を求めるためには、まず

```math
K=X^{T}X+\mathbb{R}^{n\times n}
```
を満たす行列$K$を、入力として与えられた距離行列＄D$から計算できるとします.この行列$K$は固有値分解によって

```math
K=V\varLambda V^T
```
と分解することができます.ここで$V$は固有ベクトルを並べた行列,$Λ$は固有値を対角成分に持つ行列です.これらを用いて

```math
X=\sqrt{(\varLambda)V^T}
```
とすることで、所望の$X$を計算することができます.

###混合ガウスモデル
####式説明
　混合ガウスモデル(Gaussian Mixture Model: GMM)は，正規分布を要素とした混合分布モデルです．正規分布のような単純な形をしたものでも混合モデルとしては十分な多様性を持ち，正規分布は平均と分散のみで記述できるので，入力の次元が高くなってもパラメータ数は入力次元の2乗程度であり，高次元データ処理に適しています．
入力$x∈\Re^{p}$とし，平均$\mu$，分散$\sigma$の$p$次元正規分布の密度関数を
$$
\phi(x;\mu,\sigma)=\frac{1}{\sqrt(2\pi)^p|\sigma|}exp(-\frac{1}{2}(x-\mu)^r\sigma^{-1}(x-\mu))
$$

とするとき，M個の正規分布を要素モデルとする混合正規分布の確率密度関数は

$$
p(x;θ)=\sum_{m=1}^M\pi_m\phi(x;\mu_m,\sigma_m)
$$

で表されます．ただし，パラメータ$\pi_{m}$は

$$
\pi_m≥0(m=1,…,M),\sum_{m=1}^{M}π_m =1
$$

を満たし，$\theta$はすべてのパラメータをまとめて表したもので

$$
θ=(π_m,μ_m,\sigma_m;m=1,…,M)
$$

となります．
　混合パラメータ$\pi_{m}$に関する条件から，$\pi_{m}$は要素モデル$ϕ(x;\mu_{m},\sigma_{m})$がデータ生成するために選ばれる確率を表しているとみなせます．要素モデル$m$が選ばれ，そのモデルからデータ$x$が生成される同時分布の密度関数は
 
$$
p(x,m;\theta)=\pi_m\phi(x;\mu_m,\sigma_m )
$$

で与えられます．このため混合分布は，どの要素モデルから生成されたかという情報に関して上の同時分布を周辺化したものと考えられます．このときベイズの定理を用いれば，データ$x$が観測された時，それが要素モデル$m$から生成された事後確率は

$$
p(m|x;\theta)=\frac{p(x,m;\theta)}{p(x;\theta)}=\frac{(\pi_m \phi(x;\mu_m,\sigma_m ))}{(\sum_{m=1}^{M}\pi_m \phi(x;\mu_m,\sigma_m))}
$$

によって計算されます．この事後分布を用いてデータが生成されたと考えられる要素モデルを決定し，同じ要素モデルから生成されたと考えられるデータを1つのクラスタにまとめることによってクラスタリングを行っています．
n個のデータ${x_{1},…,x_{n}$が観測されたとき，混合正規分布モデルの最尤推定量は以下の式で定義されます．

```math
\begin{align}
\hat{\theta}&=\textrm{argmax}\theta\sum_{i=1}^{n}\log p(x_i;\theta)
\\ &=\textrm{argmax}_\theta\sum_{i=1}^{n}\log\sum_{m=1}^{M}\pi_m\phi(x_i\mu_m,\sigma_m)
\end{align}
```

単純そうな式に見えるが，対数のなかに和があるため，推定量を陽に記述することができません．データ $x_i$ を生成した要素モデルが$m_i$であるとわかっていて， ${(x_{1},m_{1}),…,(x_{n},m-{n})}$ が与えられているのであれば，最尤推定量は，

```math
\begin{align}
\hat{\theta}&=\textrm{argmax}_\theta\sum_{i=1}^{n}\log p(x_i,m_i;\theta)
\\
&=\textrm{argmax}_\theta\sum_{i=1}^{n}\log\pi_{m_i}\phi(x_i,\mu_m;\sigma_m)
\end{align}
```

で与えられます．このとき$x_{i}$がどの要素モデルから生成されたか示す変数として

```math
\mu_{i}^{m} = \left\{
\begin{array}{ll}
1,m=m_i \\
 0,m\ne m_i
\end{array}
\right.
```


を導入すると最尤推定量は

$$
\hat{\pi}=\frac{\sum_{i=1}^{n}\mu_{i}^{m}}{n}
$$

```math
\hat{\mu}_m=\frac{\sum_{i=1}^{n}\mu_{i}^{m}x_i}{\sum_{i=1}^{n}\mu_{i}^{m}}
```

```math
\hat{\sigma}_m=\frac{\sum_{i=1}^{n}\mu_{i}^{m}x_ix_{i}^{T}}{\sum_{i=1}^{n}\mu_{i}^{m}}-\hat{\mu}_m\hat{\mu}_{m}^{T}
```

と陽に記述することができます．データを生成した要素モデル$m_{i}$が欠損しているが故に，混合分布の最尤推定の難しさが発生しています．
混合正規分布の混合数の推定には要素モデルの数をいくつに設定するかという問題が極めて重要です．モデル選択の手法には交差検証法，情報量基準を用いるものがありますが，本研究では情報量基準BIC(Bayesian Information Criterion)を使用しています．

####mclustモデル説明
　EIIは分散が単位行列の定数倍で，すべての要素モデルで共通な制約．
　VIIは分散が単位行列の定数倍だが，すべての要素モデルで異なる制約．
　EEIは分散が対角だが，すべての要素モデルで全て共通な制約．
　VEIは分散が対角だが，すべての要素モデルで体積は異なるが，形は共通な制約．
　EVIは分散が対角だが，すべての要素モデルで体積は同じだが，形は異なる制約．
　VVIは分散が対角だが，すべての要素モデルで全てが異なる制約．
　EEEは分散に制約なく，すべての要素モデルで全て共通な制約．
　EEVは分散に制約なく，すべての要素モデルで体積，形は同じだが，向きは異なる制約．
　VEVは分散に制約なく，すべての要素モデルで形は同じ制約．
　VVVは分散に制約なく，すべての要素モデルで体積，形，向きが異なる制約．
　本研究では比較的自由度の高いモデルのVVVを用いました．

###BIC
　ベイズ情報量規準は，統計学における情報量規準の一つであり，ベイズ定理を用いて複数のモデルの中で，あるデータが生起する事後確率に対応する周辺尤度を最大にするモデルを 選択する理論に基づいて提案されたモデルの評価基準です.最小記述長(Minimum Description Length: MDL)と呼ばれる情報量基準と形式的には同じものです.本研究ではBICを用いてモデル選択を行った．$X$を真の分布にしたがう未知データとした場合のBICの式を以下に示す．
$$
\textrm{BIC}=-2nE \biggl(\log\int p(X;\theta_l)p(\theta_l|D)d\theta_l \biggl)
$$
となります．ここでp(θ|D)はパラメータ事後確率で与えられたデータから推測されるパラメータの分布を表します.BICはモデル全体での平均的なあてはまりの良さを評価しています． BICの値は6607.3でした． 式の右辺第2項は最尤推定量に対する対数尤度を補正した量．負の対数尤度である第1項は自由度が高くなるほど小さくなるが，第2項はモデルの自由度が高くなるに従って大きくなる罰金項です．一般に，モデルを複雑にすればするほどデータによくフィットするようになる反面，未知のデータに対する汎化性能が悪化するという問題があります(過学習).モデルの複雑さに対する罰則項は過学習を起こすことがないようにモデルに制約をかけます．

##使用した環境
macOS 10.14.6 Mojave
R 3.6.2
mclust 5.4.7
shapes 1.2.5
stats 3.6.1
##ソースコード
###パッケージインストール
とりあえずpackageをインストールします．

```
install.packages("mclust")#GMMを行うパッケージ
install.packages("shapes")#procGPAを行うパッケージ
install.packages("stats")#MDSを行うパッケージ
```

###procGPA
ここではprocGPAを行います．

```
library(shapes)
options(max.print = 2147483647)
fix(datdata.dat)#座標点の入力
datamanga <-procGPA(datdata.dat)#procGPAで計算を行う
shapepca(datamanga,pcno = c(7),type = "v",mag = 3)#pcaの表示
```
正常に入力を行うと，“function”と書かれている新しいウィンドウが開くので，全てを消して何も書かれていない状態にします．
全てなくなったところに以下を書きます．本研究では特徴点数が67,次元数が2,データ数が417です．

```
Structure(c(座標データ), . Dim = as.intger(c(特徴点数, 次元数, データ数)))
```

保存するデータは stdscores です．間違えないように．
###MDS
ここにはMDSを行う際のソースコードを載せておきます．
$k$=2となってますが，これは圧縮後の次元数になっているので必要に応じて増やしたり減らしたりしてください．

```
library(stats)
dimension <- 2
Ydata <- read.csv("年代の対数表",header = F)
YMDS <- cmdscale(Ydata,k=dimension)
```

###顔のユークリッド距離を求める
規格化が正しいものと判断するために，PCAの数値からユークリッド距離を求め，顔形状の近いキャラを求めます．

```
ncluster <- 7  #meanの数
ndim <- 7  #次元数
#各主成分のmeanを行ごとに入れる
mean <- read.csv("平均値",header = T,row.names = 1)
cmean <- as.matrix(mean)　#matrix型に変換する．
is.matrix(cmean)
nchara <- 417  #キャラ数
#各キャラの座標を行ごとに入れる
pca <- read.csv("PCAの値",header = T,row.names = 1)
All<-as.matrix(pca) #matrix型に変換する．
#ユークリッド距離の計算
distmat <- matrix(NA, nrow=nchara, ncol=ncluster)
for(i in 1:nchara){
  for(j in 1:ncluster){
    tmp <- cmean[j,]-All[i,]
    distmat[i,j] <- sqrt(tmp %*% tmp)
  }
}
#distmatに行列の名前をつける
rownames(distmat) <- rownames(All)
colnames(distmat) <-rownames(mean)
```

###ユークリッド距離で求めた近い顔
ここに載っていないキャラクタは正面真顔が無かったり，標識点のパーツが無いキャラクタです．

|鬼滅の刃キャラクタ|作品名	|キャラクタ名|
|:-:|:-:|:-:|
|竈門炭治郎|逃げるは恥だが役にたつ|津崎|
|竈門禰豆子|ONEPIECE|ルフィ|
|富岡義勇|ドラゴンボール|孫悟空|
|吾妻善逸|ドラえもん|野比のび助|
|嘴平伊之助|キン肉マン|マリ|
|栗花落カナヲ|ちはやふる|綿谷新|
|甘露寺蜜璃|櫻の園|妃子|
|時透無一郎|風と木の詩|セルジュの母|
|煉獄杏寿郎|マカロニほうれん荘|和子|
|胡蝶しのぶ|はみだしっこ|エダ|
|胡蝶カナエ|あずまんが大王|ちよちゃん|
|累に刻まれた剣士|タッチ|南|
|童磨|pink|ユミの母|
|堕姫|Dr.スランプ|山吹みどり|
|猗窩座|よつばと!|小岩井|
|鬼舞辻無惨|がきデカ|もも子|
|鬼舞辻無惨(女)|わたしは真悟|さとる|
|鬼舞辻無惨(鬼)|湘南爆走族|江口|
|継国縁壱|大きく振りかぶって|百枝|
|竈門炭十郎|アイランド|ミホ|
|神崎アオイ|おいしい関係|百恵|
|珠世|スラムダンク|春子|
|愈史郎|うしおととら|うしお|
|産屋敷耀哉|サイボーグ009の5期以降|女|

似ているキャラが多く，規格化が正しいと言えます．
漫画顔形状自体を載せるわけには行かないので，それぞれで調べてくださいすいません．

###GMM
ではGMMでクラスタリングを行います．

```
library(mclust)
dat <- read.csv("主成分分析した後のデータ",header = T,row.names = 1)
gmm<- Mclust(dat)#クラスターごとに分ける
gmm$classification#クラスわけを見る
class<-gmm$classification
plot.mclustBIC(gmm$BIC) #BICのグラフを見る
dev.new()
randProj(dat,parameters = gmm$parameters,classification = class)#ランダムな二次元散布図
clPairs(dat,gmm$classification,symbols = as.character(1:5))#推定したクラスの二次元散布図
```

Mclustを使用する際にMclust(modelnames"モデル名")でモデルを選択することができます．$G$はクラス数です．$G$やモデル指定をするとそのモデルやクラス数でクラスタリングされます．指定をしなかった場合は全てのモデルとクラス数の組み合わせのなかからBICが最大値となるものでクラスタリングが行われます．
本研究のBICのグラフは以下です．

![bic.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/683895/5a1dd150-7592-854c-3ea6-12c365a0b74f.png)

VEEの9が最も高い数値となっています．ですが本研究では比較的自由度の高いVVVを用いました．(大切なことなので二回言いました．)

こちらはクラスタリングを行った二次元散布図です．7次元のデータを2次元にして表示しているので見にくいです．

![PLOT.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/683895/148c4f0b-e88b-2e23-e036-5070e03e4591.jpeg)


こちらは各主成分ごとに総当たりで散布図をプロットした場合です．

![souataricluster.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/683895/1e95d1f1-0b50-9396-3cdc-5b176d9dfafc.png)


##結果
###クラスタリングの結果
クラスタリングおよびモデル評価しました．BICで最適なクラス数を推定した結果，$K=4$になりました．
クラス1は少女漫画のヒロインが多く集まるクラス．クラス2は青年漫画に登場する中年の男性キャラクタが多く集まるクラス．クラス3は鳥山明作品や『ONE PIECE』といったジャンプ作品の主人公が集まるクラス．クラス4は少女漫画の美青年なキャラが多く集まるクラスとなりました．一部抜粋した結果を以下に示します．

####GMMによるキャラクタ顔形状のクラスタ分析の結果(K=4)
|クラス	|代表的なキャラクタ|
|:-:|:-:|
|1|鬼舞辻無惨(鬼滅の刃),幽助(幽遊白書),千早(ちはやふる)ナルト(NARUTO),ちよちゃん(あずまんが大王),妃子(櫻の国)ユミ(pink),みくり(逃げるは恥だが役に立つ),達也(タッチ)|
|2|	煉獄(鬼滅の刃),キン肉マン(キン肉マン),エレン(進撃の巨人)正助(ぶっせん),梶原一騎(空手バカ一代),無用ノ介(無用ノ介)宗方(エースをねらえ),エド(鋼の錬金術師),花男(花男)|
|3|竈門炭治郎(鬼滅の刃),禰豆子(鬼滅の刃),悟空(ドラゴンボール)アンジー(はみだしっ子),サザエ(サザエさん),カムイ(カムイ伝)飛鳥(デビルマン),ルフィ(ONEPIECE),パタリロ(パタリロ!),|
|4|廉(ストロボエッジ),冬威(To-y),ジョナサン(JOJOの奇妙な冒険)ヴィエラ(のだめカンタービレ),オスカル(ベルサイユのばら)ブラックジャック(ブラックジャック),ケンシロウ(北斗の拳),夜神月(デスノート),シャンクス(ONEPIECE)|

##考察
　GMMでクラスタリングすることによって鬼滅の刃のキャラクターをいままでの日本の漫画の累計に分類することができた．
　炭治郎，禰荳子といった主人公キャラクタはクラスタ3に分類されていた．ジャンプの鳥山明や『ONE PIECE』の系譜にある主人公キャラの特徴を反映している事がわかる.それに対して主人公ではないカナヲやしのぶといった女性や，鬼舞辻無惨や童磨はクラスタ1に分類されていた．少女漫画の系譜にあることがわかる．煉獄はクラスタ2に分類されていた．青年漫画の系譜にあることがわかる．
竈門炭治郎，竈門禰豆子，吾妻善逸，富岡義勇といった主人公キャラクタはクラス3に分類されていた．主人公級のキャラクタが『ドラゴンボール』といった鳥山明作品や，『ONE PIECE』といったジャンプ作品の系譜にある主人公キャラクタの特徴と，『サザエさん』や『ドラえもん』といった国民的伝統漫画作品の特徴を反映している事がわかる．竈門炭治郎，竈門禰豆子，吾妻善逸は『鬼滅の刃』の話中において大きく成長していくキャラクタで孫悟空やルフィのような，成長を遂げてもらうという願いからこのような顔形状になったと考えられる．
　主人公ではないキャラクタである栗花落カナヲや胡蝶しのぶといった女性や，鬼舞辻無惨や童磨といった鬼はクラス1に分類されていた．『ちはやふる』といった少女漫画や『幽遊白書』といった少女漫画の影響を受けた作品の系譜にあるキャラクタの特徴を反映していることがわかる．主人公キャラが少年漫画で，主人公でないキャラを少女漫画の特徴にすることで，様々な性別に違和感がもたれない顔形状になったと考えられる．
　煉獄杏寿郎はクラス2に分類されていた．青年漫画や，少年漫画，少女漫画に登場する師匠や先生のキャラクタの影響を受けていることがわかる．煉獄杏寿郎は『鬼滅の刃』本編でも主人公である竈門炭治郎の成長に大いに貢献しているキャラクタであるためこのような顔形状にしたのだと考えられる．

　以上のことから鬼滅の刃は『ドラゴンボール』や『ONE PIECE』といった伝統的なジャンプの少年や主人公の漫画の描き方に則りながら，『ちはやふる』等，少女漫画のキャラクタに類似したキャラを出すことで，性別を超えてジャンプの読者層以外の幅広い年代の人に受け入れられたと考えられる．

##最後に
結果には色々ありましたがこのような結果となりました．
ほぼ初めての投稿ですので色々と不親切のところがありましたがご視聴ありがとうございました！

#参考文献
[1]蒲池みゆき, 向田茂, 吉川左紀子, 加藤隆, 尾田政臣, 赤松茂. (1998). 顔・表情認知に関する心理実験のための顔画像合成システム: FUTON System. 電子情報通信学会技術研究報告. HIP, ヒューマン情報処理, 97(509), 73–80.
[2]三中信宏. (1999). 形態測定学. 古生物の形態と解析, 朝倉書店, 61–99.