SDS-2.2, Scalable Data Science

Archived YouTube video of this live unedited lab-lecture:

Introduction to Simulation

Core ideas in Monte Carlo simulation

• modular arithmetic gives pseudo-random streams that are indistiguishable from 'true' Uniformly distributed samples in integers from $\{0,1,2,...,m\}$
• by diving the integer streams from above by $m$ we get samples from $\{0/m,1/m,...,(m-1)/m\}$ and "pretend" this to be samples from the Uniform(0,1) RV
• we can use inverse distribution function of von Neumann's rejection sampler to convert samples from Uniform(0,1) RV to the following:
  • any other random variable
  • vector of random variables that could be dependent
  • or more generally other random structures:
    • random graphs and networks
    • random walks or (sensible perturbations of live traffic data on open street maps for hypothesis tests)
    • models of interacting paticle systems in ecology / chemcal physics, etc...

• https://en.wikipedia.org/wiki/Inversetransformsampling
• https://en.wikipedia.org/wiki/Rejection_sampling - will revisit below for Expoential RV

breeze.stats.distributions

Breeze also provides a fairly large number of probability distributions. These come with access to probability density function for either discrete or continuous distributions. Many distributions also have methods for giving the mean and the variance.

import breeze.stats.distributions._

val poi = new Poisson(3.0);

import breeze.stats.distributions._
poi: breeze.stats.distributions.Poisson = Poisson(3.0)

val s = poi.sample(5); // let's draw five samples - black-box

s: IndexedSeq[Int] = Vector(1, 4, 8, 0, 6)

Getting probabilities of the Poisson samples

s.map( x => poi.probabilityOf(x) ) // PMF

res36: IndexedSeq[Double] = Vector(0.14936120510359185, 0.16803135574154085, 0.008101511794681432, 0.049787068367863944, 0.05040940672246224)

val doublePoi = for(x <- poi) yield x.toDouble // meanAndVariance requires doubles, but Poisson samples over Ints

doublePoi: breeze.stats.distributions.Rand[Double] = MappedRand(Poisson(3.0),<function1>)

breeze.stats.meanAndVariance(doublePoi.samples.take(1000));

res40: breeze.stats.MeanAndVariance = MeanAndVariance(3.013999999999998,2.9787827827827824,1000)

(poi.mean, poi.variance) // population mean and variance

res41: (Double, Double) = (3.0,3.0)

Exponential random Variable

Let's focus on getting our hands direty with a common random variable.

NOTE: Below, there is a possibility of confusion for the term rate in the family of exponential distributions. Breeze parameterizes the distribution with the mean, but refers to it as the rate.

val expo = new Exponential(0.5);

expo: breeze.stats.distributions.Exponential = Exponential(0.5)

expo.rate // what is the rate parameter

res42: Double = 0.5

A characteristic of exponential distributions is its half-life, but we can compute the probability a value falls between any two numbers.

expo.probability(0, math.log(2) * expo.rate)

res43: Double = 0.5

expo.probability(math.log(2) * expo.rate, 10000.0)

res44: Double = 0.5

expo.probability(0.0, 1.5)

res45: Double = 0.950212931632136

The above result means that approximately 95% of the draws from an exponential distribution fall between 0 and thrice the mean. We could have easily computed this with the cumulative distribution as well.

1 - math.exp(-3.0) // the CDF of the Exponential RV with rate parameter 3

res47: Double = 0.950212931632136

Drawing samples from Exponential RV

val samples = expo.sample(2).sorted; // built-in black box - we will roll our own shortly in Spark

samples: IndexedSeq[Double] = Vector(0.5392432880681155, 2.0560411023149476)

expo.probability(samples(0), samples(1));

res48: Double = 0.32373622100122956

breeze.stats.meanAndVariance(expo.samples.take(10000)); // mean and variance of the sample

