- Task 2 Simulating Negative Binomial Distributions In This Task We Will Compare Two Different Algorithms For Simulat 1 (171.27 KiB) Viewed 103 times
In this task, we will compare two different algorithms for
simulating from a negative binomial distribution.
### Problem (a)
Recall that a negative binomial random variable *NB(r, p)* is the
sum of *r Geometric(p)* random variables. Use the algorithm
from Task 1 to simulate 1000 *NB(10, 0.6)* random variates.
### Code set-up
Note that we merely need to wrap the core code from Task 1 within a
for-loop. Here is the core of the code chunk, where we are
thinking of a for-loop over a variable `sims` to replicate the
single negative binomial draw. Note that this code chunk will
not run since the for-loop over `sims` is not coded, thus the
`eval=FALSE` option. **Note that this code has the
`eval=FALSE` option just to present the code without output.
Your code will not use this option.**
```{r codecore, eval=FALSE}
r=10
nbvar1 = numeric(r)
for(nbsims in 1:r){
# for-loop allows us to simulate until r successes;
# in this problem, r=10 and p=0.6
tossnum = 0
success = 0
while(success==0){
success = (runif(1)<p)
tossnum = tossnum + 1
}
nbvar1[sims] = nbvar1[sims] + tossnum
}
```
### Problem (b)
The negative binomial pmf induces the following recursion relation.
If $X \sim B(r,\ p)$, then
\[ P(X = i+1) = \frac{i(1-p)}{(i+1-r)} \cdot P(X=i).\]
Use this recursion relation to generate 1000 $NB(10,\ 0.6)$ random
variates.
### Code set-up
Below is binomial.R, the binomial simulator used in the video
lectures and found also on the class Blackboard site.
```{r binomial, echo=TRUE}
simnum = 1000
p = 0.6; r = 10 # for point of comparison with the negative
binomial, we will use r here
y=0
for(sims in 1:simnum){
pmf=(1-p)^r; cdf=pmf; # pmf and cdf
j=0;
u=runif(1) # uniform random variate
# find Binomial variate
while(u >= cdf){
pmf=((r-j)/(j+1))*(p/(1-p))*pmf # recursion relation
cdf=cdf + pmf # compute
cdf
j=j+1
}
y[sims]=j
}
```
This binomial simulator may be applied directly after changing
just three lines:
* $j=r$
* the recursion relation formula
* `pmf = p^r`
### Report the following for each of the simulations in problems
(a) and (b)
* Histogram of the variates
* Mean and standard deviation of the simulated variates
* Run time: compare computing speed between the two algorithms.
In R, can wrap your algorithm or sequence of operations as
follows to time your code.
```{r timewrapper}
x = proc.time()
# [the code you want to time here]
timer = proc.time() - x
algtime = timer[3] # algtime will store the algorithm run
time in seconds
```
### Questions
* How do the histograms compare?
*[Answer here]*
* How do the mean and standard deviation from the simulations
compare to the true mean and standard deviation of a $NB(0.6,\ 10)$
distribution?
*[Answer here]*
* How do the computing times compare? Which algorithm is
faster?
*[Answer here]*
* "Simulation flops": Which simulator do you think uses more
uniform random numbers (call to the `runif()` function)?
Why?
*[Answer here]*
## Task 2: Simulating negative binomial distributions In this task, we will compare two different algorithms for simulating from a negative binomial distribution. ### Problem (a) Recall that a negative binomial random variable *NB(r, p)* is the sum of *r Geometric(p)* random variables. to simulate 1000 *NB(10, 0.6)* random variates. Use the algorithm from Task 1 ### Code set-up Note that we merely need to wrap the core code from Task 1 within a for-loop. Here is the core of the code chunk, where we are thinking of a for-loop over a variable ‘sims' to replicate the single negative binomial draw. Note that this code chunk will not run since the for-loop over 'sims is not coded, thus the eval=FALSE' option. **Note that this code has the 'eval=FALSE' option just to present the code without output. Your code will not use this option. ** {r codecore, eval=FALSE} r=10 nbvar1 = numeric(r) for(nbsims in 1:r){ # for-loop allows us to simulate until r successes; # in this problem, r=10 and p=0.6 tossnum 0 success while(success==0){ success = (runif(1)<p) tossnum tossnum + 1 } nbvar1[sims] nbvar1[sims] + tossnum }
### Problem (b) The negative binomial pmf induces the following recursion relation. If $X \sim Bar, \ p)$, then \[ PCX = i +1) = \frac{i(1-p)}{(i+1-r)} \cdot PCX=i). ] Use this recursion relation to generate 1000 $NB(10, 0.6)$ random variates. ### Code set-up Below is binomial.R, the binomial simulator used in the video lectures and found also on the class Blackboard site. **'{r binomial, echo=TRUE} simnum = 1000 p = 0.6; r = 10 # for point of comparison with the negative binomial, we will use r here y=0 for(sims in 1: simnum) { pmf=(1- pr; cdf pmf; # pmf and cdf jo; u=runif(1) # uniform random variate # find Binomial variate while(u >= cdf){ pmf=((r-j)/(1+1))*(p/(1-p)*pmf # recursion relation cdf-cdf + pmf # compute cdf j-j+1 } y[sims] =j This binomial simulator may be applied directly after changing just three lines: * Sj=rs * the recursion relation formula * pmf = par
### Report the following for each of the simulations in problems (a) and (b) * Histogram of the variates * Mean and standard deviation of the simulated variates * Run time: compare computing speed between the two algorithms. In R, can wrap your algorithm or sequence of operations as follows to time your code. {r timewrapper) X = proc. time # [the code you want to time here] timer = proc. time algtime = timer [3] # algtime will store the algorithm run time in seconds х ### Questions * How do the histograms compare? *[Answer here]* * How do the mean and standard deviation from the simulations compare to the true mean and standard deviation of a $NB(0.6, \ 10$ distribution? * [Answer here]* * How do the computing times compare? Which algorithm is faster? *[Answer here]* * "Simulation flops": Which simulator do you think uses more uniform random numbers (call to the 'runif() function)? Why? *[Answer here]