## Abstract

The wind speed probability density function (PDF) is used in a variety of applications in meteorology, oceanography, and climatology usually as a dataset comparison tool of a function of a quantity such as momentum flux or wind power density. The wind speed PDF is also a function of measurement scale and sampling error. Thus, quantities derived from a function of the wind PDF estimated from measurements taken at different scales may yield vastly different results. This is particularly true in the assessment of wind power density and studies of model subgrid-scale processes related to surface energy fluxes. This paper presents a method of estimating the PDF of wind speed representing a specific scale, whether that is in time, space, or time–space. The concepts used have been developed in the field of nonlinear geostatistics but have rarely been applied to meteorological problems. The method uses an expansion of orthogonal polynomials that incorporates a scaling parameter whose values can be found from the variance of wind speed at the desired scale. Possible uses of this technique are for scale homogenization of model or satellite datasets used in comparison studies, investigations of subgrid-scale processes for development of parameterization schemes, or wind power density assessment.

## 1. Introduction

The probability density function (PDF) of wind speed is used in many meteorological, oceanographical, and climatological investigations. Meissner et al. (2001) used it to intercompare satellite-sensed ocean surface wind speeds to global surface wind analyses from the National Centers for Environmental Prediction–National Center for Atmospheric Research (NCEP–NCAR) Annual Reanalysis and the European Centre for Medium-Range Weather Forecasts (ECMWF) Tropical Ocean and Global Atmosphere Global Advanced Operational Surface Analysis. Wanninkhof (1992) studied the exchange of gas at the ocean surface as a function of wind speed PDF. Justus et al. (1979) use the wind speed distribution to study the intra-annual variation in wind speed across the United States.

In addition to simple product validation and the assessment of wind climate statistics, the wind speed PDF is also essential to many more physically based investigations. Holland (1973) used the PDFs of wind and temperature data to study turbulent atmospheric eddies at the ocean surface. The ocean surface energy and momentum fluxes depend nonlinearly on wind speed (Wright and Thompson 1983; Wang et al. 1998; Feely et al. 2004). This tends to create a bias in general circulation model (GCM) energy and momentum fluxes when grid-scale winds are used to parameterize them (Vickers and Esbensen 1998; Monahan 2006a, 2007). One way to capture the subgrid-scale wind speed variability is to use its PDF estimated from subgrid-scale winds. Capps and Zender (2008) did just this to establish estimates of the magnitude of this bias in one GCM.

The wind speed PDF is increasingly used in the wind power industry where it is required for the assessment of power potential at different sites (e.g., Hennessey 1977; Garcia-Bustamante et al. 2008; Li and Li 2005; Lackner et al. 2008). The wind power density (WPD) computed at a location in space using measurements taken over time is a function of the wind speed PDF *f* (*W*) (Li and Li 2005); that is,

where *W* is the wind speed (m s^{−1}) and *ρ*_{air} is the air density (kg m^{−3}).

Not always recognized is that the wind speed PDF is a function of the scale of measurements from which it is constructed. The use of time, space, or time–space averaging smoothes out the higher frequencies captured by high-resolution estimates. This results in a decrease in the variance of the estimates, one of the moments of the PDF. Obviously, the use of filters or the interpolation methods has the same effect. This will certainly affect the WPD estimates, because wind turbines require WPD estimates at a point in space. Sometimes, as in the case of wind speed estimates from numerical models, it is unclear exactly what time–space scale an individual grid value truly represents. Numerical models use a combination of parameters, some of which are spatially averaged and some are not. For a discussion of the effects of subgrid-scale processes on spatially average flux parameterization used in numerical models, refer to Mahrt and Sun (1995).

An example of how the use of different wind products representing different scales can contaminate results of a study is illustrated in the work by Monahan (2006b). In this investigation, the first four moments of the wind speed PDF were computed from the following datasets: 1) daily level 3.0 SeaWinds scatterometer data from the National Aeronautics and Space Administration (NASA) Quick Scatterometer (QuikSCAT) satellite, which has a 0.25° spatial resolution and approximately twice a day return period over a particular location; 2) 6-hourly blended Special Sensor Microwave Imager (SSM/I) satellite/model data with a nominal 1° spatial resolution; 3) 6-hourly 40-yr ECMWF Re-Analysis (ERA-40) model output with a 2.5° resolution; 4) 6-hourly NCEP–NCAR reanalysis, which has a nominal 2° spatial resolution; and 5) daily buoy observations from the Tropical Atmosphere Ocean (TAO) and Pilot Research Moored Array (PIRATA) dataset. Significant differences were observed in the wind speed skewness among all five datasets, especially in the tropics (refer to Figs. 5, 6 of Monahan 2006b). Particularly revealing is the very large value for skewness in the buoy observations (Fig. 6 of Monahan 2006b). However, as noted earlier, the PDF is a function of scale, and the point buoy measurements have the largest variance by consequence of their high time resolution and statistically should have the largest skewness value. Unfortunately, all the datasets used in this study represented different time and space scales, which made it unclear whether the differences in the estimated PDFs were due to estimation errors or scale differences in the different products. In fact, Isemer and Hasse (1991) pointed out the importance of scale in the compilation of climate statistics in their study of the different ship-based wind estimation methods.

Most important is the effect of scale differences on physical quantities derived from the PDF. Many of these quantities depend nonlinearly on the wind speed PDF. For instance, note that the WPD in (1) may change drastically depending upon whether a researcher uses 10-min-resolution (Lackner et al. 2008), hourly (Hennessey 1977; Garcia-Bustamante et al. 2008), or daily data (Liu et al. 2008), given that the WPD is a cubed function of the wind speed and wind speed is often positively skewed. The units of the WPD remain the same regardless of the measurement scale. By utilizing low-resolution wind speed measurements, very high wind speed values are not captured by the PDF. Thus, the WPD is likely to be biased low compared to that constructed using high-resolution data. The amount of bias depends upon the magnitude of the spatial correlation structure of the field. A recent example of this can be shown by the study of Liu et al. (2008). In their study, QuikSCAT level 2 (12.5 km, twice a day resolution) satellite estimates of wind were used to produce monthly maps of oceanic WPD. They assumed that wind speed *W* had a Weibull PDF, and they computed the scale *c* and shape *k* parameters using the mean *W* and standard deviation *σ* of the satellite wind estimates; that is,

