Elixir の Explorer で主成分分析（PCA）をやってみた

## はじめに

前回の記事では Explorer をとりあえず使ってみました

https://qiita.com/RyoWakabayashi/items/dfb5495b9c25eff710a2

結果として、 Python の pandas のようにデータ分析が簡単にできることが分かりました

今回はより実践的な内容として、主成分分析 = PCA をやってみたいと思います

実装したコード（Livebook）はこちら

https://github.com/RyoWakabayashi/elixir-learning/blob/main/livebooks/explorer/pca.livemd

主成分分析とは、については以下の記事を参考にしてください

https://best-biostatistics.com/summary/pca.html

## 参考にした記事

https://laid-back-scientist.com/pca-imple

https://qiita.com/maskot1977/items/082557fcda78c4cdb41f

## 対象データ

今回はこちらのデータを使います

https://archive.ics.uci.edu/ml/datasets/wine

E 資格の勉強資料などでもよく使われる、ワインのデータです

３品種のワイン 178 本について、アルコール度数や色の濃さなど 13 項目の数値を計測しています

13 項目を学習するのは大変なので、 ２ 項目くらいに圧縮したいと思います

## 準備

必要なパッケージをインストールします

- Explorer: データ探索
- Kino.VegaLite: グラフ描画
- Nx: 行列計算

```elixir
Mix.install([
  {:explorer, "~> 0.2.0"},
  {:kino_vega_lite, "~> 0.1.1"},
  {:nx, "~> 0.3.0-dev", github: "elixir-nx/nx", branch: "main", sparse: "nx"}
])
```

エイリアスを付けておきます

```elixir
alias Explorer.DataFrame
alias Explorer.Series
alias VegaLite, as: Vl
```

## データロード

Explorer はワインのデータセットを含んでいます

ダウンロードしたり整形したりしなくていいので楽です

```explorer
wine_df = Explorer.Datasets.wine()
```

テーブル表示してみましょう

```elixir
DataFrame.table(wine_df, limit: :infinity)
```

![スクリーンショット 2022-07-21 14.14.32.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/8628165f-cab7-90f5-48c9-81dff6430a84.png)

項目数が多いので表示しきれていませんが、
左端にワインの品種を示す `class` があり、各項目が続いています

品種が3種類であることを確認しておきましょう

```elixir
wine_df
|> DataFrame.distinct(columns: ["class"])
|> DataFrame.pull("class")
```

![スクリーンショット 2022-07-21 14.17.40.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/682b6a97-a23c-8b5b-3a30-1826c174fa50.png)

この後頻繁に使うので、 `class` 以外の １３項目の名前を取っておきます

```elixir
cols =
  wine_df
  |> DataFrame.names()
  |> Enum.filter(fn name -> name != "class" end)
```

![スクリーンショット 2022-07-21 14.19.43.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/324bfe77-4b8f-a8c5-5030-5ab7e6bba67f.png)

## 統計情報の取得

主成分分析には必要ありませんが、 pandas の describe を Explorer で再現したいと思います

https://pandas.pydata.org/docs/reference/api/pandas.DataFrame.describe.html

主要な統計項目を取得することで、データの傾向がある程度見えます

まず、 describe で出力する各統計項目を取得する関数を定義します

* count: データ数
* mean: 平均
* std: 標準偏差
* min: 最小値
* 25%: 1/4 分位数
* 50%: 中央値
* 75%: 3/4 分位数
* max: 最大値

```elixir
get_stats = fn df, col ->
  series = DataFrame.pull(df, col)

  [
    DataFrame.n_rows(df),
    Series.mean(series),
    Series.std(series),
    Series.min(series),
    Series.quantile(series, 0.25),
    Series.median(series),
    Series.quantile(series, 0.75),
    Series.max(series)
  ]
end

get_stats.(wine_df, "alcohol")
```

![スクリーンショット 2022-07-21 14.21.30.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/5954c4c4-f0c4-9b5d-6066-1e91de118f2d.png)

これを各列に適用します（`class` は連続値ではないので対象外）

