Python ラテン超方格法などで水準表を作成する「lhsmdu」

Python 回帰 グラフ化

'21/12/17更新：余剰なコードを削除して可読性を良くしました。
　はじめに、実験計画法（Design of Experiments,DOE)は広義には、直交表計画、中心複合計画、要因計画、ラテン超方格法、最適ラテン超方格法、パラメータスタディ(1因子実験)などがあります。実験やシミュレーションをする目的は、回帰式(回帰モデル)を作成することです。これにより、入力因子(設計変数)が出力指標(目的関数)に与える感度をわかり、設計因子としてどれに着目して制御すべきかの判断材料になります。

　本記事では、次の２つの実験計画法と２つの確率分布の計４つの手法で、主にシミュレーションや機械学習で使用するための入力因子(設計変数)を作成する雛形コードを載せました。作成した乱数データはcsvファイルで保存します。これで、フルDOEではなく実験数を抑制することや、ばらつきを考慮した解析をすることができます。
1.モンテカルロ法　…　lhsmdu.createRandomStandardUniformMatrix(dimension, sampling_num)
2.ラテン超方格法　…　lhsmdu.sample(dimension, sampling_num)
3.正規分布　…　np.random.normal(0.5, 0.16, sampling_num)
4.ベータ分布　…　np.random.beta(5, 2, sampling_num)

上記1, 2をするには、下記のように「lhsmdu」をpipでインストールします。

pip install lhsmdu

　本記事のコードの特徴は、下記①,②です。その下に設定例を示します。
① 因子数とその上下限を自由に設定できます。「入力因子名」とその指定範囲を（下限、上限）で辞書型で与えます。正規化したデータがよければ、上下限に(0, 1)を指定します。
② サンプリング数を指定できます。

if __name__ == "__main__":
    # parameter
    # 因子名とその乱数を生成する上下限を辞書で設定
    my_factors = {
        "height":(50, 200),
        "width":(0.06, 0.1),
        "density":(1e15, 9e15),
        "temp":(-50, 250)
    }

乱数手法別に生成した全ての因子についての行列散布図が次の4つ画像です。
下図は、モンテカルロによる結果です。そこそこ均等な分布で生成されてる様子です。

f:id:HK29:20191130184719p:plain

下図は、ラテン超方格法による結果です。ほぼ満遍なく均等に分布していることがわかります。

f:id:HK29:20191130184731p:plain

下図は、正規分布による結果です。中心にデータ点数が多く、山形状に分布してる様子がわかります。

f:id:HK29:20191130184744p:plain

下図は、ベータ分布による結果です。分布が右寄りになってる様子がわかります。本例では、np.random.beta(5, 2, sampling_num)の引数が5と2となっているためです。逆にすれば左寄りの分布になり、等しい値の場合は中心を軸に対称になります。

f:id:HK29:20191130184759p:plain

そして、下図のようにcsvファイルに手法別に出力します。

f:id:HK29:20191130191502p:plain

▼本プログラム

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
import lhsmdu
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

# 区間0～1の乱数データを指定した最小値～最大値の区間へ変換する関数
def my_calc(my_factors_column_names, my_random_np_array):
    for i, key in enumerate(my_factors_column_names):
        my_min, my_max = my_factors[key]
        my_random_np_array[i] = ((my_max - my_min) * my_random_np_array[i] + my_min)
    fixed_data_np = np.array(my_random_np_array)

    return fixed_data_np


# モンテカルロ法
def monte_carlo_func(my_factors, my_factors_column_names):
    monte_carlo = lhsmdu.createRandomStandardUniformMatrix(len(my_factors), sampling_num)
    print("monte_carlo_sampling")
    print(monte_carlo)    
    print(sampling_num)
    
    fixed_data_np = my_calc(my_factors_column_names, monte_carlo)
    monte_carlo_DF = pd.DataFrame(fixed_data_np.T)
    monte_carlo_DF.columns = my_factors_column_names
    monte_carlo_DF.to_csv("01_monte_carlo.csv", sep=",", index=False, encoding="utf-8")
    plot_matrix_scatter("01_monte_carlo", monte_carlo_DF, "winter")
    
    return monte_carlo_DF


