bnp

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

dpmm.r (971B)

      1 BNP_DIR = "/Users/jtobin/projects/bnp"
      2 SBP_SRC = paste(BNP_DIR, "stick-breaking-process/src/sbp.r", sep = "/")
      3 
      4 require(gtools)
      5 require(magrittr)
      6 source(SBP_SRC)
      7 
      8 mixing_model = function(n, a) {
      9   if (n <= 1) {
     10     stop("need > 1 observation.")
     11     } else {
     12     sbp(n - 1, a)
     13     }
     14   }
     15 
     16 label_model     = function(n, p) drop(rmultinom(1, size = n, prob = p))
     17 location_model  = function(k, l, r) rnorm(k, l, 1 / r)
     18 precision_model = function(k, b, w) rgamma(k, b, 1 / w)
     19 
     20 parameter_model = function(n, a) {
     21   p  = mixing_model(n, a)
     22   k  = length(p)
     23   c  = label_model(n, p)
     24   mu = location_model(k, 0, 0.1)
     25   s  = precision_model(k, 1, 1)
     26   list(c, mu, s)
     27   }
     28 
     29 data_model = function(config) {
     30   sampler = function(y, m, s) rnorm(y, m, 1 / s)
     31   mapply(sampler, config[[1]], config[[2]], config[[3]])
     32   }
     33 
     34 model = function(n, a) {
     35   clusters = parameter_model(n, a) %>% data_model
     36   nonempty = unlist(lapply(clusters, function(x) length(x) != 0))
     37   clusters[nonempty]
     38   }
     39