「Pythonを使って現実世界のデータにおけるべき乗則の検出」

「Pythonを使って現実世界のデータのべき乗則を検出する方法」

最尤法に基づくアプローチの解説とサンプルコード

これはパワーローとファットテールに関するシリーズの2番目の記事です。前回の投稿では、パワーローに関する初心者向けの紹介と、標準統計ツールでの問題点を3つ紹介しました。これらの問題を回避するための認識は役立ちますが、実際には与えられたデータがどのような分布に従っているのかは常に明確ではありません。この記事では、実世界のデータからパワーローを客観的に検出する方法と、ソーシャルメディアデータの具体的な例を共有します。

注：パワーロー分布やファットテールなどの用語に慣れていない場合は、このシリーズの最初の記事をじっくり読んでください。

パワーローは統計学の基本を覆す

前の記事では、ガウシアン分布とパワーロー分布の2種類の分布に焦点を当てました。これらの分布は、統計的な性質がまったく異なることがわかりました。具体的には、パワーローはまれなイベントによって駆動される一方、ガウシアンはそうではありません。

このまれなイベントによって駆動される性質により、パワーローはお気に入りの統計ツール（平均、標準偏差、回帰など）で3つの問題を引き起こします。結論として、データがガウシアンに似ている場合は、回帰や期待値の計算など、一般的なアプローチを心配せずに使用できます。ただし、データがよりパワーローのような形状をしている場合は、これらの手法は間違った結果や誤解を引き起こす可能性があります。

• 現代のサイバーセキュリティの脅威に対処する効果的な方法
• 最適化アルゴリズム：ニューラルネットワークの基礎101
• 『ブンブンの向こう側 産業における生成型AIの実用的な応用を探求する』

さらに、ガウシアンとパワーローの両方に似ている（それにもかかわらず性質が逆）対数正規分布と呼ばれる3番目の分布も見ました。

この曖昧さは、与えられたデータセットを分析する最良の方法を決定するための実践者に課題を提供します。これらの課題を克服するためには、データがパワーロー、対数正規分布、またはその他のタイプの分布に適合するかどうかを判断することが有利です。

対数対数アプローチ

実世界のデータにパワーローを適合させる人気のある方法は、私が「対数対数アプローチ」と呼ぶものです[1]。このアイデアは、下記に示すように、パワーローの確率密度関数（PDF）の対数をとることから派生しています。

上記の導出により、パワーローのPDFの定義を線形方程式として表すことができます。以下に示す図に示しています。

これにより、べき乗則に従うデータのヒストグラムは、直線に従うことが示唆されます。実際のところ、これはデータのヒストグラムを生成し、それを対数-対数プロットにプロットすることで確認できます[1]。さらに進んで、線形回帰を実行して分布のα値（ここでは、α = -m + 1）を推定することもできます。

ただし、このアプローチには重要な制限があります。これらは参考文献[1]で詳しく説明されており、以下にまとめられています。

• 傾き（したがってα）の推定は系統誤差の影響を受ける
• 回帰誤差の推定が難しい場合がある
• 分布がべき乗則に従っていなくてもフィットは良く見える場合がある
• フィットが確率分布の基本的な条件（例：正規化）に従っていない場合がある

最尤推定法

対数-対数法は簡単に実装できますが、その制限から考えると最適とは言えません。代わりに、より数学的に妥当なアプローチである最尤推定法を用いることができます。これは、与えられたデータに対して最適なパラメータを推測するための広く使用されている統計的な方法です。

最尤推定法は2つの主要なステップで構成されています。 ステップ1：尤度関数の取得。 ステップ2：モデルパラメータに関する尤度を最大化します。

ステップ1：尤度関数の作成

尤度は特別なタイプの確率です。単純に言えば、これは特定のモデルでのデータの確率を定量化するものです。すべての観測データについての共同確率として表すことができます[3]。パレート分布の場合、次のように表すことができます。

尤度を最大化することを少し簡単にするために、対数尤度で作業するのが一般的です（それらは同じα値によって最大化されます）。

ステップ2：尤度の最大化

(対数)尤度関数が手元にある状態で、パラメータの最適な選択を決定するタスクを最適化問題としてフレーム化できます。データに基づいて最適なα値を見つけるため、これは、l(α)のαに関する導関数をゼロに設定し、αについて解くということになります。これに関する導出は以下の通りです。

これにより、いわゆる最大尤度推定が得られます。これにより、観測されたxの値を代入してパレート分布のα値を推定することができます。

