Z Skrypty dla studentów Ekonofizyki UPGOW
Spis treści |
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
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)
