LoginSignup
26
5

More than 1 year has passed since last update.

ElixirAdvent Calendar 2022Day 2
@RyoWakabayashi(若林 諒)in株式会社オーイーシー

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

Last updated at Posted at 2022-07-21

はじめに

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

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

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

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

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

参考にした記事

対象データ

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

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

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

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

準備

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

  • Explorer: データ探索
  • Kino.VegaLite: グラフ描画
  • Nx: 行列計算
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"}
])

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

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

データロード

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

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

wine_df = Explorer.Datasets.wine()

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

DataFrame.table(wine_df, limit: :infinity)

スクリーンショット 2022-07-21 14.14.32.png

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

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

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

スクリーンショット 2022-07-21 14.17.40.png

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

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

スクリーンショット 2022-07-21 14.19.43.png

統計情報の取得

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

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

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

  • count: データ数
  • mean: 平均
  • std: 標準偏差
  • min: 最小値
  • 25%: 1/4 分位数
  • 50%: 中央値
  • 75%: 3/4 分位数
  • max: 最大値
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

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

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

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

ヒストグラム

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

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

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

全項目でやってみます

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

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

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

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

散布図

次に散布図を作ります

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 を指定することで、クラス毎に点の色を変えています

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

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

visualization (3).png

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

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

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

visualization (4).png

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

散布図行列

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

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

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

標準化

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

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

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

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

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

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

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

相関行列

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

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

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

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

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

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

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

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

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

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

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

s_{x} : x の標準偏差

s_{y} : y の標準偏差

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

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

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

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

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

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

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

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

固有値、固有ベクトル

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

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

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

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

スクリーンショット 2022-07-21 15.21.13.png

寄与率

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

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

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

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

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

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

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

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

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

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

次元圧縮

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

これを射影行列とします

eigenvecs = Nx.transpose(eigenvecs)

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

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

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

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

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

スクリーンショット 2022-07-21 15.31.02.png

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

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

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

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

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

visualization (8).png

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

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

おわりに

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

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

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

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

26
5
1

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
26
5