buhlmann.R

Nathan Esau — Mar 2, 2014, 10:44 PM

probs <- c(0.1,0.2,0.3,0.4)
bayesians <-c(1,1.5,2.5,4)
data <- c(c(2,2), c(2,3), c(3,2), c(3,3)); data <- matrix(data,nrow=4,ncol=2, byrow=TRUE)
pastObs <- c(3,3)

rref <- function(A, tol=sqrt(.Machine$double.eps),verbose=FALSE,
                 fractions=FALSE){
  ## A: coefficient matrix
  ## tol: tolerance for checking for 0 pivot
  ## verbose: if TRUE, print intermediate steps
  ## fractions: try to express nonintegers as rational numbers
  ## Written by John Fox
  if (fractions) {
    mass <- require(MASS)
    if (!mass) stop("fractions=TRUE needs MASS package")
  }
  if ((!is.matrix(A)) || (!is.numeric(A)))
    stop("argument must be a numeric matrix")
  n <- nrow(A)
  m <- ncol(A)
  for (i in 1:min(c(m, n))){
    col <- A[,i]
    col[1:n < i] <- 0
    # find maximum pivot in current column at or below current row
    which <- which.max(abs(col))
    pivot <- A[which, i]
    if (abs(pivot) <= tol) next     # check for 0 pivot
    if (which > i) A[c(i, which),] <- A[c(which, i),]  # exchange rows
    A[i,] <- A[i,]/pivot            # pivot
    row <- A[i,]
    A <- A - outer(A[,i], row)      # sweep
    A[i,] <- row                    # restore current row
    if (verbose)
      if (fractions) print(fractions(A))
    else print(round(A,round(abs(log(tol,10)))))
  }
  for (i in 1:n)
    if (max(abs(A[i,1:m])) <= tol)
      A[c(i,n),] <- A[c(n,i),] # 0 rows to bottom
  if (fractions) fractions (A)
  else round(A, round(abs(log(tol,10))))
}

buhlmann <- function(probs, bayesians, x, pastObs) {
  numrowsx = length(x[, 1]) # first column of x
  numcolsx = length(x[1 ,]) # first row of x

  Q <- numeric(0)
  # go through column-wise and get alphas coefficients
  for (k in (1:length(probs))) {
    Q = c (Q, 1)  # a0 coefficient

    # row-wise
    for (j in (1:length(x[1,]))) {
      currentRow = x[k,]
      Q = c(Q,currentRow[j])
    }
  }

  # X1col represents coefficients of ao, X2col represnets coefficients of a1...
  Q = matrix(Q,nrow=length(probs), ncol=numcolsx+1, byrow=TRUE)

  numrowsQ = length(Q[,1]) # first column of Q
  numcolsQ = length(Q[1,]) # first row of Q
  normal <- numeric(0)      # eventual matrix for normal eqns  

  for (m in (1:numcolsQ)) {
    eqns <- numeric(numcolsQ+1) # store one row of normal eqns at a time
    # gets the nth normal eqn
    for (k in (1:numrowsQ)) {
      # factor for a0 normal eqns
      factor = Q[,m][k]
      currentRow = Q[k,]
      for (j in (1:(numcolsQ+1))) {
        if(j>numcolsQ) {
          eqns[j] = eqns[j] + bayesians[k]*probs[k]*factor
        } else {
          eqns[j] = eqns[j] + currentRow[j]*probs[k]*factor
        }
      }
    }
    normal = c(normal,eqns)
  }
  # dimensions of normaleqn matrix
  #   rows: numcolsQ  
  #   cols: numcolsQ+1 (for constants)
  normal = matrix(normal,nrow=numcolsQ, ncol=numcolsQ+1, byrow=TRUE)

  # row reduce normal eqns
  normal = rref(normal)
  alphas = normal[,numcolsQ+1]

  total = 0
  for (k in (0:length(pastObs))) {
    if (k==0) {
      total = total + alphas[k+1]
    }
    else {
      total = total + alphas[k+1]*pastObs[k]
    }
  }
  return(total)
}

ans = buhlmann(probs,bayesians,data, pastObs)
print(ans)

[1] 3.88