# A Novel Bayesian Method for Calculating Circular Error Probability with Systematic-Biased Prior Information

Authors: Bowen Liu; Xiaojun Duan; Liang Yan
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

Circular Error Probability (CEP) is defined as the radius of a circle where the probability of an impact point being inside is 50%, which is also widely used as a measure of the guidance weapon systems’ precision. In order to achieve a fusion of various test information, Bayesian methods and improved Bayesian methods have been extensively studied in calculating the CEP. Nevertheless, these methods could fail when there exists unknown systematic bias in the prior information. Therefore, a novel method called Bayesian estimation based on representative points (BERP) with an optimization procedure for determining the optimal number of representative points is proposed in this paper. Explicit theoretical analyses demonstrate that the BERP outperforms the classical Bayesian methods when fusing the slightly biased prior information and also give the bound of the systematic bias for stopping using the heavily biased prior information. Moreover, when the systematic bias is within the bound, simulation results indicate that our method is credible and outperforms the classical Bayesian method in calculating the CEP of guidance weapon systems.


---

## Body

## 1. Introduction

Performance evaluation of complex equipment, such as the guidance weapon systems, is very important before application. In the assessment of impact accuracy of the guidance weapon systems, CEP is the most commonly used as a measure, which can integrate the precision and dispersion to assess the impact accuracy [1–3]. Usually, the tighter the pattern of impact point errors, the smaller the CEP we can get; i.e., there is higher impact accuracy of guidance weapon systems. Let (X,Z) be the impact point errors of the projectile, where X and Z are downrange and cross-range misses, respectively; then the CEP can be defined as follows [4, 5]:(1)PX2+Z2≤CEP2=∬x2+z2≤CEP2fx,zdxdz=0.5It is shown in the above equation that the CEP can be calculated by numerical integrations when giving the probability density function f(x,z). In practice, (X,Z) are usually assumed to follow the bivariate normal distribution, i.e., (X,Z)~N(μ,Σ),where μ and Σ are the mean and the covariance, respectively. Then the problem of calculating CEP rests on estimating μ and Σ, whose accuracy will determine the precision of the CEP calculation.In the performance evaluation of guidance weapon systems, the impact points used to calculate the CEP are collected from various tests. The realistic tests of guidance weapon systems are usually extremely expensive and time-consuming, so generally the sample size of impact points is very small. Moreover, the calculation of CEP based on the small sized data would lead to unreliable results, so it is reasonable to introduce the prior information to refine the results. Therefore, Bayesian estimation is widely used to achieve the fusion of the prior information and the realistic test information [6, 7] and hence to increase the reliability and the accuracy of the estimation of (μ,Σ), as well as the CEP. In fact, the most important prior information is provided by the substitute tests in the development processes. The performance evaluation is a sequential process, so the data collected from substitute tests are usually regarded as prior samples. In order to use the prior information more reasonably, several strategies are introduced in bringing the credibility of the prior information into Bayesian estimation, such as data compatibility tests, the information divergence, and the theory of fuzzy operators [8, 9]. In the CEP calculation, Huang comes up with a measure of credibility from physics resources of data [10], which has a good estimation accuracy if we have a clear understanding of the physical background. Duan suggested a method with all prior information normalized into one test sample [11], which can reduce the deviation when prior information got distorted, but more theoretical discussion is needed for this method when considering the information loss and fusion efficiency of the normalization.However, the substitute tests for the guidance weapon systems may be systematically biased compared with realistic tests. Since the pattern of systematic bias is unknown, we cannot estimate the bias but can only give a rough range of it. Simulation results demonstrate that the systematically biased prior information will cause serious impact on the mean estimationμ when applying classical Bayesian estimation directly. The estimation of the parameter Σ is slightly affected by the unknown systematic bias because of the same guidance system of these tests. The improved methods considering the normalization of prior information could reduce the estimate bias of μ, but this will cause much loss of prior information and inevitably leads to an unreliable calculation of CEP. One possible way to solve the problem is by choosing appropriate samples, such as representative points to generate new prior information.Representative points (RPs), also known as Principal Points [12], are a group of points that could represent a distribution with the least mean squared error (MSE). RPs’ theory is brought up in 1990s and now widely applied in clustering analyses [13], statistical simulations [14], image processing [15], and so on. In this paper, we resample samples from the systematically biased data to get the RPs which are regarded as the new prior information, expecting to reduce the estimate bias of μ substantially. This new estimation method is called Bayesian estimation based on representative points (BERP). Meanwhile, we propose an optimization procedure which balances the effects of estimate bias and information loss to determine the optimal number of RPs. In addition, two theorems are proposed to prove that the estimate bias and MSE of μ with RPs are smaller than those with the raw data. Furthermore, we also analyze the bound of the systematic bias for stopping using the heavily biased prior information. Within the bound, both the simulations and authentic experiments show that our BERP outperforms the classical Bayesian methods in estimating the parameters. In the performance evaluation of guidance weapon systems, via BERP it is better to choose RPs from the raw prior samples as new prior information when calculating the CEP. On the whole, by using the theory of RPs, our works enrich the Bayesian methods especially on the prior information fusion patterns theoretically and provide a possible way to solve the engineering problem in performance evaluation of guidance weapon systems.The rest of this paper is organized as follows. Section2 introduces the classical Bayesian estimation applied in CEP calculation. Section 3 describes some notations and preliminaries of RPs, the process of BERP, and the optimization procedure for determining the optimal number of RPs. In Section 4, we propose two theorems to compare the estimation performances by BERP and classical Bayesian estimation and analyze the bound of systematic bias for stopping using the biased prior information. Section 5 shows three numerical experiments about our new method. The conclusion is given in Section 6.