Given that the satellite spatial resolution is much lower than the point scale, estimates of the standard deviation are biased low to some degree. If the spatial decorrelation distance of the wind field is large relative to the satellite measurement resolution, then the bias will be small because of minor differences in the moments between the point and the spatially averaged data (i.e., point data being quite redundant). If the spatial correlation scale is quite small relative to the satellite measurement scale, significant bias is likely due to large differences in the second- and higher-order moments of the PDFs of measurements at the two resolutions.

The consequence of scaling on the PDF can be shown by a simple example using wind speed data from the Oklahoma Mesonetwork (Brock et al. 1995) site at Hooker, Oklahoma. Constructing hourly and daily wind speed from the basic 5-min-resolution wind speed for the year 2000 and comparing the two histograms (Fig. 1), an obvious difference in shape can be observed. The PDF for both datasets (shown by the smooth curves in Fig. 1) were computed using the Gauss–Hermite method described in appendix B. Both the hourly and daily wind speeds had a mean value of *W* = 5.14 m s^{−1}, whereas the standard deviation values were *σ* = 2.69 m s^{−1} and *σ* = 1.77 m s^{−1}, respectively. If we assume a constant air density of *ρ* = 1.0 kg m^{−3}, then the WPD estimated from (2) for hourly data is 226 W m^{−2}, whereas that using the daily data is 189 W m^{−2} or a change in WPD of approximately 20%. However, although this is an extreme example, it indicates that the scale effect on such quantities cannot be ignored. Also, different physical processes contribute to the measurement statistics differently. For example, using the same data, the PDFs of 5-min and hourly data (not shown) were almost identical, indicating a large amount of redundancy between successive 5-min measurements for this site. Thus, the temporal decorrelation time was quite large. It should also be noted that the 5-min data PDF most likely would be unrepresentative of the “instantaneous” PDF because of a smoothing of the effects of turbulence.

What is required in simple product intercomparison studies or more physically based investigations is the use of a “standard” scale for comparison over the dimensions of interest. For example, a standard time scale in the assessment of the WPD may be that determined by aerodynamic engineers for optimal loads on wind turbine blades. Investigations of the degree of data statistical redundancy would have to be an important component in any study of a standard scale.

The present study presents one method that may be used to approximate a scale-standardized PDF. The technique utilizes concepts from the field of geostatistics to develop a transform function of the wind speed PDF as a function of scale. Using this method and assuming the variance of the wind speed at a given scale is known (or can be estimated), the PDFs representing that scale may be estimated. It can be applied to one or more dimensional problems where data are measured discretely. In simple terms, the PDF from the higher-resolution estimates can be “upscaled” to approximate that of lower-resolution estimates. Thus, the PDFs can be scale corrected, albeit with the caveats discussed below. Admittedly, for wind power assessment, a technique that “downscales” wind estimates is more desirable. However, the method presented here illustrates the complexity involved in scale adjustment and will hopefully motivate researchers to produce downscaling techniques for wind speed and other variables.

The PDF of a continuous random variable *x* has units of the reciprocal of *x*. More importantly, its shape is a function of the basic measurement scale used. This is, of course, a result of the well-known central limit theorem (CLT), which states that a distribution of averages will converge toward normality with an increase in the number of independent samples used to construct the averages. One often overlooked corollary to the CLT is that distributions will also converge toward normality as the scale of the averaging domain increases. Because averages are nothing more than normalized sums, there is an obvious reduction of information when high-resolution measurements are summed to produce an average. This loss of information results from the reduction in the total variance resulting from smoothing of the small-scale fluctuations detectable only in high-resolution measurements. Because variance is one of the moments of the PDF, the PDF must therefore change with measurement scale. The use of low-, band-, or high-pass filters will also reduce the information content of the data, also affecting the PDF, and thus may mask this effect to the unaware researcher. It should be noted that the convergence toward Gaussian is not necessarily monotonic. The CLT says nothing about the rate of convergence of the PDF moments to Gaussian as the scale increases nor does it say anything about how the PDF “shape” evolves as it asymptotically convergences toward a Gaussian PDF. As mentioned above different physical processes generally have different correlation structures. As the averaging scale increases, it is possible, even likely, that the PDF at the averaging scale changes dramatically (i.e., not in a smooth fashion) as the effects of some processes, such as the turbulence, are averaged out.

Often in an attempt to homogenize the scale of different datasets, researchers will utilize simple averaging or a given interpolation routine. If the averages are constructed from data taken discretely, a certain amount of sampling error will exist in the averages, the amount of which depends upon the correlation between the pairs of data points. Sampling error adds to the total variance of the averages. Some amount of sampling error will exist, whether the averages are constructed using simple averaging or interpolation. Because the true mean within a domain is always unknown given discretely sampled data, the sample mean will almost always differ to some degree from the true mean. Thus, the estimated sample variance constructed using the sample mean, not the true mean, will be biased high, thereby biasing the second and higher moments of sample PDF. The effect of sampling error on averages over different dimensions has been extensively studied by many researchers (e.g., Rodríguez-Iturbe and Mejía 1974; Parker 1984; Morrissey et al. 1995). Thus, given the combination of variance reduction resulting from averaging or interpolation and the variance enhancement resulting from sampling error, it can be a very difficult and complex problem to rescale data accurately.

