TensorFlow प्रायिकता केस स्टडी: सहप्रसरण अनुमान

TensorFlow.org पर देखें Google Colab में चलाएं GitHub पर स्रोत देखेंनोटबुक डाउनलोड करें

मैंने इस नोटबुक को TensorFlow प्रायिकता सीखने के लिए केस स्टडी के रूप में लिखा था। जिस समस्या को मैंने हल करने के लिए चुना है, वह 2-डी माध्य 0 गाऊसी यादृच्छिक चर के नमूनों के लिए एक सहप्रसरण मैट्रिक्स का अनुमान लगा रही है। समस्या में कुछ अच्छी विशेषताएं हैं:

  • यदि हम सहसंयोजक (एक सामान्य दृष्टिकोण) से पहले एक उलटा विशार्ट का उपयोग करते हैं, तो समस्या का एक विश्लेषणात्मक समाधान है, इसलिए हम अपने परिणामों की जांच कर सकते हैं।
  • समस्या में एक विवश पैरामीटर का नमूना लेना शामिल है, जो कुछ दिलचस्प जटिलता जोड़ता है।
  • सबसे सीधा समाधान सबसे तेज़ नहीं है, इसलिए कुछ अनुकूलन कार्य करना है।

जैसे-जैसे मैं आगे बढ़ा, मैंने अपने अनुभवों को लिखने का फैसला किया। मुझे अपने सिर को TFP के बारीक बिंदुओं के चारों ओर लपेटने में थोड़ा समय लगा, इसलिए यह नोटबुक काफी सरलता से शुरू होती है और फिर धीरे-धीरे अधिक जटिल TFP सुविधाओं तक काम करती है। मैं रास्ते में बहुत सारी समस्याओं में भाग गया, और मैंने उन दोनों प्रक्रियाओं को पकड़ने की कोशिश की है जिससे मुझे उन्हें पहचानने में मदद मिली और अंततः मुझे मिले कामकाज। मैं (परीक्षण यकीन है कि अलग-अलग चरणों सही हैं के बहुत सारे सहित) विस्तार की बहुत सारी शामिल करने के लिए कोशिश की है।

TensorFlow प्रायिकता क्यों सीखें?

मैंने कुछ कारणों से अपने प्रोजेक्ट के लिए TensorFlow प्रायिकता को आकर्षक पाया:

  • TensorFlow प्रायिकता आपको प्रोटोटाइप जटिल मॉडल को एक नोटबुक में अंतःक्रियात्मक रूप से विकसित करने देती है। आप अपने कोड को छोटे टुकड़ों में तोड़ सकते हैं जिन्हें आप अंतःक्रियात्मक रूप से और यूनिट परीक्षणों के साथ परीक्षण कर सकते हैं।
  • एक बार जब आप विस्तार करने के लिए तैयार हो जाते हैं, तो आप हमारे पास मौजूद सभी बुनियादी ढांचे का लाभ उठा सकते हैं ताकि TensorFlow को कई मशीनों पर कई, अनुकूलित प्रोसेसर पर चलाया जा सके।
  • अंत में, जबकि मुझे वास्तव में स्टेन पसंद है, मुझे डीबग करना काफी मुश्किल लगता है। आपको अपने सभी मॉडलिंग कोड को एक स्टैंडअलोन भाषा में लिखना होगा जिसमें आपको अपने कोड को पोक करने, मध्यवर्ती राज्यों का निरीक्षण करने आदि के लिए बहुत कम टूल हों।

नकारात्मक पक्ष यह है कि TensorFlow प्रायिकता स्टैन और PyMC3 की तुलना में बहुत नई है, इसलिए प्रलेखन एक कार्य प्रगति पर है, और बहुत सारी कार्यक्षमता है जिसे अभी बनाया जाना बाकी है। खुशी की बात है कि मैंने टीएफपी की नींव को ठोस पाया, और इसे एक मॉड्यूलर तरीके से डिजाइन किया गया है जो इसकी कार्यक्षमता को काफी सरलता से विस्तारित करने की अनुमति देता है। इस नोटबुक में, केस स्टडी को हल करने के अलावा, मैं टीएफपी को विस्तारित करने के कुछ तरीके दिखाऊंगा।

यह किसके लिए है

मैं मान रहा हूँ कि पाठक इस नोटबुक में कुछ महत्वपूर्ण पूर्वापेक्षाएँ लेकर आ रहे हैं। तुम्हे करना चाहिए:

  • बायेसियन अनुमान की मूल बातें जानें। (आप नहीं करते हैं, वास्तव में एक अच्छा पहली पुस्तक है सांख्यिकीय पुनर्विचार )
  • एक एमसीएमसी नमूना पुस्तकालय, जैसे से कुछ परिचित हैं स्टेन / PyMC3 / बग
  • का एक ठोस समझ है NumPy (एक अच्छा परिचय है डेटा विश्लेषण के लिए अजगर )
  • साथ कम से कम गुजर परिचित पर है TensorFlow , लेकिन जरूरी नहीं कि विशेषज्ञता। ( लर्निंग TensorFlow अच्छा है, लेकिन TensorFlow के तेजी से विकास का अर्थ है सबसे किताबें थोड़ा दिनांकित किया जाएगा। स्टैनफोर्ड CS20 कोर्स भी अच्छा है।)

पहला प्रयास

यहाँ समस्या पर मेरा पहला प्रयास है। स्पोइलर: मेरा समाधान काम नहीं करता है, और चीजों को ठीक करने के लिए कई प्रयास करने जा रहे हैं! यद्यपि इस प्रक्रिया में कुछ समय लगता है, नीचे दिया गया प्रत्येक प्रयास TFP के एक नए भाग को सीखने के लिए उपयोगी रहा है।

एक नोट: TFP वर्तमान में व्युत्क्रम विशार्ट वितरण को लागू नहीं करता है (हम अंत में देखेंगे कि अपने स्वयं के व्युत्क्रम विशार्ट को कैसे रोल किया जाए), इसलिए इसके बजाय मैं एक विशर्ट का उपयोग करके एक सटीक मैट्रिक्स का अनुमान लगाने की समस्या को बदल दूंगा।

import collections
import math
import os
import time

import numpy as np
import pandas as pd
import scipy
import scipy.stats
import matplotlib.pyplot as plt

import tensorflow.compat.v2 as tf
tf.enable_v2_behavior()
import tensorflow_probability as tfp
tfd = tfp.distributions
tfb = tfp.bijectors

चरण 1: प्रेक्षणों को एक साथ प्राप्त करें

यहां मेरा डेटा सभी सिंथेटिक है, इसलिए यह वास्तविक दुनिया के उदाहरण की तुलना में थोड़ा कठिन प्रतीत होगा। हालांकि, ऐसा कोई कारण नहीं है कि आप अपना खुद का कुछ सिंथेटिक डेटा नहीं बना सकते।

युक्ति: जब आप अपने मॉडल के रूप में तय कर लिया है, तो आप कुछ पैरामीटर मान लेने और अपने चुने हुए मॉडल का उपयोग कुछ सिंथेटिक डेटा उत्पन्न करने के लिए कर सकते हैं। अपने कार्यान्वयन की एक विवेकपूर्ण जाँच के रूप में, फिर आप यह सत्यापित कर सकते हैं कि आपके अनुमानों में आपके द्वारा चुने गए मापदंडों के सही मान शामिल हैं। अपने डिबगिंग/परीक्षण चक्र को तेज़ बनाने के लिए, आप अपने मॉडल के सरलीकृत संस्करण पर विचार कर सकते हैं (उदाहरण के लिए कम आयामों या कम नमूनों का उपयोग करें)।

युक्ति: यह NumPy सरणी के रूप में अपनी टिप्पणियों के साथ काम करने के लिए सबसे आसान है। ध्यान देने वाली एक महत्वपूर्ण बात यह है कि डिफ़ॉल्ट रूप से NumPy फ्लोट64 का उपयोग करता है, जबकि TensorFlow डिफ़ॉल्ट रूप से फ्लोट32 का उपयोग करता है।

सामान्य तौर पर, TensorFlow संचालन चाहता है कि सभी तर्क एक ही प्रकार के हों, और आपको प्रकार बदलने के लिए स्पष्ट डेटा कास्टिंग करनी होगी। यदि आप फ्लोट64 अवलोकनों का उपयोग करते हैं, तो आपको बहुत सारे कास्ट ऑपरेशंस जोड़ने होंगे। इसके विपरीत, NumPy स्वचालित रूप से कास्टिंग का ख्याल रखेगा। इसलिए, यह बहुत आसान float32 में अपने Numpy डेटा कन्वर्ट करने के लिए की तुलना में यह प्रयोग float64 को TensorFlow मजबूर करने के लिए है।

कुछ पैरामीटर मान चुनें

# We're assuming 2-D data with a known true mean of (0, 0)
true_mean = np.zeros([2], dtype=np.float32)
# We'll make the 2 coordinates correlated
true_cor = np.array([[1.0, 0.9], [0.9, 1.0]], dtype=np.float32)
# And we'll give the 2 coordinates different variances
true_var = np.array([4.0, 1.0], dtype=np.float32)
# Combine the variances and correlations into a covariance matrix
true_cov = np.expand_dims(np.sqrt(true_var), axis=1).dot(
    np.expand_dims(np.sqrt(true_var), axis=1).T) * true_cor
# We'll be working with precision matrices, so we'll go ahead and compute the
# true precision matrix here
true_precision = np.linalg.inv(true_cov)

# Here's our resulting covariance matrix
print(true_cov)
# Verify that it's positive definite, since np.random.multivariate_normal
# complains about it not being positive definite for some reason.
# (Note that I'll be including a lot of sanity checking code in this notebook -
# it's a *huge* help for debugging)
print('eigenvalues: ', np.linalg.eigvals(true_cov))

[[4.  1.8]
 [1.8 1. ]]
eigenvalues:  [4.843075   0.15692513]

कुछ सिंथेटिक अवलोकन उत्पन्न करें

नोट TensorFlow संभावना सम्मेलन है कि आपके डेटा की प्रारंभिक आयाम (रों) नमूना सूचकांक प्रतिनिधित्व का उपयोग करता है, और अपने डेटा के अंतिम आयाम (रों) अपने नमूनों की आयामी स्वरूप प्रतिनिधित्व करते हैं।

यहाँ हम चाहते हैं 100 नमूनों की लंबाई का एक वेक्टर है, जिनमें से प्रत्येक 2. हम एक सरणी उत्पन्न करेंगे my_data आकार (100, 2) के साथ। my_data[i, :] है \(i\)वें नमूना है, और यह लंबाई 2 का एक वेक्टर है।

(बनाने के लिए याद रखें my_data प्रकार float32 है!)

# Set the seed so the results are reproducible.
np.random.seed(123)

# Now generate some observations of our random variable.
# (Note that I'm suppressing a bunch of spurious about the covariance matrix
# not being positive semidefinite via check_valid='ignore' because it really is
# positive definite!)
my_data = np.random.multivariate_normal(
    mean=true_mean, cov=true_cov, size=100,
    check_valid='ignore').astype(np.float32)

my_data.shape

(100, 2)

विवेक टिप्पणियों की जाँच करें

बग का एक संभावित स्रोत आपके सिंथेटिक डेटा को खराब कर रहा है! आइए कुछ आसान जांच करते हैं।

# Do a scatter plot of the observations to make sure they look like what we
# expect (higher variance on the x-axis, y values strongly correlated with x)
plt.scatter(my_data[:, 0], my_data[:, 1], alpha=0.75)
plt.show()

पीएनजी

print('mean of observations:', np.mean(my_data, axis=0))
print('true mean:', true_mean)

mean of observations: [-0.24009615 -0.16638893]
true mean: [0. 0.]

print('covariance of observations:\n', np.cov(my_data, rowvar=False))
print('true covariance:\n', true_cov)

