# Problem 1

Suppose $(x_{1}, y_{1}), (x_{2}, y_{2}), \cdots, (x_{n}, y_{n}) \in \mathbb{R} \times \{0,1\}$. Set
\begin{equation*}
\begin{split}
\sigma(t) &= \frac{1}{1 + e^{-t}}
\\ p_{i}(a, b) &= \sigma(a + b x_{i})^{y_{i}} (1 - \sigma(a + b x_{i}))^{1 -y_{i}} 
\\ l_{i}(a, b) &= -\log p_{i}(a, b)
\\ &= -y_{i} \log (\sigma(a + b x_{i})) - (1 - y_{i}) \log((1 - \sigma(a + b x_{i})))
\\ l(a, b) &= \sum_{i=1}^{n} l_{i}(a , b)
\end{split}
\end{equation*}


## (a) 

\textbf{Proof:} 
\begin{equation*}
\begin{split}
\sigma(x) & = \frac{1}{1 + e^{-x}} = (1 + e^{-x})^{-1}
\\ \sigma '(x) &= (-1) (1 + e^{-x})^{-2} (0 + e^{-x}(-1))
\\ &=\frac{e^{-x}}{(1 + e^{-x})^{2}}
\\ \sigma(x) (1 - \sigma(x)) &= \frac{1}{1+e^{-x}} \left( 1 - \frac{1}{1+e^{-x}}\right) 
\\& = \frac{1}{1+e^{-x}} \left( \frac{1+e^{-x} - 1}{1+e^{-x}}\right) 
\\& = \frac{e^{-x}}{(1 + e^{-x})^{2}}
\end{split}
\end{equation*}

Therefore, we proved that $\sigma '(x)  = \sigma(x) (1 - \sigma(x))$.


## (b)

To get the $\bigtriangledown l_{i}(a, b)$, we derive $\frac{\partial l_{i}(a,b)}{\partial a}$ and $\frac{\partial l_{i}(a,b)}{\partial b}$ respectively.

\begin{equation*}
\begin{split} 
\frac{\partial l_{i}(a,b)}{\partial a} &= -y_{i} \frac{\sigma'(a + bx_{i})}{\sigma(a + bx_{i})}\frac{\partial (a + bx_{i})}{\partial a} - (1 - y_{i})\frac{-\sigma'(a + bx_{i})}{1 - \sigma(a + bx_{i})}\frac{\partial (a + bx_{i})}{\partial a}
\\& = -y_{i}\frac{\sigma(a + bx_{i})(1 - \sigma(a + bx_{i}))}{\sigma(a + bx_{i})} + (1 - y_{i})\frac{\sigma(a + bx_{i})(1 - \sigma(a + bx_{i}))}{1 - \sigma(a + bx_{i})}
\\& = -y_{i}(1 - \sigma(a + bx_{i})) + (1 - y_{i}) \sigma(a + bx_{i})
\\& = \sigma(a + bx_{i}) - y_{i}
\end{split}
\end{equation*}

\begin{equation*}
\begin{split} 
\frac{\partial l_{i}(a,b)}{\partial b} &= -y_{i} \frac{\sigma'(a + bx_{i})}{\sigma(a + bx_{i})}\frac{\partial (a + bx_{i})}{\partial b} - (1 - y_{i})\frac{-\sigma'(a + bx_{i})}{1 - \sigma(a + bx_{i})}\frac{\partial (a + bx_{i})}{\partial b}
\\& = -y_{i}\frac{\sigma(a + bx_{i})(1 - \sigma(a + bx_{i}))}{\sigma(a + bx_{i})}x_{i} + (1 - y_{i})\frac{\sigma(a + bx_{i})(1 - \sigma(a + bx_{i}))}{1 - \sigma(a + bx_{i})}x_{i}
\\& = -y_{i}(1 - \sigma(a + bx_{i}))x_{i} + (1 - y_{i}) \sigma(a + bx_{i})x_{i}
\\& = [\sigma(a + bx_{i}) -  y_{i}]x_{i}
\end{split}
\end{equation*}

Then, \begin{equation*}
\begin{split}
\bigtriangledown l_{i}(a, b) &= \begin{pmatrix}
\frac{\partial l_{i}(a,b)}{\partial a}\\ 
\frac{\partial l_{i}(a,b)}{\partial b}
\end{pmatrix}
\\& = \begin{pmatrix}
 \sigma(a + bx_{i}) - y_{i}\\ 
[\sigma(a + bx_{i}) - y_{i}]x_{i}
\end{pmatrix}
\\& = (\sigma(a + bx_{i}) - y_{i}) \begin{pmatrix}
1\\ 
x_{i}
\end{pmatrix} 
\end{split}
\end{equation*}

And 
\begin{equation*}
\begin{split} 
\bigtriangledown ^{2} l_{i}(a, b) &= \begin{pmatrix}
\frac{\partial ^{2} l_{i}(a,b)}{\partial a^{2}} & \frac{\partial ^{2} l_{i}(a,b)}{\partial a \partial b}\\  
\frac{\partial ^{2} l_{i}(a,b)}{\partial b \partial a} & \frac{\partial ^{2} l_{i}(a,b)}{\partial b^{2}}
\end{pmatrix}
\\& = \begin{pmatrix}
 \sigma'(a + bx_{i})  &   \sigma'(a + bx_{i}) x_{i}\\
 \sigma'(a + bx_{i}) x_{i}  &   \sigma'(a + bx_{i}) x_{i}^{2}
\end{pmatrix}
\\& =  \sigma'(a + bx_{i}) \begin{pmatrix}
1 & x_{i}\\ 
x_{i} & x_{i}^{2}
\end{pmatrix}
\end{split}
\end{equation*}

Now we are going to prove $\bigtriangledown ^{2} g$ is positive-semidefinite but not positive-definite, making
$g(a, b)$ convex but not strictly convex. First, show the following proposition.

\textbf{[Proposition 1]}

Assume that $g : \mathbb{R}^{2} \mapsto \mathbb{R}$ is differentiable, for every $\mathbf{x}$ and $\mathbf{y}$, if the inequality
$$g(\mathbf{y}) \geqslant g(\mathbf{x}) +\bigtriangledown ^{2} g(\mathbf{x}) (\mathbf{y} - \mathbf{x})$$
is satisfied, then $g$ is convex.

\textbf{Proof}:

Let $\mathbf{w}$, $\mathbf{z} \in \mathbb{R}^{2}$, and $0 \leqslant \lambda \leqslant 1$. Denote $\mathbf{x} := \lambda \mathbf{w} + (1 - \lambda)\mathbf{z}$. Then the inequality implies 
\begin{equation} \label{inequality}
g(\mathbf{w}) \geqslant g(\mathbf{x}) +\bigtriangledown ^{2} g(\mathbf{x}) (\mathbf{w} - \mathbf{x})\\
g(\mathbf{z}) \geqslant g(\mathbf{x}) +\bigtriangledown ^{2} g(\mathbf{x}) (\mathbf{z} - \mathbf{x})
\end{equation}

Note that 
\begin{equation*}
\mathbf{w} - \mathbf{x} = (1 - \lambda) (\mathbf{w} - \mathbf{z})\\
\mathbf{z} - \mathbf{x} = \lambda (\mathbf{z} - \mathbf{w})
\end{equation*}

We multiply the first line in \ref{inequality} with $\lambda$ and the second line with $1 - \lambda$, and then add the two inequality in \ref{inequality} together, we obtain
\begin{equation*}
\begin{split}
\lambda g(\mathbf{w}) + (1 - \lambda) g(\mathbf{z}) &\geqslant g (\mathbf{x}) + \lambda \bigtriangledown^{\top} g(\mathbf{x})(1 - \lambda)(\mathbf{w} - \mathbf{z}) + \lambda \bigtriangledown^{\top} g(\mathbf{x})(1 - \lambda)(\mathbf{z} - \mathbf{w})
\\ &= g(\lambda \mathbf{w} + (1 - \lambda)\mathbf{z})
\end{split}
\end{equation*}
Since $\mathbf{w}$ and $\mathbf{z}$ are arbitrary, this proves $g$ is a convex function.

$\square$

\textbf{[Proposition 2]}

Now deduce that $\bigtriangledown ^{2} g$ is positive-semidefinite but not positive-definite, making
$g(a, b)$ convex but not strictly convex.

\textbf{Proof}:

Now assume the Hessian $\bigtriangledown ^{2} g$ is positive-semidefinite for all $\mathbf{x} \in \mathbb{R}^{2}$. Then $\forall \mathbf{x}$ and $\mathbf{y} \in \mathbb{R}^{2}$, the Taylor's theorem implies that
\begin{equation*}
g(\mathbf{y}) = g(\mathbf{x}) + \bigtriangledown ^{\top} g(\mathbf{x})(\mathbf{y} - \mathbf{x}) +\frac{1}{2}(\mathbf{y} - \mathbf{x})^{\top}\bigtriangledown ^{2} g(t\mathbf{x} + (1-t)\mathbf{y}) (\mathbf{y} - \mathbf{x})
\end{equation*}
for some $0 \leqslant t \leqslant 1$. Since $\bigtriangledown ^{2} g$ is everywhere positive semi-definite, the quadratic
term in this equation is always non-negative. Thus we can estimate
\begin{equation*}
g(\mathbf{y}) \geqslant g(\mathbf{x}) + \bigtriangledown ^{\top} g(\mathbf{x})(\mathbf{y} - \mathbf{x})
\end{equation*}
By \textbf{proposition 1}, $g$ is convex but not strictly convex.

$\square$

Or we can use \textbf{Theorem} directly , see Sydsaeter and Hammond (1995), Theorems 17.13 (p. 640).

\begin{theorem}
Let $z = f(\mathbf{x})$ be a $C^{2}$ function in an open convex set $S$ in $R^{n}$. If $H(\mathbf{x})$ denotes the Hessian matrix of f,then:

$f$ is convex $\Leftrightarrow$ $H(\mathbf{x})$ is positive semidefinite for all $\mathbf{x} \in S$.
\end{theorem}


## (c)

From (b), we have 
\begin{equation*}
\bigtriangledown ^{2} l_{i}(a, b) = 
\begin{pmatrix}
 \sigma'(a + bx_{i})  &   \sigma'(a + bx_{i}) x_{i}\\
 \sigma'(a + bx_{i}) x_{i}  &   \sigma'(a + bx_{i}) x_{i}^{2}
\end{pmatrix}
\end{equation*}

Based on the definition of null space, let $\bigtriangledown ^{2} l_{i}(a, b) (m_{i}, n_{i})^{\top} = (0, 0)^{\top}$, we have 