# ラテン超方格法
def latin_hypercube_func(my_factors, my_factors_column_names):
    latin_hypercube = lhsmdu.sample(len(my_factors), sampling_num)
    print("latin_hypercube_sampling")
    print(latin_hypercube)
    print(sampling_num)
    
    fixed_data_np = my_calc(my_factors_column_names, latin_hypercube)
    latin_hypercube_DF = pd.DataFrame(fixed_data_np.T)
    latin_hypercube_DF.columns = my_factors_column_names
    latin_hypercube_DF.to_csv("02_latin_hypercube.csv", sep=",", index=False, encoding="utf-8")
    plot_matrix_scatter("02_latin_hypercube", latin_hypercube_DF, "autumn")
    
    return latin_hypercube_DF


# 正規乱数
def normal_random_func(my_factors, my_factors_column_names):
    buf=[]
    for cnt in range(len(my_factors)):
        buf.append(np.random.normal(0.5, 0.16, sampling_num))
    normal_random = np.array(buf)
    print("normal_random_sampling")
    print(normal_random)
    print(sampling_num)
    fixed_data_np = my_calc(my_factors_column_names, normal_random)
    normal_random_DF = pd.DataFrame(fixed_data_np.T, columns = my_factors_column_names)
    normal_random_DF.to_csv("03_normal_random.csv", sep=",", index=False, encoding="utf-8")
    plot_matrix_scatter("03_normal_random", normal_random_DF, "gray")

    return normal_random_DF


# ベータ乱数
def beta_random_func(my_factors, my_factors_column_names):
    buf=[]
    for cnt in range(len(my_factors)):
        buf.append(np.random.beta(5, 2, sampling_num))
    beta_random = np.array(buf)
    print("beta_random_sampling")
    print(beta_random)
    print(sampling_num)
    fixed_data_np = my_calc(my_factors_column_names, beta_random)
    beta_random_DF = pd.DataFrame(fixed_data_np.T, columns = my_factors_column_names)
    beta_random_DF.to_csv("04_beta_random.csv", sep=",", index=False, encoding="utf-8")
    plot_matrix_scatter("04_beta_random", beta_random_DF, "summer")

    return beta_random_DF


# 行列散布図を作成する関数
def plot_matrix_scatter(label, DF, my_color):
    sns.set(style="ticks", font_scale=1.2, palette=my_color, color_codes=True)
    g = sns.pairplot(DF, diag_kind="hist")
    g.fig.suptitle(label) #, fontsize=12)
    g.fig.subplots_adjust(top=0.9)
    plt.savefig(label + '.png')
    plt.close()


# 乱数を生成するメイン関数
def main():
    ##### Monte Carlo Sampling (MCS)
    df1 = monte_carlo_func(my_factors, my_factors_column_names)
    
    ##### Latin Hypercube Sampling (LHS)
    df2 = latin_hypercube_func(my_factors, my_factors_column_names)

    ##### Normal random Sampling
    df3 = normal_random_func(my_factors, my_factors_column_names)

    ##### Beta random Sampling
    df4 = beta_random_func(my_factors, my_factors_column_names)


if __name__ == "__main__":
    # parameter
    # 因子名とその乱数を生成する上下限を辞書で設定
    my_factors = {
        "height":(50, 200),
        "width":(0.06, 0.1),
        "density":(1e15, 9e15),
        "temp":(-50, 250)
    }

    # (因子をリストで作成)
    my_factors_column_names = []
    for k in my_factors.keys():
        my_factors_column_names.append(k)
    print(my_factors_column_names)

    # 生成数を設定
    sampling_num = 100
    
    # call function
    main()

●関連資料
組み合わせで水準表を作成したい場合は下記リンクを参照

hk29.hatenablog.jp

下記リンクには、plotやsnsについての色名について詳細な記載があります。

matplotlib.org

以上

<広告>