Basics of Linear Regression

I will share with you about basics of Linear Regression. It's somehow entry point to statistics. However my ulterior motive is understanding of regularization. Actually, it's gonna be too long and intricate to cope with regularization here. Thus, in this article, I will handle only basics of linear regression. In terms of regularization, I will upload very soon :)

0. Simple linear regression¶

First of all, we have to wrap our head around "Simple linear regression".
Suppose x be predictor variable, y be dependent variable. Then linear regression line takes form

$$\hat{y} = \beta_0 + \beta_1x$$

Let experimental unit be $(x_i, y_i)\ (i = 1,2,\cdots,n)$, "Residential error" $S$ can be denoted as following.

$$y - \hat{y_i}$$

One of the way to obtain "best fitting" is to invoke "least squares criterion" which says "minimize the sum of the squared residual errors".
$$S = \sum_{i=1}^{n}\left\{y_i -(\beta_0 + \beta_1x)\right\}^2$$

As you can see, $S$ is gonna be quadratic function of the $\beta_0$ and $\beta_1$.
Thereby we can obtain "best fitting" by computing the followings.

$$\frac{\partial S}{\partial \beta_0} = 0$$$$\frac{\partial S}{\partial \beta_1} = 0$$

1. Multiple linear regression¶

Now let's getting into "Multiple linear regression". Let predictor variable be $(x_1, x_2, \cdots,x_d)$, dependent variable $y$ . Multiple linear regression line takes form

$$y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \cdots \beta_d x_d$$

where $( \beta_0,\beta_1, \cdots \beta_d )$ are called "regression coefficient".
Let's say predictor variable is $\vec{x}_t = (x_{t1}, x_{t2}, \cdots, x_{td})^T\ (t=1,2, \cdots, n)$, response variable is $\vec{y}_t\ (t=1,2, \cdots, n)$,

$$\hat{y_t} = \beta_0 + \beta_1x_{t1} + \beta_2x_{t2} + \cdots + \beta_dx_{td}$$

Now we'd like to denote this for $n$ experimental unit. Let $X$ be

$$X = \left(\begin{array}{cccc}1 & x_{ 11 } & \ldots & a_{ 1d } \\ 1 & x_{ 21 } & \ldots & a_{ 2d }\\ \vdots & \vdots & \ldots & \vdots\\ 1 & x_{ n1 } & \ldots & a_{ nd } \end{array}\right)$$

$\vec{\hat{y}}$ be

$$\vec{\hat{y}} = (y_1,y_2, \cdots, y_n)$$

$\vec{\beta}$ be

$$\vec{\beta} = (\beta_1, \beta_2, \cdots, \beta_d)$$

We can denote, $$\vec{\hat{y}} = X\vec{\beta}$$

Suppose $\epsilon$ is $\epsilon = (\epsilon_1, \epsilon_2, \cdots, \epsilon_n)^T$ as a "residual error". $$\vec{y} = X\vec{\beta} + \epsilon$$
Thereby,
$$\vec{y} - \vec{\hat{y}} = \epsilon$$ and
$$\epsilon_t = \vec{y_t} -\vec{\hat{y_t}}$$
Same as "simple linear regression", we will apply "least square criterion" for residual error, $$\begin{eqnarray}S &=& \sum ^n_{t=1}\epsilon_t^2 \\ &=&\sum ^n_{t=1}\left(\vec{y_t} -\vec{\hat{y_t}}\right)^2 \\ &=&\sum ^n_{t=1}\left(\vec{y} -X\beta \right)^2\\ &=& \left(\vec{y} -X\beta \right)^T\left(\vec{y} -X\beta \right) \end{eqnarray}$$

For the sake of best fitting, what we have to do is compute $$\frac{\partial S}{\partial \beta} = \vec{0}$$ $$\begin{eqnarray}S &=& \left(\vec{y} -X\beta \right)^T\left(\vec{y} -X\beta \right)\\ &=&\left(X\beta -\vec{y} \right)^T\left(X\beta - \vec{y} \right)\\\end{eqnarray}$$

