Pythonで学ぶ統計学１：記述統計(データ・グラフ）

[IN]
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import japanize_matplotlib
import os
import glob

class HorizontalDisplay:
    def __init__(self, *args):
        self.args = args
    
    def _repr_html_(self):
        template = '<div style="float: left; padding: 10px;">{0}</div>'
        return "\n".join(template.format(a._repr_html_()) for a in self.args)
        

# データの読み込み
filepath = 'data/FEH_00450171_230101142238.csv'
df = pd.read_csv(filepath,
             encoding='shift-jis',
             names=[f'{i}' for i in range(100)])

# 前処理
df = df.dropna(axis=1,thresh=10) # 欠損値が10個以上ある列を削除
df = df.dropna(axis=0,thresh=10) # 欠損値が10個以上ある行を削除
df.reset_index(drop=True,inplace=True)
df = df.fillna('') # 欠損値を空白に置換

#カラム名の設定
cols = []
for c1, c2, c3 in zip(df.loc[5], df.loc[3], df.loc[1]):
    if c2 == '':
        c1+=c2
    else:
        c1 = c1+'_'+c2
    if c3 == '':
        c1 += c3
    else:
        c1 = c1+'_'+c3
    cols.append(c1)

df = df.iloc[6:, :]
df.columns = cols

#不要なインデックスを削除
def func(x):
    if '総数' in x or '再掲' in x:
        return False
    else:
        return True
df = df[df['年齢階級'].apply(func)]

#文字列型を数値型に変換するためのDict
dict_obj2float = {'人数【人】_身長(cm)_男性': int, '平均値_身長(cm)_男性': float, '標準偏差_身長(cm)_男性': float, '平均値_体重(kg)_男性': float, '標準偏差_体重(kg)_男性': float, 
                  '人数【人】_身長(cm)_女性': int, '平均値_身長(cm)_女性': float, '標準偏差_身長(cm)_女性': float, '平均値_体重(kg)_女性': float, '標準偏差_体重(kg)_女性': float}


display(df.head())
df_2018 = df[df['時間軸(年次)']=='2018年']
df_2018 = df_2018.astype(dict_obj2float)
df_2018_main = df_2018[['時間軸(年次)', '年齢階級', '平均値_身長(cm)_男性', '平均値_身長(cm)_女性']]

mean_men, mean_women = df_2018_main['平均値_身長(cm)_男性'].mean(), df_2018_main['平均値_身長(cm)_女性'].mean()

#可視化
plt.figure(figsize=(10, 5))
plt.bar(df_2018_main['年齢階級'], df_2018_main['平均値_身長(cm)_男性'], label='男性', 
        align='edge', width=-0.5, color='royalblue')
plt.bar(df_2018_main['年齢階級'], df_2018_main['平均値_身長(cm)_女性'], label='女性',
        align='edge', width=0.5, color='tomato')

plt.plot(df_2018_main['年齢階級'], np.ones(len(df_2018_main))*mean_men, label='男性平均',
         color='blue', linestyle='--')
plt.plot(df_2018_main['年齢階級'], np.ones(len(df_2018_main))*mean_women, label='女性平均',
         color = 'red', linestyle='--')

plt.xticks(rotation=90)
plt.grid()
plt.legend()
plt.show()

[OUT]

[IN※前処理]
import os
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import re

url = 'https://www.gurazeni.com/'
df_first = pd.read_html(url)[0]
df_master = df_first.copy()

for page in range(2, 12):
    url_next = f'https://www.gurazeni.com/ranking/{page}'
    df_next = pd.read_html(url_next)[0]
    if len(df_next) == 0:
        break
    else:
        df_master = pd.concat([df_master, df_next], axis=0)

df_master.reset_index(drop=True, inplace=True)
print(df_master.shape)
display(df_master)

 #年俸 (推定)の文字列をParseして数値に変換する関数
income_str = df_master['年俸(推定)']

def parse_income(income_str):
    # 正規表現パターンによる数値抽出
    pattern = r'([0-9,]+)' #0 ~9の数字とカンマを抽出
    match = re.search(pattern, income_str) #文字列の中からパターンにマッチする最初の部分を抽出 
    
    # 正規表現にマッチしない場合は "None" を返す
    if not match:
        return None
    
    # 数値文字列から数値に変換※カンマは存在しないけど念のため
    income = int(match.group(1).replace(",", ""))
    
    # 単位を取得
    if "億" in income_str:
        income *= 10000 * 10000
    elif "万" in income_str:
        income *= 10000
    
    return income

df_income = income_str.apply(parse_income)
display(df_income)

[OUT]
(1042, 8)

[IN]
df = df_income.dropna() #欠損値 (調査中(情報募集中))を削除
 #カンマ区切りの数値を表示するための設定 
pd.options.display.float_format = '{:,.0f}'.format
df.describe()

[OUT]
count         1,040
mean     37,648,250
std      72,838,244
min       2,200,000
25%       5,685,000
50%      12,000,000
75%      36,000,000
max     900,000,000
Name: 年俸(推定), dtype: float64

[IN]
import matplotlib.pyplot as plt #正規分布を作成  
import japanize_matplotlib

def conv2yen(value):
    if value >= 100000000:
        return f"{int(value / 100000000)}億円"
    elif value >= 10000:
        return f"{int(value / 10000)}万円"
    else:
        return f'{int(value)}円'
 #ヒストグラムを作成 
num_bins = 1 + np.log2(len(df)) #スタージェスの公式によるビンの数を計算 
print(f'bins: {int(num_bins)}, 3.0倍: {int(3.0*num_bins)}')
counts, _bins, _ = plt.hist(df, bins=int(3.0*num_bins), color='red')
bins = _bins[:-1] + (_bins[1] - _bins[0]) / 2 #ビンの中心を計算 

plt.xticks([i*0.5e8 for i in range(19)], rotation=90), plt.xlabel('年棒'), plt.ylabel('人数')
plt.gca().get_xaxis().get_major_formatter().set_scientific(False)#x軸の目盛り表示を1億単位にする #アノテーションを追加 （矢印あり）
mean, median, maxdata = conv2yen(df.mean()), conv2yen(df.median()), conv2yen(df.max())
plt.annotate(f'平均値: {mean}', 
             xy=(df.mean(), 100), 
             xytext=(df.mean(), 0.5*max(counts)), 
             color='blue',
             arrowprops=dict(facecolor='blue', shrink=0.02))

plt.annotate(f'中央値: {median}',
             xy=(df.median(), 100),
             xytext=(df.median(), max(counts)),
             color='green',
             arrowprops=dict(facecolor='green', shrink=0.02))

plt.annotate(f'最大値: {maxdata}',
             xy=(df.max(), 100),
             xytext=(df.max(), 0.5*max(counts)),
             color='red',
             arrowprops=dict(facecolor='red', shrink=0.02))
plt.grid()
plt.show()

[OUT]
bins: 11, 3.0倍: 33

[IN]
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

np.set_printoptions(precision=3, suppress=True) #小数点以下3桁表示、指数表示しない
np.random.seed(123)#乱数の種を設定
nums_sample = 10 #サンプル数

target_value = 10 #真値(目標値)とする値
x = target_value + np.random.randn(nums_sample) #真値にばらつきを追加
# print(np.sort(x))
x_real = np.full_like(x, target_value) #真値を格納した配列
x_mean = np.full_like(x, x.mean()) #平均値を格納した配列
x_median = np.full_like(x, np.median(x)) #中央値を格納した配列

#可視化
plt.figure(figsize=(10, 6), facecolor='w')
plt.scatter(range(nums_sample), x, label='Data', color='blue', marker='o') #サンプルデータ
plt.plot(x_real, label='Target', color='red', linestyle='dashed') #真値
plt.plot(x_mean, label='Mean', color='green', linestyle='dashed') #平均値
plt.plot(x_median, label='Median', color='orange', linestyle='dashed') #中央値
plt.grid(), plt.legend(), plt.xlabel('サンプルデータ数'), plt.ylabel('計測値')
plt.show()


[OUT]
[ 7.573  8.494  8.914  9.133  9.421  9.571 10.283 10.997 11.266 11.651]

[IN]
print(np.median(x), pd.DataFrame(x).median())

plt.figure(figsize=(10, 6), facecolor='w')
plt.bar(np.arange(nums_sample), np.sort(x)) #サンプルデータ
plt.xticks(np.arange(nums_sample), np.arange(nums_sample)), plt.xlabel('サンプル数') #x軸の目盛り
plt.text(4.5, 10.5, '中央値', color='red') #テキスト追加
plt.arrow(5, 10.3, 0, -0.5, color='red', head_width=0.2, head_length=0.1) #矢印追加
plt.show()

[OUT]
9.4962435  9.496244

[IN]
a_1, a_2, a_3 = np.quantile(x, 0.25), np.quantile(x, 0.5), np.quantile(x, 0.75) #第1四分位数、第2四分位数、第3四分位数
print(f'第1四分位数:{a_1:.3f}, 第2四分位数:{a_2:.3f}, 第3四分位数:{a_3:.3f}')

a_1p, a_2p, a_3p = np.percentile(x, 25), np.percentile(x, 50), np.percentile(x, 75) #第1四分位数、第2四分位数、第3四分位数
print(f'第1四分位数:{a_1p:.3f}, 第2四分位数:{a_2p:.3f}, 第3四分位数:{a_3p:.3f}')

df = pd.DataFrame(x, columns=['Data'])
a_1df, a_2df, a_3df = df.quantile(0.25), df.quantile(0.5), df.quantile(0.75) #第1四分位数、第2四分位数、第3四分位数
print(f'第1四分位数:{a_1df[0]:.3f}, 第2四分位数:{a_2df[0]:.3f}, 第3四分位数:{a_3df[0]:.3f}')


[OUT] ※3つとも同じ出力
第1四分位数:8.969, 第2四分位数:9.496, 第3四分位数:10.819

[IN]
#データ作成
clothes_sizes = np.random.choice(['S','M','L','XL'], 20)
print(clothes_sizes)

#Numpy
uniqs, count = np.unique(clothes_sizes, return_counts=True) #要素の種類を取得
idx = np.argmax(count)
print(f'最頻値:{uniqs[idx]}, 数:{count[idx]}')
#Pandas
df_clothes = pd.Series(clothes_sizes)
df_clothes.mode()[0], df_clothes.value_counts()['L']

[OUT]
['L' 'M' 'L' 'L' 'S' 'L' 'L' 'M' 'XL' 'L' 'XL' 'M' 'L' 'M' 'S' 'M' 'L'
 'XL' 'M' 'S']

最頻値:L, 数:8
('L', 8)

[IN]
x = np.random.randn(100)
mean, var, std = x.mean(), x.var(), x.std()
print(f'Numpy　平均：{mean:.3f}, 分散: {var:.3f}, 標準偏差: {std:.3f}')

mean, var, std = x.mean(), x.var(ddof=1), x.std(ddof=1)
print(f'Numpy　平均：{mean:.3f}, 分散: {var:.3f}, 標準偏差: {std:.3f}')

df_x = pd.DataFrame(x, columns=['x'])
mean, var, std = df_x['x'].mean(), df_x['x'].var(), df_x['x'].std()
print(f'Pandas　平均：{mean:.3f}, 分散: {var:.3f}, 標準偏差: {std:.3f}')

[OUT]
Numpy　平均：0.027, 分散: 1.273, 標準偏差: 1.128
Numpy　平均：0.027, 分散: 1.286, 標準偏差: 1.134
Pandas　平均：0.027, 分散: 1.286, 標準偏差: 1.134

[IN]
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import japanize_matplotlib

def linear(x, a, b):
    return a * x + b

def getminmax(x):
    min, max = x.min(), x.max()
    if np.abs(min) > np.abs(max):
        return np.abs(min)
    else:
        return np.abs(max)

# データの作成
np.random.seed(1233)
x = np.random.randn(20) * 10
y = linear(x, 2, 3) 
 #x軸とy軸の最大値を取得 （軸設定）
xtick = getminmax(x)
ytick = getminmax(y)
 #グラフの描画 
fig, axs = plt.subplots(3, 2, figsize=(14, 18), facecolor='w', edgecolor='k')
axs[0][0].scatter(x, y, s=50, c='r', marker='o', label='data') #x軸とy軸を記載 
axs[0][0].plot([-xtick, xtick], [0,0], c='black', lw=1)
axs[0][0].plot([0,0], [-ytick, ytick], c='black', lw=1)
axs[0][0].set(title='y=ax+b(オリジナル)',
              xlabel='x', 
              ylabel='y')
axs[0][0].grid()
 #x軸を標準化 
x_mean, y_std = x.mean(), x.std()
x_std = (x - x_mean) / y_std #標準化 

axs[0][1].scatter(x_std, y, s=50, c='r', marker='o', label='data')
axs[0][1].plot([-xtick, xtick], [0,0], c='black', lw=1)
axs[0][1].plot([0,0], [-ytick, ytick], c='black', lw=1)
axs[0][1].set(title='y=ax+b(x軸を標準化)',
              xlabel='x(標準化)',
              ylabel='y')
axs[0][1].grid()
axs[0][1].annotate('xの平均0, 分散1へ圧縮', xy=(5, -10), xytext=(15, -10),
                   arrowprops=dict(facecolor='red', shrink=0.05))
axs[0][1].annotate('', xy=(-5, 10), xytext=(-15, 10),
                   arrowprops=dict(facecolor='red', shrink=0.05))
 #y軸を標準化 
y_mean, y_std = y.mean(), y.std()
y_std = (y - y_mean) / y_std #標準化 

axs[1][0].scatter(x, y_std, s=50, c='r', marker='o', label='data')
axs[1][0].plot([-xtick, xtick], [0,0], c='black', lw=1)
axs[1][0].plot([0,0], [-ytick, ytick], c='black', lw=1)
axs[1][0].set(title='y=ax+b(y軸を標準化)',
              xlabel='x',
              ylabel='y(標準化)')
axs[1][0].grid()
axs[1][0].annotate('yの平均0, 分散1へ圧縮', xy=(5, -10), xytext=(5, -30), 
                   arrowprops=dict(facecolor='red', shrink=0.05))
axs[1][0].annotate('yの平均0, 分散1へ圧縮', xy=(-5, 10),xytext=(-5, 30),
                   arrowprops=dict(facecolor='red', shrink=0.05))
 #y軸も標準化 
y_mean, y_std = y.mean(), y.std()
y_std = (y - y_mean) / y_std #標準化 

axs[1][1].scatter(x_std, y_std, s=50, c='r', marker='o', label='data')
axs[1][1].plot([-xtick, xtick], [0,0], c='black', lw=1)
axs[1][1].plot([0,0], [-ytick, ytick], c='black', lw=1)
axs[1][1].set(title='y=ax+b(x・y軸を標準化)',
              xlabel='x(標準化)',
              ylabel='y(標準化)')
axs[1][1].grid()

axs[1][1].annotate('xの平均0, 分散1へ圧縮', xy=(5, -10), xytext=(15, -10),
                   arrowprops=dict(facecolor='red', shrink=0.05))
axs[1][1].annotate('', xy=(-5, 10), xytext=(-15, 10),
                   arrowprops=dict(facecolor='red', shrink=0.05))
axs[1][1].annotate('yの平均0, 分散1へ圧縮', xy=(5, -10), xytext=(5, -30), 
                   arrowprops=dict(facecolor='red', shrink=0.05))
axs[1][1].annotate('yの平均0, 分散1へ圧縮', xy=(-5, 10),xytext=(-5, 30),
                   arrowprops=dict(facecolor='red', shrink=0.05))
 #x軸とy軸を標準化 (拡大)
xtick_std, ytick_std = getminmax(x_std), getminmax(y_std)

axs[2][0].scatter(x_std, y_std, s=50, c='r', marker='o', label='data')
axs[2][0].plot([-xtick_std, xtick_std], [0,0], c='black', lw=1)
axs[2][0].plot([0,0], [-ytick_std, ytick_std], c='black', lw=1)
axs[2][0].set(title='y=ax+b(x・y軸を標準化(拡大))',
              xlabel='x(標準化)',
              ylabel='y(標準化)')
axs[2][0].grid()

axs[2][0].text(0.5, 1.5, 'xyは正(x:+, y:+)', fontsize=14, color='red')
axs[2][0].text(-1.5, -1.5, 'xyは正(x:-, y:-)', fontsize=14, color='red')
axs[2][0].text(0.5, -1.5, 'xyは負(x:+, y:-)', fontsize=14, color='blue')
axs[2][0].text(-1.5, 1.5, 'xyは負(x:-, y:+)', fontsize=14, color='blue')
 #標準化したxとyの積を計算 
xy_std = x_std * y_std
xytick = getminmax(xy_std)
print(f'xyのデータ数:{len(xy_std)}')
print(f'相関係数:{xy_std.sum()/len(xy_std)}') #相関係数 =Cov(x,y)/σxσy=(1/n)Σ(x-μx)(y-μy)/σxσy=(1/n)Σ(標準化x)(標準化y)

axs[2][1].scatter(xy_std, np.zeros_like(xy_std), s=30, c='r', marker='x', label='data')
axs[2][1].plot([0, xytick], [0,0], c='black', lw=1)
axs[2][1].set(title='標準化x×標準化y',
              xlabel='x(標準化)×y(標準化)の1次元プロット')
axs[2][1].grid() #y軸を非表示 
axs[2][1].axes.yaxis.set_visible(False)


plt.show()

[OUT]
xyのデータ数:20
相関係数:1.0000000000000004

[IN]
def classify_xy(x, y):
    if x>=0 and y>=0:
        return 0
    elif x<0 and y<0:
        return 1
    elif x>=0 and y<0:
        return 2
    elif x<0 and y>=0:
        return 3
    else:
        None

x = np.linspace(-2*np.pi, 2*np.pi, 21)
y = np.sin(x)

fig, axs = plt.subplots(2, 2, figsize=(12, 12), facecolor='w', edgecolor='k')
 #Fig .1 y=sin(x)
axs[0][0].scatter(x, y, s=30, c='r', marker='o', label='data')
axs[0][0].plot(x, y, c='b', lw=0.2, ls='--') #確認用 
axs[0][0].plot([-2*np.pi, 2*np.pi], [0,0], c='black', lw=1)
axs[0][0].plot([0,0], [-1,1], c='black', lw=1)
axs[0][0].set(title='y=sin(x)', xlabel='x', ylabel='sin(x)')
axs[0][0].grid()
 #Fig .2 y=sin(x)をx軸・y軸ともに標準化
mean_x, std_x, mean_y, std_y = x.mean(), x.std(), y.mean(), y.std()
x_std, y_std = (x - mean_x) / std_x, (y - mean_y) / std_y
axs[0][1].scatter(x_std, y_std, s=30, c='r', marker='o', label='data')
axs[0][1].plot([-2*np.pi, 2*np.pi], [0,0], c='black', lw=1)
axs[0][1].plot([0,0], [-1,1], c='black', lw=1)
axs[0][1].set(title='y=sin(x)のx・yを標準化', xlabel='x(標準化)', ylabel='sin(x)(標準化)')
axs[0][1].grid()

axs[0][1].text(2.0, 1.0, 'xyは正(x:+, y:+)', fontsize=14, color='red')
axs[0][1].text(-6.0, -1.0, 'xyは正(x:-, y:-)', fontsize=14, color='red')
axs[0][1].text(2.0, -1.0, 'xyは負(x:+, y:-)', fontsize=14, color='blue')
axs[0][1].text(-6.0, 1.0, 'xyは負(x:-, y:+)', fontsize=14, color='blue')
 #Fig .3 色分けしたプロットを作成
classes = np.vectorize(classify_xy)(x_std, y_std)
df = pd.DataFrame({'x':x, 'y':y, 'x_std':x_std, 'y_std':y_std, 'class':classes})
df0, df1, df2, df3 = df[df['class']==0], df[df['class']==1], df[df['class']==2], df[df['class']==3]
axs[1][0].scatter(df0['x_std'], df0['y_std'], s=30, c='r', marker='o', label='label0:x:+, y:+')
axs[1][0].scatter(df1['x_std'], df1['y_std'], s=30, c='b', marker='^', label='label1:x:-, y:-')
axs[1][0].scatter(df2['x_std'], df2['y_std'], s=30, c='g', marker='s', label='label2:x:+, y:-')
axs[1][0].scatter(df3['x_std'], df3['y_std'], s=30, c='y', marker='x', label='label3:x:-, y:+')
axs[1][0].plot([-2*np.pi, 2*np.pi], [0,0], c='black', lw=1)
axs[1][0].plot([0,0], [-1,1], c='black', lw=1)
axs[1][0].set(title='y=sin(x)のx・yを標準化（色分けVer.）', xlabel='x(標準化)', ylabel='sin(x)(標準化)')
axs[1][0].grid(), axs[1][0].legend()
 #Fig .4 y=sin(x)をx軸・y軸ともに標準化したxとyの積を計算
df['xy'] = df['x_std'] * df['y_std']
df0, df1, df2, df3 = df[df['class']==0], df[df['class']==1], df[df['class']==2], df[df['class']==3]
axs[1][1].scatter(df0['xy'], df0['xy']*0, s=20, c='r', marker='o', label='label0:x:+, y:+')
axs[1][1].scatter(df1['xy'], df1['xy']*0, s=20, c='b', marker='^', label='label1:x:-, y:-')
axs[1][1].scatter(df2['xy'], df2['xy']*0, s=20, c='g', marker='s', label='label2:x:+, y:-')
axs[1][1].scatter(df3['xy'], df3['xy']*0, s=20, c='y', marker='x', label='label3:x:-, y:+')
axs[1][1].plot([-2, 2], [0,0], c='black', lw=1)
axs[1][1].set(title='y=sin(x)のx・yを標準化したxとyの積を計算', xlabel='x(標準化)×sin(x)(標準化)の1次元プロット')
axs[1][1].axes.yaxis.set_visible(False)
axs[1][1].grid(), axs[1][1].legend()
plt.show()

print(f'xyのデータ数:{len(xy)}')
print(f'相関係数:{df["xy"].sum()/len(df["xy"])}')  #相関係数 =Cov(x,y)/σxσy=(1/n)Σ(x-μx)(y-μy)/σxσy=(1/n)Σ(標準化x)(標準化y)

[OUT]
xyのデータ数:21
相関係数:-0.3507

[IN]
df_2018_men, df_2018_women = df_2018_main[['年齢階級', '平均値_身長(cm)_男性']], df_2018_main[['年齢階級', '平均値_身長(cm)_女性']]
display(HorizontalDisplay(df_2018_men.T, df_2018_women.T))
df_2018_main.corr()

[OUT]
　　　　　　	　平均値_身長(cm)_男性	平均値_身長(cm)_女性
平均値_身長(cm)_男性	1.000000	0.988219
平均値_身長(cm)_女性	0.988219	1.000000

[IN]
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import japanize_matplotlib

heights = np.array([167, 172, 188, 154, 180, 195, 160, 177, 185, 169, 192, 157, 174, 181, 159, 186, 175, 190, 166, 198])
print(f'MIN:{heights.min()}, MAX:{heights.max()}')
hist, bins = np.histogram(heights, bins=5)
print(f'度数:{hist}, bins(階級):{bins}')
relative_frequecy = hist/len(heights)
print(f'相対度数:{relative_frequecy}, 相対度数の合計：{relative_frequecy.sum()}')
cumulative_frequency = np.cumsum(relative_frequecy)
print(f'累積度数:{cumulative_frequency}')

bin_width = bins[1] - bins[0]
plt.bar(bins[:-1], hist, width=bin_width, align='edge', edgecolor='white')
for i, h in enumerate(hist):
    plt.text(bins[i]+bin_width/2, h, str(h), ha='center', va='bottom', fontsize=14)
 #累積度数 (Cumulative Frequency)
plt.plot(bins[:-1]+bin_width/2, cumulative_frequency, marker='o', color='white')
for i, c in enumerate(cumulative_frequency):
    plt.text(bins[i]+bin_width/2, c, str(c), ha='center', va='bottom', fontsize=12, color='white')


[OUT]
MIN:154, MAX:198
度数:[4 3 5 4 4], bins(階級):[154.  162.8 171.6 180.4 189.2 198. ]
相対度数:[0.2  0.15 0.25 0.2  0.2 ], 相対度数の合計：1.0
累積度数:[0.2  0.35 0.6  0.8  1.  ]

[numpy.histogram()のAPI確認]
numpy.histogram(a, bins=10, range=None, density=None, weights=None)

[IN]
#2018年の30-39歳の男性の平均身長 https://www.e-stat.go.jp/dbview?sid=0003224177
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
from typing import List
from scipy import stats
import os 

np.random.seed(12)

 #正規分布確認用 def norm(x, mu, sigma):
    return 1/(sigma*np.sqrt(2*np.pi))*np.exp(-(x-mu)**2/(2*sigma**2))

def showGraph(xs: List, ys: List, labels: List, verbose=False):
    fig = plt.figure(figsize=(10, 6), facecolor='w')
    ax = fig.add_subplot(111)
    
    colors = {0:'red', 1:'blue', 2:'green', 3:'orange'}
    idx = 0
     #確率密度分布     for x, y, label in zip(xs, ys, labels):
         #分布を表示         ax.plot(x, y, label=label, color=colors[idx])  
        if idx == 0:
            mean, std = 172.1, 6.1
            sigma1_low, sigma1_high = mean-std, mean+std
            sigma2_low, sigma2_high = mean-2*std, mean+2*std
            y_sigma1_low = stats.norm.pdf(sigma1_low, loc=mean, scale=std)  #正規分布             y_sigma1_high = stats.norm.pdf(sigma1_high, loc=mean, scale=std)  #正規分布             y_sigma2_low = stats.norm.pdf(sigma2_low, loc=mean, scale=std)  #正規分布             y_sigma2_high = stats.norm.pdf(sigma2_high, loc=mean, scale=std)  #正規分布             
             #標準偏差を表示             ax.plot([sigma1_low, sigma1_low], [0, y_sigma1_low],
                    color='red', linestyle='--', lw=1, label='1σ')
            ax.plot([sigma1_high, sigma1_high], [0, y_sigma1_high],
                    color='red', linestyle='--', lw=1, label='1σ')
            ax.plot([sigma2_low, sigma2_low], [0, y_sigma2_low],
                    color='red', linestyle='-.', lw=1, label='2σ')
            ax.plot([sigma2_high, sigma2_high], [0, y_sigma2_high],
                    color='red', linestyle='-.', lw=1, label='2σ')
            
        idx += 1

     #体裁調整     ax.set(xlabel='平均身長', ylabel='確率(比率)', 
           xticks=np.arange(140, 200+5, 5))
    plt.title(f'正規分布(2018年の30-39歳の男性の平均身長)', y=-0.15)
    plt.grid()
    plt.legend()
    if verbose:
        if os.path.exists('image') == False:
            os.mkdir('image')  #imageフォルダがなければ作成         plt.savefig(f'image/確率密度_{label}.png')
        
    plt.show()

 #確認用 mean, std = 172.1, 6.1 #平均値μ, 標準偏差σ
