Quality analysis in acyclic production networks

TL;DR

This paper develops statistical estimators for comparing machine quality contributions in acyclic production networks, enabling the identification of unreliable machines through pairwise and multiple comparisons.

Contribution

It introduces new estimators for means and variances of machine quality contributions and discusses their asymptotic distributions for reliable statistical testing.

Findings

Estimators enable effective comparison of machine quality metrics.

Asymptotic distributions facilitate standard statistical tests.

Method helps identify unreliable machines in production networks.

Abstract

The production network under examination consists of a number of workstations. Each workstation is a parallel configuration of machines performing the same kind of tasks on a given part. Parts move from one workstation to another and at each workstation a part is assigned randomly to a machine. We assume that the production network is acyclic, that is, a part does not return to a workstation where it previously received service. Furthermore, we assume that the quality of the end product is additive, that is, the sum of the quality contributions of the machines along the production path. The contribution of each machine is modeled by a separate random variable. Our main result is the construction of estimators that allow pairwise and multiple comparison of the means and variances of machines in the same workstation. These comparisons then may lead to the identification of unreliable machines. We also discuss the asymptotic distributions of the estimators that allow the use of standard statistical tests and decision making.

Equations

$$V_{\vec{p}}:=(V_{\vec{p}}(i,j){)}_{i\in [r],j\in [c]},$$

$$S:=(s(i,j){)}_{i\in [r],j\in [c]},s(i,j)\in \mathbb{R}.$$

$$b(\vec{p})=\sum _{j=1}^{c}S(p_j,j).$$

$$B^{(n)}:=\sum _{k=1}^{n}b({\vec{p}}_k)V_{{\vec{p}}_k},\text{leavevmode}{V}^{(n)}:=\sum _{k=1}^n{V}_{{\vec{p}}_k},n\ge 1.$$

