Z Skrypty dla studentów Ekonofizyki UPGOW
(→Liczby o zadanym rozkładzie) |
(→Liczby o zadanym rozkładzie) |
||
Linia 110: | Linia 110: | ||
X=rand(1,N); | X=rand(1,N); | ||
</source> | </source> | ||
- | jest generowane N liczb z rozkładem jednostajnym (wykorzystującym wbudowany w system generator liczb losowych, w przypadku GNU/Octave jest do Mersenne Twister). | + | jest generowane N liczb z rozkładem jednostajnym (wykorzystującym wbudowany w system generator liczb losowych, w przypadku GNU/Octave jest do [http://pl.wikipedia.org/wiki/Mersenne_Twister Mersenne Twister]). |
Wersja z 08:13, 24 mar 2010
Spis treści |
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.
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.
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)\). Poniższy skrypt jest przykladem wykorzystania transformacji zmiennej losowej do otrzymania generatora o rozkładzie eksponencjalnym z jednostajnego.
N=10000 X=rand(1,N); X=-log(X); hold off; h=0.1 xmax=16 hist(X,[-1:h:xmax],1/h) # drugi argument to norma histogramu hold on; fplot(@(x) exppdf(x,1),[1e-3,xmax],'r')
W lini
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).
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])