\begin{equation*}
\left\{\begin{matrix}
 \sigma'(a + bx_{i}) m_{i} + \sigma'(a + bx_{i}) x_{i} n_{i} = 0 & (1)\\ 
 \sigma'(a + bx_{i}) x_{i} m_{i} + \sigma'(a + bx_{i}) x_{i}^{2} n_{i} = 0 & (2) 
\end{matrix}\right.
\end{equation*}

 From (1), we have $m_{i} = -x_{i}n_{i}$. Let $n_{i} = 1$, then $m_{i} = -x_{i}$. Normalize $n$ and $m$, we have 
 \begin{equation*}
 \begin{split}
n_{i0} &= \frac{1}{\sqrt{x_{i}^{2}+1}}
\\ m_{i0}& = -\frac{x_{i}}{\sqrt{x_{i}^{2}+1}}
\end{split}
\end{equation*}

Therefore, a basis of the null space $N(\bigtriangledown ^{2} l_{i}(a, b))$ is 
\begin{equation*}
\begin{pmatrix}
-\frac{x_{i}}{\sqrt{x_{i}^{2}+1}}\\
\frac{1}{\sqrt{x_{i}^{2}+1}}
\end{pmatrix}
\end{equation*}

Moreover, we can write the null space as 
\begin{equation*}
N(\bigtriangledown ^{2} l_{i}(a, b)) = \left\{ (m_{i}, n_{i})^{\top} : m_{i} =  -\frac{x_{i}}{\sqrt{x_{i}^{2}+1}} t_{i}, n_{i} = \frac{1}{\sqrt{x_{i}^{2}+1}}t_{i}, \text{$t_{i}$ is any value and $t_{i} \in \mathbb{R} $ } \right\}
\end{equation*}


## (d)

Note that $N(\bigtriangledown ^{2} l_{i}(a, b)) = \left\{ (m_{i}, n_{i})^{\top} : m_{i} =  -\frac{x_{i}}{\sqrt{x_{i}^{2}+1}} t_{i}, n_{i} = \frac{1}{\sqrt{x_{i}^{2}+1}}t_{i}, t_{i} \in \mathbb{R} \right\}$, then,

\begin{equation*}
\begin{split}
\bigcap_{i=1}^{n} N(\bigtriangledown ^{2} l_{i}(a, b)) &= \left \{ \left \{ \begin{pmatrix}
-\frac{x_{1}}{\sqrt{x_{1}^{2}+1}}t_{1}\\
\frac{1}{\sqrt{x_{1}^{2}+1}}t_{1}
\end{pmatrix}   \right\} \bigcap \left \{ \begin{pmatrix}
-\frac{x_{2}}{\sqrt{x_{2}^{2}+1}}t_{2}\\
\frac{1}{\sqrt{x_{2}^{2}+1}}t_{2}
\end{pmatrix}   \right\} \bigcap  \cdots \bigcap
\left \{ \begin{pmatrix}
-\frac{x_{n}}{\sqrt{x_{n}^{2}+1}}t_{n}\\
\frac{1}{\sqrt{x_{n}^{2}+1}}t_{n}
\end{pmatrix}   \right\}\right\}
\\& = \left \{  \begin{pmatrix}
0\\ 
0
\end{pmatrix}\right \}
\end{split}
\end{equation*}

Since when $x_{i} \neq x_{j}$, then
\begin{equation*}
 \begin{pmatrix}
-\frac{x_{i}}{\sqrt{x_{i}^{2}+1}}\\
\frac{1}{\sqrt{x_{i}^{2}+1}}
\end{pmatrix} \neq  \begin{pmatrix}
-\frac{x_{j}}{\sqrt{x_{j}^{2}+1}}\\
\frac{1}{\sqrt{x_{j}^{2}+1}}
\end{pmatrix}
\end{equation*}


Now we are going to prove that \begin{equation*}
\begin{pmatrix}
-\frac{x_{i}}{\sqrt{x_{i}^{2}+1}}t_{i}\\
\frac{1}{\sqrt{x_{i}^{2}+1}}t_{i}
\end{pmatrix} =  \begin{pmatrix}
-\frac{x_{j}}{\sqrt{x_{j}^{2}+1}}t_{j}\\
\frac{1}{\sqrt{x_{j}^{2}+1}}t_{j}
\end{pmatrix}
\end{equation*}
only when they are $0$ vector.

If $-\frac{x_{i}}{\sqrt{x_{i}^{2}+1}}t_{i} = -\frac{x_{j}}{\sqrt{x_{j}^{2}+1}}t_{j} \neq 0$, then we have $t_{i} = \frac{x_{j}}{\sqrt{x_{j}^{2}+1}} \frac{\sqrt{x_{i}^{2}+1}}{x_{i}}t_{j}$. Plug this $t_{i}$ into $\frac{1}{\sqrt{x_{i}^{2}+1}}t_{i}$, we have $\frac{1}{\sqrt{x_{i}^{2}+1}}t_{i} = \frac{1}{\sqrt{x_{i}^{2}+1}}\frac{x_{j}}{\sqrt{x_{j}^{2}+1}} \frac{\sqrt{x_{i}^{2}+1}}{x_{i}}t_{j} = \frac{x_{j}}{x_{i}\sqrt{x_{j}^{2}+1}}t_{j} \neq \frac{1}{\sqrt{x_{j}^{2}+1}}t_{j}$. Therefore, we can not find $t_{i}$ and $t_{j}$ such that \begin{equation*}
\begin{pmatrix}
-\frac{x_{i}}{\sqrt{x_{i}^{2}+1}}t_{i}\\
\frac{1}{\sqrt{x_{i}^{2}+1}}t_{i}
\end{pmatrix} =  \begin{pmatrix}
-\frac{x_{j}}{\sqrt{x_{j}^{2}+1}}t_{j}\\
\frac{1}{\sqrt{x_{j}^{2}+1}}t_{j} 
\end{pmatrix} \neq 0
\end{equation*}

## (e)

From (b) we have 
\begin{equation*}
\bigtriangledown ^{2} l_{i}(a, b) = 
\begin{pmatrix}
 \sigma'(a + bx_{i})  &   \sigma'(a + bx_{i}) x_{i}\\
 \sigma'(a + bx_{i}) x_{i}  &   \sigma'(a + bx_{i}) x_{i}^{2}
\end{pmatrix}
\end{equation*}

Then $\forall \mathbf{s} = (s_{1}, s_{2})^{\top} \neq 0$, we have 
\begin{equation*}
\begin{split}
\mathbf{s}^{\top} \bigtriangledown ^{2} l_{i}(a, b) \mathbf{s} &= (s_{1}, s_{2}) \begin{pmatrix}
 \sigma'(a + bx_{i})  &   \sigma'(a + bx_{i}) x_{i}\\
 \sigma'(a + bx_{i}) x_{i}  &   \sigma'(a + bx_{i}) x_{i}^{2}
\end{pmatrix} \begin{pmatrix}
s1\\ 
s2
\end{pmatrix}
\\& = \left(s_{1}\sigma'(a + bx_{i})+ s_{2} \sigma'(a + bx_{i}) x_{i}
, s_{1}\sigma'(a + bx_{i}) x_{i} + s_{2}\sigma'(a + bx_{i}) x_{i}^{2} \right) \begin{pmatrix}
s1\\ 
s2
\end{pmatrix}
\\& = s_{1}^{2}\sigma'(a + bx_{i})+ s_{2}s_{1} \sigma'(a + bx_{i}) x_{i} + s_{1}s_{2}\sigma'(a + bx_{i}) x_{i} + s_{2}^{2}\sigma'(a + bx_{i}) x_{i}^{2}
\\& = s_{1}^{2}\sigma'(a + bx_{i}) + 2 s_{1}s_{2}\sigma'(a + bx_{i}) x_{i} + s_{2}^{2}\sigma'(a + bx_{i}) x_{i}^{2}
\\& = \sigma'(a + bx_{i}) (s_{1} + s_{2}x_{i})^{2}
\\& > 0
\end{split}
\end{equation*}
Thus we have $\bigtriangledown ^{2} l_{i}(a, b)$ is positive definite. Since $\bigtriangledown ^{2} l(a, b) = \sum_{i=1}^{n}\bigtriangledown ^{2} l_{i}(a, b)$, therefore, $\bigtriangledown ^{2} l(a, b)$ is also positive definite, since 
\begin{equation*}
\mathbf{s}^{\top} \bigtriangledown ^{2} l(a, b) \mathbf{s} = \sum_{i=1}^{n} \mathbf{s}^{\top} \bigtriangledown ^{2} l_{i}(a, b) \mathbf{s}>0
\end{equation*}

Then by the result in (b), $l(a,b)$ is a convex function.

## (f)

From (e), for $x_{i} \neq x_{j}$, we have $l(a, b)$ is a convex function. Suppose there are two optimal solutions $\mathbf{\omega} = (a_{1}, b_{1})^{\top}$, $\mathbf{\nu} = (a_{2}, b_{2})^{\top} \in \mathbb{R^{2}}$ s.t. 

\begin{equation} \label{l_function}
l(\mathbf{\omega}) \leqslant l(\mathbf{z}), \  l(\mathbf{\nu}) \leqslant l(\mathbf{z}), \
\forall \mathbf{z} \in \mathbb{R}^{2}
\end{equation}
Let $\mathbf{z} = \frac{\omega + \nu}{2}$, by convexity of function $l(a, b)$, we have 
\begin{equation} \label{l_inequality}
\begin{split}
l(\mathbf{z}) = l(\frac{\omega + \nu}{2}) &= l(\frac{1}{2}\omega + \frac{1}{2} \nu)
\\& \leqslant \frac{1}{2}l(\omega) + \frac{1}{2}l( \nu)
\\& \leqslant \frac{1}{2}l(\omega) +  \frac{1}{2}l(\omega)
\\& = l(\omega)
\end{split}
\end{equation}

But this contradicts with \ref{l_function}. Therefore, the minimizer of $l(a, b)$ is unique, which means the maximize likelihood estimates $\hat{a}$, $\hat{b}$ of function $-l(a, b)$ is unique. 
Note that the second inequality above is because the $\nu$ is arbitrary so we can take $\nu = \omega$.

# Problem 2

## （a）

\textbf{Proof}:

We have 
\begin{equation*}
\begin{split}
\\ \sigma '(x) &= (-1) (1 + e^{-x})^{-2} (0 + e^{-x}(-1))
\\ &=\frac{e^{-x}}{(1 + e^{-x})^{2}}
\\ & = \sigma(x) (1 - \sigma(x))
\\& > 0
\end{split}
\end{equation*}
Then $\sigma(x)$ is increasing. Take $x$ approximate to $-\infty$ and $\infty$, we can get the range of 
\begin{equation*}
\sigma(x) = \frac{1}{1 + e^{-x}} \in (0,1)
\end{equation*}

Secondly, $\log(x)$ is also increasing in $(0, \infty)$. Specifically, $\log(1) = 0$. Therefore, the range of $\log(x)$ in $(0,1)$ is $(-\infty, 0)$.

Since 
\begin{equation*}
\begin{split}
l_{i}(a, b) &= -\log p_{i}(a, b)
\\ &= -y_{i} \log (\sigma(a + b x_{i})) - (1 - y_{i}) \log(1 - \sigma(a + b x_{i}))
\end{split}
\end{equation*}

For $y_{i} = 0$, we have $l_{i}(a, b) = - \log(1 - \sigma(a + b x_{i})) \in (0, \infty)$, since $1 - \sigma(a + b x_{i}) \in (0,1)$.

For $y_{i} = 1$, we have $l_{i}(a, b) = - \log (\sigma(a + b x_{i})) \in (0, \infty)$, since $\sigma(a + b x_{i}) \in (0,1)$.

From the above, $l_{i}(a, b) > 0$ for all $(a, b) \in \mathbb{R}^{2}$.

## (b)

We have 
\begin{equation*}
\begin{split}
l_{i}(t v_{1}, t v_{2}) = -y_{i} \log (\sigma(t v_{1} + t v_{2} x_{i})) - (1 - y_{i}) \log(1 - \sigma(t v_{1} + t v_{2} x_{i}))
\end{split}
\end{equation*}

When $y_{i} = 0$, $l_{i}(t v_{1}, t v_{2}) =  - \log(1 - \sigma(t v_{1} + t v_{2} x_{i}))$, then
\begin{equation*}
\begin{split}
H_{i} &= \{(v_{1}, v_{2})^{\top} \in \mathbb{R}^{2} : \lim_{t \to \infty} l_{i}(t v_{1}, t v_{2}) = \infty \} 
\\& = \{(v_{1}, v_{2})^{\top} \in \mathbb{R}^{2} : \lim_{t \to \infty} [-\log(1 - \sigma(t v_{1} + t v_{2}x_{i}))] = \infty \} 
\\& = \{(v_{1}, v_{2})^{\top} \in \mathbb{R}^{2} : \lim_{t \to \infty} (1 - \sigma(t v_{1} + t v_{2}x_{i})) = 0 \}
\\& = \{(v_{1}, v_{2})^{\top} \in \mathbb{R}^{2} : \lim_{t \to \infty} (\sigma(t v_{1} + t v_{2}x_{i})) = 1 \} 
\\& = \{(v_{1}, v_{2})^{\top} \in \mathbb{R}^{2} : \lim_{t \to \infty} t (v_{1} +  v_{2}x_{i}) = \infty  \}
\\& = \{(v_{1}, v_{2})^{\top} \in \mathbb{R}^{2} :  v_{1} +  v_{2}x_{i} > 0 \}
\\& = \{(v_{1}, v_{2})^{\top} \in \mathbb{R}^{2} :  (v_{1}, v_{2}) \cdot (1, x_{i})^{\top} > 0  \}
\\& = H(w_{i})
\end{split}
\end{equation*}
where $w_{i} = (1, x_{i})^{\top}$.

When $y_{i} = 1$, $l_{i}(t v_{1}, t v_{2}) =  - \log(\sigma(t v_{1} + t v_{2} x_{i}))$, then
\begin{equation*}
\begin{split}
H_{i} &= \{(v_{1}, v_{2})^{\top} \in \mathbb{R}^{2} : \lim_{t \to \infty} l_{i}(t v_{1}, t v_{2}) = \infty \} 
\\& = \{(v_{1}, v_{2})^{\top} \in \mathbb{R}^{2} : \lim_{t \to \infty} [-\log(\sigma(t v_{1} + t v_{2}x_{i}))] = \infty \} 
\\& = \{(v_{1}, v_{2})^{\top} \in \mathbb{R}^{2} : \lim_{t \to \infty} \sigma(t v_{1} + t v_{2}x_{i}) = 0 \}
\\& = \{(v_{1}, v_{2})^{\top} \in \mathbb{R}^{2} : \lim_{t \to \infty} t (v_{1} +  v_{2}x_{i}) = - \infty  \}
\\& = \{(v_{1}, v_{2})^{\top} \in \mathbb{R}^{2} :  v_{1} +  v_{2}x_{i} < 0 \}
\\& = \{(v_{1}, v_{2})^{\top} \in \mathbb{R}^{2} :  - v_{1} -  v_{2}x_{i} > 0  \}
\\& = \{(v_{1}, v_{2})^{\top} \in \mathbb{R}^{2} :  (v_{1}, v_{2}) \cdot (-1, -x_{i})^{\top} > 0 \}
\\& = H(w_{i})
\end{split}
\end{equation*}
where $w_{i} = (-1, -x_{i})^{\top}$.

In conclusion, $w_{i} = (-1)^{y_{i}}\cdot(1, x_{i})^{\top}$.

## (c)


```python
from matplotlib import pyplot as plt
import pandas as pd
import numpy as np
import mpl_toolkits.axisartist as axisartist

#create blank figure
fig = plt.figure(figsize=(4, 4))
# use axisartist.Subplot to create a painting object ax
ax = axisartist.Subplot(fig, 111)  
#add this object ax into blank gigure
fig.add_axes(ax)

#use set_visibleto hidden x-axis, y-axis
ax.axis[:].set_visible(False)

#use ax.new_floating_axis to add new x-axis, y-axis
ax.axis["x"] = ax.new_floating_axis(0,0)
#draw arrow of x-axis
ax.axis["x"].set_axisline_style("->", size = 1.0)
#add y-axis and draw arrow of x-axis
ax.axis["y"] = ax.new_floating_axis(1,0)
ax.axis["y"].set_axisline_style("-|>", size = 1.0)
#set the direction of x-axis and y-axis
ax.axis["x"].set_axis_direction("top")
ax.axis["y"].set_axis_direction("right")

#generate the data with step 0.1
x1 = np.arange(0,15,0.1)
x1_v = np.arange(-15,15,0.1)
x2 = np.arange(-15,0,0.1)
x2_v = np.arange(-15,15,0.1)
x3 = np.arange(0,15,0.1)
x3_v = np.arange(-15,15,0.1)
#generate the data for v,w, and -u
y1 = 2*x1
y1_v = -0.5*x1_v
y2 = -0.3*x2
y2_v = (10/3)*x2_v
y3 = (-8/7)*x3
y3_v = (7/8)*x3_v
#set the range of x-axis and y-axis
plt.xlim(-12,12)
plt.ylim(-12, 12)
#draw the picture
plt.plot(x1,y1, c='b',label='v')
plt.plot(x2,y2, c='r', label='w')
plt.plot(x1_v,y1_v, c='b', ls='--')
plt.plot(x2_v,y2_v, c='r', ls='--')
plt.plot(x3,y3, c='y', label='u')
plt.plot(x3_v,y3_v, c='y', ls='--')
plt.legend(loc = 'upper left')
```

\textbf{Proof:}

"$\Leftarrow$ ": if $-\mathbf{u} = a \mathbf{v}+b\mathbf{w}$, for $\mathbf{v} = (v_{1}, v_{2})^{\top} \neq (0, 0)^{\top}$, $\forall \mathbf{m} \in H(\mathbf{v})$, $\mathbf{m}\cdot \mathbf{v} > 0$ implies the range of angle between $\mathbf{m}$ and $\mathbf{v}$, $\theta_{1} \in [0, \frac{\pi}{2})$. Similarly, for $\mathbf{w} = (w_{1}, w_{2})^{\top} \neq (0, 0)^{\top}$, $\forall \mathbf{m} \in H(\mathbf{w})$, $\mathbf{m}\cdot \mathbf{w} > 0$ implies the range of angle between $\mathbf{m}$ and $\mathbf{w}$,  $\theta_{2} \in [0, \frac{\pi}{2})$. For $a \mathbf{v}+b\mathbf{w} \neq (0, 0)^{\top}$, $\forall a, b > 0$, it must be located between $\mathbf{v}$ and $\mathbf{w}$.Therefore, the angle between $-\mathbf{u}$ and $\mathbf{v}$,  $-\mathbf{u}$ and $\mathbf{w}$ is big than $\frac{\pi}{2}$.Then, we can say, the range of $\mathbf{m} \in H(\mathbf{u})$ will have some overlap with the range of $\mathbf{m} \in H(\mathbf{v}) \cup H(\mathbf{w})$. In addition, the range of $\mathbf{m} \in H(\mathbf{v}) \cup H(\mathbf{w}) \cup H(\mathbf{u})$ is $\mathbb{R}^{2} - \{\mathbf{0}\}$.

Visually,for any $\mathbf{v}$ with color blue, the range of $H(\mathbf{v})$ is the area above that dotted blue line. For any $\mathbf{w}$ with color red, the range of $H(\mathbf{w})$ is the left area of that dotted red line. And for any $\mathbf{u}$ with color yellow, the range of $H(\mathbf{u})$ is the area below that dotted yellow line. Wherever $\mathbf{v}, \mathbf{w}$ and $\mathbf{u}$, the range of $H(\mathbf{v}) \cup H(\mathbf{w}) \cup H(\mathbf{u})$ is $\mathbb{R}^{2} - \{\mathbf{0}\}$.

"$\Rightarrow$ ": if $H(\mathbf{v}) \cup H(\mathbf{w}) \cup H(\mathbf{u}) = \mathbb{R}^{2} - \{\mathbf{0}\}$, the range of $H(\mathbf{v}) \cup H(\mathbf{w})$ must have some overlap with $H(\mathbf{u})$. It requires $-\mathbf{u} = a \mathbf{v}+b\mathbf{w}$.

## (d)

\textbf{Proof:}

Suppose $y_{i} = y_{k} = 0$, and $y_{j} = 1$, by (b), we have 
\begin{equation*}
H_{i} = H(\mathbf{w_{i}}), \ H_{j} = H(\mathbf{-w_{j}}), \ H_{k} = H(\mathbf{w_{k}})
\end{equation*}

where $\mathbf{w_{i}} = (1, x_{i})^{\top}$, $\mathbf{w_{j}} = (1, x_{j})^{\top}$, and $\mathbf{w_{k}} = (1, x_{k})^{\top}$. Note that $1 \leqslant i < j < k \leqslant n$, and $x_{i} < x_{j}<x_{k}$.

Solve $\mathbf{w_{j}} = a \mathbf{w_{i}} + b \mathbf{w_{k}}$, we can find 

\begin{equation*}
a = \frac{x_{j} - x_{k}}{x_{i} - x_{k}} > 0
\\
b = \frac{x_{i} - x_{j}}{x_{i} - x_{k}} >0 
\end{equation*}

Then by (c), $H_{i} \cup H_{j} \cup H_{k} = \mathbb{R}^{2} - \{ \mathbf{0} \}$ .

Similarly, we can prove when $y_{i} = y_{k} = 1$, and $y_{j} = 0$, then $H_{i} \cup H_{j} \cup H_{k} = \mathbb{R}^{2} - \{ \mathbf{0} \}$.

## (e)

\textbf{Proof:}

Let $S_{K} : = \{ (a, b)^{\top}: l(a,b) \leq K\}$. Suppose there $\exists \{ t(a, b)^{\top}: t \geq 0\} \subseteq S_{K}$, s.t., $ l(ta,tb) \leq K$. However, when $(i, j ,k)$ is an increasing sequence of indices that satisfies $\dagger$, from (b), let $i$ be such that

\begin{equation*}
l_{i}(ta, tb) = \infty, \quad t \rightarrow \infty 
\\ l(ta, tb) = \sum_{i} l_{i}(ta, tb) \geqslant  l_{i}(ta,tb)
\end{equation*}

It contradicts with the definition of $S_{K}$, therefore, $S_{K}$ contains no ray from the origin.

## (f)

To deduce $S_{K}$ is bounded, using the facts $(i)$ and $(ii)$, we only need to prove that $l(a, b)$ is a convex function. From problem 1 (e), we already proved that $l(a, b)$ is a convex function. Thus the result.

## (g)

\textbf{Proof:}

By the definition of infimum, $\forall \epsilon > 0$, there exists $(\hat{a},\hat{b})^{\top} \in S_{K}$ such that $m + \epsilon > l(\hat{a},\hat{b})$. And, by the fact that $ m =\underset{(a, b)^{\top} \in S_{K}}{\text{inf}}l(a,b)$, we have $\forall (a, b)^{\top} \in S_{K}$
\begin{equation*}
m \leqslant l(a,b)
\end{equation*}

By the arbitrariness of $\epsilon$, then there exists $(\hat{a},\hat{b})^{\top} \in S_{K}$ such that $l(\hat{a},\hat{b}) = m $.

Since in Problem 1 (e), we have known that $l(a,b)$ is a convex function, therefore, $m$ is the global minimum.

## (h)

\textbf{Proof:}

Suppose there are $(\hat{a}, \hat{b})^{\top}$ and $(a', b')^{\top}$ such that $l(\hat{a}, \hat{b}) = l(a',b') = m$. Note that 

\begin{equation*}
l(a,b) = \sum_{i}^{n}l_{i}(a,b)
\end{equation*}

Then, we have $l_{i}(\hat{a}, \hat{b}) = l_{i}(a',b')$, which implies $\sigma(\hat{a}+\hat{b}x_{i}) = \sigma(a'+b'x_{i})$, and then $\hat{a}+\hat{b}x_{i} = a'+b'x_{i}$. Therefore, \begin{equation*}
\hat{a} = a', \quad \hat{b} = b'
\end{equation*}

i.e., $(\hat{a}, \hat{b})^{\top}$  is the unique point at which $l$ takes on its minimum value.


# Problem 3

## (a)

\begin{equation*}
\begin{split}
E(S^{2}) &= E\left(\frac{(m-1)S_{X}^{2}+(n-1)S_{Y}^{2}}{m+n-2} \right)
\\& = \left(\frac{(m-1)E(S_{X}^{2})+(n-1)E(S_{Y}^{2})}{m+n-2} \right)
\\& = \left(\frac{(m-1)\sigma^{2}+(n-1)\sigma^{2}}{m+n-2} \right)
\\& = \sigma^{2}
\end{split}
\end{equation*}

Therefore, the pooled variance estimator $S^{2}$ is an unbiased estimator of $\sigma^{2}$.

## (b)

Simplify the equation,  we have
\begin{equation*}
\begin{split}
\frac{(m + n -2) S^{2}}{\sigma^{2}} & = \frac{(m + n -2)}{\sigma^{2}}\cdot \frac{(m-1)S_{X}^{2}+(n-1)S_{Y}^{2}}{m+n-2}
\\& = \frac{(m -1)S_{X}^{2}}{\sigma^{2}} +  \frac{(n -1)S_{Y}^{2}}{\sigma^{2}}
\end{split}
\end{equation*}

Since $ \frac{(m -1)S_{X}^{2}}{\sigma^{2}} \sim \chi^{2}(m-1)$, and $ \frac{(n -1)S_{Y}^{2}}{\sigma^{2}} \sim \chi^{2}(n-1)$, therefore, \begin{equation*}
\frac{(m + n -2) S^{2}}{\sigma^{2}} \sim \chi^{2}(m+n-2)
\end{equation*}

Then we have

\begin{equation*}
\begin{split}
 \text{Var}(\frac{(m + n -2) S^{2}}{\sigma^{2}}) = 2 (m+n-2)
\\ \Leftrightarrow \frac{(m+n-2)^{2}}{\sigma^{4}} \text{Var}(S^{2}) = 2 (m+n-2)
\\ \Leftrightarrow \text{Var}(S^{2}) = \frac{2 \sigma^{4}}{m+n-2}
\end{split}
\end{equation*}


# Problem 4

## (a) 


```python
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import confusion_matrix
from sklearn.metrics import classification_report
from matplotlib import pyplot as plt
import pandas as pd
import numpy as np
import random
```


```python
# input data
data = pd.read_csv("dataset_1.csv")
dataset = pd.DataFrame(data)

# rename the columns of dataset
dataset.columns = ['feature_1', 'feature_2','label']

#select subdataset by label 0 or 1
data1 = dataset[dataset['label'] == 1]
data0 = dataset[dataset['label'] == 0]
```


```python
#scatter plot
p1 = plt.scatter(data1['feature_1'], data1['feature_2'], marker = 'x', color = 'r', label = '1', s = 30)
p2 = plt.scatter(data0['feature_1'], data0['feature_2'], marker = '*', color = 'b', label = '0', s = 30)
plt.legend(loc = 'upper right')
```

## (b)

Fit the data to a logistic regression model and plot the decision boundary.

When the test data size is more than $10\%$ of the total data set, it shows that \textit{"no contour levels were found within the data range."} 


```python
X = dataset[['feature_1','feature_2']]
Y = dataset.label

X_train, X_test, y_train, y_test = train_test_split(X, Y, test_size = .1, random_state=np.random.randint(10e5))

# Use Logistic Regression model and fit the data
logreg = LogisticRegression(solver='lbfgs', multi_class='multinomial')
logreg.fit(X_train, y_train)

# Plot the decision boundary. For that, we will assign a color to each point in the mesh [x_min, x_max]x[y_min, y_max].
X_train = pd.DataFrame(X_train)
y_train = pd.DataFrame(y_train)
X_train['y_train'] = y_train
new_data = X_train
new_data0 = new_data[new_data['y_train'] == 0]
new_data1 = new_data[new_data['y_train'] == 1]

x_min, x_max = new_data[['feature_1']].min() - .5, new_data[['feature_1']].max() + .5
y_min, y_max = new_data[['feature_2']].min() - .5, new_data[['feature_2']].max() + .5
h = .02  # step size in the mesh
xx, yy = np.meshgrid(np.arange(x_min, x_max, h), np.arange(y_min, y_max, h))
Z = logreg.predict(np.c_[xx.ravel(), yy.ravel()])

# Put the result into a color plot
Z = Z.reshape(xx.shape)
plt.contour(xx, yy, Z)

# Plot also the training points
plt.scatter(new_data1['feature_1'], new_data1['feature_2'], marker = 'x', color = 'r', label = '1', s = 30)
plt.scatter(new_data0['feature_1'], new_data0['feature_2'], marker = '*', color = 'b', label = '0', s = 30)
#plt.scatter(X_train[['feature_1']], X_train[['feature_2']])
plt.xlabel('Feature 1')
plt.ylabel('Feature 2')

plt.xlim(xx.min(), xx.max())
plt.ylim(yy.min(), yy.max())
plt.xticks(())
plt.yticks(())

plt.show()
```


```python
X = dataset[['feature_1','feature_2']]
Y = dataset.label

X_train, X_test, y_train, y_test = train_test_split(X, Y, test_size = .1, random_state=np.random.randint(10e5))

# Use Logistic Regression model and fit the data
logreg = LogisticRegression(solver='lbfgs', multi_class='multinomial')
logreg.fit(X_train, y_train)

#Get the predictive value of Y
y_pred = logreg.predict(X_test)

print(classification_report(y_test, y_pred))
```

                 precision    recall  f1-score   support
    
            0.0       0.59      1.00      0.74        53
            1.0       0.00      0.00      0.00        37
    
    avg / total       0.35      0.59      0.44        90
    


    D:\Users\TUSBASA\AppData\Local\Continuum\anaconda3\lib\site-packages\sklearn\metrics\classification.py:1135: UndefinedMetricWarning: Precision and F-score are ill-defined and being set to 0.0 in labels with no predicted samples.
      'precision', 'predicted', average, warn_for)