The rescaling method for the wind speed PDF presented here represents one attempt at rescaling. It is applicable to discrete data and for upscaling only. Upscaling is useful for homogenizing the PDF of datasets representing different scales and thus allowing dataset intercomparison without major scale effects. The theory behind the methodology has not been published previously in the meteorology literature. Thus, it is hoped that the main benefit from this paper is an increased awareness of the rescaling problem in tackling a complex statistical issue rather than the presentation of a general method for rescaling. The variable, wind speed, was selected primarily because of its increasing importance in the global renewable energy issues and climate change–related energy flux parameterizations. Another reason wind speed was selected is that it is often unimodally distributed with a moderate positive skew and a small atom (i.e., low frequency of values at a given point on the real axis) at zero (Pavia and O’Brien 1986; Bauer 1996). This allows the wind speed PDF to be easily represented by an expansion of Hermite polynomials (refer to appendix A), probably the most frequently used of the classical orthogonal polynomials. An expansion of Hermite polynomials using only the first term is equivalent to Gaussian PDF (van der Marel and Franx 1993). Adding additional terms to the expansion on the order of three and higher account for deviations from the Gaussian function. Thus, for illustrative purposes alone, wind speed is an appropriate variable, because a Hermite expansion tends to fit moderately skewed PDFs well. The method can be applied to other variables with more complex PDFs, but may require other types of orthogonal functions.

The presentation relies heavily on the theory of isofactorial models developed primarily in the geostatistical community. A summary of the theoretical description of the underlying theory is given for the benefit of those unfamiliar with nonlinear geostatistical concepts. This paper focuses only on the theory and application of the basic Hermitian isofactorial model. For more in-depth descriptions of additional types of isofactorial models, readers are referred to Chiles and Delfiner (1999) and Wackernagel (2003).

The basic methodology involves the development of a PDF transform function, which takes the form of an expansion of orthonormal functions with an embedded scaling parameter. For reasons described below, this transform function is also referred to as an isofactorial model (Chiles and Delfiner 1999; Wackernagel 2003). The description of the technique begins in section 2 with an overview of the theory and the development of the general form of the isofactorial transform function. The modification of the general isofactorial model into an isofactorial change-of-scale model is given in section 3. The scaling parameter used in this model can be found from the variance of wind speed representing a given scale. A method of determining this variance from discrete point values is given in section 4. For illustrative purposes, the technique is applied to a hypothetical two-dimensional domain (e.g., that represented by a satellite wind speed estimate, whose spatial domain encompasses a few scattered surface point observations of wind speed) in section 5. In section 6, an analysis of the behavior of the isofactorial change-of-scale model is undertaken. Finally, section 7 summarizes and comments on the practical uses and limitations of the technique. A review of Hermite polynomials and their statistical properties is given in appendix A. The description of the Hermite expansion is given in appendix B. The equal probability transform and a method of determination of the coefficients needed in the expansion are described in appendix C.

## 2. General theory of isofactorial models for change of scale

### a. Overview

Researchers in the field of geostatistics have been driven by the need to develop methods to estimate the grade of ore in three-dimensional “blocks” of rock from relatively widely spaced core samples. This work was the obvious result of costly decisions of whether to mine a particular block of rock. Matheron (1976) ^{1} realized that the objective was not necessarily to estimate the block ore grade but to estimate the probability that the block’s grade is above a given threshold. This is essentially the same as finding the PDF of the block ore grade. His work and others led to the development of the field of nonlinear geostatistics, of which an integral part includes the development of isofactorial models for “change of support.”^{2} In keeping with meteorological terminology, this paper refers to this as “change of scale.”

Isofactorial models are used for modeling the bivariate distribution between two random variables in terms of their marginals. They are relevant for this study, because we seek an analytical relationship between two random variables: that is, the wind speed estimated at one scale (e.g., the point scale; i.e., *W _{p}*) with that estimated at a larger scale (i.e.,

*W*either in time, space, or space–time). Using an isofactorial model, the PDF of

*W*can be approximated from realizations of Gaussian transformed values of

*W*. The theoretical relations developed below form the basis of isofactorial models, which are used in many different environmental applications (e.g., Rivoirard 1994; Chiles and Delfiner 1999; Wackernagel 2003).

_{p}Assume the existence of two stationary real, continuous random variables, *X* and *Y*, with realizations, *x* and *y*, respectively, which come from a symmetric bivariate distribution *f* (*X*, *Y* ) [i.e., *f* (*X*, *Y* ) = *f* (*Y*, *X* )] with the respective marginals of *X* and *Y*, , and . Following Wackernagel (2003), a bivariate distribution is isofactorial if a given number of orthonormal real functions *χ _{k}* (i.e., “factors”) exist such that

and

where *T _{k}* is the noncentered covariance between the orthonormal functions

*χ*. Because the factors are normalized,

_{k}*T*takes the form of a correlation: that is, −1 ≤

_{k}*T*≤ 1 for Hermite polynomials. The factors

_{k}*χ*are uncorrelated, stationary random functions, and the integrals are taken over the range of support of the functions. The term isofactorial implies that the factors

_{k}*χ*are the same for

_{k}*x*and

*y*. The type of orthonormal functions used can be any one of the many classical orthogonal polynomials or other orthogonal functions, such as the sine and cosine functions. The general form of the isofactorial model shown in (4) is used for many applications, including disjunctive kriging (Chiles and Delfiner 1999), where the bivariate PDF, estimated from the marginals, is used to determine the conditional PDF at unsited locations. This gives the user an idea of the probability of a value surpassing a given threshold at this location. This paper uses a simplified form of the general isofactorial model to approximate the conditional PDF of a variable at one scale from the marginal PDF of a variable at a smaller scale.

It should be noted that nothing proves that the function *T _{k}* exists in nature such that a bivariate PDF can always be factored into its marginals (Chiles and Delfiner 1999). Assuming a bivariate Gaussian PDF does exist, it was proved by Kendall (1941) that it can be factored into Gaussian marginals with

*T*taking the form

_{k}*ρ*[Kendall (1941)’s derivation of the tetrachoric series], where

