%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (c) 2004. Takayuki Osogami. All Rights Reserved. %
% %
% Permission to use, copy, modify, and distribute this source code %
% and its documentation for any purpose, without fee, and without %
% a written agreement, is hereby granted, provided that the above %
% copyright notice, this paragraph and the following two paragraphs %
% appear in all copies, modifications, and distributions. %
% Created by Takayuki Osogami, Department of Computer Science, %
% Carnegie Mellon University. %
% %
% IN NO EVENT SHALL THE AUTHOR BE LIABLE TO ANY PARTY FOR ANY %
% DAMAGES, INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS %
% SOURCE CODE AND ITS DOCUMENTATION, EVEN IF THE AUTHOR HAS BEEN %
% ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. %
% %
% THE AUTHOR SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, %
% BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY %
% AND FITNESS FOR A PARTICULAR PURPOSE. THE SOURCE CODE AND %
% ACCOMPANYING DOCUMENTATION, IF ANY, PROVIDED HEREUNDER IS %
% PROVIDED "AS IS." THE AUTHOR HAS NO OBLIGATION TO PROVIDE %
% MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Usage: [tau,T] = matching3PH(momentsofG)
%
% INPUT, momentsofG, is a vector of size 3.
% An example of an input is:
%
% momentsofG = [2,6,24];
%
% OUTPUT, [tau,T], is a pair of a vector and a matrix.
% When the input is [2,6,24], the output is
%
% tau = [1,0];
%
% T = [-1 1;
% 0 -1];
%
% Yes. This is an Erlang-2 distribution with mean 2.
%
% Given first three moments, momentsofG,
% the function returns a PH distribution, (tau,T), whose
% first three moments is momentsofG.
% The output PH distribution has no mass probability at zero.
% That is the sum of all the elements of tau is 1.
%
% The function is not defined everywhere, but almost everywhere.
% When the function is not defined for the input, it perturbs
% the input and returns a PH distribution that has the first three
% moments close to the input. The error is very small for most
% practical purposes.
%
% When the output is (0,0), there was an error in the function
% or the input was illegal.
%
% author: Takayuki Osogami
% Department of Computer Science
% Carnegie Mellon University
% osogami@cs.cmu.edu
% date: June 6, 2004
function [tau,T] = matching3PH(momentsofG)
%% Maximum number of phases used in the PH distribution.
%% If the number of necessary phases exceeds this limit,
%% the input is matched by the closest PH distribution.
MAX_PHASE = 50;
EPSILON = 0.001;
%% If the input moment is zero, return error.
if( (momentsofG(1) == 0) || (momentsofG(2) == 0) )
Illegal_input = momentsofG
tau = 0;
T = 0;
return;
end
%% Calculate the normalized moments of the input, G.
m2 = momentsofG(2) / momentsofG(1)^2;
m3 = momentsofG(3) / (momentsofG(1) * momentsofG(2));
%% Any PH distribution has m2 > 1 and m3 > m2.
%% If there is no PH for the input, return error.
if( (m2 <= 1 + EPSILON^2) || (m3 <= m2 + EPSILON^2) )
Illegal_input = momentsofG
tau = 0;
T = 0;
return;
end
%% If the necessary number of phases exceeds the limit, MAX_PHASE,
%% calculate the normalized moments of the "closest" PH distribution.
smallestM2 = (MAX_PHASE/(MAX_PHASE-1) + (MAX_PHASE-1)/(MAX_PHASE-2)) / 2;
if( m2 < smallestM2 )
m2 = smallestM2;
end
smallestM3 = (MAX_PHASE+2)/(MAX_PHASE+1) * m2;
if( m3 < smallestM3 )
m3 = smallestM3;
end
%% If there is no PH distribution with the input moments,
%% calculate the normalized moments of the "closest" PH distribution.
if( (m3 < 2*m2 - 1) && (abs(m3 - 3/2 * m2) < EPSILON^2) )
m3 = (1+EPSILON) * m3;
elseif( m3 > 2*m2 - 1 )
r = m3 / (2*m2-1);
n = floor(1/(m2-1));
if( abs(1/(m2-1) - n) < EPSILON^2 )
m2 = 1 + 1 / (n*(1-EPSILON));
m3 = r * (2*m2-1);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% The matching PH is exponential %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if( (abs(m2-2) < EPSILON^2) && (abs(m3-3) < EPSILON^2) )
tau = 1;
T = - 1 / momentsofG(1);
return;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% The matching PH is 2-phase PH %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if( isPH2(m2,m3) )
[tau,T] = matching3PH2(momentsofG(1),m2,m3);
return
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% The matching PH is EC (Erlang-Coxian) %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if( m3 >= 2*m2 - 1 )
[tau,T] = matching3EC(momentsofG(1),m2,m3);
return;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% The matching PH is "p * EC + (1-p) * EXP" %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
k = floor((2*m2-m3)/(m3-m2) + EPSILON^2);
if( m3 >= ((k+1)*m2 + (k+4)) / (2*(k+2)) * m2 )
y = (2-m2) / (4*(3/2-m3/m2));
p = (2-m2)^2 / ((2-m2)^2+4*(2*m2-1-m3));
m2X = 2 * y;
m3X = 2 * m2X - 1;
[tauX,TX] = matching3EC(1,m2X,m3X);
[tau,T] = mixtureofPH(tauX,TX,1,-1/y,p);
%% normalize
mean = momentofPH(tau,T,1);
T = T * mean / momentsofG(1);
return;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% The matching PH is "EC + EXP" %%
%% EC has mass probability at zero %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if( abs(m2-2) < EPSILON^2 )
y = (m3-2*(k+3)/(k+2)) / (3-m3);
m2X = 2*(1+y);
else
A = m2 * ((m3-3) - 2*(k+3)/(k+2)*(m2-2));
B = m2 * sqrt((m3-3)^2 + 8*(k+3)/(k+2)*(m2-2)*(3/2-m3/m2));
C = 2 * (k+3)/(k+2) * (m2-2)^2;
y = (A+B) / C;
m2X = (1+y) * (m2*(1+y)-2*y);
%% This is a magic spell for Matlab, who doesn't like unnecessary precision...
%% m2X = round(m2X * 1000) / 1000;
end
m3X = (k+3)/(k+2) * m2X;
[tauX,TX] = matching3EC(1,m2X,m3X);
[tau,T] = convolutionofPH(tauX,TX,1,-1/y);
%% normalize
mean = momentofPH(tau,T,1);
T = T * mean / momentsofG(1);
return;