**Purpose.**:
To develop and validate a method of predicting visual function from retinal nerve fiber layer (RNFL) structure in glaucoma.

**Methods.**:
RNFL thickness (RNFLT) measurements from scanning laser polarimetry (SLP) and visual field (VF) sensitivity from standard automated perimetry were made available for 535 eyes from three centers. In a training dataset, structure–function relationships were characterized by using linear regression and a type of neural network: radial basis function customized under a Bayesian framework (BRBF). These two models were used in a test dataset to (1) predict sensitivity at individual VF locations from RNFLT measurements and (2) predict the spatial relationship between VF locations and positions at a peripapillary RNFLT measurement annulus. Predicted spatial relationships were compared with a published anatomic structure–function map.

**Results.**:
Compared with linear regression, BRBF yielded a nearly twofold improvement (*P* < 0.001; paired *t*-test) in performance of predicting VF sensitivity in the test dataset (mean absolute prediction error of 2.9 dB [SD 3.7] versus 4.9 dB [SD 4.0]). The predicted spatial structure–function relationship showed better agreement (*P* < 0.001; paired *t*-test) with anatomic prior knowledge when the BRBF was compared with the linear regression (median absolute angular difference of 15° vs. 62°).

**Conclusions.**:
The BRBF generates clinically useful relationships that relate topographical maps of RNFL measurement to VF locations and allows the VF sensitivity to be predicted from structural measurements. This method may allow clinicians to evaluate structural and functional measures in the same domain. It could also be generalized to use other structural measures.

^{ 1 }Despite their limitations, these techniques are currently central to the diagnosis and management of glaucoma. It would, therefore, be beneficial if structure and function measurements were directly linked in some way, allowing clinicians to corroborate decisions by considering the measurements in tandem.

*R*

^{2}, Pearson, or Spearman coefficients.

^{ 2 –5 }This approach has two major flaws: The use of summary data loses spatial information and may reduce power, and these association measures and regression models assume a particular shape of the relationship. Furthermore, these analyses fail to take account of spatial associations in the data, an integral attribute of glaucomatous loss.

^{ 6 –12 }we propose that much more could be gained by analyzing the data in its high-resolution form. For example, in scanning laser polarimetry (SLP), RNFLT estimates are yielded over an image space of several thousand pixels. These are high-dimensional data and any method linking structural measures to the 50 or so individual points in the VF ideally should take this into account. Moreover, individual points from both structure (pixel or sector values) and function (areas of VF or individual locations) are more likely to interact as groups rather than single independent measurements.

^{ 13 –16 }This class of networks differs from the ANNs that use a multilayer perceptron (MLP) approach, in that the nonlinearity of the model is embedded in a hidden layer of the network. This hidden layer consists of basis functions. The key element of the RBF is that predictions can be made in a multidimensional space that consists of mini-distributions of possible predictions, analogous to a kernel or window when interpolating data. RBF networks are particularly useful when mapping two- or three-dimensional images where interpolation is a prerequisite, with some image features being preserved while others are not necessarily mapped exactly. In this way, they seem most suitable to the problem of mapping the individual points from the structural and functional measures. Moreover, the RBF method developed for our purposes is formulated under a Bayesian probabilistic framework to tackle the considerable variability in measurements as well as to form an automatic learning process. This customization leads to what we denote a Bayesian RBF (BRBF) network.

^{ 9 }In addition, we illustrate how the method can be extended to high-resolution measurements.

*P*< 1% loss or greater, three or more contiguous points with

*P*< 5% loss or greater, or a 10 dB difference across the nasal horizontal midline at two or more adjacent points in the total deviation plot. Ocular hypertension subjects were recruited from clinic on the basis of having a raised IOP >21 mm Hg in two consecutive visits and having at least two normal VFs. Optic disc appearance was not used to categorize subjects. However, the optic disc was evaluated to exclude anatomic abnormalities such as coloboma or drusen. For all participants, one eye was randomly selected for study if both were eligible. All subjects had visual acuity of 20/40 or better, with ametropia <7 diopters, and had no other significant ocular abnormality or concomitant ophthalmic condition.

^{ 17 }In short, healthy subjects had normal VF, optic discs without structural abnormalities, IOP < 21 mm Hg, no previous history of ocular disease, and no family history of glaucoma. The patients had a glaucomatous appearance of the optic disc and a corresponding nerve fiber bundle VF defect, as described by Keltner et al.