^{k}*ρ*is the correlation between the two variables. This latter model is referred to as the bi-Gaussian model in the geostatistical literature. Since Kendall’s work, many other non-Gaussian bivariate factorization models now exist that may provide more appropriate fits to certain physical process or processes. Chiles and Delfiner (1999) describe in detail many of these models. It is assumed in this paper that wind speed is a purely “diffusive” process in that the higher-order moments vary smoothly with changes in averaging scale. This may not be the case depending upon the nature of the physical processes that dominate at the different scales of interest. Finally, it will be assumed that a bivariate Gaussian PDF exists whose marginals can be constructed from Gaussian transformed wind speed values. Although potentially restrictive, this assumption is made for illustrative purposes only. As previously mentioned, other models exist that may provide appropriate alternatives if this assumption is violated. The primary purpose of this paper is to illustrate the use of isofactorial modeling for scaling purposes in meteorology. Thus, more complex and probably more appropriate non-Gaussian models are beyond the scope of this introductory paper and will be the subject of future applications.

### b. Derivation of the bi-Gaussian isofactorial model

Because wind speed is often only moderately skewed with a relatively small number of atoms at zero, it is an excellent candidate for a Gaussian transformation (i.e., *W _{p}* ⇔

*X*). Once transformed, the bi-Gaussian model can be applied to the transformed variables. In the present development,

*X*and

*Y*are restricted to be standard bivariate Gaussian distributed (i.e., zero mean, variance equal to one, correlation equal to

*ρ*, and symmetric).

Any function *ϕ* of a variable that is square integrable with respect to a given weighting function [i.e., ∫*ϕ*(*x*)^{2}*f* (*x*)*dx* < ∞, where *f* (*X* ) is the weighting function] can be expanded into a set of orthonormal functions with coefficients *ψ*(*k*) as

Because of the orthogonality of the factors, the coefficients of the expansion can be found by

The conditional expectation of the transform function can be expanded as

and coefficients

Expression (8) now becomes

from (13); the conditional distribution of *X* given *Y* is then

and the joint distribution as a function of the marginals is

Kendall (1941) proved that, if the *f* (*x*, *y*) is the bi-Gaussian PDF, then *T _{k}* takes the form

*ρ*with

^{k}*ρ*being the correlation between

*X*and

*Y*. Expression (12) can be used to expand the conditional expectation of any square integrable function of

*X*[i.e.,

*ϕ*(

*X*)], given random realizations of

*Y*.

## 3. Approximating the lower-resolution wind speed PDF

Assume that the variable *W _{p}* represents discrete high-resolution point wind speed estimates, whereas

*W*represents wind speed averages of

*W*over some domain. For the sake of illustration, assume a homogeneous, square areal domain 𝒟

_{p}^{2}exists where the PDF of

*W*within 𝒟

^{2}is to be estimated from a limited network of point observations within 𝒟

^{2}. Cartier’s relation (Alfsen 1971; Chiles and Delfiner 1999) for a given homogeneous domain is

This simply states that the conditional expectation of randomly located point values in 𝒟^{2} is equivalent to the domain mean. Using the transform model developed in appendix C [i.e., *W _{p}* =

*ϕ*(

*X*)] a new transform function can be derived:

*W*=

*ϕ*

_{𝒟2}(

*Y*). Note that

*X*and

*Y*are two different standard Gaussian variates from the joint Gaussian PDF [i.e.,

*g*(

*X*,

*Y*)] having a correlation

*ρ*, with

*X*representing the Gaussian transform variate for

*W*(i.e.,

_{p}*X*⇔

*W*) and

_{p}*Y*the same for the transformed lower-resolution values (i.e.,

*Y*⇔

*W*). Using Cartier’s relation yields

Given that *E*[*ϕ*(*X* )|*ϕ*_{D2}(*Y* )] = *E*[*ϕ*(*X* )|*Y*], the isofactorial scaling model can be expressed as an expansion of Hermite orthonormal factors by replacing *ϕ*(*X* ) in (17) with its Hermite expansion, , where *η _{k}*(

*X*) are normalized Hermite polynomials on the order of

*k*, which are orthonormal with respect to the standard Gaussian PDF

*g*(

*X*), and

*ψ*(

*k*) are the expansion coefficients on the order of

*k*(refer to appendixes A and B). This substitution yields

From (11), it was shown that

where *T _{k}* =

*E*[

*η*(

_{k}*X*)

*η*(

_{k}*Y*)]. Combining (18) with (19) results in a Hermitian isofactorial model of the form

where the coefficients are

Practical estimates of these coefficients can be calculated from the point wind speed data using the equal probability transform (see Fig. 2) given in appendix C.

In the general situation, the question becomes, what form does the actual covariance between factors representing the point and lower-resolution values take? That is, *T _{k}* =

*E*[

*η*(

_{k}*X*)

*η*(

_{k}*Y*)]. As mentioned above, this can be a complex question, depending up the nature of the physical processes that dominate at given scales. However, this question can be simplified by assuming that the Gaussian transforms of the two wind speed variables have a bivariate Gaussian PDF, then

*T*=

_{k}*ρ*through Kendall’s work (Kendall 1941; this is shown below). Chiles and Delfiner (1999) describe other forms for

^{k}*T*that may provide a better to fit the bivariate PDF of the transform variables, assuming that they have a non-Gaussian bivariate PDF.

_{k}The rationale for using *T _{k}* =

*ρ*arises from Kendall’s tetrachoric series for a Gaussian bivariate factorization, which is

