Генерация случайного числа из логнормального распределения

Это не должно быть слишком сложно, но по какой-то причине такая выборка случайных чисел из распределения действительно сбивает меня с толку.

Я знаю, что лучшими вариантами генерации случайных чисел из дистрибутива являются библиотеки boost / C ++ 11 … к сожалению, я не могу заставить этот код компилироваться с c ++ 0x, и в любом случае я бы предпочел сохранить совместимость на сервер, который я также использую, который работает под управлением gcc 4.1.2 — древний, я знаю, не поддерживает более новый C ++. Разочарование. И как всегда, временное ограничение означает, что мне нужно сделать все возможное, чтобы быстро исправить ситуацию.

Моя следующая опция — взять экспоненту случайного числа из формул уравнения Мюллера, но я не получаю логнормальное распределение с указанными параметрами. Я не понимаю, почему это не работает.

Любая помощь будет принята с благодарностью!

void testRNG(){

int mean = 5000;
int std = 50;

ofstream out("./Output/normal_samples.out");

RunningStats normal;
for (int i=0;i<2000;++i){
double sample = randomSample(mean, std, NORMAL);//call function with box muller transformation to return a number from a normal distriubtion
out<<sample<<endl;
normal.Push(sample);//keep a running average of sampled numbers
}
cout<<"Normal Mean = "<<normal.Mean()<<endl;
cout<<"Normal Std = "<<normal.StandardDeviation()<<endl;

RunningStats lognormal;
for (int i=0;i<2000;++i){
double sample = randomSample(mean, std, LOGNORMAL);
out<<sample<<endl;
lognormal.Push(sample);
}
cout<<"Lognormal Mean = "<<lognormal.Mean()<<endl;
cout<<"Lognormal Std = "<<lognormal.StandardDeviation()<<endl;
}

Функции сэмплирования, которые я не написал, сначала переходят к case-объекту из randomSample (), а затем вызывают:

РЕДАКТИРОВАТЬ — я заметил, что на самом деле он вызывал функцию, чтобы найти логнормальные параметры. Добавлено в.

double randNormal(double mean, double stdev) {
static long numSamples = 0;
static double Z2;
if ((numSamples++ & 1) == 0) {
double Z1, U1, U2;
do { U1 = randUniform(0, 1); } while (U1 <= 0 || U1 >= 1);
do { U2 = randUniform(0, 1); } while (U1 <= 0 || U1 >= 1);
Z1 = sqrt(-2 * log(U1)) * cos(6.28318531 * U2);
Z2 = sqrt(-2 * log(U1)) * sin(6.28318531 * U2);
return mean + stdev * Z1;
} else {
return mean + stdev * Z2;
}
}

double randLognormal(double mu, double sigma) {
return exp(randNormal(mu, sigma));
}

double randLognormalMeanStdev(double mean, double stdev) {
return randLognormal( log(mean) - 0.5 * log(1 + (stdev * stdev) / (mean * mean)) , log(1 + (stdev * stdev) / (mean * mean)));
}

Таким образом, я получаю вывод:

Normal Mean = 4998.72 //I said 5000
Normal Std = 49.7054 //I said 50
Lognormal Mean = 4999.74
Lognormal Std = 0.492766 //this is the part that is not working

Чего мне не хватает, чтобы получить логнормальный стандарт как то, что я хочу?

Другие варианты также будут оценены — может быть, есть что-то еще, что я пропускаю.

Заранее спасибо!

Изменить — я понял, что я должен был дать понять, что мне нужно выбрать из логнормального распределения