## 2. Bayesian Estimation in CEP Calculation

### 2.1. Calculation of CEP

Suppose that(X,Z) follows the bivariate normal distribution N(μ,Σ), where(2)μ=μx,μzT,Σ=σx2ρσxσzρσxσzσz2The parameters σx, σz are standard deviations of impact point errors for downrange direction and cross-range direction, respectively; μx, μz are means of impact point errors for each direction; ρ (0≤ρ<1) is the correlation coefficient of X and Z. In most cases, X and Z are independent of each other. Nevertheless, even if ρ≠0, we can use the orthogonal transformation to achieve the decorrelation of X,Z. So we assume that ρ=0 in the rest of the paper. Under the assumption that ρ=0, the CEP R satisfies the equation(3)∬x2+z2⩽R212πσxσzexp⁡-x-μx22σx2-z-μz22σz2dxdz=0.5After estimating the parametersμx, μz, σx, and σz, the CEP R can be calculated by numerical integrations. More details about calculating CEP are given in [16].


### 2.2. Classical Bayesian Estimation

As described in the Introduction, the performance evaluation of guidance weapon systems is a sequential process; the distribution parameters would change when fusing test data from different stage. The realistic tests are conducted to refine the previous evaluation results based on the substitute tests. Therefore, the data collected from substitute tests are usually regarded as prior samples.Because the procedures of estimating the parameters(μx,σx) and (μz,σz) are similar, we take the downrange direction of impact point errors as an example to introduce the classical Bayesian estimation, where X follows the normal distribution N(μx,σx2). For convenience, we drop the subscripts of (μx,σx2) as (μ,σ2) and let π1(μ,σ2) be the joint prior distribution of (μ,σ2). In Bayesian theory, the conjugate prior distributions of μ and σ2 are normal distribution N(μ1,σ2/k1) and inverse Gamma distribution IGaα1,β1, respectively. As for the distribution parameters μ1, k1, α1, and β1, they are determined by (7). The probability density function of IGaα1,β1 is(4)fσ2=β1α1Γα11σ2α1+1e-β1/σ2where Γ(·) means the Gamma function. Moreover, the joint prior distribution π1(μ,σ2) is the normal-inverse Gamma distribution, so we have(5)π1μ,σ2=π1μ∣σ2·π1σ2∝Nμ1,σ2k1·IGaα1,β1Suppose the prior samples for downrange direction are X(1)=x1(1),x2(1)⋯,xn1(1), n1 is the size of the samples; let(6)X¯1=1n1∑i=1n1xi1S12=1n1∑i=1n1xi1-X¯12Then the estimates of the parameters of the joint prior distribution are(7)α1=12∑i=1n1xi1-X¯12=n1S122β1=12n1-1μ1=X¯1k1=n1Similarly, when the realistic test samples X(2)=x1(2),x2(2),⋯,xn2(2) are obtained, n2 is the size of the samples; let(8)X¯2=1n2∑i=1n2xi2S22=1n2∑i=1n2xi2-X¯22Because of the property of conjugate prior distribution, the posterior distribution of (μ,σ2) is also a normal-inverse Gamma distribution:(9)π2μ,σ2∣X¯2,S22∝Nμ2,σ2k2·IGaα2,β2where μ2, k2, α2, and β2 are the parameters of the normal-inverse Gamma distribution, and the estimates of the parameters are(10)α2=α1+n2S222+n1n2X¯2-μ122n1+n2β2=β1+n22μ2=n1μ1+n2X¯2n1+n2k2=n1+n2So the estimates of (μ,σ2) by classical Bayesian estimation are(11)μ^Bayes=μ2σ^Bayes2=α2β2-1As for the cross-range direction of impact point errors, the estimates ofμz and σz can also be calculated by (11). After obtaining the estimates of μx, μz, σx, and σz, we can calculate the CEP of impact points by (3).