Now change the training data size from $90\%$ to $10\%$.


```python
X = dataset[['feature_1','feature_2']]
Y = dataset.label

random.seed(317)

X_train, X_test, y_train, y_test = train_test_split(X, Y, test_size = .9, random_state=np.random.randint(10e5))

# Use Logistic Regression model and fit the data
logreg = LogisticRegression(solver='lbfgs', multi_class='multinomial')
logreg.fit(X_train, y_train)

#Get the predictive value of Y
y_pred = logreg.predict(X_test)

print(classification_report(y_test, y_pred))
```

                 precision    recall  f1-score   support
    
            0.0       0.67      1.00      0.80       541
            1.0       0.00      0.00      0.00       269
    
    avg / total       0.45      0.67      0.53       810
    


    D:\Users\TUSBASA\AppData\Local\Continuum\anaconda3\lib\site-packages\sklearn\metrics\classification.py:1135: UndefinedMetricWarning: Precision and F-score are ill-defined and being set to 0.0 in labels with no predicted samples.
      'precision', 'predicted', average, warn_for)



```python
# Plot the decision boundary. For that, we will assign a color to each point in the mesh [x_min, x_max]x[y_min, y_max].
X_train = pd.DataFrame(X_train)
y_train = pd.DataFrame(y_train)
X_train['y_train'] = y_train
new_data = X_train
new_data0 = new_data[new_data['y_train'] == 0]
new_data1 = new_data[new_data['y_train'] == 1]

x_min, x_max = new_data[['feature_1']].min() - .5, new_data[['feature_1']].max() + .5
y_min, y_max = new_data[['feature_2']].min() - .5, new_data[['feature_2']].max() + .5
h = .02  # step size in the mesh
xx, yy = np.meshgrid(np.arange(x_min, x_max, h), np.arange(y_min, y_max, h))
Z = logreg.predict(np.c_[xx.ravel(), yy.ravel()])

# Put the result into a color plot
Z = Z.reshape(xx.shape)
plt.contour(xx, yy, Z)

# Plot also the training points
plt.scatter(new_data1['feature_1'], new_data1['feature_2'], marker = 'x', color = 'r', label = '1', s = 30)
plt.scatter(new_data0['feature_1'], new_data0['feature_2'], marker = '*', color = 'b', label = '0', s = 30)
plt.xlabel('Feature 1')
plt.ylabel('Feature 2')

plt.xlim(xx.min(), xx.max())
plt.ylim(yy.min(), yy.max())
plt.xticks(())
plt.yticks(())

plt.show()
```

