discretized.R

#### LM 3rd, Example 9.5 ####
fx = c(0.150, 0.200, 0.250, 0.125, 0.075, 0.050, 
       0.050, 0.050, 0.025, 0.025)
fn = c(0.05,0.10,0.15,0.20,0.25,0.15,0.06,0.03,0.01)

# Generate the distribution of S
f = vector("list", 8)
f[[1]] = c(1,rep(0,80))
f[[2]] = c(0, fx, rep(0,70))
for(i in 2:8) {
  f[[i+1]] = numeric(81) # init
  for(x in i:(10*i)) {
    total = 0
    for(y in 1:min(x,10)) {
      total = total + f[[i]][x-y+1]*fx[y]
    }
    f[[i+1]][x+1] = total
  }
}

fs = numeric(81)
for(x in 0:80) {
  total = 0
  for(n in 0:8) {
    total = total + f[[n+1]][x+1] * fn[n+1]
  }
  fs[[x+1]] = total
}

# Display the data
options("scipen"=100)
distribution = data.frame(0:80, round(f[[1]],5), round(f[[2]],5), round(f[[3]],5),
                          round(f[[4]],5), round(f[[5]],5), round(f[[6]],5), 
                          round(f[[7]],5), round(f[[8]],5), round(f[[9]],5), 
                          round(fs,5))
names(distribution) = c('x', 'f*(0)', 'f*(1)', 'f*(2)', 
                        'f*(3)', 'f*(4)', 'f*(5)', 'f*(6)', 
                        'f*(7)', 'f*(8)', 'fs')

#### Method of Rounding ####
g <- function(x) dexp(x, 1/10)
G <- function(x) pexp(x, 1/10)

m = 200
h = 2 # span
f = numeric(m/h+1)
f[1] = G(h/2 - 0) # Pr(X < h/2)
for(j in 1:(m/h)) { # Pr(jh - h/2 <= X < jh + h/2)
  f[j+1] = G(j*h + h/2 - 0) - G(j*h - h/2 - 0)
}
x = seq(0, m, by=h)
cat(sum(f*x), "\n", integrate(function(x) x*g(x), 0, 200)$value, "\n", sep="")

## 9.983
## 10

#### Method of Local Moment Matching ####
g <- function(x) dexp(x, 1/10)

p = 1
klimit = 200
h = 2 # span
m = vector("list", p+1)
for(i in 1:(p+1)) m[[i]] = numeric((klimit)/h)
for(k in 0:(klimit/h)) {
  for(j in 0:p) {
    f = function(x) {
      total = 1
      for(i in 0:p) {
        if(i!=j) {
          total = total * (x - h*k - i*h)/((j - i)*h)
        }
      }
      total = total * g(x)      
    }
    m[[j+1]][k+1] = integrate(f, h*k, h*k + p*h)$value
  }
}

f = numeric( (klimit)/h+1)
for(j in 0:(p-1)) { # row 0
  f[j+1] = m[[j+1]][1] # p entries per row
}

for(i in 1:(klimit/(h*p) - 1)) { # row 1 - row klimit/(h*p)
   for(j in 0:(p-1)) { # fp+j
     if(j==0) { # f(ip + 0)
       f[i*p+j+1] = m[[p+1]][i] + m[[1]][i+1]
     } else {
       f[i*p+j+1] = m[[j+1]][i+1]
     }
   }
}
x = seq(0, klimit, by=h)
cat(sum(f*x), "\n")

## 10