## 2.1. Calculation of CEP

Suppose that(X,Z) follows the bivariate normal distribution N(μ,Σ), where(2)μ=μx,μzT,Σ=σx2ρσxσzρσxσzσz2The parameters σx, σz are standard deviations of impact point errors for downrange direction and cross-range direction, respectively; μx, μz are means of impact point errors for each direction; ρ (0≤ρ<1) is the correlation coefficient of X and Z. In most cases, X and Z are independent of each other. Nevertheless, even if ρ≠0, we can use the orthogonal transformation to achieve the decorrelation of X,Z. So we assume that ρ=0 in the rest of the paper. Under the assumption that ρ=0, the CEP R satisfies the equation(3)∬x2+z2⩽R212πσxσzexp⁡-x-μx22σx2-z-μz22σz2dxdz=0.5After estimating the parametersμx, μz, σx, and σz, the CEP R can be calculated by numerical integrations. More details about calculating CEP are given in [16].


## 2.2. Classical Bayesian Estimation

As described in the Introduction, the performance evaluation of guidance weapon systems is a sequential process; the distribution parameters would change when fusing test data from different stage. The realistic tests are conducted to refine the previous evaluation results based on the substitute tests. Therefore, the data collected from substitute tests are usually regarded as prior samples.Because the procedures of estimating the parameters(μx,σx) and (μz,σz) are similar, we take the downrange direction of impact point errors as an example to introduce the classical Bayesian estimation, where X follows the normal distribution N(μx,σx2). For convenience, we drop the subscripts of (μx,σx2) as (μ,σ2) and let π1(μ,σ2) be the joint prior distribution of (μ,σ2). In Bayesian theory, the conjugate prior distributions of μ and σ2 are normal distribution N(μ1,σ2/k1) and inverse Gamma distribution IGaα1,β1, respectively. As for the distribution parameters μ1, k1, α1, and β1, they are determined by (7). The probability density function of IGaα1,β1 is(4)fσ2=β1α1Γα11σ2α1+1e-β1/σ2where Γ(·) means the Gamma function. Moreover, the joint prior distribution π1(μ,σ2) is the normal-inverse Gamma distribution, so we have(5)π1μ,σ2=π1μ∣σ2·π1σ2∝Nμ1,σ2k1·IGaα1,β1Suppose the prior samples for downrange direction are X(1)=x1(1),x2(1)⋯,xn1(1), n1 is the size of the samples; let(6)X¯1=1n1∑i=1n1xi1S12=1n1∑i=1n1xi1-X¯12Then the estimates of the parameters of the joint prior distribution are(7)α1=12∑i=1n1xi1-X¯12=n1S122β1=12n1-1μ1=X¯1k1=n1Similarly, when the realistic test samples X(2)=x1(2),x2(2),⋯,xn2(2) are obtained, n2 is the size of the samples; let(8)X¯2=1n2∑i=1n2xi2S22=1n2∑i=1n2xi2-X¯22Because of the property of conjugate prior distribution, the posterior distribution of (μ,σ2) is also a normal-inverse Gamma distribution:(9)π2μ,σ2∣X¯2,S22∝Nμ2,σ2k2·IGaα2,β2where μ2, k2, α2, and β2 are the parameters of the normal-inverse Gamma distribution, and the estimates of the parameters are(10)α2=α1+n2S222+n1n2X¯2-μ122n1+n2β2=β1+n22μ2=n1μ1+n2X¯2n1+n2k2=n1+n2So the estimates of (μ,σ2) by classical Bayesian estimation are(11)μ^Bayes=μ2σ^Bayes2=α2β2-1As for the cross-range direction of impact point errors, the estimates ofμz and σz can also be calculated by (11). After obtaining the estimates of μx, μz, σx, and σz, we can calculate the CEP of impact points by (3).


## 3. Bayesian Estimation Based on Representative Points

In this section, we propose the novel method Bayesian estimation based on representative points and the procedure of determining the optimal number of RPs. RPs can not only optimally represent the distribution of prior information in terms of MSE principle, but also have smaller sample size compared with raw prior samples. Therefore, RPs can retain the useful information of prior samples and reduce the estimate bias ofμx, μz. If we search RPs as new prior information, we may get more accurate and reliable CEP of guidance weapon systems.


### 3.1. Methods for Searching Representative Points

In this subsection, we will give a brief introduction of RPs and methods for searching RPs. Assume thatX∈Rp is a p dimensional random vector, and the probability density function of X is f(x). Define the mean squared error for a set of points r1,r2,⋯,rk⊆Rp of the random vector X as follows:(12)MSEr1,r2,⋯,rk=Emini=1,⋯,kx-ri2=∫Rpmini=1,⋯,kx-ri2fxdxwhere · stands for L2 -norm. The vectors {r1,⋯rk}∈Rp are called k representative points of a random vector X if(13)MSEr1,r2,⋯,rk≤MSEx1,x2,⋯,xkfor all sets x1,x2,⋯,xk⊆Rp.From the definition above, it is obvious that whenk=1, the single RP is the mean of X. Searching the RPs equals doing the optimal grouping [17]; it is difficult to derive the concrete RPs theoretically even if the number of RPs is given. Flury has proved that there is no theoretical derivation of RPs when k>2 [18]. So some approximation algorithms have been proposed to search RPs including k-means methods [19], parametric k-means methods [20], and nonparametric methods [21]. The k-means methods are searching the clustering centers as the RPs. Moreover, the parametric k-means methods resample large samples from a specific distribution whose parameters are estimated by maximum likelihood, then searching the RPs from the resampled samples by k-means methods. The main idea of nonparametric methods is to build the empirical distribution function of X and resample large samples from this empirical function, after which the RPs are chosen from the resampled samples. In most cases, nonparametric methods have better performance to represent a distribution in terms of MSE than the k-means methods and parametric k-means methods. The main steps of the nonparametric method introduced in [21] are shown as follows:(i)
Step1. For original samples X(0)={x1,x2,⋯,xn}, use k-means algorithm to obtain k points P(1)={P1(1),P2(1),⋯,Pk(1)} from X(0) as an initial solution.(ii)
Step2. Use the kernel estimation method to estimate the density function of original samples, denoted as f^(x).(iii)
Step3.Use the Randomized Likelihood Sampling method to generate N(N≫n) samples Z={z1,z2,⋯,zN} from the density function f^(x) as the training data.(iv)
Step4. Based on the N training data Z and the k starting points P(1), use k-means algorithm to obtain k RPs P(2)={P1(2),P2(2)⋯,Pk(2)}.Following the four steps above, we can get the RPs from prior samples. After that, we can use the RPs as new prior information to estimate the population parameters. More details about BERP will be introduced in the next subsection.


### 3.2. The Procedure of BERP

In this subsection, we will introduce the procedure of BERP which is similar to the classical Bayesian estimation in Section2.2. We also take the downrange direction of impact point errors as an example, suppose the original prior samples are X(1)=x1(1),x2(1),⋯,xn1(1). Choose nr RPs from X(1) by the nonparametric method introduced in Section 3.1 and denote X(rp)=x1(rp),x2(rp),⋯,xnr(rp), and let(14)X¯rp=1nr∑i=1nrxirpSrp2=1nr∑i=1nrxirp-X¯rp2Similar to Section 2.2, let αrp, βrp, μrp, and krp be the parameters of the joint prior distribution. By searching RPs as new prior information, the estimates of these parameters are(15)αrp=nrSrp22βrp=nr-12μrp=X¯rpkrp=nrMoreover, let μ2r, k2r, α2r, and β2r be the parameters of posterior distribution. So the estimates of these parameters are(16)α2r=αrp+n2S222+nrn2X¯2-μrp22nr+n2β2r=βrp+n22μ2r=nrμrp+n2X¯2nr+n2k2r=nr+n2where X¯(2) and S22 are calculated by (8). So the estimates of (μ,σ2) by BERP are(17)μ^BERP=μ2rσ^BERP2=α2rβ2r-1If we know the exact number of RPs, it is easy to estimate the parametersμx, μz, σx, and σz by (17). However, it is difficult to determine the optimal number of RPs because there are no theoretical methods about this. Therefore, we propose an optimization procedure to determine the optimal number of RPs when there exists unknown systematic bias in prior samples.


### 3.3. Optimal Number of Representative Points

The approach to determining the optimal number of RPs varies with the background of practical problem. On the one hand, the RPs are closer to the original prior samples as the number grows. Therefore, the more the RPs chosen, the larger bias they may bring to the estimate ofμ, which will also be validated in Theorem 1, Section 4. On the other hand, similar to the methods of normalizing prior information, choosing small number of RPs will cause great information loss. Therefore, we consider two factors when determining the optimal number of RPs: the estimate bias and the information loss. We still take the downrange direction of the impact point errors as an example to describe the optimization procedure. Bnr and Lnr stand for the estimate bias and the information loss, respectively, where nr stands for the number of RPs. The objective function Fnr is to achieve a balance between Bnr and Lnr. So we have(18)Fnr=Bnr+LnrIn fact, it is hard to quantify the estimate biasBnr without the true values of the population parameters. But there is a principle in the performance evaluation of guidance weapon systems that every realistic test sample should be used. So we can use the mean of realistic test samples to approximate the true value of parameter μ, denoted as μtrue. The estimate of μ by BERP is denoted as μBERP. So the approximated Bnr is(19)Bnr=μBERP-μtrueAs for Lnr, the information loss could be estimated by Cox (1957) [17]:(20)Lnr=EX-ξnrX2σ2=∑i=1nrσi2σ2where σ2 is the variance of all prior samples; σi2 is the variance of samples in ith class which is classified in Step4 for searching RPs. As the number of RPs increases, the estimate bias would increase and the information loss would decrease. The optimal number of RPs is determined by the minimal value of the objective function(21)minnrμBERP-μtrue+∑i=1nrσi2σ2Algorithm 1 describes the optimization procedure of determining the optimal number of RPs.Algorithm 1:The optimization procedure.
DefineN: Predetermined maximum number of RPs nr: The current number of RPs nr∗: Optimal number of RPs Bnr: Estimate bias Lnr: Information loss for  1⩽nr⩽N  do Choosenr RPs from the prior samples and calculate Bnr and Lnr Updatenr=nr+1 end for nr∗=arg⁡min⁡(Bnr+Lnr) return  nr∗


## 3.1. Methods for Searching Representative Points

In this subsection, we will give a brief introduction of RPs and methods for searching RPs. Assume thatX∈Rp is a p dimensional random vector, and the probability density function of X is f(x). Define the mean squared error for a set of points r1,r2,⋯,rk⊆Rp of the random vector X as follows:(12)MSEr1,r2,⋯,rk=Emini=1,⋯,kx-ri2=∫Rpmini=1,⋯,kx-ri2fxdxwhere · stands for L2 -norm. The vectors {r1,⋯rk}∈Rp are called k representative points of a random vector X if(13)MSEr1,r2,⋯,rk≤MSEx1,x2,⋯,xkfor all sets x1,x2,⋯,xk⊆Rp.From the definition above, it is obvious that whenk=1, the single RP is the mean of X. Searching the RPs equals doing the optimal grouping [17]; it is difficult to derive the concrete RPs theoretically even if the number of RPs is given. Flury has proved that there is no theoretical derivation of RPs when k>2 [18]. So some approximation algorithms have been proposed to search RPs including k-means methods [19], parametric k-means methods [20], and nonparametric methods [21]. The k-means methods are searching the clustering centers as the RPs. Moreover, the parametric k-means methods resample large samples from a specific distribution whose parameters are estimated by maximum likelihood, then searching the RPs from the resampled samples by k-means methods. The main idea of nonparametric methods is to build the empirical distribution function of X and resample large samples from this empirical function, after which the RPs are chosen from the resampled samples. In most cases, nonparametric methods have better performance to represent a distribution in terms of MSE than the k-means methods and parametric k-means methods. The main steps of the nonparametric method introduced in [21] are shown as follows:(i)
Step1. For original samples X(0)={x1,x2,⋯,xn}, use k-means algorithm to obtain k points P(1)={P1(1),P2(1),⋯,Pk(1)} from X(0) as an initial solution.(ii)
Step2. Use the kernel estimation method to estimate the density function of original samples, denoted as f^(x).(iii)
Step3.Use the Randomized Likelihood Sampling method to generate N(N≫n) samples Z={z1,z2,⋯,zN} from the density function f^(x) as the training data.(iv)
Step4. Based on the N training data Z and the k starting points P(1), use k-means algorithm to obtain k RPs P(2)={P1(2),P2(2)⋯,Pk(2)}.Following the four steps above, we can get the RPs from prior samples. After that, we can use the RPs as new prior information to estimate the population parameters. More details about BERP will be introduced in the next subsection.


## 3.2. The Procedure of BERP

In this subsection, we will introduce the procedure of BERP which is similar to the classical Bayesian estimation in Section2.2. We also take the downrange direction of impact point errors as an example, suppose the original prior samples are X(1)=x1(1),x2(1),⋯,xn1(1). Choose nr RPs from X(1) by the nonparametric method introduced in Section 3.1 and denote X(rp)=x1(rp),x2(rp),⋯,xnr(rp), and let(14)X¯rp=1nr∑i=1nrxirpSrp2=1nr∑i=1nrxirp-X¯rp2Similar to Section 2.2, let αrp, βrp, μrp, and krp be the parameters of the joint prior distribution. By searching RPs as new prior information, the estimates of these parameters are(15)αrp=nrSrp22βrp=nr-12μrp=X¯rpkrp=nrMoreover, let μ2r, k2r, α2r, and β2r be the parameters of posterior distribution. So the estimates of these parameters are(16)α2r=αrp+n2S222+nrn2X¯2-μrp22nr+n2β2r=βrp+n22μ2r=nrμrp+n2X¯2nr+n2k2r=nr+n2where X¯(2) and S22 are calculated by (8). So the estimates of (μ,σ2) by BERP are(17)μ^BERP=μ2rσ^BERP2=α2rβ2r-1If we know the exact number of RPs, it is easy to estimate the parametersμx, μz, σx, and σz by (17). However, it is difficult to determine the optimal number of RPs because there are no theoretical methods about this. Therefore, we propose an optimization procedure to determine the optimal number of RPs when there exists unknown systematic bias in prior samples.


## 3.3. Optimal Number of Representative Points

The approach to determining the optimal number of RPs varies with the background of practical problem. On the one hand, the RPs are closer to the original prior samples as the number grows. Therefore, the more the RPs chosen, the larger bias they may bring to the estimate ofμ, which will also be validated in Theorem 1, Section 4. On the other hand, similar to the methods of normalizing prior information, choosing small number of RPs will cause great information loss. Therefore, we consider two factors when determining the optimal number of RPs: the estimate bias and the information loss. We still take the downrange direction of the impact point errors as an example to describe the optimization procedure. Bnr and Lnr stand for the estimate bias and the information loss, respectively, where nr stands for the number of RPs. The objective function Fnr is to achieve a balance between Bnr and Lnr. So we have(18)Fnr=Bnr+LnrIn fact, it is hard to quantify the estimate biasBnr without the true values of the population parameters. But there is a principle in the performance evaluation of guidance weapon systems that every realistic test sample should be used. So we can use the mean of realistic test samples to approximate the true value of parameter μ, denoted as μtrue. The estimate of μ by BERP is denoted as μBERP. So the approximated Bnr is(19)Bnr=μBERP-μtrueAs for Lnr, the information loss could be estimated by Cox (1957) [17]:(20)Lnr=EX-ξnrX2σ2=∑i=1nrσi2σ2where σ2 is the variance of all prior samples; σi2 is the variance of samples in ith class which is classified in Step4 for searching RPs. As the number of RPs increases, the estimate bias would increase and the information loss would decrease. The optimal number of RPs is determined by the minimal value of the objective function(21)minnrμBERP-μtrue+∑i=1nrσi2σ2Algorithm 1 describes the optimization procedure of determining the optimal number of RPs.Algorithm 1:The optimization procedure.
DefineN: Predetermined maximum number of RPs nr: The current number of RPs nr∗: Optimal number of RPs Bnr: Estimate bias Lnr: Information loss for  1⩽nr⩽N  do Choosenr RPs from the prior samples and calculate Bnr and Lnr Updatenr=nr+1 end for nr∗=arg⁡min⁡(Bnr+Lnr) return  nr∗


## 4. Theoretical Analysis of BERP

Estimate bias and MSE are the common measures to evaluate the quality of an estimator. So we will use them to analyse the theoretical performance of BERP and classical Bayesian estimation in this section. Suppose the posterior distribution of parameterθ is f(θ∣X), θ¯ is the posterior expectation of θ, and θtrue is the true value of θ. The estimate bias and MSE of the estimator θ^ are(22)Biasθ^,θtrue∣X=θ¯-θtrue(23)MSEθ^∣X=Eθ^∣Xθ^-θtrue2=Eθ^∣Xθ^-θ¯+θ¯-θtrue2=Varθ^∣X+θ¯-θtrue2=Varθ^∣X+Biasθ^,θtrue∣X2where Var(θ^∣X) is the variance of θ^.Let realistic test samples follow the normal distributionN(μ,σ2), and the prior samples follow the normal distribution N(μ+η,σ2), where η is the systematic bias. Let n1, n2 be the sample sizes of prior samples and realistic test samples, respectively. Suppose the posterior estimate of μ by classical Bayesian estimation is μ^Bayes, μ¯Bayes is its posterior expectation. The posterior estimate of μ by BERP is μ^BERP, μ¯BERP is its posterior expectation, and nr is the number of RPs. Let μtrue be the true value of parameter μ. There are two theorems to compare the estimate bias and MSE of the two estimators μ^Bayes and μ^BERP.Theorem 1.
In the case that there exists systematic bias in prior samples, one has(24)Biasμ^BERP,μtrue∣X≤Biasμ^Bayes,μtrue∣XProof.
From (11) and (17), the estimate bias of μ^BERP and μ^Bayes is(25)Biasμ^BERP,μtrue∣X=μ¯BERP-μtrue=nrηnr+n2(26)Biasμ^Bayes,μtrueX=μ¯Bayes-μtrue=n1ηn1+n2 The estimate bias of μBERP increases with the number of RPs. Since the size of prior samples is much larger than the size of RPs, we have n1≫nr. Therefore there is (27)nrnr+n2<n1n1+n2 So (28)Biasμ^BERP,μtrue∣X≤Biasμ^Bayes,μtrue∣XTheorem 2.
In the case thatη2⩾σ2/n2, one has (29)MSEμ^BERP∣X≤MSEμ^Bayes∣XProof.
The MSE ofμ^Bayes and μ^BERP are(30)MSEμ^Bayes∣X=Varμ^Bayes∣X+Biasμ^Bayes,μtrue∣X2(31)MSEμ^BERP∣X=Varμ^BERP∣X+Biasμ^BERP,μtrue∣X2where (32)Varμ^Bayes∣X=α2k2β2-1=k2·n1·S12+n2·S22+n1·n2·X¯2-X¯12k22·k2-3Varμ^BERP∣X=α2rk2rβ2r-1=k2r·nr·Srp2+n2·S22+nr·n2·X¯2-X¯rp2k2r2·k2r-3
Because prior samples and realistic test samples both follow the normal distributionN(μ,σ2) and there exists systematic bias η in prior samples, we have the approximated results: (33)S12≈S22≈Srp2≈σ2X¯2-X¯12≈X¯2-X¯rp2≈η2
Letλ1=n1+n2 and λr=nr+n2, and Δ1 is determined as follows:(34)Δ1=λ12λr2λ1-3λr-3n1λr+nrλ1-n1nrλ1+λr+n2-3+n23-3n22when MSE(μ^BERP∣X)≤MSE(μ^Bayes∣X), if and only if η2⩾σ2/n2·Δ1. Because n1≫n2, n1≫nr, we can get Δ1<1. To sum up, when η2≥σ2/n2, there is (35)MSEμ^BERP∣X≤MSEμ^Bayes∣XThe estimate bias ofμ^BERP is smaller than that of μ^Bayes when there exists unknown systematic bias in prior samples. Moreover, the MSE of μ^BERP is also smaller than that of μ^Bayes when η2≥σ2/n2. In most cases, σ2/n2 is smaller than η2. Therefore, it can be concluded that BERP has better accuracy for estimating the parameter μ of normal distribution than classical Bayesian estimation when there exists unknown systematic bias in prior samples.However, when the systematic bias is beyond a certain bound, it may be better to stop fusing the biased prior samples even if they may provide some useful information. Supposeμ^MLE is estimated by maximum likelihood estimation without using the prior information, it is obviously an unbiased estimate. The MSE of μ^MLE is(36)MSEμ^MLE∣X=Varμ^MLE∣X=σ2n2whenMSE(μ^MLE∣X)≤MSE(μ^BERP∣X), if and only if η2⩾σ2/n2·Δ2, where(37)Δ2=nr-3nr+n22nrn2+nr2n2+nr-3Moreover, compared with (34), there is Δ1<Δ2, which means that when η2⩾σ2/n2·Δ2, MSE(μ^MLE∣X) is also smaller than MSE(μ^Bayes∣X). Therefore, if η2 is larger than σ2/n2·Δ2, we should stop using the prior information; σ2/n2·Δ2 is the bound of the systematic bias.


## 5. Numerical Experiments

In this section, three numerical experiments are provided to validate that BERP can help to get more accurate and reliable calculation of CEP when the systematic bias is within the bound. The first one is to show the performance of optimization procedure when determining the optimal number of RPs; the second one is to compare the estimation accuracy of BERP with that of classical Bayesian estimation; the third one is to compare the calculation of CEP based on BERP with classical Bayesian estimation.Example 1.
In this example, we will analyze the optimization procedure for determining the optimal number of RPs. Letn1=200 be the sample size of the prior samples and n2=5 be the sample size of realistic test samples. The prior samples and the realistic test samples are generated from normal distributions N(μ+η,σ2) and N(μ,σ2), respectively, where μ=0,η2=10,σ2=50. The RPs are searched by the nonparametric method. In order to reduce the influence of random factors on simulation result, resample samples 100 times under the same circumstance and get the average results. Figure 1 shows the variation of the information loss, the estimate bias, and the objective function for different numbers of RPs, where the estimate bias and information loss are both calculated after normalization.
As displayed in Figure1, the estimate bias approximated by (19) is very close to the theoretical one, which means that using Bnr to quantify the estimate bias is reasonable. In addition, when the number of RPs increases, the information loss decreases while the estimate bias increases, so the objective function can balance the two factors well when determining the optimal number of RPs. From Figure 1(b), when the number of RPs is 9, the objective function reaches the minimal value, so the optimal number of RPs is 9 in this example.Figure 1
Performance for different representative points.
(a) (b)Example 2.
In this example, we will compare the estimation performance of the two methods. Letn1=200 be the sample size of the prior samples, n2=5 be the sample size of the realistic test samples, and N=20 be the predetermined maximum number of RPs. The prior samples and realistic test samples are generated from normal distributions N(μ+η,σ2) and N(μ,σ2), respectively, where μ=0 and η is the systematic bias in prior samples. We set different values for η in the simulations to investigate the estimation performance when estimating the parameter μ. Similar to Example 1, we resample samples 100 times under the same circumstance and get the average results in order to reduce the influence of random factors. Table 1 shows the simulation results.
From Table1, we can summarize two conclusions. (1) Within the bound of the systematic bias shown in Section 4, if there exists slight systematic bias in prior samples, μ^BERP is much more closer to the true value of the parameter μ than μ^Bayes, and MSE(μ^BERP) is also smaller than MSE(μ^Bayes) in most cases. Therefore, BERP has higher accuracy to estimate the parameter μ than classical Bayesian estimation when there exists slight systematic bias in prior information. (2) If there is no systematic bias in prior samples, μ^BERP is very close to μ^Bayes, and MSE(μ^BERP) is a little larger than MSE(μ^Bayes). There is no obvious difference between the estimation accuracies of the two methods. Moreover, without the systematic bias in prior samples, the optimal number of RPs is close to the predetermined maximum number 20. In this case, BERP is degenerated into the classical Bayesian estimation to some extent.Table 1
Comparison of the simulation result of the two methods.
σ2 η2 μ ^ B a y e s μ ^ B E R P MSE(μ^Bayes) M S E ( μ ^ B E R P ) OptNRps 40 0 0.0071 0.0419 0.1985 1.6958 19  50 0 -0.0185 -0.0720 0.2482 2.1204 19  60 0 -0.0178 -0.0796 0.2979 2.5442 19  40 5 2.1657 0.1930 4.9576 4.8515 14  50 5 2.1885 0.3104 5.0081 5.4469 13  60 5 2.1661 0.3347 5.0574 6.0391 13  40 10 3.0963 0.8185 9.7162 7.032 9  50 10 3.0684 1.0232 9.7688 7.7469 9  60 10 3.0920 1.0394 9.8186 8.4485 8Example 3.
When there is systematic bias in prior information, the CEP calculation based on BERP and classical Bayesian estimation is discussed in this example. We simulaten1=200 prior samples (X(1),Z(1))T from the bivariate normal distribution N(μ+η,Σ) and n2=5 realistic test samples (X(2),Z(2))T from the bivariate normal distribution N(μ,Σ), where η=10,10T,μ=0,0T, and Σ=diag(40,40). The parameters μx, μz, σx, and σz are estimated by BERP and classical Bayesian estimation, respectively. Based on the estimates of the population parameters, we use the numerical integration to get the CEP. Figure 2 shows the simulation results about CEP calculation.Figure 2
CEP calculated by different methods.As shown in Figure2, the true CEP is CEPtrue=7.4466, the CEP calculation based on BERP is CEPBERP=7.6532, and the CEP calculation based on classical Bayesian estimation is CEPBayes=8.4402. It is easy to conclude that the CEP calculation based on BERP is much closer to the true CEP when the systematic bias is within the bound. In addition, the CEP calculation without using the prior information is defined as CEPnoprior. If the prior information were effectively ignored, the CEPnoprior would be unreliable and unstable because of the large MSE of the estimation of μ without using prior information. Therefore, using BERP to estimate the parameters μx, μz will get more accurate calculation of CEP than classical Bayesian estimation when there is slight systematic bias in prior information. In addition, CEP calculation by BERP also outperforms the method with the prior information ignored.


## 6. Conclusions

In this paper, we have investigated the methods for performance evaluation of guidance weapon systems. Because of the small sample size of the realistic test data, we would fuse the prior information in the evaluation. However, by classical Bayesian estimation, the unknown systematic bias in prior information may cause large deviation for CEP calculation. For purpose of addressing it, a novel Bayesian method called BERP is proposed in this paper, and the corresponding optimization procedure is designed. In addition, we also give the bound of systematic bias for stopping using the heavily biased prior information.Within the bound of the systematic bias, theoretical analysis and simulation results prove that our new method has smaller estimate bias and MSE for estimating the mean of normal distribution than classical Bayesian estimation when there exists slight systematic bias in prior information. As for CEP calculation, the simulation results also validate that the CEP calculated by BERP is more accurate and reliable than the CEP calculated by classical Bayesian estimation. It can be concluded that a more accurate and reliable estimation of the CEP can be obtained via the BERP when the unknown systematic bias is within the bound.There is no obvious difference of the estimation accuracy between the two methods; BERP also has a good estimation performance when there is no systematic bias. Therefore, in order to get accurate and reliable evaluation results of guidance weapon systems, it is better to calculate the CEP via BERP than classical Bayesian estimation if the systematic bias is within the bound. In contrast, if the systematic bias is beyond the bound, we should stop fusing the biased prior information and evaluate the performance only by realistic test samples even if the sample size is small.



---

*Source: 5930109-2018-07-16.xml*