Z Skrypty dla studentów Ekonofizyki UPGOW
Spis treści |
Generacja liczb losowych
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
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.
x=rand(1,3000); X=reshape(x,[3,length(x)/3]); plot3(X(1,:), X(2,:), X(3,:), 'o')
disp('Zera funkcji, metodami dwoma') clear all % na wszelki wypadek, kasujemy wszystkie zmienne w octave % definiujemy funkcje, ktore oblicza kolejne iterajec metod: Newtona i sieczncyh function x=newton(x0,N,f,df) disp('Metoda Newtona') x=zeros(1,N); x(1)=x0; for i=2:N x(i)=x(i-1)-f(x(i-1))/df(x(i-1)); endfor endfunction function x=sieczna(x0,x1,N,f) disp('Metoda siecznych') x=zeros(1,N); x(1)=x0; x(2)=x1; for i=3:N s = (f(x(i-1)) - f(x(i-2)) ) /(x(i-1)-x(i-2)); x(i)=x(i-1)-f(x(i-1))/s; endfor endfunction % to jest funkcja pomocnicza, wykonujaca wykres function plotall(newt,siecz,f,a,b) subplot(2,1,1) plot(newt,"r*-") hold on plot(siecz,"b*-") hold off subplot(2,1,2) x=linspace(a,b,101); plot(x,f(x)) hold on plot(newt,f(newt),"r*") plot(siecz,f(siecz),"g*") hold off endfunction % teraz mozemy zdefiniowac funkcje matemetyczne, np f,g % i oblicznyc ich zera f=@(x) sign(x).*sqrt(abs(x)); df=@(x) 1/(2*sqrt(abs(x))); N=10 newt=newton(1.1,N,f,df); siecz=sieczna(1.5,0.4,N,f,df); figure(1) % figure pozwoli otworzyc kilka okien z wykresami plotall(newt,siecz,f,-2.1,2) g=@(x) sin(x)./x; dg=@(x) cos(x)./x-sin(x)./(x.^2); N=10 newt=newton(1.1,N,g,dg); siecz=sieczna(1.5,0.4,N,g); figure(2) plotall(newt,siecz,g,0.01,12) h=@(x) atan(x); dh=@(x) 1.0./(1+x.^2); N=5 newt=newton(1.1,N,h,dh); siecz=sieczna(1.5,0.4,N,h); figure(3) plotall(newt,siecz,h,-3.001,4)