Scientific Research

An Academic Publisher

Robust Variance Components Estimation in the PERG Mixed Distributions of Empirical Variances—PEROBVC Method ()

Keywords

Share and Cite:

*Open Journal of Statistics*,

**10**, 640-650. doi: 10.4236/ojs.2020.104038.

1. Introduction

A detailed presentation of Variance and Covariance Components Models can be found in Rao and Kleffe [1], as well as in Krishnaiah [2].

However, in the applications deviations appear from the assumptions under which these solutions are obtained, which also leads to negative estimates of variance components. The main source of these deviations is the contaminating distribution of observations (also presence of outliers), which, as a rule, always appears. Therefore, Variance Components (VC) solutions are always studied, which are less sensitive or insensitive to these deviations. Such properties are characteristic for robust solutions, out of which the following ones are here presented:

Maximum Likelihood Approaches to Variance Component Estimation and to Related Problems [3]; Robust Estimation of Variance Components [4], Estimation and Use of Variance Components [5], Robust Variance Estimation in Linear Regression [6], Robust Variance Estimation for Random Effects Meta-analysis [7], Robust Variance Estimation in Meta-regression with Dependent Effect Size Estimates [8], Robust Estimation of Linear Mixed Models [9].

Robust estimates of variances and of their components have been continually applied in various fields [10] - [21] (a more detailed presentation can be seen in [22]).

Regarding a GPS measurement error as a random process and modeling the process realization structurally as a 2-way hierarchical classification with random effects, Perović [23] studied a robust solution for variance components in GPS measurements—PERG2FH method. Perović [22] also robustly estimated the parameters for both distributions of normally distributed observations that form Tukey’s mixed distribution.

From the said above it is seen that robust methods for estimation of Variance Components are studied. However, studies of Variance Components estimates on the basis of empirical variances distributions, except the standard procedure, have not been undertaken. For this reason the present author has devoted himself to VC research on the basis of these distributions, so that the here presented method was invented.

In the present paper a mixed distribution composed of two distributions of empirical variances is considered, where one of them is basic and the other one is contaminating, with an essential assumption that the variances of both distributions are obtained on the basis of primary normal distributions, as well as that they have the same degrees of freedom. The present author names this distribution PERG mixed-distribution of empirical variances (PERG—from initial letters of the author’s names: Perović Gligorije), which is presented in Figure 1 and described by means of the following model

$F\left({s}^{2}\right)=\left(1-\epsilon \right){F}_{1}\left({s}^{2}\right)+\epsilon {F}_{2}\left({s}^{2}\right),\text{\hspace{0.17em}}\text{\hspace{0.17em}}\left(0<\epsilon <0.5\right)$ (1)

^{1}PEROBVC is an abbreviation of the initial letters: Perović’s Robust Variance Components Estimation Model.

where F_{1} is the basic distribution, F_{2} the contaminating distribution, both of empirical variances *s*^{2} with *f* d.f. and ε (0 < ε < 0.5) the contaminating degree.

On the basis of such a composition of the mixed distribution of empirical variances the author obtained a robust solution named Perović’s Robust Variance Components Estimation Model—PEROBVC^{1} method for variances and for observation numbers for both distributions. In this way the method eliminates the influences of gross errors and outliers on the estimates of distribution parameters.

The key difference between the PEROBVC method and the existing ones is that distribution censoring is not used, but instead on the basis of structural decomposition in two distributions, except obtaining parameter estimates of the two distributions, the question of outliers is simultaneously solved.

Consequently, the PEROBVC method has the following properties:

1) Unbiased (exact) variance estimates for the basic distribution, as the most important property,

2) Unbiased (exact) variance estimates for the contaminating distribution,

3) Estimates of observation numbers for both distributions and in this way percentage estimates concerning fractions of basic and contaminating distributions in the mixed one.

The correctness of the PEROBVC method has been verified on exact (expected) values of some quantities from the mixed PERG distribution, as well as on examples of measurements in 2D control geodetic networks.

The variance estimate for the basic distribution obtained by applying the PEROBVC method is compared to its ML estimate.

The structure of the further presentation is the following: 1) establishing the PEROBVC method, 2) way of solving the formulated problem, 3) presentation of robust ML estimation method, 4) correctness verification for PEROBVC method by comparing its solution to the ML one and 5) verification on examples of geodetic measurements in to 2D control networks.

2. Basis of PEROBVC Method

The idea of PEROBVC method has been presented in [24] (Section 32.14).

Consider a PERG mixed distribution (Figure 1) of empirical variances ${{s}^{\prime}}_{i}^{2},i=1,2,\cdots ,{n}^{\prime}$ and ${{s}^{\u2033}}_{j}^{2},j=1,2,\cdots ,{n}^{\u2033}$

$F\left({s}^{2}\right)=\left(1-\epsilon \right)F\left({{s}^{\prime}}^{2}\right)+\epsilon F\left({{s}^{\u2033}}^{2}\right),\text{\hspace{0.17em}}\text{\hspace{0.17em}}\left(0<\epsilon <0.5\right)$ (2)