However, when the size of training data is small, the result here (the decision boundary) is unconvincing.

Fit the data to a Gaussian naive Bayes model and plot the decision boundary.


```python
from sklearn.naive_bayes import GaussianNB
from sklearn.model_selection import train_test_split
from matplotlib import pyplot as plt
import pandas as pd
import numpy as np
import math
```


```python
X = dataset[['feature_1','feature_2']]
Y = dataset.label

X_train, X_test, y_train, y_test = train_test_split(X, Y, test_size = .1, random_state=np.random.randint(10e5))

#Use a Gaussian naive Bayes model to fit the data
gnb = GaussianNB()
gnb.fit(X_train, y_train)

#Get the predictive value of Y
y_pred = gnb.predict(X_test)

print(classification_report(y_test, y_pred))
```

                 precision    recall  f1-score   support
    
            0.0       0.84      0.88      0.86        58
            1.0       0.76      0.69      0.72        32
    
    avg / total       0.81      0.81      0.81        90
    



```python
# Plot the decision boundary. For that, we will assign a color to each point in the mesh [x_min, x_max]x[y_min, y_max].
X_train = pd.DataFrame(X_train)
y_train = pd.DataFrame(y_train)
X_train['y_train'] = y_train
new_data = X_train
new_data0 = new_data[new_data['y_train'] == 0]
new_data1 = new_data[new_data['y_train'] == 1]

x_min, x_max = new_data[['feature_1']].min() - .5, new_data[['feature_1']].max() + .5
y_min, y_max = new_data[['feature_2']].min() - .5, new_data[['feature_2']].max() + .5
h = .02  # step size in the mesh
xx, yy = np.meshgrid(np.arange(x_min, x_max, h), np.arange(y_min, y_max, h))
Z = gnb.predict(np.c_[xx.ravel(), yy.ravel()])

# Put the result into a color plot
Z = Z.reshape(xx.shape)
plt.contour(xx, yy, Z,colors='y', linewidths=3)

# Plot also the training points
plt.scatter(new_data1['feature_1'], new_data1['feature_2'], marker = 'x', color = 'r', label = '1', s = 30)
plt.scatter(new_data0['feature_1'], new_data0['feature_2'], marker = '*', color = 'b', label = '0', s = 30)
plt.xlabel('Feature 1')
plt.ylabel('Feature 2')

plt.xlim(xx.min(), xx.max())
plt.ylim(yy.min(), yy.max())
plt.xticks(())
plt.yticks(())

plt.show()
```

