自転車シェアリングの貸出・返却のモデル化 ― ポアソン過程・指数分布

ポアソン過程
ポアソン過程確率分布統計自転車シェアリング
0

自転車シェアリングの利用パターンを理解するために、利用者が貸出や返却するといったイベントの発生間隔がどのような分布に従うのか確かめたい。

そこで、ニューヨークで自転車シェアリング事業を展開するCitiBikeのオープンデータを用いて、簡単に調べてみたのでその結果をメモする。

スポンサーリンク

ポアソン過程と指数分布

貸出や返却のようなイベントがランダムに発生する過程は「ポアソン過程」という。
そして、ポアソン過程に従う貸出や返却イベントの間隔は「指数分布」に従う。

指数分布がどういうものかに関しては、次の記事にまとめた。

ポアソン過程がランダムに発生するイベントであることを考えると明らかであるが、指数分布は無記憶の分布と言われる。
無記憶というのは過去のデータに依存せずに、次の人が来る間隔が決定されるということである。
自転車シェアリングでは、当然借りる人は過去に借りた人の情報を知らないはずである。(知人と一緒に借りるという場合があるかもしれないが、そのような場合もまたランダム性の中に埋もれていくはずである。)

スポンサーリンク

用いたデータ

実際に貸出間隔が指数分布に従うのかCiti Bikeが公開するデータを用いて検証する。
Citi BikeはNew York Cityで自転車シェアリング事業を展開しており、2020年5月時点で14500台の自転車と約900のステーションがある。

次のようにデータを準備する。(今回は2019年9月のデータを使った。)
ファイルのパスは適宜変更してください。
また、描画ツールとしてPlotlyを使ったが好きなものに変更してください。

import numpy as np
from scipy.stats import expon
import pandas as pd
from plotly.subplots import make_subplots
import plotly.graph_objects as go

year = '2019'
month = '09'

df = pd.read_csv('../data/Citi Bike Trip Data/' + year + month +'-citibike-tripdata.csv')

# 欠損値を削除する
df = df.dropna()

df['starttime'] = pd.to_datetime(df['starttime'])
df['stoptime'] = pd.to_datetime(df['stoptime'])

一般に自転車シェアリングでは、平日 or 休日、時間帯、ステーションの位置、天気、周辺のイベントなどの影響で需要が大きく変化する。
したがって、指数分布も定常的なものではなく、非定常的なものになる。

今回の検証は非常に簡単なもので、平日 or 休日、時間帯、ステーションの位置のみを考慮する。
具体的には、平日 or 休日、1時間ごとのタイムフレーム、ステーションIDをそれぞれランダムサンプリングして貸出間隔の分布を調べる。

weekday, weekend = 0, 1
day = [weekday, weekend]
time = np.arange(0, 24)
stationID = np.unique(df['start station id'])

# 組み合わせ
N = len(day) * len(time) * len(stationID)
comb = np.array(np.meshgrid(day, time, stationID)).reshape(3, N).T

# ランダムにシャッフル
np.random.shuffle(comb)

祝日の扱い方

ちなみに、2019年9月2日が祝日(Labor Day)となっている。
この日のデータを見てみると、Trip数は少ないが需要のピークが午前と午後のラッシュ時にあるように平日に近い需要分布になっている。
この分類については今回の検証の範囲外であるため、とりあえず2019年9月2日は平日とみなして検証する。

検証

サンプリング

指数分布に従うことを確認したいのだが、データ数が多い方が確認しやすいため利用回数が100以上あったものに限定した。

# 図のタイトル
subplot_titles = []

# サンプルの抽出
sample = []
N_sample = 10

for d, t, s in comb:
    # 平日と休日で分ける
    if d == weekday:
        condition = (df['starttime'].dt.weekday < 5)
    else:
        condition = (df['starttime'].dt.weekday >= 5)

    # 条件を抽出
    condition = condition & (df['starttime'].dt.hour == t) & (df['start station id'] == s)
    if np.count_nonzero(condition) < 100:
        continue
    else:
        sample.append(condition)
        
        # タイトルの設定
        if d == weekday:
            subplot_titles.append('weekday, ' + str(t) + ', ' + str(s))
        else:
            subplot_titles.append('weekend, ' + str(t) + ', ' + str(s))
        if len(sample) == N_sample: break

サンプルの貸出間隔の累積分布を求め、指数分布で最尤推定する。
そして、その結果をプロットする。

