문제

학교에서 집까지 가는 방법은 버스 A를 타고 10분 동안 이동한 뒤 내려서 5분 동안 걸어가는 방법과 버스 B를 타고 8분 동안 이동한 뒤 내려서 1분 동안 걸어가는 방법이 있다. 버스 A의 배차간격은 3분이고 버스 B의 배차간격은 15분이라고 알려져 있다. (위의 상수는 임의의 숫자로 변경될 수 있음) SeoulBus나 버스정류장의 전광판 등 버스의 도착 예정 시간을 알 수 있는 방법이 없다고 할 때, 어떤 기준으로 버스를 타야 도착 시간의 기대값을 가장 빠르게 할 수 있을까?

풀이

이 문제는 내가 버스를 기다리기 위해 정류장에 도착한 시각부터 각 버스가 오기까지 기다리는 기대 시간이 얼마인지 계산할 수 있으면 되는 문제이다. 이는 가정에 따라 생각보다 복잡해질 수도 있다. 그 이유는 length time bias 때문이다.

Length Time Bias

만약 버스가 정확하게 15분 혹은 3분의 배차 간격을 두고 운행이 된다면, 평균적으로 내가 버스를 기다리는 시간은 각각 7.5분과 1.5분이 된다. 하지만 정확하게 시간을 지키지 않고 평균적으로 그 정도의 시간만 지키면서 조금씩 달라질 수 있다면 결과가 달라지게 된다. 왜냐하면 버스 배차 간격이 긴 경우에 내가 정류장에 도달할 확률이 그렇지 않은 경우보다 높기 때문이다.

앞서 지나간 버스와 이를 뒤따르는 버스 사이의 시간 간격을 배차 간격이라 하자. 배차 간격이 주어진 경우에 승객이 정류장에 도달하는 것은 균등 분포를 따른다. 그러므로 기대 대기 시간은 배차 간격의 절반이 된다. 하지만 승객이 정류장에 도착하는 것이 균등하다면 더 큰 배차 간격에 떨어질 확률이 높다.

In [1]:
import matplotlib.image as mpimg
img = mpimg.imread('waitt.gif')
imshow(img)  # Image from http://mahalanobis.twoday.net/stories/3486587/
Out[1]:
<matplotlib.image.AxesImage at 0x10ae565d0>

위의 이미지에서 $X_{i}$는 배차 간격, $\bar{W}$는 평균 대기 시간이다. 이는 위의 삼각형으로 이루어진 함수의 평균 값이 된다.

$\bar{W} = \frac{1}{t} \sum_{i=1}^{n} \frac{1}{2} X_{i}^{2}$

평균 대기 시간은 전체 시간인 $t$를 버스의 총 개수인 $n$으로 나눈 값이므로,

$\bar{W} = \frac{1}{2} \frac{\bar{X^{2}}}{\bar{X}}$

이다.

Poisson/Exponential인 경우

만약 버스가 도달하는 사건이 푸아송 분포를 따른다면, 그 배차 간격은 지수 분포를 따를 것이며 다음의 수식을 만족한다.

$Var(X) = \bar{X}^{2}$

그리고 분산은 제곱의 평균 - 평균의 제곱이므로,

$\bar{X^{2}} = 2 \bar{X}^{2}$

이를 활용하면,

$\bar{W} = \bar{X}$

가 된다. 사실 이런 방법을 쓰지 않더라도 지수 분포는 memory less 성질이라 하여 승객의 도착 시간과 무관하게 평균 대기 시간은 항상 같으므로 $\bar{W} = \bar{X}$ 임을 알 수 있다.

정규 분포인 경우

$Var(X) = \sigma^{2}$

위와 같은 방법을 사용하면,

$\bar{X^{2}} = \sigma^{2} + \bar{X}^{2}$

이를 똑같이 대입하면,

$\bar{W} = \frac{1}{2} \frac{\sigma^{2} + \bar{X}^{2}}{\bar{X}}$

이다.

만약 배차 간격이 정확하게 지켜진다면 분산은 0이 될 것이고 그렇다면 $\bar{W} = \frac{1}{2} \bar{X}$가 될 것이다. 이는 앞서 언급한 결과와 일치한다.

엄밀히 말하자면 정규 분포는 음의 값을 가질 수 있으므로 배차 간격을 모델링하기에 적합하지 않다. 하지만 현실에서는 버스 사이의 배차 간격을 운전수들이 조절하려고 노력하기 때문에 (표준 편차가 현실적이라면) 배차 간격이 지수 분포를 따르지 않는 경우를 모델링 할 수 있을 것이다.

