ホーム / python / その他 / プログラミング / 統計 /

統計についての勉強~プログラミングで確認~

2021/09/12

python その他 プログラミング 統計

前回母集団とサンプル、平均、分散、標準偏差について勉強しました。統計についての勉強~平均、分散、標準偏差についての説明~
今回は母集団からサンプリングし、サンプルから母集団をどのくらいの精度で推定できるか確認していきます。
無いように間違いがある場合は指摘していただけると幸いです。

母集団を生成

まずは今回使用する母集団を生成します。今回はネットで見つけたどこかの大学？の健康診断の平均身長、標準偏差を元に正規分布データを作成します。http://web.tuat.ac.jp/~health/nenpou_data_file/kensin2012.pdf


import numpy as np
import matplotlib.pyplot as plt
import collections
import pandas as pd
#Grobal変数宣言
#----------------------------------------------------------------------------------
N = 10**6 #募集団のdata数
file_name = "height_data.csv"
#----------------------------------------------------------------------------------

#母集団を生成する関数
#----------------------------------------------------------------------------------
def make_population():
    data = np.random.normal(loc = 171.4,scale = 5.69, size = N)#平均身長の引用元 http://web.tuat.ac.jp/~health/nenpou_data_file/kensin2012.pdf
    df_data = pd.DataFrame(data)#PandasでCSVファイル出力するためにデータ型を変換
    df_data = df_data.astype("int64")#データをint型に変換して出力(この後描画しやすくするため)
    df_data.to_csv(file_name,index = False)#データをCSV形式で出力
#--------------------------------------------

#データを読み込んで前処理を行う関数
#----------------------------------------------------------------------------------
def preprocess(read_file_name):
    read_data = pd.read_csv(read_file_name)#データの読み込み
    read_data = read_data.T#転置をとってこの後の処理ができるようにする
    read_data = read_data.values.tolist()[0]#配列に変換()
    count_data = dict(sorted(collections.Counter(read_data).items()))#データの頻度を求める
    data_int_std = np.std(read_data)#int型に変換した後の標準編差
    data_int_mean = np.mean(read_data)#int型に変換した後の平均
    print(data_int_std)
    print(data_int_mean)

    x = list(count_data.keys())#身長データを抜粋
    y = list(count_data.values()) #頻度データを抜粋
    return x,y
#----------------------------------------------------------------------------------

#描画する関数
#----------------------------------------------------------------------------------
def plot_data(x,y):
    plt.bar(x,y)
    plt.xlabel("身長(cm)",fontname="MS Gothic")#デフォルトだと文字化けしてしまう為、フォント指定（Windowsの向け）
    plt.ylabel("人数(人)",fontname="MS Gothic")
    plt.show()    
#----------------------------------------------------------------------------------

if __name__ == "__main__":
    #make_population()#データ生成する場合はこちらのコメントアウトを解除
    x,y = preprocess(file_name)
    plot_data(x,y)

今回使用する身長データは下記の通りです。

今回使用する母集団
この母集団の平均は170.89、標準偏差が5.69の正規分布となっております。
また、今回作成したこのソースコードを編集していきます。
最後に最終的なソースコードを記載します。
(行き当たりばったりで書きなぐったコードになってしまいました。。。クラスを使った方がすっきりしたのではないかと反省です。。。)

母集団からサンプリングする

母集団が作成できたので次はサンプルを作成していきます。
サンプルサイズをいくつにするか考えていきます。

前回の記事で分散を求める際にデータが少ない場合はn-1で割るという話をしました。
nで割った場合とn-1で割った場合の誤差を求めて、それを元にn数がいくつが良いかを考えていきます。(こちらの考えは合ってるか不明です。。。)

今回はサンプルの標準偏差を母集団に近づけたいため、標準偏差で考えていきます。
nで割って求めた標準偏差とn-1で割って求めた標準偏差の誤差は次のように表します。
\begin{equation}\frac{\sqrt{V_n}}{\sqrt{V_{n-1}}} = \sqrt{1 - \frac{1}{n}} \end{equation}

これをプロットしたものが下記の図です。

上記図よりサンプルサイズが10あれば95 %程度一致し、65あれば99.2 %程度で一致する。。。はずです。

この考察を元にサンプルサイズを10とした場合と65とした場合で母集団とどれほど一致するかを見ていきます。

サンプル数と精度の関係

サンプルサイズによる母集団の値からのバラつき

先ほどの結果より、サンプルサイズは最低65個必要だという推定をしました。
ここではサンプルサイズ10の場合と65の場合で母集団の平均/標準偏差からどれほどばらつきに違いが出るか見ていきます。

バラつきの評価として母集団の平均値/標準偏差とサンプルの平均値/標準偏差の標準偏差を求めてバラつきを見ていきます。
\begin{equation}今回求める標準偏差 = \sqrt{\frac{\sum{(x_i-\mu)}}{n}} \end{equation}

また、母集団と最大何%ずれているかを相対誤差を使って評価します。
(相対誤差についてはこちらを参照して下さい。相対誤差の計算方法と意義)

まずはサンプルサイズ10の場合です。今回はサンプル数を10として評価してます。
サンプルサイズが10の場合の平均のバラつき

