First off, evaluating some existing code on the function from a website[http://sas-and-r.blogspot.com/2010/01/example-723-monty-hall-problem.html]
numsim = 1000 ; doors = 1:3
opendoor = function(x) { #what door will monty open?
if (x[1]==x[2]) #if winner==choice (contestant chose right) then...
return(sample(doors[-c(x[1])], 1)) #monty chooses a random wrong door.
else return(doors[-c(x[1],x[2])]) #if winner dne choice ..monty chooses the other bad door..so monty is helping you by showing you a bad door.
}
swapdoor = function(x) { return(doors[-c(x[1], x[2])]) } #drop the doors corresponding to x[1] and x[2]...those were cbind(open, choice)
winner = sample(doors, numsim, replace=TRUE) #the good door is "winner"
choice = sample(doors, numsim, replace=TRUE) #contestant chooses "choice"
open = apply(cbind(winner, choice), 1, opendoor) #monty then picks a door via opendoor()
newchoice = apply(cbind(open, choice), 1, swapdoor) #contestant swaps doors via swapdoor()
head(cbind(winner,choice,open,newchoice)) #these things collected.
## winner choice open newchoice
## [1,] 3 2 1 3
## [2,] 3 3 2 1
## [3,] 2 3 1 2
## [4,] 3 1 2 3
## [5,] 3 3 2 1
## [6,] 1 1 2 3
sum(winner==choice) / numsim #% correct via not swapping
## [1] 0.331
sum(winner==newchoice) / numsim #% correct via swapping
## [1] 0.669
Now to put these ideas into a functional wrapper that meets the criteria.
Monty = function(play.stay, play.swap) {
# conditions under noswap
winner = sample(doors, play.stay, replace = TRUE)
choice = sample(doors, play.stay, replace = TRUE)
open = apply(cbind(winner, choice), 1, opendoor)
newchoice = apply(cbind(open, choice), 1, swapdoor)
# conditions under swap
winner2 = sample(doors, play.swap, replace = TRUE)
choice2 = sample(doors, play.swap, replace = TRUE)
open2 = apply(cbind(winner2, choice2), 1, opendoor)
newchoice2 = apply(cbind(open2, choice2), 1, swapdoor)
noswap = sum(winner == choice)
Pct.noswap = noswap/play.stay
var.noswap = (Pct.noswap * (1 - Pct.noswap))/(play.stay - 1)
moe.noswap = 2 * sqrt(var.noswap)
CI.noswap = c(Pct.noswap - moe.noswap, Pct.noswap, Pct.noswap + moe.noswap)
swap = sum(winner2 == newchoice2)
Pct.swap = swap/play.swap
var.swap = (Pct.swap * (1 - Pct.swap))/(play.swap - 1)
moe.swap = 2 * sqrt(var.swap)
CI.swap = c(Pct.swap - moe.swap, Pct.swap, Pct.swap + moe.swap)
return(list(Number.noswap = noswap, CI.noswap.pct = CI.noswap, Number.swap = swap,
CI.swap = CI.swap))
}
Monty(100, 100)
## $Number.noswap
## [1] 38
##
## $CI.noswap.pct
## [1] 0.2824 0.3800 0.4776
##
## $Number.swap
## [1] 61
##
## $CI.swap
## [1] 0.512 0.610 0.708
##
The proof will follow but first the code is provided in Maria Rizzo's book. (Example 3.15)
Consider a Poisson random variable X with mean \( \lambda \). If \( \lambda \) is distributed as \( \Gamma(\alpha, \beta) \) then the marginal distribution on X is a negative binomial with parameters r = \( \alpha \) and p = \( 1/(\beta+1) \)
### Example 3.15 (Poisson-Gamma mixture)
# generate a Poisson-Gamma mixture
n <- 1000
r <- 4
beta <- 3
lambda <- rgamma(n, r, beta) #lambda is random
# now supply the sample of lambda's as the Poisson mean
x <- rpois(n, lambda) #the mixture
# compare with negative binomial
mix <- tabulate(x + 1)/n
negbin <- round(dnbinom(0:max(x), r, beta/(1 + beta)), 3)
se <- sqrt(negbin * (1 - negbin)/n)
round(rbind(mix, negbin, se), 3)
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## mix 0.308 0.316 0.219 0.101 0.036 0.013 0.004 0.001 0.002
## negbin 0.316 0.316 0.198 0.099 0.043 0.017 0.006 0.002 0.001
## se 0.015 0.015 0.013 0.009 0.006 0.004 0.002 0.001 0.001
There is also code illustrating the negative binomial as a mixture of Poissons with gamma distributed parameters in the Cassell book. (Example 2.6)
Nsim = 10^4
n = 6
p = 0.3
y = rgamma(Nsim, n, rate = p/(1 - p))
x = rpois(Nsim, y)
hist(x, main = "Negative Binomial, simulated and true", freq = F,
col = "grey", breaks = 40)
lines(1:50, dnbinom(1:50, n, p), lwd = 2, col = "sienna")
The negative binomial distribution:
\[ P(X=x) = \left(\frac{x+r-1}{r-1} \right){p}^{r}{q}^{r}, x=0,1,2,... \]
The Poisson distribution for a single observation (PDF):
\[ p(x) = \frac {{e}^{-\lambda}{\lambda}^{x}}{x!}, x=0,1,2,... \]
The Gamma distribution with scale parameter alpha and shape parameter beta
\[ g(\lambda) = \frac {{\alpha }^{\beta}} {\Gamma \left(\beta \right)}{\lambda }^{\beta -1}{e}^{-\alpha \lambda }, \lambda >0 \]
Now consider the Poisson distribution with lambda as a random variable:
\[ p(x|\lambda=\lambda) = \frac {{e}^{-\lambda}{\lambda}^{x}}{x!}, x=0,1,2,... \]
So a joint density looks like:
\[ p(x|\lambda=\lambda) g(\lambda) = \frac {{e}^{-\lambda}{\lambda}^{x}}{x!} \frac {{\alpha }^{\beta}} {\Gamma \left(\beta \right)}{\lambda }^{\beta -1}{e}^{-\alpha \lambda } \]
Pull out the constants, then add the lambdas & e's
\[ \int_{0}^{\infty} \frac{{\alpha }^{\beta}}{n!*\Gamma \left(\beta \right)}{\lambda }^{n+\beta -1}{e}^{-\left(\alpha +1 \right)\lambda }d\lambda \]
Multiply by one a couple times in order to get integral that integrates to one…
\[ \frac{{\alpha }^{\beta}}{n!*\Gamma \left(\beta \right)} \frac{\Gamma \left(n+\beta \right)}{\left({\alpha +1}\right)^{n+\beta }} \int_{0}^{\infty} \frac{{\left({\alpha +1}\right)^{n+\beta }}}{\Gamma \left(n+\beta \right)}{\lambda }^{n+\beta -1}{e}^{-\left(\alpha +1 \right)\lambda} d\lambda \]
So this one's easy, integrate a density and get one, so the factored constant remains
\[ \frac{{\alpha }^{\beta}}{n!*\Gamma \left(\beta \right)} \frac{\Gamma \left(n+\beta \right)}{\left({\alpha +1}\right)^{n+\beta }} \]
Here we are just rearranging and collecting terms with the same power.
\[ \frac{\Gamma \left(n+\beta \right)}{\Gamma \left(n+1 \right)\Gamma \left(\beta \right)}{\left(\frac{\alpha }{\alpha +1} \right)}^{\beta } {\left(\frac{1}{\alpha +1} \right)}^{n} \]
…and a final step. You need to know that \( \Gamma(n+ \beta) = (n+\beta-1)! \) (for integer values at least)
\[ \left(\begin{matrix}
n+\beta -1 \\
n
\end{matrix} \right)
{\left(\frac{\alpha }{\alpha +1} \right)}^{\beta } {\left(\frac{1}{\alpha +1} \right)}^{n}
, n =0,1,2,... \]
Following Robert and Casella (Chapter 2: random variable generation):
\[ P(Y\leq x|U\leq f(Y)/ \left(Mg(Y) \right))= \]
Just following the rule: P(A|B) = P(A&B) / P(B)
\[ =\frac
{P(Y\leq x, U\leq f(Y)/\left(Mg(Y) \right))}
{P(U\leq f(Y)/\left(Mg(Y) \right)} \]
Transform words into cdf's. Y is integrated up to x and U is integrated up to that fraction.
\[ \frac{\int_{-\infty}^{x}\int_{0}^{f(y)/Mg(y)}du*g(y)*dy}{\int_{-\infty}^{\infty}\int_{0}^{f(y)/Mg(y)}du*g(y)*dy} \]
Evaluating the first integral we get to:
\[ \frac{\int_{-\infty}^{x}f(y)/Mg(y)*du*g(y)*dy}{\int_{-\infty}^{\infty}f(y)/Mg(y)*du*g(y)*dy} \]
Now everything cancels. You are dividing by and multiplying g(x). After that's cleared, there's an M in the numerator and the denominator that clear, leaving…
\[ \frac{\int_{-\infty}^{x}f(y)dy}{\int_{-\infty}^{\infty}f(y)dy} \]
Integrating a pdf over its entire support yields one in the denominator so we have the cdf of x left exactly.
\[ \int_{-\infty}^{x}f(y)dy=P(X\leq x) \]
par(mfrow = c(1, 1))
nsim = 10000
# the multimodal pdf h
h = function(x) {
x * sin(20 * x) + 1
}
constant = max(h(x))
x = seq(0, 1, by = 0.01)
plot(x, h(x), type = "l", ylim = c(0, max(h(x)))) #
# an ad hoc way to look at the accept-reject method (with poor acceptance
# rate) uses a uniform for each axis. 'throwing darts'
YY = runif(nsim, 0, max(h(x)))
XX = runif(nsim, 0, 1)
factoR = YY < h(XX)
col.points = c("grey", "blue")[as.factor(factoR)]
points(XX, YY, col = col.points, pch = 20)
frac.acc = mean(factoR)
arEa = frac.acc * max(h(x)) * 1
## the accept reject algorithm for generating from our arbitrary pdf h.
## the enveloping function f1
f1 = function(x) {
constant * dunif(x)
}
AlgoU = runif(nsim)
AlgoU2 = runif(nsim)
Ratio = h(AlgoU)/constant * f1(AlgoU)
accept.rate = sum(AlgoU2 < Ratio)/nsim
accept.rate
## [1] 0.833
rfunc = AlgoU[AlgoU2 < Ratio]
par(new = T)
hist(rfunc, freq = F, ylim = c(0, max(h(x))))
lines(x, f1(x))
## perhaps I'll do a 3d later for fun... doesnt say 2d or 3d in
## question... robert and casella page 135
h = function(x, y) {
(x * sin(20 * y) + y * sin(20 * x))^2 * cosh(sin(10 * x) * x) + (x * cos(10 *
y) - y * sin(10 * x))^2 * cosh(cos(20 * y) * y)
}
x = y = seq(-3, 3, le = 435) #defines a grid for persp
z = outer(x, y, h)
par(bg = "wheat", mar = c(1, 1, 1, 1)) #bg stands for background
persp(x, y, z, theta = 155, phi = 30, col = "green4", ltheta = -120,
shade = 0.75, border = NA, box = T)
nsim = 10
x = y = runif(nsim, min = -3, max = 3)