Compare these two results, the Gaussian naive Bayes classifier yields a better result than the logistic regression classifier. 

## (c)


```python
# input data
data = pd.read_csv("dataset_1.csv")
dataset = pd.DataFrame(data)

# rename the columns of dataset
dataset.columns = ['feature_1', 'feature_2','label']

# augment the dataset with three new features
dataset['feature_3'] = dataset['feature_1'].map(lambda x: x**2)
dataset['feature_4'] =dataset['feature_1']*dataset['feature_2']
dataset['feature_5'] = dataset['feature_2'].map(lambda x: x**2)

dataset.head
```




    <bound method NDFrame.head of      feature_1  feature_2  label  feature_3  feature_4  feature_5
    0     3.675052  -6.344161    0.0  13.506006 -23.315121  40.248381
    1    -3.364928   7.041318    0.0  11.322742 -23.693530  49.580158
    2     1.694320   1.447161    0.0   2.870721   2.451954   2.094275
    3     4.054083  -1.304888    0.0  16.435592  -5.290127   1.702734
    4    -7.742027   1.264913    0.0  59.938987  -9.792989   1.600004
    5    -6.917131   8.893598    0.0  47.846706 -61.518185  79.096085
    6     3.383288  -3.178704    0.0  11.446638 -10.754471  10.104158
    7     3.962048  -4.795083    0.0  15.697828 -18.998350  22.992818
    8    -1.897284   5.798241    0.0   3.599685 -11.000909  33.619602
    9    -5.207891   1.907284    0.0  27.122130  -9.932926   3.637731
    10    3.083910  -6.765969    0.0   9.510502 -20.865640  45.778334
    11   -5.679789   1.613964    0.0  32.260000  -9.166972   2.604878
    12   -3.296934  -1.489919    0.0  10.869772   4.912165   2.219860
    13   -1.521670  -4.450176    0.0   2.315479   6.771699  19.804064
    14   -6.247812   5.942000    0.0  39.035153 -37.124495  35.307358
    15    0.768498   4.988647    0.0   0.590589   3.833766  24.886595
    16   -7.324032   3.193076    0.0  53.641451 -23.386191  10.195734
    17   -6.304129   2.691515    0.0  39.742040 -16.967656   7.244252
    18   -5.945891   3.728386    0.0  35.353620 -22.168576  13.900861
    19   -4.688559   4.869327    0.0  21.982587 -22.830128  23.710347
    20    0.418084   7.058720    0.0   0.174795   2.951141  49.825528
    21    3.184358   1.796015    0.0  10.140139   5.719154   3.225668
    22    9.393220  -4.313703    0.0  88.232583 -40.519557  18.608030
    23   -9.429204   5.420570    0.0  88.909891 -51.111662  29.382580
    24   -1.960239   3.262893    0.0   3.842536  -6.396050  10.646473
    25    3.110441  -3.565557    0.0   9.674844 -11.090454  12.713194
    26    4.626828  -5.718313    0.0  21.407538 -26.457650  32.699101
    27   -6.646018   4.664786    0.0  44.169554 -31.002254  21.760233
    28   -2.598596  -1.835376    0.0   6.752700   4.769401   3.368606
    29   -7.721713   5.997561    0.0  59.624858 -46.311444  35.970733
    ..         ...        ...    ...        ...        ...        ...
    869   2.574460  -1.128456    1.0   6.627842  -2.905164   1.273413
    870   1.328956   0.910903    1.0   1.766123   1.210550   0.829744
    871  -4.431544   4.639129    1.0  19.638578 -20.558503  21.521519
    872  -1.945830  -0.959378    1.0   3.786253   1.866786   0.920406
    873   1.689179   0.197319    1.0   2.853326   0.333307   0.038935
    874  -1.314542   0.165947    1.0   1.728022  -0.218144   0.027538
    875  -4.232121   0.357848    1.0  17.910847  -1.514457   0.128055
    876   0.766022  -2.724408    1.0   0.586789  -2.086955   7.422398
    877  -2.454223  -1.853266    1.0   6.023210   4.548327   3.434593
    878   1.416881  -1.860400    1.0   2.007552  -2.635965   3.461087
    879  -4.048996   0.631740    1.0  16.394372  -2.557914   0.399096
    880  -0.094495  -1.355411    1.0   0.008929   0.128079   1.837139
    881  -1.991145   2.031486    1.0   3.964660  -4.044984   4.126935
    882  -0.006729   1.976014    1.0   0.000045  -0.013296   3.904632
    883   4.412876  -2.281690    1.0  19.473472 -10.068813   5.206108
    884   2.307079   1.576972    1.0   5.322615   3.638199   2.486841
    885   0.402857   0.285474    1.0   0.162294   0.115005   0.081495
    886  -0.870524   2.314665    1.0   0.757812  -2.014971   5.357676
    887   1.345819  -1.793121    1.0   1.811230  -2.413217   3.215283
    888  -1.303370   0.334438    1.0   1.698773  -0.435896   0.111849
    889   1.329108  -2.457942    1.0   1.766528  -3.266871   6.041480
    890   0.842346  -0.682016    1.0   0.709546  -0.574494   0.465146
    891   0.802667  -0.093213    1.0   0.644274  -0.074819   0.008689
    892  -0.936393  -2.931473    1.0   0.876833   2.745011   8.593532
    893  -3.040933   3.440424    1.0   9.247271 -10.462099  11.836521
    894   0.869583  -1.158494    1.0   0.756175  -1.007407   1.342108
    895   3.642802  -4.796566    1.0  13.270008 -17.472940  23.007041
    896  -1.495733   3.365957    1.0   2.237218  -5.034574  11.329667
    897   1.569534  -1.398385    1.0   2.463436  -2.194813   1.955482
    898   0.891431  -1.951263    1.0   0.794649  -1.739416   3.807427
    
    [899 rows x 6 columns]>



## (d)


```python
X = dataset[['feature_1','feature_2', 'feature_3', 'feature_4', 'feature_5']]
Y = dataset.label

X_train, X_test, y_train, y_test = train_test_split(X, Y, test_size = .1, random_state=np.random.randint(10e5))

# Use Logistic Regression model and fit the data
logreg = LogisticRegression(solver='lbfgs', multi_class='multinomial')
logreg.fit(X_train, y_train)
w2 = logreg.coef_
w1 = logreg.intercept_


#Get the predictive value of Y
y_pred = logreg.predict(X_test)

print(w2, w1)
print(classification_report(y_test, y_pred))
```

    [[-0.01758604 -0.05236831 -0.11576063 -0.15439781 -0.1463601 ]] [1.33376643]
                 precision    recall  f1-score   support
    
            0.0       0.93      0.85      0.88        59
            1.0       0.75      0.87      0.81        31
    
    avg / total       0.87      0.86      0.86        90
    


