discrete.R

nesau — Oct 30, 2014, 2:13 PM

n = 30
m = 100

# X ~ Poisson(lambda=5)
fx <- function(x) dpois(x,lambda=5)
# N ~ Geom(p=0.6)
fn <- function(x) dgeom(x,p=0.5)

f <- vector("list", (n+1))
f[[1]] = c(1, rep(0,m)) # f*0
f[[2]] = unlist(lapply(0:m, fx)) # f*1
# apply a recursive algorithm to generate f[[3]] ... f[[n+1]]
for(l in 3:(n+1)) {
  f[[l]] = numeric(m+1)
  for(x in 0:m) {
    total = 0
      for(y in 0:x) {
        total = total + fx(y)*f[[l-1]][x-y+1]
      }
    f[[l]][x+1] = total
  }
}

# S = X1 + X2 + .. + XN
fs <- numeric(n+1) 
for(x in 0:m) {
  total = 0
  for(i in 1:n) {
    total = total + f[[i+1]][x+1]*fn(i)
  }
  fs[x+1] = total
}
Fs <- unlist(lapply(0:m, function(i) { force(i); sum(fs[1:(i+1)])}))

x = 0:m
par(mfrow=c(2,2))
plot(x,fs,type='l')
plot(x,Fs,type='l')
# compare to normal approximation
Es = ( (1 - 0.5) / 0.5) * 5
Vs = 5^2 * ((1 - 0.5) / (0.5)^2) + ((1 - 0.5)/(0.5)) * 5
plot(function(x) pnorm((x - Es)/sqrt(Vs)), 0, 100, ylab="Fz")