res49: breeze.stats.MeanAndVariance = MeanAndVariance(1.987723329390203,4.095923857112381,10000)

(1 / expo.rate, 1 / (expo.rate * expo.rate)) // mean and variance of the population

res50: (Double, Double) = (2.0,4.0)

Pseudo Random Numbers in Spark

import spark.implicits._
import org.apache.spark.sql.functions._

import spark.implicits._
import org.apache.spark.sql.functions._

val df = spark.range(1000).toDF("Id") // just make a DF of 100 row indices

df: org.apache.spark.sql.DataFrame = [Id: bigint]

df.show(5)

+---+
| Id|
+---+
|  0|
|  1|
|  2|
|  3|
|  4|
+---+
only showing top 5 rows

val dfRand = df.select($"Id", rand(seed=1234567) as "rand") // add a column of random numbers in (0,1)

dfRand: org.apache.spark.sql.DataFrame = [Id: bigint, rand: double]

dfRand.show(5) // these are first 5 of the 1000 samples from the Uniform(0,1) RV

+---+------------------+
| Id|              rand|
+---+------------------+
|  0|0.2289181799234461|
|  1|0.9756456161051732|
|  2|0.7781702473664945|
|  3|0.5585984240683788|
|  4|0.8305446150005453|
+---+------------------+
only showing top 5 rows

Let's use the inverse CDF of the Exponential RV to transform these samples from the Uniform(0,1) RV into those from the Exponential RV.

val dfRand = df.select($"Id", rand(seed=1234567) as "rand") // add a column of random numbers in (0,1)
               .withColumn("one",lit(1.0))
               .withColumn("rate",lit(0.5))

dfRand: org.apache.spark.sql.DataFrame = [Id: bigint, rand: double ... 2 more fields]

dfRand.show(5)

+---+------------------+---+----+
| Id|              rand|one|rate|
+---+------------------+---+----+
|  0|0.2289181799234461|1.0| 0.5|
|  1|0.9756456161051732|1.0| 0.5|
|  2|0.7781702473664945|1.0| 0.5|
|  3|0.5585984240683788|1.0| 0.5|
|  4|0.8305446150005453|1.0| 0.5|
+---+------------------+---+----+
only showing top 5 rows

val dfExpRand = dfRand.withColumn("expo_sample", -($"one" / $"rate") * log($"one" - $"rand")) // samples from expo(rate=0.5)

dfExpRand: org.apache.spark.sql.DataFrame = [Id: bigint, rand: double ... 3 more fields]

dfExpRand.show(5)

+---+------------------+---+----+------------------+
| Id|              rand|one|rate|       expo_sample|
+---+------------------+---+----+------------------+
|  0|0.2289181799234461|1.0| 0.5| 0.519921578060948|
|  1|0.9756456161051732|1.0| 0.5| 7.430086817819349|
|  2|0.7781702473664945|1.0| 0.5|3.0116901426840474|
|  3|0.5585984240683788|1.0| 0.5|1.6356004297263063|
|  4|0.8305446150005453|1.0| 0.5|3.5503312043026414|
+---+------------------+---+----+------------------+
only showing top 5 rows

display(dfExpRand)

Truncated to 30 rows

dfExpRand.describe().show() // look sensible

+-------+-----------------+--------------------+----+----+--------------------+
|summary|               Id|                rand| one|rate|         expo_sample|
+-------+-----------------+--------------------+----+----+--------------------+
|  count|             1000|                1000|1000|1000|                1000|
|   mean|            499.5| 0.49368134205225334| 1.0| 0.5|    1.98855225198386|
| stddev|288.8194360957494|  0.2925326105055967| 0.0| 0.0|   2.077189631386989|
|    min|                0|6.881987320686012E-4| 1.0| 0.5|0.001376871299039548|
|    max|              999|  0.9999092082356841| 1.0| 0.5|  18.613883955674265|
+-------+-----------------+--------------------+----+----+--------------------+