bnp

commit 881bc344619017da98d1b4222e5ab363b000d90a
parent 22c63244aa68fcf0cc8d577dc248ca9a87d82aa5
Author: Jared Tobin <jared@jtobin.ca>
Date:   Sat, 20 Feb 2016 20:10:28 +1300

Add IBP.

Diffstat:
M
indian-buffet-process/src/ibp.r
 | 
52
++++++++++++++++++++++++++++++----------------------

1 file changed, 30 insertions(+), 22 deletions(-)
diff --git a/indian-buffet-process/src/ibp.r b/indian-buffet-process/src/ibp.r
@@ -1,10 +1,8 @@
-# FIXME handle sampled zero values
 
 ibp = function(n, a) {
-  dishes = rpois(1, 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)
     }
@@ -12,27 +10,37 @@ ibp = function(n, a) {
   }
 
 arrival = function(b, a) {
-
   config  = b$buffet
-  choices = b$choices
-
-  j          = length(choices) + 1
-  n          = nrow(config)
-  probs      = sapply(config$diners, function(n) { n / j })
-  selection  = rbinom(n, 1, prob = probs) == 1
+  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))
+    )
 
-  existing_diners = config[selection, 'diners']
-  new_diners = with(config, replace(diners, selection, diners[selection] + 1))
+  list(buffet = buffet, choices = c(existing_choices, choices))
+  }
 
-  existing_dishes = data.frame(dish = config$dish, diners = new_diners)
-  num_new_dishes  = rpois(1, a / j )
-  new_dishes = data.frame(
-      dish   = (n + 1):(n + num_new_dishes)
-    , diners = rep(1, num_new_dishes)
-    )
+ibp_matrix = function(buffet) {
+  dishes = Reduce(max, buffet$choices)
+  index  = function(j) { as.numeric(seq(dishes) %in% j) }
+  rows   = lapply(buffet$choices, index)
 
-  list(
-      buffet  = rbind(existing_dishes, new_dishes)
-    , choices = list(choices, rep(TRUE, num_new_dishes)) # FIXME
-    )
+  do.call(rbind, rows)
   }
+