HK (Hegselmann and Krause) opinion dynamics model
This example showcases
- How to do synchronous updating of Agent properties (also know as Synchronous update schedule). In a Synchronous update schedule changes made to an agent are not seen by other agents until the next step, see also Wilensky 2015, p.286).
- How to terminate the system evolution on demand according to a boolean function.
- How to terminate the system evolution according to what happened on the previous step.
Model overview
This is an implementation of a simple version of the Hegselmann and Krause (2002) model. It is a model of opinion formation with the question: which parameters' values lead to consensus, polarization or fragmentation? It models interacting groups of agents (as opposed to interacting pairs, typical in the literature) in which it is assumed that if an agent disagrees too much with the opinion of a source of influence, the source can no longer influence the agent's opinion. There is then a "bound of confidence". The model shows that the systemic configuration is heavily dependent on this parameter's value.
The model has the following components:
- A set of n Agents with opinions xᵢ in the range [0,1] as attribute
- A parameter ϵ called "bound" in (0, 0.3]
- The update rule: at each step every agent adopts the mean of the opinions which are within the confidence bound ( |xᵢ - xⱼ| ≤ ϵ).
It is also available from the Models
module as Models.hk
.
Core structures
We start by defining the Agent type and initializing the model. The Agent type has two fields so that we can implement the synchronous update.
using Agents
using Statistics: mean
mutable struct HKAgent <: AbstractAgent
id::Int
old_opinion::Float64
new_opinion::Float64
previous_opinon::Float64
end
There is a reason the agent has three fields that are "the same". The old_opinion
is used for the synchronous agent update, since we require access to a property's value at the start of the step and the end of the step. The previous_opinion
is the opinion of the agent in the previous step, as the model termination requires access to a property's value at the end of the previous step, and the end of the current step.
We could, alternatively, make the three opinions a single field with vector value.
function hk_model(; numagents = 100, ϵ = 0.2)
model = ABM(HKAgent, scheduler = fastest, properties = Dict(:ϵ => ϵ))
for i in 1:numagents
o = rand()
add_agent!(model, o, o, -1)
end
return model
end
model = hk_model()
AgentBasedModel with 100 agents of type HKAgent no space scheduler: fastest properties: Dict(:ϵ => 0.2)
Add some helper functions for the update rule. As there is a filter in the rule we implement it outside the agent_step!
method. Notice that the filter is applied to the :old_opinion
field.
function boundfilter(agent, model)
filter(
j -> abs(agent.old_opinion - j) < model.ϵ,
[a.old_opinion for a in allagents(model)],
)
end
Now we implement the agent_step!
function agent_step!(agent, model)
agent.previous_opinon = agent.old_opinion
agent.new_opinion = mean(boundfilter(agent, model))
end
and model_step!
function model_step!(model)
for a in allagents(model)
a.old_opinion = a.new_opinion
end
end
From this implementation we see that to implement synchronous scheduling we define an Agent type with old
and new
fields for attributes that are changed via the synchronous update. In agent_step!
we use the old
field then, after updating all the agents new
fields, we use the model_step!
to update the model for the next iteration.
Running the model
The parameter of interest is now :new_opinion
, so we assign it to variable adata
and pass it to the run!
method to be collected in a DataFrame.
In addition, we want to run the model only until all agents have converged to an opinion. From the documentation of step!
one can see that instead of specifying the amount of steps we can specify a function instead.
function terminate(model, s)
if any(
!isapprox(a.previous_opinon, a.new_opinion; rtol = 1e-12) for a in allagents(model)
)
return false
else
return true
end
end
step!(model, agent_step!, model_step!, terminate)
model[1]
Main.ex-hk.HKAgent(1, 0.6832036080050999, 0.6832036080050999, 0.6832036080050998)
Alright, let's wrap everything in a function and do some data collection using run!
.
function model_run(; kwargs...)
model = hk_model(; kwargs...)
agent_data, _ = run!(model, agent_step!, model_step!, terminate; adata = [:new_opinion])
return agent_data
end
data = model_run(numagents = 100)
data[(end - 19):end, :]
step | id | new_opinion | |
---|---|---|---|
Int64 | Int64 | Float64 | |
1 | 27 | 81 | 0.465129 |
2 | 27 | 82 | 0.465129 |
3 | 27 | 83 | 0.465129 |
4 | 27 | 84 | 0.465129 |
5 | 27 | 85 | 0.465129 |
6 | 27 | 86 | 0.465129 |
7 | 27 | 87 | 0.465129 |
8 | 27 | 88 | 0.465129 |
9 | 27 | 89 | 0.465129 |
10 | 27 | 90 | 0.465129 |
11 | 27 | 91 | 0.465129 |
12 | 27 | 92 | 0.465129 |
13 | 27 | 93 | 0.465129 |
14 | 27 | 94 | 0.465129 |
15 | 27 | 95 | 0.465129 |
16 | 27 | 96 | 0.465129 |
17 | 27 | 97 | 0.465129 |
18 | 27 | 98 | 0.465129 |
19 | 27 | 99 | 0.465129 |
20 | 27 | 100 | 0.465129 |
Notice that here we didn't speciy when to collect data, so this is done at every step. Instead, we could collect data only at the final step, by re-using the same function for the when
argument:
model = hk_model()
agent_data, _ = run!(
model,
agent_step!,
model_step!,
terminate;
adata = [:new_opinion],
when = terminate,
)
agent_data
step | id | new_opinion | |
---|---|---|---|
Int64 | Int64 | Float64 | |
1 | 6 | 1 | 0.681177 |
2 | 6 | 2 | 0.681177 |
3 | 6 | 3 | 0.180512 |
4 | 6 | 4 | 0.681177 |
5 | 6 | 5 | 0.681177 |
6 | 6 | 6 | 0.681177 |
7 | 6 | 7 | 0.180512 |
8 | 6 | 8 | 0.681177 |
9 | 6 | 9 | 0.180512 |
10 | 6 | 10 | 0.681177 |
11 | 6 | 11 | 0.681177 |
12 | 6 | 12 | 0.180512 |
13 | 6 | 13 | 0.681177 |
14 | 6 | 14 | 0.681177 |
15 | 6 | 15 | 0.180512 |
16 | 6 | 16 | 0.681177 |
17 | 6 | 17 | 0.681177 |
18 | 6 | 18 | 0.681177 |
19 | 6 | 19 | 0.180512 |
20 | 6 | 20 | 0.681177 |
21 | 6 | 21 | 0.180512 |
22 | 6 | 22 | 0.180512 |
23 | 6 | 23 | 0.180512 |
24 | 6 | 24 | 0.180512 |
25 | 6 | 25 | 0.681177 |
26 | 6 | 26 | 0.681177 |
27 | 6 | 27 | 0.180512 |
28 | 6 | 28 | 0.180512 |
29 | 6 | 29 | 0.681177 |
30 | 6 | 30 | 0.681177 |
31 | 6 | 31 | 0.681177 |
32 | 6 | 32 | 0.681177 |
33 | 6 | 33 | 0.180512 |
34 | 6 | 34 | 0.681177 |
35 | 6 | 35 | 0.180512 |
36 | 6 | 36 | 0.180512 |
37 | 6 | 37 | 0.681177 |
38 | 6 | 38 | 0.681177 |
39 | 6 | 39 | 0.180512 |
40 | 6 | 40 | 0.681177 |
41 | 6 | 41 | 0.681177 |
42 | 6 | 42 | 0.180512 |
43 | 6 | 43 | 0.681177 |
44 | 6 | 44 | 0.681177 |
45 | 6 | 45 | 0.681177 |
46 | 6 | 46 | 0.180512 |
47 | 6 | 47 | 0.681177 |
48 | 6 | 48 | 0.681177 |
49 | 6 | 49 | 0.681177 |
50 | 6 | 50 | 0.180512 |
51 | 6 | 51 | 0.681177 |
52 | 6 | 52 | 0.681177 |
53 | 6 | 53 | 0.180512 |
54 | 6 | 54 | 0.681177 |
55 | 6 | 55 | 0.681177 |
56 | 6 | 56 | 0.681177 |
57 | 6 | 57 | 0.180512 |
58 | 6 | 58 | 0.681177 |
59 | 6 | 59 | 0.681177 |
60 | 6 | 60 | 0.681177 |
61 | 6 | 61 | 0.681177 |
62 | 6 | 62 | 0.681177 |
63 | 6 | 63 | 0.180512 |
64 | 6 | 64 | 0.180512 |
65 | 6 | 65 | 0.180512 |
66 | 6 | 66 | 0.180512 |
67 | 6 | 67 | 0.180512 |
68 | 6 | 68 | 0.681177 |
69 | 6 | 69 | 0.180512 |
70 | 6 | 70 | 0.681177 |
71 | 6 | 71 | 0.681177 |
72 | 6 | 72 | 0.180512 |
73 | 6 | 73 | 0.681177 |
74 | 6 | 74 | 0.681177 |
75 | 6 | 75 | 0.681177 |
76 | 6 | 76 | 0.681177 |
77 | 6 | 77 | 0.681177 |
78 | 6 | 78 | 0.180512 |
79 | 6 | 79 | 0.681177 |
80 | 6 | 80 | 0.681177 |
81 | 6 | 81 | 0.180512 |
82 | 6 | 82 | 0.180512 |
83 | 6 | 83 | 0.681177 |
84 | 6 | 84 | 0.681177 |
85 | 6 | 85 | 0.681177 |
86 | 6 | 86 | 0.681177 |
87 | 6 | 87 | 0.681177 |
88 | 6 | 88 | 0.180512 |
89 | 6 | 89 | 0.180512 |
90 | 6 | 90 | 0.180512 |
91 | 6 | 91 | 0.180512 |
92 | 6 | 92 | 0.681177 |
93 | 6 | 93 | 0.180512 |
94 | 6 | 94 | 0.681177 |
95 | 6 | 95 | 0.180512 |
96 | 6 | 96 | 0.681177 |
97 | 6 | 97 | 0.180512 |
98 | 6 | 98 | 0.180512 |
99 | 6 | 99 | 0.681177 |
100 | 6 | 100 | 0.180512 |
Finally we run three scenarios, collect the data and plot it.
using Plots
plotsim(data, ϵ) = plot(
data.step,
data.new_opinion,
leg = false,
group = data.id,
title = "epsilon = $(ϵ)",
)
plt001, plt015, plt03 =
map(e -> (model_run(ϵ = e), e) |> t -> plotsim(t[1], t[2]), [0.05, 0.15, 0.3])
plot(plt001, plt015, plt03, layout = (3, 1))