Pythonで正規分布の確率密度関数を描いてみます。
データは気象庁の仙台 日平均気温の月平均値（℃）を使って正規分布の確率密度関数を描いてみます。
この記事はこんな人におすすめ。
- データの前処理の方法を知りたい。
- 正規分布の確率密度関数を描いてみたい。
手順としては、気温データをエクセルに落とし、さらにpythonのDataFrameとして読み込むところから始めます。
プログラミング無料体験はこちら↓↓↓
データの前処理
# 確率密度関を計算するための下準備
import pandas as pd
import numpy as np
# データセットの読み込み。1926年からの月ごとの仙台市の平均気温
Miyagi_temp = pd.read_excel("Miyagi_temperture.xlsx", encoding = 'shift-jis', index_col=0)
Miyagi_temp
読み込んだデータを下に表示してあります。
データを見て気が付くことが何点かありますね。
まず1926年のデータが10月ま９月までかけていること。
そして、変な記号 ］ がいくつかのコラムに入っているので、まともな分析ができません。
そこで、データの前処理をする必要があるのでdtype()関数を使い、データの型を確認します。
print(Miyagi_temp.dtypes) # dtype()関数を使ってデータの型を確認。
1月 float64
2月 float64
3月 object
4月 float64
5月 float64
6月 float64
7月 float64
8月 float64
9月 float64
10月 float64
11月 object
12月 float64
年の値 object
dtype: object
すると「３月、11月、年の値」のコラムで、データタイプがobject型になっているのがわかります。
これらのコラムの余計な記号を削除し、数値をすべてfloat型にする必要があります。
まずは自力で何とかしようと数時間格闘しましたが、結局ダメ。
そこで、スタックオーバーフローに質問を投稿したところ、metropolisさんから教えてもらったのが、以下のコード。
# objet typeの列名を抽出
cols = Miyagi_temp.select_dtypes('O').columns
# 数字と負号('-')と '.' 以外を削除してdtypeをobjectからfloat64へ変更
Miyagi_temp[cols] = (Miyagi_temp[cols].replace(r'[^\d.-]', '', regex=True).astype('float64'))
pd.set_option('display.unicode.east_asian_width', True)
print(Miyagi_temp.dtypes)
上のコードを実行するとすべてのコラムがfloat型になりました。
1月 float64
2月 float64
3月 float64
4月 float64
5月 float64
6月 float64
7月 float64
8月 float64
9月 float64
10月 float64
11月 float64
12月 float64
年の値 float64
dtype: object
では、肝心の記号が消されているかどうかを確認してみましょう。
print(Miyagi_temp.head().to_string(index=False))
1月 2月 3月 4月 5月 6月 7月 8月 9月 10月 11月 12月 年の値
NaN NaN NaN NaN NaN NaN NaN NaN NaN 11.8 6.8 1.7 6.8
-0.6 -1.5 2.5 9.0 13.6 17.6 23.4 24.2 18.5 13.9 8.3 2.1 10.9
0.4 0.1 2.9 8.8 13.7 16.9 21.1 22.4 21.6 14.2 8.8 1.4 11.0
-1.3 -0.7 2.9 8.7 12.9 16.7 23.3 24.9 18.7 14.3 8.2 5.0 11.1
-0.6 1.9 5.6 9.9 14.4 18.4 22.4 24.5 19.2 13.9 6.8 2.3 11.6
記号がしっかりと消えてくれています。
何時間もかかった問題が解決するのはうれしいものですね。
では次に、データが入っていないセルがある行の処理をしましょう。
今更数字で埋め合わせることはできず、分析結果への影響はあまりないと思いますので、この行そのものを削除してしまいます。
Miyagi_temp = Miyagi_temp.dropna(how='any') # 分析の邪魔になるので、値が入っていない行を削除する。
Miyagi_temp.head()
コードを実行すると、1926年の行が削除されました。
正規分布の確率密度関数を描く
データの前処理が終わったので、11月のデータを使って正規分布の確率密度関数を描いてみます。
平均、分散
の正規分布
に従う確率変数
の確率密度関数は確率密度関数がこれ。
この式をプログラムで表す場合に必要なのが、
相加平均 = mu =
標準偏差 = sigma =
の値。
以下のプログラムを実行して、平均と標準偏差を求めてみました。
# 相加平均と標準偏差のDataFrameを作成する。
mu = Miyagi_temp["11月"].mean() # 11月の気温の相加平均を求める。
sigma = Miyagi_temp["11月"].std()# 11月の気温の標準偏差を求める。
pd.DataFrame([[mu, sigma]], columns =['相加平均', '標準偏差']) # pd.DataFrame()でDataFrameを作成する。
| 相加平均 | 標準偏差 | |
|---|---|---|
| 0 | 8.896809 | 1.165125 |
再びのエラー発生
上のコードを実行した際にpd.DataFrameを使ってDataFrameを作ったんですが、こんなエラーが発生。
Type Error:'list' object is not callable
「またエラーかよ！！！」思いつついろいろ調べてみると、
pd.DataFrame([[mu, sigma]],
の部分でエラー発生。
このエラー、listが変数として定義されている場合や、関数に格納される場合に発生します。
そこで、mu と sigmaのデータタイプをdtype()で調べると、なぜかリストになっているんですね。
解決方法をググってみると「kernelをリスタートせよ」と書いてあるので、いったん閉じてリスタートすると問題解決しました。
なんでそうなるのかよく分からん。
気を取り直して、正規分布の確率密度関数をpythonに当てはめていきます。
def normal_dist_pdf(x, mu, sigma):
# ガウスの公式
y = (1/np.sqrt(2* np.pi*sigma**2)) * np.exp(-1 * (x-mu)**2/(2*sigma**2))
return y
次に、データ全体の傾向を見てみます。
pd.DataFrame(Miyagi_temp['11月'].describe())
プログラムを実行してみましょう。
# 11月データのヒストグラムを描く。
ax = Miyagi_temp['11月'].plot(kind = 'hist', bins=94)
# 最小値6.5 から 最大値11.8 までを対象（今回のデータの範囲）
x = np.arange(6.5, 11.8)
# 今回作った関数でラインを引く
y = normal_dist_pdf(x, mu, sigma)
ax2 = ax.twinx()
ax2.plot(x, y, color = 'r', linewidth = 5.0)
このプログラムを実行して描かれたのが、以下のプロット。
仙台の11月の平均気温は過去94年間を通じて、概ね正規分布に収まるといえますね。
正規分布の確率密度関数ですが、pythonのliblary の一つscipyを使うことで余計なコードを描かずに描くことができます。
from scipy.stats import norm # ライブラリ scipy から normをインポート
mu = Miyagi_temp["11月"].mean() # 平均を求める
sigma = Miyagi_temp["11月"].std() # 標準偏差を求める
r = np.arange(6.5, 11.8) # x 軸の範囲
pdf = norm.pdf(r, loc=mu, scale=sigma) # pdf (probability density function) を実行
# グラフによる描写
ax = Miyagi_temp["11月"].plot(kind='hist', bins=94)
ax2 = ax.twinx() ax2.plot(r, pdf, color = 'r', linewidth = 5.0)
以下がプログラムを実行したときに描かれるグラフです。
数式をプログラムした時とほぼ同じグラフが描かれました。
プログラミング無料体験はこちら↓↓↓
コメント