「Pythonで脂肪尾を数値化する4つの方法」

「脂肪尾を数値化するための4つのPythonの方法」

直感と例のコード

これは、パワーローとファットテールシリーズの3番目の記事です。 前の投稿では、実データからパワーローを検出する方法について説明しました。このテクニックは便利ですが、ファットテールは単にデータをパワーロー分布に適合させる以上のものです。この記事では、ファットテールを数量化する4つの方法を解説し、実世界のデータを分析するPythonの例コードを共有します。

注意:パワーロー分布やファットテールといった用語に馴染みがない場合は、この記事を参照してください。

ファット(猫の)テール。Canvaからの画像。

このシリーズの最初の記事では、ファットテールの概念を紹介しました。それは、希少なイベントが分布の集計統計を駆動する度合いを表しています。パレート分布を例に挙げると、たとえば売上の80%が顧客の20%から生成され(売上の50%はたったの1%の顧客から生成されます)。

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

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

towardsdatascience.com

パレート(そして一般的にはパワーロー)分布はファットテールの顕著な例を提供していますが、これはより一般的な概念であり、ガウス分布のような細い尾から非常にファットテールなテール(パレート80-20など)までのスペクトル上に存在しています。

ファットテール度合いのスペクトル。著者による画像。

このファットテール度合いの視点は、データを単にパワーロー(またはそうでない)としてラベル付けするよりも柔軟で正確な分類方法を提供してくれます。ただし、ここで疑問が生じます:どのようにファットテール度合いを定義するのでしょうか?

ファットテールを定量化する4つの方法

ファットテールの「真の」測定方法はありませんが、実際のデータのファットテール性を定量化するためのいくつかのヒューリスティックス(つまり、経験則)があります。ここでは、4つのヒューリスティックスを概念的に紹介し、その後にPythonの例コードに入ります。

ヒューリスティックス1:パワーロー尾指数

最もファットなテールはパワーロー分布に現れ、パワーローの尾指数(つまり、α)が小さいほど、テールはファットです。以下の画像で示されています。

さまざまなα値を持つパワーロー分布の例。著者による画像。

小さい尾指数はファットテールを意味するというこの観察から、我々はファットテールを定量化するためにαを使用する動機付けをしたります。実践では、与えられたデータセットにパワーロー分布を適合させ、推定されたα値を抽出することに帰結します。

これは直接的なアプローチですが、明確な制限があります。つまり、パワーローウにほとんどフィットしないデータを扱う場合には、このアプローチはうまく機能しません。

ヒューリスティック2:尖度(つまり非ガウス性)

ファットテールの対極は、スリムなテール(つまり、稀なイベントが極めて稀であるということ)です。スリムテールの代表例は、愛されるガウス分布です。ガウス分布では、平均から6シグマ離れたイベントの確率は約10億分の1です。

これは、データがどれだけ「非ガウス的」であるかを定量化することで、ファットテールの別の指標を示唆しています。これは、いくつもの非ガウシアン性の指標を考案できますが、最も一般的なものは尖度 (Kurtosis)です。以下の式で定義されています。

Definition of Kurtosis according to ref [1] and [2]. Image by author.

尖度は中心から遠く(つまりテールから)の値によって決定されます。したがって、尖度が大きいほど、テールが太くなります

この指標は、すべてのモーメントが有限である場合にはうまく機能します(参考文献[3])。ただし、尖度は定義されていない分布(例:α≦4のパレート分布)では利用できないため、多くのファットテールデータには無意味です。

ヒューリスティック3:対数正規分布のσ

このシリーズの以前の記事では、以下の確率密度関数で定義される対数正規分布について説明しました。

The probability density function of log-normal distribution [4]. Image by author.

先に述べたように、この分布は少ないσではガウス風に見えますが、高いσではパレート風に見えることがあります。これは、テールが太いかどうかを定量化する別の方法を自然に提供します。つまり、σが大きいほど、テールが太くなります

この指標は、ヒューリスティック1と同様の方法で得ることができます。すなわち、データに対して対数正規分布をフィットさせ、フィットのσ値を抽出します。これは簡単な手順ですが(ヒューリスティック1と同様に)、対数正規分布が元のデータを適切に説明しない場合にはうまく機能しません。

ヒューリスティック4:タレブのκ

前述のヒューリスティック(H)は、特定の分布を前提としています(つまり、H1 – パワーローとH3 – 対数正規分布)。しかし、実際には、私たちのデータはほとんど正確に特定の分布に従うわけではありません。

