[실전 시계열 분석] - chap04 시계열 데이터의 시뮬레이션

chap04의 목적

- 시계열 데이터 시뮬레이션 vs 다른 동류 데이터 시뮬레이션 비교

- 실제 코드 기반 시뮬레이션 예 살펴보기

- 시계열 시뮬레이션 동향 살펴보기

---

4.1 시계열 시뮬레이션의 특별한 점

동일한 시계열에서는 두 데이터가 서로 다른 시간에 일어나므로 정확하게 비교하는 것은 어렵다. 따라서 특정 시간에 발생 가능한 일을 예측하려면 시뮬레이션을 실행해봐야 한다.

4.1.1 시뮬레이션과 예측

시뮬레이션과 예측은 유사하나, 차이점이 있다.

- 정성적 관측을 예측보다 시뮬레이션에 통합하는 것이 더 쉬울 수 있다.

- 시뮬레이션은 확장 가능하므로, 예측은 시뮬레이션보다 더 신중하게 처리되어야 한다.

- 시뮬레이션은 가상이므로 예측보다 위험 부담이 적다. 창의적이고 탐구적인 자세로 설계할 수 있다.

---

4.2 코드로 보는 시뮬레이션

예제1 : 누군가의 이메일 열람 행동이 기부로 이어지는지 여부에 대한 상관관계 가설 검정

예제2 : 택시 기사의 교대 시간과 하루 동안 탑승객 빈도에 대한 택시 무리의 집단 행동 살펴보기

예제3 : 점진적으로 개별 자기요소(magnetic element)의 위치를 맞춰나가는 자성물질의 물리적 과정을 시뮬레이션

4.2.1 스스로 직접 만들어보기

위와 같이 모든 회원에게 특정 가입 연도를 무작위로 부여하고, 회원의 상태 정보는 부여된 가입연도에 따라 결정된다.

이번에는 주별로 회원의 이메일 열람 시점을 나타내는 테이블을 만든다. 한 주에 이메일 세 통을 보내는 기관의 행동을 정의하고, 이메일에 관한 회원들의 행동 패턴은 다음과 같이 정의한다.

- 이메일 열람 기록 없음

- 일정한 수준의 이메일 열람 및 참여율

- 참여 수준의 증가 또는 감소

NUM_EMAILS_SENT_WEEKLY = 3

## 서로다른 패턴을 위한 몇 가지 함수를 정의 합니다
def never_opens(period_rng): 
  return []
					
def constant_open_rate(period_rng): 
  n, p = NUM_EMAILS_SENT_WEEKLY, np.random.uniform(0, 1) 
  num_opened = np.random.binomial(n, p, len(period_rng)) 
  return num_opened
					
def increasing_open_rate(period_rng): 
  return open_rate_with_factor_change(period_rng,				
                                      np.random.uniform(1.01, 
                                      1.30))
					
def decreasing_open_rate(period_rng): 
  return open_rate_with_factor_change(period_rng,
                                      np.random.uniform(0.5, 
                                      0.99))
				
def open_rate_with_factor_change(period_rng, fac): 
    if len(period_rng) < 1 :	
       return []
    times = np.random.randint(0, len(period_rng),
                               int(0.1 * len(period_rng))) 
    num_opened = np.zeros(len(period_rng))
    for prd in range(0, len(period_rng), 2): 
        try:
            n, p = NUM_EMAILS_SENT_WEEKLY, np.random.uniform(0, 
                                                             1)
            num_opened[prd:(prd + 2)] = np.random.binomial(n, p, 
                                                           2)
            p = max(min(1, p * fac), 0)
        except:	
            num_opened[prd] = np.random.binomial(n, p, 1)
    for t in range(len(times)):
        num_opened[times[t]] = 0
    return num_opened

4가지 행동에 따른 함수 정의

- never_opens() : 이메일을 한 번도 열람하지 않은 회원

