散布図からのエラーバー図の作成方法

以下の記事の改良版です.
https://qiita.com/hajkk/items/1580157feecdf7aa562d
https://qiita.com/hajkk/items/da76559c11283b140df4

もともとの図と完成図

左がもとの散布図, 右が完成図です.

データの出典

e-statから取得してきており, 2016年度における日本の市区町村ごとの事業所数(経済センサス-基礎調査結果, 単位:個)と従業者数(経済センサス-基礎調査結果, 単位:人)が記録されています. 先ほどの図において, 横軸に事業所数, 縦軸に従業者数を図示しています.

作り方

データはあらかじめ次の形で与えられているとします.

x	y
118718	1477012
36655	532655
14312	174152
14224	172162
13145	159039

作成コードは次. おまけとして, 右図で横軸と縦軸の中央値の両対数fittingのプロセスも書いています. 傾きは1.080(28), $R^2$は0.972でした. 以上です.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm

def plot_scatter(df: pd.DataFrame):
  fig, ax = plt.subplots(figsize=(4, 4))
  ax.set_xscale('log')
  ax.set_yscale('log')
  ax.set_xlim([10**0, 10**6])
  ax.set_ylim([10**0, 10**7])
  ax.scatter(df['x'], df['y'], fc='none', ec='b', s=12)
  fig.savefig('left.png')

def plot_errorbar(df: pd.DataFrame):
  xbw = 0.1
  df['xbin'] = 10**(np.round(np.log10(df['x']) / xbw) * xbw)
  summary = df.groupby('xbin').quantile([0.25, 0.50, 0.75])\
              .reset_index().pivot(index='xbin', columns='level_1', values='y')
  fig, ax = plt.subplots(figsize=(4, 4))
  ax.set_xscale('log')
  ax.set_yscale('log')
  ax.set_xlim([10**0, 10**6])
  ax.set_ylim([10**0, 10**7])
  ax.errorbar(summary.index, summary[0.50], yerr=[summary[0.50]-summary[0.25], summary[0.75]-summary[0.50]],\
              fmt='o', c='b', mfc='none', lw=1, capsize=4, label='Errorbar: 25, 50, 75%')
  ax.legend(loc=2, fontsize=12)
  fig.savefig('right.png')

 def fit():
    fitdata = np.log(summary.reset_index())
    fitdata[~np.isfinite(fitdata)] = np.nan
    fitdata.dropna(inplace=True)
    X = sm.add_constant(fitdata['xbin'])
    y = fitdata[0.50]
    res = sm.OLS(y, X).fit()
    print(res.summary())

  fit()
  

if __name__ == '__main__':
  df = pd.read_csv('data.csv')
  df.columns = ['x', 'y']
  plot_scatter(df)
  plot_errorbar(df)