bnp

Some older Bayesian nonparametrics research.
git clone git://git.jtobin.io/bnp.git
Log | Files | Refs | README | LICENSE

generalized_ibp.r (1333B)

      1 # two-parameter IBP, where b controls sparsity.  b = 1 yields the standard IBP.
      2 generalized_ibp = function(n, a, b) {
      3 
      4   dishes = max(1, rpois(1, a))
      5   diners = data.frame(dish = seq(dishes), diners = rep(1, dishes))
      6   buffet = list(buffet = diners, choices = list(seq(dishes)))
      7 
      8   for (j in seq(n - 1)) {
      9       buffet = arrival(buffet, a, b)
     10     }
     11   buffet
     12   }
     13 
     14 arrival = function(buf, a, b) {
     15   config = buf$buffet
     16   existing_choices = buf$choices
     17   j = length(existing_choices) + 1
     18   n = nrow(config)
     19 
     20   num_new_dishes = rpois(1, a * b / (j + b - 1))
     21   new_dishes     =
     22     if (num_new_dishes > 0) {
     23       (n + 1):(n + num_new_dishes)
     24     } else {
     25       numeric(0)
     26     }
     27 
     28   probs     = sapply(config$diners, function(n) { n / (b + j - 1) })
     29   selection = rbinom(n, 1, prob = probs) == 1
     30   choices   = list(c(as.integer(row.names(config[selection,])), new_dishes))
     31   diners    = with(config, replace(diners, selection, diners[selection] + 1))
     32 
     33   buffet = data.frame(
     34       dish   = seq(n + num_new_dishes)
     35     , diners = c(diners, rep(1, num_new_dishes))
     36     )
     37 
     38   list(buffet = buffet, choices = c(existing_choices, choices))
     39   }
     40 
     41 ibp_matrix = function(buffet) {
     42   dishes = do.call(max, buffet$choices)
     43   index  = function(j) { as.numeric(seq(dishes) %in% j) }
     44   rows   = lapply(buffet$choices, index)
     45 
     46   do.call(rbind, rows)
     47   }
     48