#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回/演習]]