Rozkład normalny Gaussa - źle losuje

Chcę wylosować liczby korzystając z rozkładu normalnego gaussa. Liczby maja byc z przedzialu [0,x] gdzie x ustalam sobie w programie.

napisałem sobie już kod, ale nie działa poprawnie, tzn. losowane liczby nie sa z przedzialu [0,x] ...

#include <iostream>
 #include <stdlib.h>
 #include <math.h>
#include <time.h>
using namespace std;

#define    PI    3.14159265

 double randn(double mu, const double sigma)
 {
     const double gap = 0.01;
     const int length = (int)(2*sigma*6/gap);
     double *X = new double[length];
     double *P = new double[length];

     X[0] = -sigma*6+mu;
     for (int k=1; k<length; k++)
         X[k] = X[k-1] + gap;

     P[0] = 0.0;

     for (int k=1; k<length;k++)
         P[k] = P[k-1] + exp(-((X[k]-mu)*(X[k]-mu))/(2*sigma*sigma))/(sqrt(2*PI)*sigma);

     for (int k=0; k<length; k++)
         P[k] = P[k]/P[length-1];

     double temp_rand = (double)rand()/RAND_MAX;

     int i;

     for(i=0; i<length;i++)
     {
         if(P[i]>temp_rand)
             break;
     }

     return X[i-1] + (X[i]-X[i-1])*(temp_rand-P[i-1])/(P[i]-P[i-1]);
 }

int main()
 {
     srand((unsigned)time(0));

     for(int i=1; i<100; i++)
     {
         cout<<randn(0, 10) << endl;
     }

     return 0;
 }

Dowiedz się najpierw co to jest rozkład Gaussa. Podpowiem, że nie ma on określonych granic oraz przy takim wywołaniu prawidłowa funkcja będzie z równym prawdopodobieństwem zwracała liczby dodatnie i ujemne oraz 99,7% zwracanych liczb będzie w zakresie (-30,+30).

Na dodatek nie trzeba tego pisać, bo w C++ to już jest: http://www2.research.att.com/~bs/C++0xFAQ.html#std-random

Użyj C++11

std::linear_congruential engine;
std::normal_distribution<double> dist;
auto random = std::bind(dist, engine);

random();

wiem, ze to juz jest napisane, ale ja musze to sam napisac:) to dlaczego nauczyciel kazal nam napisac funkcje, ktora bedzie losowac liczby z przedzialu [0,x] rozkladem normalnym? czy to ja czegos nie rozumiem?

Może funkcję która przyjmuje parametr x będący odchyleniem standardowym? A może nauczyciel nie wie co to rozkład normalny? Oceń sam.

takie coś powinno w przybliżeniu dać rozkład gaussa (liczby całkowite):

int gauss(int a, int b)
{
  int n = b-a+1;
  int r = 0;
  while (n--)
    r += rand()%2;
  return r+a;
}

dla 10^6 losowań w zakresie <0,30> dostałem takie wyniki:

 0: 0
 1: 0
 2: 0
 3: 0
 4: 0
 5: 8
 6: 175
 7: 1155
 8: 3860
 9: 9521
10: 21002
11: 40023
12: 65366
13: 96426
14: 123139
15: 139234
16: 139302
17: 123141
18: 96583
19: 65319
20: 40011
21: 21043
22: 9501
23: 3853
24: 1156
25: 175
26: 7
27: 0
28: 0
29: 0
30: 0

dla double można to zrobić tak:

• załóżmy, że dysponujesz funkcją, która losuje liczby w rozkładzie równomiernym [0,x]
• sumujesz k losowań taką funkcją z zakresu [0,x/k]. nim większe k tym bardziej powinno to przypominać rozkład gaussa. (około 30 powinno starczyć)

Prawdopodobnie chodziło mu (nauczycielowi) o rozkład prostokątny (równomierny). Tutaj przyda się informacja, że rand()/(RAND_MAX+1.0) ma rozkład równomierny w zakresie <0,1). Pamiętaj o dodaniu do RAND_MAX 1.0, a nie 1, bo twój program może być nieprzenośny.

