練習問題 - 確率論（Pythonで計量経済学）

19498 ワード
scipy 確率論 計量経済学 Python 統計学 Python テキストリンク

Introduction to Econometrics with Rの練習問題をPythonで解いてみます。

2.3 練習問題 - 確率論

1. サンプリング

宝くじで$ 49 $のうち$ 6 $ つのユニークな数字が抽選されます。

指示:
今週の当選番号を引きなさい。

import math

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

from scipy import integrate, stats

np.random.seed(seed=123)
np.random.randint(low=1, high=50, size=6)

array([46,  3, 29, 35, 39, 18])

2. 確率密度関数

以下の確率密度関数 (PDF)を用いて、ランダムな変数 $X$ を考えなさい。

$$f_X(x)=\frac{x}{4}e^{-x^2/8},\quad x\geq 0.$$

指示:
上記の確率密度関数を関数 f()として定義せよ。
定義した関数が、実際に確率密度関数であることを確認せよ。

def f(x):
    return x/4*math.exp(-x**2/8)

integrate.quad(f, 0, np.inf)

(1.0, 2.1730298600934144e-09)

3. 期待値と分散

この練習問題では、前の練習問題で考えたランダムな変数 $X$ の期待値と分散を計算しなければなりません。

前の練習問題の確率密度関数 f() は動作環境で利用可能なものとします。

指示:
積分すると $X$ の期待値になる適切な関数 ex() を定義せよ。
$X$ の期待値を計算せよ。結果を expected_value に格納せよ。
積分すると $X^2$ の期待値になる適切な関数 ex2() を定義せよ。
$X$ の分散を計算せよ。結果を variance に格納せよ。

# define the function ex
def ex(x):
    return x*f(x)

# compute the expected value of X
expected_value = integrate.quad(ex, 0, np.inf)[0]

# define the function ex2
def ex2(x):
    return x**2*f(x)

# compute the variance of X
variance = integrate.quad(ex2, 0, np.inf)[0] - expected_value**2

4. 標準正規分布 I

Let $Z\sim\mathcal{N}(0, 1)$.

指示:
$\phi(3)$, つまり、$c=3$ での標準正規密度の値を計算せよ。

stats.norm.pdf(3)

0.0044318484119380075

5. 標準正規分布 II

Let $Z\sim\mathcal{N}(0, 1)$.

指示:
$P(|Z|\leq 1.64)$ を計算せよ。

# compute the probability
stats.norm.cdf(1.64) - stats.norm.cdf(-1.64)

0.8989948330517925

6. 正規分布 I

Let $Y\sim\mathcal{N}(5, 25)$.

指示:
与えられた分布の99%分位を計算せよ、すなわち、$y$ が $\Phi(\frac{y-5}{5})=0.99$ となるような $y$ を見つけよ。

# compute the 99% quantile of a normal distribution with mu = 5 and sigma^2 = 25.
stats.norm.ppf(0.99, 5, np.sqrt(25))

16.631739370204205

7. 正規分布 II

Let $Y\sim\mathcal{N}(2, 12)$.

指示:
この分布から$10$の乱数を生成せよ。

# generate 10 random numbers from the given distribution.
stats.norm.rvs(loc=2, scale=np.sqrt(12), size=10, random_state=12)

array([ 3.63847098, -0.36052849,  2.83983505, -3.89152106,  4.60896331,
       -3.31643067,  2.01776072,  1.58351913, -0.79546723, 11.9482742 ])

8. カイ二乗分布 I

Let $W\sim\chi^2_{10}$.

指示:
対応する確率密度関数をプロットせよ。x値の範囲を$[0,25]$に指定せよ。

# plot the PDF of a chi^2 random variable with df = 10

x = np.arange(0, 25)
plt.plot(x, stats.chi2.pdf(x, df=10))
plt.show()

9. カイ二乗分布 II

Let $X_1$ and $X_2$ be two independent normally distributed random variables with $\mu=0$ and $\sigma^2=15$.

指示:
$P(X_1^2+X_2^2>10)$ を計算せよ。

# compute the probability

stats.chi2.sf(10/15, df=2, loc=0, scale=1)

0.7165313105737892

10. スチューデントのt分布 I

Let $X\sim t_{10000}$ and $Z\sim\mathcal{N}(0,1)$.

指示:
両方の分布の $95%$分位を計算せよ。何か発見があるか？

# compute the 95% quantile of a t distribution with 10000 degrees of freedom
# qt(0.95, df = 10000)

print(stats.t.ppf(0.95, df = 10000))

# compute the 95% quantile of a standard normal distribution

print(stats.norm.ppf(0.95))

# both values are very close to each other. This is not surprising as for sufficient large degrees of freedom the t distribution can be approximated by the standard normal distribution.

1.6450060180692423
1.6448536269514722

11. スチューデントのt分布 II

Let $X\sim t_1$. Once the session has initialized you will see the plot of the corresponding probability density function (PDF).

指示:
この分布から$1000$の乱数を生成し、変数xに代入せよ。
xの標本平均を計算せよ。結果を説明でるか？

# generate 1000 random numbers from the given distribution. Assign them to the variable x.
x = stats.t.rvs(df = 1, size = 1000, random_state = 1)

# compute the sample mean of x.
np.mean(x)

# Although a t distribution with M = 1 is, as every other t distribution, symmetric around zero it actually has no expectation. This explains the highly non-zero value for the sample mean.

10.845661965991818

12. F分布 I

Let $Y\sim F(10, 4)$.

指示:
与えられた分布の分位関数をプロットせよ。

# plot the quantile function of the given distribution

dfn = 10
dfd = 4

x = np.linspace(stats.f.ppf(0.01, dfn, dfd),
                stats.f.ppf(0.99, dfn, dfd), 100)

plt.plot(stats.f.pdf(x = x, dfn = dfn, dfd = dfd))
plt.show()

13. F分布 II

Let $Y\sim F(4,5)$.

指示:
確率密度関数を積分して $P(1<Y<10)$ を計算せよ。

# compute the probability by integration

dfn = 4
dfd = 5

x = np.linspace(stats.f.ppf(0.01, dfn, dfd),
                stats.f.ppf(0.99, dfn, dfd), 100)

def function(x):
    return stats.f.pdf(x = x, dfn = dfn, dfd = dfd)

integrate.quad(function, 1, 10)[0]

0.4723970230052129