The fitted logistic regression model is 

\begin{equation*}
log(\frac{P(\text{label} = 1 )}{1-P(\text{label} = 1 )}) = 1.334 -0.017x_{1} -0.0524 x_{2}  -0.116 x_{1}^{2} -0.154 x_{1}x_{2} -0.146 x_{2}^{2}
\end{equation*}



```python
X = dataset[['feature_1','feature_2', 'feature_3', 'feature_4', 'feature_5']]
Y = dataset.label

X_train, X_test, y_train, y_test = train_test_split(X, Y, test_size = .1, random_state=np.random.randint(10e5))

#Use a Gaussian naive Bayes model to fit the data
gnb = GaussianNB()
gnb.fit(X_train, y_train)
theta = gnb.theta_
sigma = gnb.sigma_
#e = gnb.epsilon_

#Get the predictive value of Y
y_pred = gnb.predict(X_test)

print(theta)
print(sigma)
print(classification_report(y_test, y_pred))
```

    [[ -0.08784594  -0.09504156  21.72559258 -13.09618901  22.16127563]
     [  0.03491755  -0.1615705    4.62415904  -2.55244443   3.81907758]]
    [[ 21.71787607  22.15224314 490.50762883 377.48702541 436.19409414]
     [  4.62294021   3.79297296  38.24083686  20.88025111  29.58109512]]
                 precision    recall  f1-score   support
    
            0.0       0.96      0.84      0.90        51
            1.0       0.82      0.95      0.88        39
    
    avg / total       0.90      0.89      0.89        90
    


Let $i = 1, 2, 3, 4, 5$, indicates the i-th element of the following vector,
\begin{equation*}
\mu_{i0} = (-0.08784594, -0.09504156, 21.72559258, -13.09618901, 22.16127563 )
\\ \mu_{i1} = (0.03491755, -0.1615705, 4.62415904, -2.55244443, 3.81907758)
\\ \sigma_{i0} = (21.71787607, 22.15224314, 490.50762883, 377.48702541, 436.19409414)
\\ \sigma_{i1} = (4.62294021, 3.79297296, 38.24083686, 20.88025111, 29.58109512)
\end{equation*}

Then, the fitted Gaussian naive Bayes model is 

\begin{equation*}
P(X_{i} = x | Y = k) = \frac{1}{\sqrt{2 \pi \sigma_{ik}}}exp\left\{-\frac{1}{2}(\frac{x-\mu_{ik}}{\sigma_{ik}})^{2}\right \}
\end{equation*}

Note that here $k = 0,1$.

## (e)

Plot the plane curve:  
$y = c_{0}+c_{1}x_{1}+c_{2}x_{2}+c_{3}x_{1}^{2}+c_{4}x_{1}x_{2}+c_{5}x_{2}^{2}$


```python
from matplotlib import pyplot as plt
plt.style.use('seaborn-white')

X = dataset[['feature_1','feature_2', 'feature_3', 'feature_4', 'feature_5']]
Y = dataset.label

X_train, X_test, y_train, y_test = train_test_split(X, Y, test_size = .1, random_state=np.random.randint(10e5))

logreg = LogisticRegression(solver='lbfgs', multi_class='multinomial')
logreg.fit(X_train, y_train)
w2 = logreg.coef_
w1 = logreg.intercept_
w = np.append(w1, w2)

x_min, x_max = X_train[['feature_1']].min() - .5, X_train[['feature_1']].max() + .5
y_min, y_max = X_train[['feature_2']].min() - .5, X_train[['feature_2']].max() + .5
h = .02  # step size in the mesh
xx, yy = np.meshgrid(np.arange(x_min, x_max, h), np.arange(y_min, y_max, h))


def f(x1, x2, w):
    return w[0] + w[1]*x1 + w[2]*x2 + w[3]*x1**2 + w[4]*x1*x2 + w[5]*x2**2

zz = f(xx, yy, w)
plt.contour(xx, yy, zz,levels=[0],colors='black')

```


```python
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import confusion_matrix
from sklearn.metrics import classification_report
from matplotlib import pyplot as plt
import pandas as pd
import numpy as np

# input data
data = pd.read_csv("dataset_1.csv")
dataset = pd.DataFrame(data)

# rename the columns of dataset
dataset.columns = ['feature_1', 'feature_2','label']

# augment the dataset with three new features
dataset['feature_3'] = dataset['feature_1'].map(lambda x: x**2)
dataset['feature_4'] =dataset['feature_1']*dataset['feature_2']
dataset['feature_5'] = dataset['feature_2'].map(lambda x: x**2)

#select subdataset by label 0 or 1
data1 = dataset[dataset['label'] == 1]
data0 = dataset[dataset['label'] == 0]

#scatter plot
p1 = plt.scatter(data1['feature_1'], data1['feature_2'], marker = 'x', color = 'r', label = '1', s = 30)
p2 = plt.scatter(data0['feature_1'], data0['feature_2'], marker = '*', color = 'b', label = '0', s = 30)
plt.legend(loc = 'upper right')

#plot decision boundary
X = dataset[['feature_1','feature_2', 'feature_3', 'feature_4', 'feature_5']]
Y = dataset.label

X_train, X_test, y_train, y_test = train_test_split(X, Y, test_size = .1, random_state=np.random.randint(10e5))

logreg = LogisticRegression(solver='lbfgs', multi_class='multinomial')
logreg.fit(X_train, y_train)
w2 = logreg.coef_
w1 = logreg.intercept_
w = np.append(w1, w2)

x_min, x_max = X_train[['feature_1']].min() - .5, X_train[['feature_1']].max() + .5
y_min, y_max = X_train[['feature_2']].min() - .5, X_train[['feature_2']].max() + .5
h = .02  # step size in the mesh
xx, yy = np.meshgrid(np.arange(x_min, x_max, h), np.arange(y_min, y_max, h))


def f(x1, x2, w):
    return w[0] + w[1]*x1 + w[2]*x2 + w[3]*x1**2 + w[4]*x1*x2 + w[5]*x2**2

zz = f(xx, yy, w)
plt.contour(xx, yy, zz,levels=[0],colors='yellow',linewidths = 3)

plt.xlabel('Feature 1')
plt.ylabel('Feature 2')
plt.xticks(())
plt.yticks(())
plt.show()
```

\textbf{Comment:} The decision boundary of logistic regression is linear, therefore it can not classify the samples with the shape of "concentric elliptical". The decision boundary of a Gaussian naive Bayes model yields a better result than logistic regression model, since it is an ellipse. However, it still fails to separate the two class as much as possible. As for the plane curve with degree of two, it yields the best result. 

To get the best classification result, we need to use different model and then compare the efficiency to choose the "best" model.

## (f)


```python
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import confusion_matrix
from sklearn.metrics import classification_report
from matplotlib import pyplot as plt
import pandas as pd
import numpy as np

# input data
data = pd.read_csv("dataset_2.csv")
dataset = pd.DataFrame(data)

# rename the columns of dataset
dataset.columns = ['feature_1', 'feature_2','label']

# augment the dataset with new features
dataset['feature_3'] = dataset['feature_1'].map(lambda x: x**2)
dataset['feature_4'] =dataset['feature_1']*dataset['feature_2']
dataset['feature_5'] = dataset['feature_2'].map(lambda x: x**2)
dataset['feature_6'] = dataset['feature_1'].map(lambda x: x**3)
dataset['feature_7'] = dataset['feature_3']*dataset['feature_2']
dataset['feature_8'] = dataset['feature_5']*dataset['feature_1']
dataset['feature_9'] = dataset['feature_2'].map(lambda x: x**3)
dataset['feature_10'] = dataset['feature_3']*dataset['feature_2']
dataset['feature_11'] = dataset['feature_5']*dataset['feature_1']

#select subdataset by label 0 or 1
data1 = dataset[dataset['label'] == 1]
data0 = dataset[dataset['label'] == 0]

#scatter plot
p1 = plt.scatter(data1['feature_1'], data1['feature_2'], marker = 'x', color = 'r', label = '1', s = 30)
p2 = plt.scatter(data0['feature_1'], data0['feature_2'], marker = '*', color = 'b', label = '0', s = 30)
plt.legend(loc = 'upper right')

#plot decision boundary
X = dataset[['feature_1','feature_2','feature_3', 'feature_4', 'feature_5','feature_6','feature_7','feature_8','feature_9',
             'feature_10','feature_11']]
Y = dataset.label

X_train, X_test, y_train, y_test = train_test_split(X, Y, test_size = .1, random_state=np.random.randint(10e5))

logreg = LogisticRegression(solver='lbfgs', multi_class='multinomial')
logreg.fit(X_train, y_train)
w2 = logreg.coef_
w1 = logreg.intercept_
w = np.append(w1, w2)

x_min, x_max = X_train[['feature_1']].min() - .5, X_train[['feature_1']].max() + .5
y_min, y_max = X_train[['feature_2']].min() - .5, X_train[['feature_2']].max() + .5
h = .02  # step size in the mesh
xx, yy = np.meshgrid(np.arange(x_min, x_max, h), np.arange(y_min, y_max, h))

def f(x1, x2, w):
    return w[0] + w[1]*x1 + w[2]*x2 + w[3]*x1**2 + w[4]*x1*x2 + w[5]*x2**2 + w[6]*x1**3 +w[7]*x1**2*x2 + w[8]*x2**2*x1 + w[9]*x2**3 + w[10]*x2*x1**2 + w[11]*x1*x2**2

plt.contour(xx, yy, f(xx, yy, w),0, colors='yellow',linewidths = 3)
plt.xlabel('Feature 1')
plt.ylabel('Feature 2')
plt.xticks(())
plt.yticks(())
plt.show()
```

