MKZR:Liczby losowe
Z Skrypty dla studentów Ekonofizyki UPGOW
(Różnice między wersjami)
(→Generacja liczb losowych) |
|||
Linia 1: | Linia 1: | ||
==Generacja liczb losowych== | ==Generacja liczb losowych== | ||
+ | |||
+ | |||
+ | ==Generacja liczb losowych== | ||
+ | |||
+ | 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. | ||
+ | |||
+ | Najprostrzym przykładem jest liniowy generator kongruencyjny (ang. Linear congruential generator, LCG), który jest wyznaczony przez relację rekurencyjną: | ||
+ | |||
+ | : <math>X_{n+1} = \left( a X_n + c \right)~~\bmod~~m</math> | ||
+ | |||
+ | Jego implementacja w octave: | ||
+ | |||
+ | <source lang="matlab"> | ||
+ | function y=myran(x); | ||
+ | a=1664525; | ||
+ | b=1013904223; | ||
+ | m=2^32; | ||
+ | y=mod(a*x+b,m); | ||
+ | return; | ||
+ | end | ||
+ | </source> | ||
+ | |||
+ | przykład jego użycia do generacji pseudolosowych liczb z przedziału (0,1): | ||
+ | <source lang="matlab"> | ||
+ | x(1)=123; | ||
+ | for i=2:10; | ||
+ | x(i)=myran(x(i-1)); | ||
+ | disp(x(i)/2^32); | ||
+ | end | ||
+ | </source> | ||
+ | |||
+ | |||
+ | Generator ten ma wiele mankamentów: | ||
+ | |||
+ | * okres niskich bitów jest o wiele niższy od okresu całego generatora. | ||
+ | |||
+ | <source lang="matlab"> | ||
+ | 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 | ||
+ | </source> | ||
+ | |||
+ | * Jeżeli użyjemy LCG do generacji punktów w n-wymiarowej przestrzeni to punkty te będą leżały na <math>m^{1/n}</math> hiperpowierzchniach. | ||
+ | |||
+ | <source lang="matlab"> | ||
+ | function y=myranbad(x); | ||
+ | a=65539; | ||
+ | b=0; | ||
+ | m=2^31; | ||
+ | y=mod(a*x+b,m); | ||
+ | return; | ||
+ | end | ||
+ | </source> | ||
+ | |||
+ | <source lang="matlab"> | ||
+ | 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,'*') | ||
+ | </source> | ||
+ | |||
+ | |||
+ | "Jedynym słusznym" generatorem używanym obecnie jest [http://pl.wikipedia.org/wiki/Mersenne_Twister Mersenne Twister]. Jest on szybki i zapewnia dobre własności statystyczne generowanego ciągu liczb. W systemie GNU Octave funkcja rand używa | ||
+ | Mersenne Twister-a. | ||
+ | |||
+ | |||
<source lang="matlab"> | <source lang="matlab"> |
Wersja z 20:43, 15 mar 2010
Spis treści |
Generacja liczb losowych
Generacja liczb losowych
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.
Najprostrzym 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.
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,'*')
"Jedynym słusznym" generatorem używanym obecnie jest Mersenne Twister. Jest on szybki i zapewnia dobre własności statystyczne generowanego ciągu liczb. W systemie GNU Octave funkcja rand używa Mersenne Twister-a.
function y=myran(x); a=1664525; b=1013904223; m=2^32; y=mod(a*x+b,m); return; end
x(1)=123; for i=2:10; x(i)=myran(x(i-1)); disp(x(i)/2^32); end
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)