scikit-fdaで関数主成分分析をしてみた

2022.11.07

データアナリティクス事業本部の鈴木です。

scikit-fdaはPythonで関数データ解析を行う際にとても便利なAPIを数多く備えています。

今回はその中でも関数主成分分析(以降FPCA)のAPIを使い、scikit-fdaの使い心地とFPCAの実行結果を確認します。

scikit-fdaのドキュメントの1ページであるFunctional Principal Component Analysisの内容に沿って進めていきます。

実行環境・ノートブック

今回はGoogle Colaboratory環境で実行しました。

Pythonおよびscikit-fdaのバージョンは以下になります。

  • Python: 3.7.15
  • scikit-fda: 0.7.1

ノートブックは以下のGithubレポジトリに格納しています。

関数主成分分析について

関数データ解析の分野では、多くの場合データをある基底関数の組み合わせで表現します。基底関数としては、B-スプライン関数やフーリエ基底などが使われることが多いです。

基底関数での表現については以下のドキュメントに例が記載されています。

一番簡単な場合は、ドキュメントに記載があるように、例えば以下のように基底関数を線型結合してデータを関数化します。

基底の線型結合

このように関数データ化することで、離散的なデータを平滑化できるだけでなく、微分が計算できるようになるのが大きな特徴です。

FPCAでは共分散関数に対して固有関数となる関数を求めます。

固有関数

ただし、共分散関数は、例えば以下のように推定されるものになります。

共分散関数

固有値が大きい固有関数が、関数データの特徴をよく反映した関数とされます。

やってみる

通しての実行結果はJupyter Notebookに記載しているので、ポイントを紹介していきます。

データを関数データ化する

今回はサンプルに従い、Berkeley Growth Studyデータセットを使います。このデータセットには男性39名、女性54名の1歳から18歳までの身長の推移が記録されています。

まずskfda.datasets.fetch_growth()にてFDataGridで取得します。ある時点をグリッドの1点とし、各時点における離散的なデータを格納するインスタンスです。

FDataGridにはデータを図示する機能も付いているので、scatterで散布図として可視化してみます。サンプルではplotを使い、データ点間を補完したように表示していますが、後から関数データ化した際に違いが分かるように散布図としました。

dataset = skfda.datasets.fetch_growth()
fd = dataset['data']
y = dataset['target']
fd.scatter()

散布図

続いて、基底関数を使って関数データ化します。今回はサンプルに倣ってB-スプライン関数を7本使います。

FDataGrid.to_basisにより、指定した基底関数を使った関数データへ変換することができます。

内部的には、指定した基底関数をどのように結合すれば元のデータが近似できるのか求めるために基底関数の係数を計算しないといけないはずですが、scikit-fdaではこのAPIを使うことで簡単に計算することが可能でとても便利です。

# B-スプライン関数のうち7本を基底関数として使う。
basis = skfda.representation.basis.BSpline(n_basis=7)
basis_fd = fd.to_basis(basis)
basis_fd.plot()

B-スプライン関数を使った関数データ

basis_fdインスタンスを呼び出すことで、どのような係数になっているか確認することもできます。

basis_fd
# FDataBasis(
#     basis=BSpline(domain_range=((1.0, 18.0),), n_basis=7, order=4, knots=(1.0, 5.25, 9.5, 13.75, 18.0)),
#     coefficients=[[ 81.62184377  95.32753609 122.2487525  142.51904585 184.48664294
#       199.4720117  193.26470272]
#      [ 76.82387771  94.30494976 113.90201189 141.02004067 150.21453901
#       177.65529306 178.70330065]
# (略)
#      [ 80.07666599  96.85011432 120.00975813 145.05857602 173.61369512
#       165.28674818 169.23374687]
#      [ 75.9228653   95.30102055 117.61810881 144.37605359 174.59219397
#       166.05675746 169.64664151]],
#     dataset_name=Berkeley Growth Study,
#     argument_names=('age',),
#     coordinate_names=('height',),
#     extrapolation=None)

関数主成分分析を実行する

FPCAを実行してみます。普通のPCAと同じイメージで、n_componentsで第何主成分関数まで使うか指定できます。また、求めた主成分関数の可視化も簡単にできます。

fpca = FPCA(n_components=2)
fpca.fit(basis_fd)
fpca.components_.plot()

FPCA実行後の主成分関数2つ

FPCAPlotにより、再構成した結果を確認することができます。

FPCAPlotの結果

サンプルでは、第一主成分関数は標本の身長の推移の平均を表しており、第二主成分関数は成長の傾向を表しているとしています。

ほかの基底関数も使ってみる

使用する基底関数系の例として、フーリエ基底を使ったものも紹介されていましたので確認してみました。

dataset = fetch_growth()
fd = dataset['data']
basis = skfda.representation.basis.Fourier(n_basis=7)
basis_fd = fd.to_basis(basis)
basis_fd.plot()

関数データ(フーリエ基底)

今回のデータは周期的なデータではないので、上手くはいかないようです。実際、B-スプライン基底を使った場合と比べて、関数データ化したデータが歪んでしまっています。

FPCAで取得した第一・第二主成分関数もまた、その傾向が表れていることが分かります。

fpca = FPCA(n_components=2)
fpca.fit(basis_fd)
fpca.components_.plot()

主成分関数(フーリエ基底)

最後に、主成分関数を使って再構成してみます。こちらも同様の歪みが現れていることが分かります。

FPCAPlot(
    basis_fd.mean(),
    fpca.components_,
    factor=30,
    fig=plt.figure(figsize=(6, 2 * 4)),
    n_rows=2,
).plot()

FPCAPlotの結果(フーリエ基底)

今回のようなデータはB-スプライン基底で関数データ化した方がよさそうなことが分かります。

最後に

今回はscikit-fdaの使い心地とFPCAの実行結果を確認してみました。少し慣れが必要そうですが、APIを使うだけで、好きな基底関数を使ってとても簡単にデータの変換ができ、すごく便利でした。

scikit-fdaのドキュメントには他にもさまざまな例が用意されているので、興味がある方はぜひ試してみてください。

参考資料