Zajęcia35listopada2021 DyskretnepróbkowanieDorotaCelińska-Kopczyńska Rachunekprawdopodobieństwaistatystyka

Rachunek prawdopodobieństwa i statystyka

Dyskretne próbkowanie

Dorota Celińska-Kopczyńska

Uniwersytet Warszawski

Zajęcia 3 5 listopada 2021

Idea zajęć – co i po co będziemy robić?

I Podczas poprzednich zajęć próbkowaliśmy mając znane prawdopodobieństwa

I Wykorzystaliśmy sumy prefiksowe i wyszukiwanie binarne I Próbkowanie, tworzenie symulacji jest kluczowym elementem

warsztatu statystyka (generalizacja!)

I Niestety, kod służący do losowania jest wąskim gardłem całego programu

I Podczas dzisiejszych zajęć będziemy starać się usprawnić losowania (techniki o złożoności O(1))

Problem 2a

Rozważamy problem próbkowania mając znane

prawdopodobieństwa p1, p2, . . . , pd. Załóżmy, że mamy gotową maszynkę, który pozwala nam losować z prawdopodobieństw q₁, q₂, . . . , q_d, takich, które są bliskie p₁, p₂, . . . , p_d:

∀_{i =1,...,n}(1 − ε)qi ¬ p_i ¬ (1 + ε)q_i.

I Jak wykorzystać ten gotowy kod, żeby próbkować z prawdopodobieństw p?

I Przedyskutuj poprawność algorytmu

I Popraw rozwiązanie Zadania 1c, wykorzystując omówione podejście

Czym się zajmujemy (?)?

I Pracujemy z dwoma dyskretnymi zmiennymi losowymi: P i Q.

I Znamy rozkłady zarówno dla P jak i Q, dla Q I Maszynka (bez

dostępnych szczegółów) pozwala losować

z rozkładu Q 0

10 20 30 40

2 4 6

p

count

A

0 20 40 60

2 4 6

q

count

B

0 20 40 60

1 2 3 4 5 6 value

count

var p q C

Jak będziemy to robić?

I Wykorzystujemy fakt, że prawdopodobieństwa p są niższe niż q

I Losujemy z rozkładu Q wartość R

I Losujemy wartość zmiennej losowej X (rozkład jednostajny) I Jeśli X < p(R), wartość zwracamy, w przeciwnym przypadku

powtarzamy losowanie

Wniosek?

Wykorzystana procedura zaimplementowana z wykorzystaniem standardowej składni Pythona jest bardzo wolna... zyskaliśmy niewiele w porównaniu do rozwiązania poprzedniego zadania.

Czy nie da się tego poprawić?

Problem 2b – zadanie zaliczeniowe do 3.12.2021 23:59

Popraw rozwiązanie Zadania 2a, tak, żeby liczba wykonywanych pętli została zminimalizowana. Wykorzystaj wektoryzację. Użyj samplera do rozwiązania problemu urodzinowego (wynikiem może być np. histogram).

Biblioteka numpy

I Biblioteka numpy umożliwia korzystanie z optymalizowanej algebry tablic (tablice zachowują się jak tablice w C;

zoptymalizowano pod nie operacje)

I Na tablicach numpy możemy wykonywać obliczenia i porównania

I Możemy wygenerować w jednym kroku całą tablicę numpy liczb losowych

I Można wykorzystać tablice numpy jako indeksowanie dla innej tablicy numpy

import numpy as np X = np.array([10,3,7]) I = np.array([1,1,2,2]) print(X[I])

Porównanie rozwiązań

while:

R = rand(0, n-1) X = rand(0, M-1) if X < p[R]:

return R

R = np.random.randint(0, n-1, N) X = np.random.randint(0, M-1, N) return R[X < p[R]]

R[X < p[R]]

Podana linijka zawiera trzy kroki:

1. Indeksowanie tablicy p za pomocą wartości tablicy R 2. Porównanie parami wszystkich wartości tablic X i p[R] –

zwraca tablicę booli

3. Indeksowanie tablicy R tablicą booli (filtrowanie)

Oczekiwany rezulat

Problem 2c

W tym problemie ponownie będziemy próbkować ze znanych prawdopodobieństw p1, p2, . . . , pn, tylko tym razem bez dodatkowych założeń odnośnie p_i.

Zaimplementuj metodę Alias. Wykorzystaj wektoryzację (numpy where). Zastosuj tę metodę do problemu paradoksu urodzinowego (wynikiem działania całego programu powinien być histogram).

Metoda alias

I Niech V będzie średnią wartością z p_i, czyli M = ¹_n

I Stwórz n kubełków, każdy o pojemności V , włóż każde p_i do osobnego kubełka

I Póki istnieje kubełek A, który nie jest pełny, znajdź kubełek B, który jest przepełniony i przenieś prawdopodobieństwo z B do A, aż do zapełnienia kubełka A

I Mając zapełnione kubełki losujemy jeden z nich. Losujemy wartość X z rozkładu jednostajnego i porównujemy, którą z wartości w kubełku zwrócić

Metoda alias – przykład

I Wartości, którymi wypełniamy kubełki pochodzą z pliku I Ten algorytm zawsze zostanie ukończony

I Gdy się zakończy każdy z kubełków będzie zawierać dwie wartości (druga może być zerowa)