with ${n}^{\prime}+{n}^{\u2033}=n$ and with the same degrees of freedom, ${{f}^{\prime}}_{i}={{f}^{\u2033}}_{j}\equiv f$ . The variances are obtained from primary normally distributed observations with expectations

Figure 1. Density plot for PERG mixed distribution *f *(*s*^{2}) of empirical variances with 4 d. *f*. (*f* = 4); with
${f}_{1}\left({{s}^{\prime}}^{2}\right)$ is the basic with
${\sigma}_{1}^{2}$ and
${f}_{2}\left({{s}^{\u2033}}^{2}\right)$ is the contaminating distribution with
${\sigma}_{2}^{2}$ (here
${\sigma}_{2}^{2}=2{\sigma}_{1}^{2}$ ) and contaminating degree ε = 0.40.

$E\left[{{s}^{\prime}}_{i}^{2}\right]={\sigma}_{1}^{2},\text{\hspace{0.17em}}i=1,2,\cdots ,{n}^{\prime},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}E\left[{{s}^{\u2033}}_{j}^{2}\right]={\sigma}_{2}^{2},\text{\hspace{0.17em}}j=1,2,\cdots ,{n}^{\u2033}$ ; ${\sigma}_{2}^{2}>{\sigma}_{1}^{2}$ ,

where $F\left({{s}^{\prime}}^{2}\right)$ and $F\left({{s}^{\u2033}}^{2}\right)$ are the distribution functions of the empirical variances ${{s}^{\prime}}^{2}$ and ${{s}^{\u2033}}^{2}$ . The task is to estimate the variances ${\sigma}_{1}^{2}$ and ${\sigma}_{2}^{2}$ of both distributions on the basis of the mixed distribution with f known. In addition to the solving of this task, the method also yields estimates of the numbers of observations ${n}^{\prime}$ and ${n}^{\u2033}$ for both component distributions $F\left({{s}^{\prime}}^{2}\right)$ and $F\left({{s}^{\u2033}}^{2}\right)$ , i.e. the contaminating degree is estimated.

The point B divides the domain of ${s}^{2}$ into two sub-domains X and Y (Figure 1) and let be:

${n}^{\prime}$ _{}and
${n}^{\u2033}$ —numbers of all observations, (
${{s}^{\prime}}^{2}$ and
${{s}^{\u2033}}^{2}$ ), within X and Y;

${{n}^{\prime}}_{X}$ _{}and
${{n}^{\u2033}}_{X}$ —numbers of observations for
${{s}^{\prime}}^{2}$ and
${{s}^{\u2033}}^{2}$ within X;

${{n}^{\prime}}_{Y}$ and ${{n}^{\u2033}}_{Y}$ —numbers of observations for ${{s}^{\prime}}^{2}$ and ${{s}^{\u2033}}^{2}$ within Y,

$\begin{array}{l}{n}_{X}+{n}_{Y}=n\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{{n}^{\prime}}_{X}+{{n}^{\u2033}}_{X}={n}_{X}\\ {n}^{\prime}+{n}^{\u2033}=n\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{{n}^{\prime}}_{Y}+{{n}^{\u2033}}_{Y}={n}_{Y}\end{array}\}$ , (3)

The point B will be chosen in such a way that within the interval (0, B) the observations
${{s}^{\prime}}^{2}$ dominate, whereas within [B, ¥) the observations
${{s}^{\u2033}}^{2}$ dominate from where it follow that it must be B = B_{opt} (also see Figure 1).

Here the point B is defined according as

$B={s}_{\left({n}_{X}\right)}^{2}+d/2$ , (4)

where d is the width of the rounding interval for observations s^{2}.

For the sake of simplifying the following designations are introduced:

$\begin{array}{l}(a)\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}w={s}^{2}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{with}\text{\hspace{0.17em}}\text{\hspace{0.17em}}f\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{d}\text{.f}\text{.}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\Rightarrow \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{w}^{\prime}={{s}^{\prime}}^{2},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{w}^{\u2033}={{s}^{\u2033}}^{2}\\ (b)\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}y=f\frac{{s}^{2}}{{\sigma}^{2}}={\chi}_{f}^{2}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\Rightarrow \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{y}^{\prime}=f\frac{{{s}^{\prime}}^{2}}{{\sigma}_{1}^{2}},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{y}^{\u2033}=f\frac{{{s}^{\u2033}}^{2}}{{\sigma}_{2}^{2}}\end{array}\}$ (5)

Then the density of the random variable_{
${s}^{2}$ }-distribution can be written as

$f\left({s}^{2}\right)=f\left(w\right)={k}_{f}^{-1}{\sigma}^{-f}{f}^{f/2}{w}^{f/2-1}{\text{e}}^{-wf/2{\sigma}^{2}}=\frac{f}{{\sigma}^{2}}f\left(y\right)$ , ${k}_{f}={2}^{f/2}\Gamma \left(\frac{f}{2}\right)$ (6)

