python で正規分布の確率密度関数を描く

データ分析
Sponsored Link

Pythonで正規分布の確率密度関数を描いてみます。

データは気象庁の仙台 日平均気温の月平均値(℃)を使って正規分布の確率密度関数を描いてみます。

手順としては、気温データをエクセルに落とし、さらにpythonのDataFrameとして読み込むところから始めます。

プログラミングを勉強したい!! プログラミングに興味がある!!
そんなあなたへのおすすめ記事
無料体験あり、キャッシュバック有りプログラミングスクールおすすめはこちら

 

データの前処理

# 確率密度関を計算するための下準備
import pandas as pd
import numpy as np

Miyagi_temp = pd.read_excel("Miyagi_temperture.xlsx", encoding = 'shift-jis', index_col=0) # データセットの読み込み。1926年からの月ごとの仙台市の平均気温
Miyagi_tempよ

読み込んだデータを下に表示してあります。

データを見て気が付くことが何点かありますね。

まず1926年のデータが10月ま9月までかけていること。

そして、変な記号 ] がいくつかのコラムに入っているので、まともな分析ができません。

そこで、データの前処理をする必要があるので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

すると「3月、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)

以下がプログラムを実行したときに描かれるグラフです。

数式をプログラムした時とほぼ同じグラフが描かれました。

プログラミングを勉強したい!! プログラミングに興味がある!!
そんなあなたへのおすすめ記事
無料体験あり、キャッシュバック有りプログラミングスクールおすすめはこちら

 

関連記事↓↓↓

Pythonで確率の計算。for文で集合の中の要素を抽出
確率という言葉は「降水確率」みたいな形で見たことはあるけど、なんとなく感覚的にわかっているだけで、実際に何かというと何かよくわからないですよね。 自分も統計の本を買ってみましたが、確率の項で こんなのが出てくるだけで拒否反応が出てし...
Pythonを使って度数分布表とヒストグラムを書く
統計で一番最初に勉強するのが、度数分布表とヒストグラムです。 Pythonで度数分布表とヒストグラムを描いてみましょう。 度数分布表のためのデータは、例題で学ぶ統計入門での演習問題で使われている、過去40日間の救急車出動...
Sponsored Link

コメント

タイトルとURLをコピーしました