Simpson’s paradox

I learned about the Simpson’s paradox fairly recently, and I found it quite disturbing, not because of the mere “paradox” itself, but mainly because I felt it was something I should have known already.

In case you haven’t heard about it, one instance of the paradox is a real-world medical study for comparing the success rate of two treatments for kidney stones (from Wikipedia):


Overall, Treatment B is better because its success rate is 83%, compared to 78% of Treatment A. However, when they split the patients into 2 groups: those with small stones and those with large stones, then Treatment A is better than Treatment B in both subgroups. Paradoxical enough?

Well, it’s not. It turns out that for severe cases (large stones), doctors tend to give the more effective Treatment A, while for milder cases with small stones, they tend to give the inferior Treatment B. Therefore the sum is dominant by group 2 and group 3, while the other groups contribute little to the final sums. So the results can be interpreted more accurately as: when Treatment B is more frequently applied to less severe cases, it can appear to be more effective.

Now, knowing that Treatment and Stone size are not independent, this should not come up as a paradox. In fact, we can visualize the problem as a graphical model like this

Untitled drawing

All the numbers in the table above can be expressed as conditional probabilities like so:

  • Group 1: p\left(S=true \vert T=A, St=small\right) = 0.93
  • Group 2: p\left(S=true \vert T=B, St=small\right) = 0.87
  • Group 3: p\left(S=true \vert T=A, St=large\right) = 0.73
  • Group 4: p\left(S=true \vert T=B, St=large\right) = 0.69
  • p\left(S=true \vert T=A\right) = 0.78
  • p\left(S=true \vert T=B\right) = 0.83

For any of us who studied Probability, it is no surprise that the probabilities might turn up-side-down whenever some conditional variables are stripped out of the equations. In this particular case, since S depends on both St and T, the last 2 equations do not bring any new knowledge about S.

So what is this “paradox” about? Isn’t it nothing more than the problem of confounding/lurking variables, something that most people in Probability/Statistics already known? In this particular case, Stone size is the lurking variable that dictates both Treatment and Success, therefore the scientists who designed the experiment should have taken it into account since the beginning. It is well-known among Statistic practitioners that they must try their best to identify and eliminate the effect of any lurking variables in their experiments, or at least keep them fixed, before drawing any meaningful conclusion.

From a slightly different perspective, the paradox can be understood once we understand the human bias of drawing causal relations. Human, perhaps for the sake of survival, constantly look for causal relations and often tend to ignore rates or proportions. Once we conceived something as being causal (Treatment B gives higher success rate than Treatment A in general), which might be wrong, we continue to assume a causal relation and proceed with that assumption in mind. Obviously with this assumption, we will find the success rates for the subgroups of patients to be highly counter-intuitive, or even paradoxical.

In fact, the connection of this paradox to human intuitions is so important that Judea Pearl dedicated a whole section in his book for it. Modern Statistical textbooks and curriculum, however, don’t even mention it. Instead they will generally present the topic along with lurking/confounding variables.

Therefore, if you haven’t heard about this, it is probably for a good reason, or perhaps you are simply too young.


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


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 *, 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
        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
            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.linalg.solve(l, self.Y))
            cov_posterior = k_2star + self.noise_variance * np.eye(k_2star.shape[0]) -,
            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 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.


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


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