We will explore here the use of the Bayesian Information Criterion (BIC) for model selection, specifically to select features for a Naive Bayes model. The input consists of $d$-dimensional binary feature vectors $\mathbf{x}=\left[x_1, \ldots, x_d\right]^T$ where $x_i \in{-1,1}$ and binary class labels $y \in{-1,1}$.
A Naive Bayes model may be used, for example, to classify emails as spam $(y=1)$ or not-spam $(y=-1)$. The body of each email can be represented as a bag of words (i.e., we ignore frequency, placement, and grammar etc.). We could restrict ourselves to a dictionary of $d$ words $\left{w_1, \ldots, w_d\right}$. and coordinate $x_i$ could serve as an indicator of word $w_i: x_i=1$ if $w_i$ is present in the email and $x_i=-1$ if it is absent.
Recall that the Naive Bayes model assumes that the features are conditionally independent given the label so that
$$
P(\mathbf{x}, y)=\left[\prod_{i=1}^d P\left(x_i \mid y\right)\right] P(y)
$$
We parameterize our model by introducing separate parameters for each component:
$$
\begin{aligned}
P(\mathbf{x} \mid y, \theta) & =\left[\prod_{i=1}^d P\left(x_i \mid y, \theta_i\right)\right] \text { where } \
P\left(x_i \mid y, \theta_i\right) & =\theta_{i \mid y}^{\frac{x_i+1}{2}}\left(1-\theta_{i \mid y}^{\frac{1-x_i}{2}}\right)
\end{aligned}
$$
Here $\theta=\left[\theta_1, \ldots, \theta_d\right]^T$ and $\theta_i=\left[\theta_{i \mid 1}, \theta_{i \mid-1}\right]^T$ so that
$$
\begin{aligned}
\theta_{i \mid 1} & =P\left(x_i=1 \mid y=1\right) \
\theta_{i \mid-1} & =P\left(x_i=1 \mid y=-1\right)
\end{aligned}
$$
In class, we have already discussed how to compute the Maximum Likelihood estimates of these parameters. We will take a Bayesian approach, however, and introduce a prior probability over the parameters $\theta$ as:
$$
\begin{aligned}
P(\theta) & =\prod_{i=1}^d \prod_{y \in{1,-1}} P\left(\theta_{i \mid y}\right) \
P\left(\theta_{i \mid y}\right) & =\frac{\Gamma\left(r^{+}+r^{-}+2\right)}{\Gamma\left(r^{+}+1\right) \Gamma\left(r^{-}+1\right)} \cdot \theta_{i \mid y}^{r^{+}}\left(1-\theta_{i \mid y}\right)^{r^{-}}
\end{aligned}
$$

where $r^{+}$and $r^{-}$are hyper-parameters, common across all the $\theta_{i \mid y}$ ‘s, and, for integer $k, \Gamma(k+1)=k$ !. The hyper-parameters characterize our beliefs about the parameters prior to seeing any data. You may assume here that $r^{+}$and $r^{-}$are non-negative integers.
(a) Eqn 5 specifies the prior $P\left(\theta_{i \mid y}\right)$ as a Beta distribution. This choice of the prior is particularly convenient, as we will see now. Show that the posterior distribution $P(\theta \mid \mathcal{D})$, given $n$ training examples $\left{\left(\mathbf{x}j, y_j\right) \mid j=1, \ldots, n\right}$, has the following form: $$ P(\theta \mid \mathcal{D}) \propto L(\mathcal{D} ; \theta) P(\theta)=\prod{i=1}^d \prod_{y \in{-1,1}} \theta_{i \mid y}^{m_{i \mid y}^{+}}\left(1-\theta_{i \mid y}\right)^{m_{i \mid y}^{-}}
$$
where $L(\mathcal{D} ; \theta)$ is the likelihood function. What are $m_{i \mid y}^{+}$and $m_{i \mid y}^{-}$? Your answer should use the $\hat{n}y$ and $\hat{n}{i y}\left(x_i, y\right)$ notation used in the lectures; e.g., $\hat{n}_{k y}(1,-1)$ is the number of examples where $y=-1$ and $x_k=1$.

