'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つ画像です。
下図は、モンテカルロによる結果です。そこそこ均等な分布で生成されてる様子です。
下図は、ラテン超方格法による結果です。ほぼ満遍なく均等に分布していることがわかります。
下図は、正規分布による結果です。中心にデータ点数が多く、山形状に分布してる様子がわかります。
下図は、ベータ分布による結果です。分布が右寄りになってる様子がわかります。本例では、np.random.beta(5, 2, sampling_num)の引数が5と2となっているためです。逆にすれば左寄りの分布になり、等しい値の場合は中心を軸に対称になります。
そして、下図のようにcsvファイルに手法別に出力します。
▼本プログラム
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)
random_df = pd.DataFrame(fixed_data_np.T, columns=self.factors_column_names)
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
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
以上
<広告>
リンク
リンク
リンク
