KADD 2022 Laboratorium 6 EN

From Łukasz Graczykowski

(Difference between revisions)

Jump to: navigation, search

Revision as of 11:31, 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.

Część trzecia: generacja liczb losowych oparta na transformacji rozkładu jednorodnego (3 pkt.)

Dowolna funkcja zmiennej losowej jest zmienną losową. Powstaje więc pytanie jaka jest gęstość zmiennej losowej Y jeżeli znana jest gęstość f(x). Zakładamy, że prawdopodobieństwo g(y)dy jest równe f(x)dx, gdzie dx odpowiada wartością dy. Warunek jest spełniony dla odpowiednio małych dx. Wynika stąd, że:

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