Statistics

Monte Carlo và Markov Chain Monte Carlo

The devil is in the details, hoá ra ứng dụng Bayesian Optimization trong thực tế khá là phức tạp. Vậy nên bài này sẽ ôn lại kiến thức về phương pháp Monte Carlo và MCMC. Vẫn chưa biết đến khi nào mới có thể viết một bài tổng hợp về Bayesian Optimization.

1. Luật mạnh số lớn (Strong law of large numbers – SLLN)

Luật mạnh số lớn rất đơn giản: khi N càng lớn thì trung bình mẫu hội tụ gần như chắc chắn (almost surely converge) về kì vọng:

\displaystyle \frac{1}{N}\sum_{i=1}^N x_i \overset{a.s.}{\rightarrow} \mu khi N \rightarrow \infty

với điều kiện là các mẫu x_i phải là i.i.d.

Mặt khác, ta biết rằng nếu biến ngẫu nhiên X tuân theo phân phối xác suất f\left(x\right), và g\left(X\right) là một hàm theo X, thì kì vọng của g\left(X\right) được tính là:

\displaystyle E\left[g\left(X\right)\right] = \int g\left(x\right)f\left(x\right)dx

Vậy nên theo luật mạnh số lớn, khi số mẫu càng lớn thì giá trị trung bình mẫu của g\left(X\right) sẽ hội tụ về kì vọng:

\displaystyle \left(\frac{1}{N}\sum_{i=1}^N g\left(X_i\right)\right) \overset{a.s.}{\rightarrow} \left(\int g\left(x\right)f\left(x\right)dx\right) khi N \rightarrow \infty    (1)

Như vậy khi tích phân ở bên phải không tính được, ta có thể tính xấp xỉ nó bằng cách lấy thật nhiều mẫu X_i và tính biểu thức ở bên trái. Đây là tư tưởng cơ bản của phương pháp Monte Carlo, như sẽ trình bày sau đây.

2. Phương pháp Monte Carlo: Importance Sampling

Hầu hết các mô hình máy học có giám sát (supervised ML models) đều cố gắng mô hình hoá xác suất likelihood p\left(y\vert x; \theta\right), trong đó \theta là các tham số trong mô hình. Như trong bài trước, ta biết rằng trong cách tiếp cận Bayesian thì tham số \theta tuân theo phân phối tiền nghiệm p\left(\theta\right), và ta sẽ quan tâm đến xác suất hậu nghiệm (posterior):

\displaystyle p\left(\theta\vert x, y\right) = \frac{p\left(y\vert x, \theta\right)p\left(\theta\right)}{\int p\left(y\vert x, \theta\right)p\left(\theta\right)d\theta} = \frac{1}{Z}p\left(y\vert x, \theta\right)p\left(\theta\right)    (2)

Vấn đề là

Z = \displaystyle \int p\left(y\vert x, \theta\right)p\left(\theta\right)d\theta     (3)

là tích phân trên toàn miền xác định của \theta, và nói chung là khó chịu. Vậy nên ta không thể tính chính xác xác suất hậu nghiệm này được.

Cái mẹo chủ yếu của phương pháp Monte Carlo là ta sẽ nhân và chia công thức (3) cho một đại lượng q\left(\theta\right):

\displaystyle Z = \int p\left(y\vert x, \theta\right)p\left(\theta\right)d\theta = \int \frac{p\left(y \vert x, \theta\right)p\left(\theta\right)}{q\left(\theta\right)}q\left(\theta\right)d\theta =\int \omega\left(\theta\right) q\left(\theta\right)d\theta

Để ý rằng nếu q\left(\theta\right) là phân phối của \theta thì tích phân ở trên chính là kì vọng của \omega\left(\theta\right).

Bây giờ nếu ta lấy mẫu \theta^{\left(i\right)} \sim q\left(\theta\right), thì theo luật mạnh số lớn (1), khi N càng lớn thì trung bình mẫu hội tụ về kì vọng:

\displaystyle \left(\frac{1}{N}\sum_{i=1}^N \omega\left(\theta^{\left(i\right)}\right)\right) \overset{a.s.}{\rightarrow} \left(\int \omega\left(\theta\right) q\left(\theta\right)d\theta\right) khi N \rightarrow \infty    (4)

Với công thức này, ta có thể lấy mẫu để tính xác suất hậu nghiệm của tham số \theta trong (2).

Tuy nhiên trong thực tế ta hay phải làm prediction, tức là cho tập huấn luyện D = \left\{x_i, y_i\right\}_{i=1}^t và một mẫu input x_{t+1}, ta muốn mô hình hoá xác suất p\left(y_{t+1} \vert D, x_{t+1}\right). Ta viết như sau:

\displaystyle \begin{array}{rl} p\left(y_{t+1} \vert D, x_{t+1}\right) & \displaystyle = \int p\left(y_{t+1} \vert x_{t+1}, \theta\right) p\left(\theta \vert D \right)d\theta \\ & \displaystyle = \int p\left(y_{t+1} \vert x_{t+1}, \theta\right) \frac{p\left(D \vert \theta \right) p\left(\theta\right)}{p\left(D\right)}d\theta \\ & \displaystyle = \frac{1}{p\left(D\right)} \int p\left(y_{t+1} \vert x_{t+1}, \theta \right) p\left(D \vert \theta\right) p\left(\theta\right) d\theta \\ & = \frac{\displaystyle \int p\left(y_{t+1} \vert x_{t+1}, \theta \right) p\left(D \vert \theta\right) p\left(\theta\right) \frac{q\left(\theta\right)}{q\left(\theta\right)} d\theta}{\displaystyle \int p\left(D \vert \theta\right) p\left(\theta\right) \frac{q\left(\theta\right)}{q\left(\theta\right)} d\theta } \\ & = \frac{\displaystyle \int p\left(y_{t+1} \vert x_{t+1}, \theta \right) \omega\left(\theta\right) q\left(\theta\right) d\theta}{\displaystyle \int \omega\left(\theta\right) q\left(\theta\right) d\theta } \end{array}

Trong đó \omega\left(\theta\right) = \frac{p\left(D \vert \theta \right)p\left(\theta\right)}{q\left(\theta\right)}.

Áp dụng (4) ta có:

\displaystyle \begin{array}{rl} p\left(y_{t+1} \vert D, x_{t+1}\right) & = \frac{\displaystyle \int p\left(y_{t+1} \vert x_{t+1}, \theta \right) \omega\left(\theta\right) q\left(\theta\right) d\theta}{\displaystyle \int \omega\left(\theta\right) q\left(\theta\right) d\theta } \\ & \approx \frac{\displaystyle \frac{1}{N}\sum_i p\left(y_{t+1}\vert x_{t+1}, \theta^{\left(i\right)}\right)\omega\left(\theta^{\left(i\right)}\right)}{\displaystyle \frac{1}{N}\sum_j \omega\left(\theta^{\left(j\right)}\right)} \\ & \displaystyle \approx \sum_i \tilde{\omega}\left(\theta^{\left(i\right)}\right) p\left(y_{t+1}\vert x_{t+1}, \theta^{\left(i\right)}\right)\end{array}

với \tilde{\omega}\left(\theta^{\left(i\right)}\right) = \frac{\omega\left(\theta^{\left(i\right)}\right)}{\sum_j\omega\left(\theta^{\left(j\right)}\right)}.

q\left(\theta\right) là phân phối xác suất do ta “bịa” ra, ta sẽ chọn q\left(\theta\right) sao cho dễ lấy mẫu. Như vậy ta có thể lấy thật nhiều mẫu \theta_i, sau đó dùng công thức trên để xấp xỉ xác suất hậu nghiệm.

Điều tricky là để luật mạnh số lớn áp dụng được thì các mẫu \theta_i phải là i.i.d. Hơn nữa nếu q\left(\theta\right) càng gần p\left(D\vert \theta\right) p\left(\theta\right) thì chất lượng lấy mẫu sẽ càng cao.

3. Markov Chain Monte Carlo

Ý tưởng chính là ta sẽ “cơ cấu” một chuỗi Markov để lấy mẫu từ p\left(D\vert \theta\right) p\left(\theta\right), đồng thời bằng cách đi theo chuỗi Markov ta sẽ navigate qua tất cả các region quan trọng của phân phối này.

Hai thuật toán quan trọng của MCMC là Metropolis-Hasting và Gibbs sampling, sẽ được trình bày trong một bài khác, khi có dịp.

Textbook implementation of Gaussian Processes in Python

gp

So I was looking for a textbook implementation of Gaussian Processes, and disappointingly I couldn’t find any. Things like Spearmint or GPy are way complicated and very often implement some specific use-cases of GPs like regression and so on.

For educational purpose, here is a 10 line implementation of Gaussian Process:

<pre>import numpy as np


class Kernel(object):

    def compute(self, a, b):
        raise None


class SquaredDistanceKernel(Kernel):

    def __init__(self, kernel_param=0.1):
        self.kernel_parameter = kernel_param

    def compute(self, a, b):
        sq_dist = np.sum(a ** 2, 1).reshape(-1, 1) + np.sum(b ** 2, 1) - 2 * np.dot(a, b.T)
        return np.exp(-.5 * (1/self.kernel_parameter) * sq_dist)


