-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathsamplingmarkov
119 lines (95 loc) · 2.57 KB
/
samplingmarkov
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
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
P=[0 0.3 0.7;0.4 0 0.6;0.5 0.5 0];
lengthChain=1000;
lambda=5;
channel = zeros(1,lengthChain); % 3-state Markov chain (output vector).
timeState = zeros(1,lengthChain); % Vector of time in each state
timeState1= zeros(1,lengthChain); % Auxiliary vector of time
channel(1,1) = channel1; % First sample of the process
P1 = cumsum(P,2); % P’3chan
for i = 2:lengthChain
eventP = randint(1,1,[1 100])/100;%Event to decide the switch
if channel(1,i-1) == 1
%If former sample was channel 1
if eventP < P1(1,2)
channel(1,i) = 2;%Switch to channel 2
else
channel(1,i) = 3;%Switch to channel 3
end
elseif channel(1,i-1) == 2%If former sample was channel 2
if eventP < P1(2,1)
channel(1,i) = 1;%Switch to channel 1
else
channel(1,i) = 3;%Switch to channel 3
end
elseif channel(1,i-1) == 3
if eventP < P1(3,1)
channel(1,i) = 1;
else
channel(1,i) = 2;
end; end; end
%If former sample was channel 3
%Switch to channel 1
%Switch to channel 2
% Sojourn time
for i = 1:lengthChain
%eventW = randint(1,1,[1 100])/100;
%Event to decide the sojourn time
if channel(i) == 1
timeState1(1,i) = poissrnd(lambda) ;
%Sojourn time when Channel 1 appears
elseif channel(i) == 2
timeState1(1,i) = poissrnd(lambda) ;
%Sojourn time when Channel 2 appears
else
timeState1(1,i) = poissrnd(lambda) ;
%Sojourn time when Channel 3 appears
end; end
timeState = cumsum(timeState1);
channel_new=downsample(channel,4);
idx1=length(find(channel_new==1));
idx2=length(find(channel_new==2));
idx3=length(find(channel_new==3));
q12 = zeros(1,lengthChain/4);
for i = 1:lengthChain/4
if (channel(1,i)==1)&(channel(1,i+1)==2)
q12(1,i)=i;
end
end
q13 = zeros(1,lengthChain/4);
for i = 1:lengthChain/4
if (channel(1,i)==1)&(channel(1,i+1)==3)
q12(1,i)=i;
end
end
q21 = zeros(1,lengthChain/4);
for i = 1:lengthChain/4
if (channel(1,i)==2)&(channel(1,i+1)==1)
q12(1,i)=i;
end
end
q23 = zeros(1,lengthChain/4);
for i = 1:lengthChain/4
if (channel(1,i)==2)&(channel(1,i+1)==3)
q12(1,i)=i;
end
end
q31 = zeros(1,lengthChain/4);
for i = 1:lengthChain/4
if (channel(1,i)==3)&(channel(1,i+1)==1)
q12(1,i)=i;
end
end
q32 = zeros(1,lengthChain/4);
for i = 1:lengthChain/4
if (channel(1,i)==3)&(channel(1,i+1)==2)
q12(1,i)=i;
end
end
p12=length(q12)/idx1;
p13=length(q12)/idx1;
p21=length(q21)/idx2;
p23=length(q23)/idx2;
p31=length(q31)/idx3;
p32=length(q32)/idx3;
P_new=[0 p12 p13;p21 0 p23;p31 p32 0]
pi_new=[idx1/250 idx2/250 idx3/250]