^{ 18 }with SAP. All subjects had a visual acuity of 20/40 or better and had no other significant ocular abnormality. For all participants, one eye was randomly selected for the study.

^{ 19 }Two hundred thirty healthy subjects (230 eyes) and 76 patients (76 eyes) diagnosed with glaucomatous optic neuropathy were selected from this population under strict measurement quality criteria (described later). Only one randomly selected eye per subject was used throughout. The criteria used for defining glaucomatous VF loss was an abnormal result of the glaucoma hemifield test (GHT) plus one or more of the following VF defect classifications (top row of test points in the 24-2 pattern were excluded to reduce the effect of lid artifact), which could not be explained on the basis of nonglaucomatous ocular, or neurologic, causes: (1) at least four contiguous points on the pattern deviation plot depressed at the

*P*< 0.5% level; (2) at least two horizontal points in nasal step locations with a pattern deviation plot depressed at the

*P*< 0.5% level; and (3) advanced glaucomatous field loss (hemispheric or severe generalized field loss, with residual temporal or central islands). Glaucoma was diagnosed if glaucomatous VF loss was present, combined with matching neural rim loss at the optic disc, and gonioscopy showed no evidence of angle closure, rubeosis, or secondary glaucoma, other than pseudoexfoliation syndrome.

^{ 20 }This instrumentation estimates the thickness of the peripapillary RNFLT by measuring the summed retardation in a polarized scanning laser beam reflected from the fundus. Retardation measurements at various points around the optic disc are used to construct a thickness map of the retinal nerve fiber layer (RNFL) in micrometers. SLP with variable corneal compensation (VCC) has been shown to improve the estimates of RNFLT compared with earlier versions of the technology.

^{ 21 }All SLP images from all three centers had quality scores ≥8 and typical scan scores ≥80. Single SLP images were available from all the subjects from all the centers. The 64-sector thickness profile on the peripapillary annulus (provided by the GDx software) for each subject and the raw images (.mif files) were transferred to a secure database.

^{ 22 }and are reviewed in the 1. The BRBF model was developed and written in commercial software (MATLAB, ver. 7.4.0 R2007a; The MathWorks, Inc., Natick, MA). An executable version of this code is freely available from the authors.

*(where*

_{d}*d*is one of the 52 locations in a 24-2 HFA VF) from a series of RNFLT values denoted

**x**

*for*

_{i}*i*from 1 to

*m*. This can be expressed by the following oversimplified but illustrative equation: where

*c*is a constant offset. The symbol on

_{d}**ŷ**

*indicates that it is a prediction rather than the real measured value denoted by*

_{d}**y**

*. In this example, the equation has 64 peripapillary thickness profile values (*

_{d}*m*= 64), each with its own coefficient that quantifies the contribution of each

**x**value to the prediction. Thus, each

**y**value can be predicted by a linear combination of

**x**values. With some actual data we can find some real numbers for the

**w**terms by least-squares regression. This calculation attempts to fit an equation that minimizes the difference between the predicted and measured values. It yields an individual equation for each

**y**value that can be enumerated across all the points, to predict a complete VF from a given vector of

**x**values. This classic multiple linear regression can be adapted to select only those

**x**values that are statistically significant for the prediction of

**y**values, by using techniques such as stepwise multiple regression

^{ 23 }with the forward-selection scheme. In the linear model described by equation 1, the

**w**values (divided by their standard errors) quantify the amount of meaningful contribution made by

**x**values to predict the

**y**values. The largest absolute

**w**term (with respect to the variability in estimating the term) would indicate the

**x**value that affects the

**y**

*value the most, in the sense that a change in this*

_{d}**x**value results in the largest change in

**y**

*. Similarly, the next largest absolute*

_{d}**w**term would indicate the second most important term and so on. Equivalently, one could look at the relationship between the

**y**

*value and each*