class GaussianProcess(object):
    """
    Implements a GP with mean zero and a custom kernel
    """
    def __init__(self, kernel=SquaredDistanceKernel(), noise_variance=0.00005, x=None, y=None):
        """
        Initialize the GP with the given kernel and a noise parameter for the variance
        Optionally initialize this GP with given X and Y

        :param kernel: kernel function, has to be an instance of Kernel
        :param noise_variance:
        :param x: given input data
        :param y: given input label
        :return:
        """
        self.X = x
        self.Y = y
        self.kernel = kernel
        self.noise_variance = noise_variance
        self.cov = None if self.X is None else kernel.compute(self.X, self.X)

    def predict(self, x, y=None):
        """
        Given data in x, give the mean and covariance of the posterior predictive distribution p(f*|X*, X, f)
        If y is given, the function gives the predicts, as well as update the GP internally

        x should have size N x d1, y of size N x d2, where N is the number of samples

        :param x: the input data
        :param y: optional. If given, the GP parameters will be updated
        :return: a tuple (mu, cov, s):
            - mu: the mean of the posterior predictive distribution, of size N x d1
            - cov: the covariance matrix of the posterior predictive distribution, of size N x N
            - s: the standard deviation vector, convenient for plotting. Of size N x 1
        """
        # covariance of the new data
        k_2star = self.kernel.compute(x, x)

        if self.cov is None:
            # if there is no data in this GP, this is equivalent to the prior distribution (zero mean, unit covariance)
            mu = np.zeros(x.shape)
            cov_posterior = k_2star + (self.noise_variance * np.eye(k_2star.shape[0]))

            if y is not None:
                self.X = x
                self.Y = y
                self.cov = k_2star
        else:
            l = np.linalg.cholesky(self.cov + self.noise_variance * np.eye(self.cov.shape[0]))
            k_star = self.kernel.compute(self.X, x)
            l_div_k_star = np.linalg.solve(l, k_star)

            mu = np.dot(l_div_k_star.T, np.linalg.solve(l, self.Y))
            cov_posterior = k_2star + self.noise_variance * np.eye(k_2star.shape[0]) - np.dot(l_div_k_star.T,
                                                                                              l_div_k_star)
            if y is not None:
                self.X = np.vstack((self.X, x))
                self.Y = np.vstack((self.Y, y))
                self.cov = np.hstack((self.cov, k_star))
                self.cov = np.vstack((self.cov, np.hstack((k_star.T, k_2star))))

        return mu, cov_posterior, np.sqrt(np.diag(cov_posterior))</pre>

For educational purposes, I removed all the checks to make the code really concise.

As usual, the code can be found on github. In the same folder, you might find interative.py interesting.

In the next posts, I am going to explain all this.

Parameter Estimation techniques applied in Machine Learning

Bài này tóm tắt các kĩ thuật ước lượng tham số, vốn có nguồn gốc từ lí thuyết xác suất thống kê, đang được dùng rộng rãi trong Machine Learning. Bài này chỉ là ôn lại để chuẩn bị cho 1 series về Gaussian Process nên sẽ hơi sơ sài, không có ví dụ.

Trong thống kê, bài toán ước lượng tham số có thể tóm tắt như sau: Cho tập dữ liệu \left\{x_i\right\}_{i=1}^n, x_i \in \mathbf{R}^d, biết rằng tập dữ liệu này lấy mẫu từ phân phối f\left(x, \theta\right), với \theta là tham số.  Hãy ước lượng giá trị của tham số \theta.

Như vậy ta cần tập mẫu \left\{x_i\right\}_{i=1}^n, và biết trước dạng của hàm f\left(x, \theta\right).

1. Frequentist: Maximum Likelihood Estimation (MLE)

Theo trường phái frequentist, tham số \theta chỉ đơn giản là ẩn số (constant-valued but unknown), và MLE là một phương pháp thống kê để ước lượng ẩn số này.

Theo đó, hàm likelihood được định nghĩa là:

L\left(x_1, x_2, ..., x_n | \theta\right) = \prod_{i=1}^n f\left(x_i, \theta\right)

Ý tưởng là để \theta “phù hợp” với tập mẫu thì giá trị của L phải càng lớn càng tốt. Do đó bài toán ước lượng \theta chính là tìm \theta để L đạt cực đại:

\displaystyle \theta_{ML} = \arg\max_{\theta} \prod_{i=1}^n f\left(x_i, \theta\right)     (1)

Có nhiều cách để giải bài toán này. Nếu công thức của L là closed-form thì có thể giải rất đơn giản bằng cách tính đạo hàm và giải phương trình đạo hàm bằng 0.

Thông thường vì công thức của L là tích của nhiều đại lượng nhỏ, nên thay vì tìm cực đại của L, người ta tìm cực đại của log(L):

\displaystyle \log L = \sum_{i=1}^n\log f\left(x_i, \theta\right)

Và lời giải của MLE là nghiệm của phương trình:

\displaystyle \frac{\partial \log L}{\partial \theta} = \sum_{i=1}^n \left(\frac{1}{f\left(x_i, \theta\right)}\frac{\partial f\left(x_i, \theta\right)}{\partial \theta}\right) = 0

 Tuy nhiên trong thực tế, công thức (1) không closed-form thì bài toán được chuyển thành 1 bài toán optimization.

(more…)

Hóa ra mọi thứ đều là exponential!

Image

Hầu như tất cả các phân phối xác suất mình biết từ trước tới giờ, từ multinomial, Bernoulli cho tới Dirichlet… đều thuộc exponential family. Và dưới một góc nhìn chung như vậy, các khái niệm trở nên sáng tỏ hơn hẳn: từ sufficient statistics cho tới partition function hay conjugation…

Nhân tiện, cái hình lấy từ paper này.

Data science open class

Đây (http://cs109.org/) là lớp Data Science dạy cho sinh viên undergrad ở Harvard mấy tháng trước. Toàn bộ video lectures đã up và có thể xem miễn phí.
Xem video dạng này vẫn có vẻ thú vị hơn các lớp Data Science trên coursera.

RBM 3: Contrastive Divergence

Chuỗi bài viết về RBMs:
RBM 1: Hopfield nets and Boltzmann Machines
RBM 2: Model definition and training of RBMs
RBM 3: Contrastive Divergence
RBM 4: Contrastive Divergence in training RBMs and practical applications

Mặc dù chỉ bắt đầu “nổi tiếng” sau khi được dùng để huấn luyện Restricted Boltzmann Machines (RBMs) nhưng Contrastive Divergence thực ra đã được đề xuất từ trước đó, và bản thân thuật toán này cũng có vẻ đẹp theo cách riêng của nó. Do vậy trong bài này, ta tạm thời quên đi RBMs để xét thuật toán Contrastive Divergence từ góc độ tổng quan hơn. Cách áp dụng Contrastive Divergence trong huấn luyện RBMs sẽ được đề cập trong bài tiếp theo (và hi vọng là cuối cùng) của chuỗi bài này.

1. Contrastive Divergence là gì

Giả sử ta cần mô hình phân phối xác suất của một vector x bằng cách sử dụng một hàm có dạng f\left(x; \theta\right), trong đó \theta là vector chứa toàn bộ tham số của mô hình. Phân phối xác suất của x, kí hiệu p\left(x;\theta\right) phải có tổng bằng 1 nên:

\displaystyle p\left(x; \theta\right) = \frac{1}{Z\left(\theta\right)}f\left(x; \theta\right)       (1)

trong đó Z\left(\theta\right) gọi là partition function và được tính như sau:

\displaystyle Z\left(\theta\right) = \int f\left(x; \theta\right) dx        (2)

(more…)

How to assess the convergence of MCMC algorithms

1. Tổng quan

Markov Chain Monte Carlo (MCMC) là một lớp các thuật toán được dùng khá rộng rãi trong ML và thống kê. Trong ngữ cảnh ML thì:

  • đôi khi ta cần ước lượng một phân bố xác suất nào đó, nhưng lại không thể tính chính xác phân bố đó do công thức quá phức tạp (thông thường vì trong công thức có phép tích phân), và
  • Ta lại biết rằng một chuỗi Markov, khi chạy đủ lâu, thì sẽ hội tụ về một phân bố xác suất cố định, độc lập với trạng thái ban đầu của nó (đây là tính chất đẹp và rất dễ chứng minh, ta sẽ trở lại khi có dịp). Phân bố xác suất này gọi là stationary distribution.

Như vậy tư tưởng của MCMC là cho trước một phân bố xác suất cần ước lượng, ta sẽ xây dựng một chuỗi Markov sao cho stationary distribution của nó chính là phân bố xác suất đang xét. Xuất phát từ trạng thái ngẫu nhiên nào đó, chạy chuỗi Markov một thời gian đủ lâu để nó hội tụ về stationary distribution, thì khi đó mỗi trạng thái của chuỗi Markov có thể xem là được lấy mẫu từ phân bố xác suất đang cần ước lượng.

Ví dụ 1: Một ví dụ kinh điển là các thuật toán suy diễn của undirected graphical model (còn gọi là Markov network). Thông thường joint distribution của các biến trong Markov network được định nghĩa là

\displaystyle P\left(X_1,X_2,\ldots ,X_n\right)=\frac{1}{Z}\phi_1\left(X_{\phi_1}\right)\ldots\phi_k\left(X_{\phi_k}\right)

Trong đó Z thường được gọi là partition function, các \phi_k là potential thứ k thể hiện mối liên hệ giữa các biến nằm trong tập X_{\phi_{k}}\subset\left\{ X_{1},\ldots,X_{n}\right\}X_{\phi_k} thường được gọi là scope của \phi_k. Đặc biệt partition function Z định nghĩa là:

(more…)