logo

에이전트 기반 모델 시뮬레이션에서의 사망 📂동역학

에이전트 기반 모델 시뮬레이션에서의 사망

시뮬레이션

malthusian\_growth\_simulation2.gif

이 포스트에서는 생성된 에이전트가 사망하는 액션을 주어 거시적인 관점에서 집단의 역성장을 모방하려고 한다. 이 시뮬레이션에서 공간이나 이동에 관련된 모든 것들은 단지 시각화를 위한 것이며, 실제 목적과는 아무런 관계가 없다.

변수

  • $t$: 현재 턴을 의미한다.
  • $N(t)$: $t$ 턴에서 에이전트의 수를 나타낸다.

파라미터

  • $N_{0} \in \mathbb{N}$: 시뮬레이션이 시작할 때 에이전트의 수를 나타낸다.
  • $d \in [0,1]$: 사망률death rate로써 에이전트가 번식할 확률을 나타낸다.
  • $t_{\text{end}} \in \mathbb{N}$: 시뮬레이션이 끝나는 턴을 정한다.

액션

(모든 에이전트가, 매 턴마다)

  • 사망: d의 확률로 시뮬레이션에서 제외된다.
  • 이동: 이변수 정규분포 $N \left( \mathbb{0} , \Sigma \right)$ 에서 뽑은 벡터대로 이동한다.

코드 리뷰

Step 1. 초기화 및 에이전트 생성

julia> N0 = 500 # 초기 인구수
500

julia> d = 0.04 # 사망률
0.04

julia> max\_iteration = 180 # 시뮬레이션 기간
180

julia> gaussian2 = MvNormal([0.0; 0.0], 0.04I) # 2차원 정규분포
IsoNormal(
dim: 2
μ: [0.0, 0.0]
Σ: [0.04 0.0; 0.0 0.04]
)

julia> Random.seed!(0);

julia> Time\_evolution = [] # 인구수를 기록하기 위한 스택
Any[]

에이전트 기반 시뮬레이션 튜토리얼과 비슷하게 에이전트를 생성한다. 여기서 각 변수들은 위의 파라미터 세팅을 반영한다.

  • N0 $\leftarrow N_{0} = 5$
  • d $\leftarrow d = 0.04$
  • max\_iteration $\leftarrow t_{\text{end}} = 180$

Random.seed!(0)은 결과를 재현하기 위한 시드넘버 설정,은 결과를 재현하기 위한 시드넘버 설정, time\_evolutiond는 개체수를 기록하기 위한 스택으로 시뮬레이션 자체에는 영향을 미치지 않으나 시뮬레이션을 반복하고 분석하기 위해 중요하게 쓰이는 테크닉이다.

julia> coordinate = rand(gaussian2, N0)'
500×2 Adjoint{Float64,Array{Float64,2}}:
  0.135821    0.165683
 -0.0706015  -0.0269708
  0.117323    0.0594672
  0.0129895  -0.0218035
 -0.102842    0.314866
 -0.137781   -0.152561
  0.0794965   0.162326
 -0.0692709  -0.0375145
  ⋮
 -0.22697     0.178074
  0.036792   -0.35028
 -0.253518    0.418463
 -0.226243   -0.111246
  0.120505   -0.0111552
  0.0882067   0.0526977
  0.398744   -0.114504
 -0.215111    0.259829

julia> N = N0
500

$N$ 은 개체수를 나타내므로 우선 $N \gets N_{0}$ 와 같이 초기 개체수를 받아 주어야한다.

Step 2. 에이전트의 사망

julia> coordinate = coordinate[(rand(N) .> d),:]
481×2 Array{Float64,2}:
  0.135821    0.165683
 -0.0706015  -0.0269708
  0.117323    0.0594672
  0.0129895  -0.0218035
 -0.102842    0.314866
  0.0794965   0.162326
 -0.0692709  -0.0375145
  ⋮
 -0.253518    0.418463
 -0.226243   -0.111246
  0.120505   -0.0111552
  0.0882067   0.0526977
  0.398744   -0.114504
 -0.215111    0.259829

현재 에이전트의 수만큼 0과 1사이에서 난수를 뽑는다. 이것이 d보다 작다는 것은 d의 확률로 사망 이 일어났다는 뜻이다. 위 스크린샷에서는 19개의 에이전트가 사망해 481개의 에이전트가 남은 것을 볼 수 있다. 이는 기존 에이전트 500개 중 약 4%가 사망한 것으로, 우리의 파라미터 세팅이 거의 정확하게 반영된 것을 확인할 수 있다.

Step 3. 기록과 에이전트 이동

julia> push!(time\_evolution, N)
1-element Array{Any,1}:
 481

julia> coordinate = coordinate + rand(gaussian2, N)'
481×2 Array{Float64,2}:
  0.399141     0.492689
 -0.151015    -0.00161691
 -0.00843121   0.10119
  0.111592     0.00865223
  0.0591775    0.150437
  0.0060825    0.0383255
 -0.0652535   -0.0688566
 -0.702124    -0.485473
  ⋮
  0.437486    -0.354416
 -0.149827     0.156366
 -0.0737965   -0.168831
  0.344358    -0.136392
 -0.231089    -0.0695535
  0.446688     0.143493
  0.133667     0.156873

위의 Step 2. 에서 에이전트가 19개 줄었으므로 $N$ 역시 업데이트 해주었고, 이 때의 개체수를 기록했다. 각 에이전트의 좌표는 행렬로 표현되고 있으므로 이변수 정규분포를 $N$ 번 뽑아 쌓은 $N \times 2$ 행렬을 더하면 에이전트가 이동한 것을 구현할 수 있다. 이를 마치고 Step 2. 로 돌아간다. 계속해서 반복하다가, $t = t_{\text{end}}$ 가 되면 멈춘다.

위 과정을 반복하기만하면 포스트 상단의 움짤과 같은 현상을 볼 수 있다.

전체 코드

다음은 이 포스트에 쓰인 줄리아 코드다.

cd(@__DIR__) # 파일 저장 경로

@time using Plots
@time using Random
@time using Distributions
@time using LinearAlgebra

N0 = 500 # 초기 인구수
d = 0.04 # 사망률
max_iteration = 180 # 시뮬레이션 기간
gaussian2 = MvNormal([0.0; 0.0], 0.04I) # 2차원 정규분포

Random.seed!(0);
time_evolution = [] # 인구수를 기록하기 위한 스택
let
  coordinate = rand(gaussian2, N0)'
  N = N0

  anim = @animate for t = (0:max_iteration)/100
    plot(coordinate[:,1], coordinate[:,2], Seriestype = :scatter,
      markercolor = RGB(1.,94/255,0.), markeralpha = 0.4, markerstrokewidth	= 0.1,
      title = "t = $t", aspect_ratio = 1, size = [400,400],
      xaxis=true,yaxis=true,axis=nothing, legend = false)
      xlims!(-10.,10.)
      ylims!(-10.,10.)

    coordinate = coordinate[(rand(N) .> d),:]

    N = size(coordinate, 1)
    push!(time_evolution, N)
    coordinate = coordinate + rand(gaussian2, N)'
  end
  gif(anim, "malthusian_growth_simulation2.gif", fps = 18)
end