_{d}**x**value separately and simply calculate the correlation coefficient based on the raw data or on the ranks of the data (Spearman's ρ) and end up with a similar result. Loosely speaking, this is the approach of Gardiner et al.,

^{ 9 }who used neuroretinal rim area estimated from scanning laser ophthalmoscope measurements (Heidelberg Retina Tomograph [HRT]; Heidelberg Engineering, Heidelberg, Germany) as the surrogate structural measurement for glaucomatous damage. Particularly, in the linear model implemented for this study, the VF sensitivity was unlogged from the decibel value and the prediction from this model was converted back to the decibel scale when compared with the measured sensitivity. We will refer to this method as a classic linear model.

**x**value is independent of all the other

**x**values,

^{ 24 }whereas, in reality, the

**x**and

**y**values are topographically and physiologically related and may interact as groups. Although one could try to demarcate these groups on the basis of a physiological relationship between the

**x**values or of an anatomic map, it would be preferable for the groups to be learned from the data rather than imposing any relationship from incomplete prior knowledge. Second, this approach assumes that the relationship between

**y**and

**x**is either linear or becomes linear after some transform (typically logarithmic). In reality, the relationship between

**y**and

**x**may be more complex, with the nature of the association probably changing across the measurement range of

**y**. Put simply, at different stages of disease, the apparent relationship between

**y**and

**x**could switch from linear, to noisy, to curvilinear and occasionally could be censored because of the imprecision and range of the measurement. This notion that the association may change at different levels of functional loss has also been asserted in published studies.

^{ 25,26 }The third difficulty with the classic linear model is that outlier points exert an overly strong influence and can yield a false association.

**y**and

**x**without the limiting assumptions associated with the classic linear model just described. As an illustrative example, consider one

**x**value, say

**x**

_{1}, appearing to be correlated to a

**y**value, say

**y**

_{1}. This apparent relationship might be explained by the very strong relation of

**x**

_{1}to

**x**

_{2}, which in turn is very strongly associated with

**y**

_{1}. Thus, the apparent

**x**

_{1}and

**y**

_{1}relationship may well be much weaker or not even significant. As it will be shown in the derived structure–function relationship (see Figs. 3 4–5), this covariance in the relationship between the

**y**and

**x**values is modeled better with the RBF approach. Furthermore, the central idea of the RBF is the basis functions, each of which performs very much like a dynamic window or kernel that moves across the data, both spatially and at various stages in disease severity, identifying groups of measurements that appear to behave in a similar pattern. The non-normalized Gaussian basis function used in this study has an activation field that has a center—that is, a particular input value at which it has a maximal output. The output tails off as the input moves away from this point. In this way, those hidden basis functions that have centers similar to the input

**x**patterns will have stronger activation and will thus contribute more to the prediction of

**y**. On the other hand, those basis functions with weak activation will be isolated and will not affect the prediction. Moreover, the RBF learns the parameters from the data and this is customized for our purpose by using a BRBF. The method makes predictions in multiple dimensions by extending the standard relevance vector machine

^{ 27 }where the output is restricted to one dimension. The VF sensitivities at different locations are implicitly correlated by sharing the same basis functions and parameters of weights (see α in the 1). In the BRBF model, the original VF sensitivity in decibels was used.

*M*and the number of data points

*N*of a dataset satisfy

*M*≥

*N*, the dimensionality of the dataset can be reduced from

*M*to

*N*− 1, with minimal loss of information using principal component analysis (PCA),

^{ 28 }because the

*N*data points span at least one linear

*hyperplane*in this

*M*-dimensional space. Using this technique the 16,512 dimensional SLP image vectors were reduced to 228 dimensions for analysis and transformed back to the original SLP vectors for the purpose of visualization and evaluation of results. Using the classic linear model on a reduced SLP image vector with 228 elements and a 24-2 VF with 52 sensitivity values will still result in a prohibitive number of weights (11,856) to be fitted, which will cause significant overfitting to the noise in the data.

^{ 29,30 }It is also well known that in developing a model, the input variables and weights selected may vary across different samples.

^{ 31 }Therefore, the models were developed on the MEH and REH data alone, leaving the BMES data as a test dataset.

^{ 9 }The structure–function relationship was defined by the corresponding functional change given a subtle structural change, which is mathematically modeled by the derivative of the BRBF (1). This derivative describes the relationship between

**y**(VF sensitivity) and

**x**(RNFLT) at each VF location. The other output from this analysis consisted of point-by-point predictions of each subject's VF, as represented by the HFA grayscale (which was replicated for this purpose). These outputs were considered for (1) the classic linear model, (2) the BRBF I model (based on the 64 summary RNFLT values output from the GDx software), and (3) for the BRBF II model (based on the reduced 228 dimensional data derived from the PCA on the 16,512 individual pixel retardation values in the broad annulus centered on the ONH).

Data Set | Ocular Status | Sex (n) | Age (y) | MD (dB) | PSD (dB) | NFI |
---|---|---|---|---|---|---|

MEH | Healthy | 19 M, 11 F | 40.6 ± 15.6 (21, 75) | 1.26 ± 0.76 (0.09, 3.30) | 1.47 ± 0.28 (1.12, 2.27) | 16 ± 8 (2, 30) |

OHT | 13 M, 21 F | 60.7 ± 11.1 (21, 75) | 1.07 ± 0.82 (0.10, 3.06) | 1.35 ± 0.27 (1, 2.1) | 20 ± 10 (5, 43) | |

Glaucoma | 28 M, 15 F | 60.8 ± 13.1 (31, 84) | −4.02 ± 2.55 (−12.00, 1.15) | 4.93 ± 2.93 (1.49, 12.50) | 42 ± 18 (11, 80) | |

REH | Healthy | 23 M, 23 F | 60.4 ± 12.1 (23, 77) | 0.38 ± 0.99 (−1.55, 2.73) | 1.63 ± 0.26 (1.13, 2.30) | 21 ± 9 (2, 43) |

Glaucoma | 47 M, 29 F | 62.2 ± 10.1 (30, 82) | −9.52 ± 8.43 (−30.39, 1.25) | 8.35 ± 4.32 (1.99, 15.92) | 63 ± 21 (21, 98) | |

BMES | Healthy | 99 M, 131 F | 69.3 ± 6.5 (60, 87) | 0.21 ± 1.07 (−1.26, 3.03) | 1.53 ± 0.29 (1.12, 2.55) | 19 ± 10 (2, 49) |

Glaucoma | 25 M, 51 F | 72.0 ± 6.3 (61, 78) | −7.94 ± 6.55 (−29.67, 1.56) | 6.97 ± 3.76 (1.67, 15.56) | 65 ± 24 (27, 98) |

*P*< 0.001; paired

*t*-test) in performance (mean [SD]): BRBF I 2.9 dB (3.7 dB) and BRBF II 2.8 dB (3.8 dB). In the BRBF I and II models, the training process described in the 1 selected 49 and 73 basis functions, respectively, in the hidden layer.