理論的な基盤が整ったので、実際のデータ（私のソーシャルメディアアカウントからのもの）に適用するとどのように見えるかを見てみましょう。

例: ソーシャルメディアデータへのべき乗則の適合

fat-tailed データが一般に存在するドメインの1つがソーシャルメディアです。たとえば、作成者の一部が注目を集め、VoAGIのブログの少数が読者の大多数を占めるなどです。

ここでは、powerlaw Pythonライブラリを使用して、私のさまざまなソーシャルメディアチャンネル（つまり、VoAGI、YouTube、LinkedIn）のデータが本当にべき乗則分布に従っているかどうかを判断します。これらの例のデータとコードは、GitHubリポジトリで入手できます。

YouTube-Blog/power-laws/2-detecting-powerlaws at main · ShawhinT/YouTube-Blog

VoAGIに関するYouTubeビデオおよびブログ投稿の補完コード – YouTube-Blog/power-laws/2-detecting-powerlaws at main ·…

github.com

人工データ

現実のごちゃごちゃしたデータに最大尤度ベースの手法を適用する前に、パレート分布と対数正規分布から生成された（真の）人工データにこの手法を適用してみましょう。これによって、私たちが「真の」潜在的な分布クラスを正確に知らないデータにこの手法を使用する前に、私たちの予想を具体化します。

まず、いくつかの便利なライブラリをインポートします。

import numpy as np
import matplotlib.pyplot as plt
import powerlaw
import pandas as pd
np.random.seed(0)

次に、パレート分布と対数正規分布からデータを生成しましょう。

# パワーローデータa = 2
x_min = 1
n = 1_000
x = np.linspace(0, n, n+1)
s_pareto = (np.random.pareto(a, len(x)) + 1) * x_min
# ログ正規データm = 10
s = 1
s_lognormal = np.random.lognormal(m, s, len(x)) * s * np.sqrt(2*np.pi)

これらのデータがどのように見えるかを把握するには、ヒストグラムをプロットするのが役立ちます。ここでは、各サンプルの生の値とその対数のヒストグラムをプロットします。後者の分布は、パワーローと対数正規データを視覚的に区別しやすくします。

上記のヒストグラムからわかるように、生の値の分布は両方の分布で qualitatively 似ています。しかし、対数分布には顕著な違いがあります。つまり、対数パワーロー分布は非常に偏っており、平均に対して中央寄りではない一方、対数の対数正規分布はガウス分布に似ています。

さて、powerlawライブラリを使用して、各サンプルにパワーローを適合させ、αとx_minの値を推定しましょう。以下は、パワーローサンプルに対してその様子です。

# パワーローの適合を適用するresults = powerlaw.Fit(s_pareto)# 結果を出力print("alpha = " + str(results.power_law.alpha)) # 注: powerlawライブラリのalpha定義は、標準の定義（a_standard + 1）とは異なりますprint("x_min = " + str(results.power_law.xmin))print('p = ' + str(compute_power_law_p_val(results)))# パワーロー適合のための最適な最小値を計算# alpha = 2.9331912195958676# x_min = 1.2703447024073973# p = 0.999

適切なパラメータ値（つまり、a=3、x_min=1）を推定するため、適合度はまずまずの仕事をします。上記の alpha と x_min の値によって確認できます。上記の p 値は、適合度の質を定量化します。p 値が高いほど、適合度が良くなります（この値については、参考文献[1]のセクション4.1で詳しく説明されています）。

同様のことをログ正規分布についても行うことができます。

# ログ正規データにパワーローを適用results = powerlaw.Fit(s_lognormal)print("alpha = " + str(results.power_law.alpha)) # 注意：パワーローのalphaの定義は、標準のalphaの定義と異なります (a_powerlawlib = a_standard + 1)print("x_min = " + str(results.power_law.xmin))print('p = ' + str(compute_power_law_p_val(results)))# パワーローフィットのための最適な最小値の計算# alpha = 2.5508694755027337# x_min = 76574.4701482522# p = 0.999

ログ正規分布のサンプルもまた、パワーロー分布によく適合していることが分かります（p=0.999）。ただし、x_minの値はテール部にあります。これは一部の使用例では役に立つかもしれませんが、サンプル内のすべてのデータに最も適合する分布についてはあまりわかりません。

これを克服するために、x_minの値をサンプルの最小値に手動で設定してフィットをやり直すことができます。