さらに、これらの指標を使用して、質的に異なる分布に従う2つの変数を評価する場合には問題が生じる場合があります。たとえば、パレト様とガウス様のデータを比較するためにパワーローのテール指数を使用する場合、パワーローはガウス様のデータには適合しないため、ほとんど意味がありません。

そこで、分布に依存しないファットテールの指標の使用が動機付けられます。そのような指標の1つは、Talebによって参考文献[3]で提案されました。提案された指標(κ)は、単峰性のデータ(有限な平均を持つ)について0から1の値を取るもので、0はデータが最大限スリムテールであることを示し、1はデータが最大限ファットテールであることを示します。以下の式で定義されています。

定義 Taleb's κ メトリック [3]。著者による画像。

このメトリックは、2つのサンプル(Sₙ₀とSₙとする)を比較します。ここで、Sₙは特定の分布から抽出されたn個のサンプルの合計です。たとえば、ガウス分布を評価し、n=100とする場合、ガウス分布から100個のサンプルを抽出し、それらをすべて合計してS₁₀₀を作ります。

上記の式中のM(n)は、以下の式に従って定義された平均絶対偏差を示しています。この平均絶対偏差は平均周りの分布のバラツキを示し、通常、標準偏差よりも堅牢性があります [3][5]。

κ 方程式による平均絶対偏差の定義 [3]。著者による画像。

簡単にするために、n₀=1として、以下の式を得ることができます。

n₀=1 の κ。著者による画像。

ここでの重要な用語は、M(n)/M(1)です。ここで、M(n)はn個のサンプルの合計の平均周りの分布のバラツキを量化します(ある分布の)。

テールが薄いデータの場合、M(30)はM(1)に比べてかなり近くなるでしょう。なぜならば、データは一般的に平均に近い位置に存在しているからです。したがって、M(30)/M(1) ≈ 1 となります。

しかし、テールが太いデータの場合、M(30)はM(1)よりもはるかに大きくなります。したがって、M(30)/M(1) >> 1 となります。これは下の図に示されており、左のプロットはガウス分布の合計における分散のスケーリングを示し、右のプロットは Pareto 分布におけるスケーリングを示しています。

ガウス分布(左)および Pareto 80-20 分布(右)における M(n) と M(1) のスケーリング。y軸のラベルに注意してください。注意:ガウス分布の分散のスケーリングは、n個の分布の合計により増加します。著者による画像。

したがって、テールが太いデータの場合、κ 方程式の分母は分子よりも大きくなるため、RHSの第2項は小さくなり、結果としてκが大きくなります。

これがあなたが求めていた以上に数学的な内容であれば、次のポイントをお伝えします:大きな κ = テールが太い、小さな κ = テールが薄いです。

コード例:(現実世界の)ソーシャルメディアデータのテールが太い性質の定量化

理論的な内容が終わったので、これらのヒューリスティックスを実践でどのように使用するか見てみましょう。ここでは、このシリーズの前の記事からのデータを分析するために、それぞれのアプローチを使用します。

これらのデータは私のソーシャルメディアアカウントのもので、VoAGIでの月間のフォロワー増加数、YouTube動画ごとの収益、LinkedInでの日次インプレッションが含まれています。データとコードはGitHub リポジトリで自由に利用できます。

YouTube-Blog/power-laws/3-尾の太い分布の数量化

YouTubeのビデオとブログ記事を補完するためのコード – YouTube-Blog/power-laws/3-尾の太い分布の数量化

github.com

いくつかの便利なライブラリをインポートして始めましょう。

import matplotlib.pyplot as pltimport pandas as pdimport numpy as npimport powerlawfrom scipy.stats import kurtosis

次に、各データセットを読み込んで辞書に保存します。

filename_list = ['VoAGI-followers', 'YT-earnings', 'LI-impressions']df_dict = {}for filename in filename_list:    df = pd.read_csv('data/'+filename+'.csv')    df = df.set_index(df.columns[0]) # インデックスを設定    df_dict[filename] = df

ここで、データを見ることは常に良いアイデアです。ヒストグラムをプロットし、各データセットの上位5つのレコードを表示してみましょう。