^{ 32 }across all locations from two VF tests from each of 49 individuals have been superimposed on Figure 1b. On inspection, these limits are similar to the 90% prediction limits when using the BRBF to predict the VF from a GDx RNFLT measurement of the same individual. Note that predictions at higher sensitivities (>30 dB) tend to be slightly lower than the actual values, whereas at lower sensitivities (<20 dB), the predictions tend to be higher.

^{ 33 }This correspondence improves (

*P*< 0.001; paired

*t*-test) with the BRBF I model (median absolute angular difference of 15°) in general, and improves further (

*P*< 0.001; paired

*t*-test) still with the BRBF II model (median absolute angular difference of 12°), especially at the points around central vision and blind spot. This result is consistent with the BRBF models learning to encode the structure–function relationship during the training.

^{ 34 –40 }Most of these applications have tended to use a conventional MLP ANN. Technically, there are several advantages of the RBF over MLP ANNs. For example, with the latter, there is an input, and the distributed pattern “lights up” all hidden units, to contribute to the prediction of the output, which makes them combine and interfere with each other. This typically yields a highly nonlinear training process with mathematical difficulties that result in a slow convergence of a training procedure.

^{ 16 }Moreover, the complexity of hidden unit patterns causes difficulty when interpreting the result, because the connection among units and the hidden-unit output do not have any physical or realistic meaning: they are simply numbers that can produce a correct output. Hence, MLP ANNs are less suitable for the mapping of points in different measurement spaces, which requires a detailed understanding of the hidden layer output and other manipulation (e.g., derivatives) within the network. In this study, although an MLP provided a prediction accuracy comparable to the BRBF (data not shown), the spatial aspects of the structure–function relationship were poorly predicted. This deficiency probably resulted from complex interactions among the large number of weights in the MLP model. In contrast, RBF substitutes hidden layers in MLP with a set of basis functions, which leaves the solution of weights in a linear space. The Gaussian basis function forms a local representation in the hidden unit, each of which can be understood as a representative of similar input patterns (Fig. 2). With a given RNFLT input, only a few representatives will be activated and contribute to the VF prediction (Fig. 2). Furthermore, the RBF handled within a Bayesian probabilistic framework,

