4장 시계열 데이터의 시뮬레이션 (물리)

import numpy as np

### 환경 설정
## 물리적 배치
N           = 5 # 격자의 너비
M           = 5 # 격자의 높이

## 온도 설정
temperature = 0.5
BETA        = 1 / temperature

def initRandState(N, M):		
    block = np.random.choice([-1, 1], size = (N, M))
    return block

def costForCenterState(state, i, j, n, m):
    centerS = state[i, j]
    neighbors=[((i + 1) % n, j), ((i - 1) % n, j),
               (i, (j + 1) % m), (i, (j - 1) % m)]
               
    ## 주기적 경계 조건의 상정을 위한 % n 부분에 주목하시기 바랍니다
    ## 이 부분이 이해가 안된다면, 그냥 무시하셔도 좋습니다
    ## 이는 2D 시스템이 도너츠의 표면과 같다고 말하는 시스템의 물리적인 졔약일 뿐입니다.
    interactionE = [state[x, y] * centerS for (x, y) in neighbors]
    return np.sum(interactionE)

def magnetizationForState(state): 
   return np.sum(state)		

def mcmcAdjust(state):		
    n = state.shape[0] 					
    m = state.shape[1]
    x, y = np.random.randint(0, n), np.random.randint(0, m)
    centerS = state[x, y]
    cost = costForCenterState(state, x, y, n, m)	
    if cost < 0:
        centerS *= -1
    elif np.random.random() < np.exp(-cost * BETA):
        centerS *= -1
    state[x, y] = centerS	
    return state

def runState(state, n_steps, snapsteps = None): 
    if snapsteps is None:
        snapsteps = np.linspace(0, n_steps, 
                     num = round(n_steps / (M * N * 100)), dtype = np.int32)
    saved_states = []
    sp = 0
    magnet_hist = []
    for i in range(n_steps):				
        state = mcmcAdjust(state) 
        magnet_hist.append(magnetizationForState(state)) 
        if sp < len(snapsteps) and i == snapsteps[sp]:
            saved_states.append(np.copy(state))		
            sp += 1
    return state, saved_states, magnet_hist

### 시뮬레이션을 실행합니다
init_state = initRandState(N, M)
print(init_state)
final_state = runState(np.copy(init_state), 1000)

[[ 1  1 -1  1 -1]
 [-1 -1 -1 -1  1]
 [ 1 -1  1  1  1]
 [ 1  1 -1  1 -1]
 [ 1 -1  1 -1 -1]]

import matplotlib.pyplot as plt

## results 리스트의 개별 요소로서 각 시계열을 수집합니다
results = []
for i in range(100):
    init_state = initRandState(N, M)
    final_state, states, magnet_hist = runState(init_state, 1000)
    results.append(magnet_hist)

## 겹치는 곡선을 알아볼 수 있게끔 약간의 투명도를 줘서 각 곡선을 그립니다.
for mh in results:
    plt.plot(mh,'r', alpha=0.2)