normal_distribution_fitting_mle.kt

normal_distribution_fitting_mle.kt

import java.util.concurrent.ThreadLocalRandom
import kotlin.math.exp
import kotlin.math.ln
import kotlin.math.pow


fun main() {

    val observations = doubleArrayOf(
            1.0,
            1.0,
            1.0,
            2.0,
            2.0,
            3.0,
            2.0,
            2.0,
            1.0,
            1.0,
            1.0
    )

    var bestLikelihood = Double.MIN_VALUE
    var stdDev = 1.0
    var mean = 0.0

    repeat(100_000) {

        val meanAdj = ThreadLocalRandom.current().nextDouble(-.05, .05)
        val stdDevAdj = ThreadLocalRandom.current().nextDouble(-.05, .05)

        mean += meanAdj
        stdDev += stdDevAdj

        val likelihood = observations.asSequence().map { normalDistributionPDF(it, mean, stdDev) }
                .map { ln(it) }
                .sum()
                .let(::exp)

        if (likelihood > bestLikelihood) {
            bestLikelihood = likelihood
        } else {
            mean -= meanAdj
            stdDev -= stdDevAdj
        }
    }


    println("mean=$mean, stdDev=$stdDev")

}

# https://en.wikipedia.org/wiki/Normal_distribution
fun normalDistributionPDF(x: Double, mean: Double, stdDev: Double) =
        (1.0 / (2.0 * Math.PI * stdDev.pow(2)).pow(0.5)) * exp(-1.0 * ((x - mean).pow(2) / (2.0 * stdDev.pow(2))))