density-estim.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "13c4a8f0-c420-4131-88fc-0dda57076c17",
   "metadata": {},
   "source": [
    "# Density estimation, inference about densities, and inverse problems"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ce3b47d9-94a9-424b-8464-ee2cc33963cb",
   "metadata": {},
   "source": [
    "## Background\n",
    "\n",
    "Suppose we observe $\\{X_j\\}_{j=1}^n$ i.i.d. $F$, where $F$ has density $f$ with respect\n",
    "to Lebesgue measure on the real line.\n",
    "What can we learn about $f$ from these data?\n",
    "\n",
    "Estimating $f$ can play a role in exploratory data analysis (EDA)\n",
    "as a graphical summary of the data set.\n",
    "In some contexts, more rigorous\n",
    "estimates and inferences about $f$ and properties of $f$ such as\n",
    "its value at a point $f(x_0)$, its derivative\n",
    "at a point $f'(x_0)$, a Sobolev norm of $f$ such as\n",
    "$\\|f\\|_S^2 = \\int (f^2 + f'^2 + f''^2)dx$,\n",
    "and the number and locations of modes of $f$,\n",
    "also are interesting.\n",
    "We will look at some approaches to estimating $f$, to finding confidence regious\n",
    "for $f$, and to testing hypotheses\n",
    "about $f$.\n",
    "We will not dwell on optimality considerations.\n",
    "\n",
    "### The Histogram and the Naive Estimator\n",
    "This section follows Silverman  (1990).\n",
    "\n",
    "Let $\\hat{F}_n$ denote the empirical cdf of the data $\\{X_j\\}_{j=1}^n$:\n",
    "\\begin{equation}\n",
    "    \\hat{F}_n(x) = \\frac{1}{n} \\sum_{j=1}^n 1_{X_j \\le x}.\n",
    "\\end{equation}\n",
    "Although $\\hat{F}_n$ is often a good estimator of $F$,\n",
    "$d \\hat{F}_n/dx$ is usually not a good estimator of $f = dF/dx$.\n",
    "The derivative of the empirical cdf is a sum of point masses\n",
    "at the observations.\n",
    "It usually is not an enlightening representation of the data.\n",
    "\n",
    "Suppose we have a collection of _class intervals_ or _bins_\n",
    "$\\{ \\mathcal{I}_k = (a_k, a_{k+1}] \\}_{k=1}^K$ such that every $X_j$ is in some $\\mathcal{I}_k$.\n",
    "(Choosing the intervals to be open on the  left and closed on the right is\n",
    "arbitrary; the essential point is that they be disjoint and that their\n",
    "union include all the data.)\n",
    "Let\n",
    "\\begin{equation}\n",
    "    w_k = a_{k+1} - a_k.\n",
    "\\end{equation}\n",
    "The _histogram_ of the data using these bins is\n",
    "\\begin{equation}\n",
    "    h(x) = \\frac{1}{n}\\sum_{k=1}^K \\frac{1}{w_k} 1_{x \\in \\mathcal{I}_k}\n",
    "    \\#\\left \\{ \\{X_j\\}_{j=1}^n \\cap \\mathcal{I}_k \\right \\}.\n",
    "\\end{equation}\n",
    "The histogram is an estimate of $f$.\n",
    "Its general appearance, including the number and locations of its modes and its smoothness,\n",
    "depends strongly on the locations and widths of the bins.\n",
    "It is blocky and discontinuous.\n",
    "If the bin widths and locations are chosen well, its performance--in\n",
    "the sense of convergence to $f$ in a norm as the sample size $n$\n",
    "grows--can be reasonable.\n",
    "\n",
    "Another estimate of $f$ derives from the definition of $f$ as the\n",
    "derivative of $F$:\n",
    "\\begin{equation}\n",
    "    f(x) = \\lim_{\\epsilon \\rightarrow 0} \\frac{F(x+\\epsilon) - F(x)}{\\epsilon}\n",
    "         = \\lim_{h \\rightarrow 0} \\frac{1}{2h} \\Pr \\{ x - h < X \\le x+h \\}\n",
    "\\end{equation}\n",
    "One could imagine estimating $f$ by picking a small value of $h$ and taking\n",
    "\\begin{eqnarray}\n",
    "    \\hat{f}_h(x) & \\equiv & \\frac{1}{2h} \\left ( \\hat{F}_n(x+h) - \\hat{F}_n(x-h) \\right ) \\nonumber \\\\\n",
    "               & = & \\frac{1}{2nh} \\sum_{j=1}^n 1_{x-h < X_j \\le x+h} \\nonumber \\\\\n",
    "               & = & \\frac{1}{n}\\sum_{j=1}^n \\frac{1}{h} K \\left ( \\frac{x - X_j}{h} \\right ),\n",
    "\\end{eqnarray}\n",
    "where $K(x) = \\frac{1}{2}\\times1_{-1 < x \\le 1}$.\n",
    "This is the _naive density estimate_.\n",
    "It amounts to estimating $f(x)$ by a superposition (sum) of boxcar functions\n",
    "centered at the observations, each with width $2h$ and area $1/n$.\n",
    "This sum is also blocky and discontinuous, but it avoids one of the arbitrary\n",
    "choices in constructing a histogram: the choice of locations of the bins.\n",
    "As $h \\rightarrow 0$, the naive estimate converges weakly to the sum of point masses\n",
    "at the data; for $h >0$, the naive estimator smooths the data.\n",
    "The tuning parameter $h$ is analogous to the bin width in a histogram.\n",
    "Larger values of $h$ give smoother density estimates.\n",
    "Whether \"smoother\" means \"better\" depends on the true density $f$;\n",
    "generally, there is a tradeoff between bias and variance: increasing the smoothness\n",
    "increases the bias but decreases the variance.\n",
    "\n",
    "It follows from the fact that $\\int_{-\\infty}^\\infty K(x) dx = 1$ that\n",
    "\\begin{eqnarray}\n",
    "    \\int_{-\\infty}^\\infty \\hat{f}(x) dx\n",
    "        &=&\n",
    "        \\frac{1}{n} \\sum_{j=1}^n \\frac{1}{h} \\int_{-\\infty}^\\infty K\n",
    "        \\left ( \\frac{x - X_j}{h} \\right )\n",
    "        \\nonumber \\\\\n",
    "        &=& \\frac{1}{n} \\sum_{j=1}^n 1 = 1.\n",
    "\\end{eqnarray}\n",
    "It follows from the fact that $K(x) \\ge 0$ that $\\hat{f} \\ge 0$ for all $x$.\n",
    "Thus $\\hat{f}$ is a probability density function."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "56734378-e877-4e7d-bf71-fc2b10cb3034",
   "metadata": {},
   "source": [
    "## Kernel estimates\n",
    "\n",
    "The two properties of the boxcar just mentioned---integrating to one and nonnegativity---hold \n",
    "whenever $K(x)$ is itself a probability density function, not just when $K$ is a \n",
    "unit-area boxcar function.\n",
    "Using a smoother _kernel_ function $K$, such as a Gaussian density,\n",
    "leads to a smoother estimate $\\hat{f}_K$.\n",
    "Estimates that are linear combinations of such kernel functions centered at the\n",
    "data are called _kernel density estimates_.\n",
    "We denote the kernel density estimate with bandwidth (smoothing parameter) $h$ by\n",
    "\\begin{equation}\n",
    "    \\hat{f}_h(x) = \\frac{1}{nh} \\sum_{j=1}^n K \\left ( \\frac{x - X_j}{h} \\right ).\n",
    "\\end{equation}\n",
    "The dependence of the estimate on the kernel is not evident in the notation---the kernel\n",
    "is understood from context.\n",
    "Kernels are always chosen to integrate to one, but there can be asymptotic advantages\n",
    "to kernels that are negative in places.\n",
    "The density estimates derived using such kernels can fail to be probability densities,\n",
    "because they can be negative for some values of $x$.\n",
    "Typically, $K$ is chosen to be a symmetric probability density function.\n",
    "\n",
    "There is a large body of literature on choosing $K$ and $h$ well, where \"well\"\n",
    "means that the estimate converges asymptotically as rapidly as possible in\n",
    "some suitable norm on probability density functions.\n",
    "The most common measure of performance is the mean integrated squared error (MISE):\n",
    "\\begin{eqnarray}\n",
    "    \\mbox{\\rm MISE}(\\hat{f}) & \\equiv &  \\mathbb{E} \\int (\\hat{f}(x) - f(x))^2 dx \\nonumber \\\\\n",
    "    &=& \\int \\mathbb{E} (\\hat{f}(x) - f(x))^2 dx \\nonumber \\\\\n",
    "    &=& \\int (\\mathbb{E} \\hat{f}(x) - f(x))^2 dx + \\int \\mathrm{Var}(\\hat{f}) dx.\n",
    "\\end{eqnarray}\n",
    "The MISE is sum of the integral of the squared pointwise bias of the estimate and the\n",
    "pointwise variance of the estimate.\n",
    "For kernel estimates,\n",
    "\\begin{equation}\n",
    "    \\mathbb{E} \\hat{f}(x) = \\int \\frac{1}{h} K \\left ( \\frac{x - y}{h} \\right ) f(y) dy,\n",
    "\\end{equation}\n",
    "and\n",
    "\\begin{equation}\n",
    "    \\mathrm{Var} \\hat{f}(x) = \\int \\frac{1}{h^2} K \\left ( \\frac{x-y}{h} \\right )^2 f(y) dy -\n",
    "    \\left [ \\frac{1}{h} \\int K \\left ( \\frac{x-y}{h} \\right ) f(y) dy \\right ]^2.\n",
    "\\end{equation}\n",
    "The expected value of $\\hat{f}$ is a smoothed version of $f$, the result of\n",
    "convolving $f$ with the scaled kernel.\n",
    "If $f$ is itself very smooth, smoothing it by convolution with the scaled kernel does not change\n",
    "its value much, and the bias of the kernel estimate is small.\n",
    "But in places where $f$ varies rapidly compared with the width of the scaled kernel,\n",
    "the local bias of the kernel estimate will be large.\n",
    "Note that the bias depends on the kernel function and the scale (bandwidth) $h$, not\n",
    "on the sample size.\n",
    "\n",
    "The two previous expressions for bias and variance rarely lead to tractable \n",
    "computations, but good approximations\n",
    "are available subject to some assumptions about $K$ and $f$.\n",
    "Suppose $K$ integrates to 1, is symmetric about zero so that $\\int xK(x)dx = 0$,\n",
    "and has nonzero finite second central moment $\\int x^2 K(x) dx = k_2 \\ne 0$,\n",
    "and that $f$ has as many continuous derivatives as needed.\n",
    "Then, to second order in $h$,\n",
    "\\begin{equation}\n",
    "    \\mathrm{Bias}_h(x) \\approx \\frac{1}{2} h^2 f\"(x) k_2.\n",
    "\\end{equation}\n",
    "(See Silverman  (1990), pp. 38ff.)\n",
    "Thus\n",
    "\\begin{equation}\n",
    "    \\int \\mathrm{Bias}_h^2 (x) dx \\approx \\frac{1}{4} h^4 k_2^2 \\int (f\")^2(x) dx,\n",
    "\\end{equation}\n",
    "confirming more quantitatively that (asymptotically)\n",
    "the integrated bias depends on the smoothness\n",
    "of $f$.\n",
    "A similar Taylor series approximation shows that to first order in $h^{-1}$,\n",
    "\\begin{equation}\n",
    "    \\mathrm{Var} \\hat{f}(x) \\approx n^{-1}h^{-1} f(x) \\int K^2(u) du,\n",
    "\\end{equation}\n",
    "so\n",
    "\\begin{equation}\n",
    "    \\int \\mathrm{Var} \\hat{f}(x) dx \\approx n^{-1}h^{-1} \\int K^2(u) du.\n",
    "\\end{equation}\n",
    "Thus shows that to reduce the integrated bias, one wants a narrow kernel,\n",
    "which must have large values to satisfy $\\int K = 1$,\n",
    "while to reduce the integrated variance, one wants the kernel to have small\n",
    "values, which requires it to be broad to satisfy $\\int K = 1$.\n",
    "\n",
    "By calculus one can show that the approximate MISE is minimized by choosing the\n",
    "bandwidth to be\n",
    "\\begin{equation}\n",
    "    h^* = n^{-1/5} k_2^{-2/5} \\left [ \\int K^2(u) du \\right ]^{1/5}\n",
    "        \\left [ \\int (f\")^2(x) dx \\right ]^{-1/5},\n",
    "\\end{equation}\n",
    "which depends on the unknown density $f$.\n",
    "Note that the (approximately) optimal bandwidth for MISE decreases with $n$ as\n",
    "$n^{-1/5}$.\n",
    "For the (approximately) optimal bandwidth $h^*$,\n",
    "\\begin{equation}\n",
    "    \\mbox{\\rm MISE} \\approx n^{-4/5} \\times 1.25 C(K) \\left [ \\int (f\")^2(x) dx \\right ]^{1/5},\n",
    "\\end{equation}\n",
    "where\n",
    "\\begin{equation}\n",
    "    C(K) = k_2^{2/5} \\left [ \\int K^2(u) du \\right ]^{4/5}\n",
    "\\end{equation}\n",
    "depends only on the kernel.\n",
    "The kernel that is approximately optimal for MISE thus has the smallest\n",
    "possible value of $C(K)$ subject to the restrictions on the moments of $K$.\n",
    "If we restrict attention to kernels that are probability density functions,\n",
    "the optimal kernel is the _Epanechnikov kernel_ $K_e(u)$\n",
    "\\begin{equation}\n",
    "    K_e(u) = \\frac{3}{4 \\sqrt{5}} (  1 - u^2/5 )_+.\n",
    "\\end{equation}\n",
    "This is the positive part of a parabola.\n",
    "\n",
    "One can define the relative efficiency of other kernels compared with the Epanechnikov\n",
    "kernel as the ratio of their values of $C(K)^{5/4}$.\n",
    "Other common kernels include Tukey's Biweight (suitably normalized, this is $\\frac{15}{16} (1-u^2)_+^2$),\n",
    "a triangular kernel, the rectangular kernel of the naive estimate, and the Gaussian density.\n",
    "Table 3.1 on p. 43 of Silverman  (1990) shows that there is not much\n",
    "variation in the efficiency: the rectangular kernel is worst, with an efficiency of about\n",
    "93\\%; the efficiency of the Gaussian is about 95\\%; the efficiency of the triangular\n",
    "kernel is about 99\\%; and the efficiency of the Biweight is over 99\\%.\n",
    "Thus the choice of kernel can reflect other concerns, such as desired properties of $\\hat{f}$\n",
    "(continuity, computational complexity, and so on).\n",
    "\n",
    "Choosing $h = h(n)$ is much more of a concern for the asymptotic behavior of the density\n",
    "estimate.\n",
    "To a large extent, choosing $h$ is a black art, but there are some automatic strategies\n",
    "that behave well subject to some assumptions.\n",
    "One of the most popular is least-squares cross-validation, which is a resampling method\n",
    "related to the jackknife.\n",
    "Here is a sketch of the method, following Silverman  (1990), pp. 48ff.\n",
    "\n",
    "The integrated squared error of a density estimate $\\hat{f}$ is\n",
    "\\begin{equation}\n",
    "    \\int (\\hat{f} - f)^2 dx = \\int \\hat{f}^2 dx - 2 \\int \\hat{f}f dx + \\int f^2 dx.\n",
    "\\end{equation}\n",
    "The last term does not involve the density estimate, so it is not in our control.\n",
    "Thus it is enough to try to minimize\n",
    "\\begin{equation}\n",
    "    R(\\hat{f}) \\equiv \\int \\hat{f}^2 dx - 2 \\int \\hat{f}f dx.\n",
    "\\end{equation}\n",
    "Cross validation estimates $R(\\hat{f}_h)$ from the data, and chooses $h$ to minimize\n",
    "the estimate.\n",
    "The first term in $R(\\hat{f})$ can be calculated explicitly from $\\hat{f}$.\n",
    "Estimating the second term is the crux of the method.\n",
    "By analogy to the jackknife, define the _leave one out_ kernel density estimate\n",
    "\\begin{equation}\n",
    "    \\hat{f}_{h,(i)} (x) = \\frac{1}{(n-1)h} \\sum_{j \\ne i} K \\left ( \\frac{x - X_j}{h} \\right ) .\n",
    "\\end{equation}\n",
    "Let\n",
    "\\begin{equation}\n",
    "    M_0(h) \\equiv \\int \\hat{f}_h^2 - \\frac{2}{n} \\sum_i \\hat{f}_{h,(i)}(X_i).\n",
    "\\end{equation}\n",
    "Let's compute the expected value of $M_0(h)$.\n",
    "First note that\n",
    "\\begin{eqnarray}\n",
    "    \\mathbb{E} \\frac{1}{n} \\sum_i \\hat{f}_{h,(i)} (X_i) & = & \\mathbb{E} \\hat{f}_{h,(1)}(X_1) \\nonumber \\\\\n",
    "    & = & \\mathbb{E} \\int \\hat{f}_{h,(1)}(x) f(x) dx \\nonumber \\\\\n",
    "    & = & \\mathbb{E} \\int \\hat{f}_{h}(x) f(x) dx.\n",
    "\\end{eqnarray}\n",
    "The last step uses the fact that the expected value of the kernel density estimate\n",
    "depends on $K$ and $h$ but not on the sample size.\n",
    "Thus\n",
    "\\begin{equation}\n",
    "    \\mathbb{E} R(\\hat{f}_h) = \\mathbb{E} M_0(h),\n",
    "\\end{equation}\n",
    "and $M_0(h)$ is an unbiased estimator of the ISE of $\\hat{f}_h$, less the term $\\int f^2$,\n",
    "which does not depend on $\\hat{f}_h$.\n",
    "Provided $M_0(h)$ is close to $\\mathbb{E} M_0(h)$, choosing $h$ to minimize $M_0(h)$ should\n",
    "select a good value of $h$ for minimizing the MISE of the estimate.\n",
    "The form of $M_0(h)$ is not computationally efficient; simplifications are possible,\n",
    "especially if $K(\\cdot)$ is symmetric.\n",
    "Moreover, if we use $n$ in place of $n-1$ in the denominators, we get a\n",
    "similar score function $M_1(h)$ that is easier to compute:\n",
    "\\begin{equation}\n",
    "    M_1(h) = \\frac{1}{n^2h} \\sum_i \\sum_j K^*\\left ( \\frac{X_i - X_j}{h} \\right ) +\n",
    "        \\frac{2}{nh} K(0),\n",
    "\\end{equation}\n",
    "where $K^*(u) = (K \\star K)(u) - 2 K(u)$ (here $\\star$ denotes convolution).\n",
    "The score function $M_1(h)$ can be computed very efficiently by Fourier methods; see\n",
    " section 3.5 of Silverman  (1990).\n",
    "\n",
    "A theorem due to Charles Stone (1984) justifies asymptotically choosing $h$ by\n",
    "cross validation using the score function $M_1$.\n",
    "Stone's theorem says that, subject to minor restrictions on $K$,\n",
    "the ratio of the integrated squared error choosing $h$ by minimizing $M_1$ to\n",
    "the integrated squared error for the best choice of $h$ given the sample $\\{X_j\\}$\n",
    "converges to 1 with probability 1 as $n \\rightarrow \\infty$.\n",
    "\n",
    "Cross validation tends to fail when the data have been discretized (binned),\n",
    "because the behavior of $M_1(h)$ at small $h$ is sensitive to rounding and\n",
    "discretization.\n",
    "It can be rescued sometimes by restricting the optimization to a range of\n",
    "values of $h$ that excludes very small values.\n",
    "\n",
    "The MISE (or an estimate of it) is but one of many possible score functions\n",
    "that could be used in a cross validation scheme.\n",
    "For example, one could use the log likelihood instead, which leads to\n",
    "maximizing the score function\n",
    "\\begin{equation}\n",
    "    \\mbox{CV}(h) = \\frac{1}{n} \\sum_{i=1}^n \\log \\hat{f}_{h,(i)} (X_i).\n",
    "\\end{equation}\n",
    "It turns out that under strong restrictions on $f$ and $K$,\n",
    "$-\\mbox{CV}(h)$ is (within a constant) an unbiased estimator of the\n",
    "Kullback-Leibler distance between\n",
    "$\\hat{f}_h$ and $f$:\n",
    "\\begin{equation}\n",
    "    I(f, \\hat{f}_h) \\equiv \\int f(x) \\log \\frac{f(x)}{\\hat{f}_h(x)} dx.\n",
    "\\end{equation}\n",
    "The score function $\\mbox{CV}(h)$ is not resistant.\n",
    "\n",
    "## Kernel estimates of multivariate densities\n",
    "\n",
    "This material is drawn from Chapter 4 of Silverman (1990).\n",
    "\n",
    "Let $\\{X_j\\}_{j=1}^n$ each take values in $\\mathbf{R}^d$, $d \\ge 1$.\n",
    "Let $K: \\mathbf{R}^d \\rightarrow \\mathbf{R}$ satisfy\n",
    "\\begin{equation}\n",
    "    \\int_{\\mathbf{R}^d} K(x) dx = 1.\n",
    "\\end{equation}\n",
    "Typically, the kernel $K$ is a radially symmetric probability distribution\n",
    "such as the standard multivariate normal, or the multivariate Epanechnikov\n",
    "kernel\n",
    "\\begin{equation}\n",
    "    K_e(x) \\equiv \\frac{d+2}{2c_d}(1 - \\|x\\|^2)_+,\n",
    "\\end{equation}\n",
    "where $c_d$ is the volume of the unit sphere in $\\mathbf{R}^d$:\n",
    "\\begin{equation}\n",
    "    c_d = \\int_{\\mathbf{R}^d} 1_{\\|x\\| < 1} dx.\n",
    "\\end{equation}\n",
    "The kernels\n",
    "\\begin{equation}\n",
    "    K_2 (x) \\equiv \\frac{3}{\\pi}(1 - \\|x\\|^2)_+^2\n",
    "\\end{equation}\n",
    "and\n",
    "\\begin{equation}\n",
    "    K_3(x) \\equiv \\frac{4}{\\pi}(1 - \\|x\\|^2)_+^3\n",
    "\\end{equation}\n",
    "have more derivatives than the Epanechnikov kernels, and thus produce\n",
    "smoother density estimates; also, they are easier to compute than\n",
    "the multivariate normal density.\n",
    "\n",
    "Given a multivariate kernel function, the multivariate kernel density estimate is\n",
    "\\begin{equation}\n",
    "    \\hat{f}_h(x) \\equiv \\frac{1}{nh^d} \\sum_{j=1}^n K \\left ( \\frac{x - X_j}{h} \\right ),\n",
    "\\end{equation}\n",
    "which is directly analogous to the univariate kernel density estimate.\n",
    "The kernel density estimate is a sum of \"bumps\" centered at the observations,\n",
    "each with mass $1/n$ and a common width that depends on a tuning parameter, the\n",
    "bandwidth $h$.\n",
    "The bandwidth $h$ is \"isotropic\" in that all coordinates are scaled in the\n",
    "same way.\n",
    "If the coordinates are incommensurable (e.g. , if the variances of different\n",
    "coordinates are radically different), it can help to transform the coordinate system\n",
    "before using the estimator, for example, by transforming so that the covariance\n",
    "matrix of the observations is the identity matrix.\n",
    "The estimate can then be transformed by the inverse change of variables to get\n",
    "the density estimate in the original coordinate system.\n",
    "This corresponds to the estimate\n",
    "\\begin{equation}\n",
    "    \\hat{f}_{h,S}(x) = \\frac{1}{\\sqrt{|S|} n h^d} \\sum_{j=1}^n\n",
    "        k \\left ( \\frac{\\|x - X_j \\|_{S^{-1}}^2}{h^2} \\right ),\n",
    "\\end{equation}\n",
    "where $\\|x\\|_{S^{-1}}^2 \\equiv x^T S^{-1} x$, and $k(\\|x\\|^2) = K(x)$.\n",
    "\n",
    "Most of the treatment of univariate kernel density estimates carries over,\n",
    "_mutatis mutandis_, to the multivariate case.\n",
    "For example, there is an optimal window width for minimizing the (approximate) MISE;\n",
    "it depends on the smoothness of the underlying density (through $\\int (\\nabla^2 f)^2$)\n",
    "and on the norm of the kernel and on the second moment of the kernel.\n",
    "Stone's theorem shows that choosing the bandwidth by cross validation using the\n",
    "score function\n",
    "\\begin{equation}\n",
    "    M_1(h) = \\frac{1}{n^2 h^d} \\sum_i \\sum_j K^* \\left ( \\frac{X_i - X_j}{h} \\right )\n",
    "    + \\frac{2}{n h^d}K(0)\n",
    "\\end{equation}\n",
    "is asymptotically optimal for MISE.\n",
    "\n",
    "### {The _Curse of Dimensionality_}\n",
    "The difficulty of density estimation grows very rapidly as the dimension of the\n",
    "sample space, $d$, increases.\n",
    "For example, to get relative mean squared error at 0 to be less than 0.1\n",
    "in estimating a multivariate normal density at zero using the optimal\n",
    "kernel requires $n = 4$ for $d=1$, $n=19$ for $d=2$,\n",
    "$n=768$ for $d=5$, and $n = 842,000$ for $d=10$\n",
    "(see Table 4.2 of Silverman , 1990).\n",
    "Partly, this is because of the behavior of the volume element in high dimensional\n",
    "spaces.\n",
    "\n",
    "Consider the unit sphere in dimension $d$.\n",
    "As $d$ grows, the volume of the sphere is increasingly concentrated in\n",
    "a thin shell near radius $r=1$.\n",
    "As a result, regions of low density can contribute substantially to the\n",
    "probability in higher dimensions, and regions of high density can\n",
    "remain unsampled even for relatively large sample sizes when $d$ is large.\n",
    "This makes details of the estimate matter increasingly as $d$ grows, and\n",
    "makes it harder to estimate the density even where it is large as $d$ grows.\n",
    "\n",
    "## {Nearest neighbor estimates}\n",
    "This section follows Silverman  (1990),  section 5.2.\n",
    "We start with the $d$-dimensional case.\n",
    "For any point $t \\in \\mathbf{R}^d$, define $r_k(t)$ to be the\n",
    "Euclidean distance from $t$ to the $k$th closest datum in\n",
    "the set $\\{X_j\\}_{j=1}^n$.\n",
    "Let $V_k(t)$ be the volume in $\\mathbf{R}^d$ of a sphere of radius $r_k(t)$:\n",
    "\\begin{equation}\n",
    "    V_k(t) = c_d r_k^d(t),\n",
    "\\end{equation}\n",
    "where as before $c_d$ is the volume of the unit ball in $\\mathbf{R}^d$.\n",
    "The _nearest neighbor density estimate_ is\n",
    "\\begin{equation}\n",
    "    \\hat{f}_k(t) \\equiv \\frac{k}{n V_k(t)} = \\frac{k}{n c_d r_k^d(t)}.\n",
    "\\end{equation}\n",
    "Usually $k$ is chosen to be small compared with $n$; $k \\approx \\sqrt{n}$\n",
    "is typical in dimension $d=1$.\n",
    "Larger values of $k$ produce smoother estimates, but the smoothness\n",
    "varies locally: the effective \"window\" is narrower where the local density\n",
    "of data is higher.\n",
    "\n",
    "Why does the recipe for the nearest neighbor estimate make sense?\n",
    "If the density at $t$ is $f(t)$, then in a sample of size $n$, we would\n",
    "expect there to be about $n f(t) V_k(t)$ observations\n",
    "in a small sphere of volume $V_k(t)$ centered at $t$.\n",
    "If we set the expected number equal to the observed number and solve for\n",
    "$f$, we get the nearest neighbor estimate:\n",
    "\\begin{equation}\n",
    "    \\left \\{ n \\hat{f}(t) V_k(t) = n \\hat{f}(t) c_d r_k^d(t) = k \\right \\} \\Rightarrow\n",
    "    \\left \\{ \\hat{f}(t) = \\frac{k}{n c_d r_k^d(t)} \\right \\}\n",
    "\\end{equation}\n",
    "At the point $t$, each datum within a distance $r_k(t)$ of $t$ contributes\n",
    "$1/(n c_d r_k^d(t))$ to the density estimate---as if the density estimate\n",
    "at $t$ were a kernel estimate with the kernel equal to the indicator function\n",
    "of the unit ball in $\\mathbf{R}^d$ divided by the volume of the ball\n",
    "(so the kernel integrates to 1), with bandwidth $r_k(t)$.\n",
    "Of course, this bandwidth depends on $t$ through $r_k(t)$,\n",
    "so the nearest neighbor estimate can be thought of as a kernel density\n",
    "estimate with spatially varying kernel width.\n",
    "(The kernel width depends on the point $t$ at which the estimate is\n",
    "sought, not just on the data $\\{X_j\\}$.\n",
    "This leads to some difficulties---see below.)\n",
    "\n",
    "Nearest neighbor estimates are not smooth: although $r_k(t)$ is continuous in\n",
    "$t$, its derivative fails to exist at points where two or more data\n",
    "are at distance $r_k(t)$ from $t$.\n",
    "Moreover, the nearest neighbor estimate is not itself a density.\n",
    "Consider what happens as $\\|t\\|$ grows.\n",
    "When $\\|t\\|$ is larger than $\\max_j \\| X_j \\|$,\n",
    "$r_k(t)$ grows linearly with $\\|t\\|$,\n",
    "so the density estimate falls off like $\\|t\\|^{-d}$---which has infinite\n",
    "integral.\n",
    "This rate of decay does not depend on how the tails of the sample\n",
    "fall off.\n",
    "\n",
    "Nearest neighbor estimates can be generalized to kernels more complicated\n",
    "than indicator functions.\n",
    "The generalized nearest neighbor estimate using kernel $K$ is\n",
    "\\begin{equation}\n",
    "    \\hat{f}(t) = \\frac{1}{n r_k^d(t)} \\sum_{j=1}^n K\\left ( \\frac{t - X_j}{r_k(t)} \\right ).\n",
    "\\end{equation}\n",
    "This reduces to the simple nearest neighborhood estimate when $K$ is the indicator\n",
    "of the unit ball, scaled to have integral 1.\n",
    "The tail behavior of the generalized nearest neighbor estimate depends on\n",
    "details of the kernel $K$.\n",
    "\n",
    "Even though at any fixed point, the nearest neighbor estimate is equivalent to\n",
    "a kernel estimate, it is a different kernel estimate at each point.\n",
    "The kernel estimate is a density because it is a linear combination of\n",
    "densities, with coefficients that sum to one.\n",
    "Just one \"bump\" is centered at each datum.\n",
    "In contrast, with the nearest neighbor estimate, a different bump is centered at\n",
    "each datum in finding the estimate for different values of $t$: the bandwidth\n",
    "associated with the\n",
    "contribution of the $j$th datum is a function of $t$, not just of $X_j$.\n",
    "\n",
    "### Variable Kernel Method\n",
    "In contrast, the variable kernel method allows the bandwidth associated with\n",
    "each _datum_ to be different, but holds those bandwidths fixed as\n",
    "$t$ varies.\n",
    "Let $d_{jk}$ be the distance from the $X_j$ to its $k$th nearest neighbor; _i.e._,\n",
    "$d_{jk} = r_{k+1}(X_j)$.\n",
    "Then the variable kernel estimate is\n",
    "\\begin{equation}\n",
    "    \\hat{f} = \\frac{1}{n} \\sum_{j=1}^n\n",
    "        \\frac{K \\left ( \\frac{t - X_j}{hd_{jk}} \\right )}{h^d d_{jk}^d}.\n",
    "\\end{equation}\n",
    "As $h$ or $k$ grows, the estimate gets smoother.\n",
    "This estimate centers one \"bump\" of mass $1/n$ at each datum, but the widths of the\n",
    "bumps depend on the local density of observations through $d_{jk}^{-1}$.\n",
    "Because of this, the estimate is itself a density if the basic kernel $K$ is a density.\n",
    "When the distance to the $k$th nearest neighbor is large, the width of the bump is large.\n",
    "Using $d_{jk}$ is an attempt to adapt the bandwidth to the height of the underlying\n",
    "density.\n",
    "However, there are better estimates of the local density to use to adjust the\n",
    "bandwidth.\n",
    "Of course, one can allow the bandwidth to vary in ways other than through\n",
    "$d_{jk}$.\n",
    "\n",
    "### Adaptive Kernel estimates\n",
    "This material is drawn from Silverman (1990,  section 5.3).\n",
    "Instead of using $d_{jk}^{-1}$ as (proportional to) an estimate of the local\n",
    "density for picking the bandwidth for the kernel centered at $X_j$, one\n",
    "could use a different density estimate.\n",
    "This approach leads to adaptive kernel estimates.\n",
    "\n",
    "The idea is to make a pilot density estimate, usually highly smoothed,\n",
    "and to base the bandwidth choice for the final estimate on the pilot.\n",
    "Let $\\tilde{f}(t)$ be a pilot density estimate for which $\\min_j \\tilde{f}(X_j) > 0$.\n",
    "Let $\\alpha \\in [0, 1]$.\n",
    "Let $g$ be the geometric mean of the values $\\{\\tilde{f}(X_j)\\}_{j=1}^n$:\n",
    "\\begin{equation}\n",
    "    g = \\left ( \\prod_{j=1}^n \\tilde{f}(X_j) \\right )^{1/n} = \\exp \\left \\{ \\frac{1}{n}\n",
    "        \\sum_{j=1}^n \\log \\tilde{f}(X_j) \\right \\}.\n",
    "\\end{equation}\n",
    "Define\n",
    "\\begin{equation}\n",
    "    \\lambda_j = \\left ( \\frac{\\tilde{f}(X_j)}{g} \\right )^{-\\alpha}.\n",
    "\\end{equation}\n",
    "The adaptive kernel estimate is\n",
    "\\begin{equation}\n",
    "    \\hat{f}(t) = \\frac{1}{n} \\sum_{j=1}^n \\frac{K\n",
    "    \\left ( \\frac{t-X_j}{h \\lambda_j} \\right )}{h^d \\lambda_j^d}.\n",
    "\\end{equation}\n",
    "Note that $\\lambda_j$ plays the role of $d_{jk}$ of the variable kernel method.\n",
    "The estimate depends on a number of tuning constants: $h$, $\\alpha$, and the\n",
    "tuning constants of the method used to derive the pilot estimate.\n",
    "The overall bandwidth $h$ plays the same role as before.\n",
    "The _sensitivity parameter_ $\\alpha$ controls how much the bandwidth varies\n",
    "as the pilot estimate varies---the rapidity of variation with $\\tilde{f}$.\n",
    "For $\\alpha=0$, the method becomes the ordinary kernel estimate.\n",
    "Silverman  (1990) says that \"there are good reasons for setting $\\alpha=1/2$.\"\n",
    "(See his section 5.3.3 and reference to Abramson, 1982.)\n",
    "\n",
    "Silverman (1990)\n",
    "says that the fine details and smoothness of the pilot estimate don't matter much\n",
    "for the final estimate, and recommends using an Epanechnikov kernel estimate\n",
    "with bandwidth chosen to perform well for a standard distribution, calibrated to\n",
    "have the same variance as the sample.\n",
    "He does not advocate using cross validation or other computationally intensive\n",
    "schemes for the pilot estimate.\n",
    "\n",
    "Silverman  reports simulation studies by Breiman, Meisel and Purcell (1977),\n",
    "showing that with the bandwidth chosen optimally, the adaptive\n",
    "kernel method performs remarkably better than the fixed kernel method,\n",
    "even for tame densities such as the normal.\n",
    "\n",
    "The overall smoothing parameter $h$ for the adaptive kernel estimate\n",
    "can be chosen by least squares cross validation.\n",
    "\n",
    "### Maximum Penalized Likelihood\n",
    "This section follows section 5.4 of Silverman  (1990);\n",
    "it is connected to an approach to solving\n",
    "inverse problems, which we will discuss later in the course.\n",
    "\n",
    "Recall that the MLE of the distribution $F$ is just a sum of point masses\n",
    "at the observations, each with mass $1/n$.\n",
    "This is not very satisfactory as a density estimate because it is so rough.\n",
    "The idea of maximum penalized likelihood is to give up some likelihood in\n",
    "favor of smoothness.\n",
    "We need a functional $R$ that assigns a finite positive number to some subset\n",
    "of all density functions.\n",
    "For example, we might take\n",
    "\\begin{equation}\n",
    "    R(g) = \\int (g\")^2 dt.\n",
    "\\end{equation}\n",
    "Let $\\mathcal{F}$ be the set of probability density functions for which $R$ is defined and finite.\n",
    "For a fixed positive number $\\lambda$ (the smoothing parameter),\n",
    "the penalized log-likelihood function of the density $g$ is\n",
    "\\begin{equation}\n",
    "    \\ell_\\lambda (g) =\\sum_{j=1}^n \\log g(X_j) - \\lambda R(g).\n",
    "\\end{equation}\n",
    "The maximum penalized likelihood density estimate $\\hat{f}$ is any density in\n",
    "$\\mathcal{F}$ for which\n",
    "\\begin{equation}\n",
    "    \\ell_\\lambda(\\hat{f}) \\ge \\ell_\\lambda(g) \\;\\; \\forall g \\in \\mathcal{F}.\n",
    "\\end{equation}\n",
    "The maximum penalized log likelihood estimate is an \"optimal\" compromise\n",
    "between maximizing the likelihood and being as smooth as possible\n",
    "(in the sense of minimizing\n",
    "$R$).\n",
    "Estimators that optimize a tradeoff between data fit and simplicity\n",
    "are quite common in many settings; they are called\n",
    "_regularized_ estimates.\n",
    "The functional $R$ is called the _regularization functional_ or \n",
    "_penalty functional_.\n",
    "Measures of fidelity to the data other than the likelihood are also common.\n",
    "\n",
    "Finding the maximum penalized likelihood estimator\n",
    "can be made more tractable numerically in a variety of ways,\n",
    "depending on the choice of $R$.\n",
    "For example, imposing the constraint $g > 0$ is easier if one works with the square-root\n",
    "of the density or with the logarithm of the density, although imposing the\n",
    "other part of the constraint $g \\in \\mathcal{F}$---that the density integrates to one---is\n",
    "harder then.\n",
    "Discrete approximations\n",
    "to the density (such as truncated expansions in orthogonal sets of functions)\n",
    "also can simplify the numerics of finding an approximate maximum\n",
    "penalized likelihood estimate.\n",
    "See Silverman  (1990) for references and more detail.\n",
    "\n",
    "The penalized maximum likelihood approach, using roughness penalties\n",
    "like those described, treats the underlying density as having homogeneous\n",
    "smoothness.\n",
    "We will talk more about maximum penalized likelihood in the context of\n",
    "nonparametric regression (function estimation).\n",
    "\n",
    "\n",
    "## {Confidence sets for densities with shape restrictions; lower confidence interval\n",
    "for the number of modes}\n",
    "Reference: Hengartner and Stark (1995).\n",
    "\n",
    "## {Wavelet shrinkage}\n",
    "Reference: Kerkyacharian and Picard (1993).\n",
    "\n",
    "### Time-frequency localization: windowed Fourier transform and wavelets\n",
    "\n",
    "### Haar wavelets\n",
    "\n",
    "### Unconditional bases\n",
    "\n",
    "### {Hard and soft thresholding}\n",
    "\n",
    "## Inverse Problems\n",
    "Reference: Evans and Stark (2002).\n",
    "\n",
    "### Nonparametric regression\n",
    "\n",
    "### {Example: Abel's problem}\n",
    "Given a frictionless bowling ball with mass $m$, a stopwatch and the ability to roll the\n",
    "ball with any desired initial velocity $v_j = v_j(0)$, find\n",
    "the shape of a (2-dimensional)\n",
    "hill by rolling the ball with different velocities and measuring how long it\n",
    "takes the ball to return.\n",
    "The measurements have errors.\n",
    "You can think of this as a way to survey San Francisco on a foggy day.\n",
    "\n",
    "The _forward problem_ is to predict how long it takes the ball\n",
    "to return, if we know the shape of the hill.\n",
    "The _inverse problem_ is to use a finite set of measurements to\n",
    "learn something about the shape of the hill.\n",
    "\n",
    "Let's solve the forward problem.\n",
    "It is convenient to express the shape of the hill as the\n",
    "height $h(s)$ of the hill at an arc distance $s$ along the\n",
    "surface of the hill from where the ball is launched.\n",
    "The initial kinetic energy of the ball is\n",
    "\\begin{equation}\n",
    "    E_j = m v_j^2(0)/2.\n",
    "\\end{equation}\n",
    "As the ball ascends, its energy is conserved (the ball is frictionless),\n",
    "but it is partitioned into a kinetic component and a potential\n",
    "component.\n",
    "The potential energy component at arc distance $s$ is\n",
    "\\begin{equation}\n",
    "    E_{Pj}(s) = gmh(s),\n",
    "\\end{equation}\n",
    "so, by conservation of energy, the kinetic energy at arc\n",
    "distance $s$ is\n",
    "\\begin{equation}\n",
    "    E_{Kj}(s) = m v_j^2(0)/2 - gmh(s).\n",
    "\\end{equation}\n",
    "We can find the velocity of the ball at arc distance $s$ as follows\n",
    "\\begin{eqnarray}\n",
    "    mv_j^2(s)/2 & = & m v_j^2(0)/2 - gmh(s) \\nonumber \\\\\n",
    "    v_j^2(s) & = & v_j^2(0) - 2gh(s) \\nonumber \\\\\n",
    "    v_j(s) & = & \\sqrt{v_j^2(0) - 2gh(s)}.\n",
    "\\end{eqnarray}\n",
    "The velocity of the ball goes to zero (and the ball starts to come back) when\n",
    "$v_j^2(0) = 2gh(s)$, provided the slope of the hill does not vanish\n",
    "there (then the ball would balance and never return).\n",
    "Let $s_j$ satisfy $v_j^2(0) = 2gh(s_j)$.\n",
    "The time it takes the ball to return is equal to the time it takes\n",
    "the ball to ascend.\n",
    "The time it takes the ball to come back is thus\n",
    "\\begin{equation}\n",
    "    \\tau_j = 2\\int_{s=0}^{s_j} \\frac{ds}{\\sqrt{v_j^2(0) - 2gh(s)}}.\n",
    "\\end{equation}\n",
    "This is the solution to the forward problem.\n",
    "Each transit time $\\tau_j$ is a nonlinear functional of the\n",
    "hill profile $h(s)$.\n",
    "The inverse problem is to learn something about $h(s)$ from measurements\n",
    "\\begin{equation}\n",
    "    d_j = \\tau_j + \\epsilon_j, \\;\\; j = 1, \\ldots, n,\n",
    "\\end{equation}\n",
    "where $\\{\\epsilon_j\\}_{j=1}^n$ are stochastic errors whose joint\n",
    "distribution is assumed to be known---at least up to a parameter or two.\n",
    "(Rarely does anybody allow the joint distribution to be more\n",
    "general than a multivariate zero mean Gaussian with independent components\n",
    "whose variances are known.)\n",
    "\n",
    "It turns out that this stylized surveying problem is related to inverse\n",
    "problems in seismology, helioseismology, and stereology.\n",
    "\n",
    "### {General framework for inverse problems}\n",
    "Observe data $X$ drawn from a distribution $\\Pr_\\theta$ where\n",
    "$\\theta$ is unknown, but it is known that $\\theta \\in \\Theta$.\n",
    "Use $X$ and the constraint $\\theta \\in \\Theta$ to learn about\n",
    "$\\theta$.\n",
    "For example, we might want to estimate a parameter $g(\\theta)$.\n",
    "Assume that $\\Theta$ contains at least two points; otherwise, we\n",
    "know $\\theta$ and $g(\\theta)$ perfectly even without data.\n",
    "\n",
    "The parameter $g(\\theta)$ is _identifiable_ if\n",
    "\\begin{equation}\n",
    "    \\{ g(\\theta) \\ne g(\\eta) \\} \\Leftarrow \\{ \\Pr_\\theta \\ne \\Pr_\\eta\\},\\;\\;\n",
    "    \\forall \\theta, \\eta \\in \\Theta.\n",
    "\\end{equation}\n",
    "In most inverse problems, $\\theta$ is not identifiable.\n",
    "Little general is known about nonlinear inverse problems,\n",
    "although there are particular nonlinear inverse problems, like the\n",
    "surveying problem above, that are well understood.\n",
    "\n",
    "(The \"trick\" to solving the surveying problem is to work with\n",
    "$s(h)$ instead of $h(s)$, on the assumption that $h(s)$ is strictly\n",
    "monotonic.\n",
    "Then the forward mapping $\\tau$ is a linear functional of $s(h)$,\n",
    "but there are nonlinear constraints--$s(h)$ must also be monotonic.\n",
    "Although hills in San Francisco are not monotonic, in the seismic\n",
    "problem, there are thermodynamic arguments that the corresponding\n",
    "quantity--seismic velocity as a function of radius, divided by radius--is\n",
    "monotonic in Earth's core.\n",
    "See, e.g. \n",
    "Stark, P.B., 1992.  Inference in infinite-dimensional inverse problems:\n",
    "discretization and duality, _J. Geophys. Res._, _97_,\n",
    "14,055--14,082.)\n",
    "\n",
    "### Linear forward and inverse problems\n",
    "When the forward problem has more structure, more can be said.\n",
    "The best studied class of inverse problems are linear inverse problems.\n",
    "\n",
    "A forward problem is _linear_ if the constraint set $\\Theta$ is\n",
    "a subset of a separable Banach space $\\mathcal{T}$ and for some collection $\\{\\kappa_j\\}_{j=1}^n$\n",
    "of bounded linear functionals on $\\mathcal{T}$,\n",
    "\\begin{equation}\n",
    "    X_j = \\kappa_j \\theta + \\epsilon_j,\n",
    "\\end{equation}\n",
    "where $\\{\\epsilon_j\\}_{j=1}^n$ are random errors whose distribution does not depend on\n",
    "$\\theta$.\n",
    "Usually, such a forward problem is written\n",
    "\\begin{equation}\n",
    "    X = K\\theta + \\epsilon, \\;\\; \\theta \\in \\Theta.\n",
    "\\end{equation}\n",
    "The corresponding _linear inverse problem_ is to use the data $X$ and the constraint\n",
    "$\\theta \\in \\Theta$ to learn about $g(\\theta)$.\n",
    "In a linear inverse problem, the distribution of $X$ depends on $\\theta$ through\n",
    "$K\\theta$, so if there exist $\\theta, \\eta \\in \\Theta$ such that\n",
    "\\begin{equation}\n",
    "    K\\theta = K\\eta \\;\\mbox{ but }\\; g(\\theta) \\ne g(\\eta)\n",
    "\\end{equation}\n",
    "then $g(\\theta)$ is not identifiable.\n",
    "\n",
    "Let's simplify the setup even further---we assume that $\\mathcal{T}$ is a Hilbert space,\n",
    "that $\\Theta = \\mathcal{T}$, and that $g$ is a _linear parameter_; that is,\n",
    "\\begin{equation}\n",
    "    g(a \\theta + b \\eta) = a g(\\theta) + b g(\\eta)\n",
    "\\end{equation}\n",
    "for all $a, b \\in \\mathbf{R}$ and all $\\theta, \\eta \\in \\Theta$.\n",
    "The fundamental theorem of Backus and Gilbert says that then\n",
    "$g(\\theta)$ is identifiable if and only if $g = \\sum_{j=1}^n a_j \\kappa_j$\n",
    "for some constants $\\{a_j\\}$.\n",
    "In that case, if $\\mathbb{E} \\epsilon = 0$,\n",
    "$\\sum_j a_j X_j$ is unbiased for $g(\\theta)$, and if $e$ has covariance\n",
    "matrix $\\Sigma$, the MSE of this (linear) estimator is $a\\cdot \\Sigma \\cdot a^T$.\n",
    "See Evans and Stark (2002) for more details and proofs."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b7f7500a-50b5-433b-a0ce-ab92b21d7ed5",
   "metadata": {},
   "source": [
    "## Methods for inverse problems\n",
    "\n",
    "There is a huge number of methods for \"solving\" inverse problems, although\n",
    "what qualifies as a solution is debatable.\n",
    "These solution methods can be analyzed using traditional statistical measures\n",
    "of performance, including bias and various loss criteria.\n",
    "Perhaps the most important message is that without constraints, little\n",
    "can be said.  The issue is finding constraints justified by the science\n",
    "of the situation that still are helpful in reducing the uncertainty.\n",
    "\n",
    "+ Backus-Gilbert estimation.  Finding linear functionals close, in some sense,\n",
    "to point evaluators.\n",
    "\n",
    "+ MLE and variants (regularization, maximum penalized likelihood, method of sieves,\n",
    "singular value truncation and weighting).\n",
    "Trading off fidelity to the data and a measure of complexity or roughness.\n",
    "With suitable assumptions, can show consistency and good rates of convergence\n",
    "if the tradeoff is tuned appropriately.\n",
    "\n",
    "+ Bayes estimation.  Difficult to justify in infinite-dimensional problems,\n",
    "because prior probability distributions on infinite-dimensional spaces are\n",
    "strange---it's hard to capture constraints without injecting lots of additional\n",
    "information.\n",
    "\n",
    "+ Minimax estimation.  Interesting papers by Donoho and others on estimating\n",
    "linear functionals or the entire model in the Hilbert space case.\n",
    "Connection between deterministic optimal recovery and minimax statistical\n",
    "estimation in the case that the errors are Gaussian, $\\Theta$ is convex,\n",
    "and the parameter is a linear functional.\n",
    "\n",
    "+ Shrinkage estimation.  Shrinkage can improve MSE of estimates of\n",
    "high-dimensional means.  Can help with multiple Backus-Gilbert estimates.\n",
    "\n",
    "+ Wavelet-vaguelette shrinkage estimation. Analogue of wavelet shrinkage\n",
    "density estimation we looked at earlier.  Can outperform any linear method\n",
    "in some problems.  Papers by Donoho, Johnstone, and others.  Key idea is that\n",
    "the wavelet-vaguelette decomposition almost diagonalizes both the prior information\n",
    "and the forward problem.\n",
    "\n",
    "+ Strict bounds.  Analog of the method for confidence bounds on shape-restricted\n",
    "densities we looked at earlier in the class. Can get conservative joint confidence\n",
    "sets for arbitrarily many parameters of the model by finding upper and lower bounds\n",
    "on functionals over a set of models that satisfies the constraints and is in an\n",
    "infinite-dimensional confidence set based on fit to the data. Not generally\n",
    "optimal for standard measures of misfit to the data."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}