{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Biostat M280 Homework 5\n", "\n", "**Due May 31 @ 11:59PM**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this homework, we build a classifer for handwritten digit recognition. Following figure shows example bitmaps of handwritten digits from U.S. postal envelopes. \n", "\n", "\n", "\n", "Each digit is represented by a $32 \\times 32$ bitmap in which each element indicates one pixel with a value of white or black. Each $32 \\times 32$ bitmap is divided into blocks of $4 \\times 4$, and the number of white pixels are counted in each block. Therefore each handwritten digit is summarized by a vector $\\mathbf{x} = (x_1, \\ldots, x_{64})$ of length 64 where each element is a count between 0 and 16. \n", "\n", "We will use a model-based method by assuming a distribution on the count vector and carry out classification using probabilities. A common distribution for count vectors is the multinomial distribution. However as you will see in Q10, it is not a good model for handwritten digits. Let's work on a more flexible model for count vectors. In the Dirichlet-multinomial model, we assume the multinomial probabilities $\\mathbf{p} = (p_1,\\ldots, p_d)$ follow a Dirichlet distribution with parameter vector $\\alpha = (\\alpha_1,\\ldots, \\alpha_d)$, $\\alpha_j>0$, and density\n", "$$\n", "\\begin{eqnarray*}\n", "\t\\pi(\\mathbf{p}) = \\frac{\\Gamma(|\\alpha|)}{\\prod_{j=1}^d \\Gamma(\\alpha_j)} \\prod_{j=1}^d p_j^{\\alpha_j-1},\n", "\\end{eqnarray*} \n", "$$\n", "where $|\\alpha|=\\sum_{j=1}^d \\alpha_j$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Q1\n", "\n", "For a multivariate count vector $\\mathbf{x}=(x_1,\\ldots,x_d)$ with batch size $|\\mathbf{x}|=\\sum_{j=1}^d x_j$, show that the probability mass function for Dirichlet-multinomial distribution is\n", "$$\n", " f(\\mathbf{x} \\mid \\alpha) \n", "\t= \\int_{\\Delta_d} \\binom{|\\mathbf{x}|}{\\mathbf{x}} \\prod_{j=1}^d p_j^{x_j} \\pi(\\mathbf{p}) \\, d \\mathbf{p} \n", " = \\binom{|\\mathbf{x}|}{\\mathbf{x}} \\frac{\\prod_{j=1}^d \\Gamma(\\alpha_j+x_j)}{\\prod_{j=1}^d \\Gamma(\\alpha_j)} \\frac{\\Gamma(|\\alpha|)}{\\Gamma(|\\alpha|+|\\mathbf{x}|)}\n", "$$\n", "where $\\Delta_d$ is the unit simplex in $d$ dimensions and $|\\alpha| = \\sum_{j=1}^d \\alpha_j$.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Q2\n", "\n", "Given independent data points $\\mathbf{x}_1, \\ldots, \\mathbf{x}_n$, show that the log-likelihood is\n", "$$\n", "L(\\alpha) = \\sum_{i=1}^n \\ln \\binom{|\\mathbf{x}_i|}{\\mathbf{x}_i} + \\sum_{i=1}^n \\sum_{j=1}^d [\\ln \\Gamma(\\alpha_j + x_{ij}) - \\ln \\Gamma(\\alpha_j)] - \\sum_{i=1}^n [\\ln \\Gamma(|\\alpha|+|\\mathbf{x}_i|) - \\ln \\Gamma(|\\alpha|)].\n", "$$\n", "Is the log-likelihood a concave function?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Q3\n", "\n", "Write Julia function to compute the log-density of the Dirichlet-multinomial distribution." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "dirmult_logpdf" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "\"\"\"\n", " dirmult_logpdf(x::Vector, α::Vector)\n", " \n", "Compute the log-pdf of Dirichlet-multinomial distribution with parameter `α` \n", "at data point `x`.\n", "\"\"\"\n", "function dirmult_logpdf(x::Vector, α::Vector)\n", " # your code here\n", "end\n", "\n", "function dirmult_logpdf!(r::Vector, X::Matrix, α::Vector)\n", " for j in 1:size(X, 2)\n", " r[j] = dirmult_logpdf(X[:, j], α)\n", " end\n", " return r\n", "end\n", "\n", "\"\"\"\n", " dirmult_logpdf(X, α)\n", " \n", "Compute the log-pdf of Dirichlet-multinomial distribution with parameter `α` \n", "at each data point in `X`. Each column of `X` is one data point.\n", "\"\"\"\n", "function dirmult_logpdf(X::Matrix, α::Vector)\n", " r = zeros(size(X, 2))\n", " dirmult_logpdf!(r, X, α)\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Q4\n", "\n", "Read in [`optdigits.tra`](http://hua-zhou.github.io/teaching/biostatm280-2019spring/hw/hw5/optdigits.tra), the training set of 3823 handwritten digits. Each row contains the 64 counts of a digit and the last element (65th element) indicates what digit it is. For grading purpose, evaluate the total log-likelihood of this data at parameter values $\\alpha=(1,\\ldots,1)$ using your function in Q3." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Q5\n", "\n", "Derive the score function $\\nabla L(\\alpha)$, observed information matrix $-\\nabla^2L(\\alpha)$, and Fisher information matrix $\\mathbf{E}[-\\nabla^2L(\\alpha)]$ for the Dirichlet-multinomial distribution.\n", "\n", "Explain why Fisher scoring method is inefficient for computing MLE in this example." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Q6\n", "\n", "What structure does the observed information matrix possess that can facilitate the evaluation of the Newton direction? Is the observed information matrix always positive definite? What remedy can we take if it fails to be positive definite? (Hint: HW1 Q5.)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Q7\n", "\n", "Discuss how to choose a good starting point. Implement this as the default starting value in your function below. (Hint: Method of moment estimator may furnish a good starting point.)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Q8\n", "\n", "Write a function for finding MLE of Dirichlet-multinomial distribution given iid observations $\\mathbf{x}_1,\\ldots,\\mathbf{x}_n$, using the Newton's method. " ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "dirmult_newton" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "\"\"\"\n", " dirmult_newton(X)\n", "\n", "Find the MLE of Dirichlet-multinomial distribution using Newton's method.\n", "\n", "# Argument\n", "* `X`: an `n`-by-`d` matrix of counts; each column is one data point.\n", "\n", "# Optional argument \n", "* `alpha0`: a `d` vector of starting point (optional). \n", "* `maxiters`: the maximum allowable Newton iterations (default 100). \n", "* `tolfun`: the tolerance for relative change in objective values (default 1e-6). \n", "\n", "# Output\n", "* `maximum`: the log-likelihood at MLE. \n", "* `estimate`: the MLE. \n", "* `gradient`: the gradient at MLE. \n", "* `hessian`: the Hessian at MLE. \n", "* `se`: a `d` vector of standard errors. \n", "* `iterations`: the number of iterations performed.\n", "\"\"\"\n", "function dirmult_newton(X::Matrix; α0::Union{Vector, Nothing} = nothing, \n", " maxiters::Int = 100, tolfun::Float64 = 1e-6)\n", " \n", " # if α0==nothing, set default starting point from Q7\n", " \n", " # Newton loop\n", " for iter in 1:maxiters\n", " # evaluate gradient (score)\n", " \n", " # compute Newton's direction\n", " \n", " # line search loop\n", " for lsiter in 1:10\n", " # step halving\n", " end\n", " \n", " # check convergence criterion\n", " if abs(logl - loglold) < tolfun * (abs(loglold) + 1)\n", " break;\n", " end\n", " end\n", " \n", " # compute logl, gradient, Hessian from final iterate\n", " \n", " # output\n", " \n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Q9\n", "\n", "Read in [`optdigits.tra`](http://hua-zhou.github.io/teaching/biostatm280-2019spring/hw/hw5/optdigits.tra), the training set of 3823 handwritten digits. Find the MLE for the subset of digit 0, digit 1, ..., and digit 9 separately using your function." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Q10\n", "\n", "As $|\\alpha| \\to \\infty$ and $\\alpha / |\\alpha| \\to \\mathbf{p}$, the Dirichlet-multinomial distribution converges to a multinomial with parameter $\\mathbf{p}$. Therefore multinomial can be considered as a special Dirichlet-multinomial with $|\\alpha|=\\infty$. Perform a likelihood ratio test (LRT) whether Dirichlet-multinomial offers a better fit than multinomial for digits 0, 1, ..., 9 respectively." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Q11\n", "\n", "Now we can construct a simple Bayesian rule for handwritten digits recognition:\n", "$$\n", "\t\\mathbf{x}\t\\mapsto \\arg \\max_k \\widehat \\pi_k f(x|\\widehat \\alpha_k).\n", "$$\n", "Here we can use the proportion of digit $k$ in the training set as the prior probability $\\widehat \\pi_k$. Report the performance of your classifier on the test set of 1797 digits in [`optdigits.tes`](http://hua-zhou.github.io/teaching/biostatm280-2019spring/hw/hw5/optdigits.tes)." ] } ], "metadata": { "@webio": { "lastCommId": null, "lastKernelId": null }, "kernelspec": { "display_name": "Julia 1.1.0", "language": "julia", "name": "julia-1.1" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.1.0" }, "toc": { "colors": { "hover_highlight": "#DAA520", "running_highlight": "#FF0000", "selected_highlight": "#FFD700" }, "moveMenuLeft": true, "nav_menu": { "height": "65px", "width": "252px" }, "navigate_menu": true, "number_sections": false, "sideBar": true, "skip_h1_title": true, "threshold": 4, "toc_cell": false, "toc_section_display": "block", "toc_window_display": false, "widenNotebook": false } }, "nbformat": 4, "nbformat_minor": 2 }