• TOP
  • ブログ
  • 【論文紹介】確率分布『Modified Lognormal Power-Law Distributi...

ホットリンク 公式ブログ

データから見えてくる注目のトピックスをとりあげます。

テック

【論文紹介】確率分布『Modified Lognormal Power-Law Distribution』で実験してみた

20170221_shutterstock_555665161

開発本部研究開発グループR&Dチームのサンティです。

今回のテックブログはエンジニア向けの論文(下記参考1の文献、以下「参考1の論文」)のご紹介です。面白い確率分布を扱った論文なのですが、最近業務上でも類似の分布を見たばかりであり、強く共感したのでご紹介します。

この論文は天体物理学分野の論文であり、対数正規分布とベキ分布で成り立つ分布を持つ現象を扱っています。そのような2つの分布で成り立つ現象を扱う場合、先行研究では分布の分割条件を付け、それぞれの領域のモデルを作るのが一般的でした。 しかし、参考1の論文では上記のような事象を分割条件なしで一つの関数ですっきりと記述する数理モデルを提案しています。 提案された確率分布は対数正規分布とベキ分布のマリアージュなので「Modified Lognormal Power-Law Distribution」(略してMLP分布)と命名されました。MPL分布と対数正規分布の確率密度の比較は下記の図のとおりです。 興味深い点は、この論文に取り上げられた事象は宇宙分野だけではく、自然科学、生物学、ソーシャル科学、金融科学においてもよく起きていることです。 もしかしたら、仕事に応用できるかもしれないと思い、少し実験してみました。

まずは簡単にご紹介

参考1の論文の著者がWikipedia上にも論文の要約を記載していますが、内容の誤りもあるようで、あまり広く知られていないようです。 せっかくの綺麗な新種の分布を無駄にしないために、実験で分かったことをエンジニア向けに簡単に紹介したいと思います。

なお、本実験の記述言語はJulia言語バージョン0.5です。 必要なパッケージをインストールしたら、どなたでも下記の実験を再現できると思います。 (データは別途公開サイトから取得する必要があります。)

さあ、はじめましょう。 まずは必要なパッケージを取り込みます。

In [1]:

# 必要なパッケージを取り込む
using StatsBase
using DataFrames
using Gadfly

つぎに参考1の論文からMLP分布の乱数発生関数を記述します。

In [2]:

# MLP分布に従う乱数を発生させる関数
# N:発生乱数の個数
# α,μ,σ: MLPのパラメータ
function randMLP(N,α,μ,σ)
    z = exp(μ + σ*randn(N) - log(rand(N))/α)
end

Out[2]:

randMLP (generic function with 1 method)

早速、試しに初めてのMLP分布の乱数を1万個発生してみましょう。

In [3]:

x = randMLP(10000,2.0,8.0,0.5);
plot(x=x,Geom.density,
            Scale.x_log10, Scale.y_log10, 
            Guide.title("MLP分布で発生させた乱数の確率密度"))

Out[3]:

x 10-4 10-3 10-2 10-1 100 101 102 103 104 105 106 107 108 109 1010 1011 1012 1013 10-3.0 10-2.8 10-2.6 10-2.4 10-2.2 10-2.0 10-1.8 10-1.6 10-1.4 10-1.2 10-1.0 10-0.8 10-0.6 10-0.4 10-0.2 100.0 100.2 100.4 100.6 100.8 101.0 101.2 101.4 101.6 101.8 102.0 102.2 102.4 102.6 102.8 103.0 103.2 103.4 103.6 103.8 104.0 104.2 104.4 104.6 104.8 105.0 105.2 105.4 105.6 105.8 106.0 106.2 106.4 106.6 106.8 107.0 107.2 107.4 107.6 107.8 108.0 108.2 108.4 108.6 108.8 109.0 109.2 109.4 109.6 109.8 1010.0 1010.2 1010.4 1010.6 1010.8 1011.0 1011.2 1011.4 1011.6 1011.8 1012.0 10-5 100 105 1010 1015 10-3.0 10-2.5 10-2.0 10-1.5 10-1.0 10-0.5 100.0 100.5 101.0 101.5 102.0 102.5 103.0 103.5 104.0 104.5 105.0 105.5 106.0 106.5 107.0 107.5 108.0 108.5 109.0 109.5 1010.0 1010.5 1011.0 1011.5 1012.0 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 -1.50 -1.45 -1.40 -1.35 -1.30 -1.25 -1.20 -1.15 -1.10 -1.05 -1.00 -0.95 -0.90 -0.85 -0.80 -0.75 -0.70 -0.65 -0.60 -0.55 -0.50 -0.45 -0.40 -0.35 -0.30 -0.25 -0.20 -0.15 -0.10 -0.05 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00 1.05 1.10 1.15 1.20 1.25 1.30 1.35 1.40 1.45 1.50 1.55 1.60 1.65 1.70 1.75 1.80 1.85 1.90 1.95 2.00 2.05 2.10 2.15 2.20 2.25 2.30 2.35 2.40 2.45 2.50 2.55 2.60 2.65 2.70 2.75 2.80 2.85 2.90 2.95 3.00 -2 0 2 4 -1.5 -1.4 -1.3 -1.2 -1.1 -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.0 MLP分布で発生させた乱数の確率密度 Created with Snap 

対数スケールで確率密度をプロットしましたので見た目は当然ながらピークの左側は対数正規分布に似ていて、右側はちょっと歪んでテールが長いですね。次に相補累積分布関数をみてみましょう。

In [4]:

myecdf = ecdf(x)
function myccdf(x) 1 - myecdf(x) end
xs = 1.0e2:1.0e2:1.0e4
randdf = DataFrame(x=xs, y=map(i -> myccdf(i), xs ));
plot(randdf,x=:x,y=:y,
        Scale.x_log10,Scale.y_log10,
        Guide.title("MLP分布で発生させた乱数の相補累積分布"))

Out [4]:

x 10-0.5 100.0 100.5 101.0 101.5 102.0 102.5 103.0 103.5 104.0 104.5 105.0 105.5 106.0 106.5 100.0 100.1 100.2 100.3 100.4 100.5 100.6 100.7 100.8 100.9 101.0 101.1 101.2 101.3 101.4 101.5 101.6 101.7 101.8 101.9 102.0 102.1 102.2 102.3 102.4 102.5 102.6 102.7 102.8 102.9 103.0 103.1 103.2 103.3 103.4 103.5 103.6 103.7 103.8 103.9 104.0 104.1 104.2 104.3 104.4 104.5 104.6 104.7 104.8 104.9 105.0 105.1 105.2 105.3 105.4 105.5 105.6 105.7 105.8 105.9 106.0 100 102 104 106 100.0 100.2 100.4 100.6 100.8 101.0 101.2 101.4 101.6 101.8 102.0 102.2 102.4 102.6 102.8 103.0 103.2 103.4 103.6 103.8 104.0 104.2 104.4 104.6 104.8 105.0 105.2 105.4 105.6 105.8 106.0 10-2.25 10-2.00 10-1.75 10-1.50 10-1.25 10-1.00 10-0.75 10-0.50 10-0.25 100.00 100.25 100.50 100.75 101.00 101.25 10-2.00 10-1.95 10-1.90 10-1.85 10-1.80 10-1.75 10-1.70 10-1.65 10-1.60 10-1.55 10-1.50 10-1.45 10-1.40 10-1.35 10-1.30 10-1.25 10-1.20 10-1.15 10-1.10 10-1.05 10-1.00 10-0.95 10-0.90 10-0.85 10-0.80 10-0.75 10-0.70 10-0.65 10-0.60 10-0.55 10-0.50 10-0.45 10-0.40 10-0.35 10-0.30 10-0.25 10-0.20 10-0.15 10-0.10 10-0.05 100.00 100.05 100.10 100.15 100.20 100.25 100.30 100.35 100.40 100.45 100.50 100.55 100.60 100.65 100.70 100.75 100.80 100.85 100.90 100.95 101.00 10-2 10-1 100 101 10-2.0 10-1.9 10-1.8 10-1.7 10-1.6 10-1.5 10-1.4 10-1.3 10-1.2 10-1.1 10-1.0 10-0.9 10-0.8 10-0.7 10-0.6 10-0.5 10-0.4 10-0.3 10-0.2 10-0.1 100.0 100.1 100.2 100.3 100.4 100.5 100.6 100.7 100.8 100.9 101.0 y MLP分布で発生させた乱数の相補累積分布 Created with Snap