上図の黒い線は母集団の平均の値です。
こちらを見るとサンプルによっては平均値174付近となってしまい、相対誤差で見ると1.7 %程度ずれております(今回求めた相対誤差は符号を無視しております。)。また、こちらのデータの母集団の値との標準偏差は1.5となっております。

次は母集団の標準偏差とサンプルの標準偏差がどのくらい違いがあるか見ていきます。
(意義のあるデータかは不明ですが。。。)
サンプルサイズが10の場合の標準偏差のバラつき

こちらを見ると最大の相対誤差は40.5 %程度、標準偏差としては1.30程度となっております。


次にサンプルサイズを65とした場合のデータを見ていきます。(サンプル数は先ほどと同じ10です。)


サンプルサイズが65の場合の平均のバラつき


サンプル数を増やしてことにより、母集団の値に近づくようになりました。
相対誤差も1.4 %程度となり、サンプルサイズが10の場合の1.7 %と比較しても精度がアップしていることがわかります。

また、標準偏差も0.87となっておりサンプルサイズが10の時の1.30と比較するとバラつきも少なくなっております。


次はサンプルサイズを65にした場合の標準偏差のバラつきです10
サンプルサイズが65の場合の標準偏差のバラつき

相対誤差の最大値は14.2 %でした。サンプルサイズが10の時の40 %と比較すると大分改善されています。
また、標準偏差も0.38だったためサンプルサイズ10の時の1.30と比較するとバラつきも少なくなってます。

まとめ

母集団からサンプリングをしてサンプル数の数によってどの程度母集団からずれていくかを見ていきました。サンプル数を増やすことで母集団の値に近づいていくことがわかりました。

個人的にはサンプルサイズが10程度でも十分評価できるのではないかと思いました。(評価対象/評価項目によると思いますが。。。)

(補足)ソースコード

思いつくままにソースを書いたため、あまりきれいなコードではないですが、参考までに今回のコードを添付します。Global宣言に記載した変数を色々変えて挙動を確認してみてください。
import numpy as np
import matplotlib.pyplot as plt
import collections
import pandas as pd
import random
#Grobal変数宣言
#----------------------------------------------------------------------------------
N = 10**6 #募集団のdata数
FILE_NAME = "height_data.csv"
SAMPLE_NUM_10 = 10#サンプル数10

SAMPLE_SIZE_10 = 10#サンプルサイズ10
SAMPLE_SIZE_65 = 65#サンプルサイズ65

POPURATION_MEAN = 170.89#母集団の平均(一番初めに作成した母集団から算出)
POPURATION_STD = 5.69#母集団の標準偏差(一番初めに作成した母集団から算出)
#----------------------------------------------------------------------------------

#母集団を生成する関数
#----------------------------------------------------------------------------------
def make_population():
    data = np.random.normal(loc = 171.4,scale = 5.69, size = N)#平均身長の引用元 http://web.tuat.ac.jp/~health/nenpou_data_file/kensin2012.pdf
    df_data = pd.DataFrame(data)#PandasでCSVファイル出力するためにデータ型を変換
    df_data = df_data.astype("int64")#データをint型に変換して出力(この後の処理をしやすくするため)
    df_data.to_csv(FILE_NAME,index = False)#データをCSV形式で出力
#--------------------------------------------

#読み込んだデータを返す関数(引数：ファイル名　戻り値：読み込んだデータ)
#----------------------------------------------------------------------------------
def read_data(read_file_name):
    read_data = pd.read_csv(read_file_name)#データの読み込み
    return read_data

#---------------------------------------------------------------------------------

#データの頻度を求める関数(引数：頻度を求めるデータ　戻り値：読み込んだデータを配列にしたもの、身長データ、頻度データ)
#----------------------------------------------------------------------------------
def preprocess(read_data):
    read_data = read_data.T#転置をとってこの後の処理ができるようにする
    read_data = read_data.values.tolist()[0]#配列に変換()
    count_data = dict(sorted(collections.Counter(read_data).items()))#データの頻度を求める
    x = list(count_data.keys())#身長データを抜粋
    y = list(count_data.values()) #頻度データを抜粋
    return read_data,x,y
#----------------------------------------------------------------------------------

#平均などの統計データを求める関数(引数：平均/標準編差を求めるデータ　返り値：データの標準編差、データの平均)
#----------------------------------------------------------------------------------
def cal_statics(data):
    data_std = np.std(data)
    data_mean = np.mean(data)
    return data_std,data_mean
#----------------------------------------------------------------------------------

#描画する関数(引数：xデータ、yデータ、x軸のラベル。y軸のラベル、棒グラフか散布図かの選択)
#----------------------------------------------------------------------------------
def plot_data(x,y,x_rabel,y_rabel,mode):
    if mode == "bar":
        plt.bar(x,y)
    elif mode == "scatter":
        plt.scatter(x,y)
    plt.xlabel(x_rabel,fontname="MS Gothic")#デフォルトだと文字化けしてしまう為、フォント指定（Windowsの向け）
    plt.ylabel(y_rabel,fontname="MS Gothic")
    plt.show()    
