Najwa ˙zniejsze wnioski z wykładu:
u zastosowania Mocnego Prawa Wielkich Liczb.
u Centralne Twierdzenie Graniczne.
u Definicje estymatorów i ich najwazniejsze cechy.
Zadania do wykonania:
u W Polsce jest n = 24.6·106podatników6i ka ˙zdy z nich mo ˙ze pomyli´c si ˛e przy wypełnianiu zeznania podatkowego, przy czym warto´s´c pomyłki opisuje zmienna losowa Xi, dla którejEXi = 0 oraz σ = pD2Xi = 102PLN. Jaka jest szansa, ˙ze straty pa ´nstwa w wyniku tych bł ˛edów przekrocz ˛a jeden grosz na podatnika? A 5 groszy?
u Liczba podatników:
n=24, 6∗106
√n=4960 E(X) =0
Odchylenie standardowe:
σ= q
D2Xi =100PLN
Korzystamy z Centralnego Twierdzenia Granicznego ( n jest bar-dzo du ˙z ˛a liczb ˛a):
p(X1+X2+X3+...+Xn >xn) =p(X1+X2+...+Xn
σ
√n > xn σ
√n) =1−φ(x
√n σ ) Gdzie φ(t)to warto´s´c dystrybuanty rozkładu normalnego w
punk-cie t.
Za x podstawiamy 0,01 bo szukamy jaka jest szansa, ˙ze pa ´nstwo straci 0,01 PLN na jednego podatnika. Prawdopodobie ´nstwo to wynosi:
p=1−φ(0, 01∗4960
100 ) =1−φ(0, 496) =1−0, 69=0, 31=31%
Dla x=0,05
p=1−φ(0, 05∗4960
100 ) =1−φ(2, 48) =1−0, 9934=0, 0066=0, 66%
n o tat k i z ´c w i c z e ´n z p ro b a b i l i s t y k i 1 9/20 31
u Metoda numeryczna w j ˛ezyku java.
import j a v a . u t i l . Random ; public c l a s s BledyPodatkowe {
s t a t i c double std ev = 1 0 0 ; s t a t i c double n = 2 4 . 6 * 1 0 0 0 0 0 0 ;
public s t a t i c void main ( S t r i n g [ ] a r g s ) { Random r=new Random ( ) ;
double bladNaPodatnika ; double sum ;
i n t k = 0 ; i n t m= 0 ; i n t l i c z b a P r o b = 1 0 0 ;
f o r( i n t j = 1 ; j < l i c z b a P r o b ; j ++) { sum= 0 ;
f o r( i n t i = 0 ; i <n ; i ++) {
sum+=r . nextGaussian ( ) * stdev ; }
bladNaPodatnika =sum/n ; i f( bladNaPodatnika > 0 . 0 1 ) { k + + ; } i f( bladNaPodatnika > 0 . 0 5 ) {m+ + ; } }
System . out . p r i n t l n ( ( double ) k/ l i c z b a P r o b ) ; System . out . p r i n t l n ( ( double )m/ l i c z b a P r o b ) ; }
}
Dla 100 iteracji prawdopodobie ´nstwo strat pa ´nstwa wi ˛ekszych ni ˙z 1 grosz na podatnika wynosi 0,31, a wi ˛ekszych ni ˙z 5 groszy wynosi 0,01. Zrobienie symulacji dla wi ˛ekszych n wymagałoby wi ˛ekszej mocy obliczeniowej.
u W ubiegłym roku studenci Probabilistyki wramach prac domo-wych rzucili kostkami k6 ok. n =104razy. Ile wynosi prawdopo-dobie ´nstwo, ˙ze suma wyrzuconych oczek mie´sci si ˛e mi ˛edzy 34950 a 35050?
u Skorzystano z Centralnego Twierdzenia Granicznego:
∑ni=1xi−nµ σ
√n (28)
Xi - wynik rzutu kostk ˛a
n=104=10000 (29)
EXi=3, 5 (30)
D2= 91
6 − (3, 5)2=2, 9 (31) Poniewa ˙z:
1
6(12+22+32+42+52+62) = 1
6)1+4+9+16+25+36) = 91 (32)6
EX2=
∑
6 i=1k2P(X=k) =
∑
6 i=1k21 6 = 7
2 (33)
P 34950−n·3, 5
√2, 9·n < ∑
6i=1−n·3, 5
√2, 9·n < 35050−n·3, 5
√2, 9·n
!
= (34)
=P(−0, 29< N(0, 1) <0, 29) =F(0, 29) −F(−0, 29) = (35)
=2·0, 61−1=1, 22−1=0, 22 (36) Prawdopodobie ´nstwo wynosi 0,22.
u Metoda symulacyjna z wykorzystaniem Wolfram Mathematica 12.1.
sum = Table [ T o t a l [ RandomInteger [ 5 , 1 0 0 0 0 ] + 1 ] , 1 0 0 0 0 0 ] ; sum2 = S e l e c t [ sum , # < 35500 & ] ;
sum3 = S e l e c t [ sum2 , # > 34950 & ] ; Length [ sum3 ]
Histogram [ sum , 5 0 ]
Zaprezentowany fragment kodu kolejno:
-generuje list ˛e 100 000 sum 10 000 liczb od 1 do 6 (od 0 do 5 +1)(sum);
-przekształca list ˛e sum na list ˛e spełniaj ˛ac ˛a pierwszy warunek <
35050(sum2);
-przekształca list ˛e sum2 na list ˛e spełniaj ˛ac ˛a drugi warunek <
35500(sum3);
-zlicza ilo´s´c elementów elementów tej listy;
-dodatkowo generuje histogram dla pierwszej listy;
n o tat k i z ´c w i c z e ´n z p ro b a b i l i s t y k i 1 9/20 33
34 500 35 000 35 500
0 1000 2000 3000 4000
Rysunek 8:Histogram sum
W wersji podstawowej mo ˙zecie Pa ´nstwo policzy´c ile prób wpadło do zbioru
"r ˛ecznie".
Wynikiem Zlicze ´n w tym przypadku było 22 733 co przy podziele-niu przez liczb ˛e wszystkich przypadków daj ˛e szacowane prawdo-podobie ´nstwo zdarzenia równe 22,733%
u Narysuj zamkni ˛et ˛a krzyw ˛a mieszcz ˛ac ˛a si ˛e w kwadracie jednostko-wym. Oblicz przy pomocy metody Monte Carlo pole powierzchni ograniczone t ˛a krzyw ˛a. Spróbuj oszacowa´c dokładno´s´c oblicze ´n.
u Symulacja wykonana w j ˛ezyku Java wyznacza pole koła wpisane-go w kwadrat o boku 1:
p u b l i c s t a t i c void main ( S t r i n g [ ] a r g s ) { Random r = new Random ( ) ;
i n t a = 1 ; //dlugosc boku kwadratu o srodku w ( 1 / 2 , 1/2) i n t n = 1 0 0 0 0 0 0 ; // i l o s c punktow losowanych
i n t c t r = 0 ; // z l i c z a n i e wystapien double [ ] X = new double [ n ] ;
double [ ] Y = new double [ n ] ; f o r ( i n t i = 0 ; i <n ; i ++) {
X [ i ] = r . nextDouble ( ) ; Y [ i ] = r . nextDouble ( ) ;
i f ( Math . s q r t ( ( X [ i ]−0. 5 ) * (X[ i ]−0. 5 ) + (Y [ i ]−0. 5 ) * (Y[ i ]−0. 5 ) ) < 0 . 5 ) c t r ++;
}
System . out . p r i n t l n ( " Pole k o l a wpisanego w kwadrat o boku 1 : "
+ ( c t r /( double ) n ) ) ; }
Kilka przykładowych wyników działania programu: 0.785028, 0.785754, 0.785059, 0.785717, 0.785501. Obliczaj ˛ac pole na podsta-wie wzoru konwencjonalnego przybli ˙zaj ˛ac π = 3.14159
otrzy-Aby nie zmniejszy´c czytelno´sci kodu, pozwolono aby wchodził na margi-nes(oczywi´scie niektóre fragmenty nadal nale ˙zało przenie´s´c do nast ˛epnej linii
mujemy pole równe 0, 78540, co pozwala zauwa ˙zy´c ˙ze metoda ta jest dokładna (dla 106próbek i kwadratu o polu 1) do trzeciego miejsca po przecinku.
u Poni ˙zszy kod napisany w CERN ROOT 6.08 pozwala na oblicze-nie pola powierzchni dowolnej elipsy za pomoc ˛a metody Monte Carlo. Aby oblicza´c pole poza kwadratem jednostkowym, nale ˙zy zamieni´c definicj ˛e limitów.
void e l i p s a ( ) {
TRandom * gen = new TRandom ( ) ;
// u s t a w i e n i e parametrow p o l o s i a i b ; double osA = 0 . 3 ;
double osB = 0 . 4 ;
// u s t a w i e n i e p o l o z e n i a srodka e l p i s y ; double polozenieX = 0 . 5 ;
double polozenieY = 0 . 5 ;
// u s t a w i e n i e i l o s c i wylosowanych punktow ; c o n s t i n t rozmiar = 3 0 0 0 0 ;
// i l o s c punktow wewnatrz e l i p s y ; i n t num1 = 0 ;
// i l o s c punktow na zewnatrz e l i p s y ; i n t num2 = 0 ;
TGraph * graf1 = new TGraph ( ) ; TGraph * graf2 = new TGraph ( ) ; double x = 0 ;
double y = 0 ;
double g o r n y l i m i t X = 1 ; //polozenieX + 1 . 1 * osA ; double d o l n y l i m i t X = 0 ; //polozenieX − 1. 1 * osA ; double g o r n y l i m i t Y = 1 ; //polozenieY + 1 . 1 * osB ; double d o l n y l i m i t Y = 0 ; //polozenieY − 1. 1 * osB ; f o r ( i n t i = 0 ; i < rozmiar ; i + + ) {
x = gen−>Uniform ( dolnylimitX , g o r n y l i m i t X ) ; y = gen−>Uniform ( dolnylimitY , g o r n y l i m i t Y ) ;
i f ( TMath : : S q r t ( TMath : : Power ( ( x − polozenieX )/osA , 2 ) +
TMath : : Power ( ( y − polozenieY )/ osB , 2 ) ) <= 1 && num1 < rozmiar ) { g r a f 1−>S e t P o i n t ( num1 , x , y ) ;
num1++;
}
i f ( TMath : : S q r t ( TMath : : Power ( ( x − polozenieX )/osA , 2 ) +
TMath : : Power ( ( y − polozenieY )/ osB , 2 ) ) >1 && num2 < rozmiar ) { g r a f 2−>S e t P o i n t ( num2 , x , y ) ;
num2++;
}
n o tat k i z ´c w i c z e ´n z p ro b a b i l i s t y k i 1 9/20 35
}
TCanvas *kanwa1 =
new TCanvas ( " kanwa1 " , " P l o t t i n g Canvas 1 " , 1 5 0 , 1 0 , 9 9 0 , 9 9 0 ) ; g r a f 2−>Draw ( " PA " ) ;
g r a f 1−>Draw ( " Psame " ) ; g r a f 1−>SetMarkerColor ( 4 ) ; g r a f 2−>SetMarkerColor ( 6 ) ;
T E l l i p s e * kolo = new T E l l i p s e ( polozenieX , polozenieY , osA , osB ) ; kolo−>Draw ( " same " ) ;
kolo−> S e t F i l l S t y l e ( 3 0 0 8 ) ;
double Pole = ( 1 . 0 * num1/(num1 + num2 ) ) *
( ( g o r n y l i m i t X − d o l n y l i m i t X ) * ( gornylimitY − d o l n y l i m i t Y ) ) ; cout <<" Pole f i g u r y wynosi : "<< Pole <<endl ;
double niepewnosc = 1 . 0 / TMath : : S q r t ( num1 + num2 ) ;
cout <<"Niepewnosc pola f i g u r y wynosi : "<<niepewnosc <<endl ; }
Poni ˙zej znajduje si ˛e przedstawienie sposobu działania tego pro-gramu. Mo ˙zna zauwa ˙zy´c, ˙ze graf nie jest wypełniony całkowicie.
Poniewa ˙z nie wpływa to na wynik, postanowiono nie dodawa´c zb ˛ednych linijek kodu do programu.
Rysunek 9:Przykładowy graf
Pole figury zostało obliczone zgodnie ze wzorem Pf = na
nall ∗Pk
nato wszystkie punkty znajduj ˛ace si ˛e wewn ˛atrz figury, nall to wszystkie punkty, Pkto pole wygenerowanego kwadratu.
Niepewno´s´c pola została przyj ˛eta jako wzgl ˛edna dokładno´s´c obli-czenia całki
∆Pf = √1 nall
7Kod zapisany w Mathematice
Płaszczak- popularnonaukowa nazwa oznaczaj ˛aca istoty istniej ˛ace w prze-strzeni dwuwymiarowej. Zpolszczona nazwa stworze ´n z ksi ˛azki Flatlandia czy-li Kraina Płaszczaków autorstwa Edwina Abbotta.
Warto´sci zwrócone przez program:
Rysunek 10:przykładowe warto´sci zwracane przez program
u Wyznaczanie przy pomocy metody Monte Carlo pola koła Pk
wpisanego w kwadrat jednostkowy o polu P = 1 opiera si ˛e na poni ˙zszym równaniu:
Pk=Pk
n (37)
gdzie k to liczba punktów w okr ˛egu i na okr ˛egu, a n to liczba wszystkich punktów.
p o i n t s = RandomReal [ {−0. 5 , 0 . 5 } , { 1 0 0 0 0 , 2 } ] ; p l o t = L i s t P l o t [ p o i n t s , A s p e c t R a t i o −> 1 ] ; c i r c l e = Graphics [ C i r c l e [ { 0 , 0 } , 0 . 5 ] ] ; Show [ p l o t , c i r c l e ]
Count [Map[Norm, p o i n t s ] , _ ? ( # <= 0 . 5 & ) ] / 1 0 0 0 0 .
7
-0.4 -0.2 0.2 0.4
-0.4 -0.2 0.2 0.4
Rysunek 11:Wykres okr ˛egu o promie-niu R =0.5 wpisanego w kwadrat z wylosowanymi 10 000 punktów.
Dla n = 10000 obliczone pole wynosiło Pk = 0.7854, dla n = 100000 Pk = 0.7851, a dla n = 1000000 Pk = 0.7855. Gdy przyj-miemy π = 3.1415926535 pole tego koła obliczone z klasycz-nego wzoru z dokładno´sci ˛a pi ˛eciu miejsc po przecinku wynosi Pk =πr2=0.78539
u Pomimo, i ˙z w zadaniu mowa o krzywej zamkni ˛etej na jednostko-wej płaszczy´znie, człowiek nie płaszczak 3-ci wymiar lubi, wi ˛ec
n o tat k i z ´c w i c z e ´n z p ro b a b i l i s t y k i 1 9/20 37
postanowili´smy odwzorowa´c zadanie z u ˙zyciem metody Monte Carlo w trójwymiarowej przestrzeni.
Stworzono, wi ˛ec sfer ˛e jednostkow ˛a umieszczon ˛a w jednostkowym sze´scianie.
– Obj ˛eto´s´c sze´scianu Vc=a3=1
– Obj ˛eto´s´c sfery Vs= 43πr3= 43π(0, 5)3=0, 5235987756 Metoda Monte Carlo
Vc
Vs
≈ k N
gdzie N liczba losowo wybranych punktów w sze´scianie, za´s k liczba punktów które znalazły si ˛e wewn ˛atrz sfery.
Wizualizacja symulacji dla N=10000:
Wyniki dla symulacji przy N=10000:
k=5206 Vs ≈ k
N∗Vc= 5206
10000 ∗1=0, 5206
Wida´c ró ˙znice wynosz ˛ace 1001 realnej warto´sci obj ˛eto´sci.
Vs =0, 5235987756≈0, 5206
W teorii je´sli N→∞ to
k N ∗Vc
→Vs, spróbujmy wi ˛ec podnie´s´c N.
N=1000000 k=523473 Vs≈0, 523473
Dla tak du ˙zego N ró ˙znica obj ˛eto´sci wynosi około 10001 .
Obj ˛eto´s´c wyznaczon ˛a mo ˙zna na wykorzysta´c do oszacowania π:
π≈ 3 4
V r3 = 3
4
0, 523473
(0, 5)3 =3, 140838 u Zadanie zostało wykonane w Matlabie.
n o tat k i z ´c w i c z e ´n z p ro b a b i l i s t y k i 1 9/20 39
Wykres z lewej przedstawia 1000 punktów, a wykres z prawej 10000.
Mo ˙zna wyznaczy´c pole koła wpisanego w kwadrat jednostkowy (o polu 1) na podstawie liczby punktów, które zawarte s ˛a w tym ˙ze kształcie. Stanowi ˛a one jak ˛a´s cz ˛e´s´c całego kwadratu, zatem pole kształtu mo ˙zemy wyznaczy´c według takiego wzoru:
Pko =Pk
nin nall gdzie:
Pko - pole koła, Pk - pole kwadratu, n in - liczba punktów w kole, n all - liczba wszystkich punktów
Narysowane koło ma pole:
Pko =πr2=π· (0, 5)2=0,785 Pole koła wyznaczone na podstawie liczby punktów:
Dla 1000 punktów:
Pko =1 805
1000 =0, 805 Dla 10000 punktów:
Pko =1 7837
10000 =0, 7837 Dla 100000 punktów:
Pko =1 78521
100000 =0, 78521 Dla 1000000 punktów:
Zach ˛ecam do rozwi ˛aza ´n analitycznych, jak równie ˙z symulacji.
Pko =1 785412
1000000 =0, 7854
Wyznaczenie pola staje si ˛e dokładniejsze wraz ze wzrostem liczby punktów. W przybli ˙zeniu staje si ˛e równe warto´sci teoretycznej przy liczbie 100 000 punktów.
Skrypt w Matlabie:
%t a f u n k c j a z n a j d u j e s i e w o d d z i e l n y m s k r y p c i e f u n c t i o n h = c i r c l e ( x , y , r )
hold on
th = 0 : pi / 5 0 : 2 * pi ;
x u n i t = r * cos ( th ) + x ; y u n i t = r * sin ( th ) + y ;
h = p l o t ( xunit , yunit , ’ l i n e w i d t h ’ , 1 . 5 ) ; hold o f f
end
rng ( 0 , ’ t w i s t e r ’ ) a = −0.5;
b = 0 . 5 ;
x = ( b−a ) . * rand ( 1 0 0 0 0 , 1 ) + a ; y = ( b−a ) . * rand ( 1 0 0 0 0 , 1 ) + a ; c o u n t e r = 0 ;
hold on
s c a t t e r ( x , y , ’ . ’ ) ; h = c i r c l e ( 0 , 0 , 0 . 5 ) ; hold o f f
f o r i = 1 : length ( x )
i f( hypot ( x ( i ) , y ( i ) ) <= 0 . 5 ) c o u n t e r = c o u n t e r + 1 ; end
end
u Dodajemy 10 000 liczb, zaokr ˛aglonych z dokładno´sci ˛a do 10−m. Bł ˛edy zaokr ˛aglenia s ˛a niezale ˙znymi zmiennymi losowymi o roz-kładzie jednostajnym w przedziale(−10−m/2, 10−m/2). Wyzna-czy´c przedział, w którym z prawdopodobie ´nstwem 0.99 b ˛edzie si ˛e zawierał bł ˛ad sumy.
n o tat k i z ´c w i c z e ´n z p ro b a b i l i s t y k i 1 9/20 41
u n = 10000 liczb zaokr ˛aglonych do ok. 10−m Xk- niepewno´s´c k-tej liczby
Sn =
∑
n k=1Xk
dla P(X)=0,99 wyznaczamy przedział w którym b ˛edzie zawierał si ˛e bł ˛ad sumy
z prawdopodobie ´nstwem 0,99:
Sn
u Metoda symulacyjna z wykorzystaniem Wolfram Mathematica 12.1.
Zaprezentowany fragment kodu kolejno:
-generuje list ˛e 100 000 sum 10 000 liczb od -1/2 do 1/2 (od 0 do 1 -1/2)(sum);
-sortuj ˛e t ˛a list ˛e;
-wypisuje 500 i 99 500 element tej listy;
-dodatkowo generuje histogram dla pierwszej listy;
-100 -50 0 50 100
0 1000 2000 3000 4000 5000 6000
7000 Rysunek 12:Histogram sum
Rozkład wyników powinien by´c symetryczny wzgl ˛edem zera, dlatego bior ˛ac pod uwag ˛e wygenerowanie 100 000 wyników od-rzucamy po 500 z ka ˙zdego skaju wybieraj ˛ac element znajduj ˛acy si ˛e kraw˛edzi naszego zbioru w którym znajduje si ˛e 99% wyników.
Otrzymane w ten sposób wyniki wynosz ˛a w tym przypadku od-powiednio−74.1079 oraz 74.2388. Moduły wyników s ˛a do siebie zbli ˙zone co zgadza si ˛e z intuicj ˛a o symetrii rozkładu. Z cało´sci wynika ˙ze po zsumowaniu 10000 liczb w ten sposób nasze za-okr ˛aglenie spowoduje bł ˛ad mieszcz ˛acy si ˛e, w 99% przypadków, w przedziale od około−74.1079∗10−m do 74.2388∗10−m.
n o tat k i z ´c w i c z e ´n z p ro b a b i l i s t y k i 1 9/20 43