世の中、「平均」だけで統計を語ってる人や記事が多すぎるので、いろんな分布を自在に作ってみようと思いました。せっかくなので分散も固定して、平均50・分散10の分布を色々作ってみましたよ！

Optuna インストール

Optuna はほとんど何も考えずに最適化できるので楽で良いですね。I love Optuna!

!pip install optuna

必要なライブラリのインポート

import optuna
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt

正規分布

平均50・分散10の正規分布を作るためのObjectiveを作ってみました。え？平均50・分散10の正規分布なんて、Optuna使わなくても作れるって？その通りなんですが、他の分布を使うための練習ですよ、練習。

class Objective:
    def __init__(self, target_mean = 50, target_var = 10):
        self.best_score = 530000
        self.best_X = None
        self.target_mean = target_mean
        self.target_var = target_var

    def __call__(self, trial):
        a  = trial.suggest_uniform("a", 0, 100)
        b = trial.suggest_uniform("b", 0, 100)
        rand = a + b * np.random.randn(10000)
        score = abs(self.target_mean - rand.mean()) + abs(self.target_var - rand.var())
        if self.best_score > score:
            self.best_score = score
            self.best_X = rand
        return score

objective = Objective()
study = optuna.create_study()
study.optimize(objective, n_trials=100)

rand = objective.best_X
plt.title("mean = {0:.3f}, var={1:.3f}".format(rand.mean(), rand.var()))
sns.distplot(rand, bins=50)
plt.grid()
plt.show()

はい、厳密に一致するわけではもちろんありませんが、平均50・分散10の正規分布が作成できましたね。ちなみに、大学入試の模試などでよく使われる「偏差値」ってやつは、平均50・分散10の正規分布を仮定しているわけです。

一様分布

次は、平均50・分散10の一様分布を作成します。

class Objective:
    def __init__(self, target_mean = 50, target_var = 10):
        self.best_score = 530000
        self.best_X = None
        self.target_mean = target_mean
        self.target_var = target_var

    def __call__(self, trial):
        a  = trial.suggest_uniform("a", 0, 100)
        b = trial.suggest_uniform("b", 0, 100)
        rand = a + b * np.random.rand(10000)
        score = abs(self.target_mean - rand.mean()) + abs(self.target_var - rand.var())
        if self.best_score > score:
            self.best_score = score
            self.best_X = rand
        return score

objective = Objective()
study = optuna.create_study()
study.optimize(objective, n_trials=100)

rand = objective.best_X
plt.title("mean = {0:.3f}, var={1:.3f}".format(rand.mean(), rand.var()))
sns.distplot(rand, bins=50)
plt.grid()
plt.show()

普通、こんな分布を考える人はいないんでしょうけれども、これだって平均50・分散10の分布なわけです。

指数分布

次は、平均50・分散10の指数分布を作成します。

class Objective:
    def __init__(self, target_mean = 50, target_var = 10):
        self.best_score = 530000
        self.best_X = None
        self.target_mean = target_mean
        self.target_var = target_var

    def __call__(self, trial):
        a  = trial.suggest_uniform("a", 0, 100)
        b = trial.suggest_uniform("b", 0, 1000)
        c = trial.suggest_uniform("c", 0, 1000)
        rand = a + b * np.random.exponential(c, 10000)
        score = abs(self.target_mean - rand.mean()) + abs(self.target_var - rand.var())
        if self.best_score > score:
            self.best_score = score
            self.best_X = rand
        return score

objective = Objective()
study = optuna.create_study()
study.optimize(objective, n_trials=1000)

rand = objective.best_X
plt.title("mean = {0:.3f}, var={1:.3f}".format(rand.mean(), rand.var()))
sns.distplot(rand, bins=50)
plt.grid()
plt.show()

指数分布をフィッティングさせるのはなかなか大変ですね...

平均だけで語ってはいけないよ、と言いたい時によく使われる印象があるのが指数分布です。これは、平均値と最頻値が異なること、最大値が極端に高いことなどが特徴ですね。

β分布

次は、平均50・分散10のベータ分布を作成します。

class Objective:
    def __init__(self, target_mean = 50, target_var = 10):
        self.best_score = 530000
        self.best_X = None
        self.target_mean = target_mean
        self.target_var = target_var

    def __call__(self, trial):
        a  = trial.suggest_uniform("a", 0, 100)
        b = trial.suggest_uniform("b", 0, 100)
        rand = a + b * np.random.beta(2, 5, 10000)
        score = abs(self.target_mean - rand.mean()) + abs(self.target_var - rand.var())
        if self.best_score > score:
            self.best_score = score
            self.best_X = rand
        return score

objective = Objective()
study = optuna.create_study()
study.optimize(objective, n_trials=100)

rand = objective.best_X
plt.title("mean = {0:.3f}, var={1:.3f}".format(rand.mean(), rand.var()))
sns.distplot(rand, bins=50)
plt.grid()
plt.show()

平均値と最頻値がほとんど同じなので、これもまた分布をよく見ないと分かりにくい例ですよね。

分布は他にも二項分布とかガンマ分布とか色々あるんですが、そろそろ飽きたのでこのへんで終わりとします。

２つの正規分布

平均だけで統計を語るのはよくないよ、の別の例ですね。

class Objective:
    def __init__(self, target_mean = 50, target_var = 10):
        self.best_score = 530000
        self.best_X = None
        self.target_mean = target_mean
        self.target_var = target_var

    def __call__(self, trial):
        a  = trial.suggest_uniform("a", 0, 48)
        b = trial.suggest_uniform("b", 0, 20)
        c = trial.suggest_uniform("c", 52, 100)
        d = trial.suggest_uniform("d", 0, 20)
        rand1 = a + b * np.random.randn(10000)
        rand2 = c + d * np.random.randn(10000)
        rand = np.concatenate([rand1, rand2])
        score = abs(self.target_mean - rand.mean()) + abs(self.target_var - rand.var())
        if self.best_score > score:
            self.best_score = score
            self.best_X = rand
        return score

objective = Objective()
study = optuna.create_study()
study.optimize(objective, n_trials=1000)

rand = objective.best_X
plt.title("mean = {0:.3f}, var={1:.3f}".format(rand.mean(), rand.var()))
sns.distplot(rand, bins=50)
plt.grid()
plt.show()

平均50・分散10の分布なのに、平均値にはほとんどデータがないっていう例ですね。

２つのβ分布

同じように２つのβ分布でやってみました。

class Objective:
    def __init__(self, target_mean = 50, target_var = 10):
        self.best_score = 530000
        self.best_X = None
        self.target_mean = target_mean
        self.target_var = target_var

    def __call__(self, trial):
        a  = trial.suggest_uniform("a", 0, 50)
        b = trial.suggest_uniform("b", 0, 100)
        c = trial.suggest_uniform("c", 50, 100)
        d = trial.suggest_uniform("d", 0, 100)
        rand1 = a + b * np.random.beta(2, 5, 10000)
        rand2 = c + d * np.random.beta(5, 2, 10000)
        rand = np.concatenate([rand1, rand2])
        score = abs(self.target_mean - rand.mean()) + abs(self.target_var - rand.var())
        if self.best_score > score:
            self.best_score = score
            self.best_X = rand
        return score

objective = Objective()
study = optuna.create_study()
study.optimize(objective, n_trials=1000)

rand = objective.best_X
plt.title("mean = {0:.3f}, var={1:.3f}".format(rand.mean(), rand.var()))
sns.distplot(rand, bins=50)
plt.grid()
plt.show()