```elixir
add_stats_label = fn list ->
  [{"Stats", ["count", "mean", "std", "min", "25%", "50%", "75%", "max"]} | list]
end

cols
|> Enum.map(fn col ->
  {col, get_stats.(wine_df, col)}
end)
|> add_stats_label.()
|> DataFrame.new()
|> DataFrame.table(limit: :infinity)
```

![スクリーンショット 2022-07-21 14.24.53.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/cb3dc27a-dc3a-7e17-44d1-5d3cd3e9630b.png)

まさに pandas の describe と同じものが表示できました

## ヒストグラム

散布図行列を表示する準備として、まずヒストグラムを表示してみます

```elixir
get_values = fn df, col ->
  df
  |> DataFrame.pull(col)
  |> Series.to_list()
end

histgram = fn df, col ->
  x = get_values.(df, col)
  y = List.duplicate(1, DataFrame.n_rows(df))

  Vl.new(width: 300, height: 300, title: col)
  |> Vl.data_from_values(x: x, y: y)
  |> Vl.mark(:bar)
  |> Vl.encode_field(
    :x,
    "x",
    type: :quantitative,
    bin: %{maxbins: 20},
    title: col
  )
  |> Vl.encode_field(
    :y,
    "y",
    type: :quantitative,
    aggregate: :count
  )
end

histgram.(wine_df, "alcohol")
```

![visualization.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/6b80f47b-f7c5-bc5f-c8ab-af088f62a6b6.png)

いい感じに表示できました

全項目でやってみます

```elixir
histgram_list =
  cols
  |> Enum.map(fn col ->
    histgram.(wine_df, col)
  end)

Vl.new(width: 300, height: 300 * Enum.count(cols))
|> Vl.concat(histgram_list, :vertical)
```

![visualization (1).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/7f441eb6-7817-1284-43ce-036e7a8aa5c6.png)

続いて、クラス（ワインの品種）毎に出力してみましょう

```elixir
get_class_values = fn df, col, class ->
  df
  |> DataFrame.filter(Series.equal(df["class"], class))
  |> get_values.(col)
end

class_histgram = fn class, color ->
  Vl.new(width: 300, height: 300)
  |> Vl.mark(:bar, color: color, opacity: 0.5)
  |> Vl.encode_field(
    :x,
    "x#{class}",
    type: :quantitative,
    bin: %{maxbins: 20},
    title: "class#{class}"
  )
  |> Vl.encode_field(
    :y,
    "y",
    type: :quantitative,
    aggregate: :count
  )
end

all_class_histgram = fn df, col ->
  x1 = get_class_values.(df, col, 1)
  x2 = get_class_values.(df, col, 2)
  x3 = get_class_values.(df, col, 3)
  y = List.duplicate(0, DataFrame.n_rows(df))

  Vl.new(width: 300, height: 300, title: col)
  |> Vl.data_from_values(x1: x1, x2: x2, x3: x3, y: y)
  |> Vl.layers([
    class_histgram.(1, :blue),
    class_histgram.(2, :yellow),
    class_histgram.(3, :red)
  ])
end

all_class_histgram.(wine_df, "alcohol")

histgram_list =
  cols
  |> Enum.map(fn col ->
    all_class_histgram.(wine_df, col)
  end)

Vl.new(width: 300, height: 300 * Enum.count(cols))
|> Vl.concat(histgram_list, :vertical)
```

![visualization (2).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/15582ccc-cecb-fb39-ac14-598cc0a6b466.png)

結構クラス毎に特徴が出ている項目もありますね

## 散布図

次に散布図を作ります

```elixir
scatter = fn df, x_col, y_col ->
  x = get_values.(df, x_col)
  y = get_values.(df, y_col)
  class = get_values.(wine_df, "class")

  Vl.new(width: 300, height: 300)
  |> Vl.data_from_values(x: x, y: y, class: class)
  |> Vl.mark(:point)
  |> Vl.encode_field(:x, "x",
    type: :quantitative,
    scale: [domain: [Enum.min(x), Enum.max(x)]],
    title: x_col
  )
  |> Vl.encode_field(:y, "y",
    type: :quantitative,
    scale: [domain: [Enum.min(y), Enum.max(y)]],
    title: y_col
  )
  |> Vl.encode_field(:color, "class", type: :nominal)
end
```

