Гаусcово распределение

[m_ax]

Вечер добрый!

Никто не подскажет алгоритм для моделирования Гаусова распределения?

В настоящий момент я использую следующую схему:

Код:

double rand(const double &lx, const double &rx)
{
    return (rx - lx) * rand()/RAND_MAX + lx;
}
//-----------------------------------------------------------------------------

double RandG(const double &Mean, const double &StdDev)
{
    double s = 0, u = 0;
    do {
        u = rand(-1.0, 1.0);
        s = pow(u, 2) + pow(rand(-1.0, 1.0), 2);
    } while (s > 1);
    return sqrt(-2 * log(s) / s) * u * StdDev + Mean;
}

Мне очень важна скорость, поскольку функция RandG вызывается в цикле очень много раз,
а реализация её (см. выше) слишком дорого обходится...

Помогите кто чем может 

[Rcus]

/*Открывает конспекты по моделированию случайных величин*/
есть другая формула для моделирования нормального распределения:
если alpha1 и alpha2 - два случайных значения на интервале [0;1), то
eta1 = sqrt(-2*ln(alpha1))*cos(2*pi*alpha2);
eta2 = sqrt(-2*ln(alpha1))*sin(2*pi*alpha2);
Два независимых значения с нормальным распределением с параметрами (0,1)
Собственно именно такой алгоритм использует Boost::Random
http://www.boost.org/doc/libs/1_38_0/boost/random/normal_distribution.hpp

[m_ax]

Спасибо, просветили!

Вариант приведённый мной был позаимствован из Delphi 

[m_ax]

Понравился мне там один отрывок из кода в Boost библиотеке:

Код:

if(!_valid) {
      _r1 = eng();
      _r2 = eng();
      _cached_rho = sqrt(-result_type(2) * log(result_type(1)-_r2));
      _valid = true;
    } else {
      _valid = false;
    }

else - чтоб наверняка 

[Rcus]

Нет там все правильно. При первом вызове генерируется два значения и записывается в кеш, выставляется флаг valid, возвращается первое значение. При втором вызове флаг снимается, возвращается второе значение.

[m_ax]

Да, сорри, там всё правильно...

Кстати, оба алгоритма: это фактически две реализации метода Box'а и Muller'а.
Если реализацию предложенную мной оформить аналогично в виде класса, то по скорости данный алгоритм чуть выигрывает...

Ладно, как гласит девиз компании Borland: Не изобретайте велосипед - унаследуйте его

[Rcus]

/*shrugs*/ Я пробовал обе реализации в тесте на заполнение массива 10^6 (вычисления в double).

Код:

main@rchome:~/projects/test_norm$ gcc main.c -O2 -lm
main@rchome:~/projects/test_norm$ ./a.out
19 time - VCL
12 time - Boost

main.c:

Код

C
#include <time.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <malloc.h>
 
inline double randCust(double lx, double rx)
{
        return (rx - lx) * rand()/RAND_MAX + lx;
}
 
inline double RandG(double Mean, double StdDev)
{
    double s = 0, u = 0;
    do {
        u = randCust(-1.0, 1.0);
        s = pow(u, 2) + pow(randCust(-1.0, 1.0), 2);
    } while (s > 1);
    return sqrt(-2 * log(s) / s) * u * StdDev + Mean;
}
 
int main(int argc, char *argv[])
{
    const int cycles = 100000000;
    double *values = malloc(sizeof(double)*cycles);
    time_t t = time(NULL);
    srand(t);
    int i = 0;
    for (i = 0; i < cycles; ++i) {
        values[i] = RandG(0,1);
    }
    printf("%d time - VCL\n",(int)(time(NULL)-t));
    t = time(NULL);
    for (i = 0; i < cycles; i += 2) {
        double rnd1 = (double)rand()/RAND_MAX;
        double rnd2 = (double)rand()/RAND_MAX;
        double m1 = sqrt(-2*log(rnd1));
        double m21 = cos(2*M_PI*rnd2);
        double m22 = sin(2*M_PI*rnd2);
        values[i] = m1 * m21;
        values[i+1] = m1 * m22;
    }
    printf("%d time - Boost\n",(int)(time(NULL)-t));
    return 0;
}
 

[m_ax]

Ну разумеется   
А попробуйте вот такую реализацию:

.h

Код

C++ (Qt)
#ifndef RANDG_H
#define RANDG_H
 
class RandG
{
public:
    RandG(const double &mean = 0.0, const double &sigma = 1.0);
    RandG(const RandG &other);
    double operator()() const;
    double mean() const;
    double sigma() const;
    void setMean(const double &mean);
    void setSigma(const double &sigma);
 
private:
    static const double weight;
    double urand() const;
    double m_mean;
    double m_sigma;
    mutable double r1, r2, s, rho;
    mutable bool valid;
};
 
#endif // RANDG_H
 

.cpp

Код

C++ (Qt)
#include <cmath>
#include <cstdlib>
#include <ctime>
#include "randg.h"
 
 
const double RandG::weight = 2.0/double(RAND_MAX);
 
RandG::RandG(const double &mean, const double &sigma)
    : m_mean(mean), m_sigma(sigma), valid(false)
{
    srand(time(0));
}
//-----------------------------------------------------------------------------
 
RandG::RandG(const RandG &other)
    : m_mean(other.mean()), m_sigma(other.sigma()), valid(false)
{
    srand(time(0));
}
//-----------------------------------------------------------------------------
 
double RandG::urand() const
{
    return double(rand()) * weight - 1.0;
}
//-----------------------------------------------------------------------------
 
double RandG::mean() const
{
    return m_mean;
}
//-----------------------------------------------------------------------------
 
double RandG::sigma() const
{
    return m_sigma;
}
//-----------------------------------------------------------------------------
 
void RandG::setMean(const double &mean)
{
    m_mean = mean;
}
//-----------------------------------------------------------------------------
 
void RandG::setSigma(const double &sigma)
{
    m_sigma = sigma;
}
//-----------------------------------------------------------------------------
 
double RandG::operator()() const
{
    if (!valid)
    {
        do {
            r1 = urand();
            r2 = urand();
            s = r1 * r1 + r2 * r2;
        } while (s > 1);
        rho = sqrt(-2.0 * log(s)/s);
        valid = true;
    } else {
        valid = false;
    }
    return rho * (valid ? r1 : r2) * m_sigma + m_mean;
}
//-----------------------------------------------------------------------------
 

[Rcus]

Понятно Генерим значения, отбрасывая ненужные (enwi подсказывает: 4/pi - 1 ~= 27.32%) и взамен получаем более простую формулу без синусов и косинусов. Для стандартных генераторов работает и то хорошо

[m_ax]

В общем как то так