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

Python 回帰 グラフ化

'23/11/04更新：ひな形コードの可読性を高めるため、クラスと継承で書き直しました。
　はじめに、実験計画法（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

　本プログラムの入力パラメータは、次の①,②です。
① 入力因子数は任意であり、さらにそれら乱数を生成する上下限の値は自由に設定できます。キーは入力因子名で、バリューはタプルで（下限値, 上限値）のように設定します。
② サンプリング数を指定できます。

下記は、その入力パラメータの例です。入力因子はheight, width, density, tempの4つです。そして、例えばheightの場合のデータ生成範囲は50～200の間です。生成する乱数の数は200です。

if __name__ == "__main__":
    # 作成するランダムデータを辞書型で作成する
    # キーは因子名、バリューはタプルで生成データの下限と上限を指定する
    input_params = {
        "height": (50, 200),
        "width": (0.06, 0.1),
        "density": (1e15, 9e15),
        "temp": (-50, 250)
    }
    
    # 生成するサンプル数
    sampling_num = 200

乱数手法別に生成した全ての因子についての行列散布図が次の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
import datetime
now = datetime.datetime.now()
now = now.strftime("%y%m%d")

class RandomSampler:
    # インスタンスの生成
    def __init__(self, factors, sampling_num, save_file_name, plot_color):
        self.factors = factors # 入力データ
        self.factors_column_names = list(factors.keys()) # 因子名をリスト化
        self.sampling_num = sampling_num # 生成するサンプル数
        self.save_file_name = save_file_name # 保存するファイル名
        self.plot_color = plot_color # グラフにプロットする点の色

    # データ生成
    def generate_random_data(self, random_data):
        # データの逆正規化
        fixed_data_np = self._normalize_data(random_data)
        
        # pandasデータフレームにする
        random_df = pd.DataFrame(fixed_data_np.T, columns=self.factors_column_names)
        
        # csvファイルへ出力する
        random_df.to_csv(f"{self.save_file_name}.csv", sep=",", index=False, encoding="utf-8")
        
        # グラフへプロット
        self.plot_matrix_scatter(random_df)
        
        return random_df

    # データの逆正規化。生成されたデータは区間0～1の乱数データのため、指定した最小値～最大値の区間へ変換する
    def _normalize_data(self, data):
        for i, key in enumerate(self.factors_column_names):
            my_min, my_max = self.factors[key]
            data[i] = ((my_max - my_min) * data[i] + my_min)
        return np.array(data)

    # 行列散布図
    def plot_matrix_scatter(self, df):
        sns.set(style="ticks", font_scale=1.2, palette=self.plot_color, color_codes=True)
        g = sns.pairplot(df, diag_kind="hist")
        g.fig.suptitle(self.save_file_name)
        g.fig.subplots_adjust(top=0.9)
        plt.savefig(f"{self.save_file_name}.png")
        plt.close()

# モンテカロ法によるデータ生成
class MonteCarlo(RandomSampler):
    def __init__(self, factors, sampling_num, save_file_name, plot_color):
        super().__init__(factors, sampling_num, save_file_name, plot_color)

    def generate_samples(self):
        monte_carlo = lhsmdu.createRandomStandardUniformMatrix(len(self.factors), self.sampling_num)
        return self.generate_random_data(monte_carlo)

# ラテン超方格法によるデータ生成
class LatinHypercube(RandomSampler):
    def __init__(self, factors, sampling_num, save_file_name, plot_color):
        super().__init__(factors, sampling_num, save_file_name, plot_color)

    def generate_samples(self):
        latin_hypercube = lhsmdu.sample(len(self.factors), self.sampling_num)
        return self.generate_random_data(latin_hypercube)

# 正規乱数によるデータ生成
class NormalRandom(RandomSampler):
    def __init__(self, factors, sampling_num, save_file_name, plot_color):
        super().__init__(factors, sampling_num, save_file_name, plot_color)

    def generate_samples(self):
        buf=[]
        for cnt in range(len(self.factors)):
            buf.append(np.random.normal(0.5, 0.16, self.sampling_num))
        normal_random = np.array(buf)
        return self.generate_random_data(normal_random)

# ベータ関数によるデータ生成
class BetaRandom(RandomSampler):
    def __init__(self, factors, sampling_num, save_file_name, plot_color):
        super().__init__(factors, sampling_num, save_file_name, plot_color)

    def generate_samples(self):
        buf=[]
        for cnt in range(len(self.factors)):
            buf.append(np.random.beta(5, 2, self.sampling_num))
        beta_random = np.array(buf)
        return self.generate_random_data(beta_random)

def main():
    # モンテカルロ法
    monte_carlo = MonteCarlo(input_params, sampling_num, f'{now}_monte_carlo', 'winter')
    df1 = monte_carlo.generate_samples()
    print(df1)

    # ラテン超方格法
    latin_hypercube  = LatinHypercube(input_params, sampling_num, f'{now}_latin_hypercube', 'autumn')
    df2 = latin_hypercube.generate_samples()
    print(df2)

    # 正規乱数
    normal_random = NormalRandom(input_params, sampling_num, f'{now}_normal_random', 'gray')
    df3 = normal_random.generate_samples()
    print(df3)

    # ベータ乱数
    beta_random = BetaRandom(input_params, sampling_num, f'{now}_beta_random', 'summer')
    df4 = beta_random.generate_samples()
    print(df4)
    
if __name__ == "__main__":
    # 作成するランダムデータを辞書型で作成する
    # キーは因子名、バリューはタプルで生成データの下限と上限を指定する
    input_params = {
        "height": (50, 200),
        "width": (0.06, 0.1),
        "density": (1e15, 9e15),
        "temp": (-50, 250)
    }
    
    # 生成するサンプル数
    sampling_num = 200

    main()

（参考）格子状に水準表を作成したい場合は下記リンクを参照ください。

hk29.hatenablog.jp

以上

<広告>