Z Skrypty dla studentów Ekonofizyki UPGOW

Aktualna wersja na dzień 15:20, 23 sty 2011

Wektoryzacja obliczeń w środowsku matlab/GNU Octave

Matlab/GNU Octave są językami interpretowanymi. W znacznym stopniu ułatwia to przejście do sedna rozwiązywanego problemu z pominięciem wielu informatycznych aspektów programowania. Jednak ceną za taki komfort jest często wydajność pisanych programów. W przypadku ogólnym bezpośrednie przepisanie kodu w języku typu Fortran czy C/C++ na matlab może spowolnić jego działanie kilkaset razy. Stosując jednak języki interpretowane świadomie i wykorzystując ich możliwości w wielu przypadkach można uniknąć spowolnienia i otrzymać bardzo wydajny program.

Podstawowym elementem języka wykorzystywanym w jest pętla "for". Rozważmy następujący program:

N=15000;
x=linspace(0,2*pi,N);
tic;
for i=2:N
   x(i)=x(i)*0.1;
endfor
toc

Inicjalizujemy wektor x zawierający punkty jednorodnie położone na odcinku $(0,2 \pi)$. Chcemy pomnożyć każdą ze współrzędnych przez 0.1. Powyższy program wykonuje to za pomocą pętli for. Jest to poprawny sposób wykonania tego zadania. Instrukcje tic/toc mierzą czas wykonania wszystkich komend pomiędzy nimi. Na typowym komputerze PC wynosi on ok 0.2 sekundy.

Wykonajmy identyczną operację w sposób który umożliwia nam śwodowisko matlab/GNU Octave. Operator mnożenia jest zdefiniowany nie tylko na liczbach, ale również na wektorach. Możeme więc zapisać:

N=15000;
x=linspace(0,2*pi,N);
tic
x=0.1*x;
toc

Przy takim sposobie zapisu czas wykonania jest rzędu 0.001 sekundy.

Rozważmy teraz trochę bardziej skomplikowany problem. Majać wektor z wartosciami pewnej funkcji na odcinku chcemy policzyć pochodna korzystając ze wzroru na iloraz różnicowy. Kod napisany z użyciem pętli for może wyglądać tak:

N=15000;
x=linspace(0,2*pi,N);
z=sin(x);
h=2*pi/N;
 
tic
for i=2:N
   diffz(i)=(z(i)-z(i-1))/h;
endfor
toc
 
diffz(1)=diffz(2);
plot(x,z,x,diffz)

a czas jego wykonania będzie również rzędy 0.2 sekundy. By napisać program wektorowy warto zapoznać sie z funkcja przesuwającą index w wektorze "shift".

octave:142> shift([1,2,3],1)
ans =
 
   3   1   2

Dysponując taką funkcję można wektorowo zapisać żądaną operacje w poniższy sposób:

x=linspace(0,2*pi,N);
z=sin(x);
h=2*pi/N;
tic
for i=2:N
   diffz(i)=(z(i)-z(i-1))/h;
endfor
 
toc
diffz(1)=diffz(2);
tic 
 
diff2z=(z-shift(z,1))/h;
toc
plot(x,z,x,diff2z)

Czas wykonania takiego programu jest znacząca krótszy i wynosi ok 0.01 sekundy.

Histogramy

Ważną umiejętnością jest estymacja funkcji gęstości z pewnej (dużej) liczby danych. Proces taki jest nazywany histogramowaniem.

Histogram jest to sposób przedstawienia częstości występowania pewnej cechy w postaci wykresu słupkowego. Weźmy na przykład ciągłą zmienną losową. Podzielmy przestrzeń wartości przyjmowanych przez tą zmienna na N równych przedziałów. Przypuścmy, że dysponumemy M realizacjami zmiennej losowej. Na histogramie można umieścić liczbę realizacji, które należa do danego przedziału. Co więcej taki wykres bedzie w granicy $N,M\to \infty$ taki sam jak wykres gęstości prawdopodobieństwa nasze zmiennej losowej. By otrzymać dokładnie gęstość prawdopodobieństwa trzeba przeskalować liczbę zliczeń do gęstości. Niech $n_i, i=1..N$ będą zliczeniami w kolejnych przedziałach. Gęstość prawdopodobieństwa jest unormowana do jedności (bierzemy zmienną losową przyjmującą wartości między a i b):

$1=\int_a^b f(x) dx=\sum_{i=1}^N f(x_i) h$