where $\Gamma \left(f/2\right)$ is the gamma function of f/2. Then

$F\left({s}^{2}\right)=F\left(w\right)={\displaystyle {\int}_{0}^{w}f\left(w\right)\text{d}w}={\displaystyle {\int}_{0}^{y}f\left(y\right)\text{d}y}=F\left(y\right)$ (7)

$wf\left(w\right)\text{d}w=\frac{{\sigma}^{2}}{f}yf\left(y\right)\text{d}y\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\Rightarrow \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\displaystyle {\int}_{0}^{w}wf\left(w\right)\text{d}w}=\frac{{\sigma}^{2}}{f}{\displaystyle {\int}_{0}^{y}yf\left(y\right)\text{d}y}$ . (8)

From the equality condition for the ordinates of the two distribution densities, i.e. from the condition, one obtains optimal point B (B_{opt})

${B}_{\text{opt}}=\frac{\frac{2}{f}\mathrm{ln}\left(\frac{{n}^{\prime}}{{n}^{\u2033}}{\left(\frac{{\sigma}_{2}}{{\sigma}_{1}}\right)}^{f}\right)}{\frac{1}{{\sigma}_{1}^{2}}-\frac{1}{{\sigma}_{2}^{2}}}.$ (9)

3. PEROBVC Solution

The estimates of the variance components for both distributions are obtained from the event-probability maximum

$\left[{{D}^{\prime}}_{\left(1\right)}\wedge \cdots \wedge {{D}^{\prime}}_{\left({{n}^{\prime}}_{X}\right)}\right]\wedge \left({{s}^{\u2033}}^{2}<B\right)\wedge \left({{s}^{\prime}}^{2}>B\right)\wedge \left[{{D}^{\u2033}}_{\left({{n}^{\u2033}}_{X}+1\right)}\wedge \cdots \wedge {{D}^{\u2033}}_{\left({n}^{\u2033}\right)}\right]$

where ${{D}^{\prime}}_{\left(1\right)}={{s}^{\prime}}_{\left(1\right)}^{2}<{{s}^{\prime}}^{2}<{{s}^{\prime}}_{\left(1\right)}^{2}+\text{d}{{s}^{\prime}}^{2}={{w}^{\prime}}_{\left(1\right)}<{w}^{\prime}<{{w}^{\prime}}_{\left(1\right)}+\text{d}{w}^{\prime}$ , $\cdots $ , ${{D}^{\u2033}}_{\left({n}^{\u2033}\right)}={{s}^{\u2033}}_{\left({n}^{\u2033}\right)}^{2}<{{s}^{\u2033}}^{2}<{{s}^{\u2033}}_{\left({n}^{\u2033}\right)}^{2}+\text{d}{{s}^{\u2033}}^{2}={{w}^{\u2033}}_{\left({n}^{\u2033}\right)}<{w}^{\u2033}<{{w}^{\u2033}}_{\left({n}^{\u2033}\right)}+\text{d}{w}^{\u2033}$ , whose likelihood function, to the proportionality constant k, is

$L=k\cdot {\displaystyle {\prod}_{X}f\left({{w}^{\prime}}_{\left(i\right)}\right)\cdot p\left({w}^{\u2033}<B\right)\cdot p\left({w}^{\prime}>B\right)}\cdot {\displaystyle {\prod}_{Y}f\left({{w}^{\u2033}}_{\left(i\right)}\right)}$ ,

where: ${\prod}_{X}f\left({{w}^{\prime}}_{\left(i\right)}\right)}={\displaystyle {\prod}_{1}^{{{n}^{\prime}}_{X}}f\left({{w}^{\prime}}_{\left(i\right)}\right)$ ,.

In view of (5)-(8) one finds:

$\begin{array}{ll}{{f}^{\prime}}_{B}={k}_{f}^{-1}{{y}^{\prime}}_{B}^{f/2-1}{\text{e}}^{-{{y}^{\prime}}_{B}/2}\hfill & {{f}^{\u2033}}_{B}={k}_{f}^{-1}{{y}^{\u2033}}_{B}^{f/2-1}{\text{e}}^{-{{y}^{\u2033}}_{B}/2}\hfill \\ {{y}^{\prime}}_{B}=\frac{f\cdot B}{{\sigma}_{1}^{2}}\hfill & {{y}^{\u2033}}_{B}=\frac{f\cdot B}{{\sigma}_{2}^{2}}\hfill \\ {{F}^{\prime}}_{B}={\displaystyle {\int}_{0}^{{{y}^{\prime}}_{B}}f\left({y}^{\prime}\right)\text{d}{y}^{\prime}}\hfill & {{F}^{\u2033}}_{B}={\displaystyle {\int}_{0}^{{{y}^{\u2033}}_{B}}f\left({y}^{\u2033}\right)\text{d}{y}^{\u2033}}\hfill \\ {{A}^{\prime}}_{B}={\displaystyle {\int}_{0}^{{{y}^{\prime}}_{B}}\frac{1}{f}{y}^{\prime}f\left({y}^{\prime}\right)\text{d}{y}^{\prime}}\hfill & {{A}^{\u2033}}_{B}={\displaystyle {\int}_{0}^{{{y}^{\u2033}}_{B}}\frac{1}{f}{y}^{\u2033}f\left({y}^{\u2033}\right)\text{d}{y}^{\u2033}}\hfill \end{array}\text{}\}$ . (10)