`color` に `class` を指定することで、クラス毎に点の色を変えています

アルコールとリンゴ酸の散布図

```elixir
scatter.(wine_df, "alcohol", "malic_acid")
```

![visualization (3).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/77538d2d-90cb-6af6-5a57-48e76a5c7d14.png)

かなりバラバラで品種の分類は難しそうです

フラボノイドと色の濃さの散布図

```elixir
scatter.(wine_df, "flavanoids", "color_intensity")
```

![visualization (4).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/4a669caa-5e31-ee0f-dc24-77f40ed2ceb0.png)

クラス毎にまとまっていて、分類に使えそうな感じがします

## 散布図行列

準備ができたので、散布図行列を出してみましょう

```elixir
graphs =
  cols
  |> Enum.map(fn col_1 ->
    h_graphs =
      cols
      |> Enum.map(fn col_2 ->
        cond do
          col_1 == col_2 ->
            all_class_histgram.(wine_df, col_1)

          true ->
            scatter.(wine_df, col_1, col_2)
        end
      end)

    Vl.new(width: 300 * Enum.count(cols), height: 300)
    |> Vl.concat(h_graphs, :horizontal)
  end)

Vl.new(width: 300 * Enum.count(cols), height: 300 * Enum.count(cols))
|> Vl.concat(graphs, :vertical)
```

![visualization (5).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/0e424f13-3dd7-a87d-cf76-76716e25a054.png)

これを眺めれば、どの組み合わせが有効そうか、なんとなく見えてきます

## 標準化

この後相関行列を作るため、各項目の値を標準化しておきます

標準化とは、についてはこちらを参照してください

https://kotodori.jp/user-research/analytics/what-is-standardization/#:~:text=%E5%BD%B9%E7%AB%8B%E3%81%A1%E8%B3%87%E6%96%99%E3%83%95%E3%82%A9%E3%83%BC%E3%83%A0-,%E7%B5%B1%E8%A8%88%E3%81%AB%E3%81%8A%E3%81%91%E3%82%8B%E6%A8%99%E6%BA%96%E5%8C%96%E3%81%A8%E3%81%AF%EF%BC%9F,%E3%81%AF%E3%80%81%E6%AC%A1%E3%81%AE%E9%80%9A%E3%82%8A%E3%81%A7%E3%81%99%E3%80%82

標準化は各値について、平均を引いてから標準偏差で割ります

従って、各項目に対して以下のような関数で定義できます

```elixir
standardize = fn df, column ->
  mean =
    wine_df
    |> DataFrame.pull(column)
    |> Series.mean()

  std =
    wine_df
    |> DataFrame.pull(column)
    |> Series.std()

  df
  |> DataFrame.mutate(%{column => &Series.subtract(&1[column], mean)})
  |> DataFrame.mutate(%{column => &Series.divide(&1[column], std)})
end
```

この関数を各項目に対して適用します

```elixir
standardized_wine_df =
  cols
  |> Enum.reduce(wine_df, fn col, standardized_df ->
    standardize.(standardized_df, col)
  end)

standardized_wine_df
|> DataFrame.table(limit: :infinity)
```

![スクリーンショット 2022-07-21 14.48.24.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/631a3c31-c8c5-47d5-83f8-741d4bd26d8e.png)

## 相関行列

各項目がお互いにどれだけ関係しているか、を示す相関係数を全項目の組み合わせで出力します

前回の記事ではここを無理矢理ループして計算していましたが、
今回は Nx を使って行列計算でスマートに実装します

まず、データフレームをテンソルに変換します

```elixir
df_to_tensor = fn df ->
  df
  |> DataFrame.names()
  |> Enum.map(fn col ->
    standardized_wine_df
    |> DataFrame.pull(col)
    |> Series.to_tensor()
  end)
  |> Nx.concatenate()
  |> Nx.reshape({DataFrame.n_columns(df), DataFrame.n_rows(df)})
end

standardized_wine_tensor =
  standardized_wine_df
  |> DataFrame.select(["class"], :drop)
  |> df_to_tensor.()
  |> Nx.transpose()
```