# 図の準備
fig = make_subplots(rows=N_sample//2, cols=2, subplot_titles=subplot_titles)

for i, c in enumerate(sample):
    # 条件に合う行を取得
    sample_df = df[c]

    # 間隔を秒単位に変換
    startTime = np.sort(sample_df['starttime'])
    startInterval = (startTime[1:] - startTime[:-1]).astype(np.float64) / 1e9

    # タイムフレームの1時間を超えたものを除く
    startInterval = startInterval[startInterval < 60 * 60]

    # 累積分布を求める(正規化する)
    hist_startInterval = np.histogram(startInterval, bins=400, range=(0, 60*60))
    x = hist_startInterval[1]
    cdf_startInterval = hist_startInterval[0].cumsum() / hist_startInterval[0].sum()

    # 指数分布の最尤推定
    loc, scale = expon.fit(startInterval)
    ex_cdf = expon.cdf(np.linspace(0, 60*60, 400), loc=loc, scale=scale)

    # plotする(累積分布は1つずれるので0から始める)
    trace0 = go.Bar(x=x, y=np.append(0, cdf_startInterval))
    trace1 = go.Scatter(x=x, y=ex_cdf, mode='lines', opacity=0.7)
    fig.add_trace(trace0, row=i//2+1, col=i%2+1)
    fig.add_trace(trace1, row=i//2+1, col=i%2+1)

    # 軸の設定
    fig.update_xaxes(title_text='sec', gridwidth=0.2, showline=True, linewidth=0.3, range=(-10, 1600), linecolor='black', row=i//2+1, col=i%2+1)
    fig.update_yaxes(gridwidth=0.2, showline=True, linewidth=0.3, range=(-0.05, 1), dtick=0.2, linecolor='black', row=i//2+1, col=i%2+1)
    
fig.update_layout(width=500, height=1000, font=dict(family='Arial'), margin=dict(l=12, r=6, t=24, b=12), showlegend=False)
fig.show()

結果

次のような結果が得られる。

10通りをサンプリングしたが、いずれもよく一致しているように見える。
確かに、自転車シェアリングの貸出間隔は指数分布に従うようである。

休日の利用数が少ないため、平日だけのデータになってしまった。

時間・空間的な偏り

上述してきたような検証ではランダム性を担保するためにTrip数が多いものをピックアップした。
しかし、実際はTrip数が少ない場合の方が多いことに注意する必要がある。

下に1日の平均の利用回数の散布図を示す。平均であるため、平日は21で、休日は10で割っている。
横軸が平日、縦軸が休日になっているので、各点は時間帯とステーションで分類されている。(平日と休日の相関を示したいわけではない。)

青い点は7, 8, 17, 18時台の「ラッシュ時」であり、ダイアモンドのマーカーは利用回数がTop10%に入るステーションである。
特に平日(横軸)はラッシュ時の利用が多く、また、平日休日に関わらず需要が一部のステーションに偏っていることが分かる。

平均を赤い線で示しているが、度数分布に示したように平均以下のとても小さい領域に多くの点が収まっていることに注意されたい。

この図から言いたいのは、大半の場合が利用回数が少ないため、指数分布のランダム性が確認できるのはごく限られた一部の場合であるということである。

データ数を増やせば確認できるようになるかもしれない。
例えば、今回は2019年9月のデータを使ったが、これを1年に拡張することもできる。
しかし、季節によって利用パターンが変わったり、ステーションが追加されている場合もあることに気をつける必要もある。

また、今回は無作為に1時間ごとのタイムフレームに分割したが、夜間など利用回数の少ない時間をまとめることも有効であるように感じる。

利用回数が少ない場合はどうなるか(おまけ)

上では利用回数が100回以上の場合を取り上げた。
しかし、中央値は19回であるため、大半が指数分布に従うのか確認できない

なぜかというと、利用回数が少ないと間隔が1時間以内のタイムフレームに入らないものが多くなるからである。

下図に10回以上20回以下の利用があった場合の分布を示す。
ステップが非常に大きく、指数分布にしたがっているのか分からないものがある。

返却間隔

返却間隔も指数分布に従うことができる。
上のコードのstarttimestoptimeに、start station idに、end station idにするだけでできる。

ステーションの台数の推移 ― マルコフ過程

これまで貸出間隔と返却間隔が指数分布に従うことを見てきた。
ところでステーションの台数の推移はどのような過程になるだろうか?

結論を言うとマルコフ過程に従う。
貸出や返却が無相関だとして大きく単純化すると、マルコフ過程の有名な応用である「出生死亡過程」に従うということもできる。

貸出や返却が無相関であるというのは、返却は必ず貸出の後に起こるが、それぞれが全く関係なくランダムに発生するということである。

このモデルを使って自転車シェアリングをモデリングした研究例がある。

A Bike-sharing Optimization Framework Combining Dynamic Rebalancing and User Incentives

スポンサーリンク
H-MEMO

コメント