#----------------------------------------------------------------------------------

#hline付きで描画する場合(引数：xデータ、yデータ、x軸のラベル。y軸のラベル、棒グラフか散布図かの選択、y軸に線を引く値)
#----------------------------------------------------------------------------------
def plot_data_hline(x,y,x_rabel,y_rabel,mode,y_line):
    if mode == "bar":
        plt.bar(x,y)
    elif mode == "scatter":
        plt.scatter(x,y)
    plt.xlabel(x_rabel,fontname="MS Gothic")#デフォルトだと文字化けしてしまう為、フォント指定（Windowsの向け）
    plt.ylabel(y_rabel,fontname="MS Gothic")
    plt.hlines(y=y_line,xmin=np.min(x),xmax=np.max(x))
    plt.show()    

#----------------------------------------------------------------------------------

#サンプル数をいくつが良いか検討する関数(引数：データ個数)
#----------------------------------------------------------------------------------
def sample_consider(sample_size):
    x = np.array(range(1,sample_size+1))
    y = np.sqrt(1-(1/x))*100
    return x,y

#----------------------------------------------------------------------------------

#サンプルを作成する関数(引数：母集団データ、サンプルのサイズ　戻り値：作成したサンプル)
#----------------------------------------------------------------------------------
def sampling(raw_data,sample_size):
    random_list = []#サンプルリングする際に重複しないランダムな数を格納する配列
    sample_list = []#サンプリングした物を格納する配列
    while len(random_list)<sample_size:#母集団から重複が無いようにサンプリングをする。
        n = random.randint(0,N)
        if not n in random_list:
            random_list.append(n)
    for i in random_list:
        sample_list.append(raw_data[i])
    return sample_list
#----------------------------------------------------------------------------------

#母集団の平均・標準編差を使ってサンプルの平均・標準編差との標準編差を求める(引数：データ、母集団の平均/標準編差、母集団との標準編差)
#----------------------------------------------------------------------------------
def cal_popuration_sample_std(data,mean):
    V = 0
    max_diff = 0#母集団との誤差最大値
    for i in range(len(data)):
        V += (data[i]-mean)**2#標準編差を求める
        max_diff = max(max_diff,abs(data[i]-mean))
    V = V/len(data)
    s = np.sqrt(V)
    max_r_diff = max_diff*100/mean
    print("相対誤差最大値"+str(max_r_diff))
    return s

#----------------------------------------------------------------------------------

#サンプル数分のサンプリングを実施し平均値、標準編差を求め描画する(引数：サンプル数、母集団のデータ、サンプルサイズ)
#----------------------------------------------------------------------------------
def cal_sample_statics(num,popuration_data,sample_size):
    sample_list_std = []
    sample_list_mean = []
    x = range(num)
    for i in range(num):#引数num回サンプリングする
        sample_list = sampling(popuration_data,sample_size)#母集団データからサンプルサイズ分サンプリング
        sample_std,sample_mean = cal_statics(sample_list)
        sample_list_std.append(sample_std)
        sample_list_mean.append(sample_mean)    
    s_std = cal_popuration_sample_std(sample_list_std,POPURATION_STD)#サンプルのデータが母集団の値からどれだけばらつくかを求める(平均値)
    s_mean = cal_popuration_sample_std(sample_list_mean,POPURATION_MEAN)#サンプルのデータが母集団の値からどれだけばらつくかを求める(標準編差)
    print("標準編差の標準編差："+str(s_std))
    print("平均の標準編差："+str(s_mean))
    plot_data_hline(x,sample_list_std,"サンプル数","標準編差","scatter",POPURATION_STD)
    plot_data_hline(x,sample_list_mean,"サンプルサイズ","平均","scatter",POPURATION_MEAN)
#----------------------------------------------------------------------------------



#main関数
#----------------------------------------------------------------------------------
if __name__ == "__main__":
    # make_population()#データ生成する場合はこちらのコメントアウトを解除
    execel_data = read_data(FILE_NAME)
    int_data,data_height,data_freq, = preprocess(execel_data)
    # plot_data(data_height,data_freq,"身長(cm)","人数(人)","bar")#作成したデータをプロット(基本コメントアウト)
    # data_sample_x , data_sample_y, = sample_consider(100)#サンプル数がいくつが良いか検討用。基本コメントアウト
    # plot_data(data_sample_x,data_sample_y , "サンプル数","精度(%)","scatter")#検討したサンプルを描画する用。基本コメントアウト
    sample_list_10 = sampling(int_data,SAMPLE_NUM_10)#サンプルサイズ10のサンプルを生成
    cal_sample_statics(SAMPLE_NUM_10,int_data,SAMPLE_SIZE_10)
    cal_sample_statics(SAMPLE_NUM_10,int_data,SAMPLE_SIZE_65)
#----------------------------------------------------------------------------------





自己紹介

はじめまして 社会人になってからバイクやプログラミングなどを始めました。 プログラミングや整備の記事を書いていますが、独学なので間違った情報が多いかもしれません。 間違っている情報や改善点がありましたらコメントしていただけると幸いです。

X(旧Twitter)

フォローお願いします！

QooQ