통계학 입문 제공 Julia (1) - 표준 득점, 편차 득점 - Qiita
통계학 입문 동행 Julia (2) - 피어슨의 적률 상관 계수 - Qiita
통계학 입문의 공 Julia (3) - 순상관계수와 자기상관계수 - Qiita
통계학 입문의 공공 Julia (4) - 최소제곱법과 결정계수 - Qiita
통계학 입문 동행 Julia (5) - 베이즈 정리 - Qiita
통계학 입문 동행 Julia (6) - 드 메레 문제 = Oiita

통계학을 공부하기 위해 「통계학 입문」(도쿄 대학 출판회)을 읽기 시작했습니다. 연습 문제가 있습니다만, 이것도 공부중의 Julia를 사용해 계산해 가고 싶습니다.

「확률론 입문」(치쿠마 학예 문고, 아카시야)도 참고로 하고 있습니다. 확률론의 기초적인 사고 방식을 자세히 설명합니다. 특히 다중 시도의 설명은 상세합니다. 문고 책이므로 가지고 쉽고 읽기 쉽습니다.

조합을 계산하려면 Julia 언어의 Base Mathematics 범주의 binomial 함수를 사용합니다. Base이기 때문에 Mathematics 카테고리의 함수는 특별히 아무것도하지 않고도 사용할 수 있습니다.
Mathematics
htps : // / cs. 주아아 g. 오 rg/엔/v1/바세/마 th/

1. 드 멜레 문제

P85 질문 4.1

질문 i)

주사위를 4회 던질 때, 6번째 눈이 적어도 1회 나오는데 베팅할까, 1회도 나오지 않는 분에게 베팅할까?

세 가지 방법으로 해결해 보겠습니다.

\begin{align}
&サイコロを振りどの目が出るかを確認する試行を考える。\\
&\Omega : 全事象 \\
&E_i \quad (i=1,2,3,4): \Omega の部分集合とするとき以下のように定義する。\\
&(E_1, E_2, E_3, E_4) \equiv \{ (a_1, a_2, a_3, a_4) | a_i \in E_i  \quad 
 i=1,2,3,4\}\\
&これは4回続けてサイコロを振る事象を表している。\\
\\
&今サイコロを振った時に、6が出る事象をEとする。\\
\\
&4回続けてサイコロを振った時に少なくとも1回は6が出る事象は、\\
&4回連続して6以外が出る事象の補集合と考えられる。\\
&また4回連続して6以外が出る事象は以下のように表現できる。\\
&F^{(4)} \equiv (E^c, E^c, E^c, E^c)\\
\\
&故に\\
&P(F^{(4)}) = P(E^c) P(E^c) P(E^c) P(E^c) = P(E^c)^4=(5/6)^4\\
\\
&求める確率は 1-(5/6)^4 = 0.5177469135802468\\
\\
\\
\\
&************************************\\
\\
&別解です。\\
&補集合で考えるのではなく、もう少し正面から解いてみます。\\
\\
&4回続けてサイコロを振った時に、少なくとも1回は6が出る事象は以下のように表現される。\\
&E^{(4)} \equiv (E, \Omega, \Omega, \Omega) \cup (\Omega, E, \Omega, \Omega) \cup (\Omega, \Omega, E, \Omega) \cup (\Omega, \Omega, \Omega, E) \tag{1}\\
\\
\\
&加法の定理 : P(A \cup B) =P(A) + P(B) - P(A \cap B) \tag{2}\\
\\
&加法の定理を繰り返し使うことで、以下が導かれる\\
&P(A \cup B \cup C \cup D) = P(A) + P(B) + P(C) +P(D)\\
&\quad 
 - P(A \cap B) - P(A \cap C) - P(A \cap D) - P(B \cap C) - P(B \cap D) - P(C \cap D)\\
