KADD 2022 Laboratorium 6 EN

From Łukasz Graczykowski

(Difference between revisions)

Jump to: navigation, search

Revision as of 11:34, 4 April 2022

Exercise

Part one: linear congruent generator of pseudorandom numbers (1 pkt.)

Please write a generator of pseudorandom numbers and save generated numbers into a file.

The generator should be based on the following formula:

x[j+1] = (g*x[j] + c) mod m.

This generator is called LCG - linear congruent generator. Providing the first (beginning) number x[0] defines the whole series, which is periodic. The period depends on the parameters and under certain conditions it reaches maximum value of m. These conditions are:

c and m do not have joint divisors,
b = g-1 is a multiply of every prime number p, which is a divisor of a number m,
b is a multiply of 4 if n is also a multiply of 4.

For simplicity we can uuse c = 0, and in that case we get a multiplicative generator (MLCG).

Value of g and m should be easy to modify in the program.

The result of the macro execution should be a file with the name name.dat containing a series of generated numbers for a given set of parameters. The macro should be executed three times, resulting in three files: random1.dat, random2.dat, random3.dat, for the following parameters, respectively:

m=97 i g=23,
m=32363 i g=157,
m=147483647 i g=16807.

Second part: spectral test (1 pkt.)

Please perform a spectral test to see what is the quality of the generator. In order to do so please draw plot on the (x[n], x[n+1]) points from generated numbers. The obtained graph will show a spectral pattern of the generator - hence the name of the test.

If the points are distributed uniformly, the test can be judged good. If they show a certain periodic pattern - the test doesn't work properly. Of course, for the position of the points what maters are the parameters g and m.

In order to draw spectral test, please use TH2D.

As a result, we should have three spectral tests.

Third part: generation of pseudorandom numbers based on transformation of a uniform distribution (3 pkt.)

Each function of the random variable is also a random variable. What is the probability distribution of a random variable Y, if the probability distribution of X, f(x), is known. We assume, that the probability g(y)dy is equal to f(x)dx, where dx equals dy. The condition is of course fulfilled for (infinitely) small dx. What results from this is the following:

g(y) = dx/dy f(x)

Teraz jeżeli założymy, że gęstość prawdopodobieństwa f(x) wynosi 1 w 0<=x<=1 i f(x) = 0 dla x<= 0 i x>1 to powyższe równanie możemy zapisać w postaci:

g(y)dy = dx = dG(y),

gdzie G(y) jest dystrybuantą zmiennej losowej Y. Co po całkowaniu daje nam

x = G(y) => y = G^-1(x).

Jeśli zmienna losowa X ma rozkład jednostajny na odcinku pomiędzy 0 i 1 oraz jeśli znana jest funkcja odwrotna G^-1(x) to funkcja g(y) opisuje gęstość prawdopodobieństwa rozkładu zmiennej losowej Y.

Używając tej metody należy wygenerować 10000 liczb z rozkładu:

Dla tau = 2:

Należy wygenerować 10000 liczb z rozkładu 0 do 1 używając generatora z części pierwszej (zapisać do pliku wartości xn makrem z pierwszej części a następnie je wczytać w makrze z części drugiej).
Analitycznie (na kartce) policzyć dystrybuantę tego rozkładu, a następnie funkcję odwrotną. (1 pkt.)
Wygenerować rozkład f(x) - wrzucając wygenerowane wartości do histogramu - korzystając z: (1 pkt.)
- liczb wygenerowanych wcześniej i wczytanych z plików losowe1.dat, losowe2.dat, losowe3.dat,
- standardowego generatora ROOT'a, np. gRandom->Uniform(1) (obiekt gRandom istnieje domyślnie w uruchomionej instancji ROOT; można oczywiście również stworzyć samodzielnie obiekt TRandom - link).
Narysować na jednym wykresie histogram (odpowiednio unormowany) oraz funkcję teoretyczną f(x) (obiekt TF1). (1 pkt.)

Uwagi

Czytamy dokładnie Wykład 4 (link), zwłaszcza slajdy 6 oraz 18-25
Na samym początku, przed losowaniem, musimy samodzielnie ustawić wartość pierwszej liczby pseudolosowej x0 (tzw. ziarno, "seed"). Jeżeli chcemy, by za każdym razem liczby pseudolosowe były inne, możemy je ustawić z zegara systemowego:

  x0 = time(NULL);

Parametry histogramów z obrazków poniżej:

 TH1D *hUniform = new TH1D("hUniform","Uniform distribution",100,0,1);
 TH2D *hCorr = new TH2D("hCorr","Correlation",100,xmin,xmax,100,0,1);

Ilość losowań w części pierwszej:

 const int N = 1000000;

Wczytywanie danych z pliku:

ifstream ifile;
ifile.open("dane.dat");
double val;
while(ifile>>val)
{
  cout<<val<<endl;
}
ifile.close();

Zapisywanie danych do pliku:

ofstream ofile;
ofile.open("dane.dat");
for(int i=0;i<N;i++)
  ofile<<val<<endl;
}
ofile.close();

Wynik

Przykładowy rozkład dla parametrów:

m=97, g=23

m=2147483647, g=16807

Przykładowy wynik transformacji rozkładu jednorodnego:

Navigation