전략

버스 배차 간격이 지수 분포를 따른다면 자명하게 먼저 오는 버스를 타는 것이 가장 좋은 전략이다. 버스 B가 먼저 도착한다면 당연히 버스 B를 타는 것이 좋고, 버스 A가 먼저 도착하더라도 버스 B가 도착하기까지의 기대 대기 시간이 15분이기 때문이다.

만약 버스 배차 간격이 정규 분포를 따른다면 어떨까?

버스 B가 먼저 도착하면 당연히 버스 B를 타면 된다. 버스 A가 먼저 도착한다면 다음 버스 B가 도달하는 시간을 추론해 볼 수 있다. 내가 버스 A를 $\alpha$분 만큼 기다렸다면 버스 B가 도착하기까지 평균적으로 $\frac{\sigma^{2} + 225}{30} - \alpha$분 남았다고 생각할 수 있다. 하지만 이미 $\alpha$만큼의 시간이 지났으므로 그대로 계산할 수는 없다.

(분산이 충분히 작아서 무시할 수 있다고 가정하면) 이미 지나간 $\alpha$만큼의 시간은 일어날 확률이 0이 되므로 평균 대기 시간은 $\frac{15 + \alpha}{2}$가 되고 앞으로 남은 기대 대기 시간은 $\frac{15 + \alpha}{2} - \alpha = \frac{15 - \alpha}{2}$이다.

지금 당장 버스 A를 타면 15분 후 집에 도착하게 되고, 버스 B를 기다렸다가 타면 평균적으로 $7.5 - \frac{1}{2}\alpha + 9 = 16.5 - \frac{1}{2}\alpha$분 후에 집에 도착하게 된다. 결국 내가 3분 이상 버스 A를 기다렸다면, 버스 B가 오기를 기다렸다가 버스 B를 타고 내려가는 것이 평균적으로 집에 더 빠르게 도착하는 방법이 된다.

실험

위에서는 수식을 사용해서 문제를 풀어봤는데, 실제로 시뮬레이션을 해보면 어떨까?

In [2]:
import numpy as np
import pymc as mc
import statsmodels.api as sm
Couldn't import dot_parser, loading of dot files will not be possible.

우선 length-time bias의 효과를 살펴보자.

In [3]:
a_interval = 3.0
tau1 = 2.0

unbiased_a = mc.Normal('unbiased_a', mu=a_interval, tau=tau1)

@mc.stochastic(name='unbiased_abs_a', trace=True)
def unbiased_abs_a(value=5.0, mu=a_interval, tau=tau1):
    if value < 0:
        return np.NINF
    return mc.normal_like(value, mu, tau)

@mc.stochastic(name='biased_a', trace=True)
def biased_a(value=5.0, mu=a_interval, tau=tau1):
    if value < 0:
        return np.NINF
    return np.sum(np.log(value) + mc.normal_like(value, mu, tau))
In [4]:
model = mc.MCMC([biased_a, unbiased_a, unbiased_abs_a])
model.sample(iter=100000, burn=1000)
 [-----------------100%-----------------] 100000 of 100000 complete in 12.2 sec
In [5]:
x = model.trace('biased_a')[:]
y = model.trace('unbiased_a')[:]
z = model.trace('unbiased_abs_a')[:]
In [6]:
figure(figsize=(10, 5))
hist(x, alpha=0.5, label="Biased A", bins=50, normed=True)
hist(z, alpha=0.5, label="Unbiased A", bins=50, normed=True)
legend()
title("Comparison of Bus Arrival Time")
xlabel("Time (min)")
ylabel("Frequency")
savefig("dist_comp_bus_a.png")

확실히 bias가 있는 샘플이 더 긴 체감 시간을 느끼게 된다.

조금 더 확실하게 확인할 수 있도록 CDF를 비교해보자.

In [7]:
figure(figsize=(10, 5))

ecdf_x = sm.distributions.ECDF(x)
ecdf_z = sm.distributions.ECDF(z)

x_x = np.linspace(min(x), max(x))
y_x = ecdf_x(x_x)

x_z = np.linspace(min(z), max(z))
y_z = ecdf_z(x_z)

plot(x_x, y_x, label="Biased A")
plot(x_z, y_z, label="Unbiased A")