&\quad + P(A \cap B \cap C) + P(A \cap B \cap D) + P(A \cap C \cap D) + P(B \cap C \cap D)\\
&\quad - P(A \cap B \cap C \cap D) \tag{3}\\
\\
&A=(E, \Omega, \Omega, \Omega), \quad B= (\Omega, E, \Omega, \Omega), \\
&C=(\Omega, \Omega, E, \Omega),\quad D=(\Omega, \Omega, \Omega, E)として(3)を置き換える。\\
\\
\\
&P(E^{(4)}) =P(E, \Omega, \Omega, \Omega)+P(\Omega, E, \Omega, \Omega)+P(\Omega, \Omega, E, \Omega)+P(\Omega, \Omega, \Omega, E)\\
&-P(E, E, \Omega, \Omega)-P(E, \Omega, E, \Omega)-P(E, \Omega, \Omega, E)-P(\Omega, E, E, \Omega)\\
&-P(\Omega, E, \Omega, E)-P(\Omega, \Omega, E, E)\\
&+P(E, E, E, \Omega)+P(E, E, \Omega, E)+P(E, \Omega, E, E)+P(\Omega, E, E, E)\\
& - P(E, E, E, E)\\
&=4*(1/6) -6*(1/6)^2+4*(1/6)^3-(1/6)^4=0.5177469135802469
\\
\\
&誤差を除いて、計算結果は同じです。\\
&最初の方が簡単だけど、別解の方が全体の構造が良くわかり汎用性があります。\\
\\
\\
&************************************\\
\\
&別解２です。\\
&2項分布の考え方で解いてみます。\\
\\
&確率変数 X_i^{(4)} を以下のように定義する。\\
\\
&X_i^{(4)}(a_1,a_2,a_3,a_4) = 1 \qquad (a_i \in E)\\
&X_i^{(4)}(a_1,a_2,a_3,a_4) = 0 \qquad (a_i \notin E)\\
\\
&この時、確率変数S^{(4)}を以下のように定義する。\\
\\
&S^{(4)} \equiv X_1^{(4)} + X_2^{(4)} + X_3^{(4)} + X_4^{(4)}\\
\\
&S^{(4)}(a_1,a_2,a_3,a_4) = k は、\\
&a_1,a_2,a_3,a_4 のうちEに属するものの個数がk個であることを意味します。\\
&当然、0 \leqq k \leqq 4　が成り立ちます。\\
\\
&\{ S^{(4)}=k \} \equiv \{(a_1,a_2,a_3,a_4) \in (\Omega,\Omega,\Omega,\Omega) \quad | \quad S^{(4)}
(a_1,a_2,a_3,a_4)=k \}
\\
\\
&例えばk=3の時は、以下のように4個から3個を選び出す組み合わ数\begin{pmatrix} 4 \\3 \end{pmatrix} 通りの\\
&お互いに背反な事象に分割されます。\\
\\
&\{ S^{(4)}=3 \} = (E,E,E,E^c) \cup (E,E,E^c,E) \cup (E,E^c,E,E) \cup (E^c,E,E,E)
\\
\\
\\
&一般に\{ S^{(4)}=k \}は、4個からk個を選び出す組み合わ数\begin{pmatrix} 4 \\k \end{pmatrix} 通りの\\
&お互いに背反な事象に分割されます。\\
&それぞれの事象の確率は、k個のEと、(4-k)個のE^cですからP(E)^k(1-P(E)^{(4-k)}です。\\
&P(\{ S^{(4)}=k \}) = \begin{pmatrix} 4 \\k \end{pmatrix} P(E)^k(1-P(E))^{(4-k)}\\
\\
\\
&各 \{ S^{(4)}=k \} \quad i=1,2,3,4は明らかに排反だから、求めるものは以下のようになる。\\
&P(\{ S^{(4)}=1 \})+P(\{ S^{(4)}=2 \})+P(\{ S^{(4)}=3 \})+P(\{ S^{(4)}=4 \})\\
\\
&後はJuliaのプログラムで解いてみたいと思います。
\\
&\qquad \qquad  \qquad \qquad  \qquad \qquad  \qquad \qquad  \qquad \qquad  \qquad \qquad  \qquad \qquad  \qquad\\
\end{align}

다음의 프로그램은 보다 범용이 되고 있습니다만, n=4, k=1로 실행하면 문제와 같습니다.

function mybinomial(k, n, p)
    s = 0.0
    for i in k:n
        s += binomial(n,i)*(p)^i*(1-p)^(n-i)
    end
    s
end

mybinomial(1,4,1/6)
0.5177469135802469

2항 분포로 풀어도 같은 결과군요.

질문 ii)

주사위를 2회 동시에 24회 던질 때, (6,6)의 눈이 적어도 1회 나오는데 베팅할까, 1회도 나오지 않는 쪽에 베팅할까?

(6,6)の目以外が24回連続して出る確率は以下の計算で求まる。
(35/36)^24
0.5085961238690966

2項分布でも求まる

(6,6)の目が少なくとも1回出る確率
mybinomial(1,24,1/36)
0.49140387613090303

(6,6)の目以外が24回連続して出る確率
1-mybinomial(1,24,1/36)
0.508596123869097

2. 호이헨스 문제

P85 질문 4.2

2개의 주사위를 몇 번 던지면, 그 중 1회는 눈의 합이 12가 될 확률이 0.9를 넘는가?