You have just shown that the Beta distribution is a conjugate prior to the Binomial distribution. In other words, if the prior probability $P(\theta)$ has the form of a Beta distribution, and the likelihood $P(\mathcal{D} \mid \theta)$ has the form of a Binomial distribution, then the posterior probability $P(\theta \mid \mathcal{D})$ will be a Beta distribution. When the class labels and features are not binary but take on $k$ values, the same relationship holds between the Dirichlet and Multinomial distributions.

(b) Recall that, in the Bayesian scheme for model selection, our aim is to find the model that maximizes the marginal likelihood, $P\left(\mathcal{D} \mid \mathcal{F}_d\right)$, which is the normalization constant for the posterior:
$$
P\left(\mathcal{D} \mid \mathcal{F}_d\right)=\int L(\mathcal{D} ; \theta) P(\theta) d \theta
$$
where $\mathcal{F}_d$ denotes the model. Compute a closed-form expression for $P\left(\mathcal{D} \mid \mathcal{F}_d\right)$. (Hint: use the form of the normalization constant for the beta distribution).

(c) Let us now use our results to choose between two models $\mathcal{F}1$ and $\mathcal{F}_2$. The only difference between them is in how they use feature $i: \mathcal{F}_2$ involves $P\left(x_i \mid y\right)$ term that connects the feature to the label $y$ whereas $\mathcal{F}_1$ only has $P\left(x_i\right)$, assuming that feature $x_i$ is independent of the label. In $\mathcal{F}_2$ the distribution of $x_i$ ‘s is parameterized by two parameters $\left[\theta{i \mid 1}, \theta_{i \mid-1}\right]^T$ as described before. In contrast, $\mathcal{F}_1$ only requires a single parameter $\theta_i^{\prime}$,
$$
P\left(x_i \mid \theta_i^{\prime}\right)=\theta_i^{\prime \frac{x_i+1}{2}}\left(1-\theta_i^{\prime}\right)^{\frac{1-x_i}{2}}
$$
The prior probability over $\theta_i^{\prime}$ is the same as Eqn 6 :
$$
P\left(\theta_i^{\prime}\right) \propto \theta_i^{r^{+}}\left(1-\theta_i^{\prime}\right)^{r^{-}}
$$
The expressions of the marginal likelihood from the two models will differ only in terms of feature $i$. To compare the two models we can ignore all the other terms. From your results in part (b), extract from $P\left(\mathcal{D} \mid \mathcal{F}_2\right)$ the factors involving feature $i$. Use a similar analysis to evaluate the terms involving feature $i$ in $P\left(\mathcal{D} \mid \mathcal{F}_1\right)$. Using these, evaluate the “decision rule”, $\log \left[P\left(\mathcal{D} \mid \mathcal{F}_2\right) / P\left(\mathcal{D} \mid \mathcal{F}_1\right)\right]>0$ to include feature $x_i$.

The data and scripts for this problem are available in hw2/prob1. You can load the data using the MATLAB script load_al_data. This script should load the matrices y_noisy, y_true, X_in. The $y$ vectors are $n \times 1$ while $\mathrm{X}{-}$in is a $n \times 3$ matrix with each row corresponding to a point in $\mathcal{R}^3$. The $y{\text {true }}$ vectors correspond to the ideal $y$ values, generated directly from the “true” model (whatever it may be) without any noise. In contrast, the $y_{\text {noisy }}$ vectors are the actual, noisy observations, generated by adding Gaussian noise to the $y_{\text {true }}$ vectors. You should use $y_{n o i s y}$ for any estimation. $y_{\text {true }}$ is provided only to make it easier to evaluate the error in your predictions (simulate an infinite test data). You would not have $y_{\text {true }}$ in any real task.
(a) Write MATLAB functions theta = linear_regress $(\mathrm{y}, \mathrm{X})$ and $\mathrm{y}$ hat $=$ linear_pred(theta, X_test). Note that we are not explicitly including the offset parameter but instead rely on the feature vectors to provide a constant component. See part (b).
(b) The feature mapping can substantially affect the regression results. We will consider two possible feature mappings:
$$
\begin{aligned}
& \phi_1\left(x_1, x_2, x_3\right)=\left[1, x_1, x_2, x_3\right]^T \
& \phi_2\left(x_1, x_2, x_3\right)=\left[1, \log x_1^2, \log x_2^2, \log x_3^2\right]^T
\end{aligned}
$$
Use the provided MATLAB function feature mapping to transform the input data matrix into a matrix
$$
X=\left[\begin{array}{c}
\phi\left(\mathbf{x}1\right)^T \ \phi\left(\mathbf{x}_2\right)^T \ \cdots \ \phi\left(\mathbf{x}{\mathbf{n}}\right)^T
\end{array}\right]
$$
For example, $\mathrm{X}$ = feature_mapping ( $\mathrm{X}_{-}$in, 1 ) would get you the first feature representation. Using your completed linear regression functions, compute the mean squared prediction error for each feature mapping (2 numbers).

