---
title: "Conditional Example"
author: "Prof. Eric A. Suess"
date: "August 31, 2022"
format: pdf
---

## 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.

1. What is $P(X = k | N = n) =$ ?
2. What is $P(N=n)$?
3. Compute $P(X=k)$.

**Is the conditional probability a regression?**

## Simulation of the conditional distribution of $X|N = n$.

The imperfect particle counter

```{r}
lam <- 10
p <- 0.9
```

The conditional probability distributions for $n = 0, 1, 2,...,5$.

```{r}
len <- 5
x <- matrix(0, 1+len, 1+len) # distributions down the columns of the matrix

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

n <- c(0,1,2,3,4,5)

dimnames(x) <- list(n, n)

x
```

The sums of each conditional distribution should be 1.

```{r}
x.sum <- apply(x, 2, sum) # column sum
x.sum
```

Determine the conditional mean $E[X|N = n]$  Column means.

```{r}
x.mean <- 0

for(i in 1:len){
       x.mean <- c(x.mean, i*p)
}
x.mean
```

Plots of each probability distribution.

```{r}
for(i in 0:len){
       plot(n, x[,1+i], type="h", main = paste("n = ", i), ylab = "prob")
}
```

```{r}
library(scatterplot3d)

x

s3d.dat <- data.frame(columns = c(col(x)), rows = c(row(x)), prob = c(x))

scatterplot3d(s3d.dat, type = "h", lwd = 5, pch = " ",
              x.ticklabs = colnames(x), y.ticklabs = rownames(x),
              main = "3D barplot")
```


## Simulation, how to estimate $p$ using linear regression through the origin.

Simulate the data.

```{r}
B <- 20000

n <- rpois(B,lam)

mean(n)

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

mean(x)
```

Fit the model through the origin.

```{r}
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, col="blue", lwd = 3)        # E[X|N]
abline(v = mean(n), col="navy", lwd = 3)  # E[N] 
abline(h = mean(x), col="navy", lwd = 3)  # E[X]
```

```{r}
library(car)

scatterplot(n, x, xlab="n", ylab="x", main="E[X|N]", smooth=FALSE)
```

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