^{ 27 }as developed for our purposes, is unlike most other ANN approaches because it is independent of any subjective input parameter and thus requires no model validation on a test dataset, provided that the training dataset is representative and sufficiently large to enable modeling of the many different states of the VF. However, for this study, the trained model is still tested on a separate, independent dataset to illustrate the generalization of the model performance. One technical limitation of the current BRBF model, despite its good performance, is that it assumes that the variability in the VF measurements is largely Gaussian, which is not optimal, given that it is often skewed and heavily tailed.

^{ 32 }

^{ 33 }derived by overlaying an appropriately scaled VF grid on RNFL photographs and tracing nerve fiber bundles or defects from various VF locations to the ONH margin. Gardiner et al.

^{ 9 }produced a topographical map of the relationship between sectors of the ONH and locations in the VF by considering the linear correlation between Heidelberg retina tomography neuroretinal rim area and VF sensitivity at each point in the VF. The classic linear model is akin to that derived by Gardiner et al. except that the SLP measurements from GDx imaging are used. The maps derived from the BRBF models (Figs. 4, 5) indicated a closer concordance with the anatomically derived landmark than did the map derived from the classic linear model (Fig. 3). In addition, the predictions on point-wise VF sensitivities in the validation dataset were, on average, better with the BRBF method than with the classic linear model. A plausible explanation for the improvements afforded by the BRBF technique is that it models the spatial and quantitative structure–function relationship more precisely. The technical and statistical advantages of the BRBF model over the classic linear model support this notion. For example, the BRBF considers that the VF points, and indeed RNFLT values acquired from different discreet areas, interact as groups rather than as independent measurements. The BRBF also accommodates the covariance and nonindependence of the measurements. It is less affected by outlying observations and makes no assumption about the linearity of relationships. Of course, one difficulty in generating any type of map, driven by data or anatomic observation, is the restriction in the sampling of VF points. These are, for example, probably not optimally placed for estimating RGC density nearer the fovea.

^{ 32 }This similarity suggests that, on average, a VF predicted by the BRBF from RNFLT values has measurement noise equivalent to that found in a newly measured field. This finding is not as exciting as it may first appear, because it is well established that the measurement noise in VFs is already very high, prohibiting straightforward clinical diagnosis of glaucomatous defects and monitoring progression. Nevertheless, this finding illustrates that the range and scale of the average predictive performance of the BRBF model is much better than the classic linear model approach, which completely fails to predict the full range of VF values (Fig. 1a).

^{ 32 }on average, are very small at the normal or healthy end. At the damaged end, the retest values tend to be higher, but the median BRBF predictions are slightly higher than the retest values. This likely reflects the difficulty the prediction from RNFLT images have in identifying small, focal defects. Moreover, the floor effect in the VFs and SLP images

^{ 41,42 }and the atypical scan pattern in SLP images, which may be associated with glaucoma severity,

^{ 43 }may be additional causes of the overestimation at the lower end of the VF sensitivity. Furthermore, because the diagnosis of subjects in REH datasets includes structural criteria, the normal subjects in the training dataset may have supernormal structure and the glaucomatous subjects have greater than average structural damage. This potential bias on subject selection may distort the structure–function relationship. Therefore, the training process may be improved by including a range of glaucoma severity in subjects defined only by VF loss. Another factor that may confound the reported prediction accuracy is that the models were trained with VFs tested with both SITA and full-threshold programs, whereas they were evaluated on VFs tested only with full-threshold program. This may, in part, account for the observed tendency for the prediction to overestimate the sensitivities in the test dataset due to the higher sensitivities obtained with SITA compared with the full-threshold program.

^{ 44 }However, this effect would be small because the sensitivity difference with two programs is just 1.3 dB, on average.

^{ 44 }

^{ 45,46 }For example, loss of function without loss in structure that does not adhere to a particular structure–function model's prediction may be an indicator of a nonglaucomatous process, measurement imprecision, or some artifact in the image of the structure or function test. Our long-term goal is to use the BRBF technique to provide a relevant clinical tool that indicates concordance between the VF and the chosen surrogate measure for structural loss. For example, when a VF and structural measure derived from GDx, HRT, or optical coherence tomography (OCT) is available, a chart mapped in VF space will be provided indicating areas where the measurements are in concordance (within a certain range of accuracy and precision) and where they are not, it could provide clinically useful information about the reliability of the individual measurements or diagnostically useful information.