Furthermore lnL will be used there is:

$\mathrm{ln}L=c+{\displaystyle {\sum}_{X}\mathrm{ln}\text{\hspace{0.05em}}f\left({{w}^{\prime}}_{\left(i\right)}\right)}+\mathrm{ln}{\left({{F}^{\u2033}}_{B}\right)}^{{{n}^{\u2033}}_{X}}+\mathrm{ln}{\left(1-{{F}^{\prime}}_{B}\right)}^{{{n}^{\prime}}_{Y}}+{\displaystyle {\sum}_{Y}\mathrm{ln}f\left({{w}^{\u2033}}_{\left(i\right)}\right)}$ , (11)

where c is a constant.

According to the asymptotic theory it is:

$\sum}_{X}{{w}^{\u2033}}_{i}}\underset{{n}^{\u2033}\to \infty}{\overset{p}{\to}}{{n}^{\u2033}}_{X}{\displaystyle {\int}_{0}^{B}{w}^{\u2033}f\left({w}^{\u2033}\right)\text{d}{w}^{\u2033$ , (12a)

${\sum}_{Y}{{w}^{\prime}}_{i}}\underset{{n}^{\prime}\to \infty}{\overset{p}{\to}}{{n}^{\prime}}_{Y}{\displaystyle {\int}_{0}^{B}{w}^{\prime}f\left({w}^{\prime}\right)\text{d}{w}^{\prime}}={{n}^{\prime}}_{Y}{\sigma}_{1}^{2}\left(1-{{A}^{\prime}}_{B}\right)$ . (12b)

Bearing in mind (10)-(12) and introducing the substitutions

$\begin{array}{l}(\text{a})\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\displaystyle {\sum}_{X}{{s}^{\prime}}_{i}^{2}}={\displaystyle {\sum}_{X}{s}_{i}^{2}}-{\displaystyle {\sum}_{X}{{s}^{\u2033}}_{i}^{2}}\text{,}{\displaystyle {\sum}_{Y}{{s}^{\u2033}}_{i}^{2}}={\displaystyle {\sum}_{Y}{s}_{i}^{2}}-{\displaystyle {\sum}_{Y}{{s}^{\prime}}_{i}^{2}}\\ (\text{b})\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\displaystyle {\sum}_{X}{{s}^{\u2033}}_{i}^{2}}={n}^{\u2033}{\displaystyle {\int}_{X}{{s}^{\u2033}}^{2}f\left({{s}^{\u2033}}^{2}\right)\text{d}{{s}^{\u2033}}^{2}}={n}^{\u2033}{\sigma}_{2}^{2}{{A}^{\u2033}}_{B}\\ (\text{c})\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\displaystyle {\sum}_{Y}{{s}^{\prime}}_{i}^{2}}={n}^{\prime}{\displaystyle {\int}_{Y}{{s}^{\prime}}^{2}f\left({{s}^{\prime}}^{2}\right)\text{d}{{s}^{\prime}}^{2}}={n}^{\prime}{\sigma}_{1}^{2}\left(1-{{A}^{\prime}}_{B}\right)\end{array}\}$

from the necessary maximum condition: (a) $\partial \mathrm{ln}L/\partial {\sigma}_{1}^{2}=0$ and (b) $\partial \mathrm{ln}L/\partial {\sigma}_{2}^{2}=0$ we get the equations

$\begin{array}{l}(\text{a})\text{}\frac{-{{n}^{\prime}}_{X}}{2{\sigma}_{1}^{2}}+\frac{1}{2{\sigma}_{1}^{4}}{\displaystyle \underset{X}{\sum}{s}_{i}^{2}}-\frac{{{n}^{\u2033}}_{X}{\sigma}_{2}^{2}}{2{\sigma}_{1}^{4}}{{A}^{\u2033}}_{B}+\frac{B{{n}^{\prime}}_{Y}{{f}^{\prime}}_{B}}{{\sigma}_{1}^{4}\left(1-{{F}^{\prime}}_{B}\right)}=0\\ (\text{b})\text{}\frac{-{{n}^{\u2033}}_{Y}}{2{\sigma}_{2}^{2}}+\frac{1}{2{\sigma}_{2}^{4}}{\displaystyle \underset{Y}{\sum}{s}_{i}^{2}}-\frac{B{{n}^{\u2033}}_{X}{{f}^{\u2033}}_{B}}{{\sigma}_{2}^{4}{{F}^{\u2033}}_{B}}=0\end{array}\}$

