Генератор нормальных случайных чисел

[tyomitch]

Ни у кого под рукой нету готовой реализации на VB?

[Sebas]

ты что, магией решил заняться? К чему такое стремление?

[tyomitch]

Вон, спортировал что нашёл.
Надо в Кирпичи?

Код: Выделить всё

Option Explicit

Private Sub Form_DblClick()
Dim x As Long, y As Long, c As Long
For y = 1 To 100
    For x = 1 To 100
        '1..255: m=128, 3s=127, s=42
        c = Norm1 * 42 + 127
        If c < 0 Then c = 0: If c > 255 Then c = 255
        PSet (x, y), RGB(c, c, c)
    Next
Next

For y = 1 To 100
    For x = 101 To 200
        c = Norm2 * 42 + 127
        If c < 0 Then c = 0: If c > 255 Then c = 255
        PSet (x, y), RGB(c, c, c)
    Next
Next

For y = 101 To 200
    For x = 1 To 100
        c = Norm3 * 42 + 127
        If c < 0 Then c = 0: If c > 255 Then c = 255
        PSet (x, y), RGB(c, c, c)
    Next
Next

For y = 101 To 200
    For x = 101 To 200
        c = Norm4 * 42 + 127
        If c < 0 Then c = 0: If c > 255 Then c = 255
        PSet (x, y), RGB(c, c, c)
    Next
Next

For y = 1 To 100
    For x = 201 To 300
        c = Rnd * 256
        PSet (x, y), RGB(c, c, c)
    Next
Next
For y = 101 To 200
    For x = 201 To 300
        Dim r As Long, g As Long, b As Long
        r = Norm2 * 42 + 127
        If r < 0 Then r = 0: If r > 255 Then r = 255
        g = Norm3 * 42 + 127
        If g < 0 Then g = 0: If g > 255 Then g = 255
        b = Norm4 * 42 + 127
        If b < 0 Then b = 0: If b > 255 Then b = 255
        PSet (x, y), RGB(r, g, b)
    Next
Next
End Sub

Private Sub Form_Load()
Randomize
ScaleMode = vbPixels
Exit Sub

Dim i As Long, Start As Double
Start = Timer
For i = 1 To 1000000
    Norm1
Next
Debug.Print Timer - Start

Start = Timer
For i = 1 To 1000000
    Norm2
Next
Debug.Print Timer - Start

Start = Timer
For i = 1 To 1000000
    Norm3
Next
Debug.Print Timer - Start

Start = Timer
For i = 1 To 1000000
    Norm4
Next
Debug.Print Timer - Start
End Sub


