Linear regression with an observation distribution model

Abstract

Despite the high complexity of the real world, linear regression still plays an important role in estimating parameters to model a physical relationship between at least two variables. The precision of the estimated parameters, which can usually be considered as an indicator of the solution quality, is conventionally obtained from the inverse of the normal equations matrix for which intensive computation is required when the number of observations is large. In addition, the impacts of the distribution of the observations on parameter precision are rarely reported in the literature. In this paper, we propose a new methodology to model the distribution of observations for linear regression in order to predict the parameter precision prior to actual data collection and performing the regression. The precision analysis can be readily performed given a hypothesized data distribution. The methodology has been verified with several simulated and real datasets. The results show that the empirical and model-predicted precisions match very well, with discrepancies of up to 6% and 3.4% for simulated and real datasets, respectively. Simulations demonstrate that these differences are simply due to finite sample size. In addition, simulation also demonstrates the relative insensitivity of the method to noise in the independent regression variables that causes deviations from the data distribution function. The proposed methodology allows straightforward prediction of the parameter precision based on the distribution of the observations related to their numerical limits and geometry, which greatly simplify design procedures for various experimental setups commonly involved in geodetic surveying such as LiDAR data collection.

Appendix A

In this appendix, the unbiasedness of the estimates of the slope parameter, m, and the y-axis intercept parameter, b, estimated by the linear least-squares model incorporating the x-observation distribution function, p(x), is proven. The limits of integration [a₁, a₂] are unequal (a₁ < a₂) to keep the proof general. The distribution function is assumed to be centred at x = 0.

The expectation of the estimated slope stems from the solution to the least-squares normal equations (Eq. 9) and substituting the specific forms given by Eqs. 30 and 31

$$ \begin{gathered} {\mathbf{x}} = {\mathbf{N}}^{ - 1} {\mathbf{U}} \\ = \left( {\begin{array}{*{20}c} {\int\limits_{{a_{1} }}^{{a_{2} }} {x^{2} p\left( x \right)dx} } & 0 \\ 0 & {\int\limits_{{a_{1} }}^{{a_{2} }} {p\left( x \right)dx} } \\ \end{array} } \right)^{ - 1} \left( {\begin{array}{*{20}c} {\int\limits_{{a_{1} }}^{{a_{2} }} {x \cdot y \cdot p\left( x \right)dx} } \\ {\int\limits_{{a_{1} }}^{{a_{2} }} {y \cdot p\left( x \right)dx} } \\ \end{array} } \right) \\ \end{gathered} $$

(A1)

results in

$$ E\left\{ {\hat{m}} \right\} = E\left\{ {\frac{{\int\limits_{{a_{1} }}^{{a_{2} }} {x \cdot y \cdot p\left( x \right)dx} }}{{\int\limits_{{a_{1} }}^{{a_{2} }} {x^{2} p\left( x \right)dx} }}} \right\} $$

(A2)

Since x has been assumed to be error-free, this expression reduces to

$$ E\left\{ {\hat{m}} \right\} = \frac{{E\left\{ {\int\limits_{{a_{1} }}^{{a_{2} }} {x \cdot y \cdot p\left( x \right)dx} } \right\}}}{{\int\limits_{{a_{1} }}^{{a_{2} }} {x^{2} p\left( x \right)dx} }} $$

(A3)

The numerator is analysed by substituting the model of Eq. 1

$$ \begin{aligned} & E\left\{ {\int\limits_{{a_{1} }}^{{a_{2} }} {x \cdot y \cdot p\left( x \right)dx} } \right\} = E\left\{ {\int\limits_{{a_{1} }}^{{a_{2} }} {x\left( {m \cdot x + b - e} \right)p\left( x \right)dx} } \right\} \\ &\quad = E\left\{ {\int\limits_{{a_{1} }}^{{a_{2} }} {m \cdot x^{2} p\left( x \right)dx} } \right\} + E\left\{ {\int\limits_{{a_{1} }}^{{a_{2} }} {b \cdot x \cdot p\left( x \right)dx} } \right\} - E\left\{ {\int\limits_{{a_{1} }}^{{a_{2} }} {x \cdot e \cdot p\left( x \right)dx} } \right\} \\ & \quad = mE\left\{ {\int\limits_{{a_{1} }}^{{a_{2} }} {x^{2} p\left( x \right)dx} } \right\} + bE\left\{ {\int\limits_{{a_{1} }}^{{a_{2} }} {x \cdot p\left( x \right)dx} } \right\} - E\left\{ {e\int\limits_{{a_{1} }}^{{a_{2} }} {x \cdot p\left( x \right)dx} } \right\} \\ \end{aligned} $$

(A4)

Under the stated assumptions, the numerator reduces to

$$ \begin{gathered} E\left\{ {\int\limits_{{a_{1} }}^{{a_{2} }} {x \cdot y \cdot p\left( x \right)dx} } \right\} = m\int\limits_{{a_{1} }}^{{a_{2} }} {x^{2} p\left( x \right)dx} + b\left( 0 \right) - E\left\{ {e\left( 0 \right)} \right\} \\ = m\int\limits_{{a_{1} }}^{{a_{2} }} {x^{2} p\left( x \right)dx} \\ \end{gathered} $$

(A5)

Division of Equation A5 by the denominator of Equation A3 results in the following

$$ E\left\{ {\hat{m}} \right\} = \frac{{m\int\limits_{{a_{1} }}^{{a_{2} }} {x^{2} p\left( x \right)dx} }}{{\int\limits_{{a_{1} }}^{{a_{2} }} {x^{2} p\left( x \right)dx} }} = m $$

(A6)

Therefore, the estimated slope parameter is unbiased.

The expected value of the intercept parameter is given by

$$ E\left\{ {\hat{b}} \right\} = E\left\{ {\frac{{\int\limits_{{a_{1} }}^{{a_{2} }} {y \cdot p\left( x \right)dx} }}{{\int\limits_{{a_{1} }}^{{a_{2} }} {p\left( x \right)dx} }}} \right\} = \frac{{E\left\{ {\int\limits_{{a_{1} }}^{{a_{2} }} {y \cdot p\left( x \right)dx} } \right\}}}{{\int\limits_{{a_{1} }}^{{a_{2} }} {p\left( x \right)dx} }} $$

(A7)

Since the denominator is unity by definition (Eq. 20), it is sufficient to analyse only the numerator

$$ \begin{aligned} & E\left\{ {\int\limits_{{a_{1} }}^{{a_{2} }} {y \cdot p\left( x \right)dx} } \right\} = E\left\{ {\int\limits_{{a_{1} }}^{{a_{2} }} {\left( {m \cdot x + b - e} \right)p\left( x \right)dx} } \right\} \\ & \quad = m\int\limits_{{a_{1} }}^{{a_{2} }} {x \cdot p\left( x \right)dx} + b\int\limits_{{a_{1} }}^{{a_{2} }} {p\left( x \right)dx} - E\left\{ {e\int\limits_{{a_{1} }}^{{a_{2} }} {p\left( x \right)dx} } \right\} \\ & \quad = m\left( 0 \right) + b\left( 1 \right) - E\left\{ {e\left( 1 \right)} \right\} = b \\ \end{aligned} $$

(A8)

Therefore, the estimated y-axis intercept parameter is also unbiased.