legend()
xlabel("Time")
ylabel("CDF")
title("Comparison of Bus Arrival Time (CDF)")
savefig("dist_comp_bus_a_cdf.png")
In [8]:
b_interval = 15.0
tau = .5

unbiased_b = mc.Normal('unbiased_b', mu=b_interval, tau=tau)

@mc.stochastic(name='unbiased_abs_b', trace=True)
def unbiased_abs_b(value=5.0, mu=b_interval, tau=tau):
    if value < 0:
        return np.NINF
    return mc.normal_like(value, mu, tau)

@mc.stochastic(name='biased_b', trace=True)
def biased_b(value=5.0, mu=b_interval, tau=tau):
    if value < 0:
        return np.NINF
    return np.sum(np.log(value) + mc.normal_like(value, mu, tau))
In [9]:
model2 = mc.MCMC([biased_b, unbiased_b, unbiased_abs_b])
model2.sample(iter=100000, burn=1000)
 [-----------------100%-----------------] 100000 of 100000 complete in 12.8 sec
In [10]:
x2 = model2.trace('biased_b')[:]
y2 = model2.trace('unbiased_b')[:]
z2 = model2.trace('unbiased_abs_b')[:]
In [11]:
figure(figsize=(10, 5))
hist(x2, alpha=0.5, label="Biased B", bins=50, normed=True)
hist(z2, alpha=0.5, label="Unbiased B", bins=50, normed=True)
legend()
xlabel("Time (min)")
ylabel("Frequency")
title("Comparison of Bus Arrival Time")
savefig("dist_comp_bus_b.png")
In [12]:
figure(figsize=(10, 5))

ecdf_x2 = sm.distributions.ECDF(x2)
ecdf_z2 = sm.distributions.ECDF(z2)

x_x2 = np.linspace(min(x2), max(x2))
y_x2 = ecdf_x2(x_x2)

x_z2 = np.linspace(min(z2), max(z2))
y_z2 = ecdf_z2(x_z2)

plot(x_x2, y_x2, label="Biased B")
plot(x_z2, y_z2, label="Unbiased B")

legend()
xlabel("Time")
ylabel("CDF")
title("Comparison of Bus Arrival Time (CDF)")
savefig("dist_comp_bus_b_cdf.png")

$\tau$가 작아서 잘 보기 힘들지만 역시 우편향된 분포를 볼 수 있다.

In [13]:
@mc.deterministic
def waittime_a(sample_a=biased_a):
    return sample_a * np.random.random()
In [14]:
@mc.deterministic
def waittime_b(sample_b=biased_b):
    return sample_b * np.random.random()
In [15]:
m = mc.MCMC([biased_a, biased_b, waittime_a, waittime_b])
m.sample(iter=100000, burn=1000)
 [-----------------100%-----------------] 100000 of 100000 complete in 13.1 sec
In [16]:
wt_a = m.trace('waittime_a')[:]
wt_b = m.trace('waittime_b')[:]
In [17]:
figure(figsize=(10, 5))
hist(wt_a, alpha=0.5, bins=50, range=(0, 20), label="Bus A", normed=True)
hist(wt_b, alpha=0.5, bins=50, range=(0, 20), label="Bus B", normed=True)
title("Comparison of Bus Wait Time")
xlabel("Time (min)")
ylabel("Frequency")
legend()
savefig("dist_comp_wait_time.png")

대기 시간의 분포이다.

In [18]:
def hop_on_fast_arrival(a_wait, b_wait):    
    a_arrival = a_wait + 15
    b_arrival = b_wait + 9
    
    result = np.empty(len(a_wait))
    
    result[(b_wait <= a_wait) & (b_arrival <= a_arrival)] = 1
    result[(b_wait <= a_wait) & (b_arrival > a_arrival)] = 0
    
    result[(b_wait > a_wait) & (b_arrival >= a_arrival)] = 1
    result[(b_wait > a_wait) & (b_arrival < a_arrival)] = 0
    return result
In [19]:
result = hop_on_fast_arrival(wt_a, wt_b)
In [20]:
np.count_nonzero(result)
Out[20]:
58826
In [21]:
np.count_nonzero(result) / float(len(result))
Out[21]:
0.5942020202020202

먼저 오는 버스를 무조건 타는 전략은 약 60%의 확률로 빨리 집에 갈 수 있다.