Their solution is given as simple iterations

${\sigma}_{1,k+1}^{2}=\frac{1}{{{n}^{\prime}}_{X,k}}\left[{\displaystyle {\sum}_{X}{s}_{i}^{2}}-{{n}^{\u2033}}_{k}{\sigma}_{2,k}^{2}\text{\hspace{0.05em}}{{A}^{\u2033}}_{B,k}+2B\frac{{{n}^{\prime}}_{Y,k}{{f}^{\prime}}_{B,k}}{1-{{F}^{\prime}}_{B,k}}\right]$ , (13a)

${\sigma}_{2,k+1}^{2}=\frac{1}{{{n}^{\u2033}}_{Y,k}}\left[{\displaystyle {\sum}_{Y}{s}_{i}^{2}}-{{n}^{\prime}}_{k}{\sigma}_{1,k}^{2}\left(1-{{A}^{\prime}}_{B,k}\right)-2B\frac{{{n}^{\u2033}}_{X,k}{{f}^{\u2033}}_{B,k}}{{{F}^{\u2033}}_{B,k}}\right]$ . (13b)

where the numbers ${{n}^{\prime}}_{X}$ , ${{n}^{\prime}}_{Y}$ , ${{n}^{\u2033}}_{X}$ and ${{n}^{\u2033}}_{Y}$ are determined from the expressions

$\begin{array}{l}(a)\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{{n}^{\prime}}_{X}={n}^{\prime}{{F}^{\prime}}_{B},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{{n}^{\prime}}_{Y}={n}^{\prime}-{{n}^{\prime}}_{X}\\ (b)\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{{n}^{\u2033}}_{X}={n}^{\u2033}{{F}^{\u2033}}_{B},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{{n}^{\u2033}}_{Y}={n}^{\u2033}-{{n}^{\u2033}}_{X}\end{array}\}$ , (14)

whereas
${n}^{\prime}$ _{ }and
${n}^{\u2033}$ are determined as

$\begin{array}{l}{{F}^{\prime}}_{B}{n}^{\prime}+{{F}^{\u2033}}_{B}{n}^{\u2033}={n}_{X}\\ \left(1-{{F}^{\prime}}_{B}\right){n}^{\prime}+\left(1-{{F}^{\u2033}}_{B}\right){n}^{\u2033}={n}_{Y}\end{array}\}$ (15)

Note 1. The initial values for the variance components can be assumed, for example, as:

${\sigma}_{1,0}^{2}=0.6\cdot \left(\frac{1}{n}{\displaystyle \underset{1}{\overset{n}{\sum}}{s}_{i}^{2}}\right)$ and ${\sigma}_{2,0}^{2}=3\cdot \left(\frac{1}{n}{\displaystyle \underset{1}{\overset{n}{\sum}}{s}_{i}^{2}}\right)$ .

If the model assumptions are not satisfied, for example, if the distribution instead of a long tail has a shortened tail (from the right), then for one of the variances we obtain negative estimates. ∆

Note 2. Blunders in the observations (measurements with large gross errors) must be rejected before applying the PEROBVC method. For example, first find

${s}^{2}=\frac{1}{n}{\displaystyle {\sum}_{1}^{n}{s}_{i}^{2}}$ and then ${s}_{0.9999}^{2}=\frac{1}{f}{\chi}_{f;0.9999}^{2}\cdot {s}^{2}$ .

Those ${s}_{i}^{2}$ exceeding ${s}_{0.9999}^{2}$ are rejected, whereas this procedure is repeated once more with non-rejected results. ∆

The advantages of the PEROBVC method are:

1) Unbiased estimators for
${\sigma}_{1}^{2}$ and
${\sigma}_{2}^{2}$ , if B are close to *B*_{opt}, and

2) Minimal variances for ${\sigma}_{1}^{2}$ .

The disadvantages of the PEROBVC method are:

1) A sensitivity to the choice of the point B, (point B must be close to B_{opt}), which can result in negative estimates for either of the variances
${\sigma}_{1}^{2}$ or
${\sigma}_{2}^{2}$ , or for both, and

2) Sensitivity to the choice of the initial values for the variances, ${\sigma}_{1}^{2}$ and ${\sigma}_{2}^{2}$ , which, also, can result in negative estimates for one or both variances.

The stopping criterion for the iterative process: Let ${x}_{k}^{\text{T}}=\left[\begin{array}{cc}{\sigma}_{1,k}^{2}& {\sigma}_{2,k}^{2}\end{array}\right]$ be the vector of parameter estimates in the k-th iteration and d the difference vector of these estimates from the (k + 1)-th and k-th iterations, ${d}^{\text{T}}=\left[\begin{array}{cc}\left({\sigma}_{1,k+1}^{2}-{\sigma}_{1,k}^{2}\right)& \left({\sigma}_{2,k+1}^{2}-{\sigma}_{2,k}^{2}\right)\end{array}\right]$ , then the iterations should be stopped if