- constant_open_rate() : 매주 같은 양의 이메일을 열람한 회원

- decreasing_open_rate() : 매주 열람한 이메일의 양이 줄어드는 회원

- increasing_open_rate() : 매주 열람한 이메일의 양이 늘어나는 회원

위 함수를 정의하고, 이제 기부 행동을 모델링하는 시스템을 생각해봐야 한다.

## 기부 행동
def produce_donations(period_rng, member_behavior, num_emails, 
                      use_id, member_join_year):
    donation_amounts = np.array([0, 25, 50, 75, 100, 250, 500, 
                                 1000, 1500, 2000])
    member_has = np.random.choice(donation_amounts)
    email_fraction = num_emails / (NUM_EMAILS_SENT_WEEKLY * len(period_rng)) 
    member_gives = member_has * email_fraction
    member_gives_idx = np.where(member_gives >= donation_amounts)[0][-1]
    member_gives_idx = max(min(member_gives_idx,
                               len(donation_amounts) - 2),
                           1)
    num_times_gave = np.random.poisson(2) * (2018 - member_join_year)
    times = np.random.randint(0, len(period_rng), num_times_gave)
    dons = pd.DataFrame({'member'   : [],
                         'amount'   : [],
                         'timestamp': []})
					
    for n in range(num_times_gave):
        donation = donation_amounts[member_gives_idx + np.random.binomial(1, .3)]
        ts = str(period_rng[times[n]].start_time + random_weekly_time_delta())
        dons = dons.append(pd.DataFrame(
                  {'member' :[use_id],
                   'amount' :[donation],
                   'timestamp': [ts]}))

        if dons.shape[0] > 0:
            dons = dons[dons.amount != 0]
            ## we don't report zero donation events as this would not
            ## be recorded in a real world database 
            ## 기부액이 0인 경우에는 보고하지 않습니다.
            ## 실세계에서 다룰 때 이러한 정보는 데이터베이스에 반영되지 않습니다.
    
            return dons

- 회원자격의 기간에 따라서 기부의 전체 횟수를 정함

- 회원별 재정 상태를 만들어, 한 사람의 안정적인 재산의 양이 기부의 양에 밀접한 관련이 있다는 행동 가설에 기반

def random_weekly_time_delta():
    days_of_week = [d for d in range(7)]
    hours_of_day = [h for h in range(11, 23)] 
    minute_of_hour = [m for m in range(60)] 
    second_of_minute = [s for s in range(60)]
    return pd.Timedelta(str(np.random.choice(days_of_week)) + " days" ) + pd.Timedelta(str(np.random.choice(hours_of_day)) + " hours" ) + pd.Timedelta(str(np.random.choice(minute_of_hour)) + " minutes") + pd.Timedelt

회원의 행동은 특정 타임스탬프와 연관되어 있어, 각 회원이 기부를 한 주와 각 주 내에 기부를 한 시점을 고르는 코드 작성. 특정 주 내의 시간을 무작위로 고르기 위한 유틸리티 함수

behaviors        = [never_opens, 
                    constant_open_rate, 
                    increasing_open_rate, 
                    decreasing_open_rate] 
member_behaviors = np.random.choice(behaviors, 1000,
                                    [0.2, 0.5, 0.1, 0.2]) 
				
rng = pd.period_range('2015-02-14', '2018-06-01', freq = 'W')
emails = pd.DataFrame({'member'      : [],
                       'week'        : [],
                       'emailsOpened': []})	
donations = pd.DataFrame({'member'   : [],
                          'amount'   : [],
                          'tunestamp': []})	
					