covariance of observations:
 [[3.95307734 1.68718486]
 [1.68718486 0.94910269]]
true covariance:
 [[4.  1.8]
 [1.8 1. ]]

ठीक है, हमारे नमूने उचित दिखते हैं। अगला कदम।

चरण 2: NumPy में संभावना फ़ंक्शन को लागू करें

टीएफ प्रोबेबिलिटी में अपना एमसीएमसी सैंपलिंग करने के लिए हमें जो मुख्य चीज लिखनी होगी, वह एक लॉग संभावना फ़ंक्शन है। आम तौर पर NumPy की तुलना में TF लिखना थोड़ा मुश्किल होता है, इसलिए मुझे NumPy में प्रारंभिक कार्यान्वयन करने में मदद मिलती है। मैं 2 टुकड़े, एक डेटा संभावना समारोह के अनुरूप है में संभावना समारोह विभाजित करने के लिए जा रहा हूँ \(P(data | parameters)\) और एक पूर्व संभावना समारोह मेल खाती है कि \(P(parameters)\)।

ध्यान दें कि इन NumPy फ़ंक्शंस को सुपर अनुकूलित / वेक्टरकृत करने की आवश्यकता नहीं है क्योंकि लक्ष्य केवल परीक्षण के लिए कुछ मान उत्पन्न करना है। शुद्धता महत्वपूर्ण विचार है!

सबसे पहले हम डेटा लॉग संभावना टुकड़ा लागू करेंगे। यह बहुत सीधा है। याद रखने वाली एक बात यह है कि हम सटीक मैट्रिसेस के साथ काम करने जा रहे हैं, इसलिए हम तदनुसार पैरामीटर करेंगे।

def log_lik_data_numpy(precision, data):
  # np.linalg.inv is a really inefficient way to get the covariance matrix, but
  # remember we don't care about speed here
  cov = np.linalg.inv(precision)
  rv = scipy.stats.multivariate_normal(true_mean, cov)
  return np.sum(rv.logpdf(data))

# test case: compute the log likelihood of the data given the true parameters
log_lik_data_numpy(true_precision, my_data)

-280.81822950593767

हम वहाँ के बाद से पीछे के लिए एक विश्लेषणात्मक समाधान एक विशार्ट पूर्व परिशुद्धता मैट्रिक्स के लिए उपयोग करने के लिए (देखें जा रहे हैं संयुग्म महंतों के विकिपीडिया के काम तालिका )।

विशार्ट वितरण 2 पैरामीटर:

  • स्वतंत्रता की डिग्री की संख्या (लेबल \(\nu\) विकिपीडिया में)
  • पैमाने मैट्रिक्स (लेबल \(V\) विकिपीडिया में)

मानकों के साथ एक विशार्ट वितरण के लिए मतलब \(\nu, V\) है \(E[W] = \nu V\), और विचरण है \(\text{Var}(W_{ij}) = \nu(v_{ij}^2+v_{ii}v_{jj})\)

कुछ उपयोगी अंतर्ज्ञान: आप उत्पन्न करके एक विशार्ट नमूना उत्पन्न कर सकते हैं \(\nu\) स्वतंत्र ड्रॉ \(x_1 \ldots x_{\nu}\) मतलब 0 और सहप्रसरण के साथ एक मल्टीवेरिएट सामान्य यादृच्छिक चर से \(V\) और फिर राशि के गठन \(W = \sum_{i=1}^{\nu} x_i x_i^T\)।

आप उन्हें से विभाजित करके विशार्ट नमूने rescale तो \(\nu\), आप में से नमूना सहप्रसरण मैट्रिक्स मिल \(x_i\)। यह नमूना सहप्रसरण मैट्रिक्स की ओर जाते हैं चाहिए \(V\) रूप \(\nu\) बढ़ जाती है। जब \(\nu\) छोटा है, वहाँ नमूना सहप्रसरण मैट्रिक्स में बदलाव के बहुत सारे, के इतने छोटे मान है \(\nu\) कमजोर महंतों और के बड़े मूल्यों के अनुरूप \(\nu\) मजबूत महंतों के अनुरूप हैं। ध्यान दें कि \(\nu\) कम से कम अंतरिक्ष आप कर रहे हैं नमूना या आप विलक्षण मैट्रिक्स उत्पन्न करेंगे के आयाम के रूप में बड़े रूप में किया जाना चाहिए।

हम इस्तेमाल करेंगे \(\nu = 3\) तो हम एक कमजोर से पहले है, और हम ले लेंगे \(V = \frac{1}{\nu} I\) जो (याद है कि इसका मतलब है पहचान की ओर हमारे सहप्रसरण अनुमान खींच लेंगे \(\nu V\))।

PRIOR_DF = 3
PRIOR_SCALE = np.eye(2, dtype=np.float32) / PRIOR_DF

def log_lik_prior_numpy(precision):
  rv = scipy.stats.wishart(df=PRIOR_DF, scale=PRIOR_SCALE)
  return rv.logpdf(precision)

# test case: compute the prior for the true parameters
log_lik_prior_numpy(true_precision)

-9.103606346649766

विशार्ट वितरण में जाना जाता है मतलब के साथ एक मल्टीवेरिएट सामान्य की शुद्धता मैट्रिक्स के आकलन के लिए संयुग्म पहले है \(\mu\)।

मान लीजिए पहले विशार्ट मापदंड हैं \(\nu, V\) और कहा कि हमारे पास \(n\) हमारे मल्टीवेरिएट सामान्य, की टिप्पणियों \(x_1, \ldots, x_n\)। पीछे मापदंड हैं \(n + \nu, \left(V^{-1} + \sum_{i=1}^n (x_i-\mu)(x_i-\mu)^T \right)^{-1}\)।

n = my_data.shape[0]
nu_prior = PRIOR_DF
v_prior = PRIOR_SCALE
nu_posterior = nu_prior + n
v_posterior = np.linalg.inv(np.linalg.inv(v_prior) + my_data.T.dot(my_data))
posterior_mean = nu_posterior * v_posterior
v_post_diag = np.expand_dims(np.diag(v_posterior), axis=1)
posterior_sd = np.sqrt(nu_posterior *
                       (v_posterior ** 2.0 + v_post_diag.dot(v_post_diag.T)))

पोस्टीरियर और सच्चे मूल्यों की एक त्वरित साजिश। ध्यान दें कि पोस्टीरियर नमूना पोस्टीरियर के करीब हैं लेकिन पहचान की ओर थोड़ा सा सिकुड़ा हुआ है। यह भी ध्यान दें कि वास्तविक मान पोस्टीरियर के मोड से बहुत दूर हैं - संभवतः ऐसा इसलिए है क्योंकि पहले हमारे डेटा के लिए बहुत अच्छा मेल नहीं है। एक वास्तविक समस्या में हम संभावना एक सहप्रसरण के लिए पहले विशार्ट उलटा बढ़ाया की तरह कुछ के साथ बेहतर कर चाहते हैं (उदाहरण के लिए देखें, एंड्रयू गेल्मैन की कमेंटरी विषय पर), लेकिन फिर हम एक अच्छा विश्लेषणात्मक पीछे नहीं होगा।

sample_precision = np.linalg.inv(np.cov(my_data, rowvar=False, bias=False))
fig, axes = plt.subplots(2, 2)
fig.set_size_inches(10, 10)
for i in range(2):
  for j in range(2):
    ax = axes[i, j]
    loc = posterior_mean[i, j]
    scale = posterior_sd[i, j]
    xmin = loc - 3.0 * scale
    xmax = loc + 3.0 * scale
    x = np.linspace(xmin, xmax, 1000)
    y = scipy.stats.norm.pdf(x, loc=loc, scale=scale)
    ax.plot(x, y)
    ax.axvline(true_precision[i, j], color='red', label='True precision')
    ax.axvline(sample_precision[i, j], color='red', linestyle=':', label='Sample precision')
    ax.set_title('precision[%d, %d]' % (i, j))
plt.legend()
plt.show()

पीएनजी

चरण 3: TensorFlow में संभावना फ़ंक्शन को लागू करें

स्पोइलर: हमारा पहला प्रयास काम नहीं करेगा; हम नीचे क्यों बात करेंगे।

युक्ति: उपयोग TensorFlow उत्सुक मोड जब आपके संभावना कार्यों का विकास। तुरंत, तो आप सहभागी का उपयोग करने के बजाय डीबग कर सकते हैं सब कुछ कार्यान्वित - उत्सुक मोड TF व्यवहार अधिक NumPy की तरह बना देता है Session.run() । नोट देखें यहाँ ।

प्रारंभिक: वितरण वर्ग

TFP में वितरण वर्गों का एक संग्रह है जिसका उपयोग हम अपनी लॉग संभावनाओं को उत्पन्न करने के लिए करेंगे। ध्यान देने वाली एक बात यह है कि ये वर्ग केवल एक नमूने के बजाय दसियों नमूनों के साथ काम करते हैं - यह वेक्टरकरण और संबंधित गति के लिए अनुमति देता है।

एक वितरण नमूने के टेंसर के साथ 2 अलग-अलग तरीकों से काम कर सकता है। इन 2 तरीकों को एक ठोस उदाहरण के साथ स्पष्ट करना सबसे आसान है जिसमें एकल स्केलर पैरामीटर के साथ वितरण शामिल है। मैं इस्तेमाल करेंगे प्वासों वितरण, जो एक है rate पैरामीटर।

  • अगर हम के लिए एक एकल मूल्य के साथ एक प्वासों बनाने rate पैरामीटर, इसके लिए एक कॉल sample() विधि एक भी मान। यह मान एक कहा जाता है event , और इस मामले में घटनाओं सभी scalars हैं।
  • अगर हम के लिए मूल्यों की एक टेन्सर के साथ एक प्वासों बनाने rate पैरामीटर, इसके लिए एक कॉल sample() विधि अब एक से अधिक मान देता है, दर टेन्सर में प्रत्येक मान के लिए एक। वस्तु स्वतंत्र Poissons का एक संग्रह, अपने स्वयं के दर के साथ एक के रूप में कार्य करता है, और मूल्यों से प्रत्येक के लिए एक कॉल द्वारा दिया sample() इन Poissons में से एक से मेल खाती है। स्वतंत्र नहीं बल्कि हूबहू वितरित की घटनाओं का यह संकलन एक कहा जाता है batch ।
  • sample() विधि एक लेता है sample_shape एक खाली टपल करने के लिए पैरामीटर जो चूक। के लिए एक गैर-रिक्त मान पासिंग sample_shape नमूने में परिणाम कई बैचों लौटने। बैचों के इस संग्रह एक कहा जाता है sample ।

एक वितरण के log_prob() विधि एक तरह से है कि समानताएं कैसे में डेटा की खपत sample() यह उत्पन्न करता है। log_prob() कई के लिए नमूने के लिए संभावनाओं, यानी, घटनाओं की स्वतंत्र बैचों देता है।

  • हम अपने प्वासों वस्तु एक अदिश से बनाया गया था कि है, तो rate , प्रत्येक बैच एक अदिश है, और अगर हम नमूनों की एक टेन्सर में पारित, हम लॉग संभावनाओं का एक ही आकार के एक टेन्सर बाहर मिल जाएगा।
  • हम अपने प्वासों वस्तु आकार का एक टेन्सर से बनाया गया था कि है, तो T की rate मूल्यों, प्रत्येक बैच आकार का एक टेन्सर है T । यदि हम आकार डी, टी के नमूने के एक टेंसर में गुजरते हैं, तो हम आकार डी, टी की लॉग संभावनाओं का एक टेंसर प्राप्त करेंगे।

नीचे कुछ उदाहरण दिए गए हैं जो इन मामलों को स्पष्ट करते हैं। देखें इस नोटबुक की घटनाओं, बैचों, और आकार पर एक अधिक विस्तृत ट्यूटोरियल के लिए।

