Бокса-Мюллера генератор гауссовских случайных чисел и график его - PullRequest
0 голосов
/ 23 мая 2018

Я пытаюсь написать код на python, используя уравнения Бокса-Мюллера, но я не знаю, как начать!

Вот пример, который я пытаюсь решить:

  • Пик, наблюдаемый при 900 кэВ, показывает FWHM 2 кэВ.Используя метод гауссовой выборки, указанный ниже, сгенерируйте 15 000 отсчетов, соответствующих пику 900 кэВ, и сохраните энергии выборки.

  • Создайте и построите гистограмму с шириной ячейки 0,2 кэВ и сравните сфункция Гаусса с той же площадью пика.

  • Используя программное обеспечение для анализа данных, попробуйте гауссово подгонку к данным Монте-Карло и убедитесь, что результат достаточно близок к модели пика.

Метод выборки Гаусса-Бокса-Мюллера: [Обратите внимание, две приведенные ниже выборочные переменные y1, y2 относятся к единице распределения Гаусса (т. Е. Mu = 0, segma = 1).

y1 = (-2 ln r1)^1/2  * cos(2pi*r2)
y2 = (-2 ln r1)^1/2 * sin(2pi*r2)

   (r1, r2: random numbers)}

Есть предложения?

* Обновление *

Я получаю сообщение об ошибке:

g1 = BoxMuller (v) NameError: имя 'v' не определено

Используемый код:

import random    
import matplotlib.pyplot as plt    
import numpy as np

def BoxMuller():    
    r1 = np.random.randn(15000)*10    
    r2 = np.random.randn(15000)    
    a = 2.0 * np.pi * r2        
    v = np.sqrt( -2.0*np.log(1.0 - r1)) * np.sin(a)        
    u = np.sqrt( -2.0*np.log(1.0 - r1)) * np.cos(a)

g1 = BoxMuller(v)
g2 = BoxMuller(u)
q = 900.0 + g1*2.0
k = 900.0 + g2*2.0
plt.hist(q, k)
plt.show()

1 Ответ

0 голосов
/ 23 мая 2018

Что ж, вот простая реализация для запуска и повозки с

import math
import random

def BoxMuller():
    r1 = random.random()
    r2 = random.random()

    a  = 2.0 * math.pi * r1
    v  = math.sqrt( -2.0*math.log(1.0 - r2))

    return (v * math.sin(a), v * math.cos(a))


g1, g2 = BoxMuller()

q = 900.0 + g1*2.0
...

ОБНОВЛЕНИЕ

Очевидно, что дано FWHM, а не std.dev.Чтобы получить сигму, нужно разделить FWHM на 2*sqrt(2*log(2)) ~ 2.355.Таким образом, код выборки должен быть

FWHM = 2.0
q = 900.0 + g1 * FWHM/2.355
Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...