テールの部分が対数スケールで直線になっているのはベキ分布の出現です。 次に論文に従って確率密度関数(PDF)を記述します。

In [5]:

# probability density function 
function pdfMLP(x,α,μ,σ)
    0.5α * exp(α*μ + 0.5α^2σ^2) * x^-(1+α) * erfc(1/√2*(α*σ - (log(x)-μ)/σ))
end

Out[5]:

pdfMLP (generic function with 1 method)

続いて累積分布関数を記述し、ついでに相補累積分布関数も記述します。

In [6]:

# cumulative distribution functioncumulative distribution function
function cdfMLP(x,α,μ,σ)
    u = (log(x) - μ)/(√2*σ)
    a = 0.5*erfc(-u)
    b = 0.5*exp(α*μ+α^2*σ^2*0.5)*x^(-α)
    c = erfc(α*σ/√2 - u)
    a-b*c
end

# complementary cumulative distribution function
function ccdfMLP(x,α,μ,σ)
    1 - cdfMLP(x,α,μ,σ)
end

Out[6]:

ccdfMLP (generic function with 1 method)

ここで描画のための便利関数を記述したら、準備万端です。

In [7]:

# 描画のための便利関数
function makedf(xs,params;fn=ccdfMLP)
    df = DataFrame()
    for (i,j,k) in params
        df = vcat(df,DataFrame(x=xs,y=map(x -> fn(x,i,j,k), xs) , settings="α:$i μ:$j σ:$k"))
    end
    df
end

Out[7]:

makedf (generic function with 1 method)

次にMLP分布の3つのパラメータでいろいろな実験をして、それぞれのパラメータの性質を把握しておきましょう。 数学の得意方は論文の式からしっかりと把握しているかもしれませんが、あえてここはエンジニア的に感覚を掴めます。 まずはαの性質から調べます。

In [8]:

df = makedf(1.0:10.0:1000.0,[(0.5,2.5,0.01), (1.0,2.5,0.01), (1.5,2.5,0.01), (2.0,2.5,0.01)]);
plot(df,x=:x,y=:y,color=:settings,Scale.x_log10,Scale.y_log10,Geom.line,
            Guide.title("相補累積分布"),
            Guide.xlabel("x"))

Out[8]:

x 10-4 10-3 10-2 10-1 100 101 102 103 104 105 106 107 10-3.0 10-2.9 10-2.8 10-2.7 10-2.6 10-2.5 10-2.4 10-2.3 10-2.2 10-2.1 10-2.0 10-1.9 10-1.8 10-1.7 10-1.6 10-1.5 10-1.4 10-1.3 10-1.2 10-1.1 10-1.0 10-0.9 10-0.8 10-0.7 10-0.6 10-0.5 10-0.4 10-0.3 10-0.2 10-0.1 100.0 100.1 100.2 100.3 100.4 100.5 100.6 100.7 100.8 100.9 101.0 101.1 101.2 101.3 101.4 101.5 101.6 101.7 101.8 101.9 102.0 102.1 102.2 102.3 102.4 102.5 102.6 102.7 102.8 102.9 103.0 103.1 103.2 103.3 103.4 103.5 103.6 103.7 103.8 103.9 104.0 104.1 104.2 104.3 104.4 104.5 104.6 104.7 104.8 104.9 105.0 105.1 105.2 105.3 105.4 105.5 105.6 105.7 105.8 105.9 106.0 10-3 100 103 106 10-3.0 10-2.8 10-2.6 10-2.4 10-2.2 10-2.0 10-1.8 10-1.6 10-1.4 10-1.2 10-1.0 10-0.8 10-0.6 10-0.4 10-0.2 100.0 100.2 100.4 100.6 100.8 101.0 101.2 101.4 101.6 101.8 102.0 102.2 102.4 102.6 102.8 103.0 103.2 103.4 103.6 103.8 104.0 104.2 104.4 104.6 104.8 105.0 105.2 105.4 105.6 105.8 106.0 α:0.5 μ:2.5 σ:0.01 α:1.0 μ:2.5 σ:0.01 α:1.5 μ:2.5 σ:0.01 α:2.0 μ:2.5 σ:0.01 settings 10-9 10-8 10-7 10-6 10-5 10-4 10-3 10-2 10-1 100 101 102 103 104 105 10-8.0 10-7.8 10-7.6 10-7.4 10-7.2 10-7.0 10-6.8 10-6.6 10-6.4 10-6.2 10-6.0 10-5.8 10-5.6 10-5.4 10-5.2 10-5.0 10-4.8 10-4.6 10-4.4 10-4.2 10-4.0 10-3.8 10-3.6 10-3.4 10-3.2 10-3.0 10-2.8 10-2.6 10-2.4 10-2.2 10-2.0 10-1.8 10-1.6 10-1.4 10-1.2 10-1.0 10-0.8 10-0.6 10-0.4 10-0.2 100.0 100.2 100.4 100.6 100.8 101.0 101.2 101.4 101.6 101.8 102.0 102.2 102.4 102.6 102.8 103.0 103.2 103.4 103.6 103.8 104.0 10-10 10-5 100 105 10-8.0 10-7.5 10-7.0 10-6.5 10-6.0 10-5.5 10-5.0 10-4.5 10-4.0 10-3.5 10-3.0 10-2.5 10-2.0 10-1.5 10-1.0 10-0.5 100.0 100.5 101.0 101.5 102.0 102.5 103.0 103.5 104.0 y 相補累積分布 Created with Snap

グラフから見ると、αがベキの傾きを決めるようです。値が大きければ大きいほど斜線が急傾斜になる。 次にμの性質を調べます。

In [9]:

df = makedf(1.0:10.0:1000.0,[(1.0,2.5,0.01), (1.0,3.5,0.01), (1.0,4.5,0.01), (1.0,5.5,0.01)]);
plot(df,x=:x,y=:y,color=:settings,Scale.x_log10,Scale.y_log10,Geom.line,
            Guide.title("相補累積分布"),
            Guide.xlabel("x"))

Out [9]:

x 10-4 10-3 10-2 10-1 100 101 102 103 104 105 106 107 10-3.0 10-2.9 10-2.8 10-2.7 10-2.6 10-2.5 10-2.4 10-2.3 10-2.2 10-2.1 10-2.0 10-1.9 10-1.8 10-1.7 10-1.6 10-1.5 10-1.4 10-1.3 10-1.2 10-1.1 10-1.0 10-0.9 10-0.8 10-0.7 10-0.6 10-0.5 10-0.4 10-0.3 10-0.2 10-0.1 100.0 100.1 100.2 100.3 100.4 100.5 100.6 100.7 100.8 100.9 101.0 101.1 101.2 101.3 101.4 101.5 101.6 101.7 101.8 101.9 102.0 102.1 102.2 102.3 102.4 102.5 102.6 102.7 102.8 102.9 103.0 103.1 103.2 103.3 103.4 103.5 103.6 103.7 103.8 103.9 104.0 104.1 104.2 104.3 104.4 104.5 104.6 104.7 104.8 104.9 105.0 105.1 105.2 105.3 105.4 105.5 105.6 105.7 105.8 105.9 106.0 10-3 100 103 106 10-3.0 10-2.8 10-2.6 10-2.4 10-2.2 10-2.0 10-1.8 10-1.6 10-1.4 10-1.2 10-1.0 10-0.8 10-0.6 10-0.4 10-0.2 100.0 100.2 100.4 100.6 100.8 101.0 101.2 101.4 101.6 101.8 102.0 102.2 102.4 102.6 102.8 103.0 103.2 103.4 103.6 103.8 104.0 104.2 104.4 104.6 104.8 105.0 105.2 105.4 105.6 105.8 106.0 α:1.0 μ:2.5 σ:0.01 α:1.0 μ:3.5 σ:0.01 α:1.0 μ:4.5 σ:0.01 α:1.0 μ:5.5 σ:0.01 settings 10-4.5 10-4.0 10-3.5 10-3.0 10-2.5 10-2.0 10-1.5 10-1.0 10-0.5 100.0 100.5 101.0 101.5 102.0 102.5 10-4.0 10-3.9 10-3.8 10-3.7 10-3.6 10-3.5 10-3.4 10-3.3 10-3.2 10-3.1 10-3.0 10-2.9 10-2.8 10-2.7 10-2.6 10-2.5 10-2.4 10-2.3 10-2.2 10-2.1 10-2.0 10-1.9 10-1.8 10-1.7 10-1.6 10-1.5 10-1.4 10-1.3 10-1.2 10-1.1 10-1.0 10-0.9 10-0.8 10-0.7 10-0.6 10-0.5 10-0.4 10-0.3 10-0.2 10-0.1 100.0 100.1 100.2 100.3 100.4 100.5 100.6 100.7 100.8 100.9 101.0 101.1 101.2 101.3 101.4 101.5 101.6 101.7 101.8 101.9 102.0 10-4 10-2 100 102 10-4.0 10-3.8 10-3.6 10-3.4 10-3.2 10-3.0 10-2.8 10-2.6 10-2.4 10-2.2 10-2.0 10-1.8 10-1.6 10-1.4 10-1.2 10-1.0 10-0.8 10-0.6 10-0.4 10-0.2 100.0 100.2 100.4 100.6 100.8 101.0 101.2 101.4 101.6 101.8 102.0 y 相補累積分布 Created with Snap

グラフから見ると、μはベキに切り替える位置を決めるようです。大きければ大きいほど右の方へ移動します。 最後にσの性質を調べます。

In [10]:

df = makedf(1.0:10.0:1000.0,[(1.0,2.5,0.1), (1.0,2.5,0.6), (1.0,2.5,1.0), (1.0,2.5,2.0)]);
plot(df,x=:x,y=:y,color=:settings,Scale.x_log10,Scale.y_log10,Geom.line,
            Guide.title("相補累積分布"),
            Guide.xlabel("x"))

Out [10]:

x 10-4 10-3 10-2 10-1 100 101 102 103 104 105 106 107 10-3.0 10-2.9 10-2.8 10-2.7 10-2.6 10-2.5 10-2.4 10-2.3 10-2.2 10-2.1 10-2.0 10-1.9 10-1.8 10-1.7 10-1.6 10-1.5 10-1.4 10-1.3 10-1.2 10-1.1 10-1.0 10-0.9 10-0.8 10-0.7 10-0.6 10-0.5 10-0.4 10-0.3 10-0.2 10-0.1 100.0 100.1 100.2 100.3 100.4 100.5 100.6 100.7 100.8 100.9 101.0 101.1 101.2 101.3 101.4 101.5 101.6 101.7 101.8 101.9 102.0 102.1 102.2 102.3 102.4 102.5 102.6 102.7 102.8 102.9 103.0 103.1 103.2 103.3 103.4 103.5 103.6 103.7 103.8 103.9 104.0 104.1 104.2 104.3 104.4 104.5 104.6 104.7 104.8 104.9 105.0 105.1 105.2 105.3 105.4 105.5 105.6 105.7 105.8 105.9 106.0 10-3 100 103 106 10-3.0 10-2.8 10-2.6 10-2.4 10-2.2 10-2.0 10-1.8 10-1.6 10-1.4 10-1.2 10-1.0 10-0.8 10-0.6Out[13]: 10-0.4 10-0.2 100.0 100.2 100.4 100.6 100.8 101.0 101.2 101.4 101.6 101.8 102.0 102.2 102.4 102.6 102.8 103.0 103.2 103.4 103.6