# Problem 5

## (a)


```python
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import confusion_matrix
from sklearn.metrics import classification_report
from matplotlib import pyplot as plt
import pandas as pd
import numpy as np
```


```python
# input data
data = pd.read_csv("dataset_3.csv")
dataset = pd.DataFrame(data)

# rename the columns of dataset
dataset.columns = ['feature_1', 'feature_2','label']

#select subdataset by label 0, 1, 2, 3
data0 = dataset[dataset['label'] == 0]
data1 = dataset[dataset['label'] == 1]
data2 = dataset[dataset['label'] == 2]
data3 = dataset[dataset['label'] == 3]

#scatter plot
p0 = plt.scatter(data0['feature_1'], data0['feature_2'], marker = '*', color = 'b', label = '0', s = 30)
p1 = plt.scatter(data1['feature_1'], data1['feature_2'], marker = 'o', color = 'r', label = '1', s = 30)
p2 = plt.scatter(data2['feature_1'], data2['feature_2'], marker = 'x', color = 'm', label = '2', s = 30)
p3 = plt.scatter(data3['feature_1'], data3['feature_2'], marker = 's', color = 'g', label = '3', s = 30)
plt.legend(loc = 'lower left')
plt.xlabel('Feature 1')
plt.ylabel('Feature 2')
```

## (b)


```python
X = dataset[['feature_1','feature_2']]
Y = pd.DataFrame(dataset.label)

np.random.seed(20190313)

X_train, X_test, y_train, y_test = train_test_split(X, Y, test_size = .25, random_state=np.random.randint(10e5))

y_tr = y_train
y_te = y_test
```

## (c)


```python
# construct new training target vectors, making it binary
y_tr0 = y_tr.replace({1:4, 2:4, 3:4})
y_tr1 = y_tr.replace({0:4, 2:4, 3:4})
y_tr2 = y_tr.replace({0:4, 1:4, 3:4})
y_tr3 = y_tr.replace({0:4, 1:4, 2:4})

# construct new testing target vectors, making it binary
y_te0 = y_te.replace({1:4, 2:4, 3:4})
y_te1 = y_te.replace({0:4, 2:4, 3:4})
y_te2 = y_te.replace({0:4, 1:4, 3:4})
y_te3 = y_te.replace({0:4, 1:4, 2:4})

# fit the training data with logistic regression model
#y_tr0_new = column_or_1d(y_tr0, warn=True)
logreg0 = LogisticRegression(solver='lbfgs', multi_class='multinomial')
logreg0.fit(X_train, y_tr0.values.ravel())

logreg1 = LogisticRegression(solver='lbfgs', multi_class='multinomial')
logreg1.fit(X_train, y_tr1.values.ravel())

logreg2 = LogisticRegression(solver='lbfgs', multi_class='multinomial')
logreg2.fit(X_train, y_tr2.values.ravel())

logreg3 = LogisticRegression(solver='lbfgs', multi_class='multinomial')
logreg3.fit(X_train, y_tr3.values.ravel())
```




    LogisticRegression(C=1.0, class_weight=None, dual=False, fit_intercept=True,
              intercept_scaling=1, max_iter=100, multi_class='multinomial',
              n_jobs=1, penalty='l2', random_state=None, solver='lbfgs',
              tol=0.0001, verbose=0, warm_start=False)



## (d)


```python
# Use score method to get accuracy of model
score0 = logreg0.score(X_test, y_te0)
score1 = logreg1.score(X_test, y_te1)
score2 = logreg2.score(X_test, y_te2)
score3 = logreg3.score(X_test, y_te3)

avg_score = (score0 + score1 + score2 + score3) / 4 
print(score0, score1, score2, score3)
print(avg_score)
```

    0.97 0.75 0.72 0.99
    0.8574999999999999


For each class $y_{i} = 0, 1, 2, 3$, the accuracy of logistic regression model is $0.97\%$, $0.75\%$, $0.72\%$, $0.99\%$. The overall relative accuracy is $0.86\%$.  

## (e)


```python
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import confusion_matrix
from sklearn.metrics import classification_report
from matplotlib import pyplot as plt
import pandas as pd
import numpy as np

# input data
data = pd.read_csv("dataset_3.csv")
dataset = pd.DataFrame(data)

# rename the columns of dataset
dataset.columns = ['feature_1', 'feature_2','label']

#select subdataset by label 0 or 1
data0 = dataset[dataset['label'] == 0]
data1 = dataset[dataset['label'] == 1]
data2 = dataset[dataset['label'] == 2]
data3 = dataset[dataset['label'] == 3]

X = dataset[['feature_1','feature_2']]
Y = dataset.label

X_train, X_test, y_train, y_test = train_test_split(X, Y, test_size = .1, random_state=np.random.randint(10e5))

# Use Logistic Regression model and fit the data
logreg = LogisticRegression(solver='lbfgs', multi_class='multinomial')
logreg.fit(X_train, y_train)

#Get the predictive value of Y
y_pred = logreg.predict(X_test)

print(classification_report(y_test, y_pred))

# Plot the decision boundary. For that, we will assign a color to each point in the mesh [x_min, x_max]x[y_min, y_max].
X_train = pd.DataFrame(X_train)
y_train = pd.DataFrame(y_train)
X_train['y_train'] = y_train
new_data = X_train
new_data0 = new_data[new_data['y_train'] == 0]
new_data1 = new_data[new_data['y_train'] == 1]
new_data2 = new_data[new_data['y_train'] == 2]
new_data3 = new_data[new_data['y_train'] == 3]

x_min, x_max = new_data[['feature_1']].min() - .5, new_data[['feature_1']].max() + .5
y_min, y_max = new_data[['feature_2']].min() - .5, new_data[['feature_2']].max() + .5
h = .02  # step size in the mesh
x_min = np.asarray(x_min)
x_max = np.asarray(x_max)
y_min = np.asarray(y_min)
y_max = np.asarray(y_max)
xx, yy = np.meshgrid(np.arange(x_min, x_max, h), np.arange(y_min, y_max, h))
Z = logreg.predict(np.c_[xx.ravel(), yy.ravel()])

# Put the result into a color plot
Z = Z.reshape(xx.shape)
plt.contour(xx, yy, Z)

# Plot also the training points
p0 = plt.scatter(data0['feature_1'], data0['feature_2'], marker = '*', color = 'b', label = '0', s = 30)
p1 = plt.scatter(data1['feature_1'], data1['feature_2'], marker = 'o', color = 'r', label = '1', s = 30)
p2 = plt.scatter(data2['feature_1'], data2['feature_2'], marker = 'x', color = 'm', label = '2', s = 30)
p3 = plt.scatter(data3['feature_1'], data3['feature_2'], marker = 's', color = 'g', label = '3', s = 30)
plt.xlabel('Feature 1')
plt.ylabel('Feature 2')

plt.xlim(xx.min(), xx.max())
plt.ylim(yy.min(), yy.max())
plt.xticks(())
plt.yticks(())

plt.show()
```

01a30749e03ccd624ac1909467d39e7e34f222b4
ipynb
Jupyter Notebook
docs/Assignment 2 of STAT 641 _ Qian Wang.ipynb
Autumn-grass/QianWang
74af8b6b92c3e647e8f2cf60716c54a5e431adb1
MIT
docs/Assignment 2 of STAT 641 _ Qian Wang.ipynb
Autumn-grass/QianWang
74af8b6b92c3e647e8f2cf60716c54a5e431adb1
MIT
docs/Assignment 2 of STAT 641 _ Qian Wang.ipynb
Autumn-grass/QianWang
74af8b6b92c3e647e8f2cf60716c54a5e431adb1
MIT
Qwen/Qwen-72B
1. YES
2. YES