Название: Гаусcово распределение
Отправлено: m_ax от Апрель 03, 2009, 19:30
Вечер добрый! Никто не подскажет алгоритм для моделирования Гаусова распределения? В настоящий момент я использую следующую схему: 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 вызывается в цикле очень много раз, а реализация её (см. выше) слишком дорого обходится... Помогите кто чем может :)
Название: Re: Гаусово распределение
Отправлено: Rcus от Апрель 03, 2009, 21:07
/*Открывает конспекты по моделированию случайных величин*/
есть другая формула для моделирования нормального распределения:
если 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
Название: Re: Гаусово распределение
Отправлено: m_ax от Апрель 03, 2009, 21:21
Спасибо, просветили!
Вариант приведённый мной был позаимствован из Delphi :)
Название: Re: Гаусово распределение
Отправлено: m_ax от Апрель 04, 2009, 00:50
Понравился мне там один отрывок из кода в Boost библиотеке: if(!_valid) {
_r1 = eng();
_r2 = eng();
_cached_rho = sqrt(-result_type(2) * log(result_type(1)-_r2));
_valid = true;
} else {
_valid = false;
}
else - чтоб наверняка ;D
Название: Re: Гаусово распределение
Отправлено: Rcus от Апрель 04, 2009, 06:48
Нет там все правильно. При первом вызове генерируется два значения и записывается в кеш, выставляется флаг valid, возвращается первое значение. При втором вызове флаг снимается, возвращается второе значение.
Название: Re: Гаусово распределение
Отправлено: m_ax от Апрель 07, 2009, 18:44
Да, сорри, там всё правильно...
Кстати, оба алгоритма: это фактически две реализации метода Box'а и Muller'а.
Если реализацию предложенную мной оформить аналогично в виде класса, то по скорости данный алгоритм чуть выигрывает...
Ладно, как гласит девиз компании Borland: Не изобретайте велосипед - унаследуйте его ;)
Название: Re: Гаусово распределение
Отправлено: Rcus от Апрель 07, 2009, 19:18
/*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;
}
Название: Re: Гаусово распределение
Отправлено: m_ax от Апрель 07, 2009, 19:42
Ну разумеется ;D А попробуйте вот такую реализацию: .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;
}
//-----------------------------------------------------------------------------
Название: Re: Гаусово распределение
Отправлено: Rcus от Апрель 07, 2009, 20:50
Понятно :) Генерим значения, отбрасывая ненужные (enwi подсказывает: 4/pi - 1 ~= 27.32%) и взамен получаем более простую формулу без синусов и косинусов. Для стандартных генераторов работает и то хорошо :)
Название: Re: Гаусово распределение
Отправлено: m_ax от Апрель 07, 2009, 22:06
В общем как то так :)
|