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)

Liczby o zadanym rozkładzie

Rozkład Gaussa

Algorytm Boxa-Mullera

Źródło „http://172.16.3.1/ekonofizyka/index.php/MKZR:Liczby_losowe”