bnp

commit e636ba1ac2f0c1b5bdbc1239ec88f0997a12fe4b
parent 353f1d88410b9ba30e83f83e24ad7fa14ef6cd20
Author: Jared Tobin <jared@jtobin.ca>
Date:   Sun, 21 Feb 2016 09:05:22 +1300

Add generalized IBP.

Diffstat:
A
indian-buffet-process/src/generalized_ibp.r
 | 
48
++++++++++++++++++++++++++++++++++++++++++++++++

1 file changed, 48 insertions(+), 0 deletions(-)
diff --git a/indian-buffet-process/src/generalized_ibp.r b/indian-buffet-process/src/generalized_ibp.r
@@ -0,0 +1,48 @@
+# two-parameter IBP, where b controls sparsity.  b = 1 yields the standard IBP.
+generalized_ibp = function(n, a, b) {
+
+  dishes = max(1, rpois(1, a))
+  diners = data.frame(dish = seq(dishes), diners = rep(1, dishes))
+  buffet = list(buffet = diners, choices = list(seq(dishes)))
+
+  for (j in seq(n - 1)) {
+      buffet = arrival(buffet, a, b)
+    }
+  buffet
+  }
+
+arrival = function(buf, a, b) {
+  config = buf$buffet
+  existing_choices = buf$choices
+  j = length(existing_choices) + 1
+  n = nrow(config)
+
+  num_new_dishes = rpois(1, a * b / (j + b - 1))
+  new_dishes     =
+    if (num_new_dishes > 0) {
+      (n + 1):(n + num_new_dishes)
+    } else {
+      numeric(0)
+    }
+
+  probs     = sapply(config$diners, function(n) { n / (b + j - 1) })
+  selection = rbinom(n, 1, prob = probs) == 1
+  choices   = list(c(as.integer(row.names(config[selection,])), new_dishes))
+  diners    = with(config, replace(diners, selection, diners[selection] + 1))
+
+  buffet = data.frame(
+      dish   = seq(n + num_new_dishes)
+    , diners = c(diners, rep(1, num_new_dishes))
+    )
+
+  list(buffet = buffet, choices = c(existing_choices, choices))
+  }
+
+ibp_matrix = function(buffet) {
+  dishes = do.call(max, buffet$choices)
+  index  = function(j) { as.numeric(seq(dishes) %in% j) }
+  rows   = lapply(buffet$choices, index)
+
+  do.call(rbind, rows)
+  }
+