# R program for calculations in the beer game
# Jens Lund, March 2003

s <- 6 # Number of sides on each dice
k <- 6 # Number of dices

fac <- function(x) gamma(x+1)

# Density for the j'th highest dice among k independent dices with s sides.
ddj <- function(i,j,k,s){
  # i: point in which to take the density, i = 1,...,s
  # j: rank of dice, j=1,...,k
  # k: number of dices
  # s: sides on each dice
  polydensity <- function(yplus,yi,i,k,s){
    fac(k)/fac(yplus)/fac(yi)/fac(k-yplus-yi)*((s-i)/s)^yplus/(s^yi)*((i-1)/s)^(k-yplus-yi)
  }
  innersum <- function(yplus) {
    tmpfun <- function(yi) polydensity(yplus,yi,i,k,s)
    sum(unlist(lapply((j-yplus):(k-yplus),tmpfun)))
  }
  sum(unlist(lapply(0:(j-1),innersum))) 
}

# Mean for the j'th highest dice among k independent dices with s sides
dmeanj <- function(j,k,s){
  ddjks <- function(i) ddj(i,j,k,s)
  sum((1:s)*unlist(lapply(1:s,ddjks)))
}

# Print information matrix
# k x k matrices
# Rows count j, the rank of the dice (counting 1 from the top)
# Columns count number of remaining dices (1 at the left, k at the right)
# Matrix 1 gives the mean of the j'th highest roll among k left dices
# Matrix 2 gives the cumulative mean of the 1 to j'th highest rolls.
#  I.e. sums over rows
# Matrix 3 gives the sum of the 1 highest rolls in the remaing k rolls.
#  I.e. sums over columns to the left. Each row sums this many columns to
#  the left.
#  I.e. entry (2,4) gives the mean sum of the highest rolls
#  with 4 and 3 dices.

printinfomatrix <- function(kstart,s){
  # Initialize
  mat1 <- matrix(0,nrow=k,ncol=k)
  mat2 <- mat1
  mat3 <- mat1

  # Do means (not very efficient due to the for loops)
  for (k in 1:kstart)
    for (j in 1:k){
      mat1[j,k] <- dmeanj(j,k,s)
      if(j==1) mat2[j,k] <- mat1[j,k]
      else mat2[j,k] <- mat2[j-1,k]+mat1[j,k]
      mat3[j,k] <- sum(mat1[1,(k-j+1):k])
    }
  list(mean=mat1,jhighest=mat2,rollshighest=mat3)
}

print(printinfomatrix(k,s))