BTW ta implementacja w oryginalnym poście wygląda tragicznie. Dobry algorytm jest w części 3.4.1.C Sztuki programowania Knutha (nie chce mi się przepisywać).

No właśnie on tak dziwnie powiedział, że mają być losowane z okreslonego przedziału rozkładem normalnym.

@Zjarek: nigdzie nie mogę znaleźć książki, mógłbyś może jednak coś podpowiedzieć?:)

EndrjuNa dodatek nie trzeba tego pisać, bo w C++ to już jest: http://www2.research.att.com/~bs/C++0xFAQ.html#std-random</quote napisał(a)

coś niedorobiony i niekompletny ten przykład. przerobiłem trochę:

#include <random>
#include <iostream>
#include <map>
#include <functional>
 
using namespace std;
 
int main()
{
    random_device rd;
    default_random_engine re(rd());
    normal_distribution<double> nd(31,6);
    auto norm = bind(nd,re);
    map<int,int> Gauss;
    for (int i=0; i<1000000; i++)
        Gauss[(int)norm()]++;
    for (int i=10; i<50; i++)
    {
        cout << i <<":\t";
        int starCount = Gauss[i]/3000;
        for (int j=0; j<starCount; j++)
            cout << '*';
        cout<<'\n';
    }
}

wynik:
58550fbe726b26.76870

coś niedorobiony i niekompletny ten przykład. przerobiłem trochę:

#include <random>
#include <iostream>
#include <map>
#include <functional>
 
using namespace std;
 
int main()
{
    random_device rd;
    default_random_engine re(rd());
    normal_distribution<double> nd(31,6);
    auto norm = bind(nd,re);
    map<int,int> Gauss;
    for (int i=0; i<1000000; i++)
        Gauss[(int)norm()]++;
    for (int i=10; i<50; i++)
    {
        cout << i <<":\t";
        int starCount = Gauss[i]/3000;
        for (int j=0; j<starCount; j++)
            cout << '*';
        cout<<'\n';
    }
}

wynik:

[azarien@freebsd9:~/myprogs]$ g++47 random.cpp --std=c++11 -o random
[azarien@freebsd9:~/myprogs]$ ./random
10:        
11:        
12:        
13:        
14:        
15:        
16:        *
17:        *
18:        **
19:        ***
20:        ****
21:        ******
22:        ********
23:        **********
24:        ************
25:        **************
26:        ****************
27:        ******************
28:        ********************
29:        *********************
30:        **********************
31:        *********************
32:        *********************
33:        ********************
34:        ******************
35:        ****************
36:        **************
37:        ************
38:        **********
39:        ********
40:        ******
41:        ****
42:        ***
43:        **
44:        *
45:        *
46:        
47:        
48:        
49:

łał.

ok, działa pięknie. tylko ja dalej nie wiem, jak mam rozumiec to, ze te liczby maja byz z przedzialu? wyrzucac te, ktore nie beda? to taka troche na chama metoda :/

No i ok ... zapytałem się nauczyciela [e-mailem], o co dokładnie mu chodziło, odpisał mi tak: mamy losować liczby z przedzialu jednostajengo ... Co to moze oznaczac?

No to rozkład równomierny, nie Gaussa.

Spoko, tzn, że zwykłe: ((double) rand() / (x+1.0)) starczy?

Tak, starczy. Chyba, że chcesz mieć dużo lepszy szum to wtedy musisz użyć innego PRNG'a. Ale jeśli to tylko zadanie domowe, to absolutnie wystarczy, zwykły człowiek nie rozróżni tego od szumu.

okazało się, że to jest źle - to wcale mi nie losuje liczb z przedziału [0,x] tylko jakieś kosmiczne liczby...

Ogólne rozwiązanie:

double min = 1.0;  // dolna granica przedzialu 
double max = 10.0; // gorna granica przedzialu

// losowa liczba z przedzialu [min, max]
double random_number = (rand() / static_cast<double>(RAND_MAX)) * (max - min) + min;

http://ideone.com/MwUuq

U Ciebie min = 0 i kod się uprości.

wygląda ok, to chyba będzie to, o ile to daje rozkład jednostajny:>