PCA xác suất

Xem trên TensorFlow.org Chạy trong Google Colab Xem nguồn trên GitHub Tải xuống sổ ghi chép

Xác suất chính thành phần phân tích toàn diện (PCA) là một kỹ thuật giảm chiều để phân tích dữ liệu thông qua một không gian tiềm ẩn thấp chiều ( Tipping và Đức Giám mục 1999 ). Nó thường được sử dụng khi thiếu các giá trị trong dữ liệu hoặc để chia tỷ lệ đa chiều.

Nhập khẩu

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')

Ngươi mâu

Hãy xem xét một tập dữ liệu \(\mathbf{X} = \{\mathbf{x}_n\}\) của \(N\) điểm dữ liệu, trong đó mỗi điểm dữ liệu là \(D\)chiều, $ \ mathbf {x} _n \ in \ mathbb {R} ^ D\(. We aim to represent each \)\ mathbf {x} _n $ dưới một biến tiềm ẩn \(\mathbf{z}_n \in \mathbb{R}^K\) với kích thước thấp hơn, $ K <D\(. The set of principal axes \)\ mathbf {W} $ liên hệ các biến tiềm ẩn các dữ liệu.

Cụ thể, chúng tôi giả định rằng mỗi biến tiềm ẩn được phân phối bình thường,

\[ \begin{equation*} \mathbf{z}_n \sim N(\mathbf{0}, \mathbf{I}). \end{equation*} \]

Điểm dữ liệu tương ứng được tạo thông qua phép chiếu,

\[ \begin{equation*} \mathbf{x}_n \mid \mathbf{z}_n \sim N(\mathbf{W}\mathbf{z}_n, \sigma^2\mathbf{I}), \end{equation*} \]

nơi ma trận \(\mathbf{W}\in\mathbb{R}^{D\times K}\) được gọi là trục chính. Trong PCA xác suất, chúng ta thường quan tâm đến việc lập dự toán các trục chính \(\mathbf{W}\) và thời hạn tiếng ồn\(\sigma^2\).

PCA xác suất tổng quát PCA cổ điển. Loại bỏ biến tiềm ẩn, phân phối của mỗi điểm dữ liệu là

\[ \begin{equation*} \mathbf{x}_n \sim N(\mathbf{0}, \mathbf{W}\mathbf{W}^\top + \sigma^2\mathbf{I}). \end{equation*} \]

PCA cổ điển là trường hợp cụ thể của PCA xác suất khi hiệp phương sai của tiếng ồn trở nên cực nhỏ, \(\sigma^2 \to 0\).

Chúng tôi thiết lập mô hình của chúng tôi bên dưới. Trong phân tích của chúng tôi, chúng tôi giả định \(\sigma\) được biết đến, và thay vì điểm lập dự toán \(\mathbf{W}\) như một tham số mô hình, chúng tôi đặt một trước khi qua nó để suy ra một bản phân phối trên trục chính. Chúng tôi sẽ thể hiện các mô hình như một TFP JointDistribution, cụ thể, chúng tôi sẽ sử dụng 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)

Dữ liệu

Chúng tôi có thể sử dụng mô hình để tạo dữ liệu bằng cách lấy mẫu từ phân phối trước chung.

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)

Chúng tôi hình dung tập dữ liệu.

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()

png

Suy luận Posteriori tối đa

Trước tiên, chúng tôi tìm kiếm ước tính điểm của các biến tiềm ẩn tối đa hóa mật độ xác suất sau. Đây được gọi là tối đa một posteriori (MAP) suy luận, và được thực hiện bằng cách tính toán các giá trị của \(\mathbf{W}\) và \(\mathbf{Z}\) đó phát huy tối đa mật độ sau \(p(\mathbf{W}, \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>]

png

Chúng ta có thể sử dụng mô hình để dữ liệu mẫu cho các giá trị suy cho \(\mathbf{W}\) và \(\mathbf{Z}\), và so sánh với số liệu thực tế, chúng tôi có điều kiện trên.

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)>

png

Suy luận biến đổi

MAP có thể được sử dụng để tìm chế độ (hoặc một trong các chế độ) của phân phối sau, nhưng không cung cấp bất kỳ thông tin chi tiết nào khác về nó. Tiếp theo chúng ta sử dụng suy luận variational, nơi distribtion sau \(p(\mathbf{W}, \mathbf{Z} \mid \mathbf{X})\) là xấp xỉ bằng một bản phân phối biến phân \(q(\mathbf{W}, \mathbf{Z})\) parametrised bởi \(\boldsymbol{\lambda}\). Mục đích là để tìm các thông số biến phân \(\boldsymbol{\lambda}\) đó giảm thiểu sự phân kỳ KL giữa q và sau, \(\mathrm{KL}(q(\mathbf{W}, \mathbf{Z}) \mid\mid p(\mathbf{W}, \mathbf{Z} \mid \mathbf{X}))\), hoặc tương đương, tối đa hóa các bằng chứng thấp hơn bị ràng buộc, \(\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)>

png

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()

png

Sự nhìn nhận

Hướng dẫn này ban đầu được viết bằng Edward 1.0 ( nguồn ). Chúng tôi cảm ơn tất cả những người đóng góp để viết và sửa đổi phiên bản đó.

Người giới thiệu

[1]: Michael E. Tipping và Christopher M. Bishop. Phân tích thành phần chính xác suất. Tạp chí của Hiệp hội thống kê Hoàng gia: Series B (Phương pháp thống kê), 61 (3): 611-622, 1999.