#author("2023-10-04T17:58:16+09:00","default:cmdsadmin","cmdsadmin")
[[第1回]]
#author("2024-05-07T10:35:46+09:00","default:cmdsadmin","cmdsadmin")

* 第1回ハンズオン [#e0d61039]

#contents

** 0. 準備 [#fd9a7c5f]
- [[Google Drive:https://drive.google.com/]]にアクセスして,Google Colabを新規作成してください
- 名前を ''handson1.ipynb'' とします
- 下記のコードを貼り付けて実行してください
 # PandasとNumpyをインポート
 import pandas as pd
 import numpy as np
 
 # 日本語化Matplotlibもインポート
 import matplotlib.pyplot as plt
 !pip install japanize-matplotlib
 import japanize_matplotlib
 
** 1. Pythonで確率分布を見てみよう [#eafee07a]
*** SciPy statsライブラリ [#rc31d39d]
- Pythonには,[[SciPy:https://ja.wikipedia.org/wiki/SciPy]]という強力な数値解析ライブラリがある.その中のstatパッケージを使えば,
様々な統計関数を簡単に利用できる
-- 説明書:https://docs.scipy.org/doc/scipy/reference/stats.html

*** statsライブラリの基本的な使い方 [#nbddff14]
- まずはimportする
 from scipy import stats
- ''確率分布を定義''する
 変数名 = stats.分布(パラメータ)
- 離散型確率分布
-- ベルヌーイ分布: stats.bernoulli(p=確率)
-- 二項分布: stats.binom(n=個数, p=確率)
-- ポワソン分布: stats.poisson(mu=平均回数) 
- 連続型確率分布
-- 正規分布: stats.norm(loc=平均, scale=標準偏差)
-- 指数分布: stats.expon(scale=平均時間)
- 確率分布から''確率関数'' f_X(x)=P(X=x)の値を求める
 変数名.pmf(xの値)  #離散型の場合 probability mass function
 変数名.pdf(xの値)  #連続型の場合 probability density function
- 確率分布から''分布関数'' F_X(x)=P(X<=x))を求める
 変数名.cdf(xの値) #cumulative distribution function
- 確率分布から''分布関数'' F_X(x)=P(X>=x))を求める
 変数名.sf(xの値)  #survival function.1-変数名.cdf(x)に同じ
- 確率分布の''下側パーセント点''を求める
 変数名.ppf(下側確率)  #percent point function
- 確率分布の''上側パーセント点''を求める
 変数名.isf(上側確率)  #inverse survival function
- 確率分布に従った''乱数''を発生させる
 変数名.rvs(size=個数) #random variates

*** 正規分布を可視化してみよう [#j6d4ddda]
- scipy.statsを使って正規分布を可視化してみよう
-- 平均μ(loc)と標準偏差σ(scale)を変えて分布を見てみる

 from scipy import stats
 
 #平均100,標準偏差30の正規分布を定義
 dis1 = stats.norm(loc=100,scale=30)
 #平均130,標準偏差30の正規分布を定義
 dis2 = stats.norm(loc=130,scale=30)
 #平均100,標準偏差10の正規分布を定義
 dis3 = stats.norm(loc=100,scale=10)
 
 #xの範囲を適当に決める
 x=range(0,200)
 
 #確率密度関数をxにそってプロット
 plt.plot(x,dis1.pdf(x), label="N(100,30)")
 plt.plot(x,dis2.pdf(x), label="N(130,30)")
 plt.plot(x,dis3.pdf(x), label="N(100,10)")
 plt.grid()
 plt.legend()
 plt.show()
 
 #1000個の乱数を作成
 nums1 = dis1.rvs(size=1000)
 nums2 = dis2.rvs(size=1000)
 nums3 = dis3.rvs(size=1000)
 #ヒストグラムを作成
 plt.hist(nums1, bins=range(0,200,5), alpha=0.5)
 plt.hist(nums2, bins=range(0,200,5), alpha=0.5)
 plt.hist(nums3, bins=range(0,200,5), alpha=0.5)
 plt.grid()
 plt.show()

- Q: 正規分布N(100,30),N(130,30),N(100,10)からそれぞれ1つの数字をランダムに選んだ場合,
その数字が50以上75以下である確率を求めなさい
 #それぞれ分布関数から求める
 print(dis1.cdf(75)-dis1.cdf(50))
 print(dis2.cdf(75)-dis2.cdf(50))
 print(dis3.cdf(75)-dis3.cdf(50))

- Q: 正規分布N(100,30),N(130,30),N(100,10)について,それぞれ下側95%点を求めなさい
 #それぞれ分布関数から求める
 print(dis1.ppf(0.95))
 print(dis2.ppf(0.95))
 print(dis3.ppf(0.95))

*** 例題:神戸市で交通事故が起こる確率を調べる [#xe127cb7]

神戸市における交通事故は1か月約350件と報告されている.この時,ポワソン分布を使って,以下の問いに答えなさい.
なお,1日あたりの交通事故の平均件数をλ=350/30として計算せよ.

- RQ1: 1日あたり,ちょうど10件の交通事故が発生する確率はいくらか?
- RQ2: 1日あたり,5件以内に収まる確率はいくらか?
- RQ3: 1日あたり,10件以上~15件以内に収まる確率はいくらか?

 #ポワソン分布を作成
 kobe_acc = stats.poisson(mu=350/30)

 #概形をプロットしてみる
 x = range(1,30)
 plt.bar(x, kobe_acc.pmf(x))
 plt.grid()

 #確率関数,分布関数,パーセント点をそのまま表示するだけ
 print(f"[RQ1]: {kobe_acc.pmf(10)}")
 print(f"[RQ2]: {kobe_acc.cdf(5)}")
 print(f"[RQ3]: {kobe_acc.cdf(15) - kobe_acc.cdf(9)}")


** 2. 統計的推測をPythonでシミュレーションしてみよう [#q0b0c9ff]

*** 例題:全国17歳男子の身長の分布を推定する [#b6bcf85b]

[[ある調査:https://ai-trend.jp/basic-study/normal-distribution/example/]]によると,全国610,122人の17歳男子の平均身長は,

- 平均 μ = 170.6931cm
- 分散 σ^2=33.735 (標準偏差σ=5.8082)

の正規分布にほぼ従うという.

この時,母集団からn個の標本を取り出し,標本平均X~,不偏標本分散U^2を計算して,μ,σ^2を推定できるかやってみよう.

*** 準備 [#y54d0d48]
+ 母集団を表す正規分布を作成する
 #stats.normで全男子の母集団を生成(σ^2でなくσを使うことに注意)
 all_boys = stats.norm(loc=170.6931,scale=5.8082)

*** 平均μ(=170.6931)を推定してみよう [#sbd1430b]
+ all_boysからn=10人の標本X={x1,x2,...,x10}を取り出し,標本平均X~を計算しよう
 #nをとりあえず10として
 n=10
 #ランダムにn人取り出してみる
 sample_boys = all_boys.rvs(size=n)
 #標本平均を計算する(何回も試してみる)
 sample_boys.mean()
+ この作業を1000回繰り返した結果を表に入れてみよう
 #表を作る
 df = pd.DataFrame()
 #1000回繰り返して表に入れる(リストの内包表現)
 df[f"{n}人の標本平均"] = [all_boys.rvs(size=n).mean() for i in range(1000)]
 df
+ 標本平均の平均E(X~)を計算してみる
 df.mean()
+ nを変化させて同じ実験をする
 for n in [10,100,1000,10000]:
   df[f"{n}人の標本平均"] = [all_boys.rvs(size=n).mean() for i in range(1000)]
 df
+ 標本平均の平均E(X~)が平均的にμの近くの値になっていることがわかる(''X~がμの不偏推定量である'')
 df.mean()
+ ヒストグラムを描いてみる.nが大きくなるほどX~の分布がμに近づくことがわかる (''X~がμの一致推定量である'')
 df.plot.hist(subplots=True, bins=100)

*** σ^2(=33.735)を推定してみよう [#z6433da4]
+ μの時とほぼ同じ.numpyの標本偏分散U^2は,var(ddof=1)で計算する
 #σ用の表を作る
 df2 = pd.DataFrame()
 #nを変化させてn個の不偏標本分散を1000回取る実験をする(不偏推定量の確認)
 for n in [10,100,1000,10000]:
   df2[f"{n}人の不偏標本分散"] = [all_boys.rvs(size=n).var(ddof=1) for i in range(1000)]
 #不偏標本分散の平均E(U^2) (U^2がσ^2の不偏推定量である)
 df2.mean()
+ ヒストグラムを描いてみる(U^2がσ^2の一致推定量である)
 df2.plot.hist(subplots=True, bins=100)

** 第1回演習へ [#f6c13d5f]

- [[第1回/演習]]


トップ   編集 差分 履歴 添付 複製 名前変更 リロード   新規 一覧 検索 最終更新   ヘルプ   最終更新のRSS