In [22]:
def wait_for_b(a_wait, b_wait, k=3.0):    
    a_arrival = a_wait + 15
    b_arrival = b_wait + 9
    
    result = np.empty(len(a_wait))
    
    result[(b_wait <= a_wait) & (b_arrival <= a_arrival)] = 1
    result[(b_wait <= a_wait) & (b_arrival > a_arrival)] = 0
    
    # a가 먼저 오고 k 분이 넘게 지난 상황. b를 기다렸다 탄다.
    result[(b_wait > a_wait) & (a_wait >= k) & (b_arrival <= a_arrival)] = 1
    result[(b_wait > a_wait) & (a_wait >= k) & (b_arrival > a_arrival)] = 0

    # a가 먼저 오고 아직 k 분 못되게 기다린 상황. 바로 a를 탄다.
    result[(b_wait > a_wait) & (a_wait < k) & (b_arrival <= a_arrival)] = 0
    result[(b_wait > a_wait) & (a_wait < k) & (b_arrival > a_arrival)] = 1
    return result
In [23]:
grid = np.linspace(0, 10, 50)
v = []
for k in grid:
    d = wait_for_b(wt_a, wt_b, k)
    v.append(np.count_nonzero(d)/float(len(d)))
In [24]:
figure(figsize=(10, 5))
plot(grid, v)
xlabel("Wait Time (min)")
ylabel("Winning Rate")
title("Winning Rate for Various Wait Time")
savefig("wait_time_success_rate.png")
In [25]:
print grid[np.argmax(v)], np.max(v)
2.65306122449 0.599686868687
In [26]:
np.max(v)/float(len(result))
Out[26]:
6.0574431180491793e-06

위의 예제에서는 약 3.26분(3분 15초) 기다리는 동안 버스 A가 오지 않았다면 버스 A가 곧 오더라도 더 기다렸다가 버스 B를 타는 편이 이득이다. 하지만 그 이득은 미미하다.

In [27]:
def hop_on_fast_arrival_time(a_wait, b_wait):    
    a_arrival = a_wait + 15
    b_arrival = b_wait + 9
    
    result = np.empty(len(a_wait))
    
    result[b_wait <= a_wait] = b_arrival[b_wait <= a_wait]
    result[b_wait > a_wait] = a_arrival[b_wait > a_wait]
    return result
In [28]:
strategy1 = hop_on_fast_arrival_time(wt_a, wt_b)
In [29]:
figure(figsize=(10, 5))
hist(strategy1, bins=50)
xlabel("Time (min)")
ylabel("Frequency")
title("Home Arrival Time - FCFS")
savefig("home_arrival_time_fcfs.png")
In [30]:
print "Mean arrival time: {}".format(strategy1.mean())
print "STD of arrival time: {}".format(strategy1.std())
print "Median arrival time: {}".format(np.median(strategy1))
Mean arrival time: 15.8524538159
STD of arrival time: 2.17607477284
Median arrival time: 16.2567454086
In [31]:
def wait_for_b_time(a_wait, b_wait, k=3.0):    
    a_arrival = a_wait + 15
    b_arrival = b_wait + 9
    
    result = np.empty(len(a_wait))
    
    result[b_wait <= a_wait] = b_arrival[b_wait <= a_wait]    
    # a가 먼저 오고 k 분이 넘게 지난 상황. b를 기다렸다 탄다.
    result[(b_wait > a_wait) & (a_wait >= k)] = b_arrival[(b_wait > a_wait) & (a_wait >= k)]
    # a가 먼저 오고 아직 k 분 못되게 기다린 상황. 바로 a를 탄다.
    result[(b_wait > a_wait) & (a_wait < k)] = a_arrival[(b_wait > a_wait) & (a_wait < k)]
    return result