$\Vert d\Vert <{10}^{-q}\Vert {x}_{k}\Vert ,\text{\hspace{0.17em}}\text{\hspace{0.17em}}q\in \left\{5,6,7,8,9\right\}.$ (16)

4. Some Robust Estimation

In this case we shall use the method of Maximum Likelihood (ML).

The Robust Maximum Likelihood (ML) Estimation

The ML method is based on the assumption that in the domain X there exists only variances ${{s}^{\prime}}^{2}$ ( ${{s}^{\prime}}^{2}\equiv {s}^{2}$ ; with $E\left[{s}^{\text{2}}\right]={\sigma}^{\text{2}}$ ) of the basic distribution, whereas in the Y interval in addition to basic variances exist the variances of contaminating distributions. Therefore the point B in Figure 1 will be the censoring point.

The ML solution is obtained from the probability maximum for the event

$\left[{D}_{\left(1\right)}\wedge \cdots \wedge {D}_{\left({n}_{X}\right)}\right]\wedge \left({s}^{2}>B\right)$ ,

where ${D}_{\left(i\right)}={s}_{\left(i\right)}^{2}<{s}^{2}<{s}_{\left(i\right)}^{2}+\text{d}{s}^{2}$ , $i=1,\cdots ,{n}_{X}$ , one obtains the likelihood function (with proportionality constant k):

$L=k{\displaystyle {\prod}_{i=1}^{{n}_{X}}f\left({s}_{i}^{2}\right)}{\left[1-F\left(B\right)\right]}^{{n}_{Y}}$ .

Hence it follows

$\mathrm{ln}L\propto {\displaystyle {\sum}_{X}\mathrm{ln}f\left({s}_{i}^{2}\right)}+{n}_{Y}\mathrm{ln}\left[1-F\left(B\right)\right]$ .

From the necessary condition $\partial \mathrm{ln}L/\partial {\sigma}_{1}^{2}=0$ one obtains

$-{n}_{X}{\sigma}^{2}+{\displaystyle {\sum}_{X}{s}_{i}^{2}}+2B{n}_{Y}\frac{f\left({y}_{B}\right)}{1-F\left({y}_{B}\right)}=0$ .

Their solution is given as simple iterations

${\sigma}_{k+1}^{2}=\frac{1}{{n}_{X}}{\displaystyle {\sum}_{X}{s}_{i}^{2}}+2B\frac{{n}_{Y}}{{n}_{X}}\cdot \frac{f\left({y}_{B,k}\right)}{1-F\left({y}_{B,k}\right)}$ , (17)

where

${y}_{B,k}=\frac{f\cdot B}{{\sigma}_{k}^{2}}$ , ${F}_{B,k}={\displaystyle {\int}_{0}^{{y}_{B,k}}f\left(y\right)\text{d}y}$ , $\left(y=f\frac{{s}^{2}}{{\sigma}^{2}}={\chi}_{f}^{2},f\left(y\right)={k}_{f}^{-1}{y}^{f/2-1}{\text{e}}^{-y/2}\right)$ . (18)

The stopping criterion for the iterative process can be based on the parameter difference ${r}_{\sigma}=\left({\sigma}_{k+1}^{2}-{\sigma}_{k}^{2}\right)/{\sigma}_{k+1}^{2}$ between two iterations:

$\left|{r}_{\sigma}\right|<{10}^{-c}$ , $c\in \left\{4,5,6,7,8\right\}$ , (19)

5. Results and Discussion

Example 1. By using the functions F-RAN1 and F-RAN2 [25], one simulates 10,000, (n = 10,000), samples (series) of size 4, f = 3, with F-RAN1 5000 samples are simulated ( ${n}^{\prime}=5000$ ), ${\sigma}_{1}^{2}=1$ and the variance estimate (from these 5000 samples) ${s}_{1}^{2}=0.9997$ , with F-RAN2 also 5000 samples ( ${n}^{\u2033}=5000$ ), of the same size, with ${\sigma}_{2}^{2}=2$ and variance estimate (from these 5000 samples) ${s}_{2}^{2}=1.9886$ . Then a PERG mixed distribution (2), is formed ( $\epsilon =50\%$ ).

Since the variances
${\sigma}_{1}^{2}$ and
${\sigma}_{2}^{2}$ are known then, according to (6), one can also find, B optimal, obtaining B_{opt} = 2.772589.

In Table 1, for 4 partitions (X, Y), we present the PEROBVC estimates of the variances ${\sigma}_{1}^{2}$ and ${\sigma}_{2}^{2}$ , as well as the contamination percentage $\epsilon $ , ( $\epsilon ={n}^{\u2033}/n$ ).

From Table 1 it is easily noticed that the variance components, as well as the numbers of observations for both distributions, are well estimated, i.e. their values are very close to the real ones. It is also noticed that the values of the ML estimates for ${\sigma}^{2}$ are smaller than those of the PEROBVC ones for ${\sigma}_{1}^{2}$ (where it is ${\sigma}^{2}\equiv {\sigma}_{1}^{2}$ ). ▲

