LINE Solver
MATLAB API documentation
Loading...
Searching...
No Matches
Coxian.m
1classdef Coxian < Markovian
2 % The coxian statistical distribution
3 %
4 % Copyright (c) 2012-2026, Imperial College London
5 % All rights reserved.
6
7 methods
8 function self = Coxian(varargin)
9 % SELF = COXIAN(VARARGIN)
10 % Constructs a Coxian distribution from phase rates and
11 % completion probabilities, with entry probability 1 on the
12 % first phase
13 self@Markovian('Coxian',2);
14 if length(varargin)==2
15 mu = varargin{1};
16 phi = varargin{2};
17 if abs(phi(end)-1)>GlobalConstants.CoarseTol && isfinite(phi(end))
18 line_error(mfilename,sprintf('The completion probability in the last Cox state must be 1.0 but it is %0.1f', phi(end)));
19 end
20 setParam(self, 1, 'mu', mu);
21 setParam(self, 2, 'phi', phi); % completion probability in phase 1
22 elseif length(varargin)==3
23 mu1 = varargin{1};
24 mu2 = varargin{2};
25 phi1 = varargin{3};
26 setParam(self, 1, 'lambda0', mu1);
27 setParam(self, 2, 'lambda1', mu2);
28 setParam(self, 3, 'phi0', phi1); % completion probability in phase 1
29 else
30 line_error(mfilename,'Coxian accepts at most 3 parameters.');
31 end
32 jline_mu = java.util.ArrayList();
33 jline_phi = java.util.ArrayList();
34 if length(self.params) == 3
35 jline_mu.add(self.getParam(1).paramValue);
36 jline_mu.add(self.getParam(2).paramValue);
37 jline_phi.add(self.getParam(3).paramValue);
38 else
39 mu = self.getParam(1).paramValue;
40 phi = self.getParam(2).paramValue;
41 for i = 1:length(mu)
42 jline_mu.add(mu(i));
43 end
44
45 for i = 1:length(phi)
46 jline_phi.add(phi(i));
47 end
48 end
49 self.obj = jline.lang.processes.Coxian(jline_mu, jline_phi);
50 if length(self.params) == 2
51 mu = getMu(self);
52 phi = getPhi(self);
53 ph = {diag(-mu)+diag(mu(1:end-1).*(1-phi(1:end-1)),1),[phi.*mu,zeros(length(mu),length(mu)-1)]};
54 else
55 mu1 = self.getParam(1).paramValue;
56 mu2 = self.getParam(2).paramValue;
57 phi1 = self.getParam(3).paramValue;
58 ph={[-mu1,(1-phi1)*mu1;0,-mu2],[phi1*mu1,0;mu2,0]};
59 end
60 self.process = ph;
61 self.nPhases = length(ph{1});
62 end
63 end
64
65 methods
66 function phases = getNumberOfPhases(self)
67 % PHASES = GETNUMBEROFPHASES()
68 % Return number of phases in the distribution
69 if length(self.params) == 2
70 phases = length(self.getParam(1).paramValue);
71 else
72 phases = 2;
73 end
74 end
75
76 function ex = getMean(self)
77 % EX = GETMEAN()
78 % Get distribution mean
79 if length(self.params) == 2
80 mu = getMu(self);
81 phi = getPhi(self);
82 ex = map_mean({diag(-mu)+diag(mu(1:end-1).*(1-phi(1:end-1)),1),[phi.*mu,zeros(length(mu),length(mu)-1)]});
83 else
84 % Get distribution mean
85 mu1 = self.getParam(1).paramValue;
86 mu2 = self.getParam(2).paramValue;
87 phi1 = self.getParam(3).paramValue;
88 ex = 1/mu1 + (1-phi1)/mu2;
89 end
90 end
91
92 function SCV = getSCV(self)
93 % SCV = GETSCV()
94 % Get the squared coefficient of variation of the distribution (SCV = variance / mean^2)
95 if length(self.params) == 2
96 mu = getMu(self);
97 phi = getPhi(self);
98 SCV = map_scv({diag(-mu)+diag(mu(1:end-1).*(1-phi(1:end-1)),1),[phi.*mu,zeros(length(mu),length(mu)-1)]});
99 else
100 mu1 = self.getParam(1).paramValue;
101 mu2 = self.getParam(2).paramValue;
102 phi1 = self.getParam(3).paramValue;
103 mean = 1/mu1 + (1-phi1)/mu2;
104 var = ((2*mu2*(mu1 - mu1*phi1))/(mu1 + mu2 - mu1*phi1) + (2*mu1*mu2*phi1)/(mu1 + mu2 - mu1*phi1))/(mu1*mu1*((mu2*(mu1 - mu1*phi1))/(mu1 + mu2 - mu1*phi1) + (mu1*mu2*phi1)/(mu1 + mu2 - mu1*phi1))) - (1/mu1 - (phi1 - 1)/mu2)*(1/mu1 - (phi1 - 1)/mu2) - (((phi1 - 1)/(mu2*mu2) + (phi1 - 1)/(mu1*mu2))*((2*mu2*(mu1 - mu1*phi1))/(mu1 + mu2 - mu1*phi1) + (2*mu1*mu2*phi1)/(mu1 + mu2 - mu1*phi1)))/((mu2*(mu1 - mu1*phi1))/(mu1 + mu2 - mu1*phi1) + (mu1*mu2*phi1)/(mu1 + mu2 - mu1*phi1));
105 SCV = var / mean^2;
106 end
107 end
108
109 function mu = getMu(self)
110 % MU = GETMU()
111 % Get vector of rates
112 if length(self.params) == 2
113 mu = self.getParam(1).paramValue(:);
114 else
115 mu1 = self.getParam(1).paramValue;
116 mu2 = self.getParam(2).paramValue;
117 mu = [mu1;mu2];
118 end
119 end
120
121 function phi = getPhi(self)
122 % PHI = GETPHI()
123 % Get vector of completion probabilities
124 if length(self.params) == 2
125 phi = self.getParam(2).paramValue(:);
126 else
127 phi1 = self.getParam(3).paramValue;
128 phi = [phi1;1.0];
129 end
130 end
131
132 end
133
134 methods(Static)
135 function cx = fitCentral(MEAN, VAR, SKEW)
136 % CX = FITCENTRAL(MEAN, VAR, SKEW)
137 cx = Cox2.fitCentral(MEAN, VAR, SKEW);
138 SCV = VAR/MEAN^2;
139 if abs(1-map_scv(cx.getProcess)/SCV) > 0.01
140 cx = Coxian.fitMeanAndSCV(MEAN, SCV);
141 end
142 cx.immediate = MEAN < GlobalConstants.CoarseTol;
143 end
144
145 function [cx,mu,phi] = fitMeanAndSCV(MEAN, SCV)
146 % [CX,MU,PHI] = FITMEANANDSCV(MEAN, SCV)
147 % Fit a Coxian distribution with given mean and squared coefficient of variation (SCV=variance/mean^2)
148 if SCV >= 1-GlobalConstants.CoarseTol && SCV <= 1+GlobalConstants.CoarseTol
149 n = 1;
150 mu = 1/MEAN;
151 phi = 1;
152 elseif SCV > 0.5+GlobalConstants.CoarseTol && SCV<1-GlobalConstants.CoarseTol
153 phi = 0.0;
154 n = 2;
155 mu = zeros(n,1);
156 mu(1) = 2/MEAN/(1+sqrt(1+2*(SCV-1)));
157 mu(2) = 2/MEAN/(1-sqrt(1+2*(SCV-1)));
158 elseif SCV <= 0.5+GlobalConstants.CoarseTol
159 n = ceil(1/SCV);
160 lambda = n/MEAN;
161 mu = lambda*ones(n,1);
162 phi = zeros(n,1);
163 else % SCV > 1+GlobalConstants.CoarseTol
164 n = 2;
165 %transform hyperexp into coxian
166 mu = zeros(n,1);
167 mu(1) = 2/MEAN;
168 mu(2) = mu(1)/( 2*SCV );
169 phi = zeros(n,1);
170 phi(1) = 1-mu(2)/mu(1);
171 phi(2) = 1;
172 end
173 phi(n) = 1;
174 cx = Coxian(mu,phi);
175 cx.immediate = MEAN < GlobalConstants.CoarseTol;
176 end
177
178 end
179
180end
181
Generated by doxygen 1.9.8