''from comp.lang.c FAQ, question 13.20 (http://c-faq.com/lib/gaussian.html)

'Central Limit Theorem
Function Norm1() As Double
    Norm1 = Rnd + Rnd + Rnd + Rnd + Rnd + Rnd + Rnd + Rnd + Rnd + Rnd + Rnd + Rnd - 6
End Function

'Abramowitz, Stegun, Box, Muller
Function Norm2() As Double
Const PI = 3.14159265358979
Static U As Double, V As Double
Static phase As Boolean

    If phase Then
        Norm2 = U * Cos(V)
    Else
        U = Sqr(-2 * Log((Rnd + 0.0000001) / 1.0000001))
        V = 2 * PI * Rnd
        Norm2 = U * Sin(V)
    End If

    phase = Not phase
End Function

'Knuth, 3.4.1.C.1
Function Norm3() As Double
Static V1 As Double, V2 As Double, S As Double
Static phase As Boolean
   
    If phase Then
        Norm3 = V2 * S
    Else
        Do
            V1 = 2 * Rnd - 1
            V2 = 2 * Rnd - 1
            S = V1 * V1 + V2 * V2
        Loop Until (S > 0) And (S < 1)
        S = Sqr(-2 * Log(S) / S)
        Norm3 = V1 * S
    End If

    phase = Not phase
End Function


'http://www.enlight.ru/demo/faq/smth.phtml?query=alg_rnd
Function Norm4() As Double
Const C1 = 0.029899776, C2 = 0.008355968, C3 = 0.076542912, C4 = 0.252408784, C5 = 3.949846138
Dim S As Double: S = Norm1 / 4
Dim r As Double: r = S * S
    Norm4 = ((((C1 * r + C2) * r + C3) * r + C4) * r + C5) * S
End Function


[uhm]

Уук...

1ый метод - лажа полная. ЦПТ говорит, вообще-то, о пределе при n стремящемся к бесконечности, а не к 12

Про остальные - интересно было бы прочитать теоретическую основу. А самый прозрачный метод, на мой взгляд, ты пропустил:

http://c-faq.com/lib/gaussrand.luben.html

(берется обратная функция плотности нормального распределения от равномерно распределенной случайной величины).

[tyomitch]

1ый метод - лажа полная. ЦПТ говорит, вообще-то, о пределе при n стремящемся к бесконечности, а не к 12

Ну, бесконечную сумму я сосчитать ведь не могу? Сколько смог, столько и сосчитал.

Про остальные - интересно было бы прочитать теоретическую основу.

Третий описан в Кнуте. Двум других я объяснения не знаю сам; но они встречаются во многих местах, и, кроме того, результаты у всех четырёх получаются визуально схожие.

(берется обратная функция плотности нормального распределения от равномерно распределенной случайной величины).

А какая у этого способа теоретическая основа? Обычно-то ведь обращают функцию распределения, а не плотность.

[uhm]

12 ближе к 1, чем к +бесконечности

Шутка, конечно, но в общем лучше этот метод не использовать.

А с функцией плотности я как-то погорячился, обычно обращают действительно функцию распределения.

[tyomitch]

12 ближе к 1, чем к +бесконечности

Ну хорошо, а какое число тогда ближе к +бесконечности, чем к 1?

[uhm]

Вопрос, вообще-то, вполне осмысленный (в контексте данной задачи, разумеется) Я уже просто всю теорию помню достаточно смутно, но насколько я понимаю, вполне можно оценить, насколько сильно распределение, полученное сложением n равномерно распределенных величин, отличается от нормального. Я к тому, что при n=12 погрешность будет весьма велика...

Эта штука (как, кажется, было написано в источнике) хреново работает с хвостами. При n=12, вероятность того, что полученная величина будет больше 6 - строго 0. Согласно экселевской функции распределения, она где-то в районе 0,000000001. В принципе, кстати, не так уж и много. Anyway, если интересно моделировать выбросы, эта штука не подходит.

[tyomitch]

А самый прозрачный метод, на мой взгляд, ты пропустил:

http://c-faq.com/lib/gaussrand.luben.html

Дошли руки проверить: на глаз видно, что распределение левое
Любена фтопку!




Кроме того, посмотрел характеристики моих четырёх распределений.

Код: Выделить всё

Norm1:         1,938
Mean:         -0,0004
Variance:      1,0006
Skewness:      0,0059
Kurtosis:     -0,0848

Norm2:         1,1406
Mean:         -0,0001
Variance:      0,9989
Skewness:      0,0002
Kurtosis:     -0,0081

Norm3:         1,3288
Mean:          0
Variance:      1,0016
Skewness:      0,0014
Kurtosis:      0,0111

Norm4:         2,532
Mean:         -0,0007
Variance:      1,0004
Skewness:      0,0005
Kurtosis:      0,0266


Norm2 и самый быстрый, и самый точный.

[uhm]

Ок, значит, его и будем использовать

Ты эксцесс считал с вычитанием 3?

[Денис Победря]

А что значит нормальных?

[uhm]

Это название класса распределений. Их также называют распределениями Гаусса. Более подробно - в учебниках по теории вероятности и мат. статистике.

[tyomitch]

Ок, значит, его и будем использовать

Ты эксцесс считал с вычитанием 3?

Да, конечно.

--------------------
Ещё я сегодня пересдавал теорию вероятности, и не пересдал: не смог решить задачу. Причём до сих пор не могу её решить.

"Произведено 1000 независимых испытаний с тремя равновероятными исходами А, Б и В. Какой коэффициент корреляции между числом исходов А и числом исходов Б?"

Я пробовал много разных подходов, и каждый раз получал новый неправильный результат

[GSerg]

Ноль?

[tyomitch]

Нет, минус половина.

Но интересно само решение.

[Amed]

А как это объяснить физически?
Исходы равновероятны, P(A) = P(Б) = P(В) = 1/3, как я понимаю?

[tyomitch]

Сами исходы равновероятны, но количества исходов разных видов не независимы. А именно, их сумма задана жёстко.

[uhm]

Э-э, я правильно понимаю, что если величин было бы две, то корреляция была бы (-1) ?

[tyomitch]

uhm, совершенно правильно.

[tyomitch]

Решил набить руку с набором формул в ТеХе; решение той задачи, если кому интересно, тут. Исходник тут.

[tyomitch]

Вчера пересдавал ещё раз, и в этот раз пересдал успешно

Что более интересно, препод показал мне короткое решение той задачи.
Вот оно целиком:

ξ+η+ζ = 1000
0 = D(ξ+η+ζ) = 3D + 6cov
cov = -½D; r = -½