Example 2. In the 2D geodetic control network of object “TUZLA” [26] from the measurements of 328 (n = 328), horizontal angles, with a theodolite Wild T3, in 6 quick locks, in stable and unstable conditions, a mixed distribution of empirical variances (in ["^{2}]), with f = 5 d.f., is obtained (Table 2).

In Table 3, for 7 partitions (X, Y), we present the estimates of the variances

Table 1. Results of PEROBVC and ML estimations in simulated PERG mixed distribution of empirical variances s^{2} with
$n=10000$ ,
${n}^{\prime}=5000$ ,
${n}^{\u2033}=5000$ (
$\epsilon =50\%$ ),
${\sigma}_{1}^{2}=1$ ,
${\sigma}_{2}^{2}=2$ , f = 3 d.f., and.

Table 2. Ordered grouped increasing sequence of 328 (n = 328) empirical variances, with f = 5 d.f., of horizontal angle measurements with a theodolite Wild T3 in the geodetic 2D control network of object “TUZLA” in epochs 1979, 1981 and 1984.

Table 3. Results of PEROBVC and ML estimates of variances for PERG mixed distribution composed of 328 (n = 328), empirical variances of horizontal angles in the 2D control network of object “TUZLA”, in epochs 1979, 1981 and 1984; f = 5 d.f., and.

and, as well as the contamination percentage ().

According to an a priori analysis for the horizontal-angle measurements in stable conditions, following Činklović [27], one obtains ["^{2}], resulting in a very good agreement with 1.8109 ["^{2}] obtained by using PEROBVC method based on the measurements.

Besides, on the basis of the results presented in Table 3 one concludes about a very high stability of the estimator for the variance component.

In Example 2 the average value from the ML estimates for, (), is also smaller than the average value of the PEROBVC estimates for (1.6753 < 1.8109).

It should be noted that for the partition the numbers of observations and are well estimated, although B = 1.8 is significantly smaller than B_{opt} = 8.18. ▲

Note 3. In Examples 1-2 the ML estimates for are smaller than the corresponding PEROBVC estimates for, (), because higher values of empirical variances are in the end of the distribution so that they do not enter the ML estimate for. ∆

On the basis of applications of the PEROBVC method in the treatment of real geodetic measurements, Example 2, the present author concludes that the variance estimator of the basic distribution, has a good stability.

Generally, the PEROBVC method offers good estimates, especially for, which is illustrated by Examples 1-2.

6. Conclusions

On the basis of the obtained results in Examples 1-2 we can conclude the following:

1) On the basis of exact (expected) values from Example 1 the validity of the PEROBVC method in the variance estimation, as well as the estimates of numbers of observations, for both distributions in the PERG mixed distribution of empirical variances is confirmed. Here the variance estimates for both distributions, basic and contaminating ones, are correct; i.e. their values are exact.

2) On the basis of realistic measurements for the horizontal angles in the geodetic 2D control network from Examples 2 good (satisfactory) parameter estimates for both distributions are also confirmed.

3) In Example 2 the variance estimate for the basic distribution is confirmed through the result of the a priori analysis.

Acknowledgements

Dr ZORICA Cvetković from the Astronomical Observatory in Belgrade is the author of the FORTRAN programmes who performed the calculations in the examples. The OBJEKTIV GEO D.O.O. Company from Belgrade (Serbia) has financially supported the publishing of the article.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

[1] | Rao, C.R. and Kleffe, J. (1988) Estimation of Variance Components and Applications. North-Holland, Amsterdam. |

[2] | Krishnaiah, P.R., (1984) Analysis of Variance. North-Holland, Amsterdam. |

[3] |
Harville, D.A. (1977) Maximum Likelihood Approaches to Variance Component Estimation and to Related Problems. Journal of the American Statistical Association, 72, 320-338. https://doi.org/10.1080/01621459.1977.10480998 |

[4] |
Fellner, W.H. (1986) Robust Estimation of Variance Components. Technometrics, 28, 51-60. https://doi.org/10.1080/00401706.1986.10488097 |

[5] |
Robinson, D.L. (1987) Estimation and Use of Variance Components. Journal of the Royal Statistical Society. Series D (The Statistician), 36, 3-14. https://doi.org/10.2307/2988267 |

[6] |
Herbert, J.H. and Kott, Ph.S. (1988) Robust Variance Estimation in Linear Regression. Journal of Applied Statistics, 15, 341-345. https://doi.org/10.1080/02664768800000044 |

[7] |
Sidik, K. and Jonkman, J.N. (2006) Robust Variance Estimation for Random Effects Meta-Analysis. Computational Statistics & Data Analysis, 50, 3681-3701. https://doi.org/10.1016/j.csda.2005.07.019 |

