%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    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] = coxian2(mean,secondmoment,thirdmoment);
%
% Given the first three moments of a nonnegative distribution, G,
% coxian2 returns a 2-phase PH distribution, (tau,T) that best approximates G.

function [tau,T] = coxian2(mean,secondmoment,thirdmoment);

error = 1e-6;

if( ((2*mean^2 - secondmoment)/secondmoment < error) & ((secondmoment - 2*mean^2)/secondmoment < error) & ((6*mean^3 - thirdmoment)/thirdmoment < error) & ((thirdmoment - 6*mean^3)/thirdmoment < error) ) 
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% approximating by exponential distribution %%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 tau = 1;
 T = -1/mean;

else

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% approximating by 2-phase PH distribution %%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 %% normalized moments %%
 m2 = secondmoment / mean^2;
 m3 = thirdmoment / (mean * secondmoment);

 %% if a given distribution cannot be well-represented 
 %% by a 2-phase PH distribution,
 %% move the given distribution to the "closest" distribution 
 %% that can be well-represented by a 2-phase PH distriution.
 if( m2 < 3/2 )
  m2 = 3/2;
  m3 = 2;
 elseif( m2 > 2)
  if( m3 <= 3*m2/2 )
   m3 = 1.01*3*m2/2;
  end
 else
  if( m3 > 6*(m2-1)/m2 )
   m3 = 6*(m2-1)/m2;
   m3 = 2*m2 - 1;
  elseif( m3 < (9*m2-12+3*(2-m2)*sqrt(2*(2-m2))) / m2 )
   m3 = (9*m2-12+3*(2-m2)*sqrt(2*(2-m2))) / m2;
   m3 = 2*m2 - 1;
  end
 end

 %% finding the 2-phase PH distribution %%
 u = (6-2*m3) / (3*m2-2*m3);
 v = (12-6*m2) / (m2*(3*m2-2*m3));
 if( v < 0 )
     WarningV = 1
 end

 sqrtD = sqrt(u^2-4*v);
 if( sqrtD < 0 )
     WarningD = 1
 end

 mu1 = (u+sqrtD)/2;
 mu2 = (u-sqrtD)/2;
 p = mu2*(mu1-1)/mu1;
 if( (p<0) | (1<p) )
     WarningP = 1
 end

 MU1 = mu1 / mean;
 MU2 = mu2 / mean;

 tau = [1 0];
 T = [-MU1 p*MU1;
      0    -MU2];

end