# x_minを修正してフィットにすべてのデータを含めるresults = powerlaw.Fit(s_lognormal, xmin=np.min(s_lognormal))print("alpha = " + str(results.power_law.alpha))print("x_min = " + str(results.power_law.xmin))# alpha = 1.3087955873576855# x_min = 2201.318351239509

.Fit()メソッドは、自動的にログ正規分布の推定値も生成します。

print("mu = " + str(results.lognormal.mu))print("sigma = " + str(results.lognormal.sigma))# mu = 10.933481999687547# sigma = 0.9834599169175509

推定されたログ正規分布のパラメータ値は、実際の値（mu=10、sigma=1）に近く、フィットは再びうまくいきました！

ただし、x_minを修正することで、品質指標pを失いました（何らかの理由で、x_minが指定された場合に値が生成されません）。これが疑問を投げかけることです。どちらの分布パラメータを選べばよいのでしょうか？パワーローまたはログ正規分布ですか？

この質問に答えるために、対数尤度比（R）を使用してパワーローのフィットを他の候補分布と比較することができます。正のRはパワーローがより適合していることを示し、負のRは代替分布がより適合していることを示します。さらに、各比較には有意性値（p）が得られます。これは以下のコードブロックで示されています。

distribution_list = ['lognormal', 'exponential', 'truncated_power_law', \      'stretched_exponential', 'lognormal_positive']    for distribution in distribution_list:    R, p = results.distribution_compare('power_law', distribution)    print("power law vs " + distribution +           ": R = " + str(np.round(R,3)) +           ", p = " + str(np.round(p,3)))# power law vs lognormal: R = -776.987, p = 0.0# power law vs exponential: R = -737.24, p = 0.0# power law vs truncated_power_law: R = -419.958, p = 0.0# power law vs stretched_exponential: R = -737.289, p = 0.0# power law vs lognormal_positive: R = -776.987, p = 0.0

上記のように、すべての代替分布はログ正規サンプルのすべてのデータを含めた場合にパワーローに対して優先されます。また、尤度比に基づくと、lognormalとlognormal_positiveのフィットが最も良く機能します。

実世界のデータ

真の分布がわからないデータにpowerlawライブラリを適用してみましょう。

ここでは、次のデータを分析します。私のVoAGIプロフィールの月間フォロワー獲得数、私のすべてのYouTubeビデオの収益、および過去1年間の私のLinkedIn投稿の日ごとのインプレッション。

まず、ヒストグラムをプロットします。

これらのプロットから私が把握した2つのことです。1つは、これまで見たパワーローヒストグラムよりも、これらの3つはすべて対数正規ヒストグラムに似ているということです。 2つは、VoAGIとYouTubeの分布がまばらであることであり、強力な結論を導くための十分なデータが不足している可能性があることを意味しています。

次に、x_minを各サンプルの最小値として、3つの分布全てにパワーローフィットを適用します。その結果は以下に出力されます。

どの分布が最適であるかを判断するために、パワーローフィットといくつかの代替分布を比較することができます。以下にその結果が示されています。

有意性のカットオフポイントであるp<0.1を使用して、以下の結論を導くことができます。VoAGIのフォロワーとLinkedInの表示回数は、対数正規分布に最も適合し、YouTubeの収益をパワーローが最も表しています。

もちろん、ここでのVoAGIフォロワーやYouTubeの収益データは限られています（N<100）ので、それらのデータから得られる結論は慎重に扱う必要があります。

結論

パワーロー分布に従うデータに対しては、多くの標準的な統計ツールが使えません。したがって、実データでのパワーローの検出は、誤った分析や誤解を防ぐのに役立ちます。

ただし、パワーローはより一般的な現象であるファットテールの極端な例です。このシリーズの次の記事では、この作業をさらに進めて、任意のデータセットに対してファットテールを定量化するための4つの便利なヒューリスティックスを紹介します。

パレート、パワーロー、ファットテール

統計学では教えてくれないこと

towardsdatascience.com

リソース

連絡先：私のウェブサイト | 電話で予約する | 何でも聞いてください

ソーシャル：YouTube 🎥 | LinkedIn | Twitter

サポート：コーヒーを買ってください ☕️

データエンタープライズ

データ領域の起業家コミュニティ。👉 Discordに参加する！

VoAGI.com

[1] arXiv:0706.1062 [physics.data-an]

[2] arXiv:2001.10488 [stat.OT]

[3] https://en.wikipedia.org/wiki/Likelihood_function

[4] https://en.wikipedia.org/wiki/Pareto_distribution

We will continue to update VoAGI; if you have any questions or suggestions, please contact us!