^{k}Because *g*(*x*|*y*) = *g*(*x*, *y*)/*g*(*y*), we have the conditional PDF

and the conditional expectation,

Thus, assuming a bivariate Gaussian PDF for the Gaussian transformed speed, Eq. (24) can be equated to Eq. (19), allowing *T _{k}* =

*ρ*to be substituted into (20), arriving at the transform function for the lower-resolution wind speed values,

^{k}The expected value of (25) is

and the variance of the transform is computed as

using the relations given in appendix A. Thus, we have

The importance of (28) is that the parameter *ρ*, which can be thought of as a “scaling parameter,” can be found from the variance of *W* by simply inverting the relation. A method of determining the *W* variance from point observations is given later. Comparing the expression for the variance of the transform function for given in appendix B [i.e., (B4)] with (28), the only difference between the two relations is the term *ρ*^{2k} in (28).

## 4. The variance of averaged wind speed from discrete point observations

From (28), we have a relationship between the variance of *W* and the scaling parameter *ρ*. It is known from the CLT that the variance of *W* must be less than or equal^{3} to *W _{p}*. The variance of

*W*can be estimated from

*W*values using the method below. Once the variance of

_{p}*W*is obtained, then (28) can be inverted to determine the scaling parameter

*ρ*.

Given a network or time series of *m* discretely taken wind measurements within a one-, two-, or three-dimensional averaging domain (i.e., 𝒟^{1,2or3}), simple averaging of these measurements generally produce high biased estimates of the variance of the resulting averages resulting from sampling error. Given a sample of *m* dependent measurements, the expression for the variance of *W* (Rodríguez-Iturbe and Mejía 1974; Parker 1984; Morrissey et al. 1995; Jones et al. 1997) becomes

where *ρ*_{D1,2,3} is the expected value of the correlation between all points within the averaging domain. For a given one-, two-, or three-dimensional averaging domain, this value is computed using

where *l* is the separation in time, space, or space–time between all pairs of vector locations within 𝒟^{1,2,or3} and *f* (*l*) is the frequency function of all random distances (or time–space separation) within 𝒟^{1,2,or3}. This expression is simply the expected value of the point correlation within the averaging domain. The largest separation distance, time period, or time–space separation measure within the domain is *d*. Analytical expressions for *f* (*l*) are given for one-dimensional lengths and two-dimensional rectangular shapes in appendix D. The frequency function for three-dimensional volumes is given by Zilinskas (2003). Taking the limit as the number of measurements approach infinity [i.e., sampled in a continuous manner (*m* → ∞)], from (29) we obtain the variance of *W* without contamination by sampling error,

Thus, within a stationary and homogenous domain, the variance of *W* can be computed from the product of the variance of *W _{p}* and the domain-averaged correlation, which is used to determine the scaling parameter

*ρ*in (28).

## 5. A worked example

At this point, it would be useful to demonstrate the steps in applying the method above to estimate a low-resolution wind speed PDF from point measurements taken within a hypothetical spatial domain. A realistic analogy to this example would be perhaps the comparison of satellite surface-estimated wind speed to point surface wind speed observations collocated in space and time with the satellite estimates, such as that done by Monahan (2006a). The satellite estimate’s highest spatial resolution may be the NASA QuikSCAT level 3 winds available on a 25 km × 25 km grid, which is to be compared with point wind speed from in situ buoy measurements. In this example, it is assumed that a particular satellite 25 km × 25 km grid square contains *m* = 5 randomly distributed buoy wind measuring sites. Also, suppose that the buoy winds are spatially averaged to obtain grid estimates for comparison to the satellite estimates. The tool for comparison will be the PDF (as in Monahan 2006a). It is supposed that there are 200 buoy wind speed values from the five sites, each measured simultaneously in time with the satellite estimates. This provides *m* × 200 → *n* = 1000 random point data values of wind speed measurements *W _{p}* (nominally in units of m s

^{−1}). Thus, the domain has an area of 625 km

^{2}, with sides equal to 25 km. Because our domain is hypothetical, the buoy data were randomly generated from a Weibull distribution with a shape parameter

*α*equal to 2 and a scale parameter

*β*equal to 6 (Fig. 3). The selected parameter values were chosen to simulate a typical wind distribution with a moderate amount of skewness. Assume that the values drawn from the Weibull distribution represent realizations of the random variable

*W*(i.e.,

_{p}*w*).

The mean, variance, and skewness from the Weibull PDF with the given parameters {*α* → 2, *β* → 6} are

In practice, one would first apply the inverse Gaussian transform function to sorted point measurements, *W _{p}* =

*w*(

*i*),

*i*= 1, 2, 3, … ,

*n*,

*w*(1) <

*w*(2), … ,

*w*(

*n*), to determine the corresponding standard Gaussian values,

*X*=

*x*(

*i*),

*i*= 1, 2, 3, … ,

*n*(as shown in appendix C). The equal probability concept is illustrated graphically in Fig. 1 using the empirical cumulative distribution function (CDF) and the standard Gaussian CDF. A value of

*k*= 15 was initially selected for the number of coefficients to be used in the expansion. The method shown in appendix C was used to compute the coefficients, which were then used to build the transform function, (C5). The transformed values were compared to the simulated

*W*values using a Q–Q plot to determine if an appropriate number of coefficients was selected to adequately model the point PDF. If the fit proved insufficient, then the number of coefficients

_{p}*k*could be increased until a reasonable fit was found. Note that 1000 standard Gaussian random values were used in transform function to produce 1000 transformed values. The resulting Q–Q plot (Fig. 4) shows a very good fit and, thus

*k*= 15 was deemed sufficiently large for this example. A Q–Q plot is a good way to determine if a sufficient number of terms have been included to produce an adequate model. For example, the Q–Q plot in Fig. 5 resulted from incorporating only the first three terms in the expansion. Obviously, the tails of the distribution are inadequately represented by the expansion.

The next step would be to use the point wind data to estimate the spatial correlation function. Because this worked example is hypothetical, the following exponential function was arbitrarily selected to represent an isotropic spatial correlation function with an *e*-folding distance of 25 km; that is,

The expected value of the correlation was determined by numerically integrating Eq. (30) using the frequency function in appendix D [i.e., (D3)] over all separation distances *l* within the 25 km × 25 km domain. This value turned out to be *ρ*_{D2} = 0.61, which when substituted in Eq. (31) and using the computed point variance of the *n* values produced a domain variance of .

Given this value, expression (28) was used to determine the value of the scaling parameter *ρ*. A plot of versus *ρ* is shown in Fig. 6. By inverting Eq. (28), a scaling parameter of *ρ* = 0.79 was found. This value corresponds with that estimated graphically using Fig. 6. Note that the *ρ* value corresponding to the point variance of *W _{p}* (i.e., 7.37 m

^{2}s

^{−2}) is approximately 1.0, as it should be because (25) will then reduce to the expression for the point transform function.

The scaling parameter can now be inserted into Eq. (25) to produce any number of simulated satellite values. Again 1000 standard Gaussian random numbers (*Y* = *y*(*i*), *i* = 1, 2, 3, … , 1000) were used in Eq. (25) to produce 1000 simulated values [i.e., *W* = *w*(*i*), *i* = 1, 2, 3, … , 1000]. The resulting histogram of these values is shown in Fig. 7. Table 1 gives the statistics for *W _{p}* and the transformed

*W*values. The mean of the simulated

*W*values is approximately equal to that of the

*W*values, and the

_{p}*W*variance is equal to the analytically determined value . In Fig. 7 note that the

*W*values have a distribution whose shape is considerably more Gaussian in appearance than that for the

*W*values (Fig. 7). The skewness computed from the 1000 transformed values decreased from 0.60 to 0.37 (Table 1). The lack of any significant difference between the first two moments of the transformed values and their theoretically expected values indicated that the number of Gaussian values used in the transform function was large enough to avoid excessive random noise in the resulting

_{p}*W*histogram and the estimated moments.

We now have a set of values whose distribution approximates the true *W* histogram. It should be noted that we use the word “approximates” only in the sense that the moments greater than the second behave as expected from the CLT, given a single diffusive physical process. We could now apply any number of parametric or nonparameter methods to estimate the *W* PDF from its histogram. In keeping with the theme of the paper, a nonparameter method using Hermite polynomials will be used.

The nonparametric method of choice for estimating the PDF using orthogonal functions is the Gauss–Hermite expansion (van der Marel and Franx 1993). It has been demonstrated that this expansion has much better convergence properties than that of the Gram–Charlier series A expansion (Blinnikov and Moessner 1998). The Gauss–Hermite expansion requires a slightly different form of the Hermite polynomial; that is,

and is expressed by

where *x* represents the scaled values of *W*; *N* + 1 is the number of coefficients selected; and *μ* and *σ* are its mean and standard deviation of *W*, respectively. The Hermite polynomials defined in (36) have the orthogonal relation

The coefficients in (37) are determined from

Because *f* (*w*) is a PDF, the coefficients Ψ(*k*) are equal to *E*[*Hp _{k}*(

*x*)

*g*(

*x*)], which can be estimated from scaled data,

*x*(

*i*) = [

*w*(

*i*) −

*μ*]/

*σ*,

*i*= 1, 2, 3, … ,

*n*,

*n*= 1000, using

Equation (37) produces a nonparametric estimate of the *W* PDF, which is plotted together with the *W* histogram in Fig. 7. For comparison, an expansion was also fit to the *W _{p}* values and plotted, together with its histogram, in Fig. 8.

When using the above nonparameter technique for fitting a PDF, a number of factors need to be considered. The amount of smoothing of the resulting PDF is controlled by the number of coefficients used. The lower the number of coefficients, the smoother the function becomes. A small number of randomly selected Gaussian values used in the transform function results in a very noisy histogram. So, it is recommended to use a very large number of values and try different numbers of expansion coefficients to find a PDF that provides a reasonable fit to the histogram. More details can be found about nonparametric density estimation in using kernels in Silverman (1998) and polynomial orthogonal functions in Efromovich (1999).

As can be observed in Figs. 7 and 8, the transform function produced a near-Gaussian PDF and a shape that is generally consistent with the CLT given diffusive processes. Notice in Fig. 7 that the function produces small negative values for the PDF in the right tail. This is an inherent problem with orthogonal polynomial fits. The negative values are often small and can be set to zero as one option (Efromovich 1999). Integrating the function over the range of support (−∞, ∞) in Fig. 6 produced a value of 1.0, so the negative values had minimal effect of the PDF estimate in this study. It should be noted that negative density values resulting from polynomial fits can have serious effects on the Gaussian transform, because the cumulative distribution function is no longer monotonically increasing. This will negate the required bijective capability of the transform. Thus, it is imperative that an expansion producing negative density values be corrected in some way. One may follow the corrections suggested by Efromovich (1999). We have found that adjusting the number of terms in the expansion works in many situations. Nevertheless, numerical integration of the resulting function in Fig. 7 produced values for the first three moments that were extremely close to those shown in Table 1. Another experiment using a *ρ* value of 0.9 produced a PDF with a skewness in between the PDFs for *W* and *W _{p}*, which is also consistent with the CLT for diffusive processes (not shown).

## 6. Analysis of the transform equation

One intriguing question pertaining to the method presented in this paper is how the isofactorial transform function behaves with decreasing variance [or equivalently, to decreasing values of the scaling parameter *ρ*; refer to the function (28)].

To examine this behavior closer, the coefficients computed in the worked example given in section 5 are used to observe the variation in skewness with different values of the scaling parameter *ρ*. The skewness of the transformed values can be expressed in terms of the expansion shown in (25) as a function of *ρ* as

where *μ*_{3} is the third central moment; *nk* is the coefficient cut-off limit for the expansion, selected as 15 for this experiment; and *n* is the number of standard Gaussian random numbers used (*n* = 1000). Given that *ρ* is an increasing function of variance (Fig. 6), the variance and *ρ* are de facto measures of changing scale (all things else being equal). Thus, a decreasing value of |*ρ*| ≤ 1 is consistent with an increase in the scale of the domain. By observing the variation in the skewness as a function of *ρ*, the convergence of the PDF toward Gaussian can be studied. With the coefficients, *ψ*(*k*) computed in section 5, the function (41) is plotted against *ρ* in Fig. 9. Notice that the relationship between *ρ* and skewness, although close, is not linear. It should be noted that a *ρ* value of zero implies a Gaussian PDF, because skewness approaches zero as *ρ* → 0 (assuming that the higher-order moments also approach zero).

It is of interest to observe how the transformed values change with changing *ρ* values. Figure 10 is a contour plot of the transform function (25) with different *ρ* values as a function of standard Gaussian variates ranging from −3 to 3. A standard Gaussian PDF is superimposed on this plot for reference purposes. If a horizontal line is taken across the figure at any given value of the scaling parameter (e.g., lines A and B in Fig. 10), this provides an idea of the shape of the transformed PDF by observing the change in the values along a given line. Thus, moving this line up and down shows how the transform PDF changes with different values of the scaling parameter. A convergence toward normality can easily be observed with decreasing values of *ρ* as the contour lines become increasingly more symmetric about the Gaussian mean of zero. Given a *ρ* value of 0.2 and Gaussian values ranging from 3 to −3, there is a very small probability of obtaining a transform value larger than 7.0 m s^{−1} or less than 4.0 m s^{−1}. The mean, 5.3 m s^{−1}, is nearly halfway between these two values, indicating near symmetry. For a *ρ* value of 0.8, the respective transformed values are 0.5 and 13.0 m s^{−1}. The differences between these transformed values and the mean are −4.8 and 7.7 m s^{−1}, respectively. As the scale decreases (i.e., *ρ* increases), the probability of obtaining larger values increases dramatically. It is highly unlikely to obtain values over 7.0 m s^{−1} with *ρ* less than 2.0. This is, of course, a result of the positive skew in the parent PDF associated with the *W _{p}* variable. Note that the expected value of transformed value (i.e., ∼5.3 m s

^{−1}) is independent of scale as expected by Cartier’s relation.

On the left side of the plot, the values also converge with increasing scale, albeit at a slower rate. Given that wind speed is always greater than or equal to zero, a convergence rate of equal magnitude on both sides of the mean would not be consistent with the CLT and would result in a shape other than Gaussian. Thus, the behavior of the transform function with changing scale is consistent with the CLT and by extension what is usually observed in nature.

## 7. Summary

The primary purpose of this paper is to introduce the isofactorial method of estimating a PDF from point wind measurements averaged over a domain given a certain amount of sampling error. A comprehensive treatment of isofactorial models is beyond the scope of this paper. The reader is referred to Chiles and Delfiner (1999) for more details and other applications of isofactorial models.

As mentioned several times above, the development and application of the bi-Gaussian isofactorial model for scaling is not without assumptions for proper application. First, it is assumed that not only are the marginals of the Gaussian transforms of *X* = *G*^{−1}(*F*(*W _{p}*)) and

*Y*=

*G*

^{−1}(

*F*(

_{ν}*W*)) Gaussian [where

*F*(

_{ν}*W*) is the CDF of

*W*], but their bivariate distribution is Gaussian as well. The first assumption is easier to check than the second. For example, it might reasonably be assumed that 0.1-s wind speed measurements are dominated by the effects of turbulence, whereas daily wind speed has the turbulence effects averaged out and are influenced primarily by synoptic-scale systems, seasonal effects, etc. Assuming the statistics on the turbulence scale are considerably different from those at the daily scale, the use of

*T*=

_{k}*ρ*will still produce a bi-Gaussian PDF from the marginals of the transformed data, but the actual joint PDF may be quite different. The consequences will be that the use of

^{k}*T*=

_{k}*ρ*in (25) would be invalid and the simulated

^{k}*W*values would have a PDF different to some degree from reality. Different forms for

*T*could be tried, and the resulting scatterplot could be compared with observations.

_{k}In addition, the influence of different physical processes at different scales affect the way the PDF convergences with averaging as it asymptomatically approaches Gaussian. It is quite possible for the unimodel PDF to become bimodel then back to unimodel as different physical effects are being averaged out. It should also be reemphasized that the main problem at hand, determining the lower-scale PDF, actually has no exact solution. The true PDF is a function of an infinite number of mixed moments.

Given these caveats, one may question the worth of isofactorial scaling models. It turns out that, in most situations, where the scales are not vastly different and the physical processes at these scales have statistical properties, which vary smoothly with averaging, an isofactorial scaling model will provide an adequate approximation of the PDF at the lower scale. The first two moments will be exact and the higher-order moments will approach zero. This behavior is commonly seen in nature with increasing scale with the results shown in Fig. 1 being one such example.

For applications of the bi-Gaussian isofactorial model, the transformation of non-Gaussian data values to Gaussian values should be done only when the data values have a slightly skewed distribution with no large number of atoms at a given value. Thus, the technique given in this paper is not recommended for application to high-resolution rainfall data, which are highly skewed and have many values at zero. In this situation, it may be reasonable to transform observed rainfall into values with a gamma PDF and apply a transformation using Laguerre polynomials, which utilized a gamma PDF for its weighting function (Chiles and Delfiner 1999; Wackernagel 2003). Although the example shown in this paper dealt with the PDF of averages of point values within a spatial domain, it can also be applied to a temporal series whose measurements are taken discretely in time as well as space–time averaging domains, such as those found in box averages of monthly ship data (Morrissey and Greene 2008).

An interesting question remains as to how well the method reproduces the higher-ordered moments of a PDF. The answer obviously depends upon the data and its mixed moments, which are usually unknown. However, it is known that increasing scale produces PDFs tending toward Gaussian, given the caveats discussed earlier. Without any additional information, this method at least produces a first-order approximation of a lower-resolution PDF that is consistent with the CLT’s contention that a PDF converges asymptotically toward Gaussian with increasing scale, albeit not necessarily monotonically. Further research is needed to find ways to incorporate the mixed moments into the analytical transform function in a similar manner that the second moment is currently utilized.

The method presented here may prove useful to those wishing homogenize different wind datasets so that they all represent the same scale, especially when using the wind speed PDF for comparison. Given that the WPD is a strong nonlinear function of the PDF, this study demonstrates the need to incorporate scaling considerations into its estimation.

## Acknowledgments

The authors would like to acknowledge Grant 105128600 from the Oklahoma EDGE (Economic Development Generating Excellence) initiative for funding this work. Also, thanks goes to the Oklahoma Mesonetwork for providing the wind speed data from Hooker, Oklahoma. A special thanks goes to Ethan Cook, who provided needed “thinking” on Gauss–Hermite expansion.

## REFERENCES

**,**

**,**

**,**

**,**

_{2}fluxes in the equatorial Pacific Ocean.

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

### APPENDIX A

#### Hermite Polynomials

Geostatisticians generally define Hermite polynomials as a function of the standard Gaussian PDF,

Hermite polynomials of different order *k* can be produced using a recurrence relation,

They are orthogonal with respect to the standard Gaussian PDF, *g*(*x*), *x* ∈ *R*, in the sense

Orthonormal Hermite polynomials are defined by *η _{k}*(

*x*) = [

*H*(

_{k}*x*)]/

*k*!, where

*k*! is called the “norm,” and thus

Noting that *η*_{0}(*x*) = 1, the first few moments of Hermite polynomials are

### APPENDIX B

#### Hermite Expansions

Any function of a standard Gaussian random variable *ϕ*(*X* ) can be expanded into a series of Hermite polynomials if its second moment is finite: that is, *E*[*ϕ*^{2}(*x*)] < ∞. Using

with coefficients

from (B2), it can be observed that

The first raw moment and second central moment of the expansion are

When expanding a function using (B1), the order of the expansion *k* is usually truncated at a point where the individual contributions of additional terms become negligible. This is usually chosen at a *k* value where 10 < *k* ≤ 30.

### APPENDIX C

#### The Equal Probability Transform and the Determination of the Hermite Coefficients

The Hermitian isofactorial model requires standard Gaussian variates for input. A transformation of the non-Gaussian distributed wind speed values into standard Gaussian variates (i.e., *X* ) is necessary so that the change-of-scale isofactorial model can be applied. Following Journel and Huijbregts (1978), this is done using an equal probability transform. Assume a model exists where

where *w*(*i*), *i* = 1, 2, 3, … , *n* is the available set of *n* point wind speed measurements taken at random locations within, say, a grid-square domain 𝒟^{2} (e.g., 5 point sites and 200 temporal observations producing *n* = 5 × 200 = 1000 values). Let *x*(*i*), *i* = 1, 2, 3, … , *n* be the set of *n* standard Gaussian values corresponding in probability with the wind speed values [i.e., *P*(*W _{p}* <

*w*(

*i*)) =

*P*(

*X*<

*x*(

*i*)) for each

*i*]. The function

*ϕ*(

*X*) is bijective in the sense that there is a one-to-one correspondence in probability between

*w*(

*i*) and

*x*(

*i*). Let

*F*(

*w*) =

*P*(

*W*<

_{p}*w*) be the cumulative PDF (i.e., CDF) of point wind speed and

*G*(

*x*) =

*P*(

*X*<

*x*) be the standard Gaussian CDF. Assuming the inverse of

*F*(

*w*) exists, an equal probability transformation function can then be constructed as

and conversely as

By sorting the wind speed observations *w*(*i*), *i* = 1, 2, 3, … , *n* in ascending order, an estimate of the empirical cumulative distribution function of *W _{p}* can be obtained using

where 1_{Wp<w(i)} means that, for *W _{p}* <

*w*(

*i*), 1

_{Wp<w(i)}= 1 and 0 otherwise. The standard Gaussian value

*x*(

*i*) corresponding to a given sorted

*w*(

*i*) value can be found by mapping the inverse error function to the sorted

*w*(

*i*) values or using the approximation given by Abramowitz and Stegun (1972, p. 933). This can also be done by using a simple graphical method (Fig. 2).

The Gaussian transform function is of the form

with coefficients equal to

Given the *n* pairs of point wind speed and standard Gaussian values, {*w*(*i*), *x*(*i*), *i* = 1, 2, 3, … , *n*}, the coefficients can be approximated as an integral of a sum of the *n* elements (for a more detailed derivation, see Ortiz et al. 2005); that is,

with *g*(*x*(1)) = *g*(−∞) = 0 and *g*(*x*(*n*)) = *g*(∞) = 0.

Substituting these coefficients into (C5) produces a transform function where standard random Gaussian values can be used to produce values with the same PDF as the point wind speed data. By applying the transform to a large number of standard random Gaussian values, a pool of *w* = *ϕ*(*x*) samples is available, where a nonparametric technique can be applied to estimate the *W _{p}* PDF.

### APPENDIX D

#### The Frequency Functions between Two Points in One- and Two-Dimensional Domains

From Rozanov (1997), the cumulative distribution of the time separation *t* between any two instances in time within a given time interval *L* equals

If we take the derivative of (D1) with respect to *t*, we arrive at the PDF of *t* within *L* as

The distribution of the distance between two points chosen at random in a plane convex region as given by Ghosh (1951) is given by first noting that the domain 𝒟^{2} is a rectangle of area *A* with sides *A*_{1} and *A*_{2}. The PDF of the distance *l* between two points in 𝒟^{2} is

## Footnotes

*Corresponding author address:* Mark L. Morrissey, School of Meteorology, University of Oklahoma, 120 David L. Boren Blvd., Suite 5900, Norman, OK 73072. Email: mmorriss@ou.edu

^{1}

Georges Matheron is known as the father of geostatistics. As an engineer and mathematician, he generated a plethora of papers, mostly in French, during the 1960s and 1970s in this field. For a concise compilation of the results of his work and the work of others in this field, refer to Chiles and Delfiner (1999).

^{2}

This is terminology from the geostatistical literature, where “support” refers to “scale.”

^{3}

If the correlation everywhere in a given domain is 1.0, then the *W* and *W _{p}* variances will equal each other.