%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ------------------------------------------------------------------------
% Simulation Example *****************************************************
% ------------------------------------------------------------------------
% Submitted:      Dr. Eric Suess,Statistics 6501,Winter 2002
% Submitted Date: January 1st, 2002
% Accepted Date:  January 9th, 2002
% ------------------------------------------------------------------------
% This is a first draft of the handout "Particle Detector" in LaTeX.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\documentclass[12]{report}

\usepackage{latexsym}
\usepackage[centertags]{amsmath}
\usepackage{amsfonts}
\usepackage{newlfont}
%\usepackage{graphicx,fullpage,subfigure,pstricks}

\setlength{\oddsidemargin}{-0.4cm}
\setlength{\evensidemargin}{0cm}
\setlength{\topmargin}{-0.75cm}
\setlength{\textheight}{8.9in}
\setlength{\textwidth}{6.5in}
\renewcommand{\baselinestretch}{1}
\setlength{\parskip}{3mm}
\setlength{\parindent}{0mm} % indented paragraphs just don't
%\pagestyle{plain}

\begin{document}

\pagestyle{myheadings}
\markright{Statistics 6501, Winter 2009}

\large
\begin{center}
\textbf{Conditional Example}
\end{center}

\normalsize
\textbf{Particle Counter}

A particle counter is imperfect and independently detects each 
incoming particle with probability $p$.  If the distribution of the 
number of incoming particles in a unit of time is a Poisson 
distribution with parameter $\lambda$, what is the distribution of 
the number of counted particles?

Let $N = \#$ of incoming particles and $X = \#$ counted.

\begin{enumerate}
\item What is $P(X = k | N = n) =$ ?
\item What is $P(N=n)$?
\item Compute $P(X=k)$.
\end{enumerate}

\normalsize
\textbf{Is the conditional probability a regression?}

\begin{verbatim}
### imperfect particle counter

lam = 10
p = 0.9

# probability distribution

len = 5
z = matrix(0,1+len,1+len)

for(i in 0:len){
       for(j in 0:i){
               z[1+j,1+i] = dbinom(j,size=i,prob=p)
       }
}

z

z.sum = apply(z,2,sum)
z.sum

x.mean=0
for(i in 1:len){
       x.mean = c(x.mean,i*p)
}
x.mean

for(i in 0:len){
       X11()
       plot(z[,1+i],type="h")
}


# simulation, how to estimate p using linear regression through the
# origin

B = 1000000

n = rpois(B,lam)

x = rbinom(n=B,size=n,prob=p)

plot(n,x,main="Counts detected vs Counts emitted, with E[X|N]")

x.fit = lm(x ~ 0+n)
summary(x.fit)
abline(x.fit)

# Note that the estimated regression slope is very close to the true p

\end{verbatim}


\end{document}
% ------------------------------------------------------------------------
