公平な不平等、不公平な平等

長さ1000cmの棒を一様乱数で1000個に切り分けたら、どんな長さの破片ができるか？

１cmの破片が1000個できるのかな、と想像しがちなところですが、実際にはこうなります。

これはパソコンでシミュレーションした結果のヒストグラムです。
確かに破片１個の“平均”は１cmなのですが、数で言えば平均以下の小さい破片の方がずっと多く、
その一方で、ごく小数の極端に長い破片があります。
最も数が多いのは、最も短い0.0〜0.18cmで、ここに160個以上の破片が含まれています。
反対に、最も長いものは 8cm以上、次いで7cm台、6cm台にも、ごく少数の破片があります。

以上は、ほんの数行のPythonスクリプトで確かめることができます。

# 長い棒を一様乱数で切ったら、破片の分布はどうなる？

import numpy as np
import matplotlib.pyplot as plt

rds = np.sort(np.random.rand(1000) * 1000)   # 1000個の乱数をソート
sample = [ rds[i+1] - rds[i] for i in range(rds.size-1) ] # 乱数の間隔を取得する
sample[-1] = rds[0] + (1000 - rds[-1])   # 最後は末尾と先頭をつなげる

plt.hist( sample, bins=50 )
plt.show()

たとえPytonを知らなくても、エクセルで1000個の乱数を作って試すこともできます。

手間を厭わなければ、実験して確かめることもできます。 >> [id:rikunora:20091213]

これがタイトルに掲げた「公平な不平等」です。
乱数は公平です。しかし、公平な乱数で分配した結果は、不平等なのです。
“公平＝平等”という思い込みは、必ずしも正しくありません。

似たようなことを、交換によって確かめてみましょう。

・1000人が当初1.0ずつの財産を持ち、お互いにランダムに財産を交換する。
・交換は、ランダムに選んだ２人がお互いの財産を出し合い、それを一様乱数によって振り分ける。
・交換を10万回繰り返す。

結果はこうなりました。
赤い線は「指数分布」と呼ばれている曲線で、理論上はこうなる、という形です。

さらに、交換のルールを少し変えてみましょう。

・1000人が当初1.0ずつの財産を持ち、お互いにランダムに財産を交換する。
・交換は、ランダムに選んだ２人が少ない方を越えない財産を出し合い、それを一様乱数によって振り分ける。
・交換を10万回繰り返す。

先ほどとの違いは「少ない方を越えない」という点で、たとえ少ない方が負けても全財産を失わないようにとの配慮からです。

結果は極端で、ごく一握りの勝ち組以外、大半はほとんど０になります。
なぜこうなるかと言うと、いったん財産が０近くになると、そこから抜け出すのが極めて困難だからです。

これでは余りにも勝ち負けがはっきりしているので、ハンディを付けましょう。

・交換は、ランダムに選んだ２人の財産の２乗の和が一定になるように乱数で振り分ける。

どういうことかと言うと、金持ちは持てる財産の大きさに比例してハンディを負え、というルールです。

先ほどより、だいぶ平等に近づきました。
※ 赤い線は「ガンマ分布」と呼ばれている曲線です。
※ なんとなく当てはまりそうですが、理論的にこれが正解なのかどうか、私にはよく分かりません。

もっともっとハンディを付けたら、どうなるか。

・交換は、ランダムに選んだ２人の財産の３乗の和が一定になるように乱数で振り分ける。

金持ちは持てる財産の２乗に比例してハンディを負えという、かなり金持ちに厳しいルールです。

さらに平等に近づいてきました。

反対に、金持ちが有利になるようなハンディを付けてみたら、どうなるか。

・交換は、ランダムに選んだ２人の財産の平方根の和が一定になるように乱数で振り分ける。

予想通り、かなり不平等な結果となりました。

以上の結果をまとめて描くと、こうなります。

このグラフは、分配ルールのハンディ乗数を0.5(平方根)〜4.0まで、0.5刻みに変えた結果を重ねて描いたものです。
（人数は8000人に増やしています。）
ハンディの大きさに応じて、結果が格差から平等に変わる様子が見て取れることと思います。

昔から言い古されてきたことなのですが、自由とは格差社会であり、かといって出る杭を打つ社会に生まれた天才は不幸です。
この事実は今も昔も変わりませんが、今が昔と違うところは、分配ルールによって平等が調整できる姿を、誰もがパソコン１つで試せるようになったことです。

