Juliaで学ぶ確率変数(7) - 指数分布（連続型）

[Juliaで学ぶ確率変数(1) - 確率変数の定義 - Qiita](https://qiita.com/sand/items/afb5976ad28120d216c4)
[Juliaで学ぶ確率変数(11) - まとめ - Qiita](https://qiita.com/sand/items/e4b32f7cfd74dea40298)



　確率変数を勉強中ですが、**「確率統計」（森北出版）**は数学的に明確な定義がしっかり書かれているので、これを中心に勉強しています。あわせて**「統計学入門」（東京大学出版会）**と**「確率論入門」（ちくま学芸文庫、赤攝也）**も併読しています。

　本記事は、それらの教科書を読みながら、実際に例題や問題をJuliaで解いていく試みです。Juliaの連続型確率変数のライブラリのドキュメントです。　==>[Distributions/Univariate/ContinuousDistributions](
https://juliastats.github.io/Distributions.jl/latest/univariate.html#Continuous-Distributions-1)


#１．指数分布 Ex(λ)

指数分布 Ex(λ)は、**『連続的な待ち時間分布』と呼ばれ、期間1/λあたりに1回起きると期待される事象が1回起きるまでの時間の分布**、という意味を持つとされています。**幾何分布**と同じく待ち時間分布という性質を持っています。機械の部品などの故障時間の分布として使われることが多い。

[ガンマ分布のはなし - 統計学といくつかのよしなしごと](http://yanbow221.hatenablog.com/entry/2017/09/06/211456)

```math
\begin{align}
\\
\\
&\lambdaを正の定数とする。\\
\\
&X: \Omega \rightarrow 非負の実数値\\
\\
&確率密度関数f(x)が以下のように書けるとき、Xは指数分布Ex(\lambda)に従うという。\\
\\
&f(x) =
\left\{
\begin{array}{ll}
\lambda e^{-\lambda x} & (x \geqq 0 \;のとき) \\
0 & (x<0\;のとき)
\end{array}
\right.\\
\\
\\
&ここで、以下が成り立つ。\\
&p(\Omega)=\int_{-\infty}^{\infty} f(x) dx = \int_{0}^{\infty} \lambda e^{-\lambda x} dx = 1\\
&E[X] = \int_{0}^{\infty} xf(x) dx = \frac {1} {\lambda}\\
&V[X] = \int_{0}^{\infty} (x - E[X])^2 f(x) dx = \frac {1} {\lambda^2}\\
&F(x) = p(\{X \leqq x\}) 
= \int_{-\infty}^{x} f(t) dt 
= \int_{0}^{x} \lambda e^{-\lambda t} dt = 1-e^{-\lambda x}\\
\\
&\qquad \qquad  \qquad \qquad  \qquad \qquad  \qquad \qquad \qquad  \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad  \qquad  \\
\end{align}
```




##１－１．p(Ω)の証明


解析学の練習にp(Ω)とE[X]、F(x)の証明を追ってみたいと思います。「確率統計」（森北出版）の証明は全く分からなかったので、「改訂版 解析学大要」（養賢堂）で解析学を復習しつつ試みました。


```math
\begin{align}

\\
&p(\Omega)=\int_{-\infty}^{\infty} f(x) dx = \int_{0}^{\infty} \lambda e^{-\lambda x} dx\\
\\
&上の定義から以下を導きます。\\
&p(\Omega)  = 1\\
\\
\\
&まず \; x=\varphi(t) = -\frac{t}{\lambda} \; とすると、t = -\lambda x となります。\\
&\varphi(0)=0,\; \varphi(-\infty)=\infty \; であるから、置換積分により、\\
&p(\Omega) = \int_{0}^{\infty} \lambda e^{-\lambda x} dx 
= \lambda \int_{0}^{-\infty} e^{t} \frac{dx}{dt} dt 
= \lambda \int_{0}^{-\infty} e^{t} \frac{-1}{\lambda} dt  \\
&= -\int_{0}^{-\infty} e^{t} dt 
= -[e^t]_0^{-\infty}\\
\\
\\
&一般に、\lim_{x \to -\infty} e^x = 0,\; e^0=1が成り立つから、\\
\\
&p(\Omega) = -(0-1)=1  \qquad \qquad  Q.E.D.\\
\\
&\qquad \qquad  \qquad \qquad  \qquad \qquad  \qquad \qquad \qquad  \qquad \qquad \qquad \qquad \qquad \qquad \qquad  \\
\end{align}
```



##１－２． E[X]の証明

```math
\begin{align}

\\
&E[X]=\int_{0}^{\infty} xf(x) dx = \int_{0}^{\infty} x \lambda e^{-\lambda x} dx\\
\\
&上の定義から以下を導きます。\\
&E[X]  = \frac {1} {\lambda}\\
\\
\\
&まず \; x=\varphi(t) = \frac{t}{\lambda} \; とすると、t = \lambda x となります。\\
&\varphi(0)=0,\; \varphi(\infty)=\infty \; であるから、置換積分により、\\
&E[X] = \int_{0}^{\infty} x \lambda e^{-\lambda x} dx = \int_{0}^{\infty} t e^{-t} \frac{dx}{dt} dt =\int_{0}^{\infty} t e^{-t} \frac{1}{\lambda} dt  = \frac{1}{\lambda} \int_{0}^{\infty} t e^{-t} dt \\
\\
&また\frac{d\;(-e^{-t})}{dt}=e^{-t}だから、部分積分により、\\
&\lambda \; E[X] = \int_{0}^{\infty} t e^{-t} dt
= \int_{0}^{\infty} t \frac{d\;(-e^{-t})}{dt} dt\\
&= [t(-e^{-t})]_0^\infty - \int_{0}^{\infty} (-e^{-t})dt \\
&= [t(-e^{-t})]_0^\infty - [e^{-t}]_0^\infty\\
\\
&一般に、\lim_{x \to -\infty} e^x = 0,\; \lim_{x \to \infty} xe^{-x} = 0 ,\; e^0=1が成り立つから、\\
&\lambda \; E[X] = 0 - (-1) = 1
\\
&よって \; E[X]=\frac {1} {\lambda}  \qquad \qquad  Q.E.D.\\
\\
&\qquad \qquad  \qquad \qquad  \qquad \qquad  \qquad \qquad \qquad  \qquad \qquad \qquad \qquad \qquad \qquad \qquad  \\
\end{align}
```

##１－３．F(x)の証明



```math
\begin{align}

\\
&F(x) = p(\{X \leqq x\}) 
= \int_{-\infty}^{x} f(t) dt 
= \int_{0}^{x} \lambda e^{-\lambda t} dt\\
\\
&上の定義から以下を導きます。\\
&F(x)  = 1-e^{-\lambda x}\\
\\
\\
&「１－１．P(\Omega)の証明」と同じようにして以下が成り立つ。\\
\\
&まず \; t=\varphi(s) = -\frac{s}{\lambda} \; とすると、s = -\lambda t となります。\\
&\varphi(0)=0,\; \varphi(-\lambda x)=x \; であるから、置換積分により、\\
&F(x) = \int_{0}^{x} \lambda e^{-\lambda t} dt 
= \lambda \int_{0}^{-\lambda x} e^{s} \frac{dt}{ds} ds 
= \lambda \int_{0}^{-\lambda x} e^{s} \frac{-1}{\lambda} ds  \\
&= -\int_{0}^{-\lambda x} e^{s} ds 
= -[e^s]_0^{-\lambda x}=1-e^{-\lambda x}\\
\\
&Q.E.D.\\
\\
&\qquad \qquad  \qquad \qquad  \qquad \qquad  \qquad \qquad \qquad  \qquad \qquad \qquad \qquad \qquad \qquad \qquad \\
\end{align}
```

計算するとき便利でした。==>[「数列の極限」](http://www.sist.ac.jp/~kanakubo/research/hosoku/kyokugen.html)





##１－４．Juliaで解く例題(1)

Juliaで分布を扱うための基礎知識は以下のページを参照してください。
[「Juliaで学ぶ確率変数(1) - 確率変数の定義」の「４．Juliaで確率分布を扱う」](https://qiita.com/sand/items/afb5976ad28120d216c4)



以下のサイトの例題を考えてみます。

[練習問題（15. いろいろな確率分布3）- 統計WEB]
https://bellcurve.jp/statistics/course/8282.html

**--- 問題 ---**


**Kさんのもとには、日々たくさんのお問い合わせメールが届く。お問い合わせメールに対する返信メールの作成時間は1通あたり平均15分の指数分布に従うと仮定できる。**

１．作成時間が30分以上となる確率を求めよ。
２．また、作成時間が5分以上10分以下となる確率を求めよ。


**--- 解答 ---**

```math
\begin{align}
&まず\lambdaを求めます。指数分布の期待値の公式より、E[X]=\frac{1}{\lambda}=15なので、\\
&\lambda=\frac{1}{15}と求められます。\\
\\
\\
&メールの作成時間が30分以下となる確率は、分布F(30)と表現できます。\\
&F(30) = 1-e^{-\lambda 30} = 1-e^{-\frac{1}{15} 30}= 1-e^{-2}\\\
\\
&作成時間が30分以上となる確率は、余事象の確率で、1-F(30)、となります。
\\
\\
\\
&作成時間が5分以上10分以下となる確率は分布で、F(10) - F(5)、となります。\\
&F(10) - F(5) = (1-e^{-\lambda 10}) - (1-e^{-\lambda 5}) = e^{-\frac{1}{15} 5}-e^{-\frac{1}{15} 10}\\
\\
\\
&これ以上の具体的な計算は、Juliaに任せます。\\
&Q.E.D.\\
\\
&\qquad \qquad  \qquad \qquad  \qquad \qquad  \qquad \qquad \qquad  \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \\
\\

\end{align}
```


**--- Juliaの計算 ---**

JuliaのパッケージDistributionsでは、**確率密度関数**(Probability Density Function = PDF)の評価は**pdf**関数でした。**分布（累積分布関数）** (Cumulative Distribution Function = CDF)の評価は**cdf**関数を使います。


```:Distributions.cdf
Distributions.cdf — Method.
cdf(d::UnivariateDistribution, x::Real)
Evaluate the cumulative probability at x.
```


Juliaでは正規分布の確率密度の定義はDistributions.Exponentialを使います。\lambdaが与えられれば、**Exponential(1/lambda)**で定義します。Exponential(lambda)でないことに注意してください。


```julia
using Plots
using Distributions
using StatPlots

d=Exponential(15)
plot(d,fill=(0, .5,:orange))

println("F(30) = ", cdf(d,30))
println("1 - F(30) = ", 1-cdf(d,30))
println("F(10) - F(5) = ", cdf(d,10) -cdf(d,5))

plot(d,fill=(0, .5,:orange))
```

出力は以下のようになります。

```
F(30) = 0.8646647167633873          # メールの作成時間が30分以下となる確率
1 - F(30) = 0.1353352832366127      # 作成時間が30分以上となる確率
F(10) - F(5) = 0.20311419154119725  # 作成時間が5分以上10分以下となる確率
```


![image.png](https://qiita-image-store.s3.amazonaws.com/0/61915/fd8ec83b-037b-dd3d-7995-1a86e0dad20b.png)


**確率密度関数と分布**の関係をグラフで表します。

```
d=Exponential(15)
scatter(d, leg=false)         # 確率密度関数
bar!(d, func=cdf, alpha=0.3)  # 分布
```


![image.png](https://qiita-image-store.s3.amazonaws.com/0/61915/cff12116-191a-a079-0059-bc39d3403903.png)




##１－５．Juliaで解く例題(2)

「確率統計」（森北出版）にある例題です。

**--- 問題 ---**

**ある機械がx時間以内に故障する確率F(x)が以下の式で与えられているとする。**
$$F(x)=1-e^{-\frac{x}{100}} \qquad (x>0)$$

１．機械が故障するまでの時間X従う分布と確率密度関数を求めよ。
２．機械が故障するまでの時間の平均E[X]と分散V[X]を求めよ。
３．この機械が故障するまでの時間が200時間を超える確率pを求めよ。


**--- 解答 ---**


```math
\begin{align}
&１．まず分布の定義と問題の条件より、\\
&F(x) = p(\{X \leqq x\})= 1-e^{-\frac{x}{100}} \qquad (x>0)\\
\\
&ゆえに確率密度関数f(x)は以下のようになる。\\
&f(x) = \frac{d \; F(x)}{dx} = \frac{1}{100} e^{-\frac{x}{100}}\\
\\
&つまり\lambda = \frac{1}{100}として、Xは指数分布E_x(\frac{1}{100})に従うと言える。\\
\\
\\
\\
&２．\lambda = \frac{1}{100}だから、\\
&E[X] = \frac {1} {\lambda} = 100\\
&V[X] = \frac {1} {\lambda^2} = 10000\\
\\
\\
\\
&３．200時間までに機械が故障する確率はF(200)で表される。\\
&求める確率は余事象の確率なので、p=1-F(200)となる。\\
\\
\\
&これ以上の具体的な計算は、Juliaに任せます。\\
&Q.E.D.\\
\\
&\qquad \qquad  \qquad \qquad  \qquad \qquad  \qquad \qquad \qquad  \qquad \qquad \qquad \qquad \qquad \qquad  \qquad \qquad \\
\\

\end{align}
```


**--- Juliaの計算 ---**



```julia
using Plots
using Distributions
using StatPlots

d=Exponential(100)
plot(d,fill=(0, .5,:orange))

println("F(200) = ", cdf(d,200))
println("1 - F(200) = ", 1-cdf(d,200))

plot(d,fill=(0, .5,:orange),xlabel = "time")
```

出力は以下のようになります。

```
F(200) = 0.8646647167633873         # 200時間までに機械が故障する確率
1 - F(200) = 0.1353352832366127     # 故障するまでの時間が200時間を超える確率
```

![image.png](https://qiita-image-store.s3.amazonaws.com/0/61915/6bcfdafe-7925-88ca-6bb0-a1d6416d81d8.png)



**確率密度関数と分布**の関係をグラフで表します。

```
d=Exponential(100)
scatter(d, leg=false)         # 確率密度関数
bar!(d, func=cdf, alpha=0.3)  # 分布
```

![image.png](https://qiita-image-store.s3.amazonaws.com/0/61915/8473d622-9ce8-615d-5564-fd93ccb97d71.png)


今回は以上です。