# case 1: get log probabilities for a vector of iid draws from a single
# normal distribution
norm1 = tfd.Normal(loc=0., scale=1.)
probs1 = norm1.log_prob(tf.constant([1., 0.5, 0.]))

# case 2: get log probabilities for a vector of independent draws from
# multiple normal distributions with different parameters.  Note the vector
# values for loc and scale in the Normal constructor.
norm2 = tfd.Normal(loc=[0., 2., 4.], scale=[1., 1., 1.])
probs2 = norm2.log_prob(tf.constant([1., 0.5, 0.]))

print('iid draws from a single normal:', probs1.numpy())
print('draws from a batch of normals:', probs2.numpy())

iid draws from a single normal: [-1.4189385 -1.0439385 -0.9189385]
draws from a batch of normals: [-1.4189385 -2.0439386 -8.918939 ]

डेटा लॉग संभावना

पहले हम डेटा लॉग संभावना फ़ंक्शन को लागू करेंगे।

VALIDATE_ARGS = True
ALLOW_NAN_STATS = False

NumPy मामले से एक महत्वपूर्ण अंतर यह है कि हमारे TensorFlow संभावना फ़ंक्शन को केवल एकल मैट्रिक्स के बजाय सटीक मैट्रिक्स के वैक्टर को संभालने की आवश्यकता होगी। जब हम कई श्रृंखलाओं से नमूना लेते हैं तो पैरामीटर के वेक्टर का उपयोग किया जाएगा।

हम एक वितरण ऑब्जेक्ट बनाएंगे जो सटीक मैट्रिसेस के एक बैच के साथ काम करता है (अर्थात प्रति श्रृंखला एक मैट्रिक्स)।

हमारे डेटा की लॉग संभावनाओं की गणना करते समय, हमें अपने डेटा को उसी तरीके से दोहराने की आवश्यकता होगी जैसे हमारे पैरामीटर ताकि प्रति बैच चर एक प्रति हो। हमारे दोहराए गए डेटा का आकार इस प्रकार होना चाहिए:

[sample shape, batch shape, event shape]

हमारे मामले में, घटना का आकार 2 है (चूंकि हम 2-डी गाऊसी के साथ काम कर रहे हैं)। नमूना आकार 100 है, क्योंकि हमारे पास 100 नमूने हैं। बैच आकार केवल उन सटीक मैट्रिक्स की संख्या होगी जिनके साथ हम काम कर रहे हैं। हर बार जब हम संभाव्यता फ़ंक्शन कहते हैं, तो डेटा को दोहराना बेकार है, इसलिए हम डेटा को पहले से दोहराएंगे और प्रतिकृति संस्करण में पास करेंगे।

ध्यान दें कि यह एक अकुशल दिया गया है: MultivariateNormalFullCovariance कुछ विकल्प के लिए महंगा है कि हम अंत में अनुकूलन अनुभाग में के बारे में बात करेंगे।

def log_lik_data(precisions, replicated_data):
  n = tf.shape(precisions)[0]  # number of precision matrices
  # We're estimating a precision matrix; we have to invert to get log
  # probabilities.  Cholesky inversion should be relatively efficient,
  # but as we'll see later, it's even better if we can avoid doing the Cholesky
  # decomposition altogether.
  precisions_cholesky = tf.linalg.cholesky(precisions)
  covariances = tf.linalg.cholesky_solve(
      precisions_cholesky, tf.linalg.eye(2, batch_shape=[n]))
  rv_data = tfd.MultivariateNormalFullCovariance(
      loc=tf.zeros([n, 2]),
      covariance_matrix=covariances,
      validate_args=VALIDATE_ARGS,
      allow_nan_stats=ALLOW_NAN_STATS)

  return tf.reduce_sum(rv_data.log_prob(replicated_data), axis=0)

# For our test, we'll use a tensor of 2 precision matrices.
# We'll need to replicate our data for the likelihood function.
# Remember, TFP wants the data to be structured so that the sample dimensions
# are first (100 here), then the batch dimensions (2 here because we have 2
# precision matrices), then the event dimensions (2 because we have 2-D
# Gaussian data).  We'll need to add a middle dimension for the batch using
# expand_dims, and then we'll need to create 2 replicates in this new dimension
# using tile.
n = 2
replicated_data = np.tile(np.expand_dims(my_data, axis=1), reps=[1, 2, 1])
print(replicated_data.shape)

(100, 2, 2)

सुझाव: एक बात मैं अत्यंत उपयोगी मेरी TensorFlow कार्यों का थोड़ा विवेक चेकों लिख रहा है होना करने के लिए मिल गया है। TF में वैश्वीकरण को गड़बड़ाना वास्तव में आसान है, इसलिए TF आउटपुट को सत्यापित करने के लिए NumPy के आसपास सरल कार्य करना एक शानदार तरीका है। इन्हें छोटे इकाई परीक्षणों के रूप में सोचें।

# check against the numpy implementation
precisions = np.stack([np.eye(2, dtype=np.float32), true_precision])
n = precisions.shape[0]
lik_tf = log_lik_data(precisions, replicated_data=replicated_data).numpy()

for i in range(n):
  print(i)
  print('numpy:', log_lik_data_numpy(precisions[i], my_data))
  print('tensorflow:', lik_tf[i])

0
numpy: -430.71218815801365
tensorflow: -430.71207
1
numpy: -280.81822950593767
tensorflow: -280.8182

पूर्व लॉग संभावना

पूर्व आसान है क्योंकि हमें डेटा प्रतिकृति के बारे में चिंता करने की ज़रूरत नहीं है।

@tf.function(autograph=False)
def log_lik_prior(precisions):
  rv_precision = tfd.WishartTriL(
      df=PRIOR_DF,
      scale_tril=tf.linalg.cholesky(PRIOR_SCALE),
      validate_args=VALIDATE_ARGS,
      allow_nan_stats=ALLOW_NAN_STATS)
  return rv_precision.log_prob(precisions)

# check against the numpy implementation
precisions = np.stack([np.eye(2, dtype=np.float32), true_precision])
n = precisions.shape[0]
lik_tf = log_lik_prior(precisions).numpy()

for i in range(n):
  print(i)
  print('numpy:', log_lik_prior_numpy(precisions[i]))
  print('tensorflow:', lik_tf[i])

0
numpy: -2.2351873809649625
tensorflow: -2.2351875
1
numpy: -9.103606346649766
tensorflow: -9.103608

संयुक्त लॉग संभावना फ़ंक्शन बनाएं

ऊपर दिया गया डेटा लॉग संभावना फ़ंक्शन हमारी टिप्पणियों पर निर्भर करता है, लेकिन सैंपलर के पास वे नहीं होंगे। हम एक [क्लोजर] (https://en.wikipedia.org/wiki/Closure_(computer_programming) का उपयोग करके वैश्विक चर का उपयोग किए बिना निर्भरता से छुटकारा पा सकते हैं। क्लोजर में एक बाहरी फ़ंक्शन शामिल होता है जो एक ऐसे वातावरण का निर्माण करता है जिसमें एक के लिए आवश्यक चर होते हैं आंतरिक कार्य।

def get_log_lik(data, n_chains=1):
  # The data argument that is passed in will be available to the inner function
  # below so it doesn't have to be passed in as a parameter.
  replicated_data = np.tile(np.expand_dims(data, axis=1), reps=[1, n_chains, 1])

  @tf.function(autograph=False)
  def _log_lik(precision):
    return log_lik_data(precision, replicated_data) + log_lik_prior(precision)

  return _log_lik

चरण 4: नमूना

ठीक है, नमूना लेने का समय! चीजों को सरल रखने के लिए, हम केवल 1 श्रृंखला का उपयोग करेंगे और हम प्रारंभिक बिंदु के रूप में पहचान मैट्रिक्स का उपयोग करेंगे। हम बाद में और सावधानी से काम करेंगे।

दोबारा, यह काम नहीं करेगा - हमें एक अपवाद मिलेगा।

@tf.function(autograph=False)
def sample():
  tf.random.set_seed(123)
  init_precision = tf.expand_dims(tf.eye(2), axis=0)

  # Use expand_dims because we want to pass in a tensor of starting values
  log_lik_fn = get_log_lik(my_data, n_chains=1)

  # we'll just do a few steps here
  num_results = 10
  num_burnin_steps = 10
  states = tfp.mcmc.sample_chain(
     num_results=num_results,
     num_burnin_steps=num_burnin_steps,
     current_state=[
         init_precision,
     ],
     kernel=tfp.mcmc.HamiltonianMonteCarlo(
         target_log_prob_fn=log_lik_fn,
         step_size=0.1,
         num_leapfrog_steps=3),
     trace_fn=None,
     seed=123)
  return states

try:
  states = sample()
except Exception as e:
  # shorten the giant stack trace
  lines = str(e).split('\n')
  print('\n'.join(lines[:5]+['...']+lines[-3:]))

Cholesky decomposition was not successful. The input might not be valid.
     [[{ {node mcmc_sample_chain/trace_scan/while/body/_79/smart_for_loop/while/body/_371/mh_one_step/hmc_kernel_one_step/leapfrog_integrate/while/body/_537/leapfrog_integrate_one_step/maybe_call_fn_and_grads/value_and_gradients/StatefulPartitionedCall/Cholesky} }]] [Op:__inference_sample_2849]

Function call stack:
sample
...
Function call stack:
sample

समस्या की पहचान

InvalidArgumentError (see above for traceback): Cholesky decomposition was not successful. The input might not be valid. यह सुपर मददगार नहीं है। आइए देखें कि क्या हुआ इसके बारे में हम और जान सकते हैं।

  • हम प्रत्येक चरण के लिए मापदंडों का प्रिंट आउट लेंगे ताकि हम वह मान देख सकें जिसके लिए चीजें विफल होती हैं
  • विशिष्ट समस्याओं से बचाव के लिए हम कुछ अभिकथन जोड़ेंगे।

दावे मुश्किल हैं क्योंकि वे TensorFlow संचालन हैं, और हमें ध्यान रखना होगा कि वे निष्पादित हो जाएं और ग्राफ़ से अनुकूलित न हों। यह की कीमत पढ़ने इस अवलोकन TensorFlow डीबगिंग के लिए यदि आप TF दावे से परिचित नहीं हैं। आप स्पष्ट रूप से उपयोग करते हुए निष्पादित करने के लिए मजबूर कर सकते हैं दावे tf.control_dependencies (नीचे दिए गए कोड में टिप्पणी देखें)।

TensorFlow के मूल Print समारोह दावे के रूप में ही व्यवहार किया है - यह एक ऑपरेशन है, और आप यह सुनिश्चित करें कि यह कार्यान्वित कुछ देखभाल करने के लिए की जरूरत है। Print जब हम एक नोटबुक में काम कर रहे हैं अतिरिक्त सिर दर्द का कारण बनता है: इसके उत्पादन के लिए भेजा जाता stderr , और stderr सेल में प्रदर्शित नहीं है। हम यहाँ एक चाल का उपयोग करेंगे: उपयोग करने के बजाय tf.Print , हम के माध्यम से हमारे अपने TensorFlow प्रिंट कार्रवाई पैदा हो जाएगी tf.pyfunc । दावे के साथ, हमें यह सुनिश्चित करना होगा कि हमारी विधि निष्पादित हो।

def get_log_lik_verbose(data, n_chains=1):
  # The data argument that is passed in will be available to the inner function
  # below so it doesn't have to be passed in as a parameter.
  replicated_data = np.tile(np.expand_dims(data, axis=1), reps=[1, n_chains, 1])

  def _log_lik(precisions):
    # An internal method we'll make into a TensorFlow operation via tf.py_func
    def _print_precisions(precisions):
      print('precisions:\n', precisions)
      return False  # operations must return something!
    # Turn our method into a TensorFlow operation
    print_op = tf.compat.v1.py_func(_print_precisions, [precisions], tf.bool)

    # Assertions are also operations, and some care needs to be taken to ensure
    # that they're executed
    assert_op = tf.assert_equal(
        precisions, tf.linalg.matrix_transpose(precisions),
        message='not symmetrical', summarize=4, name='symmetry_check')

    # The control_dependencies statement forces its arguments to be executed
    # before subsequent operations
    with tf.control_dependencies([print_op, assert_op]):
      return (log_lik_data(precisions, replicated_data) +
              log_lik_prior(precisions))

  return _log_lik

@tf.function(autograph=False)
def sample():
  tf.random.set_seed(123)
  init_precision = tf.eye(2)[tf.newaxis, ...]
  log_lik_fn = get_log_lik_verbose(my_data)
  # we'll just do a few steps here
  num_results = 10
  num_burnin_steps = 10
  states = tfp.mcmc.sample_chain(
      num_results=num_results,
      num_burnin_steps=num_burnin_steps,
      current_state=[
          init_precision,
      ],
      kernel=tfp.mcmc.HamiltonianMonteCarlo(
          target_log_prob_fn=log_lik_fn,
          step_size=0.1,
          num_leapfrog_steps=3),
      trace_fn=None,
      seed=123)

try:
  states = sample()
except Exception as e:
  # shorten the giant stack trace
  lines = str(e).split('\n')
  print('\n'.join(lines[:5]+['...']+lines[-3:]))

precisions:
 [[[1. 0.]
  [0. 1.]]]
precisions:
 [[[1. 0.]
  [0. 1.]]]
precisions:
 [[[ 0.24315196 -0.2761638 ]
  [-0.33882257  0.8622    ]]]
 assertion failed: [not symmetrical] [Condition x == y did not hold element-wise:] [x (leapfrog_integrate_one_step/add:0) = ] [[[0.243151963 -0.276163787][-0.338822573 0.8622]]] [y (leapfrog_integrate_one_step/maybe_call_fn_and_grads/value_and_gradients/matrix_transpose/transpose:0) = ] [[[0.243151963 -0.338822573][-0.276163787 0.8622]]]
     [[{ {node mcmc_sample_chain/trace_scan/while/body/_96/smart_for_loop/while/body/_381/mh_one_step/hmc_kernel_one_step/leapfrog_integrate/while/body/_503/leapfrog_integrate_one_step/maybe_call_fn_and_grads/value_and_gradients/symmetry_check_1/Assert/AssertGuard/else/_577/Assert} }]] [Op:__inference_sample_4837]

Function call stack:
sample
...
Function call stack:
sample

यह विफल क्यों है

नमूनाकर्ता द्वारा आजमाया गया पहला नया पैरामीटर मान एक विषम मैट्रिक्स है। इससे चोल्स्की अपघटन विफल हो जाता है, क्योंकि यह केवल सममित (और सकारात्मक निश्चित) मैट्रिक्स के लिए परिभाषित है।

यहां समस्या यह है कि हमारी रुचि का पैरामीटर एक सटीक मैट्रिक्स है, और सटीक मैट्रिक्स वास्तविक, सममित और सकारात्मक निश्चित होना चाहिए। सैंपलर को इस बाधा के बारे में कुछ भी पता नहीं है (संभावित रूप से ग्रेडिएंट के माध्यम से छोड़कर), इसलिए यह पूरी तरह से संभव है कि सैंपलर एक अमान्य मान का प्रस्ताव करेगा, जिससे अपवाद हो सकता है, खासकर यदि चरण का आकार बड़ा हो।

हैमिल्टनियन मोंटे कार्लो सैंपलर के साथ, हम बहुत छोटे चरण आकार का उपयोग करके समस्या को हल करने में सक्षम हो सकते हैं, क्योंकि ग्रेडिएंट को मापदंडों को अमान्य क्षेत्रों से दूर रखना चाहिए, लेकिन छोटे चरण आकार का मतलब धीमा अभिसरण है। मेट्रोपोलिस-हेस्टिंग्स नमूना के साथ, जो ग्रेडियेंट के बारे में कुछ भी नहीं जानता, हम बर्बाद हो गए हैं।

संस्करण 2: अप्रतिबंधित मापदंडों के लिए पुनरावृति करना

उपरोक्त समस्या का एक सीधा समाधान है: हम अपने मॉडल को इस तरह से पुन: व्यवस्थित कर सकते हैं कि नए मापदंडों में अब ये बाधाएँ न हों। ऐसा करने के लिए TFP टूल - बायजेक्टर - का एक उपयोगी सेट प्रदान करता है।

बायजेक्टर के साथ पुनर्मूल्यांकन

हमारा सटीक मैट्रिक्स वास्तविक और सममित होना चाहिए; हम एक वैकल्पिक पैरामीटरकरण चाहते हैं जिसमें ये बाधाएं न हों। एक प्रारंभिक बिंदु सटीक मैट्रिक्स का एक चोल्स्की गुणनखंड है। चोल्स्की कारक अभी भी विवश हैं - वे निचले त्रिकोणीय हैं, और उनके विकर्ण तत्व सकारात्मक होने चाहिए। हालांकि, अगर हम चोल्स्की कारक के विकर्णों का लॉग लेते हैं, तो लॉग अब सकारात्मक होने के लिए बाध्य नहीं हैं, और फिर यदि हम निचले त्रिकोणीय हिस्से को 1-डी वेक्टर में समतल करते हैं, तो हमारे पास अब कम त्रिकोणीय बाधा नहीं है . हमारे मामले में परिणाम बिना किसी बाधा के लंबाई 3 वेक्टर होगा।

( स्टेन मैनुअल परिवर्तनों का उपयोग कर मानकों के आधार पर बाधाओं के विभिन्न प्रकार दूर करने के लिए पर एक बड़ा अध्याय है।)

इस पुनरावर्तन का हमारे डेटा लॉग संभावना फ़ंक्शन पर बहुत कम प्रभाव पड़ता है - हमें बस अपने परिवर्तन को उल्टा करना होगा ताकि हम सटीक मैट्रिक्स वापस प्राप्त कर सकें - लेकिन पूर्व पर प्रभाव अधिक जटिल है। हमने निर्दिष्ट किया है कि दिए गए सटीक मैट्रिक्स की संभावना विशार्ट वितरण द्वारा दी गई है; हमारे रूपांतरित मैट्रिक्स की प्रायिकता क्या है?

याद रखें कि अगर हम एक monotonic समारोह लागू \(g\) एक 1-डी यादृच्छिक चर को \(X\), \(Y = g(X)\), के लिए घनत्व \(Y\) द्वारा दिया जाता है

\[ f_Y(y) = | \frac{d}{dy}(g^{-1}(y)) | f_X(g^{-1}(y)) \]

के व्युत्पन्न \(g^{-1}\) अवधि तरीका है कि के लिए खातों \(g\) स्थानीय संस्करणों बदल जाता है। उच्च आयामी यादृच्छिक परिवर्तनीय के लिए, सुधारात्मक कारक के Jacobian का निर्धारक का निरपेक्ष मान है \(g^{-1}\) (देखें यहाँ )।

हमें अपने लॉग पूर्व संभावना फ़ंक्शन में व्युत्क्रम परिवर्तन का एक जैकोबियन जोड़ना होगा। सौभाग्य से, TFP के Bijector वर्ग हमारे लिए इस की देखभाल कर सकते हैं।

Bijector वर्ग उलटी, चिकनी प्रायिकता घनत्व कार्यों में चर को बदलने के लिए इस्तेमाल किया कार्यों का प्रतिनिधित्व करने के लिए किया जाता है। Bijectors सब एक है forward() विधि है कि प्रदर्शन एक को बदलने, एक inverse() विधि है कि उलट यह, और forward_log_det_jacobian() और inverse_log_det_jacobian() तरीकों कि Jacobian सुधार हम की जरूरत है जब हम एक पीडीएफ reparaterize प्रदान करते हैं।

TFP उपयोगी bijectors है कि हम के माध्यम से रचना के माध्यम से गठजोड़ कर सकते हैं का एक संग्रह प्रदान करता है Chain ऑपरेटर काफी जटिल रूपांतरण के रूप में। हमारे मामले में, हम निम्नलिखित 3 बायजेक्टर लिखेंगे (श्रृंखला में संचालन दाएं से बाएं किया जाता है):

  1. हमारे परिवर्तन का पहला चरण सटीक मैट्रिक्स पर चोल्स्की फ़ैक्टराइज़ेशन करना है। उसके लिए कोई बिजेक्टर वर्ग नहीं है; हालांकि, CholeskyOuterProduct bijector 2 Cholesky कारकों में से उत्पाद लेता है। हम चाहते हैं कि आपरेशन का उपयोग करने का उलटा उपयोग कर सकते हैं Invert ऑपरेटर।
  2. अगला कदम चोल्स्की कारक के विकर्ण तत्वों का लॉग लेना है। हम के माध्यम से यह पूरा TransformDiagonal bijector और का प्रतिलोम Exp bijector।
  3. अंत में हम एक सदिश का प्रतिलोम का उपयोग करने के मैट्रिक्स के निचले त्रिकोणीय भाग समतल FillTriangular bijector।
# Our transform has 3 stages that we chain together via composition:
precision_to_unconstrained = tfb.Chain([
    # step 3: flatten the lower triangular portion of the matrix
    tfb.Invert(tfb.FillTriangular(validate_args=VALIDATE_ARGS)),
    # step 2: take the log of the diagonals    
    tfb.TransformDiagonal(tfb.Invert(tfb.Exp(validate_args=VALIDATE_ARGS))),
    # step 1: decompose the precision matrix into its Cholesky factors
    tfb.Invert(tfb.CholeskyOuterProduct(validate_args=VALIDATE_ARGS)),
])

# sanity checks
m = tf.constant([[1., 2.], [2., 8.]])
m_fwd = precision_to_unconstrained.forward(m)
m_inv = precision_to_unconstrained.inverse(m_fwd)

# bijectors handle tensors of values, too!
m2 = tf.stack([m, tf.eye(2)])
m2_fwd = precision_to_unconstrained.forward(m2)
m2_inv = precision_to_unconstrained.inverse(m2_fwd)

print('single input:')
print('m:\n', m.numpy())
print('precision_to_unconstrained(m):\n', m_fwd.numpy())
print('inverse(precision_to_unconstrained(m)):\n', m_inv.numpy())
print()

print('tensor of inputs:')
print('m2:\n', m2.numpy())
print('precision_to_unconstrained(m2):\n', m2_fwd.numpy())
print('inverse(precision_to_unconstrained(m2)):\n', m2_inv.numpy())

single input:
m:
 [[1. 2.]
 [2. 8.]]
precision_to_unconstrained(m):
 [0.6931472 2.        0.       ]
inverse(precision_to_unconstrained(m)):
 [[1. 2.]
 [2. 8.]]

tensor of inputs:
m2:
 [[[1. 2.]
  [2. 8.]]

 [[1. 0.]
  [0. 1.]]]
precision_to_unconstrained(m2):
 [[0.6931472 2.        0.       ]
 [0.        0.        0.       ]]
inverse(precision_to_unconstrained(m2)):
 [[[1. 2.]
  [2. 8.]]

 [[1. 0.]
  [0. 1.]]]

TransformedDistribution वर्ग के लिए आवश्यक Jacobian सुधार एक वितरण करने के लिए एक bijector लागू करने और बनाने की प्रक्रिया को स्वचालित log_prob() हमारा नया पूर्व बन जाता है:

def log_lik_prior_transformed(transformed_precisions):
  rv_precision = tfd.TransformedDistribution(
      tfd.WishartTriL(
          df=PRIOR_DF,
          scale_tril=tf.linalg.cholesky(PRIOR_SCALE),
          validate_args=VALIDATE_ARGS,
          allow_nan_stats=ALLOW_NAN_STATS),
      bijector=precision_to_unconstrained,
      validate_args=VALIDATE_ARGS)
  return rv_precision.log_prob(transformed_precisions)

# Check against the numpy implementation.  Note that when comparing, we need
# to add in the Jacobian correction.
precisions = np.stack([np.eye(2, dtype=np.float32), true_precision])
transformed_precisions = precision_to_unconstrained.forward(precisions)
lik_tf = log_lik_prior_transformed(transformed_precisions).numpy()
corrections = precision_to_unconstrained.inverse_log_det_jacobian(
    transformed_precisions, event_ndims=1).numpy()
n = precisions.shape[0]

for i in range(n):
  print(i)
  print('numpy:', log_lik_prior_numpy(precisions[i]) + corrections[i])
  print('tensorflow:', lik_tf[i])

0
numpy: -0.8488930160357633
tensorflow: -0.84889317
1
numpy: -7.305657151741624
tensorflow: -7.305659

हमें अपने डेटा लॉग संभावना के लिए केवल परिवर्तन को उल्टा करने की आवश्यकता है:

precision = precision_to_unconstrained.inverse(transformed_precision)

चूंकि हम वास्तव में सटीक मैट्रिक्स के चोल्स्की कारककरण चाहते हैं, इसलिए यहां केवल आंशिक उलटा करना अधिक कुशल होगा। हालांकि, हम बाद के लिए अनुकूलन छोड़ देंगे और पाठक के लिए एक अभ्यास के रूप में आंशिक उलटा छोड़ देंगे।

def log_lik_data_transformed(transformed_precisions, replicated_data):
  # We recover the precision matrix by inverting our bijector.  This is
  # inefficient since we really want the Cholesky decomposition of the
  # precision matrix, and the bijector has that in hand during the inversion,
  # but we'll worry about efficiency later.
  n = tf.shape(transformed_precisions)[0]
  precisions = precision_to_unconstrained.inverse(transformed_precisions)
  precisions_cholesky = tf.linalg.cholesky(precisions)
  covariances = tf.linalg.cholesky_solve(
      precisions_cholesky, tf.linalg.eye(2, batch_shape=[n]))
  rv_data = tfd.MultivariateNormalFullCovariance(
      loc=tf.zeros([n, 2]),
      covariance_matrix=covariances,
      validate_args=VALIDATE_ARGS,
      allow_nan_stats=ALLOW_NAN_STATS)

  return tf.reduce_sum(rv_data.log_prob(replicated_data), axis=0)

# sanity check
precisions = np.stack([np.eye(2, dtype=np.float32), true_precision])
transformed_precisions = precision_to_unconstrained.forward(precisions)
lik_tf = log_lik_data_transformed(
    transformed_precisions, replicated_data).numpy()

for i in range(precisions.shape[0]):
  print(i)
  print('numpy:', log_lik_data_numpy(precisions[i], my_data))
  print('tensorflow:', lik_tf[i])

0
numpy: -430.71218815801365
tensorflow: -430.71207
1
numpy: -280.81822950593767
tensorflow: -280.8182

फिर से हम अपने नए कार्यों को बंद कर देते हैं।

def get_log_lik_transformed(data, n_chains=1):
  # The data argument that is passed in will be available to the inner function
  # below so it doesn't have to be passed in as a parameter.
  replicated_data = np.tile(np.expand_dims(data, axis=1), reps=[1, n_chains, 1])

  @tf.function(autograph=False)
  def _log_lik_transformed(transformed_precisions):
    return (log_lik_data_transformed(transformed_precisions, replicated_data) +
            log_lik_prior_transformed(transformed_precisions))

  return _log_lik_transformed

# make sure everything runs
log_lik_fn = get_log_lik_transformed(my_data)
m = tf.eye(2)[tf.newaxis, ...]
lik = log_lik_fn(precision_to_unconstrained.forward(m)).numpy()
print(lik)

[-431.5611]

सैम्पलिंग

अब जब हमें अमान्य पैरामीटर मानों के कारण हमारे सैंपलर के उड़ने की चिंता करने की आवश्यकता नहीं है, तो चलिए कुछ वास्तविक नमूने उत्पन्न करते हैं।

नमूना हमारे मापदंडों के अप्रतिबंधित संस्करण के साथ काम करता है, इसलिए हमें अपने प्रारंभिक मूल्य को इसके अप्रतिबंधित संस्करण में बदलने की आवश्यकता है। हम जो नमूने उत्पन्न करते हैं, वे भी सभी अपने अप्रतिबंधित रूप में होंगे, इसलिए हमें उन्हें वापस बदलने की आवश्यकता है। बिजेक्टर वेक्टरकृत होते हैं, इसलिए ऐसा करना आसान है।

# We'll choose a proper random initial value this time
np.random.seed(123)
initial_value_cholesky = np.array(
    [[0.5 + np.random.uniform(), 0.0],
     [-0.5 + np.random.uniform(), 0.5 + np.random.uniform()]],
    dtype=np.float32)
initial_value =  initial_value_cholesky.dot(
  initial_value_cholesky.T)[np.newaxis, ...]

# The sampler works with unconstrained values, so we'll transform our initial
# value
initial_value_transformed = precision_to_unconstrained.forward(
  initial_value).numpy()

# Sample!
@tf.function(autograph=False)
def sample():
  tf.random.set_seed(123)
  log_lik_fn = get_log_lik_transformed(my_data, n_chains=1)

  num_results = 1000
  num_burnin_steps = 1000

  states, is_accepted = tfp.mcmc.sample_chain(
      num_results=num_results,
      num_burnin_steps=num_burnin_steps,
      current_state=[
          initial_value_transformed,
      ],
      kernel=tfp.mcmc.HamiltonianMonteCarlo(
          target_log_prob_fn=log_lik_fn,
          step_size=0.1,
          num_leapfrog_steps=3),
      trace_fn=lambda _, pkr: pkr.is_accepted,
      seed=123)
  # transform samples back to their constrained form
  precision_samples = [precision_to_unconstrained.inverse(s) for s in states]
  return states, precision_samples, is_accepted

states, precision_samples, is_accepted = sample()

आइए हमारे सैंपलर के आउटपुट के माध्य की तुलना विश्लेषणात्मक पश्च माध्य से करें!

print('True posterior mean:\n', posterior_mean)
print('Sample mean:\n', np.mean(np.reshape(precision_samples, [-1, 2, 2]), axis=0))

True posterior mean:
 [[ 0.9641779 -1.6534661]
 [-1.6534661  3.8683164]]
Sample mean:
 [[ 1.4315274  -0.25587553]
 [-0.25587553  0.5740424 ]]

हम बहुत दूर हैं! आइए जानें क्यों। आइए पहले हमारे नमूने देखें।

np.reshape(precision_samples, [-1, 2, 2])

array([[[ 1.4315385, -0.2558777],
        [-0.2558777,  0.5740494]],

       [[ 1.4315385, -0.2558777],
        [-0.2558777,  0.5740494]],

       [[ 1.4315385, -0.2558777],
        [-0.2558777,  0.5740494]],

       ...,

       [[ 1.4315385, -0.2558777],
        [-0.2558777,  0.5740494]],

       [[ 1.4315385, -0.2558777],
        [-0.2558777,  0.5740494]],

       [[ 1.4315385, -0.2558777],
        [-0.2558777,  0.5740494]]], dtype=float32)

उह ओह - ऐसा लगता है कि उन सभी का मूल्य समान है। आइए जानें क्यों।

kernel_results_ चर नाम टपल कि प्रत्येक राज्य में नमूना बारे में जानकारी देता है। is_accepted क्षेत्र यहाँ मुख्य बिंदु है।

# Look at the acceptance for the last 100 samples
print(np.squeeze(is_accepted)[-100:])
print('Fraction of samples accepted:', np.mean(np.squeeze(is_accepted)))

[False False False False False False False False False False False False
 False False False False False False False False False False False False
 False False False False False False False False False False False False
 False False False False False False False False False False False False
 False False False False False False False False False False False False
 False False False False False False False False False False False False
 False False False False False False False False False False False False
 False False False False False False False False False False False False
 False False False False]
Fraction of samples accepted: 0.0

हमारे सभी नमूने खारिज कर दिए गए थे! संभवत: हमारे कदमों का आकार बहुत बड़ा था। मैं चुना stepsize=0.1 विशुद्ध रूप से मनमाने ढंग से।

संस्करण 3: अनुकूली चरण आकार के साथ नमूना लेना

चूंकि चरण आकार की मेरी मनमानी पसंद के साथ नमूनाकरण विफल रहा, हमारे पास कुछ एजेंडा आइटम हैं:

  1. एक अनुकूली चरण आकार लागू करें, और
  2. कुछ अभिसरण जाँच करें।

वहाँ में कुछ अच्छा नमूना कोड है tensorflow_probability/python/mcmc/hmc.py अनुकूली कदम आकार लागू करने के लिए। मैंने इसे नीचे अनुकूलित किया है।

नोट एक अलग नहीं है कि sess.run() प्रत्येक चरण के लिए बयान। यह डिबगिंग के लिए वास्तव में सहायक है, क्योंकि यह हमें जरूरत पड़ने पर आसानी से कुछ प्रति-चरण निदान जोड़ने की अनुमति देता है। उदाहरण के लिए, हम वृद्धिशील प्रगति, प्रत्येक चरण का समय आदि दिखा सकते हैं।

युक्ति: अपने नमूना ऊपर गड़बड़ करने के लिए एक जाहिरा तौर पर आम तरीका पाश में अपने ग्राफ बढ़ने के लिए है। (सत्र चलने से पहले ग्राफ़ को अंतिम रूप देने का कारण केवल ऐसी समस्याओं को रोकना है।) यदि आप अंतिम रूप () का उपयोग नहीं कर रहे हैं, तो एक उपयोगी डिबगिंग जांच यदि आपका कोड क्रॉल में धीमा हो जाता है तो ग्राफ़ का प्रिंट आउट लेना है के माध्यम से हर कदम पर आकार len(mygraph.get_operations()) - अगर लंबाई बढ़ जाती है, तो आप शायद कुछ बुरा कर रहे हैं।

हम यहां 3 इंडिपेंडेंट चेन चलाने जा रहे हैं। जंजीरों के बीच कुछ तुलना करने से हमें अभिसरण की जाँच करने में मदद मिलेगी।

# The number of chains is determined by the shape of the initial values.
# Here we'll generate 3 chains, so we'll need a tensor of 3 initial values.
N_CHAINS = 3

np.random.seed(123)

initial_values = []
for i in range(N_CHAINS):
  initial_value_cholesky = np.array(
      [[0.5 + np.random.uniform(), 0.0],
       [-0.5 + np.random.uniform(), 0.5 + np.random.uniform()]],
      dtype=np.float32)
  initial_values.append(initial_value_cholesky.dot(initial_value_cholesky.T))
initial_values = np.stack(initial_values)

initial_values_transformed = precision_to_unconstrained.forward(
  initial_values).numpy()

@tf.function(autograph=False)
def sample():
  tf.random.set_seed(123)
  log_lik_fn = get_log_lik_transformed(my_data)

  # Tuning acceptance rates:
  dtype = np.float32
  num_burnin_iter = 3000
  num_warmup_iter = int(0.8 * num_burnin_iter) 
  num_chain_iter = 2500

  # Set the target average acceptance ratio for the HMC as suggested by
  # Beskos et al. (2013):
  # https://projecteuclid.org/download/pdfview_1/euclid.bj/1383661192
  target_accept_rate = 0.651

  # Initialize the HMC sampler.
  hmc = tfp.mcmc.HamiltonianMonteCarlo(
      target_log_prob_fn=log_lik_fn,
      step_size=0.01,
      num_leapfrog_steps=3)

  # Adapt the step size using standard adaptive MCMC procedure. See Section 4.2
  # of Andrieu and Thoms (2008):
  # http://www4.ncsu.edu/~rsmith/MA797V_S12/Andrieu08_AdaptiveMCMC_Tutorial.pdf
  adapted_kernel = tfp.mcmc.SimpleStepSizeAdaptation(
      inner_kernel=hmc,
      num_adaptation_steps=num_warmup_iter,
      target_accept_prob=target_accept_rate)

  states, is_accepted = tfp.mcmc.sample_chain(
      num_results=num_chain_iter,
      num_burnin_steps=num_burnin_iter,
      current_state=initial_values_transformed,
      kernel=adapted_kernel,
      trace_fn=lambda _, pkr: pkr.inner_results.is_accepted,
      parallel_iterations=1)
  # transform samples back to their constrained form
  precision_samples = precision_to_unconstrained.inverse(states)
  return states, precision_samples, is_accepted

states, precision_samples, is_accepted = sample()

एक त्वरित जांच: हमारे नमूने के दौरान हमारी स्वीकृति दर 0.651 के हमारे लक्ष्य के करीब है।

print(np.mean(is_accepted))

0.6190666666666667

इससे भी बेहतर, हमारे नमूना माध्य और मानक विचलन विश्लेषणात्मक समाधान से हम जो अपेक्षा करते हैं, उसके करीब हैं।

precision_samples_reshaped = np.reshape(precision_samples, [-1, 2, 2])
print('True posterior mean:\n', posterior_mean)
print('Mean of samples:\n', np.mean(precision_samples_reshaped, axis=0))

True posterior mean:
 [[ 0.9641779 -1.6534661]
 [-1.6534661  3.8683164]]
Mean of samples:
 [[ 0.96426415 -1.6519215 ]
 [-1.6519215   3.8614824 ]]

print('True posterior standard deviation:\n', posterior_sd)
print('Standard deviation of samples:\n', np.std(precision_samples_reshaped, axis=0))

True posterior standard deviation:
 [[0.13435492 0.25050813]
 [0.25050813 0.53903675]]
Standard deviation of samples:
 [[0.13622096 0.25235635]
 [0.25235635 0.5394968 ]]

अभिसरण की जाँच

सामान्य तौर पर हमारे पास जांच के लिए कोई विश्लेषणात्मक समाधान नहीं होगा, इसलिए हमें यह सुनिश्चित करने की आवश्यकता होगी कि नमूना अभिसरण हो गया है। एक मानक की जांच Gelman-रुबिन है \(\hat{R}\) आंकड़ा, जो कई नमूने चेन की आवश्यकता है। \(\hat{R}\) उपाय हैं जो विचरण (साधन के) श्रृंखलाओं के बीच करने के लिए डिग्री से अधिक है, तो चेन हूबहू वितरित किए गए एक उम्मीद होती है क्या। के मान \(\hat{R}\) 1 के करीब अनुमानित अभिसरण इंगित करने के लिए उपयोग किया जाता है। देखें स्रोत जानकारी के लिए।

r_hat = tfp.mcmc.potential_scale_reduction(precision_samples).numpy()
print(r_hat)

[[1.0038308 1.0005717]
 [1.0005717 1.0006068]]

मॉडल आलोचना

यदि हमारे पास विश्लेषणात्मक समाधान नहीं होता, तो यह कुछ वास्तविक मॉडल आलोचना करने का समय होता।

हमारी जमीनी सच्चाई (लाल रंग में) के सापेक्ष नमूना घटकों के कुछ त्वरित हिस्टोग्राम यहां दिए गए हैं। ध्यान दें कि नमूने नमूना सटीक मैट्रिक्स मानों से पहले पहचान मैट्रिक्स की ओर सिकुड़ गए हैं।

fig, axes = plt.subplots(2, 2, sharey=True)
fig.set_size_inches(8, 8)
for i in range(2):
  for j in range(2):
    ax = axes[i, j]
    ax.hist(precision_samples_reshaped[:, i, j])
    ax.axvline(true_precision[i, j], color='red',
               label='True precision')
    ax.axvline(sample_precision[i, j], color='red', linestyle=':',
               label='Sample precision')
    ax.set_title('precision[%d, %d]' % (i, j))
plt.tight_layout()
plt.legend()
plt.show()

पीएनजी

सटीक घटकों के जोड़े के कुछ स्कैटरप्लॉट बताते हैं कि पश्च की सहसंबंध संरचना के कारण, वास्तविक पश्च मान उतने असंभाव्य नहीं हैं जितने वे ऊपर के हाशिए से दिखाई देते हैं।

fig, axes = plt.subplots(4, 4)
fig.set_size_inches(12, 12)
for i1 in range(2):
  for j1 in range(2):
    index1 = 2 * i1 + j1
    for i2 in range(2):
      for j2 in range(2):
        index2 = 2 * i2 + j2
        ax = axes[index1, index2]
        ax.scatter(precision_samples_reshaped[:, i1, j1],
                   precision_samples_reshaped[:, i2, j2], alpha=0.1)
        ax.axvline(true_precision[i1, j1], color='red')
        ax.axhline(true_precision[i2, j2], color='red')
        ax.axvline(sample_precision[i1, j1], color='red', linestyle=':')
        ax.axhline(sample_precision[i2, j2], color='red', linestyle=':')
        ax.set_title('(%d, %d) vs (%d, %d)' % (i1, j1, i2, j2))
plt.tight_layout()
plt.show()

पीएनजी

संस्करण 4: विवश मापदंडों का सरल नमूना

बायजेक्टर्स ने सटीक मैट्रिक्स के नमूने को सीधा कर दिया, लेकिन अप्रतिबंधित प्रतिनिधित्व से और उसके लिए उचित मात्रा में मैन्युअल रूप से परिवर्तित किया गया था। एक आसान तरीका है!

ट्रांसफॉर्मेड ट्रांज़िशन कर्नेल

TransformedTransitionKernel इस प्रक्रिया को सरल। यह आपके नमूने को लपेटता है और सभी रूपांतरणों को संभालता है। यह एक तर्क के रूप में उन बायजेक्टरों की एक सूची लेता है जो अप्रतिबंधित पैरामीटर मानों को विवश लोगों के लिए मैप करते हैं। तो यहाँ हम का प्रतिलोम जरूरत precision_to_unconstrained bijector हम ऊपर का इस्तेमाल किया। हम सिर्फ इस्तेमाल कर सकते हैं tfb.Invert(precision_to_unconstrained) , लेकिन यह प्रतिलोम के प्रतिलोम के ले जा रहा शामिल होगा (TensorFlow सरल करने के लिए स्मार्ट पर्याप्त नहीं है tf.Invert(tf.Invert()) के लिए tf.Identity()) , तो बजाय हम बस एक नया बायजेक्टर लिखूंगा।

विवश बिजेक्टर

# The bijector we need for the TransformedTransitionKernel is the inverse of
# the one we used above
unconstrained_to_precision = tfb.Chain([
    # step 3: take the product of Cholesky factors
    tfb.CholeskyOuterProduct(validate_args=VALIDATE_ARGS),
    # step 2: exponentiate the diagonals    
    tfb.TransformDiagonal(tfb.Exp(validate_args=VALIDATE_ARGS)),
    # step 1: map a vector to a lower triangular matrix
    tfb.FillTriangular(validate_args=VALIDATE_ARGS),
])

# quick sanity check
m = [[1., 2.], [2., 8.]]
m_inv = unconstrained_to_precision.inverse(m).numpy()
m_fwd = unconstrained_to_precision.forward(m_inv).numpy()

print('m:\n', m)
print('unconstrained_to_precision.inverse(m):\n', m_inv)
print('forward(unconstrained_to_precision.inverse(m)):\n', m_fwd)

m:
 [[1.0, 2.0], [2.0, 8.0]]
unconstrained_to_precision.inverse(m):
 [0.6931472 2.        0.       ]
forward(unconstrained_to_precision.inverse(m)):
 [[1. 2.]
 [2. 8.]]

TransformedTransitionKernel के साथ नमूना लेना

साथ TransformedTransitionKernel , हम अब हमारे मानकों के मैनुअल परिवर्तनों क्या करना है। हमारे प्रारंभिक मूल्य और हमारे नमूने सभी सटीक मैट्रिक्स हैं; हमें बस अपने अनियंत्रित बायजेक्टर को कर्नेल में पास करना है और यह सभी परिवर्तनों का ध्यान रखता है।

@tf.function(autograph=False)
def sample():
  tf.random.set_seed(123)
  log_lik_fn = get_log_lik(my_data)

  # Tuning acceptance rates:
  dtype = np.float32
  num_burnin_iter = 3000
  num_warmup_iter = int(0.8 * num_burnin_iter) 
  num_chain_iter = 2500

  # Set the target average acceptance ratio for the HMC as suggested by
  # Beskos et al. (2013):
  # https://projecteuclid.org/download/pdfview_1/euclid.bj/1383661192
  target_accept_rate = 0.651

  # Initialize the HMC sampler.
  hmc = tfp.mcmc.HamiltonianMonteCarlo(
      target_log_prob_fn=log_lik_fn,
      step_size=0.01,
      num_leapfrog_steps=3)

  ttk = tfp.mcmc.TransformedTransitionKernel(
      inner_kernel=hmc, bijector=unconstrained_to_precision)

  # Adapt the step size using standard adaptive MCMC procedure. See Section 4.2
  # of Andrieu and Thoms (2008):
  # http://www4.ncsu.edu/~rsmith/MA797V_S12/Andrieu08_AdaptiveMCMC_Tutorial.pdf
  adapted_kernel = tfp.mcmc.SimpleStepSizeAdaptation(
      inner_kernel=ttk,
      num_adaptation_steps=num_warmup_iter,
      target_accept_prob=target_accept_rate)

  states = tfp.mcmc.sample_chain(
      num_results=num_chain_iter,
      num_burnin_steps=num_burnin_iter,
      current_state=initial_values,
      kernel=adapted_kernel,
      trace_fn=None,
      parallel_iterations=1)
  # transform samples back to their constrained form
  return states

precision_samples  = sample()

अभिसरण की जाँच

\(\hat{R}\) अभिसरण की जांच दिखता अच्छा!

r_hat = tfp.mcmc.potential_scale_reduction(precision_samples).numpy()
print(r_hat)

[[1.0013582 1.0019467]
 [1.0019467 1.0011805]]

विश्लेषणात्मक पोस्टीरियर के खिलाफ तुलना

आइए फिर से विश्लेषणात्मक पोस्टीरियर के खिलाफ जाँच करें।

# The output samples have shape [n_steps, n_chains, 2, 2]
# Flatten them to [n_steps * n_chains, 2, 2] via reshape:
precision_samples_reshaped = np.reshape(precision_samples, [-1, 2, 2])

print('True posterior mean:\n', posterior_mean)
print('Mean of samples:\n', np.mean(precision_samples_reshaped, axis=0))

True posterior mean:
 [[ 0.9641779 -1.6534661]
 [-1.6534661  3.8683164]]
Mean of samples:
 [[ 0.96687526 -1.6552585 ]
 [-1.6552585   3.867676  ]]

print('True posterior standard deviation:\n', posterior_sd)
print('Standard deviation of samples:\n', np.std(precision_samples_reshaped, axis=0))

True posterior standard deviation:
 [[0.13435492 0.25050813]
 [0.25050813 0.53903675]]
Standard deviation of samples:
 [[0.13329624 0.24913791]
 [0.24913791 0.53983927]]

अनुकूलन

अब जबकि हमारे पास चीजें एंड-टू-एंड चल रही हैं, आइए अधिक अनुकूलित संस्करण करते हैं। इस उदाहरण के लिए गति बहुत अधिक मायने नहीं रखती है, लेकिन एक बार जब मैट्रिसेस बड़े हो जाते हैं, तो कुछ अनुकूलन एक बड़ा अंतर लाएंगे।

एक बड़ा गति सुधार जो हम कर सकते हैं, वह है चोल्स्की अपघटन के संदर्भ में पुन: निर्धारित करना। इसका कारण यह है कि हमारे डेटा संभावना फ़ंक्शन के लिए कॉन्वर्सिस और सटीक मैट्रिसेस दोनों की आवश्यकता होती है। मैट्रिक्स उलट महंगा (है\(O(n^3)\) एक के लिए \(n \times n\) मैट्रिक्स), और अगर हम या तो सहप्रसरण या परिशुद्धता मैट्रिक्स के मामले में parameterize, हम अन्य प्राप्त करने के लिए एक व्युत्क्रम करने की ज़रूरत है।

आपको याद दिला दें, जैसा कि एक वास्तविक, सकारात्मक-निश्चित, सममित मैट्रिक्स \(M\) प्रपत्र का एक उत्पाद में विघटित किया जा सकता है \(M = L L^T\) जहां मैट्रिक्स \(L\) कम त्रिकोणीय है और सकारात्मक विकर्णों है। की Cholesky अपघटन को देखते हुए \(M\), हम और अधिक कुशलता से दोनों प्राप्त कर सकते हैं \(M\) (कम से उत्पाद और एक ऊपरी त्रिकोणीय मैट्रिक्स) और \(M^{-1}\) (बैक-प्रतिस्थापन के माध्यम से)। चोल्स्की गुणनखंड स्वयं गणना करने के लिए सस्ता नहीं है, लेकिन यदि हम चोल्स्की कारकों के संदर्भ में मापन करते हैं, तो हमें केवल प्रारंभिक पैरामीटर मानों के छोलेक्सी गुणनखंड की गणना करने की आवश्यकता है।

सहप्रसरण मैट्रिक्स के चोल्स्की अपघटन का उपयोग करना

TFP मल्टीवेरिएट सामान्य वितरण, का एक संस्करण है MultivariateNormalTriL , कि सहप्रसरण मैट्रिक्स के Cholesky कारक के मामले में parameterized है। इसलिए यदि हम सहप्रसरण मैट्रिक्स के चोल्स्की कारक के संदर्भ में मानकीकरण करते हैं, तो हम डेटा लॉग संभावना की कुशलता से गणना कर सकते हैं। चुनौती समान दक्षता के साथ पूर्व लॉग संभावना की गणना करने में है।

यदि हमारे पास प्रतिलोम विशार्ट वितरण का एक संस्करण होता जो नमूनों के चोल्स्की कारकों के साथ काम करता, तो हम पूरी तरह तैयार हो जाते। काश, हम नहीं। (हालांकि, टीम कोड सबमिशन का स्वागत करेगी!) एक विकल्प के रूप में, हम विशार्ट वितरण के एक संस्करण का उपयोग कर सकते हैं जो नमूने के चोल्स्की कारकों के साथ-साथ बायजेक्टर की एक श्रृंखला के साथ काम करता है।

फिलहाल, चीजों को वास्तव में कुशल बनाने के लिए हमारे पास कुछ स्टॉक बायजेक्टर नहीं हैं, लेकिन मैं इस प्रक्रिया को एक अभ्यास और टीएफपी के बायजेक्टर की शक्ति का एक उपयोगी उदाहरण के रूप में दिखाना चाहता हूं।

एक विशार्ट वितरण जो चोल्स्की कारकों पर संचालित होता है

Wishart वितरण एक उपयोगी झंडा है, input_output_cholesky कि निर्दिष्ट करता है कि इनपुट और आउटपुट मैट्रिक्स Cholesky कारकों होना चाहिए,। पूर्ण मैट्रिक्स की तुलना में चोल्स्की कारकों के साथ काम करना अधिक कुशल और संख्यात्मक रूप से फायदेमंद है, यही कारण है कि यह वांछनीय है। ध्वज के शब्दों के बारे में एक महत्वपूर्ण बात: यह केवल एक संकेत है कि इनपुट और वितरण के लिए उत्पादन का प्रतिनिधित्व बदलना चाहिए है - यह वितरण, जो करने के लिए एक Jacobian सुधार शामिल होगा की एक पूरी reparameterization का संकेत नहीं है log_prob() समारोह। हम वास्तव में यह पूर्ण पुनर्मूल्यांकन करना चाहते हैं, इसलिए हम अपना वितरण स्वयं बनाएंगे।

# An optimized Wishart distribution that has been transformed to operate on
# Cholesky factors instead of full matrices.  Note that we gain a modest
# additional speedup by specifying the Cholesky factor of the scale matrix
# (i.e. by passing in the scale_tril parameter instead of scale).

class CholeskyWishart(tfd.TransformedDistribution):
  """Wishart distribution reparameterized to use Cholesky factors."""
  def __init__(self,
      df,
      scale_tril,
      validate_args=False,
      allow_nan_stats=True,
      name='CholeskyWishart'):
    # Wishart has a bunch of methods that we want to support but not
    # implement.  We'll subclass TransformedDistribution here to take care of
    # those.  We'll override the few for which speed is critical and implement
    # them with a separate Wishart for which input_output_cholesky=True
    super(CholeskyWishart, self).__init__(
        distribution=tfd.WishartTriL(
            df=df,
            scale_tril=scale_tril,
            input_output_cholesky=False,
            validate_args=validate_args,
            allow_nan_stats=allow_nan_stats),
        bijector=tfb.Invert(tfb.CholeskyOuterProduct()),
        validate_args=validate_args,
        name=name
    )
    # Here's the Cholesky distribution we'll use for log_prob() and sample()
    self.cholesky = tfd.WishartTriL(
        df=df,
        scale_tril=scale_tril,
        input_output_cholesky=True,
        validate_args=validate_args,
        allow_nan_stats=allow_nan_stats)

  def _log_prob(self, x):
    return (self.cholesky.log_prob(x) +
            self.bijector.inverse_log_det_jacobian(x, event_ndims=2))

  def _sample_n(self, n, seed=None):
    return self.cholesky._sample_n(n, seed)

# some checks
PRIOR_SCALE_CHOLESKY = np.linalg.cholesky(PRIOR_SCALE)

@tf.function(autograph=False)
def compute_log_prob(m):
  w_transformed = tfd.TransformedDistribution(
      tfd.WishartTriL(df=PRIOR_DF, scale_tril=PRIOR_SCALE_CHOLESKY),
      bijector=tfb.Invert(tfb.CholeskyOuterProduct()))
  w_optimized = CholeskyWishart(
      df=PRIOR_DF, scale_tril=PRIOR_SCALE_CHOLESKY)
  log_prob_transformed = w_transformed.log_prob(m)
  log_prob_optimized = w_optimized.log_prob(m)
  return log_prob_transformed, log_prob_optimized

for matrix in [np.eye(2, dtype=np.float32),
               np.array([[1., 0.], [2., 8.]], dtype=np.float32)]:
  log_prob_transformed, log_prob_optimized = [
      t.numpy() for t in compute_log_prob(matrix)]
  print('Transformed Wishart:', log_prob_transformed)
  print('Optimized Wishart', log_prob_optimized)

Transformed Wishart: -0.84889317
Optimized Wishart -0.84889317
Transformed Wishart: -99.269455
Optimized Wishart -99.269455

उलटा विशार्ट वितरण का निर्माण

हम अपने सहप्रसरण मैट्रिक्स है \(C\) में विघटित \(C = L L^T\) जहां \(L\) कम त्रिकोणीय है और एक सकारात्मक विकर्ण है। हम की संभावना जानना चाहते \(L\) यह देखते हुए कि \(C \sim W^{-1}(\nu, V)\) जहां \(W^{-1}\) विशार्ट वितरण उल्टा होता है।

उलटा विशार्ट वितरण संपत्ति है कि अगर \(C \sim W^{-1}(\nu, V)\), तो सटीक मैट्रिक्स \(C^{-1} \sim W(\nu, V^{-1})\)। इसलिए हम की संभावना प्राप्त कर सकते हैं \(L\) एक के माध्यम से TransformedDistribution कि विशार्ट वितरण और एक bijector कि उसके व्युत्क्रम की एक Cholesky कारक के लिए सटीक मैट्रिक्स के Cholesky कारक नक्शे पैरामीटर के रूप में लेता है।

की Cholesky कारक से प्राप्त करने के लिए एक सरल (लेकिन सुपर कुशल नहीं) जिस तरह से \(C^{-1}\) को \(L\) वापस सुलझाने द्वारा Cholesky कारक को उलटने के लिए, तो इन उल्टे कारकों से सहप्रसरण मैट्रिक्स के गठन, और फिर एक Cholesky गुणन कर रही है .

की Cholesky अपघटन करते \(C^{-1} = M M^T\)। \(M\) तो हम का उपयोग कर इसे उलटने के कर सकते हैं, कम त्रिकोणीय है MatrixInverseTriL bijector।

बनाने \(C\) से \(M^{-1}\) : एक छोटे से मुश्किल है \(C = (M M^T)^{-1} = M^{-T}M^{-1} = M^{-T} (M^{-T})^T\)। \(M\) , कम त्रिकोणीय है तो \(M^{-1}\) भी कम त्रिकोणीय हो जाएगा, और \(M^{-T}\) ऊपरी त्रिकोणीय हो जाएगा। CholeskyOuterProduct() bijector केवल, कम त्रिकोणीय मैट्रिक्स के साथ काम करता तो हम इसका इस्तेमाल नहीं कर सकते हैं बनाने के लिए \(C\) से \(M^{-T}\)। हमारा वर्कअराउंड बायजेक्टर की एक श्रृंखला है जो एक मैट्रिक्स की पंक्तियों और स्तंभों को क्रमित करता है।

सौभाग्य से इस तर्क में समझाया गया है CholeskyToInvCholesky bijector!

सभी टुकड़ों को मिलाकर

# verify that the bijector works
m = np.array([[1., 0.], [2., 8.]], dtype=np.float32)
c_inv = m.dot(m.T)
c = np.linalg.inv(c_inv)
c_chol = np.linalg.cholesky(c)
wishart_cholesky_to_iw_cholesky = tfb.CholeskyToInvCholesky()
w_fwd = wishart_cholesky_to_iw_cholesky.forward(m).numpy()

print('numpy =\n', c_chol)
print('bijector =\n', w_fwd)

numpy =
 [[ 1.0307764   0.        ]
 [-0.03031695  0.12126781]]
bijector =
 [[ 1.0307764   0.        ]
 [-0.03031695  0.12126781]]

हमारा अंतिम वितरण

चोल्स्की कारकों पर काम कर रहा हमारा व्युत्क्रम विशार्ट इस प्रकार है:

inverse_wishart_cholesky = tfd.TransformedDistribution(
    distribution=CholeskyWishart(
        df=PRIOR_DF,
        scale_tril=np.linalg.cholesky(np.linalg.inv(PRIOR_SCALE))),
    bijector=tfb.CholeskyToInvCholesky())

हमें हमारा उलटा विशार्ट मिल गया है, लेकिन यह थोड़ा धीमा है क्योंकि हमें बायजेक्टर में चोल्स्की अपघटन करना है। आइए सटीक मैट्रिक्स पैरामीटराइजेशन पर वापस जाएं और देखें कि हम अनुकूलन के लिए वहां क्या कर सकते हैं।

अंतिम (!) संस्करण: सटीक मैट्रिक्स के चोल्स्की अपघटन का उपयोग करना

एक वैकल्पिक दृष्टिकोण सटीक मैट्रिक्स के चोल्स्की कारकों के साथ काम करना है। यहां पूर्व संभावना फ़ंक्शन की गणना करना आसान है, लेकिन डेटा लॉग संभावना फ़ंक्शन अधिक काम लेता है क्योंकि टीएफपी में बहुभिन्नरूपी सामान्य का एक संस्करण नहीं है जो कि सटीक द्वारा पैरामीटर किया गया है।

अनुकूलित पूर्व लॉग संभावना

हम का उपयोग CholeskyWishart वितरण हम ऊपर बनाया पहले के निर्माण के लिए।

# Our new prior.
PRIOR_SCALE_CHOLESKY = np.linalg.cholesky(PRIOR_SCALE)

def log_lik_prior_cholesky(precisions_cholesky):
  rv_precision = CholeskyWishart(
      df=PRIOR_DF,
      scale_tril=PRIOR_SCALE_CHOLESKY,
      validate_args=VALIDATE_ARGS,
      allow_nan_stats=ALLOW_NAN_STATS)
  return rv_precision.log_prob(precisions_cholesky)

# Check against the slower TF implementation and the NumPy implementation.
# Note that when comparing to NumPy, we need to add in the Jacobian correction.
precisions = [np.eye(2, dtype=np.float32),
              true_precision]
precisions_cholesky = np.stack([np.linalg.cholesky(m) for m in precisions])
precisions = np.stack(precisions)
lik_tf = log_lik_prior_cholesky(precisions_cholesky).numpy()
lik_tf_slow = tfd.TransformedDistribution(
    distribution=tfd.WishartTriL(
        df=PRIOR_DF, scale_tril=tf.linalg.cholesky(PRIOR_SCALE)),
    bijector=tfb.Invert(tfb.CholeskyOuterProduct())).log_prob(
    precisions_cholesky).numpy()
corrections = tfb.Invert(tfb.CholeskyOuterProduct()).inverse_log_det_jacobian(
    precisions_cholesky, event_ndims=2).numpy()
n = precisions.shape[0]

for i in range(n):
  print(i)
  print('numpy:', log_lik_prior_numpy(precisions[i]) + corrections[i])
  print('tensorflow slow:', lik_tf_slow[i])
  print('tensorflow fast:', lik_tf[i])

0
numpy: -0.8488930160357633
tensorflow slow: -0.84889317
tensorflow fast: -0.84889317
1
numpy: -7.442875031036973
tensorflow slow: -7.442877
tensorflow fast: -7.442876

अनुकूलित डेटा लॉग संभावना

We can use TFP's bijectors to build our own version of the multivariate normal. Here is the key idea:

Suppose I have a column vector \(X\) whose elements are iid samples of \(N(0, 1)\). We have \(\text{mean}(X) = 0\) and \(\text{cov}(X) = I\)

Now let \(Y = A X + b\). We have \(\text{mean}(Y) = b\) and \(\text{cov}(Y) = A A^T\)

Hence we can make vectors with mean \(b\) and covariance \(C\) using the affine transform \(Ax+b\) to vectors of iid standard Normal samples provided \(A A^T = C\). The Cholesky decomposition of \(C\) has the desired property. However, there are other solutions.

Let \(P = C^{-1}\) and let the Cholesky decomposition of \(P\) be \(B\), ie \(B B^T = P\). Now

\(P^{-1} = (B B^T)^{-1} = B^{-T} B^{-1} = B^{-T} (B^{-T})^T\)

So another way to get our desired mean and covariance is to use the affine transform \(Y=B^{-T}X + b\).

Our approach (courtesy of this notebook ):

  1. Use tfd.Independent() to combine a batch of 1-D Normal random variables into a single multi-dimensional random variable. The reinterpreted_batch_ndims parameter for Independent() specifies the number of batch dimensions that should be reinterpreted as event dimensions. In our case we create a 1-D batch of length 2 that we transform into a 1-D event of length 2, so reinterpreted_batch_ndims=1 .
  2. Apply a bijector to add the desired covariance: tfb.Invert(tfb.ScaleMatvecTriL(scale_tril=precision_cholesky, adjoint=True)) . Note that above we're multiplying our iid normal random variables by the transpose of the inverse of the Cholesky factor of the precision matrix \((B^{-T}X)\). The tfb.Invert takes care of inverting \(B\), and the adjoint=True flag performs the transpose.
  3. Apply a bijector to add the desired offset: tfb.Shift(shift=shift) Note that we have to do the shift as a separate step from the initial inverted affine transform because otherwise the inverted scale is applied to the shift (since the inverse of \(y=Ax+b\) is \(x=A^{-1}y - A^{-1}b\)).
class MVNPrecisionCholesky(tfd.TransformedDistribution):
  """Multivariate normal parameterized by loc and Cholesky precision matrix."""

  def __init__(self, loc, precision_cholesky, name=None):
    super(MVNPrecisionCholesky, self).__init__(
        distribution=tfd.Independent(
            tfd.Normal(loc=tf.zeros_like(loc),
                       scale=tf.ones_like(loc)),
            reinterpreted_batch_ndims=1),
        bijector=tfb.Chain([
            tfb.Shift(shift=loc),
            tfb.Invert(tfb.ScaleMatvecTriL(scale_tril=precision_cholesky,
                                  adjoint=True)),
        ]),
        name=name)

@tf.function(autograph=False)
def log_lik_data_cholesky(precisions_cholesky, replicated_data):
  n = tf.shape(precisions_cholesky)[0]  # number of precision matrices
  rv_data = MVNPrecisionCholesky(
      loc=tf.zeros([n, 2]),
      precision_cholesky=precisions_cholesky)
  return tf.reduce_sum(rv_data.log_prob(replicated_data), axis=0)

# check against the numpy implementation
true_precision_cholesky = np.linalg.cholesky(true_precision)
precisions = [np.eye(2, dtype=np.float32), true_precision]
precisions_cholesky = np.stack([np.linalg.cholesky(m) for m in precisions])
precisions = np.stack(precisions)
n = precisions_cholesky.shape[0]
replicated_data = np.tile(np.expand_dims(my_data, axis=1), reps=[1, 2, 1])
lik_tf = log_lik_data_cholesky(precisions_cholesky, replicated_data).numpy()

for i in range(n):
  print(i)
  print('numpy:', log_lik_data_numpy(precisions[i], my_data))
  print('tensorflow:', lik_tf[i])

0
numpy: -430.71218815801365
tensorflow: -430.71207
1
numpy: -280.81822950593767
tensorflow: -280.81824

Combined log likelihood function

Now we combine our prior and data log likelihood functions in a closure.

def get_log_lik_cholesky(data, n_chains=1):
  # The data argument that is passed in will be available to the inner function
  # below so it doesn't have to be passed in as a parameter.
  replicated_data = np.tile(np.expand_dims(data, axis=1), reps=[1, n_chains, 1])

  @tf.function(autograph=False)
  def _log_lik_cholesky(precisions_cholesky):
    return (log_lik_data_cholesky(precisions_cholesky, replicated_data) +
            log_lik_prior_cholesky(precisions_cholesky))

  return _log_lik_cholesky

Constraining bijector

Our samples are constrained to be valid Cholesky factors, which means they must be lower triangular matrices with positive diagonals. The TransformedTransitionKernel needs a bijector that maps unconstrained tensors to/from tensors with our desired constraints. We've removed the Cholesky decomposition from the bijector's inverse, which speeds things up.

unconstrained_to_precision_cholesky = tfb.Chain([
    # step 2: exponentiate the diagonals    
    tfb.TransformDiagonal(tfb.Exp(validate_args=VALIDATE_ARGS)),
    # step 1: expand the vector to a lower triangular matrix
    tfb.FillTriangular(validate_args=VALIDATE_ARGS),
])

# some checks
inv = unconstrained_to_precision_cholesky.inverse(precisions_cholesky).numpy()
fwd = unconstrained_to_precision_cholesky.forward(inv).numpy()
print('precisions_cholesky:\n', precisions_cholesky)
print('\ninv:\n', inv)
print('\nfwd(inv):\n', fwd)

precisions_cholesky:
 [[[ 1.         0.       ]
  [ 0.         1.       ]]

 [[ 1.1470785  0.       ]
  [-2.0647411  1.0000004]]]

inv:
 [[ 0.0000000e+00  0.0000000e+00  0.0000000e+00]
 [ 3.5762781e-07 -2.0647411e+00  1.3721828e-01]]

fwd(inv):
 [[[ 1.         0.       ]
  [ 0.         1.       ]]

 [[ 1.1470785  0.       ]
  [-2.0647411  1.0000004]]]

Initial values

We generate a tensor of initial values. We're working with Cholesky factors, so we generate some Cholesky factor initial values.

# The number of chains is determined by the shape of the initial values.
# Here we'll generate 3 chains, so we'll need a tensor of 3 initial values.
N_CHAINS = 3

np.random.seed(123)

initial_values_cholesky = []
for i in range(N_CHAINS):
  initial_values_cholesky.append(np.array(
      [[0.5 + np.random.uniform(), 0.0],
       [-0.5 + np.random.uniform(), 0.5 + np.random.uniform()]],
      dtype=np.float32))
initial_values_cholesky = np.stack(initial_values_cholesky)

Sampling

We sample N_CHAINS chains using the TransformedTransitionKernel .

@tf.function(autograph=False)
def sample():
  tf.random.set_seed(123)
  log_lik_fn = get_log_lik_cholesky(my_data)

  # Tuning acceptance rates:
  dtype = np.float32
  num_burnin_iter = 3000
  num_warmup_iter = int(0.8 * num_burnin_iter) 
  num_chain_iter = 2500

  # Set the target average acceptance ratio for the HMC as suggested by
  # Beskos et al. (2013):
  # https://projecteuclid.org/download/pdfview_1/euclid.bj/1383661192
  target_accept_rate = 0.651

  # Initialize the HMC sampler.
  hmc = tfp.mcmc.HamiltonianMonteCarlo(
      target_log_prob_fn=log_lik_fn,
      step_size=0.01,
      num_leapfrog_steps=3)

  ttk = tfp.mcmc.TransformedTransitionKernel(
      inner_kernel=hmc, bijector=unconstrained_to_precision_cholesky)

  # Adapt the step size using standard adaptive MCMC procedure. See Section 4.2
  # of Andrieu and Thoms (2008):
  # http://www4.ncsu.edu/~rsmith/MA797V_S12/Andrieu08_AdaptiveMCMC_Tutorial.pdf
  adapted_kernel = tfp.mcmc.SimpleStepSizeAdaptation(
      inner_kernel=ttk,
      num_adaptation_steps=num_warmup_iter,
      target_accept_prob=target_accept_rate)

  states = tfp.mcmc.sample_chain(
      num_results=num_chain_iter,
      num_burnin_steps=num_burnin_iter,
      current_state=initial_values,
      kernel=adapted_kernel,
      trace_fn=None,
      parallel_iterations=1)
  # transform samples back to their constrained form
  samples = tf.linalg.matmul(states, states, transpose_b=True)
  return samples

precision_samples = sample()

Convergence check

A quick convergence check looks good:

r_hat = tfp.mcmc.potential_scale_reduction(precision_samples).numpy()
print(r_hat)

[[1.0013583 1.0019467]
 [1.0019467 1.0011804]]

Comparing results to the analytic posterior

# The output samples have shape [n_steps, n_chains, 2, 2]
# Flatten them to [n_steps * n_chains, 2, 2] via reshape:
precision_samples_reshaped = np.reshape(precision_samples, newshape=[-1, 2, 2])

And again, the sample means and standard deviations match those of the analytic posterior.

print('True posterior mean:\n', posterior_mean)
print('Mean of samples:\n', np.mean(precision_samples_reshaped, axis=0))

True posterior mean:
 [[ 0.9641779 -1.6534661]
 [-1.6534661  3.8683164]]
Mean of samples:
 [[ 0.9668749 -1.6552604]
 [-1.6552604  3.8676758]]

print('True posterior standard deviation:\n', posterior_sd)
print('Standard deviation of samples:\n', np.std(precision_samples_reshaped, axis=0))

True posterior standard deviation:
 [[0.13435492 0.25050813]
 [0.25050813 0.53903675]]
Standard deviation of samples:
 [[0.13329637 0.24913797]
 [0.24913797 0.53983945]]

Ok, all done! We've got our optimized sampler working.