![スクリーンショット 2022-07-21 14.52.58.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/0824a5dd-b499-d67d-c6de-b44ba37fef37.png)

データフレームと同じ形のテンソルになったことが分かると思います

分散共分散行列は標準化したテンソルを A とすると、以下のように表せます

```math
Cov = \frac{1}{n}　A^T \cdot A
```

```elixir
covariance_tensor =
  standardized_wine_tensor
  |> Nx.transpose()
  |> Nx.dot(standardized_wine_tensor)
  |> Nx.divide(DataFrame.n_rows(standardized_wine_df))
```

x と y の相関係数は以下の数式で表します

```math
r_{xy} = \frac{cov_{xy}}{s_x s_y}
```

```math
cov_{xy} : x と y の共分散
```

```math
s_{x} : x の標準偏差
```

```math
s_{y} : y の標準偏差
```

標準化済のため、各項目の標準偏差は 1 になっています
従って、共分散 = 相関係数になっています
そのため、 分散共分散行列 = 相関行列になります

取得した相関行列をテンソルからデータフレームに変換してテーブル表示してみます

```elixir
add_cols_label = fn list, cols_ ->
  [{"x", cols_} | list]
end

covariance_df =
  cols
  |> Stream.with_index()
  |> Enum.map(fn {col, index} ->
    {col, Nx.to_flat_list(covariance_tensor[index])}
  end)
  |> add_cols_label.(cols)
  |> DataFrame.new()

covariance_df
|> DataFrame.table(limit: :infinity)
```

![スクリーンショット 2022-07-21 15.15.04.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/f4cc7be8-5b78-7ccb-6187-e2206932118c.png)

これでも見にくいので、ヒートマップにしてみます

```elixir
heatmap =
  cols
  |> Stream.with_index()
  |> Enum.map(fn {col_1, index_1} ->
    cols
    |> Stream.with_index()
    |> Enum.map(fn {col_2, index_2} ->
      covariance =
        covariance_tensor[index_1][index_2]
        |> Nx.to_number()
      %{
        x: col_1,
        y: col_2,
        covariance: covariance
      }
    end)
  end)
  |> List.flatten()

Vl.new(width: 200, height: 200)
|> Vl.data_from_values(heatmap)
|> Vl.mark(:rect)
|> Vl.encode_field(:x, "x", type: :nominal)
|> Vl.encode_field(:y, "y", type: :nominal)
|> Vl.encode_field(
  :fill, "covariance",
  type: :quantitative,
  scale: [
    domain: [-1, 1],
    scheme: :blueorange
  ]
)
```

![visualization (6).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/1945e85d-e122-5470-c13d-394da3f544a9.png)

オレンジに近いところは正の相関、青に近いところは負の相関を持っています

対角線は同じ項目なので濃いオレンジです

項目間の関連具合が何となく見えますね

## 固有値、固有ベクトル

PCA のための固有値、固有ベクトルを取得します

これに関しては Nx がそのものずばりの関数を提供していたので、1行で済みました

しかも、固有値の大きい順に並べてくれているので、こちらで並べ替える必要がありません

```elixir
{eigenvals, eigenvecs} = Nx.LinAlg.eigh(covariance_tensor)
```

![スクリーンショット 2022-07-21 15.21.13.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/61309c05-1ea3-8783-a050-786238b2e61f.png)

## 寄与率

各固有ベクトル = 主成分得点がどれくらい情報量を持っているのか、寄与率を出してみます

```elixir
get_contribution_rate = fn index ->
  Nx.to_number(eigenvals[index - 1]) / Nx.size(eigenvals)
end

contribution_rate_list =
  1..Nx.size(eigenvals)
  |> Enum.to_list()
  |> Enum.map(fn index ->
    get_contribution_rate.(index)
  end)
```