$$T^{(n)}:=(T^{(n)}(i,j){)}_{i\in [r],j\in [c]},\text{mbox}where{T}^{(n)}(i,j):=\{\begin{array}{ll}\frac{B^{(n)}(i,j)}{V^{(n)}(i,j)},& \text{if }{V}^{(n)}(i,j)\ne 0;\\ \text{0, }& \text{otherwise}.\end{array}$$

$${\Sigma }^{(n)}(i,j):=\{\begin{array}{ll}\frac{1}{V^{(n)}(i,j)}{\sum }_{k=1}^n{\left(b({\vec{p}}_k)V_{{\vec{p}}_k}(i,j)-T^{(n)}(i,j)\right)}^2,& \text{if }{V}^{(n)}(i,j)\ne 0;\\ \text{0, }& \text{otherwise}.\end{array}$$

$$D_{(i,j)}^{(n)}:=\{\vec{p}\in D^{(n)}:\vec{p}[j]=i\}$$

$$T^{(n)}(i,j)-T^{(n)}(i^{\prime },j)\underset{n\to \infty }{\overset{a.s.}{\to }}\mathbb{E}[S(i,j)]-\mathbb{E}[S(i^{\prime },j)]$$

$${\Sigma }^{(n)}(i,j)-{\Sigma }^{(n)}(i^{\prime },j)\underset{n\to \infty }{\overset{a.s.}{\to }}\mathbb{V}[S(i,j)]-\mathbb{V}[S(i^{\prime },j)].$$

$$T^{(n)}(i,j)$$

$$T^{(n)}(i^{\prime },j)=\frac{1}{|D_{(i^{\prime },j)}^{(n)}|}\sum _{\vec{p}\in D_{(i^{\prime },j)}^{(n)}}b(\vec{p})\underset{n\to \infty }{\overset{a.s.}{\to }}\mathbb{E}\left[b(\vec{\mathcal{P}})|\vec{\mathcal{P}}[j]=i^{\prime }\right].$$

$$\mathbb{E}\left[b(\vec{\mathcal{P}})|\vec{\mathcal{P}}[j]=i\right]-\mathbb{E}\left[b(\vec{\mathcal{P}})|\vec{\mathcal{P}}[j]=i^{\prime }\right]$$

$${\Sigma }^{(n)}(i,j)$$

$$\mathbb{E}\left[b(\vec{\mathcal{P}}{)}^2|\vec{\mathcal{P}}[j]=i\right]-\mathbb{E}{\left[b(\vec{\mathcal{P}})|\vec{\mathcal{P}}[j]=i\right]}^2=A_j+\mathbb{V}\left[S(i,j)\right]$$

$${\Sigma }^{(n)}(i,j)-{\Sigma }^{(n)}(i^{\prime },j)\underset{n\to \infty }{\overset{a.s.}{\to }}\mathbb{V}\left[S(i,j)\right]-\mathbb{V}\left[S(i^{\prime },j)\right].$$

$$\sum _{\vec{p}\in D_{(i,j)}^{(n)}}b(\vec{p})=\sum _{k=1}^{n}b({\vec{\mathcal{P}}}_k)1\{{\vec{\mathcal{P}}}_k[j]=i\}$$

$$\mathbb{P}(Y\le y)=\sum _{k=1}^{\infty }\mathbb{P}(b({\vec{\mathcal{P}}}_k)\le y|K=k)\mathbb{P}(K=k)=\mathbb{P}(b({\vec{\mathcal{P}}}_1)\le y|{\vec{\mathcal{P}}}_1[j]=i).$$

$$\frac{1}{|D_{(i,j)}^{(n)}|}\sum _{k=1}^{|D_{(i,j)}^{(n)}|}{Y}_k,$$

$$\sqrt{|D_{(i,j)}^{(n)}|}\left(T^{(n)}(i,j)-{\mu }_{i,j}\right)\underset{n\to \infty }{\overset{\mathcal{D}}{\to }}\mathcal{N}(0,{\sigma }_{i,j}^2),$$

$${\mu }_{i,j}:=\mathbb{E}\left[b(\vec{\mathcal{P}})|\vec{\mathcal{P}}[j]=i\right]\text{mbox}and{\sigma }_{i,j}^2:=\mathbb{E}\left[b(\vec{\mathcal{P}}{)}^2|\vec{\mathcal{P}}[j]=i\right]-\mathbb{E}{\left[b(\vec{\mathcal{P}})|\vec{\mathcal{P}}[j]=i\right]}^2.$$

$${\Sigma }^{(n)}(i,j)=\frac{1}{|D_{(i,j)}|}\sum _{\vec{p}\in D_{(i,j)}^{(n)}}{\left(b(\vec{p})-\frac{1}{|D_{(i,j)}^{(n)}|}\sum _{\vec{p}\in D_{(i,j)}^{(n)}}b(\vec{p})\right)}^2.$$

$${\Sigma }^{(n)}(i,j)=\frac{1}{|D_{(i,j)}|}\sum _{\vec{p}\in D_{(i,j)}^{(n)}}b(\vec{p}{)}^2-{\left(T^{(n)}(i,j)\right)}^2.$$

$$\sqrt{D_{(i,j)}^{(n)}}\left({\Sigma }^{(n)}(i,j)-{\sigma }_{i,j}^2\right)$$

$$\frac{\sqrt{|D_{(i,j)}^{(n)}|}}{{\sigma }_{i,j}}(T^{(n)}(i,j)-{\mu }_{i,j})\sim \mathcal{N}(0,1)\text{mbox}and\frac{\sqrt{|D_{i^{\prime },j}^{(n)}|}}{{\sigma }_{i^{\prime },j}}(T^{(n)}(i^{\prime },j)-{\mu }_{i^{\prime },j})\sim \mathcal{N}(0,1).$$

$${\left(\frac{|D_{i,j}^{(n)}|\cdot |D_{i^{\prime },j}^{(n)}|}{|D_{i,j}^{(n)}|+|D_{i^{\prime },j}^{(n)}|}\right)}^{1/2}\frac{T^{(n)}(i,j)-T^{(n)}(i^{\prime },j)}{\sigma }\sim \mathcal{N}(0,1).$$

$$\hat{\sigma }={\left(\frac{1}{|D_{i,j}^{(n)}|+|D_{i^{\prime },j}^{(n)}|}\left(|D_{i,j}^{(n)}|{\Sigma }^{(n)}(i,j)+|D_{i^{\prime },j}^{(n)}|{\Sigma }^{(n)}(i^{\prime },j)\right)\right)}^{1/2},$$

$${\left(\frac{|D_{i,j}^{(n)}|\cdot |D_{i^{\prime },j}^{(n)}|}{|D_{i,j}^{(n)}|+|D_{i^{\prime },j}^{(n)}|}\right)}^{1/2}\frac{T^{(n)}(i,j)-T^{(n)}(i^{\prime },j)}{\hat{\sigma }}$$

$$X_k=T^{(n)}(k,j)-T^{(n)}(r_j,j)k\in \{1,\dots ,r_j-1\}.$$

$$\mathbb{E}[S]:=\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 2& 0& 0\\ 0& 0& \ast & 0\\ 0& \ast & \ast & 0\end{array}\right]\text{mbox}and\mathbb{V}[S]:=\left[\begin{array}{cccc}1& 1& 4& 1\\ 1& 1& 1& 1\\ 1& 1& \ast & 1\\ 1& \ast & \ast & 1\end{array}\right],$$

$$T=\left[\begin{array}{cccc}0.60& -0.03& 0.75& 0.70\\ 0.25& 2.06& 0.40& 0.89\\ 0.74& -0.10& \ast & 0.38\\ 0.69& \ast & \ast & 0.14\end{array}\right],\text{hfill}\Sigma =\left[\begin{array}{cccc}6.46& 5.62& 8.92& 7.12\\ 8.76& 6.13& 5.38& 4.89\\ 4.97& 6.58& \ast & 8.41\\ 7.37& \ast & \ast & 8.14\end{array}\right].$$

Full text

Our main result is the construction of estimators that allow pairwise and multiple comparison of the means and variances of machines in the same workstation. These comparisons then may lead to the identification of unreliable machines. We also discuss the asymptotic distributions of the estimators that allow the use of standard statistical tests and decision making.

Keywords: direct acyclic graphs, production networks, quality estimation, anomaly detection, variability

AMS MSC 2010: 90B30, 90B15, 62M02, 62M05

1 Introduction

In order to maintain competitiveness, industrial manufacturers have to pursue new sources to enhance process agility. In parallel production networks, the quality of the end product depends on various intermediate production steps. Therefore, a recent challenge is to achieve a level of visibility into the production flows that allow to optimize throughput by guaranteeing at the same time given quality standards. The challenge to maintain a visibility across all parallel workflows and to identify eventual sources of errors becomes particularly difficult in situations where the qualities of single machines are not observable.

In this note we discuss a new approach allowing to compare the impact of different machines of the same workstation, i.e. performing the same task in parallel, on the quality of the product. We assume that the quality of the product is only observable at the end of the process and that the qualities of the nodes along the path are latent variables. See Figure 1 for an illustration. Although the main interest is usually in the process mean, the size of the process variability is often crucial. Our approach allows to model the process mean and process variability in parallel under very general conditions. In particular, the only condition on the production network is that it is acyclic, i.e. a part does not return to a workstation where it previously received service. These kind of networks are modeled in mathematics and computer science as directed acyclic graphs (DAGs).

The concept of DAGs is becoming increasingly important having various applications in different fields of science and engineering. For instance DAGs have many applications not only in modeling production networks, e.g. [8], but also in process networks and distributed computing, e.g. [6], scheduling for systems of tasks, e.g. [17], DAG-networks for deep learning, e.g. [18], and DAG-based alternatives of the blockchain technology, e.g. [14]. Moreover, directed acyclic graphs (DAGs) are more and more often used to represent causal relationships among random variables in complex systems, e.g. [13]. We speak in general about production networks, but due to the wide applications of DAGs we want to emphasize here that our approach can be applied to any kind of acyclic network where effects along a path are additive and where measurements can only be performed at the end of the process.

Parameter estimation is not only important for the decision making process, but it is also an essential prerequisite in order to obtain meaningful simulations. However, complex production networks depend on a high number of parameters and their estimation is challenging. Since simulations in manufacturing plays an increasingly important role, e.g. [10, 12], we believe that our approach is an interesting contribution to this topic as well.

Previous work

Estimation for the mean can also be conducted using a linear regression with categorical covariates, e.g. [11]. However, since in a linear regression homoscedasticity, i.e. the homogeneity of the variances, is a crucial condition, this approach naturally does not allow to compare differences of the variances. In fact, linear regression models are often plagued by different variabilities or heteroskedasticity. We refer to [9, Chapter 4] for an overview on how to detect and control heteroskedasticity. In contrast to these problems, our approach naturally allows different variabilities between variables and is able to detect differences in variability of values of a given variable.

Multifactor experimental designs, [5, 16] are also alternatives to estimate the mean differences but also rely on homoscedasticity. They are mostly used in the context of statistically planned experiments, which consists of a few experimental runs to obtain data on the product characteristics. If the number of observations for each setting (or path in our notation) is sufficiently high and under further conditions described in [5, Section 4] these methods allow a comparison of the variances, too. While this may offer a feasible, however not direct, way to identify differences in variability if the number of machines is small, it seems not practical in more complex networks.

There is also a connection to critical paths analysis, see e.g. [3, 15]. While these methods allow to find critical paths in acyclic networks they are not suited to compare nor estimate differences in mean and variances of given tasks.

2 The model

A directed acyclic graph (DAG) is a finite directed graph with no directed cycles. It consists of a finite vertex set $V$ and a finite set of directed edges $E=\{(v,w):v,w\in V,v\neq w\}$. In our setting the DAG contains two special vertices: a source $s$ and a sink $t$. We are interested in the paths from the source to the sink in this graph. We denote a path $\vec{p}$ in the DAG as $\vec{p}=(p_{0},p_{1},\ldots,p_{c},p_{c+1})$ where $p_{0}=s$ and $p_{c+1}=t$ and $(p_{i},p_{i+1})\in E$. We define $\vec{p}[j]:=p_{j}$. We refer to Figure 2 for an illustration and to [1] for more details on directed graphs.

We assume that at each step $1\leq i\leq c$ the path $\vec{p}$ has $r_{i}\leq r$ different choices and the nodes in each column are always numerated starting with $1$. The possible choices of a path can therefore be modeled through an $r\times c$ matrix. More precisely, given a path $\vec{p}$, we associate an $r\times c$ binary matrix $V_{\vec{p}}$ that has $1$’s only in the nodes visited by the path:

[TABLE]

where we denote $[k]:=\{1,2,\ldots,k\}$ for an integer $k$. We call $V_{\vec{p}}$ the indicator matrix of the path $\vec{p}$.

Each path contains exactly one node of each column. The aim of this paper is to study differences among nodes of the same columns. We think of nodes in the same columns as different possibilities for a given task, as different persons performing the same job, as different machines in the same workstation or as variations of the same kind of treatment.

The given data consists of the list of paths $\{\vec{p}_{i}\}_{i=1,\ldots,n}$ in the DAG and the list of outputs $\{b(\vec{p}_{i})\}_{i=1,\ldots,n}$ for each path.

We consider the quality matrix $S$, which is a random matrix of size $r\times c$ with real entries

[TABLE]

We model the paths with a random vector $\vec{\mathcal{P}}:=(P_{1},\ldots,P_{c})$ where the components $P_{i}$ are random variables over the set $[r]$.

Throughout the paper we work under the following standing assumptions:

Assumption 1**.**

We assume that:

all entries of $S$ have finite second moments; 2. 2.

the random variables $S(i,j),i\in[r],\,j\in[c]$ are (jointly) independent; 3. 3.

the paths $\vec{\mathcal{P}}_{1},\vec{\mathcal{P}}_{2},\ldots$ are chosen independently and uniformly.

Note that, we do not assume the entries of $S$ to be identically distributed nor having the same variance.

Let $\vec{p}=(p_{0},p_{1},\ldots,p_{c},p_{c+1})$ be a realization of $\vec{\mathcal{P}}$, where $p_{0}=s$ and $p_{c+1}=t$ almost surely. Then, the quality of the construction path $\vec{p}$ is defined as

[TABLE]

We also can think of $b(\vec{p})$ as the quality (or error) cumulated along the path $\vec{p}$.

Let us make precise where the randomness enters in our model. We choose a random path $\vec{\mathcal{P}}$ and random matrix $S$. The corresponding probability measure is denoted by $\mathbb{P}$. The random choice of $\vec{\mathcal{P}}$ and $S$ induces a random variable $b(\vec{\mathcal{P}})$ and allows to generate a sequence of i.i.d. random variables $(\vec{\mathcal{P}}_{1},b(\vec{\mathcal{P}}_{1})),(\vec{\mathcal{P}}_{2},b(\vec{\mathcal{P}}_{2})),\ldots$.

The goal of our study is to give estimates on the law of $S$ by observing the paths $\vec{\mathcal{P}}$ and its cumulated qualities $b(\vec{\mathcal{P}})$. Note that, $(\vec{\mathcal{P}}_{n},b(\vec{\mathcal{P}}_{n}))_{n\in\mathbb{N}}$ is in general not a sufficient statistic for $S$, i.e. we can not recover the distribution of $S$ by only observing realizations of $(\vec{\mathcal{P}},b(\vec{\mathcal{P}}))$, as we see in the next remark

Remark 1**.**

Let us consider the case $r=1$ and $c=2$. Let $S(1,1)\sim\mathcal{N}(0,1)$ and $S(1,2)\sim\mathcal{N}(1,1)$ and define $\tilde{S}(1,1):=S(1,1)+1$ and $\tilde{S}(1,2):=S(1,2)-1$. Then for any given path $\vec{p}$ we have that $\sum_{j=1}^{2}S(p_{j},j)=\sum_{j=1}^{2}\tilde{S}(p_{j},j)$. Hence, the statistic $(\vec{p}_{n},b(\vec{p}_{n}))_{n\in\mathbb{N}}$ does not allow us to distinguish between $S$ and $\tilde{S}$.

Example 2** (Binary errors).**

The matrix $S$ consists of independent Bernoulli random variables $S(i,j)$. The value $1$ of this Bernoulli may encode a defect and hence $b(\vec{p})$ counts the number of defects of the end product.

Example 3** (Gaussian quality).**

The matrix $S$ consists of independent Gaussian random variables $S(i,j)$. The end quality $b(\vec{p})$ is then distributed as a mixture of Gaussian random variables.

Given a sequence of realizations $(\vec{p}_{k})_{k\in[n]}$ of $\vec{\mathcal{P}}$, we define the following matrices that at are the core of our analysis:

[TABLE]

The value $B^{(n)}(i,j)$ is the sum of all cumulated qualities of paths containing node $(i,j)$, whereas $V^{(n)}(i,j)$ just counts the number of times node $(i,j)$ was used. We define the sample mean matrix as the sample mean quality matrix :

[TABLE]

The corresponding sample variance matrix $\Sigma^{(n)}$ is defined by

[TABLE]

3 Results

Denote $D^{(n)}=(\vec{\mathcal{P}}_{i})_{i\in[n]}$ the multi-set111We use a multi-set since we need to keep track of the multiplicity of the paths. or sequence of all paths up to time $n$ and let

[TABLE]

be the multi-set of all paths up to index $n$ that go through the node $(i,j)$; we note that we can recover $V^{(n)}(i,j)$ through $D^{(n)}_{(i,j)}$ by $V^{(n)}(i,j)=|D^{(n)}_{(i,j)}|$.

In general it is not possible to estimate the mean quality matrix $\mathbb{E}[S]$, see Remark 1. However, it is possible to identify nodes with higher or lower quality mean or variance in each column.

Theorem 4**.**

Let $(i,j),(i^{\prime},j)\in[r]\times[c]$, then

[TABLE]

and

[TABLE]

Proof.

Using twice the law of large numbers and the continuous mapping theorem we obtain

[TABLE]

In the same way

[TABLE]

Using the assumption that the paths are chosen uniformly and the definition of $b(\vec{\mathcal{P}})$, we obtain the first part of the theorem from

[TABLE]

For the second part of the theorem, we use the law of large numbers and the continuous mapping theorem to get that

[TABLE]

Using the definition of $b(\vec{p})$ we deduce with elementary calculations that

[TABLE]

where $A_{j}$ is a quantity that only depend on the column $j$. Applying this identity for $(i,j)$ and $(i^{\prime},j)$ we obtain that

[TABLE]

∎

3.1 Asymptotic distribution

The sum

[TABLE]

can be interpreted as the sum of random variables appearing in an acception-rejection method. More precisely, we start with $k=1$ and consider $(\vec{\mathcal{P}}_{k},b(\vec{\mathcal{P}}_{k}))$. If $\vec{\mathcal{P}}_{k}[j]=i$ we set $Y:=b(\vec{\mathcal{P}}_{k})$ and stop, otherwise we increase $k$ and repeat until $\vec{\mathcal{P}}_{K}[j]=i$ for the first $K$. Now, for $y\in\mathbb{R}$,

[TABLE]

In other words, the distribution of $Y$ equals the distribution of $b(\vec{\mathcal{P}}_{1})$ conditioned on $\vec{\mathcal{P}}_{1}[j]=i$. Iterating this acception-rejection method we see that $|D^{(n)}_{(i,j)}|$ describes the number of acceptions using $(\vec{\mathcal{P}}_{k},b(\vec{\mathcal{P}}_{k}))$, $1\leq k\leq n$. Hence the estimator $T^{(n)}$ has the same distribution as

[TABLE]

where $Y_{k},k\in\mathbb{N},$ is a sequence of i.i.d. random variables distributed as $b(\vec{\mathcal{P}}_{1})$ conditioned on $\vec{\mathcal{P}}_{1}[j]=i$. Finally, Anscombe’s theorem, [7, Theorem 1.3.1], implies that

[TABLE]

where

[TABLE]

For the variance we assume that the entries of $S$ have finite forth moments. The estimator for the variance is

[TABLE]

Since the distribution of $\Sigma^{(n)}(i,j)$ does not change if we replace $b(\vec{p})$ by $b(\vec{p})-\mu_{i,j}$ we can assume that $\mu_{i,j}=0$. Moreover,

[TABLE]

We have that $T^{(n)}(i,j)$ converges almost surely to [math] and $\sqrt{D_{(i,j)}^{(n)}}T^{(n)}(i,j)$ converges in distribution to $\mathcal{N}(0,\sigma_{i,j}^{2})$. Slutzky’s theorem implies that $\sqrt{D_{(i,j)}^{(n)}}T^{(n)}(i,j)^{2}$ converges in distribution to [math] and therefore also in probability to [math]. Using Slutzky’s theorem and Anscombe’s theorem as above we obtain that

[TABLE]

with $\tau_{i,j}^{2}=\mathbb{E}\left[(b(\vec{\mathcal{P}})-\mu_{i,j})^{4}\,|\,\vec{\mathcal{P}}[j]=i\right]-\mathbb{E}\left[(b(\vec{\mathcal{P}})-\mu_{i,j})^{2}\,|\,\vec{\mathcal{P}}[j]=i\right]^{2}.$

3.2 Pairwise and multiple comparison

The asymptotic distributions of the estimators $T^{(n)}$ and $\Sigma^{(n)}$ justify that pairwise comparison of the mean quality and the variance of quality can be done using standard tests, e.g. $t$-test, Fisher tests and Bartlett tests. We sketch only the construction of a test statistic to compare the means of two nodes under the hypothesis that they have the same variance. However, it is standard to extent this result to the case where the variances are not equal and to the estimators of the variance differences. More details on possible applications are given in Section 4. We have for $n$ sufficiently large

[TABLE]

Under the hypothesis that $\mu_{i,j}=\mu_{i^{\prime},j}$ and $\sigma=\sigma_{i,j}=\sigma_{i^{\prime},j}$ one obtains that for $n$ sufficiently large

[TABLE]

Replacing $\sigma$ by its estimator

[TABLE]

we find that

[TABLE]

can be approximated by a standard normal distribution for $n$ sufficiently large.

The identification of nodes with bad quality boils down to multiple comparisons. For a given column $j$ with $r_{j}$ nodes we define a vector $X=(X_{1},\ldots,X_{r_{j}-1})$ of test statistics

[TABLE]

The vector $X$ satisfies (asymptotically) the positive regression dependency. Therefore, we suggest the Benjamini-Yekutieli method, see [2], to control the false discovery rate.

4 Examples

4.1 Wafer production

Our study was motivated by a root-cause analysis in the wafer fabrication. Wafer fabrication is in general a procedure of many repeated sequential processes. For instance, a simplified illustration consists of $12$ subsequent fabrication steps, see [4], where intermediate measurement of qualities are not a feasible. In our concrete examples we treated up to $30$ different steps and more than $90$ machines. Unfortunately, since our industrial partner insists on the fulfillment of an NDA we are not allowed to publish any more information about the project. Probably for the same reasons, it was impossible for us to find public available data on other industrial projects.

4.2 Simulations

We consider the DAG network as given in Figure 2 and consider Gaussian qualities as described in Example 3. The matrix $S$ consists in this case of independent Gaussian random variables. In this simulation we consider the distribution that is characterized by

[TABLE]

where the $*$ are placeholders for the machines that do not exist. We simulate $n=200$ observations and obtain the following results for the estimators $T$ and $\Sigma$ (rounded to two decimals):

[TABLE]

While the difference in mean in the second column seems to be obvious, the difference of the variances in the third column might be overlooked and differences in the variance in the first and last columns could be suspected.

We perform pairwise $t$-tests for each column with Benjamini-Yekutieli adjustment. While in all but the second column no statistically significant difference is detected, the difference in mean of the second machine in column $2$ is detected with a $p$-value of $9.8\cdot 10^{-6}$. The Bartlett test does not find any differences in the variances in columns $1,2,$ and $4$. However, the difference of variance in column $3$ seems to be statistically significant with a $p$-value of $0.01$. The number of observations of $200$ is rather small compared to the number of different paths, which is $96$. This explains the fact that the matrix $\Sigma$ does not reflect the correct differences of variances; in particular we have in the third column that $8.92-5.38=3.54>3=4-1$. In order to demonstrate the convergence of the estimators we calculate the estimators for $n=10000$:

[TABLE]

Acknowledgment

The authors wish to thank Alessandro Chiancone, Herwig Friedl, Jérôme Depauw, and Marc Peigné for stimulating discussing during this project. A. G. acknowledges financial support from the Austrian Science Fund project FWF P29355- N35. Grateful acknowledgement is made for hospitality from TU-Graz where the research was carried out during visits of S. M.