for idx in range(yearJoined.shape[0]):
    ## randomly generate the date when a member would have joined
    ## 회원이 가입한 시기를 무작위로 생성합니다
    join_date = pd.Timestamp(yearJoined.iloc[idx].yearJoined) + pd.Timedelta(str(np.random.randint(0, 365)) + ' days')
    join_date = min(join_date, pd.Timestamp('2018-06-01')).to_period(freq='W')

    ## member should not have action timestamps before joining
    ## 가입 전 어떤 행동에 대한 타임스탬프도 가져서는 안됩니다
    member_rng = rng[rng > join_date] 

    if len(member_rng) < 1:
        continue

    info = member_behaviors[idx](member_rng) 
    if len(info) == len(member_rng):
        emails = emails.append(pd.DataFrame(
           {'member': [idx] * len(info),
            'week': [str(r.start_time) for r in member_rng],
            'emailsOpened': info})) 
        donations = donations.append(
           produce_donations(member_rng, member_behaviors[idx], 
                                sum(info), idx, join_date.year))

2015년에서 2018년으로 시간이 흐르멩 따라 기부 및 이메일 열람 횟수가 증가한느 것처럼 보인다. 하지만 누적 회원수가 증가했기 때문에 열람된 이메일 수가 증가했다는 사실은 정확하지 않다. 회원의 자격이 무기한으로 유지된다는 가정이 내재되어 있어, 해당 방식을 수정할 필요가 있음.

 

위 시계열을 만드는 데 주의해야 할 점

- 시계열 개수 결정

- 시간에 따라 모델링해야 하는 추세 결정 : 안정, 증가, 감소하는 이메일 열람률 / 안정적인 행동 패턴 만들기

- 회원 가입 이전에는 이메일 열람, 기부하지 않도록 설계

- 데이터가 미래로 누수되지 않도록 설계

---

4.2.2 스스로 실행하는 시뮬레이션 세계 구축

서로 다른 시간에 교대근무가 예정된 택시 무리의 전체 행동 방식 시뮬레이션

import numpy as np 	
		
def taxi_id_number(num_taxis): 
   arr = np.arange(num_taxis)
   np.random.shuffle(arr)
   for i in range(num_taxis):
       yield arr[i]

ids = taxi_id_number(10) 
print(next(ids))
print(next(ids))
print(next(ids)) 

먼저 제네레이터에 대해 이해해야 하는데, 각 개체가 각자의 상태를 독립적으로 보관하는 1회용 객체를 말한다. 출력은 7->2->5로 독립적이면서 각자 자기만의 변수를 신경 쓰는 여러 객체를 만들고 시을 때 유용하다.

def shift_info():
   start_times_and_freqs = [(0, 8), (8, 30), (16, 15)]
   indices               = np.arange(len(start_times_and_freqs)) 
   while True:
       idx = np.random.choice(indices,p=[0.25,0.5,0.25]) 
       start = start_times_and_freqs[idx]
       yield (start[0], start[0] + 7.5, start[1]) 

하루 중 서로 다른 세 개의 교대 시간을 표현하고, 각 시간대에 택시가 할당될 서로 다른 확률을 부여한다. 하루의 서로 다른 교대 시간은 서로 다른 평균 운행 횟수를 가진다.

def taxi_process(taxi_id_generator, shift_info_generator):
   taxi_id = next(taxi_id_generator)
   shift_start, shift_end, shift_mean_trips = next(shift_info_generator)
   actual_trips = round(np.random.normal(loc   = shift_mean_trips,
                                         scale = 2))
   average_trip_time = 6.5 / shift_mean_trips * 60
   # convert mean trip time to minutes
   # 평균 운행 시간을 분 단위로 변환합니다
   between_events_time = 1.0 / (shift_mean_trips - 1) * 60
   # this is an efficient city where cabs are seldom unused
   # 이 도시는 매우 효율적이어서 모든 택시가 거의 항상 사용됩니다
   time = shift_start
   yield TimePoint(taxi_id, 'start shift', time)
   deltaT = np.random.poisson(between_events_time) / 60
   time += deltaT
   for i in range(actual_trips):
       yield TimePoint(taxi_id, 'pick up ', time)
       deltaT = np.random.poisson(average_trip_time) / 60
       time += deltaT
       yield TimePoint(taxi_id, 'drop off ', time)
       deltaT = np.random.poisson(between_events_time) / 60
       time += deltaT
   deltaT = np.random.poisson(between_events_time) / 60
   time += deltaT
   yield TimePoint(taxi_id, 'end shift ', time) 

