% This is a simple LaTex sample document that gives a submission format
%   for IEEE PAMI-TC conference submissions.  Use at your own risk.

% Make two column format for LaTex 2e.\
\documentclass[12pt,twocolumn]{article} %,twocolumn
\usepackage{times,amsmath,amsfonts}

% Use following instead for LaTex 2.09 (may need some other mods as well).
%\documentstyle[times,twocolumn]{article}
\usepackage[dvips]{graphicx,graphics}
% Set dimensions of columns, gap between columns, and paragraph indent
\setlength{\textheight}{10in} \setlength{\textwidth}{7in}
%\setlength{\columnsep}{0.3125in} \setlength{\topmargin}{0in}
\setlength{\headheight}{0in} \setlength{\headsep}{-1in}
\setlength{\parindent}{1pc}
\setlength{\oddsidemargin}{-.5in}  % Centers text.
\setlength{\evensidemargin}{-.5in}

% Add the period after section numbers.  Adjust spacing.
\newcommand{\Section}[1]{\vspace{-8pt}\section{\hskip -1em.~~#1}\vspace{-3pt}}
\newcommand{\SubSection}[1]{\vspace{-3pt}\subsection{\hskip -1em.~~#1}
        \vspace{-3pt}}
\newcommand{\bqn}{\begin{eqnarray}}
\newcommand{\eqn}{\end{eqnarray}}
\newcommand {\diff}[1] {\frac{\partial}{\partial #1}}
\newcommand{\jacob}[3]{\frac{\partial^2 #3}{\partial #1 \partial #2}}
\newcommand{\der}[2]{\frac{\partial #2}{\partial #1}}
\begin{document}

% Make title bold and 14 pt font (Latex default is non-bold, 16pt)
\title{Stat 471: Lecture 11\\
Importance Sampling.}
% For single author (just remove % characters)
\author{Moo K. Chung\\
mchung@stat.wisc.edu}
% For two authors (default example)
\maketitle \thispagestyle{empty}

\begin{enumerate}
\item Estimate
$$p=\int_{-\infty}^{-2} \frac{e^{-x^2/2}}{\sqrt{2\pi}} \; dx.$$
We use the formula $$\mathbb{E} g(X) = \int g(x)f(x) \; dx$$ with
$g(x)$ to be an index function such that $g(x)= 1$ if $x \in
(-\infty,-2)$ and 0 otherwise. Let $X \sim N(0,1)$. Then $\hat p =
\sum_{i=1}^ng(X_i)/n$.

\begin{verbatim}
n=1000000
X=snorm(n);
index=find(X<-2);
g=zeros(n,1);
g(index)=1;
>>[mean(g) var(g)]
ans =
    0.0228    0.0223
>>normcdf(-2)
ans =
    0.0228
\end{verbatim}
\item An alternate approach to estimate $p$ would be to use the
method of {\em importance sampling}
$$\mathbb{E} g(X)= \int g(x)\frac{f(x)}{h(x)} h(x)\; dx =
\mathbb{E} \Big[g(X)\frac{f(X)}{h(X)}\Big]$$ for some distribution
$h$. A more proper name for this technique would be {\em weighted
sampling}. Then we estimate the integral with $\hat p=\sum_{i=1}^n
w(X_i)g(X_i)/n$ with the importance sampling weights $w(X_i)
=f(X_i)/h(X_i)$.

\item The aim of importance sampling is to sample more frequently
from values by weighting more so that the contribution to the
integral is the greatest (We will probably talk about {\em
importance resampling} in the context of bootstrap). How do we
choose $h$?

We choose $h$ such that ${\bf Var} (\hat p)$ is minimal. For the
importance sampling estimator $\hat p$,
$${\bf Var}(\hat p) = {\bf Var}[g(X_i)f(X_i)/h(X_i)]/n$$


Now note that if $h(x)=cg(x)f(x)$ for some constant $c$, the
variance vanishes. Hence this choice of $h$ must be the minimizer.
Now the only requirement is that $h$ need to be a density. It is
satisfied if
$$h(x)=\frac{|g(x)|f(x)}{\int |g(y)|f(y) \; dy}.$$
See C.P. Robert and G. Casella, Monte Carlo Statistical Method for
Jensen's inequality based proof. But this is somewhat useless!
We do not know what the value of $\int |g(y)|f(y) \; dy$ is! You
can't apply importance sampling blindly. 
In next lecture we will
study how to choose $h$.

\item See Rubinstein, R.Y. (1981) Simulation and the Monte Carlo
Method for details on the importance sampling.
\end{enumerate}

\end{document}

The rule of thumb: choose $h$ such that $f/h$ is bounded and has a
thicker tail than $f$. Further $|g|f/h$ should be almost constant
with finite variance.

 \item We estimate $p$ based on the importance sampling
technique. we can show that $F_Y(y) =
\frac{1}{\pi}\arctan(\frac{y}{c}) + \frac{1}{2}$ and $F_Y^{-1}(y)
= -c\cot(\pi y)$. Based on this we can generate Cauchy
 Let $g$ to be the same index function as before and $f
\sim N(0,1)$.

 For instance, Let
$h(x) = \frac{c}{\pi}\frac{1}{c^2 + x^2}$. This is a general
Cauchy distribution with parameter $c$ so we may denote $h \sim
Cauchy(c)$, Let's see if this blind choice may give a good
estimate.

\begin{verbatim}
cauchypdf=inline('c./(pi*(c^2+x.^2))')
snormpdf=inline('exp(-x.^2/2)/sqrt(2*pi)') f=snormpdf(X);
h=cauchypdf(1,X); z=g.*f./h;

[mean(z) var(z)] ans =

    0.0123    0.0076

\end{verbatim}