# jrj-d / cardano  0.2.3

Models random variables and bayesian nets using monads in Scala

# cardano

Models random variables and bayesian nets using monads in Scala. Inference is performed by sampling.

## Installation

Include this dependency in your `sbt` configuration:

``````"org.jrj-d" %% "cardano" % "0.2.3"
``````

## Usage

`cardano` models random variable as monads. This should feel familiar to most Scala programmers since the `collection` library is built this way.

The base trait is `Stochastic[A]` which defines a random variable of concrete type `A`.

```scala> def die: Stochastic[Int] = Stochastic.discreteUniform(6).map(_ + 1)
die: cardano.Stochastic[Int] = ...```

The above line creates a uniform discrete random variable from 0 to 5 (`Stochastic.discreteUniform(6)`). Applying `map(_ + 1)` on it effectively creates a fair die with values from 1 to 6.

Let's test our die by computing its expectation and the probability of yielding a 1.

```scala> die.expectation
res0: Double = 3.554

scala> die.map(_ == 1).expectation
res1: Double = 0.156

scala> die.map(_ == 1).expectation(10000)
res2: Double = 0.1667```

The number of samples defaults to 1000 when no argument is given.

The `filter` operation allows us to compute conditional probabilities. The expectation of a die given that the value is less than 3:

```scala> die.filter(_ < 3).expectation
res16: Double = 1.504```

To model the sum of two dice, we will combine two dice using `flatMap`:

```scala> val sumDice = die.flatMap(v => die.map(_ + v))
sumDice: cardano.Stochastic[Int] = ...

scala> sumDice.expectation
res21: Double = 7.062

scala> sumDice.std
res23: Double = 2.40265041152474```

That's it! Along with `Stochastic.constant(...)` that allows us to build constant random variables, all the standard monadic Scala operations are here.

Let's combine all of this in a solution to the Monty Hall problem. Thanks to Scala syntactic sugar, we can define random variables in for loops:

```scala> :paste
// Entering paste mode (ctrl-D to finish)

val keepStrategy = for(
winningDoor <- Stochastic.discreteUniform(3);
firstChosenDoor <- Stochastic.discreteUniform(3);
revealedDoor <- Stochastic.discreteUniform(3)
if revealedDoor != firstChosenDoor && revealedDoor != winningDoor
) yield winningDoor == firstChosenDoor

val changeStrategy = for(
winningDoor <- Stochastic.discreteUniform(3);
firstChosenDoor <- Stochastic.discreteUniform(3);
revealedDoor <- Stochastic.discreteUniform(3)
if revealedDoor != firstChosenDoor && revealedDoor != winningDoor
) yield (winningDoor != firstChosenDoor && winningDoor != revealedDoor)

// Exiting paste mode, now interpreting.

keepStrategy: cardano.Stochastic[Boolean] = ...
changeStrategy: cardano.Stochastic[Boolean] = ...

scala> keepStrategy.expectation
res65: Double = 0.311

scala> changeStrategy.expectation
res66: Double = 0.652```

Indeed, if the player changes door after the host reveals a loosing door, they have a probability of winning of 2/3, 1/3 otherwise.

### Duplicating random variables

What if we want to study the rolls of 5 dice? Chaining `flatMap` would be tedious.

Method `repeat` of `Stochastic[A]` allows us to duplicate a random variable. The duplicates are identically distributed and independent.

`repeat` takes as input a function of signature `(=> A) => F[A]` for some higher-kinded type `F`. Let's create a random variable of a sequence of 5 dice:

```scala> val fiveDice = die.repeat(Seq.fill[Int](5))
fiveDice: cardano.Stochastic[Seq[Int]] = ...

scala> fiveDice.sample
res24: Seq[Int] = List(6, 2, 6, 4, 5)

scala> fiveDice.sample
res25: Seq[Int] = List(1, 2, 2, 6, 5)

scala> fiveDice.map(s => s.distinct.size == s.size).expectation
res40: Double = 0.0933```

The last line gives an estimate of the probability that all dice have a different value.

`repeat` allows us more exotic constructs like an infinite sequence of die rolls:

```scala> val dice = die.repeat(Stream.continually[Int])
dice: cardano.Stochastic[scala.collection.immutable.Stream[Int]] = ...

scala> dice.sample.take(10).toList
res46: List[Int] = List(5, 2, 5, 6, 1, 5, 2, 6, 3, 3)

scala> dice.sample.take(10).toList
res47: List[Int] = List(4, 4, 2, 4, 2, 6, 4, 5, 4, 5)```

This allows us to estimate the number of dice required to reach a given sum, say 30:

```scala> val runningSum = dice.map(_.scan(0)(_ + _))
runningSum: cardano.Stochastic[scala.collection.immutable.Stream[Int]] = ...

scala> runningSum.sample.take(10).toList
res48: List[Int] = List(0, 6, 10, 16, 21, 22, 27, 32, 37, 42)

scala> val nbDice = runningSum.map(_.takeWhile(_ < 30).length)
nbDice: cardano.Stochastic[Int] = ...

scala> nbDice.expectation
res49: Double = 9.014

scala> nbDice.std
res51: Double = 1.4758045263516417```

### Markov chain

Method `markov` of `Stochastic` creates a Markov chain of the first order. It takes as input a function indicating the relationship between two consecutive terms.

Here is a model of the Californian weather (90% chance of a sunny day after a sunny day, 50% after a rainy day):

```scala> :paste
// Entering paste mode (ctrl-D to finish)

sealed trait Weather
case object Sunny extends Weather
case object Rainy extends Weather

def boolToWeather(b: Boolean): Weather = if(b) Sunny else Rainy

def transition(weather: Weather): Stochastic[Weather] = weather match {
case Sunny => Stochastic.coin(0.9).map(boolToWeather)
case Rainy => Stochastic.coin(0.5).map(boolToWeather)
}

val weatherChain: Stochastic[Stream[Weather]] = Stochastic.constant[Weather](Sunny).markov(transition)

// Exiting paste mode, now interpreting.

defined trait Weather
defined object Sunny
defined object Rainy
boolToWeather: (b: Boolean)Weather
transition: (weather: Weather)cardano.Stochastic[Weather]
weatherChain: cardano.Stochastic[Stream[Weather]] = ...

scala> weatherChain.sample.take(10).toList
res54: List[Weather] = List(Sunny, Sunny, Rainy, Sunny, Sunny, Sunny, Sunny, Sunny, Sunny, Rainy)```

### Bayesian inference

Sampling is often used to compute posteriors in intractable cases. For the sake of brevity, we will work here with a toy example that is in fact tractable.

We want to determine the bias of a coin by repeatedly flipping it. The prior is a Beta distribution with a rather flat belief in dealing with a fair coin (2, 2).

```scala> :paste
// Entering paste mode (ctrl-D to finish)

val unknownCoin = Stochastic.coin(0.6)

val prior = Stochastic.beta(2, 2)

val observations = Seq.fill(100)(unknownCoin.sample)

val posterior = Stochastic.posterior(prior, observations) { (bias, observation) =>
if(observation) bias else 1 - bias
}

// Exiting paste mode, now interpreting.

unknownCoin: cardano.Stochastic[Boolean] = ...
prior: cardano.Stochastic[Double] = ...
observations: Seq[Boolean] = ...

scala> posterior.expectation
res55: Double = 0.6051810933353833

scala> posterior.std
res56: Double = 0.0464702905875049```

See also `posteriorByLog` where the log-likelihood is provided by the user.

### MCMC sampling

`cardano` includes the Metropolis-Hastings and Metropolis sampling procedures, along with a facility to build maximum entropy distributions. This can be especially useful when you only have access to the unnormalized density of the distribution.

Suppose we don't have `Stochastic.beta` but we want to create a Beta random variable. Metropolis-Hastings is a sampling procedure that builds a Markov chain whose equilibrium distribution is the target distribution. It takes as input as a transition function and an unnormalized target density.

The transition function consists in adding a small uniform noise to the current value and capping it to [0,1]. The unnormalized target distribution can be found on Wikipedia.

```scala> :paste
// Entering paste mode (ctrl-D to finish)

def logDensity(a: Double, b: Double, value: Double) = (a - 1) * math.log(value) + (b - 1) * math.log(1 - value)

def cap(value: Double) = math.min(1.0, math.max(0.0, value))

val metropolisBeta = Stochastic.metropolis(Stochastic.constant(0.5))(logDensity(20, 10, _)) { lastValue =>
Stochastic.continuousUniform.map(noise => cap(lastValue + 0.1 * noise - 0.05))
}

val exactBeta = Stochastic.beta(20, 10)

// Exiting paste mode, now interpreting.

logDensity: (a: Double, b: Double, value: Double)Double
cap: (value: Double)Double
metropolisBeta: cardano.Stochastic[Double] = ...
exactBeta: cardano.Stochastic[Double] = ...

scala> metropolisBeta.expectation
res61: Double = 0.6743298870740473

scala> exactBeta.expectation
res62: Double = 0.6646562222468744

scala> metropolisBeta.std
res63: Double = 0.08492853676864531

scala> exactBeta.std
res64: Double = 0.0850264034989098```

Refer to the scaladoc to set the burn-in period and the interval between two consecutive samples.

Note that the Bayesian inference capabilities of `cardano` exposed through methods `posterior` and `posteriorByLog` builds upon those MCMC sampling techniques.