bnp

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

fmm_generative.r (1036B)


      1 source('fmm_utils.r')
      2 
      3 mixing_model    = function(k, a) drop(rdirichlet(1, (rep(a, k))))
      4 label_model     = function(n, p) drop(rmultinom(1, size = n, prob = p))
      5 location_model  = function(k, l, r) rnorm(k, l, sqrt(1 / r))
      6 precision_model = function(k, b, w) rgamma(k, b, 1 / w)
      7 
      8 parameter_model = function(k, n) {
      9   p  = mixing_model(k, 1)
     10   c  = label_model(n, p)
     11   mu = location_model(k, 0, 0.1)
     12   s  = precision_model(k, 1, 1)
     13   list(c, mu, s)
     14 }
     15 
     16 data_model = function(config) {
     17   sampler = function(y, m, s) rnorm(y, m, 1 / s)
     18   mapply(sampler, config[[1]], config[[2]], config[[3]])
     19 }
     20 
     21 model = function(k, n) {
     22   params = parameter_model(k, n)
     23   data_model(params)
     24 }
     25 
     26 lmodel = function(y, p, m, s) {
     27   score      = function(pr, mu, prec) { pr * dnorm(y, mu, sqrt(1 / prec)) }
     28   by_cluster = mapply(score, p, m, s)
     29   totalled   = apply(by_cluster, MARGIN = 1, sum)
     30 
     31   # NOTE (jtobin): adjusted for numerical stability
     32   small    = 1.379783e-316
     33   adjusted = totalled
     34   adjusted[which(adjusted == 0)] = small
     35   sum(log(adjusted))
     36 }
     37