^{ 30 }We had access to three large, independent datasets, each collected at one of three clinical centers. The inclusion criteria for glaucomatous and healthy subjects were generally consistent across the three samples. However, as the purpose of this study was not to determine diagnostic performance, the precise definitions for glaucoma were less important. In fact, the mixture of data can be viewed as an advantage in the study design. However, further testing on different datasets, especially where realistic estimates of measurement precision have been performed (from test–retest measurement), is still required and is being undertaken in our laboratory and elsewhere.

**y**from a structural one

**x**:

**y**= (

**x**). This calculation sets a mathematical framework for the question: what would be the corresponding functional change (Δ

**y**) for the structural change Δ

**x**: Δ

**y**= (

**x**+ Δ

**x**) − (

**x**)? Because we are interested in subtle structural change to model the slow progression of RNFL damage, we assume that Δ

**x**becomes very small and tends toward 0. The structure–function relationship is then defined with the general equation where is differentiable with regard to

**x**, and ∇

**is the gradient of at**

_{x}**x**. Since we are sampling subjects rather than considering the whole population, the final term in this equation must be expressed as a statistical expectation of ∇

**or, for simplicity, the mean of ∇**

_{x}**.**

_{x}^{27}where the output is originally only one-dimensional. The extension was similar to the model derived by Thayananthan et al.

^{47}In particular, VF sensitivity is assumed to be where

**y**

*is the*

_{d}^{n}*d*th element in the measurement of the

*n*th subject,

**x**

*is the RNFLT measurement (64-point profile or SLP image) vector,*

^{n}**w**

*is a weight vector, ε*

_{d}*is an additive 0-mean Gaussian noise*

_{d}*N*(0, δ

_{d}^{2}) with variance δ

_{d}^{2}, and the radial basis function vector φ(

**x**

*) is defined to be*

^{n}*M*+1 dimensional for

*M*bases: φ(

**x**

*) = [1, φ*

^{n}_{1}(

**x**

*), φ*

^{n}_{2}(

**x**

*),…, φ*

^{n}*(*

_{M}**x**

*)]*

^{n}^{T}, where each element is a radial basis function with center θ

*and an isotropic covariance If all weight vectors*

_{m}**w**

*are organized into a matrix*

_{d}**W**columnwise, then, using Bayesian methodology, we assign a prior over

**W**where Ψ is a diagonal matrix whose elements are α

_{1}

^{−1}, α

_{2}

^{−1},…, α

_{M}^{−1}on diagonal, and 0 otherwise. Each α

_{m}^{−1}represents the average variance of the weights for the

*m*th basis.

**w**

*is independent from any other weight vectors, given*

_{d}**Y**,

**α**, and

**β**. Consequently where the covariance matrix and mean for the posterior distribution of

**w**

*are where Φ(*

_{d}**X**) is a matrix: Φ(

**X**) = [φ(

**x**

^{1}),…, φ(

**x**

*)]*

^{N}^{T}.

*p*(

**α**,

**β**|

**Y**) ∝

*p*(

**Y**|

**α**,

**β**)

*p*(

**α**)

*p*(

**β**) where With no prior knowledge on

**α**and

**β**, these two parameters are assumed to be “uniformly” distributed so

*p*(

**α**) and

*p*(

**β**) have little impact on

*p*(

**α**,

**β**|

**Y**). When an approximation similar to that of RVM is used,

^{27}the hyperparameters are optimized by setting the derivative of equation 10 to 0 where μ

_{wd}

^{m}is the

*m*th element of the vector μ

_{wd}, and Σ

_{wd}

^{mm}is the diagonal element at

*m*th row and column in Σ

_{wd}.

*are initialized to contain all*

_{m}**x**

*. As with an RVM, many of α*

^{n}*become infinite during the training process so the corresponding radial bases are removed accordingly. The radial basis parameter η*

_{m}*is optimized by gradient descent.*

_{m}*φ*

_{m}*(*

_{m}**x**

*) (*

^{n}**x**

*− θ*

^{n}*), according to equation 4.*

_{m}