不平等とは、地震や台風のようにコントロール不能な災害ではなく、人がコントロールできる問題です。
もし効率を求める組織だったなら、結果としての不平等より、チャンスとしての公平を敷くべきかもしれません。
あるいは調和を求める社会だったなら、結果としての平等を重んじ、方法としての不公平を受け入れるべきかもしれません。
ひょっとすると、全体最適化のためには適切なセグメンテーション、クラス分けや階級化が必要なのかもしれません。
いずれにせよ、目的に叶ったルールは数字の上で選択可能であり、
たとえそのモデル化が不完全だったとしても、感情にまかせた言葉をぶつけ合うよりずっと合理的だと思うのです。

* なぜ統計学では釣り鐘型の分布が使われ、物理現象では右肩下がりの分布が使われるのか
>> [id:rikunora:20170321]

# 交換のルールを変えてみたら、分布はどのように変わるのか

import numpy as np
import random
import scipy.optimize
import matplotlib.pyplot as plt

class ExchgRule:
    
    # ２つの数の合計をランダムに分配する
    def exchg(self, a, b):
        rd = np.random.rand()
        s = a + b
        p = (s * rd)
        q = (s * (1.0 - rd))
        return ( p, q )
    
    # 小さい方の数と等量（双方が出せるだけの金額）をランダムに分配する
    def exchg_min(self, a, b):
        rd = np.random.rand()
        mn = min( a, b )
        mx = max( a, b )
        p = 2 * mn * rd
        q = 2 * mn * (1-rd)
        return ( p, q + mx - mn )
    
    # ２つの数のｎ乗の和が一定になるように分配する
    def exchg_pwn(self, a, b, n):
        rd = np.random.rand()
        s = np.power(a, n) + np.power(b, n)
        p = np.power(s * rd, 1/n )
        q = np.power(s * (1.0 - rd), 1/n )
        ratio = (a+b) / (p+q)   # 合計が一定となるように標準化
        return( ratio * p, ratio * q )
    
    def run(self):
        N_SAMPLE  =  1000    # 粒子数
        N_EXCHG  = 100000    # 交換回数
        
        # フィットさせたい関数、指数分布
        def exfunc(x, a, b):
            return a * np.exp( - x * b )
            # TypeError: only size-1 arrays can be converted to Python scalars
            # mathパッケージのlogやexpを用いるとエラーが出ます。
            # numpyパッケージのlogやexpを用いればオッケーです。
        
        # フィットさせたい関数、正規分布
        def nmfunc(x, a, b, c):
            return a * np.exp( - (x-c)**2 * b )
        
        # フィットさせたい関数、ガンマ分布っぽいもの
        def gmfunc(x, a, b, c):
            return a * np.power(x, b) * np.exp( - x * c )
        
        # いろんな分布からスタートしてみる
        sample = np.random.rand(N_SAMPLE) # 一様分布
        # sample = np.random.randn(N_SAMPLE) # 標準正規分布
        # sample = np.random.normal( 10, 2, N_SAMPLE )    # 正規分布、平均をずらした
        # sample = np.random.exponential( scale=1.0, size=N_SAMPLE )    # 指数分布
        
        sample = np.abs( sample )   # 絶対値に直す
        
        for i in range(N_EXCHG):
            a, b = random.sample( range(N_SAMPLE), 2 ) # ランダムに２つの数を選ぶ
            
            # p, q = self.exchg( sample[a], sample[b] )
            # p, q = self.exchg_min( sample[a], sample[b] )  # 双方が出せるだけを分配
            p, q = self.exchg_pwn( sample[a], sample[b], 2 )  # ｎ乗だったらどうなる
            
            sample[a] = p
            sample[b] = q
        
        ret = plt.hist( sample, bins=50 )   # ヒストグラムを描く
        
        # 曲線あてはめを試みる
        fit_func = exfunc   # exfunc # gmfunc   # 関数名が直接代入できるって便利.
        hist_x = ret[1][:-1] # ヒストグラムの結果は返り値に入っている
        hist_y = ret[0]
        param, cov = scipy.optimize.curve_fit( fit_func, hist_x, hist_y )
        print( param )
        fit_y = fit_func( hist_x, *param )
        
        plt.plot( hist_x, fit_y, '-', color="red")  # 曲線を描く
        
        plt.show()
        
        print( "ave={:.05f}, std={:.05f}, {:.05f}〜{:.05f}".format( \
            np.mean(sample), np.std(sample), np.min(sample), np.max(sample)) )
    
if __name__ == '__main__':
    me = ExchgRule()
    me.run()