Update 2023-09

2023-09-23 09:41:04 +03:00

3.6 MiB

Raw Permalink Blame History

In [1]:

@file:Repository("https://repo.kotlin.link")
@file:DependsOn("space.kscience:plotlykt-jupyter:0.5.0")

In [2]:

Plotly.jupyter.notebook()

Plotly notebook integration switch into the legacy mode.

Центральная предельная теорема¶

Центральная предельная теорема в вльной формулировке гласит, что сумма случайных величин, имеющих распределение с конечным средним и дисперсией, распределена по нормальному распределению со средним, равным сумме средних компонент и дисперсией, равно сумме дисперсий.

На доказательство этой теоремы в общем виде обычно уходит львиная доля семестра в курсе теории вероятности. В данной задаче мы попробуем проверить эту теорему на практике.

Для этого мы возьмем некоторое распределение и посчитаем численно распределение суммы.

Для решения задачи нам потребуется:

Определить функцию плотности вероятности распределения
Создать функцию для расчета распределения суммы двух случайных величин. Распределение суммы случайных величин вычисляется как свертка распределений (см. ниже).
Последовательно посчитатать распределение суммы нескольких величин a + a + a + a = (((a + a) + a) +a)
Сравнить результат с нормальным распределением.

Для решения задачи может понадобится представление функции в виде набора точек при помощи процедуры линейной интерпояции.

Свертка¶

Сверткой) двух функций называется следующее интегральное преобразование: $$ h(x) = f(x) \otimes g(x) = \int {f(x-y) g(y) dy } $$ Интегрирование ведется по всей области определения функций. Очевидно, что переменные под интегрированием можно менять местами. Для вычисления интеграла можно использовать метод прямоугольников или метод трапеций.

Интерполяция¶

Для помощи в решении задачи можно использовать готовую реализацию линейной интерполяции.

In [3]:

typealias Function = (Double) -> Double

fun pointsFromRange(from: Double, to: Double, numPoints: Int = 100): List<Double>{
    require( to > from)
    require(numPoints > 2)
    val step = (to - from)/(numPoints - 1)
    require(step>0)
    return List(numPoints){i -> from + step * i}
}

fun Function.interpolate(from: Double, to: Double, numPoints: Int = 100): Function{
    //compute xs
    val xs = pointsFromRange(from, to, numPoints) 
    //compute function values
    val ys = xs.map{ invoke(it) }
    
    // return an anonimous function
    return { x ->
        // check interpolation region. Zerop outside the region
        if(x < xs.first() || x > xs.last()) 0.0
        // find section number
        val num: Int = xs.indexOfFirst { it > x }
        // num >=1
        if(num <= 0){
            0.0
        } else{
            //return the result as last expression
            ys[num - 1] + (ys[num] - ys[num - 1])/(xs[num] - xs[num - 1])*(x - xs[num - 1])
        }
    }
}

Решение¶

Проверка¶

In [4]:

fun l2RelativeDiff(res: Function, order: Int, baseDispersion: Double): Double{
    val numPointForDiff = 100
    val dispersion = baseDispersion * order
    val sigma = sqrt(dispersion)
    val xs = ArrayList<Double>()
    for(i in 0 until numPointForDiff){
        xs.add(-3*sigma + i.toDouble()/(numPointForDiff-1)*6*sigma)
    }

    fun normal(x: Double): Double{
        return 1.0/sqrt(2* PI*dispersion)*exp(-x.pow(2)/2.0/dispersion)
    }

    return xs.sumOf{ x-> 
        (res(x)/normal(x) - 1.0).pow(2)
    }
}

In [5]:

fun integrate(f: Function, from: Double, to: Double, numPoints: Int = 100): Double{
    val xs = pointsFromRange(from, to, numPoints)
    return xs.sumOf { f(it) } * (to - from) / numPoints
}

fun convolve(f: Function, g: Function, from: Double, to: Double, numPoints: Int = 100): Function{
    val res: Function = { x ->
        val integrand: Function = { y -> f(x - y) * g(y)}
        integrate(integrand, from, to, numPoints)
    }

    //We need interpolation here to limit complexity
    //With interpolation it linear, without it is exponential
    return res.interpolate(from, to, numPoints)
}

In [6]:

val f: Function = {x ->  if(x in (-1.0..1.0)) 0.5  else 0.0}

val from = -5.0
val to = 5.0

val num = 300

val first: Function = convolve(f, f, from, to, num)
val second: Function = convolve(first, f, from, to, num)
val third: Function = convolve(second, f, from, to, num)
val fourth: Function = convolve(third, f, from, to, num)
val fifth: Function = convolve(fourth, f, from, to, num)

In [7]:

val xs = pointsFromRange(from, to, num)

Plotly.plot { 

    scatter { 
        name = "initial"
        x.numbers = xs
        y.numbers = xs.map(f)
    }

    scatter { 
        name = "first"
        x.numbers = xs
        y.numbers = xs.map(first)
    }

    scatter { 
        name = "second"
        x.numbers = xs
        y.numbers = xs.map(second)
    }

    scatter { 
        name = "fifth"
        x.numbers = xs
        y.numbers = xs.map(fifth)
    }

    val d = 6.0/3.0


    scatter { 
        name = "normal"
        x.numbers = xs
        y.numbers = xs.map{1.0/sqrt(2*PI*d) * exp(- it.pow(2)/d/2) }
    }
 }

</html>

In [8]:

l2RelativeDiff(fifth, 5,1.0)

53.72510082115637

In [ ]:

3.6 MiB Raw Permalink Blame History Unescape Escape