 TensorFlow.orgで表示 |  GoogleColabで実行 |  GitHubでソースを表示 |  ノートブックをダウンロード |
確率的主成分分析（PCA）は、低次元の潜在空間を介してデータを分析する次元削減手法です（ Tipping and Bishop 1999 ）。これは、データに欠落値がある場合、または多次元尺度構成法でよく使用されます。
輸入
import functools
import warnings
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
import tensorflow.compat.v2 as tf
import tensorflow_probability as tfp
from tensorflow_probability import bijectors as tfb
from tensorflow_probability import distributions as tfd
tf.enable_v2_behavior()
plt.style.use("ggplot")
warnings.filterwarnings('ignore')
モデル
$ N $データポイントのデータセット$ \ mathbf {X} = {\ mathbf {x} _n} $について考えます。ここで、各データポイントは$ D $次元であり、$ \ mathbf {x} _n \ in \ mathbb { R} ^ D $。各$ \ mathbf {x} _n $を、潜在変数$ \ mathbf {z} _n \ in \ mathbb {R} ^ K $の下で、より低い次元、$ K <D $で表すことを目指しています。主軸のセット$ \ mathbf {W} $は、潜在変数をデータに関連付けます。
具体的には、各潜在変数が正規分布していると仮定します。
対応するデータポイントは、投影を介して生成されます。
ここで、行列$ \ mathbf {W} \ in \ mathbb {R} ^ {D \ times K} $は主軸として知られています。確率的PCAでは、通常、主軸$ \ mathbf {W} $とノイズ項$ \ sigma ^ 2 $を推定することに関心があります。
確率的PCAは、古典的なPCAを一般化します。潜在変数をマージナル化すると、各データポイントの分布は次のようになります。
古典的PCAは、ノイズの共分散が非常に小さくなり、$ \ sigma ^ 2 \ to 0 $になる確率的PCAの特定のケースです。
以下にモデルを設定します。この分析では、$ \ sigma $が既知であると想定し、モデルパラメーターとして$ \ mathbf {W} $を点推定する代わりに、主軸上の分布を推測するために、その上に事前分布を配置します。モデルをTFPJointDistributionとして表現します。具体的には、 JointDistributionCoroutineAutoBatchedを使用します。
def probabilistic_pca(data_dim, latent_dim, num_datapoints, stddv_datapoints):
w = yield tfd.Normal(loc=tf.zeros([data_dim, latent_dim]),
scale=2.0 * tf.ones([data_dim, latent_dim]),
name="w")
z = yield tfd.Normal(loc=tf.zeros([latent_dim, num_datapoints]),
scale=tf.ones([latent_dim, num_datapoints]),
name="z")
x = yield tfd.Normal(loc=tf.matmul(w, z),
scale=stddv_datapoints,
name="x")
num_datapoints = 5000
data_dim = 2
latent_dim = 1
stddv_datapoints = 0.5
concrete_ppca_model = functools.partial(probabilistic_pca,
data_dim=data_dim,
latent_dim=latent_dim,
num_datapoints=num_datapoints,
stddv_datapoints=stddv_datapoints)
model = tfd.JointDistributionCoroutineAutoBatched(concrete_ppca_model)
データ
モデルを使用して、共同事前分布からサンプリングすることでデータを生成できます。
actual_w, actual_z, x_train = model.sample()
print("Principal axes:")
print(actual_w)
Principal axes: tf.Tensor( [[ 2.2801023] [-1.1619819]], shape=(2, 1), dtype=float32)
データセットを視覚化します。
plt.scatter(x_train[0, :], x_train[1, :], color='blue', alpha=0.1)
plt.axis([-20, 20, -20, 20])
plt.title("Data set")
plt.show()
最大事後推論
最初に、事後確率密度を最大化する潜在変数の点推定を検索します。これは最大事後（MAP）推論として知られており、事後密度$ p（\ mathbf {W}、\を最大化する$ \ mathbf {W} $と$ \ mathbf {Z} $の値を計算することによって行われます。 mathbf {Z} \ mid \ mathbf {X}）\ propto p（\ mathbf {W}、\ mathbf {Z}、\ mathbf {X}）$。
w = tf.Variable(tf.random.normal([data_dim, latent_dim]))
z = tf.Variable(tf.random.normal([latent_dim, num_datapoints]))
target_log_prob_fn = lambda w, z: model.log_prob((w, z, x_train))
losses = tfp.math.minimize(
lambda: -target_log_prob_fn(w, z),
optimizer=tf.optimizers.Adam(learning_rate=0.05),
num_steps=200)
plt.plot(losses)
[<matplotlib.lines.Line2D at 0x7f19897a42e8>]
モデルを使用して、$ \ mathbf {W} $と$ \ mathbf {Z} $の推定値のデータをサンプリングし、条件付けした実際のデータセットと比較できます。
print("MAP-estimated axes:")
print(w)
_, _, x_generated = model.sample(value=(w, z, None))
plt.scatter(x_train[0, :], x_train[1, :], color='blue', alpha=0.1, label='Actual data')
plt.scatter(x_generated[0, :], x_generated[1, :], color='red', alpha=0.1, label='Simulated data (MAP)')
plt.legend()
plt.axis([-20, 20, -20, 20])
plt.show()
MAP-estimated axes:
<tf.Variable 'Variable:0' shape=(2, 1) dtype=float32, numpy=
array([[ 2.9135954],
[-1.4826864]], dtype=float32)>
変分推論
MAPは、事後分布のモード（またはモードの1つ）を見つけるために使用できますが、それに関する他の洞察は提供しません。次に、変分推論を使用します。ここで、事後分布$ p（\ mathbf {W}、\ mathbf {Z} \ mid \ mathbf {X}）$は、変分分布$ q（\ mathbf {W}、\ mathbfを使用して近似されます。 {Z}）$は$ \ boldsymbol {\ lambda} $によってパラメータ化されます。目的は、qと後部の間のKL発散を最小化する変分パラメーター$ \ boldsymbol {\ lambda} $を見つけることです。$ \ mathrm {KL}（q（\ mathbf {W}、\ mathbf {Z}）\ mid \ mid p（\ mathbf {W}、\ mathbf {Z} \ mid \ mathbf {X}））$、または同等に、証拠の下限を最大化する$ \ mathbb {E} _ {q（\ mathbf {W }、\ mathbf {Z}; \ boldsymbol {\ lambda}）} \ left [\ log p（\ mathbf {W}、\ mathbf {Z}、\ mathbf {X}）-\ log q（\ mathbf {W }、\ mathbf {Z}; \ boldsymbol {\ lambda}）\ right] $。
qw_mean = tf.Variable(tf.random.normal([data_dim, latent_dim]))
qz_mean = tf.Variable(tf.random.normal([latent_dim, num_datapoints]))
qw_stddv = tfp.util.TransformedVariable(1e-4 * tf.ones([data_dim, latent_dim]),
bijector=tfb.Softplus())
qz_stddv = tfp.util.TransformedVariable(
1e-4 * tf.ones([latent_dim, num_datapoints]),
bijector=tfb.Softplus())
def factored_normal_variational_model():
qw = yield tfd.Normal(loc=qw_mean, scale=qw_stddv, name="qw")
qz = yield tfd.Normal(loc=qz_mean, scale=qz_stddv, name="qz")
surrogate_posterior = tfd.JointDistributionCoroutineAutoBatched(
factored_normal_variational_model)
losses = tfp.vi.fit_surrogate_posterior(
target_log_prob_fn,
surrogate_posterior=surrogate_posterior,
optimizer=tf.optimizers.Adam(learning_rate=0.05),
num_steps=200)
print("Inferred axes:")
print(qw_mean)
print("Standard Deviation:")
print(qw_stddv)
plt.plot(losses)
plt.show()
Inferred axes:
<tf.Variable 'Variable:0' shape=(2, 1) dtype=float32, numpy=
array([[ 2.4168603],
[-1.2236133]], dtype=float32)>
Standard Deviation:
<TransformedVariable: dtype=float32, shape=[2, 1], fn="softplus", numpy=
array([[0.0042499 ],
[0.00598824]], dtype=float32)>
posterior_samples = surrogate_posterior.sample(50)
_, _, x_generated = model.sample(value=(posterior_samples))
# It's a pain to plot all 5000 points for each of our 50 posterior samples, so
# let's subsample to get the gist of the distribution.
x_generated = tf.reshape(tf.transpose(x_generated, [1, 0, 2]), (2, -1))[:, ::47]
plt.scatter(x_train[0, :], x_train[1, :], color='blue', alpha=0.1, label='Actual data')
plt.scatter(x_generated[0, :], x_generated[1, :], color='red', alpha=0.1, label='Simulated data (VI)')
plt.legend()
plt.axis([-20, 20, -20, 20])
plt.show()
謝辞
このチュートリアルは、もともとEdward 1.0（ソース）で書かれました。そのバージョンの作成と改訂に貢献してくれたすべての人に感謝します。
参考文献
[1]：Michael E.TippingとChristopherM.Bishop。確率的主成分分析。 Journal of the Royal Statistics Society：Series B（Statistical Methodology） 、61（3）：611-622、1999。