![スクリーンショット 2022-07-21 15.22.58.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/bed54083-eb19-83c4-e0bb-27b25dc3c941.png)

そして、その累積値 = 累積寄与率を計算します

```elixir
cumulative_contribution_rate_list =
  1..Nx.size(eigenvals)
  |> Enum.map(fn index ->
    contribution_rate_list
    |> Enum.slice(0, index)
    |> Enum.sum()
  end)
```

![スクリーンショット 2022-07-21 15.23.56.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/117f81b6-17a7-1dcc-2ec9-931e0c9a2bb5.png)

寄与率、累積寄与率をテーブル表示してみましょう

```elixir
contribution_rate_df =
  %{
    "PC" => 1..Nx.size(eigenvals),
    "contribution_rate" => contribution_rate_list,
    "cumulative_contribution_rate" => cumulative_contribution_rate_list
  }
  |> DataFrame.new()

contribution_rate_df
|> DataFrame.table(limit: :infinity)
```

![スクリーンショット 2022-07-21 15.25.56.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/2e2d7023-eaaf-adf5-985d-e2c7b4adebe0.png)

第２主成分の時点で情報量の半数以上を占めていることが分かります

これをグラフ化してみましょう

```elixir
Vl.new(width: 300, height: 300)
|> Vl.data_from_values(
  x: 0..Nx.size(eigenvals),
  y: [ 0 | cumulative_contribution_rate_list]
)
|> Vl.mark(:line)
|> Vl.encode_field(
  :x, "x",
  type: :quantitative,
  title: "PC"
)
|> Vl.encode_field(
  :y, "y",
  type: :quantitative,
  title: "Cumulative Contribution Rate"
)
```

![visualization (7).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/5e2fc681-4f83-6106-16e1-5882b482e3a0.png)

今回は第２主成分までを使って次元圧縮してみましょう

## 次元圧縮

第２主成分得点まで=固有ベクトルの1番目と2番目を取り出してくっつけます

これを射影行列とします

```elixir
eigenvecs = Nx.transpose(eigenvecs)

w1 = eigenvecs[0]
w2 = eigenvecs[1]

IO.inspect(w1)
IO.inspect(w2)

w = Nx.stack([w1, w2], axis: 1)
```

標準化したデータと射影行列の内積によって、第１主成分と第２主成分を取得します

```elixir
pca_wine_tensor =
  standardized_wine_tensor
  |> Nx.dot(w)
  |> Nx.transpose()
```

![スクリーンショット 2022-07-21 15.31.02.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/4ef59d0c-fdc1-f3ae-7794-2cb1ee1a51c4.png)

主成分とクラスでデータフレームを作り、テーブル表示してみます

```elixir
classes = get_values.(wine_df, "class")

pca_wine_df =
  %{
    class: classes,
    pc1: Nx.to_flat_list(pca_wine_tensor[0]),
    pc2: Nx.to_flat_list(pca_wine_tensor[1])
  }
  |> DataFrame.new()

pca_wine_df
|> DataFrame.table(limit: :infinity)
```

![スクリーンショット 2022-07-21 15.32.09.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/e39457c3-34fc-ee5d-06ed-ea52cdfcd095.png)

今まで１３項目もあったせいで横方法に収まっていなかったテーブルが、
２項目まで次元圧縮されたおかげでスッキリ見えます

最後に第１主成分と第２主成分の散布図を表示します

```elixir
scatter.(pca_wine_df, "pc1", "pc2")
```

![visualization (8).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/22415df3-cce6-9ac3-93e5-dcaf1d1ce120.png)

明確にクラス毎に分かれてまとまっていることが分かります

この主成分を使えば、効率的に機械学習できそうです

## おわりに

データの抽出や集計、テーブル表示をしたいときは Explorer
グラフ表示したいときは Kino.VegaLite
行列計算したいときは Nx

という感じに行ったり来たりすることで色々なデータ分析ができました

Python だと pandas matplotlib numpy をそれぞれ使い分ける感じですね

Explorer と Nx があれば大概の統計はできると思うので、
まだまだ色々と使えそうです