gdzie h jest szerokością jednego przedziału $h=(x_i-x_{i-1})/(N-1)$. Wynika z tego, że

$\sum_{i=1}^N f(x_i) =1/h$,

podczas gdy liczby zliczeń są unormowane do pewnej wartości $\sum_{i=1}^N n_i$.

Widać z tego, że formuła łącząca liczby zliczeń z wartościami gęstości wewnątrz odcinków wyraża się:

$f_i=f(x_i)=\frac{1}{h}\frac{n_i}{\sum_{i=1}^N n_i}$.

W systemie GNU Octave wbudowana funkcja hist umożliwia podanie drugiego argumentu, który jest czynnikiem normalizacyjnym sumy wszystkich zliczeń i w rozważanym przypadku powinien od wynosić 1/h.

Histogramy 1D

Jako przykład wykonajmy histogram 10000 liczb o rozkładzie normalnym i porównajmy go ze wzrorem analitycznym na funkcję Gaussa. Najprostszy skrypt w GNU Octave realizujący to zadanie to:

xmax=5;
h=0.1;
mydata=normrnd (0,1,10000,1);
hist(mydata,[-xmax:h:xmax],1/h);

Zauważmy, że:

Nie została zadana liczba przedziałow, ale szerokość. Liczba przedziałow w tym przypadku to:

octave:38> length( [xmin:h:xmax] )
ans =  101

Pomimo, ze zmienna losowa przyjmuje wartości na całej osi rzeczywistej, to histogram wykonaliśmy na skończonym przedziale. Wybranie wartości x a przedziału -5..5 jest w przypadku rozkładu Gaussa dobrym wyborem. Zauważmy jest stosunek funkcji gęstości w punktach 0 i 5 wynosi

octave:42> normpdf(0,0,1)/normpdf(5,0,1)
ans =  2.6834e+05

więc dla 10000 punktów możemy się spodziewać, że zaden z nim nie wyjdzie poza przedział -5...5.

Proszę sprawdzić jak będzie wyglądał histogram jeśli weżmiemy xmin=-1 i xmax=1

Funkcja hist została wywołana bez argumentów zwrotnych i zgodnie z dokumentacja narysowała ona wykres. Można postąpić inaczej:

xmax=5;
h=0.1;
mydata=normrnd (0,1,10000,1);
[NN,XX]=hist(mydata,[-xmax:h:xmax],1/h);

Podając argumenty zwrotne funkcja hist zamiast rysować wykres, zwraca wektor liczb zliczeń i wartość zmiennej w przedziałach. Można to wykorzystać:

plot(XX,NN,"-",XX,normpdf(XX,0,1),"r-")
bar(XX,NN)

Jak działa funkcja hist?

Poniższy kod ilustruje prostą implementacje funkcji histogramującej w GNU/Octave. Nie jest ona efektywna, gdyż wykorzystuje pętlę po wszystkich danych wejściowych (kod nie jest wektorowy).

xmax=5;
xmin=-5;
h=.1;
XX=[xmin:h:xmax];
mydata=normrnd (0,1,10000,1);
NN=zeros(1,101);
for i=1:length(mydata)
	idx=  ceil ( (mydata(i)-xmin)/(xmax-xmin) * length(NN)  ) ;
	if idx>0 && idx<=length(NN)
 		NN(idx)++;
	endif
endfor
plot(XX,NN)

Histogramy 2D

Estymacja gęstości dla dwów nielależnych zmiennych losowych Gaussowskich za pomocą funkcji hist2d.

Podobnie jak w przypadku jednej zmiennej losowej, histogram można utworzyć dla wielu zmiennych losowych. Do celów wizualizacyjnych nadaje się jednak tylko przypadek dwuwymiarowy. W sysemie GNU/Octave nie ma bezpośrednio wbudowanej funkcji histogramującej dwie zmienne losowe, ale jest pakiet rozszerzający plot (dostępny z [1]). Po zciągnieciu pliku instalujemy go za pomocą

pkg install 'plot-1.0.8.tar.gz'

Z jego pomoca narysujemy histogram dwoch niezależnych gausowskich zmienych losowych:

hist2d(normrnd(0,1,250000,2),60);

Podobnie jak w przypadku funkcji hist można dane zapisać do tablic by np. użyć innej procedury rysującej:

[nn,x,y]=hist2d(normrnd(0,1,134100,2),60);
imagesc(nn)

Estymacja gęstości dla dwów nielależnych zmiennych losowych Gaussowskich za pomocą funkcji hist2d, dane przedstawione za pomocą imagesc.