SMOOTHING L O C A L L Y R E G U L A R PROCESSES BY BAYESIAN NONPARAMETRIC METHODS, WITH APPLICATIONS ACID RAIN DATA TO ANALYSIS by HON B.Sc, A THESIS WAI M A CONCORDIA UNIVERSITY, 1984 S U B M I T T E D IN P A R T I A L F U L F I L L M E N T THE REQUIREMENTS FOR T H E DEGREE OF MASTER OF SCIENCE in T H E F A C U L T Y O F G R A D U A T E STUDIES Department of Statistics We accept this thesis as confirming to the required standard THE UNIVERSITY O F BRITISH September (S) COLUMBIA 1986 H O N W A I M A , 1986 OF In p r e s e n t i n g t h i s t h e s i s i n p a r t i a l f u l f i l m e n t of requirements f o r an advanced degree a t the the University of B r i t i s h Columbia, I agree t h a t the L i b r a r y s h a l l make it f r e e l y a v a i l a b l e f o r reference and study. I further agree t h a t p e r m i s s i o n f o r e x t e n s i v e copying o f t h i s t h e s i s f o r s c h o l a r l y purposes may be granted by the head of department or by h i s or her r e p r e s e n t a t i v e s . my It is understood t h a t copying or p u b l i c a t i o n of t h i s t h e s i s f o r f i n a n c i a l gain s h a l l not be allowed without my permission. Department of The U n i v e r s i t y of B r i t i s h Columbia 1956 Main M a l l Vancouver, Canada V6T 1Y3 DE-6 C3/81') written ii ABSTRACT We consider the problem of recovering an unknown smooth function from the data using the Bayesian nonparametric approach proposed by Weerahandi and Zidek (1985). Selected nonparametric smoothing methods are reviewed and compared with this new method. At each value of the independent variable, the smooth function is assumed to be expandable in a Taylor series to the pth order. Two methods, cross-validation and "backfitting" are used to estimate the a priori unspecified hyperparameters. Moreover, a data-based procedure is introduced to select the appropriate order p. Finally, an analysis of an acid-rain, wet-deposition time series is included to indicate the efficacy of the proposed methods. iii T A B L E OF CONTENTS ABSTRACT T A B L E OF CONTENTS LIST OF TABLES LIST OF FIGURES ACKNOWLEDGEMENTS 1. INTRODUCTION ii iii v vi viii 1 2. SURVEY OF SELECTED NONPARAMETRIC SMOOTHING METHODS 3 2.1 T H E KERNEL METHOD 3 2.2 T H E SPLINE SMOOTHING METHOD 6 2.3 T H E EXPONENTIAL SMOOTHING METHOD 9 2.4 SOME ROBUST SMOOTHING METHODS 14 2.5 T H E MULTIVARIATE EXTENSION 16 2.6 DISCUSSION 19 3. THEORY OF BAYESIAN NONPARAMETRIC APPROACH FOR SMOOTHING LOCALLY REGULAR PROCESSES 20 3.1 T H E GENERAL THEORY 20 iv 3.2 ESTIMATING T H E "SMOOTHING PARAMETER" C FOR LOCALLY REGULAR STRUCTURAL MODELS 24 (I) CROSS-VALIDATION M E T H O D 24 (II) 27 "BACKFITTING" M E T H O D (III) CONCLUSIONS 3.3 CHOOSING T H E "OPTIMAL" ORDER P 4. APPLICATIONS 4.1 INTRODUCTION 29 29 32 32 4.2 COMPARISON OF T H E CROSS-VALIDATION AND T H E BACKFITTING METHODS 32 4.3 RESULTS OF CHOOSING T H E "OPTIMAL" ORDER P 35 5. CONCLUSIONS 38 BIBLIOGRAPHY 40 APPENDIX 47 V LIST OF TABLES T A B L E I: Estimated values of XV(c) M.A. pll data T A B L E II: Estimated values of C for the pH77 data p T A B L E III: Estimated values of C for a locally constant fit to the pH77 and 33a P for the M.A. pH data 33b 34b T A B L E IV: Estimated values of PSE(p) vation ahead) for the pH77 data (predict one obser35a T A B L E V: Estimated values of PSE(p) of 10 observations ahead) for the pH77 data (predict the average 35b T A B L E VI: Estimated values of PSE(p) for the Lobster data predict one observation ahead (top) and predict the average of 10 observations ahead (bottom) 36a T A B L E VII: Estimated values of PSE(p) for the sinusoid (y = 3cos rc - 5sin re 4-6) data predict one observation ahead (top) and predict the average of 10 observations ahead (bottom) 36d vi LIST OF FIGURES FIG. 1 (A-C): A comparison of the locally polynomial fits of different order p to the pH77 data (using the cross-validation method) 33c FIG. 2 (A-C): A comparison of the locally polynomial fits of different order p to the pH77 data (using the backfitting method) 33d FIG. 3: A comparison of the locally constant fits (top) and quadratic fits (bottom) using the cross-validation method (solid line) and the backfitting method (dotted line) for the pH77 data , 34a FIG. 4 (A-C): A comparison of the locally polynomial fits of different order p to the M . A . pH data (using the cross-validation method) 34c FIG. 5 (A-C): A comparison of the locally polynomial fits of different order p to the M.A. pH data (using the backfitting method) 34d FIG. 6: A comparison of the locally constant fits (top) and quartic fits (bottom) using the cross-validation method (solid line) and the backfitting method (dotted line) for the M.A. pH data 34e FIG. 7: A comparison of the credibility intervals for locally quadratic (solid line) and cubic (broken line) fits for the pH77 (top) and M.A. pH (bottom) data (using the cross-validation method) 34f FIG. 8A: A comparison of the locally constant (broken line), linear (solid line), and quadratic (dotted line) fits for the Lobster data (using the backfitting method) 36b FIG. 8B: Locally linear fit with the credibility intervals for the Lobster data (using the backfitting method) 36c vii FIG. 9A: A comparison of the locally linear (broken line), 5th order polynomial (solid line) and 10th order polynomial (dotted line) fits for the sinusoid (y = 3cos x — 5sin x + 6) data (using the cross-validation method) 37a. FIG. 9B: A comparison of the locally 5th order polynomial (solid line), 7th order polynomial (broken line) and 9th order polynomial (dotted line) fits for the sinusoid (y = 3cos x — 5s'm x + 6) data (using the backfitting method) 37b FIG. 10: Locally constant fit with the credibility intervals for the difference of the M.A. pH data, (using the backfitting method) 37c viii ACKNOWLEDGEMENTS I would like to thank Dr. Harry Joe for his guidance and assistance in the producing of this thesis. I am indebted to Dr. James Zidek for his helpful suggestions, encouragement, and for his careful reading of this thesis. I also appreciate the assistance of Brian Leroux in translating the article of Collomb (1981). The finanical support of the Department of Statistics of the University of British Columbia is gratefully acknowledged. 1. INTRODUCTION. This thesis is concerned with a Bayesian nonparametric method for estimating an unknown smooth function given some observed data with noise. Consider the model F(x ) = 5 ( x , ) + e,-, t i = l,...,n where S() is an unknown smooth function and the random variables e,- are uncorrelated with zero mean and constant variance a . 2 In recent years, many smoothing methods have been developed to estimate S nonparametrically from the observations on Y(x.), i= l,...,n. A survey of selected nonparametric smoothing methods is presented in section 2. Included are the kernel method, the spline smoothing method, the exponential smoothing method and some robust smoothing methods. Weerahandi and Zidek (1985) propose a Bayesian nonparametric method to tackle such problems. The method is flexible and easy to use. Its derivation is technically elementary and can be extended to the analysis of space-time series and the regression problem. S is assumed to be locally regular, that is expandable in a Taylor series about x = x n + i local regularity is reflected in the size of p. to the pth term. The degree of The objects of inference are an estimate of S(x) and a corresponding credibility interval, where x may lie outside the range of the data. The general theory is described in section 3 and a brief comparison of their approach to those mentioned in section 2 is also included. This thesis is a supplement to the work of Weerahandi and Zidek (ibid). Besides using the cross-validation method to estimate the a priori unspecified hy1 perparameters, a new "backfitting" method is proposed as an alternative. The asymptotic behaviour of the cross-validated sum of squared errors is also examined. p. A data-based procedure is developed to choose the appropriate order Section 3 gives detailed discussions of these proposed methods. Finally, in section 4, all the results are tested on the data of an acid-rain wet-deposition time series. The calculations were performed by using Fortran subroutines which are included in the Appendix. 2 2. SURVEY OF SELECTED NONPARAMETRIC SMOOTHING METHODS. Summary : Suppose we are given observations [ ( x , - , t = l,...,n] satisfying the model Y[xi) = s{x ) + e,-, i = 1,n. { (2.1.1) where the {e,} are independent and identically distributed (i.i.d.) random variables with zero mean and constant variance a . The abscissae {x,} are assumed 2 to be known without error, and to satisfy xi < z 5: ••• 5: x . 2 n Here s() is a fixed but unknown smooth function. We consider the problem of recovering s() when only discrete, noisy measurements of it are available. In recent years, many nonparametric smoothing methods have been developed to tackle this problem. The extensive literature in this area makes a detailed review of each paper infeasible. We aim instead for a brief summary of the principal nonparametric smoothing methods including the kernel method, the spline smoothing method and the exponential smoothing method which is used quite often in the analysis of time-series data. some robust smoothing methods are also considered. Finally, This survey will help to put the developments of section 3 into perspective. 2.1 T h e Kernel ( Moving Window ) Method Watson (1964) and Nadaraya (1964) independently and simultaneously proposed the kernel estimator of s(-) : , _ S[X) ~ Y{XJ)K((X - Xj)/h) E?= K((x-x )/h) 1 3 i where the kernal K(x) sup is a certain density function satisfying the conditions : |-KXz)| < °°> — oo<l<oo and h = h(n) J—oo x-+:too ( the bandwidth of the kernel ), a sequence of positive numbers converging to zero which specifies the degree of smoothness, of the estimator. This is analogous to the kernel density estimator introduced by Rosenblatt (1956). He proposed the kernel estimator, which can be written as where K is the kernel and xi,...,x are assumed to be i.i.d. with the unknown n density /(•). Referring back to the Watson-Nadaraya estimator, Watson (ibid), using some Monte Carlo experiments, concluded that it is h and not K which is important in the estimate. Schuster (1972) obtained asymptotic normality as a correction of a result of the same type given by Nadaraya. A related estimator was introduced by Priestly and Chao (1972) : where K is the kernel, satisfying the conditions : K(u) > K {u) 0, for all u dvu 2 < oo, and K(u) — oo Stone (1977) considered a general form of estimator, n t'=l 4 — K( — u). where v} (z) — tu ,(x, xi,x ), nt n 1 < i < n , depends more heavily on those x,- n close to x. This formulation includes the Watson-Nadaraya and the Priestly-Chao estimators as special cases. Stone (ibid) also derived the sufficient conditions on {w } for the consistency of s(x). n Clark (1977) proposed a convolution-based estimator, s(i) =/_>'M^) "' < where { for x < Fi Y xi, for x > x . n n He claimed that such estimators have some advantages over the Watson-Nadaraya and the Priestly-Chao estimators in small samples. Firstly, s(x) is smoother. Sec- ondly, the Watson-Nadaraya and the Priestly-Chao estimators are not necessary invariant under linear transformations and they do not have the shape-preserving properties of s(x). It is well known [ see Table 1 of Rosenblatt (1971) ] that the choice of the kernel function K, is of relatively unimportance. computational convenience. Usually, a kernel is chosen for While the choice of kernel does not seem to be an issue, at least for reasonably moderate sample sizes, the choice of bandwidth is quite a different matter. In practice, the estimate is often formed by computing the estimate for several h's and selecting the estimate with the most pleasing appearance. H the h is too small, the estimate will appear to be too wiggly and if too large, the estimate will begin to be too flat. 5 h from the data by a minor modification of the Clark (1977) estimated cross-validation technique described by Wahba and Wold (1975). A simulation experiment demonstrated that this approach was satisfactory in determining the appropriate degree of smoothness. Similarly, Hardle and Marron (1985) estimated h of the Watson-Nadaraya estimator by minimizing n CV{h) = n- 52(Y{x )-i (z ))* l i i i t'=i where s,- is the cross-validated estimator given by , S l { x ) E y(«jW(«-gy)A) w - T ^ K{{x-x )lh) Ji j • j They also showed that such an h is asymptotically optimal with respect to the distances : n MSE ISE (Averaged Squared Error) = n~ ^^(s(£,) — s(z,)) , 2 l (Integrated Squared Error) = j (s(x) — s(x)) f(x) dx 2 where f(x) is the marginal density of x, and CMISE (Conditional Mean Integrated Squared Error) = E[ ISE \ xi,...,x n ]. We may be interested in constructing confidence intervals for s(x) or confidence bands for s(). However, besides some theoretical results mentioned in Bickel and Rosenblatt (1973) and Collomb (1981), no practical method has been reported. 2.2 T h e Spline S m o o t h i n g M e t h o d 6 For the model in equation (2.1.1), the spline smoothing estimate s of s is a defined to be the curve that minimizes (2.2.1) over the class of all functions s for which s and s' are absolutely continuous and s" is square-integrable. The idea of (2.2.1) is to trade off fidelity to the data (small sum of the squared residuals) against smoothness (low value of integrated square second derivative). The smoothing parameter a determines the rate of exchange between these quantities and hence determines the amount by which the data are smoothed to produce the estimate s . As a —> oo, s a a becomes linear, and the limiting function SQQ is the least squares straight line passing through the data. As a —• 0 replicated, until, in the limit, SQ passes through the data; So is the natural cubic spline interpolator of the data which is discussed below. The idea underlying this estimate can be traced back to Whittaker (1923), who proposed a discrete version using finite differences in place of derivatives in (2.2.1). Schoenberg (1964) proposed the replacement of the finite differences in Whittaker's objective function by derivatives. Reinsch (1967) showed that s is a cubic spline with the following properties a : (i) it is a cubic polynomial in each interval (z,-, z, i), + (ii) at each x,-, the curve and its first two derivatives are continuous, but there may be a discontinuity in the third derivative, (iii) in each of the ranges (—co,xi) and (x„,co) the second derivative is zero, so that s is linear outside the range of the data. a The method of cross-validation ( see, for example, Stone (1974) ) has been suggested for choosing the smoothing parameter a from the data. 7 Let s~'(x ) be t the solution to problem (2.2.1) with the ith data point , (x,-, F(x,)), omitted. Of course, s~*(xi) may be regarded as an estimator of Y(x{) and a as appropriately chosen if s~'(x ) is a good estimator of Y(xi). To measure the goodness of fit, t we choose the average squared error, that is , n XVSC(af=1n- 1 J^x,-) - C(*.)) (2.2.2) 2 t-i which is called the cross-validation score. s (xi) = a Let A(a) be a matrix such that Y^Aij{<x)Y{xj). 3 Craven and Wahba (1979) showed that (2.2.2) has the easier computational form XVSC(a) = n - ' ± ^ ' ) - S f ~ f . (2.2.3) They also suggested the use of a related criterion, called generalized crossvalidation, obtain from (2.2.3) by replacing Aa(a) by its average value, n traceA(a), say, n~ tvA(a). -1 1 GXVSC{a) where RRS(a) This gives the score = n- RSS{a)/(l l - n^tr^a)) , (2.2.4) 2 is the residual sum of squares J27=i0^( i)~ a(z»)) x s 2 Craven and Wahba gave fuller details and justification of the preceding derivation. Their arguments suggest that generalized cross-validation should yield approximately the optimal value for a, in the sense of minimizing the average squared errors n _ 1 J27=i ( a ( « ) — ( » ' ) ) s x s x 2 ( A Fortran implementation can be found in the IMSL Library (1980). ) To use the formula (2.2.4) one still needs to find the trA(a). (1984) proposed a way to derive such an approximation. Silverman Let f(x) be an esti- mate of the local density of z,-, calculated using the fast algorithm of Silverman 8 (1982) on the range [a, 6] just containing the design points. Let J a It can then be shown that n « 2 + ^ ( l + c a(t-1.5) ) . trA{a) 4 0 (2.2.5) - 1 t'=3 Substituting (2.2.5) into (2.2.4) gives a score function called the asymptotic generalized cross-validation (AGXV) score which can then be minimized to give an automatic choice of smoothing parameter. A simulation study is carried out in Silverman (1984) and the practical performance of AGXV turns out to be even better then that of generalized cross-validation in that, for non-uniformly spaced data, the proportion of cases giving a bad value of the smoothing parameter is substantially reduced. In constructing confidence intervals, Wahba (1983) exploited a Bayesian interpretation of the smoothing spline, that is a Bayes estimate with respect to a certain zero mean Gaussian prior. for the s(x,). She was able to develop credibility intervals However, it should be emphasized that these are not simultaneous credibility interval estimates. Some extensions and modifications of Wahba's work can be found in Silverman (1985). 2.3 T h e Exponential Smoothing Method Exponential smoothing method was developed in the early 1950s. then it has become a popular method of forecasting. formulations are relatively simple. tional effort are required. Since It is because the model Also, only limited data storage and computa- Perhaps the most important reason for the popularity 9 of exponential smoothing is the surprising accuracy that can be obtained with minimal effort in model identification. Makridakis et al. (1982) reported that exponential smoothing method proved to be at least as accurate for forecasting as more complex approaches such as the Box-Jenkins autoregressive integrated moving average (ARIMA) technique. In essence exponential smoothing method applies an unequal set of weights to past data. These weights decay in an exponential manner from the most recent data value to the most distant value. It is also assumed that the {x,} are equally spaced and that F(x ) = F, = Y is a time series. t t The main use is to forecast F i given Fi,...,F although this method can also be used to smooth t + Yi,...,Y n t ( for example, see Fig. 1 in Brown and Meyer (1961) ). However, it should be pointed out that this method depends just on the previous data. So as a smoothing method, it is probably not as good as the kernel and the spline smoothing methods which use the data in a neighbourhood of a point. Here are some commonly used model formulations: (i) Simple exponential smoothing (no trend in the series) The general equation for the forecast value can be written as F t + l = aY + {1 - a)F t (2.3.1) t where Ft is the forecast at time t, Yt is the observed value of the series at time t and a is a smoothing constant between 0 and 1 A large value of a, say .9, gives little smoothing in the forecast, where a small value of a, say .1, gives considerable smoothing. Most often a is chosen by trial-and-error. optimal value of a may be chosen by minimizing the MSE 10 An of the forecast. The initial value (as of time 0), Fo, can be estimated simply using the mean of the data. Alternatively, backcasting is suggested to obtain FQ. It is done by inverting the data series and starting the estimation procedure from the latest (most recent) value and finishing with the first (oldest) value. Doing this will provide forecasts or parameter estimates for the beginning of the data, and these can be used as initial values when the data are forecast in the usual sequence, that is, from the beginning to the end. Trigg and Leach (1967) suggest an adaptive response rate approach for the automatic specification of a. It allows a to vary from one period to another so that the forecast can respond to the pattern of the data as closely as possible. The adaptive response rate approach is based on (2.3.1), but a varies and is calculated according to a = \E /M \ t+l t t where E t ( smooth error ) = f5e + (1 - (3)E -i, t t Mt ( smooth absolute error ) = 0\e \ + (1 — (3)Mt-i, t e = Y - F , and t t t 0 is a constant usually set at 0.1 or 0.2. We can let the initial values E\ = 0 and Mi equal an estimate of the mean absolute deviation, that is, The value of a will fluctuate between 0 ( perfect forecasts ) and 1 ( extreme t bias in the forecasts ). 11 (ii) Winters' (1960) linear exponential smoothing (linear trend in the series) A smoothed estimate of the trend is given by T = 0(S - 5 _i) + (1 - P)T t t t t (2.3.2) U where St is the equivalent of the simple exponential smoothed value and (3 is a smoothing constant, analogous to a. obtained by backcasting. squares. The initial values So and T can be Q Alternatively, they are found by using ordinary least It is done by solving the equation for a straight line to obtain the intercept and the slope, and using these as initial parameter values. The basic principle underlying (2.3.2) is the same as that of (2.3.1). The most recent trend, (S — S -i), t T -i, t t is weighted by (1 — 6). is weighted by B and the last smoothed trend, Then this smoothed trend is combined with the standard smoothing equation to obtain S = aY + (l-a)(St- +T _ ). t t l t (2.3.3) 1 In order to use these smoothed series values, S, and the smoothed trend component, T , to prepare a forecast, the trend component must be added to the basic smoothed value, for the number of periods m ahead to be forecast. The forecast is given by Ft+m = S + mT . t (iii) Winters' t (1960) linear and seasonal exponential smoothing (linear trend and seasonality in the series) This model is based on three equations, each of which smooths 12 a factor associated with one of the three components of the pattern - randomness, linear, and seasonal. These three equations are as follows : S = otYt/h-L t + (1 - <*){St-i + T = B(St-St-i) t I = Y /S t 1 t t + T -!), t + {l-B)T - , t and l (l- )I _ , 1 t L where 5 is a smoothed value of the deseasonalized series, T is a smoothed value of the trend, / is a smoothed value of the seasonal factor, and L is the length of seasonality (e.g. number of months or quarters in a year). The initial value of the smoothed seasonal factor / can be found by the classical time series decomposition methods (for example, see Makridakis and Wheelwright (1978).) The equation for / is comparable to a seasonality index. membered that Yt includes any randomness in the series. It must be re- In order to smooth this randomness, the equation for / weights the newly computed seasonal factor (Yt/St) with 7 and the most recent seasonal number corresponding to the same season, It-L> with (1 — 7). The equation for Tt is given in (2.3.2). equation for 5 , Y is divided by the seasonal factor h-Lt t In the This is done to deseasonalize ( eliminate the seasonal fluctuations of ) F . h-L is used because t It cannot be calculated until S is known. The forecast is given by t Ft+m = {St + mT )It-L t + m- For models (ii) and (iii), the best values of a, 3 and 7 are found by trying various combinations of values and choosing that set of values that minimizes MSE of the forecast. Only exponential smoothing models with simple patterns are described in the 13 above discussion. Smoothing models for other complicated patterns are desribed by Makridakis and Wheelwright (1978), Makridakis et al. (1982) and Gardner (1985). 2.4 Some robust smoothing methods For the model in (2.1.1), Cleveland (1979) proposed the use of a robust weighted regression procedure for smoothing scatterplots. Firstly, for each i compute the estimates, /?y(x,), j = 0, l,...,d of the parameters in a polynomial regression of degree d of Y(xk) on Xk, which minimize n ] T w {xi)(Y{x ) k k - 0 o - 0 i x k ~ ... - 6 x d k d ) 2 k-l where w (xi) k = W ( / i , ( x f c — x,)) and hi is the rth smallest number among _1 \xi — xy|, for j = l , . . . , n . W is some given weight function and r = fn where / is the parameter used to determine the amount of smoothing. Thus, F(x,), the fitted value of the regression at x,-, is X2y /3j(x,)x,~'. =0 Let 6 k be some robustness weights. Then compute a new K(x,) for each t by fitting a dth. degree polynomial using weighted least squares with weight 5fctUfc(x ) at (xfc,Y"(xfc)). t Now repeat this procedure successively to compute new { by the procedure indicated below, for example ) and Y[xi), final Y(xi) are robust locally weighted regression fitted values. t times; the Lastly, linear interpolation of these fitted values are used to obtain the robust smooth curve. Cleveland chose 6 = i?(efc/6s), where B, the bisquare weight function is dek fined by 14 B(x) " [ X ) tk = Y(xk) I < " ^ ~ \0 I ' < *' for |z| > 1, X F O R — Y(xk), and s is the median of the for choosing /, d, t and W in practice. 1 He also gave some guidelines ( Cleveland's procedure is available in S [ Becker and Chambers (1984) ] through the command "lowess". ) Friedman (1984) constructed a variable span smoother based on locally linear fits. It is very fast to compute and the value of the parameter that controls the amount of smoothing is automatically optimized locally ( through cross-validation ) , allowing it to adapt to the response over the range of abscissa values. In the context of robust spline smoothing, we can replace the sum of squared errors in the objective function (2.2.1) by a different function of the errors, to give (2.4.1) Here p is a convex and less rapidly increasing function with bounded derivative. Huber (1979) chose for \x\ < c, for \x\ > c. The constant c regulates the degree of robustness and good choices for c are between 1 or 2 times the standard Y(xi). deviation of the individual observation Lenth (1977) and Huber (1979) use iterative reweighted procedures to mininize (2.4.1). least squares Utreras (1981) gives a numerical procedure based on the Newton minimization algorithm for computing the minimizer of (2.4.1). 15 2.5 T h e Multivariate Extension What we have discussed so far are the nonparametric smoothing methods of an univariate function s ( ) . If s ( ) is now a function of a vector variable x, then the extension of the above smoothing methods is straightforward. Here is a brief description: (i) the kernel method Using the Watson-Nadaraya estimator as an example, now h = € (hi,h ) p R+, the collection of p-tuples of positive numbers; x = (ti,...,t ) p the kernel i f is a certain density function on R . p G R p and We can use the approach of Hardle and Marron (1985) described in section 2.1 to estimate h. Bierens (1983) proposed a new class of kernel estimator which is constructed through a Sample Moments Integrating Normal Kernel (SMINK) estimator. SMINK The estimator is itself a density and satisfies the condition that it's first and second moment integrals equal the first and second sample moments, respectively. For p = 2, the kernel estimator of s() is given by „, .. ., f+"Yfi(Y,x\hi)dy /2W12) where fi(Y,x\hi) and ^ ( x l ^ ) are SMINK estimators of the density f(Y, x) and / ( z ) , respectively, and ^ 1 and h? are kernel bandwidth parameters. 2 procedure is used to estimate hi and A tv/o-step and satisfactory performance of s(x) is revealed in a numerial example. 16 (ii) the spline smoothing For variable the model x, where method in (2.1.1), x = s is now a function of a p-dimensional vector T h e smoothness penalty function (ti,...,t ). p f in s"(u) du 2 (2.2.1) is replaced by In the particular case, p = m = 2, the smoothness penalty function becomes T , . f°° 2d s f°° ,d s 2 d s, 2 2 M J x We wish to find the multivariate spline smoothing estimate, s , of or which min- a imizes 1 ~ n J2(X(xi) - s(xi)) (2.5.1) + aJ (s). 2 m For p = m = 2, (2.5.1) becomes ^E(H^)-s(x,)) + aJ ( ). 2 (2.5.2) 2 S t= l T h e minimizer of (2.5.2) is one of the "thin plate splines", J (s) so called because is the bending energy of a thin plate. 2 W a h b a and Wendelberger generalized cross-validation (1980) and Wahba (1980) have shown how to use ( see Craven and Wahba (1979) ) to choose oc from noisy data, and have given a computational has also given computational the s a which minimizes procedures. (2.5.1), Wahba (1984a,b) do provide solutions. 17 algorithm for p — 2. In the general and Wendelberger case, Utreras (1979) i.e., estimating (1980) and Wahba It may be noteworthy that in the mining industry, a class of two dimensional smoothing methods known as "kriging" are closely related to the thin-plate splines. Comparisons of these two smoothing methods are found in several sources. Interested readers are referred to Dubrule (1983,1984) and Watson (1984). (iii) the exponential smoothing method Simple exponential smoothing method has been generalized to multivariate forcasting by Jones (1966) and Enns et al. (1982). The generalization is straightforward. We replace the scalars in equation (2.3.1) with: F t+1 = aY + (1 - a)F . t t If there are p series, the dimensions of F F + i and Y are px 1. The smoothing t) t t matrix a has dimension p x p. Jones (ibid) proposes a recursive algorithm to estimate the smoothing matrix a. Enns et al. (ibid) introduce a class of multiple exponential smoothing models as an alternative to simple univariate exponential smoothing and Trigg and Leach (1967) suggest adaptive models. Such models are adaptive in that the smoothing matrix a, which is a generalization of the smoothing constant of univariate models, changes from period-to-period. The model parameters, including the full variance-covariance structure of the multiple time series as well as the smoothing matrix a are estimated by maximum likelihood. Harvey (1984) further simplified the computations of Enns et al. (ibid) so that one can use univariate smoothing models to forecast related time series. More- over, Harvey's results hold for smoothing models containing polynomial trends 18 and seasonal components as well. (iv) the robust smoothing methods To our knowledge, no generalization to the multivariate analysis has been reported up to now. 2.6 Discussion It should be emphasized that in the above survey, we have been guided in our choice of material by concerns related to the applicability of the methodology rather than theoretical properties of different kinds of smoothers. For discussion related to the latter and comprehensive lists of references on this subject, please refer to Stone (1977), Wegman (1980), Collomb (1981), Wegman and Wright (1983) and Silverman (1985). From this survey, we can discover some shortcomings of the described smoothing methods. For instance, the spline smoothing and the robust smoothing methods are used for smoothing rather than forecasting. The exponential smoothing method is mainly used for forecasting but not for smoothing. Also, explicit form of confidence intervals of the smooth curve estimate are hard to derive. How- ever, Weerahandi'and Zidek (1985) provide an alternative approach to solve these problems. The general theory is presented in details in the following section. 19 3. THEORY OF BAYESIAN NONPARAMETRIC SMOOTHING LOCALLY REGULAR 3.1 APPROACH FOR PROCESSES. T h e General Theory Weerahandi and Zidek (1985) ( abbreviated hereafter by WZ ) propose a Bayesian nonparametric approach for the smoothing of stochastic processes. The stochastic processes of concern are of the form R = S + N where S is a smooth function and N is an independent noise process. In this thesis we take the domain to be the real line and may think of the stochastic process as a time series. Assume R has been observed at a sequence of time-values, ti, and let r,- = R(U) i = l , . . . , n . The objects of inference are an estimate, B = S(t i), 0 n+ the interpolant or predictand, and a corresponding credibility interval. S is assumed to be locally regular, that is expandable in a Taylor series to the pth term about t = t i. A local parametrization of S is thereby achieved where the n+ parameters are Bi = D S{t ), l n+l D denoting the differentiation operator. An a priori structural model for the data is (3.1.1). R = X/3 + e Here R = {R ...,R ) , Ri = R{U), 3 = {B , ...,B ) T u n Q trix, is given by X = (1, X\,X ), p [*»-<»+i] //0 tn i} '/j. + y 3 7Y(i,-), r} = (r/i,rj ), n joining and t n + T i X, an n x (p + 1) ma- 1 is an n-vector of l's, Xj J = 1,-,P, e = 7? + /Y, N T = (N ...,N ), u n T = ([ti — Ni = and rn is the remainder of the Taylor expansion of £(£,), T that is tji = [U - t p n + l i. ] p + i D p + 1 S ( 9 i ) / ( p + 1)! where 0,- is a point in the interval The degree of local regularity is reflected in the size of p. Although R is observed, 3 in equation (3.1.1) and e are not. 20 The sec- ond assumption i n W Z ' s approach is that the expansion errors a n d all other a priori uncertainity about tion. R, B a n d e have a joint multivariate n o r m a l distribu- Furthermore, the 17,'s are assumed to be independent w i t h mean zero a n d the variances of the 77,'s are of the form uncertainity about the size of D p + 1 of 5 \ti 2 — £ +i| 2p+2 n where S 2 represents S/(p + 1)! i n the remainder term of the Taylor expansion of S. In order to estimate Bo = S ( t equation (3.1.1) holds with n + l Cov(e) = a H, x Cov(p) results = S. where <r = Var(JV,) is the variance 2 of the noise, H = d i a g { ( l + c\t - t and ) , the interpolant or predictand, assume that n + i\ 2 2 p + 2 ) , ( 1 + c\t - i n n + 1 | 2 p + 2 )}, c = <5 /<r , 2 2 T h e n the inference problem can be solved by using standard ( see Lindley and Smith distribution of U given V. J R | S , a ,H (1972) Let "U\V" denote the conditional Then ~ N{XP°, 2 ). XY,X + a H) T 2 and p\R,Z,o- ,H~ 2 N(F(12~ P 1 + 0 a~ X H- R),F) 2 T 1 where F = (S- +a- X H~ X)- . 1 2 T 1 1 T h e n using a diffuse prior for P, i.e., letting E (X H~ X) T l -1 X H~ R, T is noteworthy degree of smoothing. 1 —• 0, E(P\R) converges to the generalized least squares estimator of p with a pos- l teriori covariance matrix It - that cr (X H~ X)~ . 2 c serves T l x as the parameter controlling the appropriate If c is too small, f will appear to be too smooth, a n d if 21 too large, f will be too rough. Stone (1974) ). One approach to estimate c is cross-validation ( Successively, t +i = U is chosen, and then r is dropped from n t the sample and interpolated using r,- = 0o{i), obtained by the above generalized least squares procedure. Lastly, c = c is chosen to be the value which minimizes the cross-validated sum of squared errors, i Alternatively, maximum likelihood estimation may be used to find c. However, in the acid rain example ( see section 4 ), c = 0 was obtained, so this method is ineffective and hence rejected. A "backfitting" method decribed below can also be used to estimate c. It is shown to be computationally efficient and it coheres to the general theory. A detailed description of this method is given in section 3.2 (II). For t n+i = ti, the a priori local structural model is Ri = a,+<7Z,-, i = l,...,n where ct,- = S(ti) and Zi has mean 0, variance 1. Given c = c, or,- and Var(or,) = c,tr can be estimated by the above generalized least squares procedure. 2 E(Ri — a,) 2 = Var(dj) +cr =<7 (l + c,) when the noise is regarded as independent of the smooth. 2 2 From thi3, an unbiased estimated estimator of cr is given by 2 a — n~ ^2i (R. — <*»')/(l + t ) 2 l Thus 2 c Consequently, since a posteriori & has a normal distribution with mean a, the approximate 95% credibility intervals for a are given by d ± 1.965.JE'.(a:). The above methods were tested on the data from an acid rain, wet-deposition time series and very encouraging results have been obtained. In the case con- sidered in WZ, only locally constant, linear and quadratic (i.e., p < 2) structural models are considered for illustration. 22 In general, the WZ-approach has some advantages over the other nonparametric smoothing methods decribed in section 2. Firstly, the methods mentioned in section 2 share the shortcoming that the explicit form of the estimates of s(), especially for extrapolants, and confidence bands are hard to derive. For exam- ple, Silverman (1985) points out that in spline smoothing only implicit solution can be obtained except for large samples when an approximate explicit expression may be derived. In contrast, the WZ-approach does yield explicit estimates of s() with corresponding explicit credibility bands. Secondly, the prior that Wahba (1983) uses is essentially unique so that the applicability of the spline smoothing is severely limited. The exponential smoothing method suffers from the lack of an objective procedure for model identification. In contrast, in the WZ-approach, every aspect of the analysis is guided by the local-structural model and the Bayesian paradigm itself. Prior knowledge may be incorporated directly into the analysis and "borrowing from strength" is feasible. The "backfitting" method for estimating the smoothing parameter c is simple and computationally fast compared with other procedures. Furthermore, the WZ-approach can be extended to the analysis of multiple, multivariate spacetime series and the regression problem. Weerahandi and Zidek (1986) have successfully extended the WZ-approach to the analysis of multiple independent time series. generated. In the latter case, exact credibility bands for growth curves are Their work on spatio-temporal analysis is in progress and some promising unpublished preliminary results have been obtained. 23 3.2 Estimating the "smoothing parameter" c for locally regular struc- tural models. (I) Cross-validation method In the case considered in WZ, c is found by minimizing XV(c) which depends on a 2 only through c. The results of WZ suggest that quite satisfactory esti- mates of c are obtained in the locally constant, linear and quadratic structural models. So, it is tempting to try the same approach to estimate c = c for the other pth order locally regular structural models when S is assumed to have p-f 1 derivatives. It is then necessary to find out if c is really the global minimum, rather than a local minimum of the function XV(c). Two ways which help to determine this are: (i) Plot XV(c) against c and examine the shape of XV'(c) for large c. case, we need to know the asymptotic behaviour of XV(c) In this as c-> co. Let h < ... < t, n t\ < ... < t , n d;(i)( ) = c XV(oo) = Xj = tj — ti for all j = l , . . . , n excluding j = t , the cross-validated estimate of value of an accent on YI w XV(c) at or at ti, c = oo, and ^ denote summation over all j — l , . . . , n excluding j — i. In the case of a locally constant fit ( i.e. p = 0 ), oc^(c) can be obtained by 24 the generalized least squares procedure described in section 3.1. It is given by: E'd + cx?)- • K , 0 ( C )= It follows that , , (i){°) ~* E t3 j r —2 x - a def . , (t)(°o) as a 2 C -> OO. Therefore XF(oo)«X;Cr -d (oo)) 2 t (l) «=i ^ ' Similarly, in the case of a locally linearfit( i.e., p = 1 ), let ^ ( 1 + cxJ), j 1 then <*(«) () c It follows that a ( )^ , cj T! J T! ~ x A x 2 - ( E ' Z t 3 ) 2 — = a s c->oo. Therefore X V (bo) « ^-i - a (co)^ (t) By using the same technique, we can find XV(co) for higher order locally regular structural models. 25 (ii) Examine the asymptotic behaviour of XV(c) for large c. If XV(c) is less than XV(oo) for large c, then we will have increased confidence that there is not a local minimum at infinity. In the locally constant case: r.--a <)(c) = r , - ( for large E'Cc-^x-d-c-ixT )) 2 c = rv E , for large — ' Let 1 E - = 53 ' 2 Zy 53 ' ^ ^ i 2 = 6 »' y^>^.—4 'V /* 1 = a n d v^/ -2 = «A . *y e Then ^ - a(.)(c) = r, - a,(6- - —)(1 + —) t r,- — a,-6,- — -a,(6,e,- — a*,-) c r- - or(,)(oo) + — c t where /,- = -a,-(&,e,- - dij. 26 for large c c. So XV (c) = £(•>••-«(»•) ( )) c 2 (3.2.1). X:((r,-a(oo)) + f ) £)(r- - a-)(oo)) + -/,(r,- - a-)(oo)) )t( 2 t f( f( Differentiating (3.2.1) with respect to c, we obtain for large c. Similarly, by using the same technique, we can estimate XV (c)' for higher order locally regular structural models. (II) "Backfitting" method Recall from equation (3.1.1), that if S has p + 1 derivatives, then the a priori model is of the form Ri = Po + Six, + B x 2 where B = S{t ), 0 n+1 B = t />S(i n+1 2 + ... + /3 x + 3 p p p + 1 ) , x = (ti-t y/i\, p + i i = n+i { x + e (3.2.2) 1, ...,p + 1 and e = t] + N as defined in section 3.1. Recall that <5 represents the uncertainty about the size of D S/(p+1)1 2 in the p+1 remainder term of the Taylor expansion of S i.e. 6 is the a priori variance of 2 D S/(p p+l + 1)!. Hence t +i = ti may be chosen successively , and <5 may be 2 n 27 estimated by 8 2 = sample variance of (/? sample variance of ((P+1)!) = p+1 /(p + 1)!) (0 +i) p ' 2 We also note that c = <5 /<r , so the value of c for the order p, c , say, is given 2 2 p by b a 2 2 _ sample variance of (P +i) P er ((p+l)!) = 2 2 ( 3 ("-i)- Efai(4 i(0-jgp i) i = + * ((p+l)!) a where /? i(z) is the coefficient of x p p+ a 2 ' 2 3 ) a + a , 0 +i is the sample mean of / ? p p+1 ' s , and can be found by the generalized least squares procedure described in section 3.1. Thus, if (3.2.2) holds, by the same argument, we can use (3.2.3) to estimate the values of cy for any order j = 0, ...,p. However, to use the formula (3.2.3) we need c p + i- validatory approach mentioned in section 3.1 to estimate c obtain estimates of /?,-, <5 and cr . 2 assuming p is large. the size of D S/(p p+1 large p. 2 We may use the crossp + i , from which we can Alternatively, we may simply put c p + 1 = 0 The reason is that at a certain stage, the uncertainty about + 1)! will be reduced approximately to zero for sufficiently So, it yields the result 6 —* 0, hence c 2 In summary, if (3.2.2) holds and c i p+ the lower orders j < p by using (3.2.3). "backfitting" procedure. 28 p + i —• 0. is fixed, we may estimate the cy of We refer to such an approach as the (Ill) Conclusions Comparing cross-validation efficient methods for the (I) (II), following since we do not estimate and we prefer the reasons. need to Firstly, it compute coheres resulting curves with are the theory c p + i is proposed i n section comparatively smoother. In to is computationally more f,- for each p and c, cy for all j < p within one run once approach "backfitting" approach fixed. 3.1 and we Secondly, this and we expect consequence, can the the smoothers obtained for different orders of local regularity become comparable. 3.3 Choosing the "optimal" order p. For c p many practical purposes, it is sufficient to choose p subjectively. is determined, say by either of the methods decribed in section 3.2, we may plot the resulting curves and then choose the one which "looks best". in p. Once many contexts it is desirable to Hence, one of our main value of p in the sense that have an automatic objectives in this thesis f, estimated minimizes some measure of accuracy. error PSE procedure for choosing is to find the " o p t i m a l " by the pth order locally regular fit, In our case, we use the predictive squared (to be defined later). A fit giving a m i n i m a l PSE prediction However, yields the model of the data which gives the best (in the least square sense) of each data point. m i n i m u m PSE gives the parameters In other words, the which maximize the internal consistency of the data set with respect to the estimated model. 1. F i x the span of observations m ( m < n 2. F i x the value of p 4-1 and set c p+i — 0. 29 ). O u r procedure is as follows: 3. Use the first m observations and apply the "backfitting" procedure described in section 3.2 to estimate 4. cy, , for j = 0,...,p. = f i(cy) for each j < p. m + 1 m+ Shift to the next m observations i.e. r 2 , r estimate f 6. ... CI, Once cy is found, apply the generalized least squares procedure described in section 3.1 to estimate r 5. CQ, m+2 m + i and repeat steps 3 and 4 to ( c y ) for each ; < p. Repeat step 5 for successive m observations until r (cy) for each j < p is n found. 7. Add the sum of squared deviations from steps 4 to 6 for each j < p. This is the predictive squared error, . PSE(J) n—m ^ £ n — m 8. Compare the PSE(j) minimizes the (r m+i - f m+t (cy)) , 2 ; = 0,p. — for each j < p to find the optimal j < p which PSE. In the procedure decribed above, we are predicting one observation ahead each time. It is perhaps more natural to predict the the average of w observations ahead since S(t i), n+ the quantity of interest, is the smooth, i.e. stripped of noise at time t = t +\. Our intention can be easily realized by a n slight modification of the above procedure. (i) In step 4, use the cy to find f the process m + 1 It is done as follows: (cy),r m + U ) ( c y ) for each j < p. Let 1 «'(cy) = — $3 Wfc+i-i(cy). fc=i w f (ii) In step 6, repeat step 5 for successive m observations until 7v(cy) [w' = 30 n — m — w + 1), for each j < p are found. (iii) In step 7, for each < p, replace PSE(j) by 1 w i=i where 1 r, 1• w to M k=l Finally, follow the same steps 1 to 8 in the original procedure with the above modifications. In section 4 , both types of predictions are considered for illus- trations. 31 4. APPLICATION. 4.1 Introduction In this section, the methods described in the last section are tested on the data used in WZ. The data were obtained from one of the nine stations of the MAP3S/PCN monitoring network which is found in the Northeastern part of the United States and has provided event-based chemical measurements of v/et deposition events since about 1976. A detailed description of the data set can be found in Gentleman, Zidek and Olsen (1985). our analyses are confined to two subcases: For illustrative purposes, (i) the field pH values which were measured at a station located in Pennsylvania State University during the precipitation events of 1977 ( 80 observations were recorded ), and (ii) the monthly average pH-values over the 1976-82 period at the same station ( 71 observations were recorded ). M.A. pH data. Hereafter (i) is referred to as the pH77 data and (ii) as the Throughout this section time is measured in days starting with January 1, 1976 which is regarded as the origin. Recall that in equation (3.1.1), c = 5 /cr . 2 P / ? i / ( p + 1)!, is extremely small. p + 10~ , 2 So u;,- = t = l , . . . , n , in estimating c p 2 However, 5 , the variance of 2 — t i n+ numerically. obtained are enlarged by a factor of 1 0 4j+4 is rescaled by a factor In this way, all the for all j = 0, ...,p. {cy} The results reported below are obtained by the methods described in section 3; all the calculations were performed by using Fortran subroutines which are included in the Appendix. 4.2 Comparison of the cross-validation and the backfitting methods. 32 Before making any comparisons of these two methods, we examined the asymptotic behaviour of XV(c) for large c, to confirm that the cross-validatory estimate of c p being found is really the global minimum of XV(c). Only the locally constant fit is considered as an illustration. Observe in Table I that the XV(c) values increase as c tends to infinity for both the pH77 and the M.A. pH data. Also, for large c, XV(c)' « 0.2738c" and 0.000372c 2 the M.A. pH data respectively. of XV(c), -2 Therefore, we are sure that the global minimum c , is finite. p In applications, c p , the cross-validatory value of c , is found by using an p IMSL (1982) optimization subroutine to minimize XV(c). of c P for the pH77 and The estimated values for the pH77 data are presented in Table II and the corresponding locally polynomial fits are displayed in Figure l(A-C). Figure 1 shows that the locally quadratic fit is smoother than the other polynomial fits in the sense that it cannot make the quick turns so it tends to round off the corners of the dips of the data. This may due to the fact that c is comparatively small so a greater 2 degree of smoothness is revealed in the resulting curve. As we have pointed out in section 3.1 6 is decreasing as p is increasing, so c should be decreasing 2 too. p However, as p increases, c become relatively large and the resulting curves p look quite wiggly. Nevertheless, cio = 0 is obtained. It is noteworthy that unlike what happens in regression, the fit here is found to be fairly robust with respect to the degree of the polynomial within the range of the data. When Cio = 0, we can apply the backfitting procedure descibed in section 3.2 to obtain c p , the backfitted values of c , p p = 0, ...,9. The results are summa- rized in Table n and the corresponding locally polynomial fits are portrayed in 33 TABLE I : Estimated v a l u e s of XV(c) f o r a l o c a l l y constant f i t to the pH77 and M.A. pH data. c XV(c) (pH77) XV(c) (M.A. pH) 0 6.0485 3.5161 3.2 * 4.5944 3.0132 3.8339 3.2376 1 . OE+3 4.0742 3.2611 1.OE+4 4.5993 3.2642 1.0E+5 4.7845 3.2644 1.0E+6 4.8083 3.2646 1.OE+7 4.8108 3.2647 1.OE+8 4.8110 3.2648 1.0E+9 4.8111 3.2648 1.0E+10 4.8111 3.2648 4.8112 3.2648 128.0 o o * ** ** c r o s s - v a l i d a t o r y estimate of c f o r the M.A. pH data c r o s s - v a l i d a t o r y estimate of c f o r the pH77 data 33a TABLE I I : E s t i m a t e d Order A * Cp (p) v a l u e s o f Cp f o r t h e pH77 d a t a . A rt Cp * Cp ** 0 128 45 1 1750 4034 2 1028 75721 3 350000 279977 4 1715000 452421 5 2486000 171752 6 616000 471 19 7 46300 1 989 8 4700 1 07 9 400 .002 10 0 : Cp o b t a i n e d by t h e c r o s s - v a l i d a t i o n * Cp : Cp o b t a i n e d by t h e b a c k f i t t i n g 33b method method FIG. B, p=2 p=3 p=4 1 : A comparison of the l o c a l l y polynomial f i t s of d i f f e r e n t order p t o the pH77 data, ( u s i n g the c r o s s - v a l i d a t i o n method) (solid line) (broken l i n e ) (dotted l i n e ) „ p=4 ( s o l i d l i n e ) p=7 (broken l i n e ) p=1u(dotted l i n e ) 300.0 350.0 «a.o 33c 450.0 TIME 500 0 (DATI 550.0 600.0 650.0 700.0 800.0 FIG. 2 : A comparison of the l o c a l l y polynomial f i t s of d i f f e r e n t order p to the p H 7 7 data, (using the b a c k f i t t i n g method) B. p=2 p=3 p=4 C. (solid line) (broken l i n e ) (dotted l i n e ) p=4 ( s o l i d l i n e ) p=7 (broken l i n e ) p=9 ( d o t t e d l i n e ) 300.0 350.0 4)0.0 33d 450.0 TIME 500.0 lOflTI 550.0 600.0 700.0 750.0 8 0 0 0 Figure 2(A-C). In general, the estimated values of c p than the corresponding c fits. are comparatively smaller except in the cases of locally linear and quadratic p The locally constant and quadratic fits obtained by these two methods are compared in Figure 3. difference between c p The resulting curves do not change drastically unless the and c p is greater than a factor of 10. As a final note about Figures 1 and 2, the locally polynomial fits emerge from the data at either extremity with a marked upward or downward trend as the order p increases. This is inevitable; as WZ point out , for large the a,'s become the least squares polynomial estimates. For for the M.A. pH data, the estimated values of c p and c p are given in Table III and the corresponding locally polynomial fits are displayed in Figure 4(A-C) and 5(A-C). Since the estimated c 's are relatively smaller, the correp sponding locally polynomial fits obtained by the backfitting method are comparatively smoother than those obtained by the cross-validation method. In Figure 6 are contrasted the results of the locally constant and quartic fits obtained by these two methods. However, little difference is observed. The credibility intervals for the locally quadratic and cubic fits are portrayed in Figure 7. We must be careful in interpreting those credibility intervals. bands in Figure 7 are not simultaneous interval estimates. merely indicate the pointwise credibility interval for 5 ( f i ) n + The Those in Figure 7 at each t i. n+ is clear that the bands are blowing up at either end of the data. consistent with the unsatisfactory predictions revealed in Figures 1 and 2. It This is Using the backfitting method, similar results are obtained though they are not presented here. 34 FIG 3 : A c o m p a r i s o n o f t h e l o c a l l y c o n s t a n t f i t s ( t o p ) and q u a d r a t i c f i t s (bottom) u s i n g t h e c r o s s - v a l i d a t i o n m e t h o d ( s o l i d l i n e ) a n d t h e b a c k f i t t i n g method ( d o t t e d l i n e ) f o r t h e pH77 d a t a . 300.0 350.0 100.0 150.0 500.0 TIME (DAT) I 550.0 34a 600.0 650.0 700.0 750.0 800 TABLE I I I : E s t i m a t e d O r d e r (p) v a l u e s o f Cp f o r t h e M.A. pH data. A Cp * Cp ** 0 3.20629 .94601 1 1.65739 .39114 2 . 17786 .02716 3 .06700 .00179 4 .08250 .00018 5 0 0 6 0 0 7 0 0 8 0 0 9 0 0 10 0 A * Cp : Cp o b t a i n e d by t h e c r o s s - v a l i d a t i o n * Cp : Cp o b t a i n e d by t h e b a c k f i t t i n g 34b method method FIG. B. 4: A comparison of the l o c a l l y p o l y n o m i a l f i t s of d i f f e r e n t o r d e r p t o t h e M.A. pH d a t a , ( u s i n g t h e c r o s s - v a l i d a t i o n method) p=0 p=1 p=2 (solid line) § (broken l i n e ) (dotted line) p=2 p=3 p=4 (solid line) s (broken l i n e ) (dotted line) C. J £ p=4 ( s o l i d l i n e ) p=7 ( b r o k e n l i n e ) p=lO(dotted l i n e ) 16.3 4J.35 «-2 34c 94.03 TIME IDflYl 119.9 197.45 249.15 (X10 27 F FIG. 5: A c o m p a r i s o n of t h e l o c a l l y p o l y n o m i a l d i f f e r e n t o r d e r p t o t h e M.A. pH d a t a , b a c k f i t t i n g method) f i t s of ( u s i n g the A. p=0 ( s o l i d l i n e ) p=1 ( b r o k e n l i n e ) •p=2 ( d o t t e d l i n e ) 16.3 p=2 p=3 p=4 p=4 p=7 p=9 —1 66.2 1 9*03 TIME IDPT) 1 119.9 143.73 171.S 197.43 223.3 (XlO 1 I (solid line) (broken l i n e ) (dotted l i n e ) 16.3 C. 42.33 42.33 68.2 (solid line) (broken l i n e ) (dotted l i n e ) 34d 94 0 3 TIME IDRT) IIS.9 145.73 171.6 197.45 223.3 249.13 273.0 IX10 1 ) FIG 6: A c o m p a r i s o n o f t h e l o c a l l y c o n s t a n t f i t s ( t o p ) and q u a r t i c f i t s (bottom) u s i n g t h e c r o s s - v a l i d a t i o n method ( s o l i d l i n e ) a n d t h e b a c k f i t t i n g method ( d o t t e d l i n e ) f o r t h e M.A. pH d a t a . 16.5 16.5 42.35 42.35 68.2 68.2 94.05 119.9 TIME (DAY) T r 94.05 119.9 TIME (OfiY) 34e T 145.75 145.75 171.6 171.6 197.45 197.45 T " 223.3 i 223.3 -1 249.15 r 249.15 1 275.0 1 (XlO ) 275.0 IXIO I 1 FIG 7: 16.5 A comparison of the c r e d i b i l i t y i n t e r v a l s for the l o c a l l y q u a d r a t i c ( s o l i d l i n e ) and cubic (broken l i n e ) f i t s f o r t h e pH77 ( t o p ) a n d t h e M.A. pH (bottom) d a t a , ( u s i n g t h e c r o s s - v a l i d a t i o n method) 42.35 68.2 94.05 119.9 145.73 TIME lORY) 171.6 197.45 223.3 249 15 2/5 IXIO 1 34f Some differences between our approach and classical parametric methods, e.g., least squares, are obvious. O u r approach is more flexible since we can adjust the smoothing parameter c subjectively to get c = 0 gives the least squares polynomial a least squares cannot be recaptured. a "nice looking" fit. T h e choice For example, when p = line passing through the data; Moreover, it is easy curve. however, the underlying interesting feature is shown in F i g . 4 of W Z where and fits when the deleted. are contrasted the least middle of the O n e more squares pH77 linear data are T h e width of the confidence band for the least squares linear fit reaches a m i n i m u m where data are missing. However, the width of the credibility interval for the locally linear fit increases in the gap in accordance with our 4.3 pattern to generate the locally polynomial fit of any order p with corresponding pointwise credibility intervals. locally linear 1, we get expectations. R e s u l t s of choosing t h e "optimal" order p. Applying PSE(p) results the described in section 3.3, the estimated values for the pH77 data are obtained and presented in Tables I V and V . clearly indicate stant fit attains used. procedure that the optimal order is p = the m i n i m u m PSE 0 since the of The locally con- for all of the spans of observations, m , we Because the data are not monotone and the noise is significant, the best fit is simply obtained by taking the weighted average of the observations in the span. A t the same time, we have compared our results with the simplest smoother, the m o v i n g mean, the span. where f, T h e comparisons is estimated are reported 35 by the average of the observations in in Tables I V and V . Although the TABLE I V : E s t i m a t e d v a l u e s o f P S E ( p ) f o r t h e pH77 d a t a , ( p r e d i c t one o b s e r v a t i o n a h e a d ) m 15 20 30 40 50 60 70 p 0 .0738 .0773 .0639 .0615 .0718 .0851 .0393 1 .1298 . 1 305 . 1 009 .0906 .0795 .0986 .0492 2 .3804 .3484 .2300 . 1 657 .1470 . 1 540 .0586 3 5.089 3.719 .5573 .4074 .3671 .2573 .0533 4 39.08 19.21 1 .764 1 .330 .7226 .3866 .0709 5 381 .6 1 02.8 7.905 5.760 1 .795 .471 9 .1126 6 1768 .8 461 .8 18.21 22.68 1 .994 .5999 . 1 643 7 63212 3159.3 46.04 27.64 2.818 .7892 .2380 8 2769231 81 920 378.2 3 1.22 4.956 .9524 . 2700 9 221 31 30632 68.36 5.792 1 .021 .4001 . 1 977 moving mean .0673 .0762 .0848 .1026 .0940 .1064 .0678 35a TABLE V: E s t i m a t e d v a l u e s o f P S E ( p ) f o r t h e pH77 d a t a . ( p r e d i c t t h e a v e r a g e o f 10 o b s e r v a t i o n s a h e a d ) m 15 20 30 40 P 50 60 0 .0335 .0328 .0289 ' .0311 .0387 .01 23 1 . 1 250 . 1 004 .0693 .0746 .0822 .091 9 2 1 .491 1 . 1 05 .7047 .5908 .5588 .7217 3 38.67 18.06 9.541 8.317 6.819 8.919 4 1176.4 430.9 137.4 113.7 63.56 64.48 5 39996 9323.5 1938.8 1408.3 537.2 333 .4 6 1.4E+6 2.5E+5 3751 0 4108.5 1504.1 7 7.9E+7 5.4E+6 6.3E+5 2.2E+5 18543 3038.3 8 1 .2E + 9 7.7E+8 9.5E+6 1.2E+6 75790 5949.6 9 3.5E+9 2.4E+10 2.4E+6 25596 1841.4 .0428 .0507 .0582 .0497 .0342 moving mean 35b 1 9741 362 . 1 .01 94 moving mean has a simple model formulation, it is not capable of giving a better prediction than the locally constant fit. For the M.A. pH data, there is smaller variation in the pH values over the range of time, and again the locally constant fit is best though the results are not shown here. We have tested the efficacy of this procedure on two other data sets. Data set 1: We apply the procedure to data on the annual catch of lobsters by Maine fishermen taken from Morrison (1983) ( abbreviated hereafter by Lobster data ). The data consists of the number of registered lobsters traps used by Maine fishermen over the years 1932-75. The correlation coefficient between these two variables is close to 0.9 and this strongly suggests a linear dependence of the number of traps on time. This is supported by the results given in Table VI. The locally constant, linear and quadratic fits are compared in Figure 8A. These three curves are quite similar to each other even though the locally linear fit minimizes the PSE. are portrayed in Figure 8 B . The credibility intervals of the locally linear fit As we have pointed out above such band3 are only pointwise interval estimates. Data set 2: One hundred equally spaced points were generated from the sinusoid y = 2 cos x — 5 sin z + 6 on the interval [0,10]. We apply the same procedure to the first 80 points and then compare the prediction of the last 20 points by using different locally polynomial fits. As the results show in Table VII, PSE is minimized by different values of j < p depending on the size of span m. For interpolation, all the locally polynomial fits appear to pass through the data points equally 36 TABLE V I : E s t i m a t e d v a l u e s o f P S E ( p ) f o r t h e L o b s t e r d a t a p r e d i c t one o b s e r v a t i o n a h e a d ( t o p ) a n d p r e d i c t t h e a v e r a g e o f 10 o b s e r v a t i o n s a h e a d ( b o t t o m ) . m 20 30 P 0 3.4 E + 4 5.8 E + 4 1 3.3 E + 4 5.6 E + 4 2 8.1 E + 4 1 .3 E + 5 3 2.0 E + 5 3.1 E + 5 4 4.9 E + 5 7.1 E + 5 5 1. 1E + 6 1 .5 E + 6 6 2.6 E + 6 2.7 E + 6 7 6.4 E + 6 3.6 E + 6 8 1 .7 E + 7 3.3 E + 6 9 2.8 E + 6 1 .6 E + 6 1 .7 E + 5 3.5 E + 5 moving mean m 20 30 P 0 7.7 E + 4 1 .4 E + 5 1 1 .3 E + 4 2.2 E + 4 2 6.7 E + 4 5.1 E + 4 3 4.3 E + 5 2.3 E + 5 4 4.0 E + 6 1 .6 E + 6 5 4.7 E + 7 1 .2 E + 7 6 6.1 E + 8 9.5 E + 7 7 9.4 E + 9 1 .4 E + 9 8 1.7 E +11 2.5 E + 10 9 4.9 E +12 4.8 E + 10 1 .5 E + 5 3.3 E + 5 moving mean 36a FIG 1932.0 8A: 1936.3 A comparison of the l o c a l l y constant (broken l i n e ) , l i n e a r ( s o l i d l i n e ) and q u a d r a t i c ( d o t t e d l i n e ) f i t s f o r t h e L o b s t e r d a t a , ( u s i n g t h e b a c k f i t t i n g method) 1940.6 1944.9 1949.2 TIME (YERR) 1953.5 1957.8 1962.1 1966.4 1970.7 1975 FIG 8B: L o c a l l y l i n e a r f i t with the c r e d i b i l i t y i n t e r v a l s f o r the Lobster data, (using the b a c k f i t t i n g method) TABLE V I I : E s t i m a t e d v a l u e s o f P S E ( p ) f o r t h e s i n u s o i d (y = 3 c o s x - 5 s i n x + 6 ) — p r e d i c t one o b s e r v a t i o n a h e a d ( t o p ) and p r e d i c t t h e a v e r a g e o f 10 o b s e r v a t i o n s a h e a d ( b o t t o m ) . m 20 40 60 P 0 7.8 1 1 .5 E - 2 1 .8 E - 2 1 .6 E - 2 2 2.4 E - 4 5.2 E - 4 6.4 3 3.3 E - 6 6.6 1 .0 E - 5 4 4.3 E - 8 1 .7 E - 7 2.4 5 2.0 E - 9 3.2 E - 9 1 .0 E - 8 6 4.5 E - 9 1 .5 E - 9 1 .0 E - 9 7 1 .4 E - 8 1 .3 E - 9 6.0 8 4.2 1 .3 E - 9 4.0 E -10 9 1 .5 E - 8 4.0 1 .2 E - 9 m E - 1 E - 8 9.8 20 E - 1 E - 6 E -10 1. 1 E 40 P 1 .2 E + 1 1 .2 E + 1 1 4.6 E 6.3 2 6.8 E E - 4 E - 7 E -10 60 0 0 0 1 .2 E + 1 0 9.3 E E - 1 1.7 E - 0 6.3 E - 1 3 1.0 E - 1 3.0 E - 1 9.0 E - 1 4 6.5 3.7 E - 2 1 .5 E - 2 5 5.6 E - 4 3.4 E - 3 1 .6 E - 2 6 4.0 2.4 E - 4 1 .8 E - 4 7 2.9 E - 4 5.0 1 . 1E - 4 8 1.9 E - 2 3.5 E - 4 9 1 .6 E + 6 9.0 E - 3 E - 5 36d E - 5 E - 5 3.0 0 E - 5 1.0 E - 5 well since no noise exists in the data. 9(A-B). These situations are illustrated in Figures Figure 9 B indicates that locally the 9th order polynomial gives the best prediction of the last 20 data points. From these two examples, we can see that the proposed procedure is capable of selecting the locally polynomial fit which captures the underlying pattern of the data. To conclude this section, it should be added that we have examined the behaviour of the yearly differences of the M.A. pH data over the range of time, D(t i) n+ = R(t i) n+ — R(t i n+ existed over the years. — 1 year ). This is done to see if any seasonality Figure 10 is a plot of the locally constant fit with credibility intervals. No apparent seasonality is observed in the picture. However, we can see that D(t ) n+l is increasing and decreasing in alternating years. a yearly cycle of the change of pH is clearly shown. Of course, related subject knowledge would be needed to explain such a phenomenon. 37 Such FIG 9A: A c o m p a r i s o n of t h e l o c a l l y l i n e a r ( b r o k e n l i n e ) , 5 t h o r d e r p o l y n o m i a l ( s o l i d l i n e ) and 1 0 t h o r d e r polynomial (dotted l i n e ) f i t s f o r the s i n u s o i d (y = 3 c o s x - 5 s i n x + 6) d a t a , ( u s i n g t h e c r o s s - v a l i d a t i o n method) • 0.0 0.99 1.98 2.97 X 3.96 4.95 5.94 6.93 7.92 8.91 9. FIG i 0.0 1 0.99 9B: A c o m p a r i s o n of t h e l o c a l l y 5 t h o r d e r p o l y n o m i a l ( s o l i d l i n e ) , 7th order p o l y n o m i a l (broken l i n e ) and 9 t h o r d e r p o l y n o m i a l ( d o t t e d l i n e ) f i t s f o r t h e s i n u s o i d (y = 3 c o s x - 5 s i n x + 6) d a t a , ( u s i n g t h e b a c k f i t t i n g method) 1 1.96 1 2.97 X 1 3.96 1 4.95 1 5.94 1 6.93 1 7.92 1 8.91 1 9.9 FIG 10: L o c a l l y c o n s t a n t f i t w i t h the c r e d i b i l i t y f o r the d i f f e r e n c e of the M.A. pH d a t a , ( u s i n g t h e b a c k f i t t i n g method) intervals 5. CONCLUSIONS. In section 2, we survey some of the popular nonparametric smoothing methods. A brief comparison of these methods with the WZ-approach is made. As mentioned in section 3.2, behaviour of XV(c), found is finite. through the examination of the asymptotic we can be sure that the cross-validatory estimate of c being Unfortunately, the cross-validatory estimates of c are relatively large and the resulting curves are wiggly. The problem is that we allow the data "speak freely for itself but neglect the formulated model. As an alternative, the backfitting method is more coherent with the theory in section 3.1 and the resulting curves are comparatively smoother in the sense that le?s turns are generated. Also, it is computationally fast and we can estimate cy for all j < p within one run once c p+1 is fixed. In section 4, the efficacy of a data-based procedure was demonstrated for choosing the "optimal" order p. examples too. Satisfactory results are achieved in two other It would be interesting to repeat these tests on data generated by a smooth function plus noise to see if the smooth function would be recovered. In practice, we can arbitrarily set c i p+ = 0 for seme high order p and then apply the procedure to choose the appropriate order of locally polynomial fit. Once the order is determined, we can use the backfitting estimate of c to generate the interpolants or predictands with corresponding pointwise credibility intervals. Nevertheless, as we have pointed out above the fit is fairly robust with respect to the order of the polynomial, a slight difference of the order p would not change the appearance of the resulting curve very much. 38 On the whole, the proposed approach isflexibleand simple to implement. In general, it is easy to generate the locally polynomial fit of any order if it is needed and appropriate degree of smoothness can be obtained through the backfitting method. Moreover, this approach does provide explicit estimates which are much different in form to those derived by classical methods. In conclusion, we suggest that further research be carried out, including 1. a proper Bayesian analysis for controlling the extrapolation of S(t), and 2. the development of the theory for multivariate space-time series and the regression problem. 39 BIBLIOGRAPHY Becker, R. A. and Chambers, J. M . (1984): S: An Interactive for Data Analysis and Graphics. Environment California: Wadsworth. Bickel, P. J. and Rosenblatt, M . (1973): density function estimates. Ann. Statist. On some global measures of the 1 1071-1095. Bierens, H . J. (1983): Sample moments integrating normal kernel estimators of multivariate density and regression functions. Sankhya B 45 160-192. Brown, R. G . and Meyer, R. F. (1961): The fundamental theorem of exponential smoothing. Operations Research 9 673-685. Clark, R. M . (1977): Nonparametric estimation of a smooth regression function. J. R. Statist. Soc. B 39 107-113. Cleveland, W. S. (1979): scatterplots. J. Amer. Statist. Collomb, G . (1981): Robust locally weighted regression and smoothing Assoc. 74 829-836. Estimation nonparametrique de la regression: revue bibliographique. Inter. Statist. Review 49 75-93. Craven, P. and Wahba, G. (1979): Smoothing noisy data with spline functions. Numer. Math. 31 377-403. Dubrule, O. (1983): Two methods with different objectives: kriging. Math. Geology 15 245-257. 40 splines and Dubrule, O. (1984): Comparing splines and kriging. Computers and Geo- sciences 10 327-338. Enns, P. G., Machak, J. A., Spivey, W. A., and Wrobleski, W. J. (1982): Forecasting applications of an adaptive multiple exponential smoothing model. Management Science 28 1035-1044. Friedman, J. H. (1984): A variable span smoother. Stanford University LCS Technical Report No. 5. Gardner, E. S. (1985): Exponential smoothing: the state of art. J. Forecasting 4 1-28. Gentleman, R., Zidek, J. and Olsen, R. (1985): precipitation monitoring data base. MAP3S/PCN acid rain University of British Columbia Statistics Department Technical Report No. 19. Hardle, VV. and Marron, J. S. (1985): Optimal bandwidth election in nonparametric regression function estimation. Ann. Statist. 13 1465-1481. Harvey, A. C. (1984): nential smoothing model. Analysis and generalization of a multivariate expoLSE Econometrics Programme Discussion Paper No. A.41. Huber, P. J. (1979): Robust smoothing. In Rob ustness in Statistics. (R. L. Launer and G. N. wilkinson, eds.) New York: Academic Press. IMSL (1980): Internationa] Matiematicai and Statistics Libraries, Inc., Sub41 routine 1 C S S C V . IMSL (1982): International routine Mathematical and Statistics Libraries, Inc., Sub- ZXLSF. Jones, R. H. (1Q66): Exponential smoothing for multivariate times series. J. R. Statist. Soc. B 34 385-392. Lenth, R. V. (1977): Robust splines. Comm. in Statist.-A 6 847-854. Lindley, D. V. and Smith, A. F. M . (1972): Bayes estimates for the linear model. J. R. Statist. Soc. B 34 1-18. Makridakis, S., Andersen, A., Carbone, R., Fildes, R., Eibcn, M . , Newton, J., Lewandowski, R., Parzen, R. and Winkler, R. (1982): The accuracy of extrapolation (time series) methods: results of a forcasting competition. J. of Forecasting 1 111-153. Makridakis, S. and Wheelwright, S. C. (1978): Interactive variate and Multivariate Methods Morrison, D.F. (1983): Applied Forecasting: Uni- (2nd ed.). San Francisco: Holden-Day. Linear Statistical Methods. Prentice-Hall. Nadaraya, E. A. (1964): On estimating regression. Theory Probab. Appl. 9 141-142. Priestley, M . B. and Chao, M . T . (1972): Nonparametric function fitting. J. R. Statist. Soc. B 34 385-392. 42 Reinsch, C. H . (1967): Smoothing by spline functions. Numer. Math. 10 177-183. Rosenblatt, M . (1956): Remarks on some nonparametric estimates of a density function. Ann. Math. Statist. 27 832-837. Rosenblatt, M . (1971): Curve estimates. Ann. Math. Schoenberg, I. J. (1964): Proc. Nat. Acad. Statist. 42 1815-1842. Spline functions and the problem of graduation. Sci. USA 52 947-950. Schuster, E. F. (1972): Joint asymptotic distribution of the estimated regression function at a finite number of distinct points. Ann. Math. Statist. 43 84-88. Silverman, B. V . (1982): transform. Appl. Statist. Kernel density estimation using the fast Fourier 31 93-99. Silverman, B. V . (1984): A fast and efficient cross-validation method for smoothing parameter choice in spline regression. J. Amer. Statist. Assoc. 79 584-589. Silverman, B. V . (1985): Some aspects of the spline smoothing approach to nonparametric regression curve fitting. J. R. Statist. Soc. B 47 1-21. Stone, C. J. (1977): Consistent nonparametric regression. 595-645. 43 Ann. Statist. 5 Stone, M . (1974): Cross-validatory choice and assessment of statistical predictions. J. R. Statist. Soc. B 36 111-133. Trigg, D. W. and Leach, A. G . (1967): Exponential smoothing with an adaptive response rate. Opl. Res. Q. 1 8 53-63. UBC MATRIX (1982): A guide to solving matrix problems. University of British Columbia Computing. Documentation. Utreras, F. I. (1979): Cross-validation techniques for smoothing spline functions in one or two dimensions. In Smoothing Techniques For Curve Estimation. (T. Gasser and M . Rosenblatt, eds.) Lecture Notes in Mathematics, No. 757. Springer-Verlag, Berlin. Utreras, F. I. (1981): On computing robust splines and applications. J. Sci. Stat. Comput. SIAM 2 153-163. Wahba, G. (1980): Spline bases, regularization, and generalised crosscross-validation for solving approximation problems with large quantities of noisy data. In Approximation Theory III (E. W. Cheney ed.) New York: Academic Press. Wahba, G. (1983): Bayesian "confidence intervals" for the cross-validation smoothing spline. J. R. Statist. Wahba, G. (1984a): Soc. B 45 133-150. Cross-validated spline methods for the estimation of multivariate functions from data on functionals. In Statistics: An Appraisal. A. David and H. T . David eds.) The Iowa State University Press. 44 (H. Wahba, G. (1984b): Surface fitting with scattered noisy data on Euclidean d-space and on the sphere. Rocky Mountain Journal of Mathematics 14 281-299. Wahba, G. and Wendelberger, J. (1980): Some new mathematical methods for variational objective analysis using splines and cross-validation. Monthly Review Weather 108 1122-1143. Wahba, G. and Wold, S. (1975): A completely automatic french curve: fitting spline functions by cross-validation. Comm. in Statist. 4 1-17. Watson, G . S. (1964): Smooth regression analysis. Sankhya A 26 359-372. Watson, G. S. (1984): Smoothing and interpolation by kriging and with splines. Math. Geology 16 601-615. Weerahandi, S. and Zidek, J. (1985): by Bayesian nonparametric methods. Smoothing locally smooth processes University of British Columbia Statistics Department Technical Report No. 26. Weerahandi, S. and Zidek, J. (1986): Analysis of multiple time series by Bayesian and Empirical Bayesian nonparametric methods. University of British Columbia Statistics Department Technical Report No. 33. Wegman, E. J. (1980): Two approaches to nonparametric regression: splines and isotonic inference. In Recent Developments Analysis (K. in Statistical Inference and Data Matusita ed.). North-Holland. Wegman, E . J. and Wright, I. W. (1983): Splines in statistics. 45 J. Amer. Statist. Assoc. 78 351-365. Whittaker, E. T . (1923): On a new method of graduation. Proc. Math. Edinburgh Soc. 41 63-75. Winters, P. R. (1960): averages. Management Sci. Forscasting sales by exponentially weighted moving 324-342. 46 APPENDIX This appendix contains Fortran subprograms for the computations of (i) the smoothing parameter c by cross-validation, (ii) the smoothing parameter c by "backfitting", (iii) the smooth curve estimates with corresponding credibility intervals, and (iv) the values of P S E . 47 C C C C C C C C C C c c c c c c c c c c c c c c c c c c c c C C C C C DOUBLE PRECISION C C C C C C F(C) T H I S I S A R O U T I N E FOR COMPUTING T H E SMOOTHING PARAMETER C BY U S I N G T H E C R O S S - V A L I D A T I O N METHOD THROUGH T H E IMSL M I N I M I Z I N G R O U T I N E Z X L S F WHICH CAN BE R E P L A C E D BY ANOTHER M I N I M I Z I N G ROUTINE. T H E C A L L I N G PROGRAM MUST I N C L U D E T H E S T A T E M E N T "CALL ZXLSF(F,X,STEP,BOUND,XACC,MAXFN,IER)". D O U B L E P R E C I S I O N I S U S E D TO IMPROVE T H E A C C U R A C Y OF E S T I M A T I O N . INPUT V A R I A B L E S A R E : ARRAY OF O B S E R V E D V A L U E S ( OREDERED I N T I M E ) ARRAY OF T I M E - V A L U E S ( ORDERED I N T I M E ) D I M . OF ARRAYS Y AND L T H E ORDER OF T H E L O C A L L Y P O L Y N O M I A L F I T P L U S 1 ; IN T H I S C A S E , N P = 1 , T H A T I S L O C A L L Y C O N S T A N T F I T . NP2 = 2 x NP NP 1 = NP2 - 1 SC THE SCALING FACTOR PROD = WORKING ARRAY OF D I M . NP2 ATY = WORKING ARRAY OF D I M . NP ARRAYS ATA AND A T A I OF D I M . NP x N P , AND WORKING ARRAY I PERM OF D I M . NP2 A R E U S E D I N T H E S U B R O U T I N E I N V ( A L L T H E ABOVE V A R I A B L E S MUST B E D E F I N E D I N T H E C A L L I N G PROGRAM ) Y L M NP c SUM OUTPUT V A R I A B L E S A R E : E S T I M A T E OF T H E SMOOTHING P A R A M E T E R = C R O S S - V A L I D A T E D SUM OF SQUARED ERRORS S U B R O U T I N E S C A L L E D BY F A R E : I N V ( C A N B E FOUND I N T H E D O C U M E N T A T I O N : U B C M A T R I X ) . I N V CAN* B E R E P L A C E D BY ANOTHER M A T R I X I N V E R S E ROUTINE WHERE A T A I S I N P U T AND A T A I I S O U T P U T . R E A L * 8 A T A ( 1 ,1 *L(80),IPERM(2 ) C C C C FUNCTION ),ATAI(1 ,1 ),ATY(1 ),PROD(2 ),Y(80), T H E D I M E N S I O N S OF T H E A B O V E ARRAYS MUST B E S P E C I F I E D ; IN THIS C A S E , N P = 1 , THAT IS L O C A L L Y CONSTANT F I T . REAL * 8 COMMON SC,Z,C,SUM,TEMP,V,W,P,YY,BB Y,L,M,NP,SC,NP1,NP2 T H E A B O V E "COMMON" S T A T E M E N T I S U S E D I N ORDER T O P A S S T H E PARAMETERS FROM T H E C A L L I N G R O U T I N E TO T H I S R O U T I N E . T H E C A L L I N G PROGRAM MUST I N C L U D E T H E SAME "COMMON" S T A T E M E N T . Z=DABS(C) 48 5 10 12 14 20 C C C C C C 22 COMPUTE T H E R E G R E S S I O N COEFFICIENTS CALL I N V ( N P , N P , A T A , I PERM,NP,ATAI,DET,JEXP,COND) DO 30 1 = 1 , N P BB=BB+ATAI(1,I)*ATY(I) BB I S T H E C R O S S - V A L I D A T O R Y 30 25 C C C C C C SUM=O.DO DO 25 J = 1 , M BB=0.D0 DO 5 J J = 1 , N P ATY(JJ)=0.D0 DO 10 1 1 = 1 , N P 1 PR0D(II)=0.D0 NT=L(J) DO 20 1 = 1 , M •IF(I . E Q . J ) GO TO 20 P=SC*(L(I)-NT) W=1.D0/(1.D0+Z*(P**NP2)) V=W DO 12 1 1 = 1 , N P 1 PR0D(II)=PR0D(II)+V V=V*P CONTINUE V=W YY=Y(I) DO 14 J J = 1 , N P ATY(JJ)=ATY(JJ)+V*YY V=V*P CONTINUE CONTINUE DO 22 1 1 = 1 , N P DO 22 J J = 1 , N P ATA(11,JJ)=PROD(11+JJ-1 ) CONTINUE 99 ESTIMATE OF T H E SMOOTH CONTINUE SUM=SUM+(Y(J)-BB)**2 CONTINUE WRITE(6,99) Z,SUM F0RMAT(E18.7,2X,E15.8) Z I S T H E A B S O L U T E V A L U E OF C AND SUM I S T H E C R O S S V A L I D A T E D SUM OF SQUARED E R R O R S . THEY A R E OUTPUT SO T H A T T H E V A L U E OF C WHICH M I N I M I Z E S SUM C A N BE DETERMINED. F = SUM RETURN END 49 * C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C SUBROUTINE BFIT(Y,L,B,ATA,ATY,PROD,ATAI,I CC,TA,M,NP,SC,C,NP2,N) PERM,TT,TTT T H I S I S A R O U T I N E FOR COMPUTING T H E SMOOTHING PARAMETER C BY U S I N G T H E B A C K F I T T I N G M E T H O D . INPUT VARIABLES A R E : ARRAY OF O B S E R V E D V A L U E S (ORDERED I N T I M E ) ARRAY OF T I M E - V A L U E S (ORDERED I N T I M E ) D I M . OF ARRAYS Y AND L = T H E ORDER OF T H E L O C A L L Y P L O Y N O M I A L F I T PLUS 1 ; I N T H I S C A S E , NP= 1 1 SC THE SCALING FACTOR C T H E SMOOTHING PARAMETER OF T H E 1OTH ORDER L O C A L L Y POLYNOMIAL F I T ; IN T H I S C A S E , C=0.0D0 TA ARRAY OF T H E D E R I V A T I V E S OF T H E SMOOTH OF D I M . M AND N NP2 = 2 x NP N = NP - 1 PROD = WORKING ARRAY OF D I M . NP2 B , T T , T T T A R E WORKING ARRAYS OF D I M . NP A T A AND A T A I OF D I M . NP x N P , A T Y OF D I M . NP AND I P E R M OF D I M . NP2 A R E ARRAYS U S E D I N T H E S U B R O U T I N E S I N V AND DGMATV ( A L L T H E A B O V E V A R I A B L E S MUST B E D E F I N E D I N T H E C A L L I N G PROGRAM. ) Y L M NP CC(K)= OUTPUT V A R I A B L E S A R E : SMOOTHING PARAMETER OF T H E K T H ORDER POLYNOMIAL F I T , K = 1,...,N LOCALLY S U B R O U T I N E S C A L L E D BY B F I T A R E : I N V AND DGMATV ( C A N B E FOUND I N T H E D O C U M E N T A T I O N : U B C M A T R I X ) . I N V CAN B E R E P L A C E D BY ANOTHER M A T R I X I N V E R S E R O U T I N E WHERE A T A I S I N P U T AND A T A I I S O U T P U T . DGMATV CAN B E R E P L A C E D BY ANOTHER M A T R I X M U L T I P I C A T I O N R O U T I N E WHERE A T A I AND A T Y A R E INPUTS AND B I S O U T P U T . IMPLICIT R E A L * 8 ( A - H , 0 - Z ) DIMENSION Y ( M ) , L ( M ) , B ( N P ) , A T A ( N P , N P ) , A T Y ( N P ) *,PROD(NP2),ATAI(NP,NP),IPERM(NP2),TT(N),TTT(N),CC(N), *,TA(M,N) NP1=NP2-1 DD=0.0D0 C C C C COMPUTE T H E S A M P L E V A R I A N C E OF T H E N O I S E T H E D E R I V A T I V E S OF T H E SMOOTH T A 2 DO 2 1 = 1 , N TT(I)=0.0D0 TTT(I)=0.0D0 50 DD AND 5 10 12 14 19 22 21 25 C C C DO 25 J = 1 , M DO 5 J J = 1 , N P ATY(JJ)=0.D0 DO 10 1 1 = 1 , N P 1 PROD(II)=0.D0 NT=L(J) DO 19 1 = 1 , M P=SC*(L(I)-NT) W=1.D0/(1.D0+C*(P**NP2)) V=W DO 12 1 1 = 1 , N P 1 PROD(II)=PR0D(II)+V V=V*P CONTINUE V=W YY=Y(I) DO 14 J J = 1 , N P ATY(JJ)=ATY(JJ)+V*YY V=V*P CONTINUE CONTINUE DO 22 1 1 = 1 , N P DO 22 J J = 1 , N P ATA(11,JJ)=PROD(11+JJ- 1) CONTINUE CALL I N V ( N P , N P , A T A , I PERM,NP,ATAI CALL DGMATV(ATAI,ATY,B,NP,NP,NP) TEMP=ATAI(1,1) DD=DD+(Y(J)-B(1))**2/(1.DO+TEMP) DO 21 K = 1 , N TA(J,K)=B(K+1) TT(K)=TT(K)+TA(J,K) CONTINUE DD=DD/(M-NP) COMPUTE T H E SMOOTHING PARAMETERS 44 66 55 77 88 DO 44 1 = 1 , N TT(I)=TT(I)/M DO 55 1 = 1 , N DO 66 K = 1 , M TTT(I)=TTT(I)+(TA(K,I)~TT(I))**2 CONTINUE DO 77 1 = 1 , N TTT(I)=TTT(I)/(M~1) DO 88 K = 1 , N CC(K)=TTT(K)/(DD) CONTINUE RETURN END 51 ,DET,JEXP,COND) CC(K), K = 1,...,N C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C SUBROUTINE S M O O T H ( Y , L , M M , M , N P , S C , C , K K , U L , V L , * P R O D , A T Y , T E M , A , A T A , A T A I , I P E R M , N P 1 ,NP2) T H I S I S T H E R O U T I N E FOR COMPUTING T H E I N T E R P O L A N T S AND T H E P R E D I C T A N D S OF T H E SMOOTH CURVE WITH T H E CORRESPONDING C R E D I B I L I T Y I N T E R V A L S . INPUT VARIABLES A R E : ARRAY OF O B S E R V E D V A L U E S (ORDERED I N T I M E ) ARRAY OF T I M E - V A L U E S (ORDERED I N T I M E ) = T H E T O T A L NUMBER OF T H E I N T E R P O L A N T S AND T H E P R E D I C T A N D S OF T H E SMOOTH CURVE M T O T A L NUMBER OF O B S E R V A T I O N S OF Y NP = T H E ORDER OF T H E L O C A L L Y POLYNOMIAL F I T PLUS 1 SC = THE SCALING FACTOR C T H E C R O S S - V A L I D A T O R Y OR T H E B A C K F I T T I N G SMOOTHING PARAMETER NP2 = 2 x NP NP1 = NP2 - 1 KK = T H E NUMBER OF P R E D I C T A N D S BEYOND T H E E I T H E R S I D E OF T H E RANGE OF DATA PROD= WORKING ARRAY OF D I M . NP2 A T Y = WORKING ARRAY OF D I M . N P 1 T E M = WORKING ARRAY OF D I M . MM ARRAYS A T A AND A T A I OF D I M . NP AND NP AND WORKING ARRAY I PERM OF D I M NP2 A R E USED I N T H E S U B R O U T I N E I N V ( A L L T H E A B O V E V A R I A B L E S MUST B E D E F I N E D I N T H E C A L L I N G PROGRAM. ) Y L MM A UL VL SE OUTPUT VARIABLES A R E : = ARRAY OF D I M . MM OF T H E I N T E R P O L A N T S OR T H E P R E D I C T A N D S OF T H E SMOOTH CURVE = ARRAY OF D I M . MM OF T H E U P P E R L I M I T S OF T H E C R E D I B I L I T Y INTERVALS = ARRAY OF D I M . MM OF T H E LOWER L I M I T S OF T H E C R E D I B I L I T Y INTERVALS = STANDARD ERROR OF T H E I N T E R P O L A N T S OR T H E PREDICTANDS S U B R O U T I N E S C A L L E D BY SMOOTH A R E : I N V ( C A N B E FOUND I N T H E D O C U M E N T A T I O N : U B C M A T R I X ) . I N V CAN B E R E P L A C E D BY ANOTHER M A T R I X I N V E R S E R O U T I N E WHERE A T A I S INPUT AND A T A I I S O U T P U T . IMPLICIT R E A L * 8 ( A - H , 0 - Z ) DIMENSION Y ( M M ) , L ( M M ) , A T A ( N P , N P ) , A T Y ( N P ) , P R 0 D ( N P 2 ) *,ATAI(NP,NP),IPERM(NP2),A(MM),TEM(MM),UL(MM),VL(MM) N=NP-1 K1=KK+1 K2=KK+M DD=0.D0 DO 2 I = 1 , K K 52 L(I)=5*(1-1)+300 L(I+92)=5*(I-1)+730 COMPUTE T H E SAMPLE V A R I A N C E OF T H E N O I S E DD DO 25 J = K 1 , K 2 BB=0.D0 DO 5 J J = 1 , N P ATY(JJ)=0.D0 DO 10 1 1 = 1 , N P 1 PR0D(II)=0.D0 NT=L(J) DO 19 I = K 1 , K 2 P=SC*(L(I)-NT) W=1.D0/(1.D0+C*(P**NP2)) V=W DO 12 1 1 = 1 , N P 1 PR0D(II)=PR0D(II)+V V=V*P CONTINUE V=W YY=Y(I) DO 14 J J = 1 , N P ATY(JJ)=ATY(JJ)+V*YY V=V*P CONTINUE CONTINUE DO 22 1 1 = 1 , N P DO 22 J J = 1 , N P ATA(II,JJ)=PROD(II+JJ-1) CONTINUE CALL I N V ( N P , N P , A T A , I PERM,NP,ATAI,DET,JEXP,COND) DO 30 1 = 1 , N P BB=BB+ATAI(1,I)*ATY(I) CONTINUE TEMP=ATAI(1,1) DD=DD+(Y(J)-BB)**2/(1.DO+TEMP) CONTINUE DD=DD/(M-NP) COMPUTE^THE STANDARD ERROR S E DO 29 J = l , M M BB=0.D0 DO 6 J J = 1 , N P ATY(JJ)=0.D0 DO 16 1 1 = 1 , N P 1 PROD(II)=O.DO NT=L(J) DO 49 I = K 1 , K 2 P=SC*(L(I)-NT) W=1.D0/(1.D0+C*(P**NP2)) 53 v=w 27 74 49 82 39 29 C C C DO 27 I I = 1 , N P 1 PR0D(II)=PROD(II)+V V=V*P CONTINUE V=W YY=Y(I) DO 74 J J = 1 , N P ATY(JJ)=ATY(JJ)+V*YY ,V=V*P CONTINUE CONTINUE DO 82 1 1 = 1 , N P DO 82 J J = 1 , N P ATA(II,JJ)=PR0D(II+JJ-1) CONTINUE CALL I N V ( N P , N P , A T A , I PERM,NP,ATAI,DET,JEXP,COND) DO 39 1 = 1 , N P BB=BB+ATAI(1,1)*ATY(l) CONTINUE TEM(J)=ATAI(1,1) A(J)=BB CONTINUE DO 101 1=1,MM SE=DSQRT(DD*TEM(I)) COMPUTE T H E C R E D I B I L I T Y 101 INTERVALS UL(I)=A(I)+1.96*SE VL(I)=A(I)-1.96*SE CONTINUE RETURN END 54 SUBROUTINE c c c r c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c r* L- c c c C C THIS PSE(SC,N1,N2,C,M,N4,NP,NP1,NP2,N,Y,L,CC, ATA,ATY,ATAI,PROD,B,I PERM,TT,TTT,TA,VAR, R, AA) I S A R O U T I N E FOR COMPUTING T H E V A L U E S P R E D I C T ONE O B S E R V A T I O N A H E A D . OF PSE INPUT VARIABLES A R E : ARRAY OF O B S E R V E D V A L U E S (ORDERED I N T I M E ) ARRAY OF T I M E - V A L U E S (ORDERED I N T I M E ) D I M . OF T H E ARRAYS Y AND L T H E ORDER OF T H E L O C A L L Y P O L Y N O M I A L F I T P L U S 1 T H E SMOOTHING PARAMETER OF T H E 1OTH ORDER L O C A L L Y POLYNOMIAL F I T ; IN T H I S C A S E , C=0.0D0 SC THE SCALING FACTOR N1,N2 = T H E S T A R T I N G L I M I T S OF T H E SPAN OF O B S E R V A T I O N S N NP - 1 NP2 2 X NP NP 1 NP2 - 1 M - N2 N4 B ARRAY OF T H E D E R I V A T I V E S OF T H E SMOOTH CC ARRAY OF T H E SMOOTHING PARAMETERS OF DIFFERENT ORDER L O C A L L Y POLYNOMIAL F I T S AA ARRAY OF T H E E S T I M A T E S OF F U T U R E OBSERVATIONS VAR WORKING ARRAY OF D I M . M PROD = WORKING ARRAY OF D I M . NP2 TA WORKING ARRAY OF D I M . N2 AND N T T , T T T = WORKING ARRAY OF D I M . N ARRAYS A T A AND A T A I OF D I M . NP x N P , ARRAY A T Y OF D I M . NP AND ARRAY IPERM OF D I M . NP2 A R E U S E D I N T H E R O U T I N E S I N V AND DGMATV ( CAN BE FOUND I N T H E D O C U M E N T A T I O N : UBC M A T R I X ) . I N V CAN BE R E P L A C E D BY ANOTHER M A T R I X I N V E R S E ROUTINE WHERE A T A I S INPUT AND ATAI I S O U T P U T . DGMATV CAN BE R E P L A C E D BY ANOTHER M A T R I X M U L T I P I C A T I O N ROUTINE WHERE A T A I AND A T Y ARE I N P U T S AND B I S O U T P U T . ( A L L T H E A B O V E V A R I A B L E S MUST BE D E F I N E D I N T H E C A L L I N G PROGRAM ) Y L M NP C RU) OUTPUT VARIABLES A R E : = ' E S T I M A T E D V A L U E OF P S E OF T H E POLYNOMIAL F I T , I = 1,...,N ITH ORDER LOCALLY S U B R O U T I N E S C A L L E D BY P S E A R E : BACK AND P R E D ( A T T A C H E D TO T H I S ROUTINE) IMPLICIT REAL*8(A-H,0-Z) DIMENSION Y ( M ) , L ( M ) , A A ( N , N 4 ) , C C ( N , N 4 ) , A T A ( N P , N P ) *,ATY(NP),ATAI(NP,NP),PROD(NP2),B(NP),IPERM(NP2),TT(N4), *TTT(N4),TA(N2,N),VAR(M),R(N) COMPUTE T H E BACKFITTING VALUES 55 OF C AND T H E SAMPLE C C V A R I A N C E OF T H E N O I S E CALL BACK(SC,N1,N2,C,M,N4,NP,NP1,NP2,N,Y,L,CC,ATA, *ATY,ATAI,PROD,B,IPERM,TT,TTT,TA,VAR) C C C COMPUTE T H E P R E D I C T I V E VALUES OF F U T U R E OBSERVATIONS CALL PRED(SC,NP,NP1,NP2,N,M,N1,N2,N4,Y,L,CC,ATA,ATY, * A T A I , P R O D , I PERM,AA) C C C COMPUTE T H E E S T I M A T E D 101 100 C C 3 2 5 10 VALUES OF P S E DO 100 1 = 1 , N R(I)=0.0 DO 101 J = 1 , N 4 R(I)=R(I)+(Y(N2+J)-AA(l,J))**2 CONTINUE R(I)=R(I)/N4 CONTINUE RETURN END SUBROUTINE BACK(SC,N1,N2,C,M,N4,NP,NP1,NP2,N,Y,L,CC, *ATA,ATY,ATAI,PROD,B,I PERM,TT,TTT,TA,VAR) IMPLICIT R E A L * 8 ( A - H , 0 - Z ) DIMENSION Y ( M ) , L ( M ) , B ( N P ) , A T A ( N P , N P ) , P R O D ( N P 2 ) , *ATAI(NP,NP),IPERM(NP2),TT(N),TTT(N),CC(N,N4),TA(N2,N), *ATY(NP),VAR(M) K1 =N1 K2=N2 DO 18 N N = 1 , N 4 DD=0.D0 DO 2 1 = 1 , N DO 3 J = 1 , N 2 TA(J,I)=0.D0 CONTINUE TT(I)=0.DO TTT(I)=0.D0 CONTINUE DO 25 J = R 1 , K 2 DO 5 J J = 1 , N P ATY(JJ)=0.D0 DO 10 1 1 = 1 , N P 1 PROD(lI)=0.D0 NT=L(J) DO 19 I = K 1 , K 2 P=SC*(L(I)-NT) W=1.D0/(1.D0+C*(P**NP2)) V=W DO 12 1 1 = 1 , N P 1 PROD(lI)=PROD(lI)+V 56 12 14 19 22 21 25 44 66 55 47 88 18 C C V=V*P CONTINUE V=W YY=Y(I) DO 14 J J = 1 , N P ATY(JJ)=ATY(JJ)+V*YY V=V*P CONTINUE CONTINUE DO 22 1 1 = 1 , N P • DO 22 J J = 1 , N P ATA(II,JJ)=PROD(II+JJ-1) CONTINUE CALL I N V ( N P , N P , A T A , I PERM,NP,ATAI,DET,JEXP,COND) CALL D G M A T V ( A T A I , A T Y , B , N P , N P , N P ) TEMP=ATAI(1,1) DD=DD+(Y(J)-B(1))**2/(1.DO+TEMP) DO 21 K = 1 , N TA(J,K)=B(K+1) TT(K)=TT(K)+TA(J,K) CONTINUE DD=DD/(N2-NP) DO 44 1 = 1 , N TT(I)=TT(I)/N2 DO 55 1 = 1 , N DO 66 K = K 1 , K 2 TTT(I)=TTT(I)+(TA(K,I)-TT(I))**2 CONTINUE DO 47 1 = 1 , N TTT(I)=TTT(I)/(N2-1) DO 88 K = 1 , N CC(K,NN)=TTT(K)/DD CONTINUE VAR(NN)=DD K1=K1+1 K2=K2+1 CONTINUE RETURN END SUBROUTINE PRED(SC,NP,NP1,NP2,N,M,N1,N2,N4,Y,L,CC,ATA * , A T Y , A T A I , P R O D , I PERM,AA) IMPLICIT REAL*8(A-H,0-Z) DIMENSION Y ( M ) , L ( M ) , A T A ( N P , N P ) , P R O D ( N P 2 ) , A A ( N , N 4 ) *,ATAI(NP,NP),IPERM(NP2),CC(N,N4),ATY(NP) K1=N1 K2=N2 DO 71 I P = 1 , N 4 DO 72 J P = 1 , N NQ=JP NQ2=2*JP 57 74 75 77 78 76 79 80 72 71 NQ1=NQ2-1 BB=O.DO DO 74 J J = 1 , N Q ATY(JJ)=0.D0 DO 75 11=1,NQ1 PR0D(II)=0.D0 KT=L(N2+IP) DO 76 J = K 1 , K 2 P=SC*(L(J)-KT) W=1.D0/(1.D0+CC(JP,IP)*(P**NQ2) ) V=W DO 77 11=1,NQ1 PR0D(II)=PR0D(II)+V V=V*P CONTINUE V=W YY=Y(J) DO 78 J J = 1 , N Q ATY(JJ)=ATY(JJ)+V*YY V=V*P CONTINUE CONTINUE DO 79 1 1 = 1 , N Q DO 79 J J = 1 , N Q ATA(II,JJ)=PROD(II+JJ-1) CONTINUE CALL I N V ( N Q , N P , A T A , I PERM,NP,ATAI,DET,JEXP,COND) DO 80 I = 1 , N Q BB=BB+ATAI(1,I)*ATY(I) CONTINUE AA(JP,IP)=BB CONTINUE K1=K1+1 K2=K2+1 CONTINUE RETURN END 58
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Smoothing locally regular processes by Bayesian nonparametric...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Smoothing locally regular processes by Bayesian nonparametric methods, with applications to acid rain.. Ma, Hon Wai 1986-12-31
pdf
Page Metadata
Item Metadata
Title | Smoothing locally regular processes by Bayesian nonparametric methods, with applications to acid rain data analysis |
Creator |
Ma, Hon Wai |
Publisher | University of British Columbia |
Date | 1986 |
Date Issued | 2010-06-27T16:53:36Z |
Description | We consider the problem of recovering an unknown smooth function from the data using the Bayesian nonparametric approach proposed by Weerahandi and Zidek (1985). Selected nonparametric smoothing methods are reviewed and compared with this new method. At each value of the independent variable, the smooth function is assumed to be expandable in a Taylor series to the pth order. Two methods, cross-validation and "backfitting" are used to estimate the a priori unspecified hyperparameters. Moreover, a data-based procedure is introduced to select the appropriate order p. Finally, an analysis of an acid-rain, wet-deposition time series is included to indicate the efficacy of the proposed methods. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Collection |
Retrospective Theses and Dissertations, 1919-2007 |
Series | UBC Retrospective Theses Digitization Project |
Date Available | 2010-06-27 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0096763 |
URI | http://hdl.handle.net/2429/26004 |
Degree |
Master of Science - MSc |
Program |
Statistics |
Affiliation |
Science, Faculty of Statistics, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- [if-you-see-this-DO-NOT-CLICK]
- UBC_1986_A6_7 M32.pdf [ 3.71MB ]
- [if-you-see-this-DO-NOT-CLICK]
- [if-you-see-this-DO-NOT-CLICK]
- Metadata
- JSON: 1.0096763.json
- JSON-LD: 1.0096763+ld.json
- RDF/XML (Pretty): 1.0096763.xml
- RDF/JSON: 1.0096763+rdf.json
- Turtle: 1.0096763+rdf-turtle.txt
- N-Triples: 1.0096763+rdf-ntriples.txt
- Original Record: 1.0096763 +original-record.json
- Full Text
- 1.0096763.txt
- Citation
- 1.0096763.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Country | Views | Downloads |
---|---|---|
China | 14 | 30 |
United States | 7 | 1 |
Japan | 4 | 0 |
Germany | 2 | 10 |
Russia | 1 | 0 |
City | Views | Downloads |
---|---|---|
Shenzhen | 10 | 30 |
Tokyo | 4 | 0 |
Beijing | 3 | 0 |
Ashburn | 2 | 0 |
Wilmington | 2 | 0 |
Bamberg | 2 | 0 |
Philadelphia | 2 | 0 |
Sunnyvale | 1 | 0 |
Saint Petersburg | 1 | 0 |
Hangzhou | 1 | 0 |
{[{ mDataHeader[type] }]} | {[{ month[type] }]} | {[{ tData[type] }]} |
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0096763/manifest