$$\begin{eqnarray}\frac{\partial S}{\partial \beta}&=&2X^TX\beta -2X^Ty\\ &=& -2X^T\left(\vec{y}-X\vec{\beta}\right) = \vec{0}\end{eqnarray}$$$$\therefore\ \vec{\beta} = \left(X^TX\right)^{-1}X^T\vec{y}$$

However sometimes we can't find $\beta$ due to $\left(X^TX\right)^{-1}$ doesn't exist which is equivalent to $\left(X^TX\right)^{-1}$ is not regular matrix. In that situation, we can apply "regularization". I will write an article about "regularization" very soon.

2. Coefficient of determination¶

After creating "linear regression model", you might wanna assess how fit your model is. "Coefficient of determination" is the quotient of the variances of the fitted value and observed values of dependent variable. Let $S_y$ and $S_{\epsilon}$ be,
$$S_y = \frac{1}{n}\sum^{n}_{t=1}\left(y_t - \bar{y}\right)^2$$ where $\bar{y}$ is mean of y. $$S_\epsilon = \frac{1}{n}\sum^{n}_{t=1}\left(y_t - \hat{y_t}\right)^2$$ and $S_r$ is
$$S_r = S_y - S_\epsilon$$
Then "coefficient determination " $R^2$ can be denoted as folows,

$$R^2 = \frac{S_r}{S_y} = 1 - \frac{S_\epsilon}{S_y}$$

It's tribial that $R$ is $0\leqq R^2 \leqq 1$, and as bigger the coefficient determinant is, linear regression line fits well.

Poisson Mixture Model ③

This article is the continuation from Poisson Mixture Model②. In this article, I'm gonna implement poisson mixture model in practice. I applied gibbs sampling to approximate posterior distribution of poisson mixture model.

0. Implementation of Gibbs sampling¶

Fist of all, we prepare sample data or observed data used in Poisson Mixture Model②.

In [48]:

# Required Import modulejj
from scipy.stats import poisson
from scipy.stats import gamma
from scipy.stats import dirichlet
from scipy.stats import multinomial
import numpy as np
import matplotlib.pyplot as plt
% matplotlib inline
import collections
import pandas as pd 
from matplotlib import cm

In [2]:

# Sample from two different poisson distribution.
# whose parameter is 15, 30 respectively.
sample_15 = poisson.rvs(15,size=200)
sample_30 = poisson.rvs(30,size=300)
sample_freq = np.concatenate([sample_15,sample_30])

We discussed how to avoid overflow when it comes to logsumexp in Poisson Mixture Model②. The following is the simple implementation.

In [3]:

def log_sum_exp(exponent):
    """
    Argument :
    assume exponent is matrix which is shapen 
    (n (number of x), k (number of cluster))
    
    Return :
    The value of "const" to be used normalization.
    It's gonna be n (number of x) dimentional vector.
    """
    max_val = np.max(exponent, axis=1).reshape(-1,1)
    return np.log(np.sum(np.exp(exponent - max_val),
                         axis=1).reshape((-1,1))) + max_val

We also implement function in order to sample from posterior distribution according to gibbs sampling. In this implementatio, tentatively hyper parameter were set 1 for a and b, $(\frac{1}{2},\frac{1}{2})$ for $\alpha$.

In [4]:

def gibbs_mix_poi(X, k_num, iter_num):
    """
    parameter:
    -------------------------
    Assume X be numpy
    """
    # Hyper parameter
    a = 1
    b = 1
    alpha = np.ones((k_num,1))
    
    # Initialize lambda
    lmd = np.ones((k_num,1))
    
    # Initialize pi 
    pi = np.ones((k_num,1)) / k_num
    
    # Something to store result of gibbs sampling
    sampled_s = np.empty((iter_num,X.shape[0],k_num))
    sampled_lmd = np.empty((iter_num,k_num))
    sampled_pi = np.empty((iter_num,k_num))
    
    for i in range(iter_num):
        # Compute eta
        exponent = np.dot(X, np.log(lmd).reshape(1,-1)) 
                                        - lmd.reshape(1,-1) + np.log(pi).reshape(1,-1)
        const = - log_sum_exp(exponent)
        eta = np.exp(exponent + const)
        #Sample lambda
        S = np.array([ multinomial.rvs(n=1,p=temp_pi) for temp_pi in eta])
        sampled_s[i]=S

        #Sample lambda
        hat_a = np.dot(X.T,S).reshape(-1,1) + a
        hat_b = np.sum(S,axis=0).reshape(-1,1) + b
        lmd =np.array([ gamma.rvs(a=tempa,scale=1/tempb) for
                                       tempa, tempb in zip(hat_a,hat_b) ])

        sampled_lmd[i] = lmd
        
        # Sample pi
        hat_alpha = np.sum(S,axis=0).reshape(-1,1) + alpha
        pi = dirichlet.rvs(alpha=hat_alpha.reshape(-1),size=1).reshape(-1,1)
        sampled_pi[i] = pi.reshape(-1)
    
    return sampled_s, sampled_lmd, sampled_pi

