bnp

commit 714f125cb6e56d6d42cbaed512a27fa774003bb4
parent ca36a3f2f866e4876afe3dbe029776a2edb9c6f2
Author: Jared Tobin <jared@jtobin.ca>
Date:   Sat, 20 Feb 2016 20:11:13 +1300

Merge pull request #1 from jtobin/jt-ibp

Add generative Indian Buffet Process.
Diffstat:
A
indian-buffet-process/src/ibp.r
 | 
46
++++++++++++++++++++++++++++++++++++++++++++++

1 file changed, 46 insertions(+), 0 deletions(-)
diff --git a/indian-buffet-process/src/ibp.r b/indian-buffet-process/src/ibp.r
@@ -0,0 +1,46 @@
+
+ibp = function(n, a) {
+  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)
+    }
+  buffet
+  }
+
+arrival = function(b, a) {
+  config  = b$buffet
+  existing_choices = b$choices
+  j = length(existing_choices) + 1
+  n = nrow(config)
+
+  num_new_dishes = rpois(1, a / j)
+  new_dishes     =
+      if (num_new_dishes > 0) {
+        (n + 1):(n + num_new_dishes)
+      } else {
+        numeric(0)
+      }
+
+  probs     = sapply(config$diners, function(n) { n / j })
+  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 = Reduce(max, buffet$choices)
+  index  = function(j) { as.numeric(seq(dishes) %in% j) }
+  rows   = lapply(buffet$choices, index)
+
+  do.call(rbind, rows)
+  }
+