Set the precision and initialize butools (load all packages)

In [1]:

```
%precision %g
%run "~/github/butools/Python/BuToolsInit.py"
```

First the global *butools.verbose* flag is set to True to obtain more messages from the functions.

In [2]:

```
butools.verbose = True
```

The BuTools PH package offers tools for both phase-type (PH) and matrix-exponential (ME) distributions. Some functions expect a PH representation for the input (MomentsFromPH, PdfFromPH, CdfFromPH, etc.), and some others expect an ME representation (MomentsFromME, PdfFromME, CdfFromME, etc.).

If the global flag called *butools.checkInput* is set to True, these functions enforce the proper representation of their input parameters (the corresponding tolerance is given by the global *butools.checkPrecision* variable).

The following example calls the checking functions on a distribution.

In [3]:

```
alpha = ml.matrix([[0.2, 0.3, 0.5]])
A = ml.matrix([[-1, 0, 0],[0, -3, 1],[0, -1, -3]])
```

In [4]:

```
CheckPHRepresentation(alpha,A)
```

Out[4]:

In [5]:

```
CheckMERepresentation(alpha,A)
```

Out[5]:

Calling a function that needs a PH representation will fail:

In [6]:

```
MomentsFromPH(alpha,A)
```

However, the ME counterpart of this function runs properly:

In [7]:

```
MomentsFromME(alpha,A)
```

Out[7]:

In [8]:

```
alpha, A = PH2From3Moments ([1,3,20])
print("alpha=", alpha)
print("A=", A)
```

It can be checked that the moments of the result are the same as the target moments.

In [9]:

```
MomentsFromPH(alpha, A)
```

Out[9]:

The extension of the previous matching procedure: matching five moments with a PH(3)

In [10]:

```
alpha, A = PH3From5Moments ([0.9,2.5,20, 500, 22000])
print("alpha=", alpha)
print("A=", A)
```

In [11]:

```
MomentsFromPH(alpha, A)
```

Out[11]:

*graphviz* tool from AT&T is installed, the resulting PH distribution can be visualized by *ImageFromPH*.

In [12]:

```
ImageFromPH(alpha,A)
```

Out[12]:

**APHFrom3Moments** is a flexible procedure. It is able to matching *any* three moments with an appropriately large APH (acyclic PH). Let us first try to match moments [1,1.28,8] with PH2From3Moments. It fails.

In [14]:

```
alpha, A = PH2From3Moments ([1,1.28,8])
```

In [15]:

```
alpha, A = APHFrom3Moments ([1,1.28,8])
print("alpha=", alpha)
print("A=", A)
```

In [16]:

```
MomentsFromME(alpha, A)
```

Out[16]:

In [17]:

```
def feasibleBounds (n):
ly = []
uy = []
x = []
for m2 in np.linspace(APH2ndMomentLowerBound(1,n), 1.8, 500):
# convert to normalized moments
ml = NormMomsFromMoms ([1, m2, APH3rdMomentLowerBound(1,m2,n)])
mu = NormMomsFromMoms ([1, m2, APH3rdMomentUpperBound(1,m2,n)])
# record bounds
ly.append(ml[2])
uy.append(min(mu[2],15))
x.append(ml[1])
return (x, ly, uy)
plt.ylim((1.2,3))
plt.xlim((1.1,1.8))
plt.fill_between (*feasibleBounds(6),color="blue")
plt.fill_between (*feasibleBounds(5),color="yellow")
plt.fill_between (*feasibleBounds(4),color="green")
plt.fill_between (*feasibleBounds(3),color="red")
plt.fill_between (*feasibleBounds(2),color="orange")
plt.xlabel("$n_2(=m_2/m_1^2)$")
plt.ylabel("$n_3(=m_3/m_1/m_2)$");
```

**MEFromMoments** returns a vector-matrix pair matching any number of moments. Unfortunatly, while the moments are mathcing indeed according to the moment formulas, the density function can be negative, thus *MeFromMoments* can return invalid ME distributions.

