Baysian simulation for the German tank problem

tank.simul.R

### Baysian estimate simulation of the "German tank problem"
### Coded by everdark, since 2013-01-30



## Simple simulation
set.seed(12345)
j <- 1000
range <- c(14:50)
result <- list()
temp <- numeric(j)
for (k in range) {
    for (i in 1:j) {
        temp[i] <- max(sample(1:k, 4))
    }
    result[[k]] <- temp
}

result <- 
    unlist(
        lapply(result, function(x) prop.table(table(x))["14"])
    )
result[is.na(result)] <- 0
plot(result)
sum(result*seq(1, length(result), by=1)) / sum(result)



## Single node computing
set.seed(12345)
j <- 1000
paraMax <- 100
sampMax <- 14
sampsize <- 4
range <- c(sampMax:paraMax)
temp <- numeric(j)
result.1 <- list()
outcome <- numeric(max(range))
system.time(
    for (k in range) {
        for (i in 1:j) temp[i] <- max(sample(1:k, sampsize))
        result.1[[k]] <- temp
        result.2 <- 
            unlist(
                lapply(result.1, function(x) prop.table(table(x))[as.character(sampMax)])
            )
        result.2[is.na(result.2)] <- 0
        outcome[k] <- sum(result.2*seq(1, length(result.2), by=1)) / sum(result.2)
    }
)
plot(outcome)



## Parallel computing (4-way)
library(snow)     # "Simple Network Of Workstation"
library(rlecuyer) # parallel random number generator

tank.simul.2 <- function(kchunks, RANGE, times, nsample) {
    
    j <- times
    paraMax <- max(RANGE)
    sampMax <- min(RANGE)
    sampsize <- nsample
    result.1 <- data.frame(shell=rep(-1, j))
    temp <- numeric(j)
    for (k in kchunks) {
        for (i in 1:j) {
            temp[i] <- max(sample(1:k, sampsize))
        }
        result.1 <- cbind(result.1, temp)
        colnames(result.1)[ncol(result.1)] <- as.character(k)
    }
    subset(result.1, select=-c(shell))
    
}

cluster <- makeCluster(type="SOCK", 
                       c("localhost", "localhost", "localhost", "localhost"))
clusterSetupRNGstream(cluster, seed=rep(12345,6))
range <- c(14:500)
kchunks <- clusterSplit(cluster, range)

system.time(
    result.1 <- clusterApply(cluster, kchunks, 
                             tank.simul.2, RANGE=range, times=5000, nsample=4)
)
stopCluster(cluster)

result.1 <- as.data.frame(result.1)
result.2 <- 
    apply(
        result.1, 2, 
        function(x) prop.table(table(x))[as.character(min(range))]
    )
result.2[is.na(result.2)] <- 0

outcome <- numeric(length(range))
for (end in min(range):max(range)) {
    i <- end - min(range) + 1
    outcome[i] <- sum(result.2[1:i]*seq(min(range), end, by=1)) / sum(result.2[1:i])
}
plot(outcome)
abline(h=19.5, col="red")
outcome[length(outcome)]
final <- outcome[length(outcome)]
min(which(outcome==final))



## Asymptotic simulation
cluster <- makeCluster(type="SOCK", 
                       c("localhost", "localhost", "localhost", "localhost"))
clusterSetupRNGstream(cluster, seed=rep(12345,6))
range <- c(14:1000)
kchunks <- clusterSplit(cluster, range)
iterator <- seq(1000, 10000, by=250)
outcome <- numeric(length(iterator))
i <- 0
for (t in iterator) {
    
    i <- i + 1
    
    result.1 <- 
        as.data.frame(
            clusterApply(
                cluster, kchunks, tank.simul.2, RANGE=range, times=t, nsample=4
            )
        )
    result.2 <- 
        apply(
            result.1, 2, 
            function(x) prop.table(table(x))[as.character(min(range))]
        )
    result.2[is.na(result.2)] <- 0
    outcome[i] <- sum(result.2*seq(min(range), max(range), by=1)) / sum(result.2)
}
stopCluster(cluster)

plot(iterator, outcome)
Bestimate <- (min(range)-1) * (3) / (2)
abline(h=Bestimate, col="red")
mean(outcome)