Now we are ready to sample with gibbs sampling. Here we're gonna sample with five hundred iterations.

In [6]:

# Sample with 500 iteration
sample_s,sample_lmd,sample_pi = gibbs_mix_poi(
                                                      sample_freq.reshape(-1,1),2,500)

Let's check the result of sample :) First of all, we'll see $\lambda$. As you may know, Expectation of poisson distribution is $\lambda$. As we prepared observed data from two different poisson distribution whose parameter is 15 and 30. Hence we expect something close to 15 and 30.

In [7]:

# Check sampled lambda 
df_lam = pd.DataFrame(sample_lmd, columns=['lambda1','lambda2'])

# Plot lambda
for lam in df_lam.columns:
    df_lam[lam].plot.hist(bins=20)
    plt.title('Sample of {}'.format(lam),fontsize=16)
    plt.text(19,10,'Mean of lambda = \n{}'.format(df_lam[lam].mean()))
    plt.xlim((10,50))
    plt.show()

As you can see expectation of $\lambda$ is close to what it should be :)
Next we're gonna check the sample of $\pi$.

In [8]:

df_pi = pd.DataFrame(sample_pi,columns=['pi_1','pi_2'])

# plot histogram of sample of pi
for pi in df_pi.columns:
    df_pi[pi].plot.hist(bins=20)
    plt.title('Sample of {}'.format(pi),fontsize=16)
    plt.text(0.22,10,'Mean of pi =\n {}'.format(df_pi[pi].mean()))
    plt.xlim((0.2,0.8))
    plt.show()

We observed data with ratio of 4 to 6. The outcome is 3.9 to 6.1. It catches the feature of samples well . :)

Checking the outcome of clustering is a little bit tricky. I implemented with folloing steps.

1. Classify sample data into multiple bin. In this example, I set 41.
2. Compute average of which cluster the sample belongs to.
3. Plot histgram and color according to the number obtained in 2.

In [36]:

# Create binn
max_samp = np.max(sample_freq)
min_samp = np.min(sample_freq)
binn = np.linspace(min_samp,max_samp,41)

# Compute assignment of sample data
assign = np.digitize(sample_freq,binn)

sample_s_t = sample_s.transpose(1,0,2)

ratio_list = np.empty((binn.shape[0],2))

# Compute ratio of cluster for each bin
for i, assign_num in enumerate(np.unique(assign)):
    
    sum_per_iter = np.sum(sample_s_t[assign == assign_num],axis=1)
    sum_iter = np.sum(sum_per_iter,axis=1)
    ratio = sum_per_iter/sum_iter[:,np.newaxis]
    ratio_list[i] = np.mean(ratio,axis=0)

Following is the outcome of clustering of poisson mixture model. You can tell the sample around the 20 has light color since they are ambiguous about which poisson distribution yeilded. where as bigger and smaller number has darker color. I believe it go with your expectation :)

In [62]:

# Plot result of clustering
plt.figure(figsize=(8,4))
_,_, patches = plt.hist(sample_freq,bins=binn)
cm_map = cm.get_cmap('coolwarm')

# Create colors for each bin
colors = np.array([cm_map(ratio) for ratio in ratio_list[:,0]])

# Set color for each bin
for patch, color in zip(patches,colors):
    patch.set_fc(color)
    
plt.title('Result of clustering',fontsize=16)
plt.show()

Poisson Mixture Model ②

This article is continuation from Poisson mixture model①. In this article, we will cope with applying gibbs sampling to posterior distribution of Poisson Mixture Model.

