# 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

The scaladoc is available here.

### Random variables as monads

`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.