-
Notifications
You must be signed in to change notification settings - Fork 0
/
example_sensitivity.m
55 lines (42 loc) · 1.3 KB
/
example_sensitivity.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
function example_sensitivity(C,Tm)
% QSN package example
%
% Sensitivity to pertubations
% SPDX-FileCopyrightText: Copyright (C) 2011-2019, 2022 Frank C Langbein <frank@langbein.org>, Cardiff University
% SPDX-FileCopyrightText: Copyright (C) 2011-2019, 2022 SM Shermer <lw1660@gmail.com>, Swansea University
% SPDX-License-Identifier: AGPL-3.0-or-later
figure(3);
T = 0:0.5:Tm+5;
T = Tm-5:0.5:Tm+5;
N = 11;
init = 1;
target = 5;
pert = 0.01;
%C = [ 10.4516 11.3318 11.6275 11.3318 10.4516 229.1179 67.8795 55.7833 55.7833 67.8795 229.1179 ];
for r = 1:100
J = 1 * (1 + pert * (rand(1,N) * 2 - 1));
chain = qsn.QSN('ring', J, 'XX', C);
P = chain.prob(T);
for in = 1:N
if in == init
p = arrayfun(@(t) P{t}(in,target),1:size(P,2));
plot (T,p,'b');
elseif in ~= target
p = arrayfun(@(t) P{t}(in,target),1:size(P,2));
plot (T,p,'m');
end
hold on;
end
end
J = 1 * (1 + pert * ones(1,N));
chain = qsn.QSN('ring', J, 'XX', C);
P = chain.prob(T);
p = arrayfun(@(t) P{t}(init,target),1:size(P,2));
plot (T,p,'r');
J = 1 * (1 - pert * ones(1,N));
chain = qsn.QSN('ring', J, 'XX', C);
P = chain.prob(T);
p = arrayfun(@(t) P{t}(init,target),1:size(P,2));
plot (T,p,'g');
hold off;
end