(c) The selection of points to query in an active learning framework might depend on the feature representation. We will use the same selection criterion as in the lectures, the expected squared error in the parameters, proportional to $\operatorname{Tr}\left[\left(X^T X\right)^{-1}\right]$. Write a MATLAB function $\mathrm{idx}=\operatorname{active_ learn}(\mathrm{X}, \mathrm{k} 1, \mathrm{k} 2)$. Your function should assume that the top $k_1$ rows in $X$ have been queried and your goal is to sequentially find the indices of the next $k_2$ points to query. The final set of $k_1+k_2$ indices should be returned in idx. The latter may contain repeated entries. For each feature mapping, and $k_1=5$ and $k_2=10$, compute the set of points that should be queried (i.e., $\mathrm{X}(:, \mathrm{idx})$ ). For each set of points, use the feature mapping $\phi_2$ to perform regression and compute the resulting mean squared prediction errors (MSE) over the entire data set (again, using $\phi_2$ ).

(d) Let us repeat the steps of part (c) with randomly selected additional points to query. We have provided a MATLAB function $i d x=\operatorname{randomly}(\operatorname{select}(\mathrm{X}, \mathrm{k} 1, \mathrm{k} 2)$ which is essentially the same as active_learn except that it selects the $k_2$ points uniformly at random from $X$. Repeat the regression steps as in previous part, and compute the resulting mean squared prediction error again. To get a reasonable comparison you should repeat this process 50 times, and use the median MSE. Compare the resulting errors with the active learning strategies. What conclusions can you draw?

(e) Let us now compare the two sets of points chosen by active learning due to the different feature representations. We have provided a function plot_points(X,idx_r,idx_b) which will plot each row of $X$ as a point in $\mathbf{R}^3$. The points indexed by $i d x _r$ will be circled in red and those marked by idx_b will be circled (larger) in blue (some of the points indexed by idx_r and idx_b might be common). Plot the original data points using the indexes of the actively selected points based on the two feature representations. Also plot the same indexes using $\mathrm{X}$ from the second feature representation with its first constant column removed. In class, we saw an example where the active learning strategy chose points at the extrema of the available space. Can you see evidence of this in the two plots?

In sequential active learning we are free to choose the next training input $x_{n+1}$, here within $[1,4]$, for which we will then receive the corresponding noisy target $y_{n+1}$, sampled from the underlying model. Suppose we already have $\left{\left(x_1, y_1\right),\left(x_2, y_2\right), \ldots,\left(x_n, y_n\right)\right}$ and are trying to figure out which $x_{n+1}$ to select. The goal is to choose $x_{n+1}$ so as to help minimize the variance of the predictions $f\left(x ; \hat{w}_n\right)=\hat{w}_n x$, where $\hat{w}_n$ is the maximum likelihood estimate of the parameter $w$ based on the first $n$ training examples.
(a) What is the variance of $f\left(x ; \hat{w}_n\right)$ due to the noise in the training outputs as a function of $x, n$, and $\sigma^2$ given fixed (already chosen) inputs $x_1, \ldots, x_n$ ?

(b) Which $x_{n+1}$ would we choose (within $[1,4]$ ) if we were to next select $x$ with the maximum variance of $f\left(x ; \hat{w}_n\right)$ ?

(T/F -2 points) Since the variance of $f\left(x ; \hat{w}_n\right)$ only depends on $x$, $n$, and $\sigma^2$, we could equally well select the next point at random from $[1,4]$ and obtain the same reduction in the maximum variance.