動機

個人所得や都市の人口など、確率変数が大きな値と小さな値を取るときに対数変換を経由して確率密度関数を求めることができる。

対数変換による確率密度関数の変換公式

変数変換$y = \log x$にともなう確率変数$X$から$Y$への変換により、$X$の確率密度関数(PDF) $p(x)$は$Y$の確率密度関数$q(y)$に次のように変換される。

p(x) = q(y) /x

導出

これは次のように導出される。累積分布関数を

P(X \leq x) = \int_{-\infty}^{x} p(t)dt

P(Y \leq y) = \int_{-\infty}^{y} q(s)ds

で定義する。変数変換$y = \log x, Y = \log X$により

P(X \leq x) = P(Y \leq \log x)

なので

q(y) = \frac{dP(Y \leq y)}{dy} = \frac{dP(X \leq x)}{dx} \frac{dx}{dy} = p(x) x

となる。

Python code

次の図を出力するpythonコードは以下。パレート分布(べき分布)から100000万点サンプルしてPDFを計算した。(青)対数空間上で等間隔にbinを区切りリニアスケールで直接PDFを計算したもの、(赤)各点を対数変換して対数空間上でPDFを計算してリニアスケールに戻したもの、そして(黒)解析解を描画した。

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

def pdfs(series, bins, ax=None, scale='linear', kwargs={}):
  cnts = pd.cut(series, bins).value_counts()
  p = cnts / series.shape[0]
  vals = p / np.diff(bins)
  if scale=='linear':
    mb = np.array([ (bins[i]+bins[i+1])/2. for i in np.arange(bins.shape[0]-1) ])
  elif scale=='log':
    mb = np.array([ np.sqrt(bins[i]*bins[i+1]) for i in np.arange(bins.shape[0]-1) ])
  else:
    raise ValueError("scale should be either 'linear' or 'log'")
  if ax!=None:
    ax.plot(mb, vals, **kwargs)
  return mb, vals, cnts, p

def analytical_pdf_pareto(a, xmin, x):
  return a/xmin * x**(-(a+1))

if __name__ == '__main__':

  a, xmin = 1.0, 1.0
  sample_size = 100000
  
  series = (np.random.pareto(a, sample_size) + 1) * xmin
  ln_transformed_series = np.log(series)

  ### calculate bins 
  ln_bin_width = 0.1

  ln_transformed_bins = np.arange(ln_transformed_series.min(), ln_transformed_series.max(), bin_width)
  bins = np.exp(ln_transformed_bins)

  ### calculate PDF of original series
  original_bin, original_pdf, _, __ = pdfs(series = series, bins = bins)
  
  ### calculate PDF of log-transformed series
  ln_bin, ln_pdf, _, __ = pdfs(series = ln_transformed_series, bins = ln_transformed_bins)

  ### plot PDF
  fig, ax = plt.subplots()
  ax.set_xlabel(r'$x$')
  ax.set_ylabel(r'$p(x)$')
  ax.set_xscale('log')
  ax.set_yscale('log')
  
  ax.plot(original_bin, original_pdf, 
  c='b', marker='o', label='PDF of original series')
  
  ax.plot(original_bin, ln_pdf / original_bin, 
  c='r', marker='x', label='PDF via log-transformation')
  
  ax.plot(original_bin, analytical_pdf_pareto(a, xmin, original_bin),
  c='k', ls='--', label='PDF Analytical')

  ax.legend(loc='lower left')
  fig.savefig('probability_distribution.png')