print(f'平均値の確率：標準正規分布:{norm(0, 0, 1):.2f}, 本試験:{norm(mean, mean, std):.3f}')  #標準正規分布の確認 x_population = np.arange(140, 200, 1)  #x軸の値を計算 y_population= stats.norm.pdf(x_population, loc=172.1, scale=6.1)  #確率密度関数 showGraph([x_population], [y_population], ['標準正規分布'])   

[OUT]
平均値の確率：標準正規分布:0.40, 本試験:0.065

[IN]
 #母集団の作成 mean, std = 172.1, 6.1 #平均値μ, 標準偏差σ
nums_p = int(1e6)  #乱数値用に大きなサンプルサイズnを指定 datas_p = np.random.normal(mean, std, size=nums_p)
freqs1, _x1 = np.histogram(datas_p, bins=30)
freqs2, _x2 = np.histogram(datas_p, bins=100)
freqs3, _x3 = np.histogram(datas_p, bins=1000)

#相対度数の計算（確率密度ではない）
freq1_rel = freqs1 / nums_p
freq2_rel = freqs2 / nums_p 
freq3_rel = freqs3 / nums_p 

x1, x2, x3 = (_x1[1:] + _x1[:-1]) / 2, (_x2[1:] + _x2[:-1]) / 2, (_x3[1:] + _x3[:-1]) / 2  #x軸の値を計算 print(f'形状確認：{freqs1.shape}, {freqs2.shape}, {freqs3.shape}, {freq1_rel.shape}, {freq2_rel.shape}, {freq3_rel.shape}')
print(f'合計値の確認 : {freq1_rel.sum():.3f}, {freq2_rel.sum():.3f}, {freq3_rel.sum():.3f}')
xs, ys, labels = [x_population, x1, x2, x3], [y_population, freq1_rel, freq2_rel, freq3_rel], ['標準正規分布', 'bins30', 'bins100', 'bins1000']
showGraph(xs, ys, labels, verbose=True)

[OUT]

[IN]
 #母集団の作成 mean, std = 172.1, 6.1 #平均値μ, 標準偏差σ
nums_p = int(1e6)  #乱数値用に大きなサンプルサイズnを指定  #標本1 datas_p = np.random.normal(mean, std, size=nums_p)
p_dense1, _x1 = np.histogram(datas_p, bins=30, density=True)
 #標本2 freqs2, _x2 = np.histogram(datas_p, bins=100)
bin_width = _x2[1] - _x2[0]  #階級幅 p_dense2 = freqs2 / (nums_p * bin_width)  #確率 
x1, x2 = (_x1[1:] + _x1[:-1]) / 2, (_x2[1:] + _x2[:-1]) / 2#x軸の値を計算
print(f'形状確認：{p_dense1.shape}, {p_dense2.shape}, {freq2_rel.shape},')
print(f'合計値の確認 : {freq2_rel.sum():.3f}, {p_dense2.sum():.3f}')
xs, ys, labels = [x_population, x1, x2], [y_population, p_dense1, p_dense2], ['標準正規分布', 'bins30', 'bins100']
showGraph(xs, ys, labels, verbose=True)

[OUT]