---
title: "DCv 5.12-- Problém společných životů"
output:
  pdf_document: default
  html_notebook: default
---

# První verze

Lidi umírají náhodně, pokud je počet přeživších jiný, než $a$, tak výsledek zahodíme. 
Neboli simuluje úlohu přesně podle zadání. 

```{r}
m = 10
a = 4
p = 0.2

N = 10^5
suma = sumb = sumd = sumf = 0
suma = sumb = sumd = sumf = 0
cnta = 0

M = c(rep(1,a),rep(0,m-a))   # libovolná fixovaná množina, využíváme toho, že a <= m

for (i in 1:N) {
  x = rbinom(m,1,p)   # přeživší ženy
  y = rbinom(m,1,p)   # přeživší muži  
  if (sum(x)+sum(y)==a) {
    cnta = cnta + 1
    sumd = sumd + sum(x&y)
    suma = suma + x[1]
    sumb = sum(sumb, (x[1])&(y[1]))
    sumf = sumf + all(x==M)
  }
}
```

Všechy pravděpobonosti/střední hodnoty uvedeme napřed nasamplované, pak vypočtenou hodnotu (pro srovnání).
Pravděpodobnost v části (a):
```{r}
suma/cnta
a/(2*m)
```
Pravděpodobnost v části (b):
```{r}
sumb/cnta
a*(a-1)/(2*m)/(2*m-1)
```

Střední hodnota v části (d):
```{r}
sumd/cnta
a*(a-1)/(2*(2*m-1))
```
m = 10
a = 4

N = 10^5
suma = sumb = sumd = sumf = 0
Pravděpodobnost P(L=M | A=a):
```{r}
sumf/cnta
1/choose(2*m,a)
```
Pravděpodobnost P(L=M). Tady (i v minulé části) si můžeme všimnout, že výsledek se podstatně odlišuje od vypočteného. 
Je to tím, že pravděpodobnost je příliš malá na to, kolik pokusů jsme udělali. Pro $N=10^6$ už se přesnost vylepší. 

```{r}
sumf/N
p^a*(1-p)^(2*m-a)
```


# Druhá verze

Vybereme náhodnou podmnožinu velikosti $a$. 
Musíme si rozmyslet, že je to totéž -- ale zato pak máme efektivnější výpočet. 
(Pravděpodobně, neměřil jsem.)
Zde $p$ nepotřebujeme. Nemůžeme z principu ověřit $P(L=M)$, všechny 
pravděpodobnosti a střední hodnoty budou podmíněny jevem $A=a$. 

```{r}
m = 10
a = 4

N = 10^5
suma = sumb = sumd = sumf = 0

people = c(rep(1,a),rep(0,2*m-a))   # máme a živých (1), zbytek mrtvých (0)
for (i in 1:N) {
  shuffle <- sample(people)   # náhodná permutace
  x = shuffle[1:m]
  y = shuffle[(m+1):(2*m)]  
  stopifnot(sum(x)+sum(y)==a)  # jenom pro jistotu, mělo by být vždy splněno
  
  sumd = sumd + sum(x&y)
  suma = suma + x[1]
  sumb = sum(sumb, (x[1])&(y[1]))
  sumf = sumf + all(shuffle==people)  
}

print(suma/N)
```


Všechy pravděpobonosti/střední hodnoty uvedeme napřed nasamplované, pak vypočtenou hodnotu (pro srovnání).
Pravděpodobnost v části (a):
```{r}
suma/N
a/(2*m)
```
Pravděpodobnost v části (b):
```{r}
sumb/N
a*(a-1)/(2*m)/(2*m-1)
```

Střední hodnota v části (d):
```{r}
sumd/N
a*(a-1)/(2*(2*m-1))
```

Pravděpodobnost P(L=M | A=a):
```{r}
sumf/N
1/choose(2*m,a)
```