<commit_before># Monadic probabilistic programming in Scala with Rainier

```scala
set libraryDependencies += "com.stripe" %% "rainier-core" % "0.1.0"
set libraryDependencies += "com.cibo" %% "evilplot" % "0.2.0"
set resolvers += Resolver.bintrayRepo("cibotech", "public")
set scalaVersion := "2.12.4"
console
```

```scala
import com.stripe.rainier.compute._
import com.stripe.rainier.core._
import com.stripe.rainier.sampler._
import com.stripe.rainier.repl._

case class Bernoulli(p: Real) extends Distribution[Int] {

  def logDensity(b: Int): Real = {
    p.log * b + (Real.one - p).log * (1 - b)
  }

  val generator = Generator.from { (r, n) =>
    val pd = n.toDouble(p)
    val u = r.standardUniform
    if (u < pd) 1 else 0
  }

}

// first simulate some data from a logistic regression model
val r = new scala.util.Random(0)
val N = 1000
val beta0 = 0.1
val beta1 = 0.3
val x = (1 to N) map { i =>
  3.0 * r.nextGaussian
}
val theta = x map { xi =>
  beta0 + beta1 * xi
}
def expit(x: Double): Double = 1.0 / (1.0 + math.exp(-x))
val p = theta map expit
val y = p map (pi => if (r.nextDouble < pi) 1 else 0)
println(y.take(10))
println(x.take(10))

// now build and fit model
val model = for {
  beta0 <- Normal(0, 5).param
  beta1 <- Normal(0, 5).param
  _ <- Predictor.from{x: Double =>
      {
        val theta = beta0 + beta1 * x
        val p = Real(1.0) / (Real(1.0) + (Real(0.0) - theta).exp)
        //Bernoulli(p)
        Binomial(p,1)
      }
    }.fit(x zip y)
} yield (beta0, beta1)

//implicit val rng = RNG.default
implicit val rng = ScalaRNG(3)
val its = 10000
val thin = 5
//val out = model.sample(Walkers(100), 10000, its)
val out = model.sample(HMC(5), 10000, its*thin, thin)
//val out = model.toStream(HMC(5), 1000).take(its).map(_()).toArray
println(out.take(10))
```

```scala
import com.cibo.evilplot.plot._
import com.cibo.evilplot.colors._
import com.cibo.evilplot.plot.aesthetics.DefaultTheme._<commit_msg>bernoulli rather than binomial...<commit_after># Monadic probabilistic programming in Scala with Rainier

```scala
set libraryDependencies += "com.stripe" %% "rainier-core" % "0.1.0"
set libraryDependencies += "com.cibo" %% "evilplot" % "0.2.0"
set resolvers += Resolver.bintrayRepo("cibotech", "public")
set scalaVersion := "2.12.4"
console
```

```scala
import com.stripe.rainier.compute._
import com.stripe.rainier.core._
import com.stripe.rainier.sampler._
import com.stripe.rainier.repl._

case class Bernoulli(p: Real) extends Distribution[Int] {

  def logDensity(b: Int): Real = {
    p.log * b + (Real.one - p).log * (1 - b)
  }

  val generator = Generator.from { (r, n) =>
    val pd = n.toDouble(p)
    val u = r.standardUniform
    if (u < pd) 1 else 0
  }

}

// first simulate some data from a logistic regression model
val r = new scala.util.Random(0)
val N = 1000
val beta0 = 0.1
val beta1 = 0.3
val x = (1 to N) map { i =>
  3.0 * r.nextGaussian
}
val theta = x map { xi =>
  beta0 + beta1 * xi
}
def expit(x: Double): Double = 1.0 / (1.0 + math.exp(-x))
val p = theta map expit
val y = p map (pi => if (r.nextDouble < pi) 1 else 0)
println(y.take(10))
println(x.take(10))

// now build and fit model
val model = for {
  beta0 <- Normal(0, 5).param
  beta1 <- Normal(0, 5).param
  _ <- Predictor.from{x: Double =>
      {
        val theta = beta0 + beta1 * x
        val p = Real(1.0) / (Real(1.0) + (Real(0.0) - theta).exp)
        //Binomial(p,1)
        Bernoulli(p)
      }
    }.fit(x zip y)
} yield (beta0, beta1)

//implicit val rng = RNG.default
implicit val rng = ScalaRNG(3)
val its = 10000
val thin = 5
//val out = model.sample(Walkers(100), 10000, its)
val out = model.sample(HMC(5), 10000, its*thin, thin)
//val out = model.toStream(HMC(5), 1000).take(its).map(_()).toArray
println(out.take(10))
```

```scala
import com.cibo.evilplot.plot._
import com.cibo.evilplot.colors._
import com.cibo.evilplot.plot.aesthetics.DefaultTheme._
{'fraction_non_alphanumeric': 0.1620583717357911, 'fraction_numerical': 0.043010752688172046, 'mean_word_length': 3.1608093716719914, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1841854', 'n_tokens_mistral': 1667, 'n_tokens_neox': 1491, 'n_words': 412}