 عرض على TensorFlow.org |  تشغيل في Google Colab |  عرض المصدر على جيثب |  تحميل دفتر |
في هذا colab ، نستكشف انحدار عملية Gaussian باستخدام TensorFlow و TensorFlow Probability. نقوم بإنشاء بعض الملاحظات الصاخبة من بعض الوظائف المعروفة ونلائم نماذج GP لتلك البيانات. ثم نقوم بأخذ عينات من GP الخلفي ورسم قيم الدالة التي تم أخذ العينات عليها عبر الشبكات في مجالاتها.
خلفية
دعونا \(\mathcal{X}\) أن يكون أي مجموعة. عملية التمويه (GP) هي عبارة عن مجموعة من المتغيرات العشوائية فهرستها من قبل \(\mathcal{X}\) بحيث إذا\(\{X_1, \ldots, X_n\} \subset \mathcal{X}\) هو أي مجموعة فرعية محدودة، وهامشية كثافة\(p(X_1 = x_1, \ldots, X_n = x_n)\) هو متعدد المتغيرات التمويه. يتم تحديد أي توزيع غاوسي تمامًا من خلال اللحظات المركزية الأولى والثانية (المتوسط والتغاير) ، ولا تعد توزيعات GP استثناءً. يمكننا تحديد GP تماما من حيث معدلها ظيفة \(\mu : \mathcal{X} \to \mathbb{R}\) والتغاير وظيفة\(k : \mathcal{X} \times \mathcal{X} \to \mathbb{R}\). يتم تغليف معظم القوة التعبيرية لـ GP في اختيار وظيفة التغاير. لأسباب مختلفة، ويشار إلى وظيفة التغاير أيضا بوصفها وظيفة النواة. هو مطلوب منها فقط أن يكون متماثل وإيجابية-محددة (انظر الفصل 4 من راسموسن ويليامز ). أدناه نستفيد من نواة التغاير الرباعي الأس. شكله
\[ k(x, x') := \sigma^2 \exp \left( \frac{\|x - x'\|^2}{\lambda^2} \right) \]
حيث \(\sigma^2\) يسمى "السعة" و \(\lambda\) مقياس الطول. يمكن تحديد معلمات kernel عبر إجراء تحسين الاحتمالية القصوى.
تتكون عينة كاملة من GP وظيفة قيمتها الحقيقية على كامل مساحة\(\mathcal{X}\) وهي في الواقع غير عملي لتحقيق. غالبًا ما يختار المرء مجموعة من النقاط التي يلاحظ عندها عينة ويرسم قيم الوظيفة في هذه النقاط. يتم تحقيق ذلك عن طريق أخذ عينات من غاوسي مناسب (محدود الأبعاد) متعدد المتغيرات.
لاحظ أنه وفقًا للتعريف أعلاه ، فإن أي توزيع غاوسي متعدد المتغيرات بأبعاد محدودة هو أيضًا عملية غاوسية. عادة، عندما يشير أحد إلى GP، فمن ضمني بأن مجموعة المؤشر بعض \(\mathbb{R}^n\)وسوف نقوم في الواقع هذا الافتراض هنا.
من التطبيقات الشائعة لعمليات Gaussian في التعلم الآلي انحدار عملية Gaussian. والفكرة هي أن نتمنى لتقدير وظيفة غير معروفة نظرا الملاحظات صاخبة \(\{y_1, \ldots, y_N\}\) وظيفة في عدد محدود من نقاط \(\{x_1, \ldots x_N\}.\) لنا أن نتصور عملية توليدي
\[ \begin{align} f \sim \: & \textsf{GaussianProcess}\left( \text{mean_fn}=\mu(x), \text{covariance_fn}=k(x, x')\right) \\ y_i \sim \: & \textsf{Normal}\left( \text{loc}=f(x_i), \text{scale}=\sigma\right), i = 1, \ldots, N \end{align} \]
كما هو مذكور أعلاه ، من المستحيل حساب دالة العينة ، لأننا سنطلب قيمتها بعدد لا نهائي من النقاط. بدلاً من ذلك ، يعتبر المرء عينة محدودة من غاوسي متعدد المتغيرات.
\[ \begin{gather} \begin{bmatrix} f(x_1) \\ \vdots \\ f(x_N) \end{bmatrix} \sim \textsf{MultivariateNormal} \left( \: \text{loc}= \begin{bmatrix} \mu(x_1) \\ \vdots \\ \mu(x_N) \end{bmatrix} \:,\: \text{scale}= \begin{bmatrix} k(x_1, x_1) & \cdots & k(x_1, x_N) \\ \vdots & \ddots & \vdots \\ k(x_N, x_1) & \cdots & k(x_N, x_N) \\ \end{bmatrix}^{1/2} \: \right) \end{gather} \\ y_i \sim \textsf{Normal} \left( \text{loc}=f(x_i), \text{scale}=\sigma \right) \]
ملاحظة الأس \(\frac{1}{2}\) على مصفوفة التغاير: هذا يدل على التحلل Cholesky. يعد تضمين Cholesky أمرًا ضروريًا لأن MVN هو توزيع عائلي على نطاق الموقع. مما يؤسف له أن التحلل Cholesky مكلفة حسابيا، مع \(O(N^3)\) الوقت و \(O(N^2)\) الفضاء. يركز الكثير من أدبيات GP على التعامل مع هذا الأس الصغير الذي يبدو غير ضار.
من الشائع اعتبار دالة المتوسط السابقة ثابتة ، وغالبًا ما تكون صفرًا. أيضا ، بعض الاصطلاحات التوضيحية مريحة. واحد كثيرا ما يكتب \(\mathbf{f}\) لناقلات محدود من القيم وظيفة عينات. ويستخدم عدد من الرموز مثيرة للاهتمام لمصفوفة التغاير الناتجة عن تطبيق \(k\) إلى أزواج من المدخلات. التالية (Quiñonero-كانديلا، 2005) ، نلاحظ أن عناصر المصفوفة هي التغايرات من القيم وظيفة في بعينها نقاط الإدخال. وهكذا يمكننا أن دلالة على مصفوفة التغاير كما \(K_{AB}\)حيث \(A\) و \(B\) بعض المؤشرات من مجموعة من القيم وظيفة على طول أبعاد مصفوفة معينة.
على سبيل المثال، فإن البيانات الملحوظ \((\mathbf{x}, \mathbf{y})\) مع القيم وظيفة الكامنة ضمنية \(\mathbf{f}\)، نستطيع كتابة
\[ K_{\mathbf{f},\mathbf{f} } = \begin{bmatrix} k(x_1, x_1) & \cdots & k(x_1, x_N) \\ \vdots & \ddots & \vdots \\ k(x_N, x_1) & \cdots & k(x_N, x_N) \\ \end{bmatrix} \]
وبالمثل ، يمكننا مزج مجموعات من المدخلات ، كما في
\[ K_{\mathbf{f},*} = \begin{bmatrix} k(x_1, x^*_1) & \cdots & k(x_1, x^*_T) \\ \vdots & \ddots & \vdots \\ k(x_N, x^*_1) & \cdots & k(x_N, x^*_T) \\ \end{bmatrix} \]
حيث أننا نفترض أن هناك \(N\) مدخلات التدريب، و \(T\) مدخلات الاختبار. يمكن بعد ذلك كتابة العملية التوليدية المذكورة أعلاه بشكل مضغوط
\[ \begin{align} \mathbf{f} \sim \: & \textsf{MultivariateNormal} \left( \text{loc}=\mathbf{0}, \text{scale}=K_{\mathbf{f},\mathbf{f} }^{1/2} \right) \\ y_i \sim \: & \textsf{Normal} \left( \text{loc}=f_i, \text{scale}=\sigma \right), i = 1, \ldots, N \end{align} \]
عملية أخذ العينات في السطر الأول ينتج مجموعة محدودة من \(N\) القيم وظيفة من التمويه المتعدد - ليست وظيفة كاملة كما هو الحال في ما سبق GP رسم التدوين. السطر الثاني يصف مجموعة من \(N\) يستمد من Gaussians حيد المتغير تركز على مختلف القيم وظيفة، مع ثابت الضوضاء المراقبة \(\sigma^2\).
مع وجود النموذج التوليدي أعلاه في مكانه الصحيح ، يمكننا المضي قدمًا في النظر في مشكلة الاستدلال اللاحق. ينتج عن هذا توزيع لاحق على قيم الوظيفة في مجموعة جديدة من نقاط الاختبار ، بشرط البيانات الصاخبة المرصودة من العملية أعلاه.
مع تدوين المشار إليه أعلاه، يمكننا مضغوط كتابة توزيع التنبؤي الخلفي على مستقبل (صاخبة) ملاحظات مشروطة على المدخلات وبيانات التدريب المقابلة على النحو التالي (لمزيد من التفاصيل، انظر §2.2 من راسموسن ويليامز ).
\[ \mathbf{y}^* \mid \mathbf{x}^*, \mathbf{x}, \mathbf{y} \sim \textsf{Normal} \left( \text{loc}=\mathbf{\mu}^*, \text{scale}=(\Sigma^*)^{1/2} \right), \]
أين
\[ \mathbf{\mu}^* = K_{*,\mathbf{f} }\left(K_{\mathbf{f},\mathbf{f} } + \sigma^2 I \right)^{-1} \mathbf{y} \]
و
\[ \Sigma^* = K_{*,*} - K_{*,\mathbf{f} } \left(K_{\mathbf{f},\mathbf{f} } + \sigma^2 I \right)^{-1} K_{\mathbf{f},*} \]
الواردات
import time
import numpy as np
import matplotlib.pyplot as plt
import tensorflow.compat.v2 as tf
import tensorflow_probability as tfp
tfb = tfp.bijectors
tfd = tfp.distributions
tfk = tfp.math.psd_kernels
tf.enable_v2_behavior()
from mpl_toolkits.mplot3d import Axes3D
%pylab inline
# Configure plot defaults
plt.rcParams['axes.facecolor'] = 'white'
plt.rcParams['grid.color'] = '#666666'
%config InlineBackend.figure_format = 'png'
Populating the interactive namespace from numpy and matplotlib
مثال: انحدار GP الدقيق على البيانات الجيبية الصاخبة
هنا نقوم بتوليد بيانات التدريب من الجيب الصاخب ، ثم أخذ عينات من مجموعة من المنحنيات من الجزء الخلفي لنموذج الانحدار GP. نحن نستخدم آدم لتحسين hyperparameters نواة (لتقليل احتمال السجل السلبي للبيانات تحت قبل). نرسم منحنى التدريب ، متبوعًا بالوظيفة الحقيقية والعينات الخلفية.
def sinusoid(x):
return np.sin(3 * np.pi * x[..., 0])
def generate_1d_data(num_training_points, observation_noise_variance):
"""Generate noisy sinusoidal observations at a random set of points.
Returns:
observation_index_points, observations
"""
index_points_ = np.random.uniform(-1., 1., (num_training_points, 1))
index_points_ = index_points_.astype(np.float64)
# y = f(x) + noise
observations_ = (sinusoid(index_points_) +
np.random.normal(loc=0,
scale=np.sqrt(observation_noise_variance),
size=(num_training_points)))
return index_points_, observations_
# Generate training data with a known noise level (we'll later try to recover
# this value from the data).
NUM_TRAINING_POINTS = 100
observation_index_points_, observations_ = generate_1d_data(
num_training_points=NUM_TRAINING_POINTS,
observation_noise_variance=.1)
سوف نضع مقدمو الاديره على hyperparameters النواة، وكتابة التوزيع المشترك للhyperparameters والبيانات المرصودة باستخدام tfd.JointDistributionNamed .
def build_gp(amplitude, length_scale, observation_noise_variance):
"""Defines the conditional dist. of GP outputs, given kernel parameters."""
# Create the covariance kernel, which will be shared between the prior (which we
# use for maximum likelihood training) and the posterior (which we use for
# posterior predictive sampling)
kernel = tfk.ExponentiatedQuadratic(amplitude, length_scale)
# Create the GP prior distribution, which we will use to train the model
# parameters.
return tfd.GaussianProcess(
kernel=kernel,
index_points=observation_index_points_,
observation_noise_variance=observation_noise_variance)
gp_joint_model = tfd.JointDistributionNamed({
'amplitude': tfd.LogNormal(loc=0., scale=np.float64(1.)),
'length_scale': tfd.LogNormal(loc=0., scale=np.float64(1.)),
'observation_noise_variance': tfd.LogNormal(loc=0., scale=np.float64(1.)),
'observations': build_gp,
})
يمكننا التحقق من صحة تنفيذنا عن طريق التحقق من أنه يمكننا أخذ عينات من العينة السابقة ، وحساب كثافة اللوغاريتمات للعينة.
x = gp_joint_model.sample()
lp = gp_joint_model.log_prob(x)
print("sampled {}".format(x))
print("log_prob of sample: {}".format(lp))
sampled {'observation_noise_variance': <tf.Tensor: shape=(), dtype=float64, numpy=2.067952217184325>, 'length_scale': <tf.Tensor: shape=(), dtype=float64, numpy=1.154435715487831>, 'amplitude': <tf.Tensor: shape=(), dtype=float64, numpy=5.383850737703549>, 'observations': <tf.Tensor: shape=(100,), dtype=float64, numpy=
array([-2.37070577, -2.05363838, -0.95152824, 3.73509388, -0.2912646 ,
0.46112342, -1.98018513, -2.10295857, -1.33589756, -2.23027226,
-2.25081374, -0.89450835, -2.54196452, 1.46621647, 2.32016193,
5.82399989, 2.27241034, -0.67523432, -1.89150197, -1.39834474,
-2.33954116, 0.7785609 , -1.42763627, -0.57389025, -0.18226098,
-3.45098732, 0.27986652, -3.64532398, -1.28635204, -2.42362875,
0.01107288, -2.53222176, -2.0886136 , -5.54047694, -2.18389607,
-1.11665628, -3.07161217, -2.06070336, -0.84464262, 1.29238438,
-0.64973999, -2.63805504, -3.93317576, 0.65546645, 2.24721181,
-0.73403676, 5.31628298, -1.2208384 , 4.77782252, -1.42978168,
-3.3089274 , 3.25370494, 3.02117591, -1.54862932, -1.07360811,
1.2004856 , -4.3017773 , -4.95787789, -1.95245901, -2.15960839,
-3.78592731, -1.74096185, 3.54891595, 0.56294143, 1.15288455,
-0.77323696, 2.34430694, -1.05302007, -0.7514684 , -0.98321063,
-3.01300144, -3.00033274, 0.44200837, 0.45060886, -1.84497318,
-1.89616746, -2.15647664, -2.65672581, -3.65493379, 1.70923375,
-3.88695218, -0.05151283, 4.51906677, -2.28117003, 3.03032793,
-1.47713194, -0.35625273, 3.73501587, -2.09328047, -0.60665614,
-0.78177188, -0.67298545, 2.97436033, -0.29407932, 2.98482427,
-1.54951178, 2.79206821, 4.2225733 , 2.56265198, 2.80373284])>}
log_prob of sample: -194.96442183797524
الآن دعنا نقوم بالتحسين للعثور على قيم المعلمات ذات الاحتمال اللاحق الأعلى. سنحدد متغيرًا لكل معلمة ، ونجعل قيمها موجبة.
# Create the trainable model parameters, which we'll subsequently optimize.
# Note that we constrain them to be strictly positive.
constrain_positive = tfb.Shift(np.finfo(np.float64).tiny)(tfb.Exp())
amplitude_var = tfp.util.TransformedVariable(
initial_value=1.,
bijector=constrain_positive,
name='amplitude',
dtype=np.float64)
length_scale_var = tfp.util.TransformedVariable(
initial_value=1.,
bijector=constrain_positive,
name='length_scale',
dtype=np.float64)
observation_noise_variance_var = tfp.util.TransformedVariable(
initial_value=1.,
bijector=constrain_positive,
name='observation_noise_variance_var',
dtype=np.float64)
trainable_variables = [v.trainable_variables[0] for v in
[amplitude_var,
length_scale_var,
observation_noise_variance_var]]
لشرط نموذج على البيانات المرصودة لدينا، فإننا سوف تحديد target_log_prob وظيفة، والتي تأخذ (لا يزال يتعين الاستدلال) hyperparameters النواة.
def target_log_prob(amplitude, length_scale, observation_noise_variance):
return gp_joint_model.log_prob({
'amplitude': amplitude,
'length_scale': length_scale,
'observation_noise_variance': observation_noise_variance,
'observations': observations_
})
# Now we optimize the model parameters.
num_iters = 1000
optimizer = tf.optimizers.Adam(learning_rate=.01)
# Use `tf.function` to trace the loss for more efficient evaluation.
@tf.function(autograph=False, jit_compile=False)
def train_model():
with tf.GradientTape() as tape:
loss = -target_log_prob(amplitude_var, length_scale_var,
observation_noise_variance_var)
grads = tape.gradient(loss, trainable_variables)
optimizer.apply_gradients(zip(grads, trainable_variables))
return loss
# Store the likelihood values during training, so we can plot the progress
lls_ = np.zeros(num_iters, np.float64)
for i in range(num_iters):
loss = train_model()
lls_[i] = loss
print('Trained parameters:')
print('amplitude: {}'.format(amplitude_var._value().numpy()))
print('length_scale: {}'.format(length_scale_var._value().numpy()))
print('observation_noise_variance: {}'.format(observation_noise_variance_var._value().numpy()))
Trained parameters: amplitude: 0.9176153445125278 length_scale: 0.18444082442910079 observation_noise_variance: 0.0880273312850989
# Plot the loss evolution
plt.figure(figsize=(12, 4))
plt.plot(lls_)
plt.xlabel("Training iteration")
plt.ylabel("Log marginal likelihood")
plt.show()
# Having trained the model, we'd like to sample from the posterior conditioned
# on observations. We'd like the samples to be at points other than the training
# inputs.
predictive_index_points_ = np.linspace(-1.2, 1.2, 200, dtype=np.float64)
# Reshape to [200, 1] -- 1 is the dimensionality of the feature space.
predictive_index_points_ = predictive_index_points_[..., np.newaxis]
optimized_kernel = tfk.ExponentiatedQuadratic(amplitude_var, length_scale_var)
gprm = tfd.GaussianProcessRegressionModel(
kernel=optimized_kernel,
index_points=predictive_index_points_,
observation_index_points=observation_index_points_,
observations=observations_,
observation_noise_variance=observation_noise_variance_var,
predictive_noise_variance=0.)
# Create op to draw 50 independent samples, each of which is a *joint* draw
# from the posterior at the predictive_index_points_. Since we have 200 input
# locations as defined above, this posterior distribution over corresponding
# function values is a 200-dimensional multivariate Gaussian distribution!
num_samples = 50
samples = gprm.sample(num_samples)
# Plot the true function, observations, and posterior samples.
plt.figure(figsize=(12, 4))
plt.plot(predictive_index_points_, sinusoid(predictive_index_points_),
label='True fn')
plt.scatter(observation_index_points_[:, 0], observations_,
label='Observations')
for i in range(num_samples):
plt.plot(predictive_index_points_, samples[i, :], c='r', alpha=.1,
label='Posterior Sample' if i == 0 else None)
leg = plt.legend(loc='upper right')
for lh in leg.legendHandles:
lh.set_alpha(1)
plt.xlabel(r"Index points ($\mathbb{R}^1$)")
plt.ylabel("Observation space")
plt.show()
تهميش المعلمات الفوقية مع مؤسسة حمد الطبية
بدلاً من تحسين المعلمات الفائقة ، دعنا نحاول دمجها مع Hamiltonian Monte Carlo. سنقوم أولاً بتعريف وتشغيل أداة أخذ العينات للرسم تقريبًا من التوزيع اللاحق عبر المعلمات الفوقية للنواة ، بالنظر إلى الملاحظات.
num_results = 100
num_burnin_steps = 50
sampler = tfp.mcmc.TransformedTransitionKernel(
tfp.mcmc.NoUTurnSampler(
target_log_prob_fn=target_log_prob,
step_size=tf.cast(0.1, tf.float64)),
bijector=[constrain_positive, constrain_positive, constrain_positive])
adaptive_sampler = tfp.mcmc.DualAveragingStepSizeAdaptation(
inner_kernel=sampler,
num_adaptation_steps=int(0.8 * num_burnin_steps),
target_accept_prob=tf.cast(0.75, tf.float64))
initial_state = [tf.cast(x, tf.float64) for x in [1., 1., 1.]]
# Speed up sampling by tracing with `tf.function`.
@tf.function(autograph=False, jit_compile=False)
def do_sampling():
return tfp.mcmc.sample_chain(
kernel=adaptive_sampler,
current_state=initial_state,
num_results=num_results,
num_burnin_steps=num_burnin_steps,
trace_fn=lambda current_state, kernel_results: kernel_results)
t0 = time.time()
samples, kernel_results = do_sampling()
t1 = time.time()
print("Inference ran in {:.2f}s.".format(t1-t0))
Inference ran in 9.00s.
دعونا نتحقق من سلامة جهاز أخذ العينات من خلال فحص آثار المعامل الفائق.
(amplitude_samples,
length_scale_samples,
observation_noise_variance_samples) = samples
f = plt.figure(figsize=[15, 3])
for i, s in enumerate(samples):
ax = f.add_subplot(1, len(samples) + 1, i + 1)
ax.plot(s)
الآن بدلا من بناء GP واحد مع hyperparameters الأمثل، ونحن بناء على توزيع التنبؤي الخلفي على شكل مزيج من الأطباء، التي يحددها عينة من توزيع الخلفي على hyperparameters لكل منهما. يتكامل هذا تقريبًا مع المعلمات الخلفية عبر أخذ عينات مونت كارلو لحساب التوزيع التنبئي الهامشي في المواقع غير المرصودة.
# The sampled hyperparams have a leading batch dimension, `[num_results, ...]`,
# so they construct a *batch* of kernels.
batch_of_posterior_kernels = tfk.ExponentiatedQuadratic(
amplitude_samples, length_scale_samples)
# The batch of kernels creates a batch of GP predictive models, one for each
# posterior sample.
batch_gprm = tfd.GaussianProcessRegressionModel(
kernel=batch_of_posterior_kernels,
index_points=predictive_index_points_,
observation_index_points=observation_index_points_,
observations=observations_,
observation_noise_variance=observation_noise_variance_samples,
predictive_noise_variance=0.)
# To construct the marginal predictive distribution, we average with uniform
# weight over the posterior samples.
predictive_gprm = tfd.MixtureSameFamily(
mixture_distribution=tfd.Categorical(logits=tf.zeros([num_results])),
components_distribution=batch_gprm)
num_samples = 50
samples = predictive_gprm.sample(num_samples)
# Plot the true function, observations, and posterior samples.
plt.figure(figsize=(12, 4))
plt.plot(predictive_index_points_, sinusoid(predictive_index_points_),
label='True fn')
plt.scatter(observation_index_points_[:, 0], observations_,
label='Observations')
for i in range(num_samples):
plt.plot(predictive_index_points_, samples[i, :], c='r', alpha=.1,
label='Posterior Sample' if i == 0 else None)
leg = plt.legend(loc='upper right')
for lh in leg.legendHandles:
lh.set_alpha(1)
plt.xlabel(r"Index points ($\mathbb{R}^1$)")
plt.ylabel("Observation space")
plt.show()
على الرغم من أن الاختلافات طفيفة في هذه الحالة ، بشكل عام ، فإننا نتوقع أن يتم تعميم التوزيع التنبئي اللاحق بشكل أفضل (يعطي احتمالية أعلى للبيانات المعلقة) من مجرد استخدام المعلمات الأكثر احتمالية كما فعلنا أعلاه.