for filename in filename_list:    df = df_dict[filename]        # ヒストグラムをプロットする(関数はGitHubのノートブックに定義されています)    plot_histograms(df.iloc[:,0][df.iloc[:,0]>0], filename, filename.split('-')[1])    plt.savefig("images/"+filename+"_histograms.png")        # 上位5つのレコードを表示する    print("Top 5 Records by Percentage")    print((df.iloc[:,0]/df.iloc[:,0].sum()).sort_values(ascending=False)[:5])    print("")
月次のVoAGIフォロワーのヒストグラム。画像は著者によるものです。
YouTubeビデオの収益のヒストグラム。注:以前の記事と異なる点に気付く場合は、データに迷えるレコードがあるためです(だから見ることの重要性です😅)。画像は著者によるものです。
LinkedInの印象数の日次ヒストグラム。画像は著者によるものです。

上記のヒストグラムに基づくと、各データセットはいくらか尾部が太いように見えます。パーセンテージ上位5つのレコードを見て、他の角度から見てみましょう。

各データセットの上位5つのパーセンテージレコード。画像は著者によるものです。

この表示から、VoAGIフォロワーは最も尾部が太く、フォロワーの約60%がわずか2か月で得られました。YouTubeの収益も強く尾部が太く、収益の約60%がわずか4つのビデオから得られています。LinkedInの印象数は最も尾部が太くないようです。

データを見るだけでは尾部の太さを qualitatively な感じることができますが、この良くなはったデータを quantitatively に評価してみましょう。

ヒューリスティック1:べき乗法の尾部指数

各データセットごとにαを取得するために、前の記事で行ったようにpowerlawライブラリを使用することができます。以下のコードブロックで、各データセットに対して適合を実行し、パラメータ推定値をforループで表示します。

filename_listの各ファイル名に対して:    df = df_dict[filename]    # パワーロー適合を実行    results = powerlaw.Fit(df.iloc[:,0])    # 結果を表示    print("")    print(filename)    print("-"*len(filename))    print("パワーロー適合")    print("a = " + str(results.power_law.alpha-1))    print("xmin = " + str(results.power_law.xmin))    print("")
パワーロー適合結果。画像は著者提供。

上記の結果は、VoAGIのフォロワーが最もfat-tailedであり、次にYouTubeの収益とLinkedInのインプレッションが続くという、私たちの質的評価と一致しています(ここで、小さなαはよりfat-tailedを意味します)。

ヒューリスティク2:尖度

オフシェルフの実装を使用して尖度を計算する簡単な方法があります。ここでは、Scipyを使用して結果を前と同様の方法で表示します。

filename_listの各ファイル名に対して:    df = df_dict[filename]    # 結果を表示    print(filename)    print("-"*len(filename))    print("尖度 = " + str(kurtosis(df.iloc[:,0], fisher=True)))    print("")
各データセットの尖度値。画像は著者提供。

尖度はヒュリスティック1とは異なる結果を示しています。この尺度によるfat-tailednessの順位付けは、LinkedIn > VoAGI > YouTubeです。

ただし、これらの結果は慎重に扱う必要があります。上記のパワーロー適合の結果で見たように、すべての3つのデータセットがα < 4のパワーローに適合し、つまり尖度は無限大です。したがって、計算は値を返しますが、これらの数値には疑念を抱くのが賢明です。

ヒューリスティック3:対数正規分布のσ

ヒューリスティック1と同様に、powerlawライブラリを使用してσの推定値を取得することもできます。以下はその様子です。

filename_listの各ファイル名に対して:    df = df_dict[filename]    # パワーロー適合を実行    results = powerlaw.Fit(df.iloc[:,0])    # 結果を表示    print("")    print(filename)    print("-"*len(filename))    print("対数正規分布適合")    print("mu = " + str(results.lognormal.mu))    print("sigma = " + str(results.lognormal.sigma))    print("")
対数正規分布の適合結果。画像は著者提供。

上記のσの値を見ると、すべての適合結果がデータがfat-tailedであることを示しており、VoAGIのフォロワーやLinkedInのインプレッションのσの推定値が似ています。一方、YouTubeの収益は、はるかに大きなσ値を持っており、(はるかに)fat-tailedを意味しています。

ただし、疑わしい点の1つは、適合結果で負のμを推定していることであり、これは対数正規分布の適合がデータを十分に説明できない可能性を示しています。

ヒューリスティック4:タレブのκ

                                                                  

                                                                                                                                                            

                                                                                                                                        

                                                                                                                                                                                   

                                                                                                                                                         

                                                                                                                                                                                                                                    