두 개의 제네레이터를 통해 각 택시 ID번호, 교대 시작 시간, 해당 시간에 대한 평균 운행 횟수를 결정한다. 이에 따라 각 택시는 자신만의 시간표와 특정 운햇 횟수에 따라 홀로 여행을 떠나고, 외부의 next()호출로 제네레이터에 접근하는 클라이언트에게 그 결과를 내어준다.

## python 							
from dataclasses import dataclass
 							
@dataclass
class TimePoint:
   taxi_id: int
   name: str		
   time: float
				
   def __lt__(self, other):			
       return self.time < other.time 	

위는 택시 제네레이터가 생산하는 TimePoint 객체의 형태이다. 

import queue
			
class Simulator:
   def __init__(self, num_taxis):
       self._time_points = queue.PriorityQueue()
       taxi_id_generator = taxi_id_number(num_taxis)
       shift_info_generator = shift_info()
       self._taxis = [taxi_process(taxi_id_generator,
                                   shift_info_generator) for
                                            i in range(num_taxis)]
       self._prepare_run()

   def _prepare_run(self): 
       for t in self._taxis:
           while True: 
               try:
                   e = next(t)
                   self._time_points.put(e)
               except:
                   break
					
   def run(self): 
       sim_time = 0
       while sim_time < 24:
           if self._time_points.empty():
               break
           p = self._time_points.get()
           sim_time = p.time
           print(p) 	

필요한 택시 개수에 따른 택시 제네레이터를 생성하고, 반복적으로 접근하여 각각 반환한 timepoint가 유효하다면 우선순위 큐에 삽입한다. 우선순위 큐에 쌓인 TimePoint들은 시간의 순서대로 내보내질 수 있다.

---

4.2.3 물리적인 시뮬레이션

이징 모델을 통해 MCMC(마르코프 연쇄 몬테카를로) 시뮬레이션 구현

 - 격자 형태 시스템 환경 세팅

- 시작 블록을 무작위로 초기화

- 인접 상태에 비례하여 중앙 정렬 상태 에너지 계산

- 전체 블록의 자화(magnetization)

위는 무작위로 생성된 강자성 물질으 초기 상태이다. 확률적으로 완벽한 5:5 체커판 상태는 얻기 어렵다. 초기 상태에서 어떤 패턴을 발견했더라도 실제로는 패턴이 아니지만 사람의 뇌가 그렇게 인식했을 가능성이 높다.

1000번째 단계의 스냅샷으로, 지배적인 상태가 초기와 비교해 역전됐다.

각 격자 지점이 무작위로 초기화된 경우에서조차, 시스템이 저온에서 자화 상태가 될 수 있는 가능성을 보여주는 결과이다. 

---

4.3 시뮬레이션에 대한 마지막 조언

시뮬레이션으로부터 얻은 양적 측정지표와 결합된 가상적인 예를 통해 데이터에 대한 인지력을 확장해나갈 수 있다.

4.3.1 통계적인 시뮬레이션

확률적인 역동성을 이미 알고 있을 때, 몇 가지 모르는 파라미터를 추정하거나 서로 다른 가정이 파라미터 추정 과정에 주는 영향을 알아보고자 한다면 유용할 수 있다.

4.3.2 딥러닝 시뮬레이션

딥러닝 시뮬레이션의 장점은 비선형적 역동성을 잡아낼 수 있다는 것이다. 한편 시스템 역동성의 근본 원리를 전혀 이해하지 못하게 될 수도 있다는 단점이 있다.