0. Approximation with gibbs sampling¶

We've already create model of Poisson mixture model and sampled from it.

$$x \sim p(x | s, \lambda, \pi)$$

Now what we wanna know is posterior distribution $p(s,\lambda,\pi|X)$, when data $X = (x_1,x_2, \cdots, x_N)$ is observed. Nonetheless posterior distribution is itractable. Thereby approximation is required. On this article we're gonna utilize gibbs samping. If you wanna wrapp your head around gibbs sampling, you can check this link. Gentle introduction to Gibbs Sampling.
In mixture model, it's well known that sampling by sperating latent variable and parameter will give you simple distribution enough to sample from it. Therefore according to the following steps, gibbs sampling will be executed.
$$s^i \sim p(s|X,\lambda_{i-1},\pi_{i-1})$$ $$\lambda^i, \pi^i \sim p(\lambda,\pi | X,s^{i}) \tag{1}$$

1. Compute posterior distribution to sample $s$¶

From above $(1)$, what we know so as to draw sample is posterior distribution of $p(s\ |\ X, \lambda, \pi)$ and $p(\pi,\lambda\ |\ s, X)$. Let's start off with $p(s\ |\ X, \lambda, \pi)$ !

$$\begin{eqnarray}p(s\ |\ X, \lambda, \pi) &=& \frac{p(s, X, \lambda, \pi)}{p(X, \lambda, \pi)}\\ \\ &=& \frac{p(X|s,\lambda)p(s|\pi)p(\lambda)p(\pi)}{p(X, \lambda, \pi)}\\ \\ &=& \frac{p(X|s,\lambda)p(s|\pi)p(\lambda)p(\pi)}{p(X, |\ \lambda, \pi)p(\lambda)p(\pi)}\\ \\ &\propto& {p(X|s,\lambda)p(s|\pi)}\\ \\ &=& \prod^{N}_{n=1}{p(x_n|s_n,\lambda)p(s_n|\pi)}\\\ \end{eqnarray}$$$$\begin{eqnarray}\log p(x_n|s_n,\lambda) &=& \sum^{K}_{k} s_{nk}(x_n\log \lambda_k - \log x_n ! - \lambda_k) \end{eqnarray}$$

Since $\log x_n ! $ is already abserved. And $s_{nk}$ is one hot vector. We can treat it as constant.Therefore
$$\begin{eqnarray}\log p(x_n|s_n,\lambda) &=& \sum^{K}_{k} s_{nk}(x_n\log \lambda_k - \lambda_k) +\ const\end{eqnarray}$$ From model we set, $p(s_n|\pi)$ is, $$\log p(s_n|\pi) = \sum_k^Ks_{nk}\log \pi_k$$

Therefore, $$\begin{eqnarray}\log \{p(x_n|s_n,\lambda)p(s_n|\pi)\} = \sum^{K}_{k} s_{nk}(x_n\log \lambda_k - \lambda_k + \log \pi_k) +\ const\end{eqnarray}$$

Since $ \sum^{K}_{k} s_{nk} = 1$, we can tell $p(x_n|s_n,\lambda)p(s_n|\pi)$ is categorical distribution as following, $$s_n \sim Cat(s_n|\eta_n)$$ $$ s. t. $$ $$\eta_{nk} \propto exp\{x_n\log \lambda_k - \lambda_k + \log \pi_k \}\ s.\ t. \sum_k^K \eta_{nk} =1$$

2. logsumexp¶

In this section, I will share slick technique to implement above categorical distribution.
First of all, Since $\eta_{nk} \propto exp\{x_n\log \lambda_k - \lambda_k + \log \pi_k \}$, It seems that we have to compute this equation as following,

$$\eta_{nk} = \frac{exp\{x_n\log \lambda_k - \lambda_k + \log \pi_k \}}{\sum_k^K exp\{x_n\log \lambda_k - \lambda_k + \log \pi_k \}}\tag{2}$$

However, since it's exponential, the value $exp$ can be enormous value. At the time threre is possibility overflow in computation. 　Thus we might wanna think altenative mehod to compute $\eta$.　Let's say there is normalization constant "$const$". Then we can denote $\eta$ as following,

