MKZR:Liczby losowe

Z Skrypty dla studentów Ekonofizyki UPGOW

(Różnice między wersjami)
(Generatory wysokiej jakości Mersenne Twister)
(Liczby o zadanym rozkładzie)
Linia 89: Linia 89:
==Liczby o zadanym rozkładzie==
==Liczby o zadanym rozkładzie==
 +
 +
Transformacja gęstości prawdopodobieństwa.
===Rozkład Gaussa===
===Rozkład Gaussa===
 +
 +
 +
 +
 +
<source lang="matlab">
 +
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])
 +
</source>
 +
 +
 +
====Algorytm Boxa-Mullera====
====Algorytm Boxa-Mullera====

Wersja z 09:13, 22 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
    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.

Przestrzeń trzech kolejnych wywołań generatora Mersenne Twister, nie jest widoczna struktura hiperpłaszczyzn.
x=rand(1,3000);
X=reshape(x,[3,length(x)/3]);
plot3(X(1,:),
      X(2,:),
      X(3,:),
      'o')

Liczby o zadanym rozkładzie

Transformacja gęstości prawdopodobieństwa.

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])



Algorytm Boxa-Mullera