r_lab_zad3

2693 days ago by aniag

%latex \section{Zadanie 3} W zadaniu 2-gim przećwiczyliście prosty "rejection sampling": punkt był akceptowany, gdy leżał w określonym zbiorze (wewnątrz okręgu), odrzucany wpp. Ogólniejsza metoda "rejection sampling" daje większe możliwości, a w skrócie polega na losowaniu pod wykresem gęstości rozkładu, z którego losować już umiemy. Mianowicie, jeśli \begin{enumerate} \item potrafimy losować z rozkładu o gęstości $g(x)$ \item istnieje $M>0$ t.że $\forall{x}\ Mg(x)>=f(x)$ \end{enumerate} wtedy losowanie z rozkładu z gęstością $g(x)$ punktów $x_1,\ldots,x_n$, następnie akceptacja każdego z nich \\ z prawdopodobieństwem sukcesu $$p_{\textrm{accept}(x_i)} = \frac{f(x_i)}{M \cdot g(x_i)}$$ (odrzucenie wpp.) produkuje próbkę z rozkładu z gęstością $f(x)$. 
       
%html <big><a href="http://www.quantum-physics.polytechnique.fr/physix/wiki/index.php/Rejection_sampling">Tutaj</a> można znaleźć wyjaśnienie z rysunkiem metody rejection sampling.</big> 
       
Tutaj można znaleźć wyjaśnienie z rysunkiem metody rejection sampling.
%latex cd. zadania 3 Napisz w {\bf R} funkcję {\bf rundernorm{\em (n, f, M)}}, która wylosuje próbkę wielkości {\em n} z rozkładu o gęstości zadanej zdefiniowaną poniżej funkcją {\bf f}, przy założeniu, że $\forall{x}\ f(x)\le M\cdot\textrm{Norm}(0,1)$. \noindent Przy rozwiązywaniu wzoruj się na rozwiązaniu Zadania 2. \noindent Do testowania użyj funkcji {\bf f} zdefiniowanej poniżej w {\bf R}: 
       
%r #Definicja f(x) poprzez definicje klasy gestosci f(a)(x) #Nie przejmuj sie skomplikowaniem, zobacz wykres ponizej #Zastosowanie ifelse pozwala ewaluowac funkcje na wektorach wartosci f_klasa <- function(a=1) function(x) ifelse(x < -a,0,ifelse(x>a,0,a-abs(x)))/a^2 #/a^2 gwarantuje całkowalnosc do 1 #Jako przykladowej uzyjemy f(a=2)(x) f <- f_klasa(2) 
       
%r #Funkcja f: curve(f,-3,3) #Jedna z metod na rysowanie funkcji png_output=dev.off() 
       
%r #Plot razem z krzywą Gaussa x<-seq(-5,5,0.1) #Generuje siatke punktow od -5 do 5 z odleglosciami 0.1 png(width=1000,height=500); par(mfcol=c(1,2)) #Uklad dwoch wykresow 1x2 plot(x,f(x),pch="+",type="b"); lines(x,dnorm(x)) #"lines" dorysowuje linie nowego wykresu do istniejacego wykresu plot(x,f(x),pch="+",type="b"); lines(x,2*dnorm(x)) png_output=dev.off() #Widzimy, ze M=2 jest wystarczajace dla naszego f 
       
%r #...jeszcze tylko sprawdzenie, ze f jest gestoscia, przez proste calkowanie: d<-0.001 x<-seq(-5,5,d) sum(f(x)*d) 
       
[1] 1
%r #Rozwiazanie rundernorm <- function(n, f, M){ #TODO } 
       
%r #Przykladowy prosty test poprawnosci v<-rundernorm(100000,f,2) length(v) summary(v) hist(v,breaks=100) png_output=dev.off() 
       
[1] 100000
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
-1.996000 -0.582800  0.006660  0.003689  0.594400  1.995000 
#Jakie inne testy proponujesz? Jaka wlasnosc powinno miec f(x)+f(x-a) w przedziale [0,a]?