$$\begin{eqnarray}\eta &=& exp\{x_n\log \lambda_k - \lambda_k + \log \pi_k +\ const\} \\ &=& exp\{const\}exp\{x_n\log \lambda_k - \lambda_k + \log \pi_k\}\end{eqnarray}\tag{3}$$

Then from equation $(2)$, $exp\{const\}$ can be captured as below, $$\begin{eqnarray}exp\{const\} &=& \sum_k^K exp\{x_n\log \lambda_k - \lambda_k + \log \pi_k \}\\ const &=& \log \sum_k^K exp\{x_n\log \lambda_k - \lambda_k + \log \pi_k \}\end{eqnarray}\tag{4}$$ So, now we can apply "logsumexp" technique. "logsumexp" technique is folowing.
$$\begin{eqnarray}\log\sum^K_k \exp(x_k) &=& \log\sum^K_k \exp(x_k - x_{max})\exp(x_{max}) \\ &=& \log\left\{\sum^K_k \exp(x_k - x_{max})\right\}+ x_{max}\end{eqnarray}$$

That's how we can avoid overflow by calculating constant denoted in $(4)$ instead of $(2)$.
On the side note, for better understanding, the following is exactly what we're doing. let's say there is vector $(e^1, e^2, e^3, e^4)$. So as to make summatin is 1,

$$(e^{1-log(e^1 + e^2 + e^3+ e^4)},e^{2-log(e^1 + e^2 + e^3+ e^4)},e^{3-log(e^1 + e^2 + e^3+ e^4)},e^{4-log(e^1 + e^2 + e^3+ e^4}))$$

Instead of $$(\frac{e^{1}}{e^1 + e^2 + e^3+ e^4},\frac{e^{2}}{e^1 + e^2 + e^3+ e^4},\frac{e^3}{e^1 + e^2 + e^3+ e^4},\frac{e^4}{e^1 + e^2 + e^3+ e^4}))$$

3. Compute posterior distribution to sample $\lambda$ and $\pi$¶

Let's think about gibbs sampling of $\lambda$ and $\pi$,

$$\begin{eqnarray}p(\pi, \lambda| X,s) &\propto& p(\pi)p(\lambda)p(X|s,\lambda)p(s|\pi)\end{eqnarray}\tag{5}$$

From aboeve equation $(5)$ when it comes to $\lambda$ and $\pi$, we can sample seperately. We are off to good start. Let's look at $\pi$ first.

$$\begin{eqnarray}p(\pi)p(S|\pi) &=& p(\pi)\prod^N_np(s_n|\pi) \\ \log (p(\pi)p(S|\pi))&=& \sum^K_k\left\{(\alpha_{k-1} + \sum^N_ns_{nk})\log \pi_k\right\}\end{eqnarray}$$

Therefore $\pi$ is drawn from Dirichlet Distribution. $$\pi_i \sim Dir(\pi|\hat{\alpha})$$ where $$\hat{\alpha_k} = \alpha_k + \sum^N_ns_{nk}$$

Next, let's get into gibbs sampling of $\lambda$. $$\begin{eqnarray}p(\lambda)p(X|s,\lambda) &=& \prod ^K_k \lambda^{\alpha -1}_k e^{-b\lambda_k} + \prod^N_n\prod^K_k \left\{\frac{\lambda^{x_n}_{k}}{x_n!}e^{-\lambda_k}\right\}^{S_{nk}}\\ \log(p(\lambda)p(X|s,\lambda) ) &=& \sum^K_k\left\{((a-1)+\sum^N_ns_{nk}x_n)\log\lambda_k - (b + \sum^N_n s_{nk}\lambda_k)\right\}\end{eqnarray}$$ Therefore $\lambda_i$ is drawn from gammbda distribution. $$\lambda_i \sim p(\lambda|\hat{a},\hat{b})$$ where $$\hat{a} = a+\sum^N_n S_{nk}x_n,\ \ \hat{b} = b + \sum^N_nS_{nk}$$

So far we've prepared all required posterior distribution. We're ready to implement Poisson Mixture Model. However it's already quite long enoiugh to be an article. Hence I'm gonna write implementation on another blog :)