[8] |
Hedges, L.V., Tipton, E. and Johnson, M.C. (2010) Robust Variance Estimation in Meta-Regression with Dependent Effect Size Estimates. Research Synthesis Methods, 1, 39-65. https://doi.org/10.1002/jrsm.5 |

[9] | Koller, M. (2013) Robust Estimation of Linear Mixed Models. PhD Thesis, ETH, Zürich. |

[10] |
Welsh, A.H. and Richardson, A.M. (1997) 13 Approaches to Robust Estimations of Mixed Models. In: Handbook of Statistics, Vol. 15, Elsevier Science, Amsterdam, 343-384. https://doi.org/10.1016/S0169-7161(97)15015-5 |

[11] | Amiri-Simkooei, A.R. (2007) Least-Squares Variance Component Estimation: Theory and GPS Applications. Ph.D. Thesis, Delft University of Technology, Delft. |

[12] | Bartlett, J. (2014) The Robust Sandwich Variance Estimator for Linear Regression (Using R). The Stats Geek, February 14, 2014. |

[13] |
Arango-Castillo, L. and Takahara, G. (2018) Robust Estimation of the Sample Mean Variance for Gaussian Processes with Long-Range Dependence. 2017 IEEE Global Conference on Signal and Information Processing, Montreal, 14-16 November 2017, 201-205. https://doi.org/10.1109/GlobalSIP.2017.8308632 |

[14] |
Chen, Y. and Jackson, D. (1995) Robust Estimation of Mean and Variance in Fisheries. Transactions of the American Fisheries Society, 124, 401-412. https://doi.org/10.1577/1548-8659(1995)124<0401:REOMAV>2.3.CO;2 |

[15] |
Grossi, L. and Laurini, F. (2011) Robust Estimation of Efficient Mean-Variance Frontiers. Advances in Data Analysis and Classification, 5, 3-22. https://doi.org/10.1007/s11634-010-0082-3 |

[16] |
Li, Y., Li, J., Qi, J. and Chen, L. (2019) Robust Cubature Kalman Filter for Dynamic State Estimation of Synchronous Machines under Unknown Measurement Noise Statistics. IEEE Access, 7, 29139-29148. https://doi.org/10.1109/ACCESS.2019.2900228 |

[17] |
Singh, A. and Nocerino, J. (2002) Robust Estimation of Mean and Variance Using Environmental Data Sets with Below Detection Limit Observations. Chemometrics and Intelligent Laboratory Systems, 60, 69-86. https://doi.org/10.1016/S0169-7439(01)00186-1 |

[18] |
Vilà-Valls, J., Wei, Q., Closas, P. and Fernández-Prades, C. (2014) Robust Gaussian Sum Filtering with Unknown Noise statistics: Application to Target Tracking. 2014 IEEE Workshop on Statistical Signal Processing, Gold Coast, 29 June-2 July 2014, 416-419. https://doi.org/10.1109/SSP.2014.6884664 |

[19] |
Wang, Ch., Wang, Y., Lin, Z., Yuille, A.L. and Gao, W. (2014) Robust Estimation of 3D Human Poses from a Single Image. CBMM Memo No. 013. https://cbmm.mit.edu/sites/default/files/publications/CBMM-Memo-013.pdf |

[20] |
Zhou, W.-Z., Bose, K., Fan, J. and Liu, H. (2018) A New Perspective on Robust M-Estimation: Finite Sample Theory and Applications to Ddependence-Adjusted Multiple Testing. The Annals of Statistics, 46, 1904-1931. https://doi.org/10.1214/17-AOS1606 |

[21] |
Zhao, J. and Mili, L. (2018) A Framework for Robust Hybrid State Estimation with Unknown Measurement Noise Statistics. IEEE Transactions on Industrial Informatics, 14, 1866-1875. https://doi.org/10.1109/TII.2017.2764800 |

[22] |
Perovic, G. (2019) Robust Estimation of the Normal-Distribution Parameters by Use of Structural Partitioning-Perobls D Method. American Journal of Computational Mathematics, 9, 302-316. https://doi.org/10.4236/ajcm.2019.94022 |

[23] | Perovic, G. (2016) Variance Components Analysis in GPS Measurements. Imperial Journal of Interdisciplinary Research, 2, 357-367. |

[24] | Perovic, G. (2005) Least Squares. Author, Belgrade. |

[25] | William, H.P., et al. (1986) Numerical Recipes. Cambridge Univ. Press, Cambridge. |

[26] | “Tuzla” Object (1985) Elaborates of Geodetic Measurements in Epochs 1979, 1981 and 1984 for the Purposes of Ground—Deformation Determining. Salt Mine, Tuzla. |

[27] | Cinklovic, N. (1978) Analysis and a Priori Accuracy Estimate for Precise Geodetic Measurements, Monograph. Faculty of Civil Engineering, Belgrade. |

Copyright © 2020 by authors and Scientific Research Publishing Inc.

This work and the related PDF file are licensed under a Creative Commons Attribution 4.0 International License.