본서의 해답과는 별회입니다. 이항 분포를 사용하여 다음과 같이 확인할 수 있습니다.

mybinomial(1,big(81),1/36)
8.97903928782527159074445136...e-01

mybinomial(1,big(82),1/36)
9.00739930760790253942017745...e-01

즉 82회 던지면 된다.
big를 사용하지 않으면 OverflowError가 됩니다. big 님이군요.

3. 간호사 문제

P85 질문 4.3

30명의 간호사를 15명씩의 2조의 교대 근무으로 나누는 나누는 방법은 무엇인가?
조합의 문제군요.

binomial(30,15)
155117520

4. 생일 문제

P85 질문 4.4

r인 집단중에 같은 생일의 사람이 적어도 2명 있을 확률을 생각하는 문제입니다.

아무래도 제대로 오지 않기 때문에. . . . . 조금 문제의 문장을 나름대로 씹고 있습니다. 확률에 대한 시도의 의미를 명확히하고 싶었기 때문에 본질은 변하지 않습니다.

\begin{align}
\\
&r人の集団を考える。\\
&壺の中に365個の玉があり、それぞれ1,2,3,..,365の番号が振られている。\\
\\
&神様はr人の誕生日を以下のrステップの手順で決める。\\
&(1)神様が壺から1個取り出し、1人目の誕生日を決める。印をつけて玉は戻す。\\
&(2)神様が壺から1個取り出し、2人目の誕生日を決める。印をつけて玉は戻す。\\
&.....\\
&(r)神様が壺から1個取り出し、r人目の誕生日を決める。印をつけて玉は戻す。\\
\\
&同じ誕生日の人が一組もいないことは、\\
&壺の中に、2つ以上の印ついているボールが存在しないことと同じ意味です。(*)\\
&そのような確率を、各ステップにおいて求めてみましょう。\\
\\
&(1)ステップにおいて印のついていない玉を選択する確率は1です。\\
&(2)ステップにおいて印のついていない玉を選択する確率は(365-1)/365です。\\
&(3)ステップにおいて印のついていない玉を選択する確率は(365-2)/365です。\\
&.....\\
&(r)ステップにおいて印のついていない玉を選択する確率は(365-r+1)/365です。\\
\\
&つまり2つ以上の印ついているボールが存在しない確率は以下の式になります。\\
\\
&P = ((365-1)/365)((365-2)/365).....((365-r+1)/365)
\\
&\qquad \qquad  \qquad \qquad  \qquad \qquad  \qquad \qquad  \qquad \qquad  \qquad \qquad  \qquad \qquad  \qquad \qquad \qquad  \qquad\\
\end{align}

확률을 구하는 공식을 알았으므로 Julia에서 프로그래밍합니다. r을 매개변수로 하여 1에서 60까지 변경하고 각 확률을 계산합니다. 마지막으로 그래프로 그립니다.

n=60
rs=zeros(n)
ps=zeros(n)

function mybirth(r)
    p = 1.0
    for i in 1:r
        p *= ((365-i+1)/365)
    end
    1-p
end

function mybirth_r()
    for r in 1:n
        p = mybirth(r)
        println("r=",r," p=",p)
        rs[r] = r
        ps[r] = p 
    end
    1
end

mybirth_r()

using Plots
plot(rs,ps,seriestype=:scatter)

출력 결과는 다음과 같습니다.

r=1 p=0.0
r=2 p=0.002739726027397249
...
r=5 p=0.02713557369979347
...
r=10 p=0.11694817771107768
...
r=15 p=0.25290131976368646
...
r=20 p=0.41143838358058027
r=21 p=0.443688335165206
r=22 p=0.4756953076625503
r=23 p=0.5072972343239857
r=24 p=0.538344257914529
r=25 p=0.568699703969464
...
r=30 p=0.7063162427192688
...
r=35 p=0.8143832388747153
...
r=40 p=0.891231809817949
...
r=50 p=0.9703735795779884
...
r=60 p=0.994122660865348

60명이면 같은 생일 쌍이 99% 이상의 확률로 존재하게 됩니다.
r을 x축, 확률을 y축으로 하여 그린 그래프는 다음과 같습니다.



또한 이 문제는 다음 사이트에서 자세히 설명합니다.

왜 "23명 있으면 같은 생일 사람이 있을 확률은 50%"인가

생일이 일치할 확률은 얼마입니까? - 통계학 입문

같은 생일 2인조가 있을 확률에 대해 - 고등학교 수학의 아름다운 이야기

이번은 이상입니다.