Since I couldn’t find an off-the-shelf Python implementation for computing κ (I didn’t look very hard), this computation requires a few extra steps. Namely, we need to define 3 helper functions, as shown below.

def mean_abs_deviation(S):    """        Computation of mean absolute deviation of an input sample S    """    M = np.mean(np.abs(S - np.mean(S)))    return Mdef generate_n_sample(X,n):    """        Function to generate n random samples of size len(X) from an array X    """    # initialize sample    S_n=0        for i in range(n):        # ramdomly sample len(X) observations from X and add it to the sample        S_n = S_n + X[np.random.randint(len(X), size=int(np.round(len(X))))]        return S_ndef kappa(X,n):    """        Taleb's kappa metric from n0=1 as described here: https://arxiv.org/abs/1802.05495        Note: K_1n = kappa(1,n) = 2 - ((log(n)-log(1))/log(M_n/M_1)), where M_n denotes the mean absolute deviation of the sum of n random samples    """    S_1 = X    S_n = generate_n_sample(X,n)        M_1 = mean_abs_deviation(S_1)    M_n = mean_abs_deviation(S_n)        K_1n = 2 - (np.log(n)/np.log(M_n/M_1))    return K_1n

The first function, mean_abs_deviation(), computes the mean absolute deviation as defined earlier.

Next, we need a way to generate and sum n samples from our empirical data. Here, I take a naive approach and randomly sample an input array (X) n times and sum the samples together.

Finally, I bring together mean_abs_deviation(S) and generate_n_sample(X,n) to implement the κ calculation defined before and compute it for each dataset.

n = 100 # number of samples to include in kappa calculationfor filename in filename_list:    df = df_dict[filename]    # print results    print(filename)    print("-"*len(filename))    print("kappa_1n = " + str(kappa(df.iloc[:,0].to_numpy(), n)))    print("")
κ(1,100) values for each dataset. Image by author.

The results above give us yet another story. However, given the implicit randomness of this calculation (recall generate_n_sample() definition) and the fact we’re dealing with fat tails, point estimates (i.e. just running the computation once) cannot be trusted.

Accordingly, I run the same calculation 1000x and print the mean κ(1,100) for each dataset.

num_runs = 1_000kappa_dict = {}for filename in filename_list:    df = df_dict[filename]    kappa_list = []    for i in range(num_runs):        kappa_list.append(kappa(df.iloc[:,0].to_numpy(), n))        kappa_dict[filename] = np.array(kappa_list)    print(filename)    print("-"*len(filename))    print("mean kappa_1n = " + str(np.mean(kappa_dict[filename])))    print("")
Mean κ(1,100) values from 1000 runs for each dataset. Image by author.

These more stable results indicate VoAGI followers are the most fat-tailed, followed by LinkedIn Impressions and YouTube earnings.

Note: One can compare these values to Table III in ref [3] to better understand each κ value. Namely, these values are comparable to a Pareto distribution with α between 2 and 3.

Although each heuristic told a slightly different story, all signs point toward VoAGI followers gained being the most fat-tailed of the 3 datasets.

結論

バイナリラベルを使ってデータをfat-tailed(またはそうでない)と分類することは誘惑されるかもしれませんが、fat-tailednessはスペクトル上に存在します。ここでは、fat-tailedデータを数量化するための4つのヒューリスティックスを紹介しました。

それぞれのアプローチには制限がありますが、実データのfat-tailednessを比較するための定量的な手法を提供しています。

👉 パワーローズとfat tailsについて詳しく知るイントロダクション | パワーローズフィット

リソース

コネクト私のウェブサイト | 電話予約 | 私に何でも聞いてください

ソーシャルメディアYouTube 🎥 | LinkedIn | Twitter

サポートコーヒーを買って応援する ☕️

The Data Entrepreneurs

データエンジニアのためのコミュニティ。👉 Discordに参加してください!

VoAGI.com

[1] Scipy Kurtosis

[2] Scipy Moment

[3] arXiv:1802.05495 [stat.ME]

[4] https://en.wikipedia.org/wiki/Log-normal_distribution

[5] Pham-Gia, T., & Hung, T. (2001). The mean and median absolute deviations. Mathematical and Computer Modelling, 34(7–8), 921–936. https://doi.org/10.1016/S0895-7177(01)00109-1

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

Share:

Was this article helpful?

93 out of 132 found this helpful

Discover more