In [32]:
strategy2 = wait_for_b_time(wt_a, wt_b, k=3.15)
In [33]:
figure(figsize=(10, 5))
hist(strategy2, bins=50)
xlabel("Time (min)")
ylabel("Frequency")
title("Home Arrival Time - WFB")
savefig("home_arrival_time_wfb.png")
In [34]:
print "Mean arrival time: {}".format(strategy2.mean())
print "STD of arrival time: {}".format(strategy2.std())
print "Median arrival time: {}".format(np.median(strategy2))
Mean arrival time: 15.8404631999
STD of arrival time: 2.31102149876
Median arrival time: 16.1980195395
In [35]:
def wait_twice(a_wait, b_wait, k=3.0, l=3.0, seed=42):  
    if k < 0:
        k = 0
    if l < 0:
        l = 0
    # 빠른 버스가 먼저 오면 탄다.
    # 느린 버스가 먼저 오면 지금 기다린 시간 k에 따라 다른 전략을 취한다
    # 그 후에도 한 번 더 느린 버스가 오면 기다린 시간 l에 따라 다른 전략을 취한다
    a_arrival = a_wait + 15
    b_arrival = b_wait + 9
    
    result = np.empty(len(a_wait))
    
    # b가 먼저 오면 b를 탄다.
    result[b_wait <= a_wait] = b_arrival[b_wait <= a_wait]
    # a가 먼저 오고 아직 k 분 못되게 기다린 상황. 바로 a를 탄다.
    result[(b_wait > a_wait) & (a_wait < k)] = a_arrival[(b_wait > a_wait) & (a_wait < k)]
    
    # a가 먼저 오고 k 분이 넘게 지난 상황. b를 일단 기다려 본다.
    cond = (b_wait > a_wait) & (a_wait >= k)
    size = np.count_nonzero(cond)
    # 두 번째 버스는 unbiased distribution에서 뽑아야 함
    next_buses_a = abs(np.random.normal(a_interval, np.sqrt(1/tau), size=size))
    # 두 번째 버스 a가 오기 전 총 대기 시간
    wait_total = next_buses_a + a_wait[cond]
    
    tmp = np.empty(size)
    b_slices = b_arrival[cond]
        
    # 다음 버스가 b이면 b를 탄다.
    tmp[wait_total >= b_wait[cond]] = b_slices[wait_total >= b_wait[cond]]
    # 다음 버스가 또 a이면 고민을 한다.
    # 추가 대기 시간이 l 이상이면 b를 기다렸다 탄다.
    tmp[(wait_total < b_wait[cond]) & (next_buses_a >= l)] = b_slices[(wait_total < b_wait[cond]) & (next_buses_a >= l)]
    # 추가 대기 시간이 l 미만이면 그냥 a를 탄다.
    tmp[(wait_total < b_wait[cond]) & (next_buses_a < l)] = wait_total[(wait_total < b_wait[cond]) & (next_buses_a < l)] + 9

    result[cond] = tmp
    return result
In [36]:
strategy3 = wait_twice(wt_a, wt_b, k=0.0, l=6.0)
In [37]:
print "Mean arrival time: {}".format(strategy3.mean())
print "STD of arrival time: {}".format(strategy3.std())
print "Median arrival time: {}".format(np.median(strategy3))
Mean arrival time: 12.8341853053
STD of arrival time: 1.96126497269
Median arrival time: 12.8047337845
In [38]:
figure(figsize=(10, 5))
hist(strategy3, bins=50)
xlabel("Time (min)")
ylabel("Frequency")
title("Home Arrival Time - Optimal")
savefig("home_arrival_time_optimal.png")
In [39]:
grid = np.linspace(0, 20, 50)
v = []
for i in grid:
    for j in grid:
        v.append(wait_twice(wt_a, wt_b, k=i, l=j).mean())
In [40]:
figure(figsize=(10, 5))
plot(v)
xlabel("Patched Time")
ylabel("Arrival Time")
title("Comparison of Different Wait Time (2-D)")
savefig("home_arrival_time_comp_2d.png")
In [41]:
figure(figsize=(10, 5))
for i in range(0, 3):
    plot(grid, v[i*len(grid):(i+1)*len(grid)], label="{} min".format(i))
legend()
xlabel("Time (min)")
ylabel("Arrival Time")
title("Comparison of Different Wait Time (1-D)")
savefig("home_arrival_time_comp_1d.png")
In [42]:
wait_twice(wt_a, wt_b, k=0, l=10).mean()
Out[42]:
12.812273737434559

전략3

  • 버스 B(빠르지만 가끔 오는 버스)가 오면 무조건 먼저 탄다.
  • 버스 A(느리지만 자주 오는 버스)가 오면 현재 대시 시간에 따라 버스를 타거나 다음 버스를 기다린다.
    • 위의 경우 다음에 다시 버스 A가 오면 현재 추가 대기 시간에 따라 버스를 타거나 다음 버스를 기다린다.

위의 결과를 보면 첫 버스 A는 무조건 보내는 편이 평균 대기 시간이 더 짧아진다.