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

[hp114frost]

Добрый день. Из предложенных блок-схем было решено выбрать правую. Привожу полный код функций по генерации распределения. Генерируемые значения практически адекватны, единственный недостаток - не сильно разнятся от раза к разу, и пожалуй малы, относительно того, что генерирует MatLab (например в matlab с достаточной частотой получается 7, у меня 3,при lambda =4, но это случайные значения, понимаю допуски). Для генерации использовал srand(), просьба проверить, можно ли так делать, как приведено в моем коде и посоветовать более изящное решение, если оно имеется (безусловно имеется). Большое спасибо за помощь и оперативные реакции, я благодарен всем участникам дискуссии. Надеюсь приведенный код поможет кому-либо еще, хотя он безусловно требует доработки.
#include <cstdlib>
#include <iostream>
#include <math.h>
#include <time.h>

using namespace std;

inline double closed_interval_rand(double x0, double x1)
{
return x0 + (x1 - x0) * rand() / ((double) RAND_MAX);
}

double srandd()
{
double pass;
srand(time(0));
pass = closed_interval_rand(0, 1);
printf(" %f",pass);
return pass;
}


void poissrnd()
{
double mu=7;
double p;
p=exp(-mu);
double x=1;
double a=1;
double r;
r=srandd();
mark1:
a = a*r;
if (a<p){
printf(" %f",x);
}
else{
x=x+1;
goto mark1;
}

}

int main(int argc, char *argv[])
{
poissrnd();
system("PAUSE");
return EXIT_SUCCESS;
}

[type_Oleg]

Я сам попробовал. Правда, на Паскале, я Си - не очень знаю
Второй вариант - что-то не то получается.
Я генерировал по сотне значений, и потом делал статистическую оценку параметров. То есть оценка мат.ожидания M(X) - выборочное среднее, оценка дисперсии D(X) - исправленная выборочная дисперсия.
Для распределения Пуассона должно быть: M(X)=D(X)=lambda.

Для второго варианта - не то. M(X) получается в несколько раз больше, а D(X) - вообще на порядки.
А вот первый - нормально. При lambda=7 и выборках в 100 значений, и M(X) и D(X) - 6,2... 7,5 .

function Poison(lam:Double):Integer;
var x:Integer;
  p,r:Double;
begin
 p:=Exp(-lam);
 x:=0;
 r:=Random-p;   // Random возвр. от 0 до 1
 while r>=0 do
  begin
   x:=x+1;
   p:=p*lam/x;
   r:=r-p;
  end;
 Result:=x;
end;
// и где-то там ...
 Randomize;  // это типа srand

PS. К вас, кстати ошибка : double x=1;
Надо double x=0;

[type_Oleg]

Вот, функция на Си, по первому алгоритму. Работает, проверял, подставляя значения вместо closed_interval_rand.

int poison(double lam)
{ double p=exp(-lam);
  double x=0;
  double r=closed_interval_rand(0.0,1.0)-p;
  while (r>=0.0){
     x+=1;
     p=p*lam/x;
     r-=p;      }
  return x;     
    }

Распределение Пуассона - дискретное, целочисленное, поэтому int.

[hp114frost]

Спасибо большое, извиняюсь за оффтоп, но не могу не написать благодарность в явном виде!
Да, пост про мат ожидание заставил меня задуматься, вы абсолютно правы!
Как здорово, что есть еще люди, которые не ленятся давать хорошие ответы и у которых хватает на это знаний!