In [18]:

```
v, H = MEFromMoments ([0.9,2.5,20, 500, 22000])
print("v=", v)
print("H=", H)
```

In [19]:

```
MomentsFromME(v, H)
```

Out[19]:

**PHFromME** looks for a solution having the same size. Warning: it is not always possible to find a non-markovian representation of the same size. It can happen that only a larger Markovian representation exists.

In [20]:

```
beta, B = PHFromME (v, H)
print("beta=", beta)
print("B=", B)
```

*PHFromME* terminates successfully. To show that $(\beta,B)$ and $(v,H)$ define the same distribution, we plot it and conclude that they match.

In [21]:

```
x = np.linspace(0, 2, 100)
y1 = PdfFromME (v, H, x)
y2 = PdfFromPH (beta, B, x)
plt.plot(x,y1,x,y2);
```

*SimilarityMatrix* function of the *RepTrans* package.

In [22]:

```
T = SimilarityMatrix(B,H)
T
```

Out[22]:

Let us now check if $\beta,B$ and $v,H$ are equivalent:

In [23]:

```
la.norm(T.I*B*T - H)
```

Out[23]:

In [24]:

```
la.norm(beta*T - v)
```

Out[24]:

Now let us take an other example $v,H$, where no Markovian representation of the same size exists.

In [25]:

```
v = ml.matrix([[0.2, 0.3, 0.5]])
H = ml.matrix([[-1, 0, 0],[0, -3, 1],[0, -1, -3]])
beta, B = PHFromME (v, H)
print("beta=", beta)
print("B=", B)
```

**MonocyclicPHFromME**. This function is a very strong tool. It is able to convert *any* ME distribution (ok, only those that do not touch the 0 axis apart from point 0) to a PH distribution.

In [26]:

```
beta, B = MonocyclicPHFromME (v, H)
print("beta=", beta)
print("B=", B)
```

Altough it is larger (9 states instead of 3), but at least Markovian.

In [27]:

```
CheckPHRepresentation(beta,B)
```

Out[27]:

Let us check the cdfs now to see if they match.

In [28]:

```
x = np.linspace(0, 2, 100)
y1 = CdfFromME (v, H, x)
y2 = CdfFromPH (beta, B, x)
plt.plot(x,y1,x,y2);
```

*MonocyclicPHFromME* is not interesting, there is *CheckMEPositiveDensity* available to decide if the conversion to PH is possible or not, i.e. if
the given vector-matrix pair has a non-negative density or not.

In [29]:

```
CheckMEPositiveDensity(v,H)
```

Out[29]:

In [30]:

```
gammac, Gc = MinimalRepFromME(beta,B,"cont")
print("controllability order = ", Gc.shape[0])
gammao, Go = MinimalRepFromME(beta,B,"obs",1e-4)
print("observability order = ", Go.shape[0])
gamma, G = MinimalRepFromME(beta,B,"obscont") # this is the default, "obscont" can be omitted
print("minimal order = ", G.shape[0])
print("Minimal representation:")
print("gamma=", gamma)
print("G=", G)
```

**MEOrderFromMoments**. Based on a set of moments it determines the minimal order of ME distribution that has those moments. For instance, it is often easy to obtain moments from several queueing models: this method can check if there is an underlying PH or ME distribution or not.

In [31]:

```
MEOrderFromMoments(MomentsFromPH(beta,B))
```

Out[31]:

In [32]:

```
alpha,A = RandomPH(4, 1.2, 8)
print("alpha=",alpha)
print("A=",A)
```

*SamplesFromPH* function returns a vector of random samples from a PH distribution. The last parameter is the number of samples requested. The samples can be used in simulations, for instance.

In [33]:

```
tr = SamplesFromPH(alpha,A,1000000);
```

In [34]:

```
np.mean(tr)
```

Out[34]: