Electronic Journal: Southwest Journal of Pure and Applied Mathematics Internet: http://rattler.cameron.edu/swjpam.html ISSN 1083-0464 Vol. No. 01, 1995, pp. 37-45. Submitted: June 8, 1995. Published: October 26, 1995.

A PROBABILISTIC APPROACH TO THE ESTIMATION OF THE DISTRIBUTION

OF THE TARGET LOCATION POINT

M.A. Tabatabai

Abstract. A probabilistic approach is used to estimate the probability distribution of the target location point. The results can be used in missile-testing and target acquisition fire support model simulation. The target may be a mobile multiple launcher rocket, an antiaircraft gun or a ship in mid-ocean. In each case, there is a time lag between the detection of the target and the arrival of the projectile.

(1991) A.M.S. (MOS) Subject Classification Codes. 60, 62.

Key Words and Phrases. Probability density function, Gamma function.

I. Introduction

The target location problem frequently appears in missile-testing and artillery field programs. In this article, interest focuses on the probability distribution of the target location point for a multiple launcher rocket and a Howitzer. Our main concern is to find the probability density function of a mobile MLR location which is deliberately moving so as to confound prediction of its position. Let us now describe the situation by calling the two parties Red and Blue. The mission of the Red multiple launcher rocket is to strike by firing and move as quickly as possible in a specific direction with an average speed R. The reaction of the Blue commander is to predict the target location and launch a round of artillery at an anticipated target location.

Let the initial target location point (the point at which the Red initially fires) be represented by a bivariate normal random vector (X, Y), where X denotes the horizontal coordinate, Y denotes the vertical coordinate and X and Y are assumed to be independent random variables with a mean vector

(μ_{X}, μ_{Y}) = (0, 0)

and a standard deviation vector

(σ_{X}, σ_{Y}) = (100, 100)

where the mean and standard deviation for both X and Y are in meters.

Let T denote the Blue artillery reaction time. This random variable is the total response time from the Red multiple launch rocket firing to the first rounds of the Blue artillery landing on target. The distribution of the random variable T will be derived in Section 2. The Red displacement time using daylight values ranges from 1.5 minutes to 2 minutes.

Typeset by AMS-TEX

Department of Mathematical Sciences, Cameron University, Lawton, OK 73505

E-mail Address: mohammad@tiger.cameron.edu

\begin{matrix}c \\ ◯ \end{matrix}1995 Cameron University

II. Distribution of Reaction Time

A sample of 103 response times from the Red firing to the first rounds of the Blue artillery landing on target is included in the analysis. The data was provided by the Fire Support Center at Cameron University. The goal of this section is to find the distribution of Blue reaction time.

Before making a specificaiton of the distribution of response time, it is a good practice to look at the data in a graphical form in order to see if they could have come from one of a known family of distributions.

One of the simplest ways of looking at the data is to make a stem-and-leaf plot. Figure 1 displays the stem-and-leaf plot for Blue reaction time (in minutes). In this display we immediately see the asymmetry of the main part of the data. The distribution is positively skewed and the figure indicates the presence of some extreme outliers. In Figure 2, we have the Boxplot. It provides a visual impression of several important aspects of the underlying distribution of the Blue reaction time. From a Boxplot we can pick out the following features of a data set:

Location Skewness Spread Tail length Outlying data points

Stem-and-leaf of Blue Re N=103 Leaf Unit=0.10

2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 | \begin{matrix}1246777789999 \\ 00001112222233333444555566666666666777777888999 \\ 001111222334444446668899 \\ 00014 \\ 2 \\ 24578 \\ 79 \\ 49 \\ 44 \\ 1 \\ 7 \end{matrix}

Figure 1. Stem-and-leaf Plot of Blue Reaction Time

The mean reaction time for this data set is 4.554 minutes with a standard deviation of 2.396 minutes. The data ranges from a low of 2.160 minutes to a high of 16.770 minutes. The first, second and third quartiles are 3.300, 3.740, and 4.630 respectively and the trimmed mean is 4.238 minutes.

−−−−−
−−−−I + I −−−−−− ∗ ∗ ∗ 0 0 0 0 0 0
−−−−−
−−−−−−−−+−−−−−−−−−+−−−−−−−−−+−−−−−−−−−+−−−−−−−−−+−−−−−−−−−
3 . 0  6 . 0  9 . 0  1 2 . 0  1 5 . 0

Figure 2. Boxplot of Blue Reaction Time

Let T be the Blue reaction time and M be the Red movement time. More specifically, M is the difference between Blue reaction time and Red displacement time. If we use 2 minutes for the Red displacement time, then we can write the random variable M as M = T − 2. Figure 3 displays the stem-and-leaf graph for Red movement time.

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 | \begin{matrix}1246777789999 \\ 00001112222233333444555566666666666777777888999 \\ 001111222334444446668899 \\ 00014 \\ 2 \\ 24578 \\ 79 \\ 49 \\ 44 \\ 1 \\ 7 \end{matrix}

Figure 3. Stem-and-leaf Plot for Red Movement Time

To have a more reliable estimate for the parameters of underlying distrubtion of random variable M, we examine the Boxplot for the Blue reaction time. Figure 4 shows the Boxplot for M.

−−−−−  0
−−−I + I −−−−−− ∗ ∗ ∗ 0 0 0 0 0 0 0
−−−−−
−−+−−−−−−−−−+−−−−−−−−−+−−−−−−−−−+−−−−−−−−−+−−−−−−−−−−+−−−
0 . 0  3 . 0  6 . 0  9 . 0  1 2 . 0  1 5 . 0

Figure 4. Boxplot of Red Movement Time

The small hollow circles in the Boxplot indicate the presence of eight extreme outliers. There is a widespread awareness of the dangers posed by the occurrence of extreme outliers, which may be a result of keypunch errors, recording or transmission errors, misplaced decimal points, exceptional phenomena such as a lack of adequate training in firing rockets, or members of a different population slipping into the sample. Outliers occur frequently in real data. When the data contain outliers, one approach is to remove the outliers and then calculate your sample statistics. A second approach is to use robust and resistant estimators to estimate the population parameters. Robust estimators limit the influence of extreme observations and are insensitive to outliers. For more details about robust statistical techniques, see Hoaglin, Mosteller, and Tukey (1983); Roussaauw and Leroy (1987); Tabatabai and Argyros (1993). In this article we first remove the extreme outliers and then calculate some statistics.

Table 1 shows a summary of some relevant estimates of the parameters of Red movement time.

N  MEAN  MEADIAN TRMEAN STDEV  SEMEAN
95  1.960  1.690  1.838  1.124  0.115
MIN  MAX  Q1  Q3
0.160  5.850  1.270  2.410

Table 1. Summary of Red Movement Time

Let us now turn our attention to the theory underlying estimation of the density function of M. Let the range be partitioned into k class intervals I_{j}, j = 1, 2, …, k and I_{j} = (a_{j−1}, a_{j}), j = 1, 2, …, k where a_{0} < a_{1} < ⋯ < a_{k}. Here we assume that a_{j} − a_{j−1} = 1 minute for j = 1, 2, …, k. Define

f̂_{M}(m) = \{ \begin{matrix}f_{j}/[n(a_{j} − a_{j−1})] & for a_{j−1} < m ≤ a_{j}, j = 1, 2, …, k \\ 0, & elsewhere \end{matrix}

where f_{j} is the frequency of the class interval I_{j} and n = ∑_{j=1}^{k} f_{j} is the total number of observations. The

function f̂_{M} defined on (−∞, ∞) is the relative frequency histogram for M. We note that

(1)  f̂_{M}(m) ≥ 0 for all m ∈ (−∞, ∞),

(2)  ∫_{−∞}^{∞} f̂_{M}(m) dm = ∫_{a_{0}}^{a_{k}} f̂_{M}(m) dm = ∑_{j=1}^{k} ∫_{a_{j−1}}^{a_{j}} \frac{f_{j}}{n(a_{j} − a_{j−1})} dm = ∑_{j=1}^{k} \frac{f_{j}}{n} = 1.

Hence f̂_{M} defines a density function and can be used as an estimate of the underlying density function f_{M}. In particular, an estimate of the probability of an event A is given by

P̂(A) = ∫_{A} f̂_{M}(m)dm.

Moreover, P̂(A) is an unbiased and consistent estimate of P (A).

For our data on Red movement time, the estimate of probability distribution of M is given by

\underline{m}  \underline{f̂_{M}(m)}  \underline{m}  \underline{f̂_{M}(m)}
.5  13/103  8.5  2/103
1.5  47/103  9.5  2/103
2.5  24/103  10.5  1/103
3.5  5/103  11.5  0/103
4.5  1/103  12.5  0/103
5.5  5/103  13.5  0/103
6.5  0/103  14.5  1/103
7.5  2/103

where m_{j} = (a_{j−1} + a_{j})/2 is the class mark of class interval I_{j}. The mean of f̂_{M}(m) is given by

μ̂_{M} = ∫_{a_{0}}^{a_{k}} mf̂_{M}(m) dm = ∑_{j=1}^{k} ∫_{a_{j−1}}^{a_{j}} \frac{mf_{j}}{n(a_{j} − a_{j−1})} dm = ∑_{j=1}^{k} \frac{f_{j}}{n} (\frac{a_{j} + a_{j−1}}{2} ) = ∑_{j=1}^{k} \frac{f_{j}m_{j}}{n} .

So far we have answered th following questions: What does the underlying probability density function f_{M} look like? How do we estimate f_{M}? The next question that one needs to answer is: Does f_{M} look like a known or specified distribution? The probability density function of M is Gamma distribution with parameters α and β. We have

f_{M}(m) = \{ \begin{matrix}\frac{m^{α−1} exp(−m/β)}{Γ(α)β^{α}}, & m > 0 \\ 0, & m ≤ 0 \end{matrix}

where Γ(α) = ∫_{0}^{∞}x^{α−1}e^{−x} dx.

The mean μ_{M} = αβ and standard deviation σ_{M} = β\sqrt{α}. To estimate the parameters α and β, we use the method of moments technique. Hence, the moment estimators α̂ and β̂ are

α̂ = \frac{(sample mean)^{2}}{sample variance} and β̂ = \frac{sample mean}{α̂}.

If we use Table 1, then the moment estimates are

α̂ = \frac{(1.960)^{2}}{(1.124)^{2}} = 3.0419 and β̂ = \frac{1.960}{3.0419} = 0.6443.

Now, we use the chi-square goodness-of-fit test to test the following hypothesis?

H_{0} : f_{M}(m) = \frac{m^{3.0419−1} exp(−m/0.6443)}{Γ(3.0419)(0.6443)^{3.0419}}

vs

H_{1} : f_{M}(m) ≠ \frac{m^{3.0419−1} exp(−m/0.6443)}{Γ(3.0419)(0.6443)^{3.0419}}.

The chi-square test statistic is 4.9153. The associated degrees of freedom is 3 and the P − value is 0.1781. Therefore, we have to conclude that the Gamma distribution with parameters α = 3.0419 and β = 0.6443 does provide an excellent fit for the data. Similar conclusions were reached when we used the Kolmogorov-Smirnov test.

The maximum likelihood estimates of α and β are solutions of the following equations

ln α − (\frac{d}{dα}Γ(α)) /Γ(α) = ln \overline{m} − \frac{1}{n} ∑_{j=1}^{n} ln m_{j} and β = \frac{\overline{m}}{α}

where \overline{m} is the sample mean.

These type of equations are not easily solvable and it is necessary to resort to numerical methods. Let F_{T} (t) denote the distribution function of Blue reaction time, then we have

F_{T} (t) = P (T ≤ t) = P (M + 2 ≤ t)

= P (M ≤ t − 2) = ∫_{−∞}^{t−2} f_{M}(m) dm

= ∫_{−∞}^{t−2} \frac{m^{α−1} exp(−m/β)}{Γ(α)β^{α}} dm.

Then, the probability density function f_{T} (t) of Blue reaction time is given by

f_{T} (t) = \frac{d}{dt}F_{T} (t) = \{ \begin{matrix}\frac{(t − 2)^{α−1} exp(−(t − 2)/β)}{Γ(α)β^{α}}, & t > 2 \\ 0, & elsewhere \end{matrix}

where the values of α and β can be estimated by 3.0419 and 0.6443 respectively. Then mean for Blue reaction time is estimated by μ̂_{T} = 3.960 and the standard deviation is estimated by σ̂_{T} = 1.124.

III. Distribution of Target Location Point

Let the random vector (D_{1}, D_{2}) be the target location point and be defined by

D_{1} = X + RM cos θ
D_{2} = Y + RM sin θ

where R is the average speed of Red, M is the Red movement time, θ is the Red angle of movement and X, Y are independent normal random variables with μ_{X} = μ_{Y} = 0 meters and σ_{X} = σ_{Y} = 100 meters.

Our goal is to find the joint probability density function of (D_{1}, D_{2}). We begin this section by first finding the probability distribution of RM.

Let W = RM and F_{W} (w) denote the distribution function of the random variable W. Then

F_{W} (w) = P (W ≤ w) = P (RM ≤ w) = P (M ≤ \frac{1}{R}w)

= ∫_{−∞}^{\frac{1}{R} w} f_{M}(m) dm = ∫_{−∞}^{\frac{1}{R} w} \frac{m^{α−1} exp(−m/β)}{Γ(α)β^{α}} dm.

Then the probability density function f_{W} (w) is given by

f_{W} (w) = \frac{d}{dW} F_{W} (w) = \{ \begin{matrix}\frac{W^{α−1} exp(−w/(Rβ))}{Γ(α)(Rβ)^{α}}, & W > 0 \\ 0, & elsewhere. \end{matrix}

Thus, we can write

D_{1} = X + W cos θ
D_{2} = Y + W sin θ

where the random variable θ is uniformly distributed on the interval [0, 2π) with the probability density function g(θ) defined by

g(θ) = \{ \begin{matrix}\frac{1}{2π} & 0 ≤ θ < 2π \\ 0, & elsewhere \end{matrix}

and the random variables X, Y, θ, and W are independent reandom variables with joint probability density function given by

f(x, y, θ, w) = f_{X}(x) f_{Y} (y) g(θ) f_{W} (w)

= \{ \begin{matrix}\frac{w^{α−1} exp [ − w/(Rβ) − \frac{1}{2}(x^{2}/100^{2} + y^{2}/100^{2})]}{40000π^{2}Γ(α)(Rβ)^{α}}, & \begin{matrix}−∞ < & y < ∞ \\ −∞ < & x < ∞ \\ 0 ≤ & θ < 2π \\ 0 < & w < ∞ \end{matrix} \\ 0, & elsewhere \end{matrix}

Let

\begin{aligned}D_{1} & = g_{1}(X, Y, θ, W) = X + W cos θ, \\ D_{2} & = g_{2}(X, Y, θ, W) = Y + W sin θ, \\ D_{3} & = g_{3}(X, Y, θ, W) = θ, \\ D_{4} & = g_{4}(X, Y, θ, W) = W, \end{aligned}

A = {(x, y, θ, w) : −∞ < x < ∞, −∞ < y < ∞, 0 ≤ θ < 2π, 0 < w < ∞},
B = {(D_{1}, D_{2}, D_{3}, D_{4}) : −∞ < D_{1} < ∞, −∞ < D_{2} < ∞, 0 ≤ D_{3} < 2π, 0 < D_{4} < ∞}.

Define the transformation g : A ⟶_{onto}^{1−1} B as

g = (g_{1}(x, y, θ, w), g_{2}(x, y, θ, w), g_{3}(x, y, θ, w), g_{4}(x, y, θ, w))
= (D_{1}, D_{2}, D_{3}, D_{4})

then the inverse transformation is g^{−1} : B ⟶_{onto}^{1−1} A which can be written as

g^{−1} = (g_{1}^{−1}(D_{1}, D_{2}, D_{3}, D_{4}), g_{2}^{−1}(D_{1}, D_{2}, D_{3}, D_{4}), g_{3}^{−1}(D_{1}, D_{2}, D_{3}, D_{4}), g_{4}^{−1}(D_{1}, D_{2}, D_{3}, D_{4}))
= (x, y, θ, w)

where

X = D_{1} − D_{4} cos D_{3}
Y = D_{2} − D_{4} sin D_{3}
θ = D_{3}
W = D_{4}.

Let

J (\frac{x, y, θ, w}{D_{1}, D_{2}, D_{3}, D_{4}} ) = | \begin{matrix}\frac{∂x}{∂D_{1}} & · & · & · & \frac{∂x}{∂D_{4}} \\ \frac{∂y}{∂D_{1}} & · & · & · & \frac{∂y}{∂D_{4}} \\ \frac{∂θ}{∂D_{1}} & · & · & · & \frac{∂θ}{∂D_{4}} \\ \frac{∂w}{∂D_{1}} & · & · & · & \frac{∂w}{∂D_{4}} \end{matrix} |

be the Jacobian of inverse transformation g^{−1}. Then

J (\frac{x, y, θ, w}{D_{1}, D_{2}, D_{3}, D_{4}} ) = | \begin{matrix}1 & 0 & D_{4} sin D_{3} & − cos D_{3} \\ 0 & 1 & −D_{4} cos D_{3} & − sin D_{3} \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{matrix} | = 1.

The joint probability density function of random variables D_{1}, D_{2}, D_{3}, D_{4} is denoted by h(D_{1}, D_{2}, D_{3}, D_{4}) and is given by

h(D_{1}, D_{2}, D_{3}, D_{4}) = |J| · f (g_{1}^{−1}, g_{2}^{−1}, g_{3}^{−1}, g_{4}^{−1})

= \{ \begin{matrix}\frac{D_{4}^{α−1} exp \{−D_{4}/(Rβ) − [D_{1}^{2} + D_{2}^{2} + D_{4}^{2} − 2D_{4}(D_{1} cos D_{3} + D_{2} sin D_{3})] /20000\}}{40000π^{2}Γ(α)(Rβ)^{α}}, & (D_{1}, D_{2}, D_{3}, D_{4}) ∈ B \\ 0, & elsewhere. \end{matrix}

Finally, the distribution of the target location point is denoted by h(D_{1}, D_{2}) and is given by h(D_{1}, D_{2})

= \{ \frac{exp \{−(D_{1}^{2} + D_{2}^{2}) /20000\}}{40000π^{2}(Rβ)^{α}Γ(α)} ∫_{0}^{2π}∫_{0}^{∞}exp \{−[D_{4}^{2} − 2D_{4}(D_{1} cos D_{3} + D_{2} sin D_{3}) + \frac{20000}{Rβ}D_{4}] /20000\} dD_{4}dD_{3}, \begin{matrix}−∞ < & D_{1} < ∞ \\ −∞ < & D_{1} < ∞ \end{matrix} 0, elsewhere.

IV. Some Special Cases

Case 1. Suppose the initial point (X, Y) is predetermined and is equal to X = a, and Y = b. Then

D_{1} = a + W cos θ
D_{2} = b + W sin θ

Define the transformation S : A_{1} ⟶_{onto}^{1−1} B_{1} as

S = (S_{1}(w, θ), S_{2}(w, θ)) = (D_{1}, D_{2})

where

A_{1} = {(w, θ) : 0 < w < ∞, 0 ≤ θ < 2π}

and

B_{1} = {(D_{1}, D_{2}) : −∞ < D_{1} < ∞, −∞ < D_{2} < ∞}.

Then the inverse transformation S^{−1} : B_{1} → A_{1} is a one-to-one and onto function defined by

S^{−1} = (S_{1}^{−1}(D_{1}, D_{2}), S_{2}^{−1}(D_{1}, D_{2})) = (W, θ)

where

W = [(D_{1} − a)^{2} + (D_{2} − b)^{2}]^{\frac{1}{2}}

and

θ = arctan (\frac{D_{2} − b}{D_{1} − a}) .

The Jacobian of inverse transformation is

J (\frac{W, θ}{D_{1}, D_{2}} ) = | \begin{matrix}\frac{∂W}{∂D_{1}} & \frac{∂W}{∂D_{2}} \\ \frac{∂θ}{∂D_{1}} & \frac{∂θ}{∂D_{2}} \end{matrix} | = \frac{1}{[(D_{1} − a)^{2} + (D_{2} − b)^{2}]^{\frac{1}{2}}} .

The joint probability density function of the target location point conditioned upon knowledge of (X, Y) = (a, b) is denoted by h_{1}(D_{1}, D_{2}) and is given by

h_{1}(D_{1}, D_{2}) = |J (\frac{W, θ}{D_{1}, D_{2}} )| · f_{W} (S_{1}^{−1}) · g(S_{2}^{−1})

= \{ \begin{matrix}\frac{[(D_{1} − a)^{2} + (D_{2} − b)^{2}]^{\frac{α}{2} −1}}{2π(Rβ)^{α}Γ(α)} · exp \{− [(D_{1} − a)^{2} + (D_{2} − b)^{2}]^{\frac{1}{2}} /(Rβ)\} , & if (D_{1}, D_{2}) ∈ B_{1} \\ 0, & elsewhere. \end{matrix}

Case 2. If θ = 0, then the joint probability density function of D_{1} and D_{2} is denoted by h_{0}(D_{1}, D_{2}) and is equal to

h_{0}(D_{1}, D_{2})

= \{ \begin{matrix}\frac{exp \{−[D_{1}^{2} + D_{2}^{2}]\}}{20000π(Rβ)^{α}Γ(α)} ∫_{0}^{∞} D_{4}^{α−1} exp \{− [D_{4}^{2} − 2D_{1}D_{4} + \frac{20000}{Rβ}D_{4}] /20000\} dD_{4}, & \begin{matrix}−∞ < & D_{1} < ∞ \\ −∞ < & D_{2} < ∞ \end{matrix} \\ 0, & elsewhere. \end{matrix}

In general, if θ = C, then the joint probability density function of D_{1} and D_{2} is denoted by h_{C}(D_{1}, D_{2}) and is given by

h_{C}(D_{1}, D_{2})

= \{ \frac{exp \{−[D_{1}^{2} + D_{2}^{2}]\}}{20000π(Rβ)^{α}Γ(α)} ∫_{0}^{∞} D_{4}^{α−1} exp \{− [D_{4}^{2} − 2D_{4}(D_{1} cos C + D_{2} sin C) + \frac{20000}{Rβ}D_{4}] /20000\} dD_{4}, \begin{matrix}−∞ < & D_{1} < ∞ \\ −∞ < & D_{2} < ∞ \end{matrix} 0, elsewhere.

References

1. Hoaglin, D.C., Mosteller, F., Tukey, J., Understanding Robust and Exploratory Data Analysis, John Wiley & Sons, Inc. (1983).

2. Rousseeuw, P.J., and Leroy, A.M., Robust Regression and Outlier Detection, John Wiley & Sons, Inc. (1987).

3. Tabatabai, M.A., and Argyros, I.K., Robust Estimation and Testing for General Nonlinear Regression Models, Applied Math and Comp 58/1 (1993), 85–101.