Z Skrypty dla studentów Ekonofizyki UPGOW

Generacja liczb losowych: liczby o rozkładzie jednostajnym

Liniowy generator kongruencyjny

Generator liczb pseudolosowych to procedura, generująca deterministycznie ciąg bitów, który pod pewnymi względami jest nieodróżnialny od ciągu uzyskanego z prawdziwie losowego źródła.

Najprostszym przykładem jest liniowy generator kongruencyjny (ang. Linear congruential generator, LCG), który jest wyznaczony przez relację rekurencyjną:

\(X_{n+1} = \left( a X_n + c \right)~~\bmod~~m\)

Jego implementacja w octave:

function y=myran(x); 
  a=1664525;
  b=1013904223;
  m=2^32;
  y=mod(a*x+b,m);
  return; 
end

przykład jego użycia do generacji pseudolosowych liczb z przedziału (0,1):

x(1)=123;
for i=2:10; 
	x(i)=myran(x(i-1));
	disp(x(i)/2^32);
end

Generator ten ma wiele mankamentów:

• okres niskich bitów jest o wiele niższy od okresu całego generatora.

x(1)=1234;
for bitn=1:10
	for i=2:30; 
		x(i)=myran(x(i-1)); 
   	printf("Bit %d: %d\n",	 bitn,bitget(x(i),bitn) ); 
	end
	a=input('cont?','s');
	if (a=='q') 
		break;
	end
end

• Jeżeli użyjemy LCG do generacji punktów w n-wymiarowej przestrzeni to punkty te będą leżały na \(m^{1/n}\) hiperpowierzchniach. W przypadku idealnego generatora jednostajnego, rozkład punktów jest izotropowy i nie powinien zawierać żadnych regularności.
  

  

  Hiperpłaszczyzny w przestrzeni trzech kolejnych wywołań generatora LCG

function y=myranbad(x); 
  a=65539;
  b=0;
  m=2^31;
  y=mod(a*x+b,m);
  return; 
end

x(1)=1234;
for i=2:3000; 
	x(i)=myranbad(x(i-1)); 
end
X=reshape(x,[3,length(x)/3]);
randmax=2^31*1.0
plot3(X(1,:)/randmax,
  X(2,:)/randmax,
  X(3,:)/randmax,
  '*')

Generatory wysokiej jakości Mersenne Twister

Jednym z lepszych generatorów używanym obecnie jest Mersenne Twister. Jest on szybki i zapewnia dobre własności statystyczne generowanego ciągu liczb. Jego wadą jest stosunkowo duża liczba instrukcji z których się składa, co ma znaczenie w przypadku jego implementacji na architekturach wbudowanych a nie odgrywa większej roli na klasychnych komputerach PC. W systemie GNU Octave funkcja rand używa właśnie generatora Mersenne Twister.

Przestrzeń trzech kolejnych wywołań generatora Mersenne Twister, nie jest widoczna struktura hiperpłaszczyzn.

Poniższy kod przedstawia graficznie kolejne trzy liczby wylosowane za pomocą generatora Mersenne Twister. Wykres trzech kolejnych wygenerowanych liczb wygląda bardziej jednorodnie od poprzedniego.

x=rand(1,3000);
X=reshape(x,[3,length(x)/3]);
plot3(X(1,:),
      X(2,:),
      X(3,:),
      'o')

Liczby o zadanym rozkładzie

Dysponując generatorem liczb pseudolosowych o rozkładzie jednostajnym z wartościami w (0,1) możemy otrzymać generator o dowolnym rozkładzie dokonując z transformacji gęstości prawdopodobieństwa. Jeśli \(u\) jest zmienną o rozkładzie jednostajmym o docelowy rozkład ma dystrybuantę \(F_{\xi}(x)=\int_{-\infty}^x f(x)\) to zmienna y: \(y=F_{\xi}^{-1}(u)\) będzię miała rozkład \(\displaystyle f(y)\). Metodę tą można wykorzystać do otrzymania generatora o rozkładzie eksponencjalnym \(e^{-x}\) dysponując rozkładem jednostajnym. W takim przypadku dystrybuanta jest również funkcją \(F(x)=e^{-x}\) i jej funkcja odwrotną jest \(F^{-1}(x)=-log(x)\). Możemy to sprawdzić wykonując eksperyment numeryczny polegający na wysymulowaniu dużej liczby próbek z generatora jednostajnego i następnie zrobieniu histogramu tych wielkości.

Histogram 10000 liczb losowych po transformacji y=-ln(x). Czerwona linia to funkcja gęstości rozkładu eksponencjalnego \(e^{-x}\)

N=10000;
X=rand(1,N);
X=-log(X);
hold off;
h=0.1
xmax=16
hist(X,[-1:h:xmax],1/h) 
hold on;
fplot(@(x) exppdf(x,1),[1e-3,xmax],'r')

• w linii

X=rand(1,N);

jest generowane N liczb z rozkładem jednostajnym (wykorzystującym wbudowany w system generator liczb losowych, w przypadku GNU/Octave jest do Mersenne Twister).

• linia:

hist(X,[-1:h:xmax],1/h)

generuje histogram z danych w tabeli X. Proszę zauważyć ze drugi argument jest normą tego histogramu, która zgodnie z dokumentacja (help hist) jest sumą wartości wszystkich słupków. Ponieważ chcemy porównać ten histogram z gęstością to mamy:

\(1=\int_0^\infty f(x) dx=\sum_{i=1}^N f(x_i) h\) z czego nam wynika, że suma wysokości słupków gęstości unormowanej do jedynki wynosi:

\(\sum_{i=1}^N f(x_i) =1/h\)

• w linii

fplot(@(x) exppdf(x,1),[1e-3,xmax],'r')

wykres rysujemy od 0.001 zamiast 0.0. Jest to spowodowane tym, że funkcja exppdf nie jest ciągła w zerze i procedura rysująca fplot będzię usiłowała zagęścić niepotrzebnie punkty wykresu w nieciągłości.

Rozkład Gaussa

N=10000
n=3;
for i=1:N;
	gclt(i)=(sum(rand(1,n))-n*0.5)/(0.3*sqrt(n));
end;
hold off;
hist(gclt,[-5:0.1:5],10)
hold off;
hist(gclt,[-5:0.1:5],10)
hold on;
fplot(@(x) normpdf(x),[-5,5])

Algorytm Boxa-Mullera