Z Skrypty dla studentów Ekonofizyki UPGOW
(→Rozkład Gaussa) |
|||
Linia 1: | Linia 1: | ||
- | |||
==Generacja liczb losowych: liczby o rozkładzie jednostajnym == | ==Generacja liczb losowych: liczby o rozkładzie jednostajnym == | ||
Linia 129: | Linia 128: | ||
</source> | </source> | ||
wykres rysujemy od 0.001(1e-3) zamiast 0.0. Jest to spowodowane tym, że funkcja exppdf nie jest ciągła w zerze i procedura rysująca fplot będzie usiłowała niepotrzebnie zagęścić punkty wykresu w nieciągłości w przypadku wpisania 0.0 w dolny zakres. | wykres rysujemy od 0.001(1e-3) zamiast 0.0. Jest to spowodowane tym, że funkcja exppdf nie jest ciągła w zerze i procedura rysująca fplot będzie usiłowała niepotrzebnie zagęścić punkty wykresu w nieciągłości w przypadku wpisania 0.0 w dolny zakres. | ||
+ | |||
+ | ===Transformacja zmiennej losowej=== | ||
+ | |||
+ | |||
===Rozkład Gaussa=== | ===Rozkład Gaussa=== |
Wersja z 10:28, 8 kwi 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)\). 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.
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(1e-3) zamiast 0.0. Jest to spowodowane tym, że funkcja exppdf nie jest ciągła w zerze i procedura rysująca fplot będzie usiłowała niepotrzebnie zagęścić punkty wykresu w nieciągłości w przypadku wpisania 0.0 w dolny zakres.
Transformacja zmiennej losowej
Rozkład Gaussa
Ważnym rozkładem, który niezwykle często się pojawia w różnego rodzaju symulacjach numerycznych jest rozkład Gaussa. Zastosowanie metody transformacji to jego otrzymywania wymaga znajomości odwrot
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])