to Q NBS TECHNICAL NOTE

U.S. DEPARTMENT OF COMMERCE/ National Bureau of Standards

Statistical Analysis o

Extreme Winds

QC

100

U5753

no ,861

1975
NATIONAL BUREAU OF STANDARDS

The National Bureau of Standards 1 was established by an act of Congress March 3, 1901. The Bureau's overall goal is to strengthen and advance the Nation's science and technology

and facilitate their effective application for public benefit. To this end, the Bureau conducts

research and provides: (1) a basis for the Nation's physical measurement system, (2) scientific and technological services for industry and government, (3) a technical basis for equity in trade, and (4) technical services to promote public safety. The Bureau consists of the Institute for Basic Standards, the Institute for Materials Research, the Institute for Applied Technology,

the Institute for Computer Sciences and Technology, and the Office for Information Programs.

THE INSTITUTE FOR BASIC STANDARDS provides the central basis within the United

States of a complete and consistent system of physical measurement; coordinates that system

with measurement systems of other nations; and furnishes essential services leading to accurate and uniform physical measurements throughout the Nation's scientific community, industry, and commerce. The Institute consists of a Center for Radiation Research, an Office of Measurement Services and the following divisions: Applied Mathematics — Electricity — Mechanics — Heat — Optical Physicc — Nuclear

Sciences " — Applied Radiation - — Quantum Electronics " — Electromagnetics 3 — Time

and Frequency * — Laboratory Astrophysics 3 — Cryogenics ". THE INSTITUTE FOR MATERIALS RESEARCH conducts materials research leading to improved methods of measurement, standards, and data on the properties of well-characterized

materials needed by industry, commerce, educational institutions, and Government; provides

advisory and research services to other Government agencies; and develops, produces, and

distributes standard reference materials. The Institute consists of the Office of Standard Reference Materials and the following divisions: Analytical Chemistry — Polymers — Metallurgy — Inorganic Materials — Reactor

Radiation — Physical Chemistry.

THE INSTITUTE FOR APPLIED TECHNOLOGY provides technical services to promote

the use of available technology and to facilitate technological innovation in industry and

Government; cooperates with public and private organizations leading to the development of

technological standards (including mandatory safety standards), codes and methods of test; and provides technical advice and services to Government agencies upon request. The Institute consists of a Center for Building Technology and the following divisions and offices: Engineering and Product Standards — Weights and Measures — Invention and Innovation — Product Evaluation Technology — Electronic Technology — Technical Analysis — Measurement Engineering — Structures, Materials, and Life Safety 4 — Building

Environment 4 — Technical Evaluation and Application 4 — Fire Technology.

THE INSTITUTE FOR COMPUTER SCIENCES AND TECHNOLOGY conducts research and provides technical services designed to aid Government agencies in improving cost effec- tiveness in the conduct of their programs through the selection, acquisition, and effective

utilization of automatic data processing equipment; and serves as the principal focus within

the executive branch for the development of Federal standards for automatic data processing

equipment, techniques, and computer languages. The Institute consists of the following

divisions: Computer Services — Systems and Software — Computer Systems Engineering — Information Technology.

THE OFFICE FOR INFORMATION PROGRAMS promotes optimum dissemination and

accessibility of scientific information generated within NBS and other agencies of the Federal Government; promotes the development of the National Standard Reference Data System and

a system of information analysis centers dealing with the broader aspects of the National Measurement System; provides appropriate services to ensure that the NBS staff has optimum

accessibility to the scientific information of the world. The Office consists of the following

organizational units: Office of Standard Reference Data — Office of Information Activities — Office of Technical

Publications — Library — Office of International Relations.

1 Headquarters and Laboratories at Gaithersburg, Maryland, unless otherwise noted; mailing address Washington, D.C. 20234.

2 Part of the Center for Radiation Research.

3 Located at Boulder, Colorado 80302.

4 Part of the Center for Building Technology.
Statistical Analysis of Extreme Winds

^-LetKvco.l MNo-Ve } uo gg&

Emil Simiu

Institute for Applied Technology

National Bureau of Standards

Washington, D.C. 20234

and

James J. Filliben

Institute for Basic Standards

National Bureau of Standards

Washington, D.C. 20234

*z%

»< °r

'«, \

in

U.S. DEPARTMENT OF COMMERCE, Rogers C. B. Morton, Secretary

\\ £ NATIONAL BUREAU OF STANDARDS, Richard W. Roberts, Director

Issued June 1975
Library of Congress Catalog Card Number: 75-600028

National Bureau of Standards Technical Note 868

Nat. Bur. Stand. (U.S.), Tech. Note 868, 52 pages (June 1975) CODEN: NBTNAE

U.S. GOVERNMENT PRINTING OFFICE

WASHINGTON: 1975

For sale by the Superintendent of Documents, U.S. Government Printing Office, Washington, D.G 20402 (Order by SD Catalog No. C13.46:868). Price $1.50 (Add 25 percent additional for other than U.S. mailing).
TABLE OF CONTENTS

Page

Abstract and Key Words 1

SI Conversion Units 1

List of Symbols 2

1. Introduction 3

2. Probability Distributions of the Largest Yearly Wind Speed 4

3. Description of Procedure for Estimating CDF, its Parameters, and the N-year Winds. ... 5

4. Summary and Interpretation of Results 8

5. Conclusions 11

6. References 12

Appendix: Computer Program Listing. Sample Input and Output 20

iii

STATISTICAL ANALYSIS OF EXTREME WINDS

Emil Simiu

James J. Filliben

With a view to assessing the validity of current probabilistic approaches to the

definition of design wind speeds, a study was undertaken of extreme wind speeds

based on records taken at 21 U.S. weather stations. For the purpose of analyzing

extreme value data, a computer program was developed which is described herein.

The following results were obtained: (1) the assumption that a single probability

distribution is universally applicable to all extreme wind data sets in a given

type of climate was not confirmed, and (2) predictions of 100-year wind speeds

based on overlapping 20-year sets of data taken at the same station differed

between themselves by as much as 100%. Similar predictions for 1000-year winds

differed by as much as a few hundred percent.^ Since wind pressures are proportional

to the square of the wind speeds, errors of such magnitude are unacceptably high

for structural design purposes. It is therefore suggested that while, in principle,

probabilistic methods provide the most rational approach to specifying design wind

speeds, it is of the utmost importance that the possible errors inherent in this

approach be carefully taken into account.

Key words: Building codes; extreme value distributions; hurricanes; probability

distribution functions; reliability; risk; statistical analysis; storms; structural

engineering; wind loads; wind speeds.

SI CONVERSION UNITS

In view of the present practice in building technology in the United States and in

publications of the National Oceanic and Atmospheric Administration (NOAA) , common U.S.

units of measurements have been used in this paper. However, in recognition of the position

of the United States as a signatory to the General Conference on Weights and Measures, which

gave official status to the metric SI system of units in 1960, conversion factors are given

as follows:

Length

1 inch (in) = 0.0254* meter (m)

1 foot (ft) = 0.3048* meter (m)

1 mile (U.S. Statute) = 1.609344 x 10

3 meter (m)

Velocity

1 mile per hour (mph) = 4.470400 x 10 meters per second (m/s)

= 1.609344 kilometers per hour (km/hr)

* exactly

i
LIST OF SYMBOLS

CDF Cumulative distribution function

D, D' Probability distributions

F(v) Mixed distribution given by Eq. 6

F Specified value of the cumulative distribution function

o

f (v),F (v) Extreme type I and type II cumulative distribution function

E T

F (v) ,F (v) Extreme type II cumulative distribution function for extratropical and tropical storm winds, respectively

F (x) Extreme type II cumulative distribution function with tail length parameter y Y

G (p) Percentage point function given by Eq. 5, corresponding to distribution of

Y random variable X

M.(D) Median of i-th ordered observation from a sample of size n from a distribution D|

M(D) £M.(D)/n

n Sample size (number of observations)

N Mean recurrence interval, in years

p Value of cumulative distribution function

'E' FT

v, Vmax

oo VN

,VN

V X

X.

1 X

y, Y

'opt

M

a

Probability of largest yearly wind being produced by an extratropical and by

a tropical storm, respectively

PPF Percentage point function

r Probability plot correlation coefficient defined by Eq. 7

SD(X) Standard deviation of variable X

s Standard deviation of the observed annual extreme wind speeds

vv Largest yearly wind speed

Mean, maximum value of the observed annual extreme wind speeds, respectively

Extreme wind speed corresponding to a N-year mean recurrence interval obtained

using a distribution with Y

= Y » Y= °°» respectively

Random variable, value taken on by X.

i-th ordered observation

IX. /n

x

Tail length parameter, optimal value of v, respectively

Location parameter

Probability distribution scale parameter

2
The Statistical Analysis of Extreme Winds

1. INTRODUCTION

In modern building codes and standards [1,3] basic design wind speeds are specified in

probabilistic terms. At any given station at which wind records over a number of years are available, a random variable may be defined, which consists of the largest yearly wind speed

for every year of record. Using these records the cumulative distribution function (CDF) of

this random variable may, at least in theory, be estimated to characterize the probabilistic

behavior of the largest yearly wind speeds. The basic design wind speed is then defined as the speed corresponding to a specified value F of the CDF or, equivalently (in view of the

relation N = 1/(1-F ), in which N = mean recurrence interval), as the speed corresponding to a specified mean recurrence interval. For example, the American National Standard A58.1 [1]

specifies that a basic design wind speed corresponding to a 50-year mean recurrence interval

(i.e., to a value F of the CDF equal to 0.98, or to a probability of exceedance of the

basic wind speed in any one year equal to 0.02) be used in designing all permanent structures, except structures with an unusually high degree of hazard to life and property in case of failure, for which a 100-year mean recurrence interval (F = 0.96) must be used, and

structures having no human occupants or where there is negligible risk to human life, for

which a 25-year mean recurrence (F = 0.96) may be used. A wind speed corresponding to an N-year recurrence interval is commonly referred to as the N-year wind.

The mean recurrence intervals specified by building codes, rather than being based on a formal risk analysis—which is in practice not feasible in the present state of the art —are

selected in such a manner as to yield basic wind speeds which, by professional consensus,

[are judged to be adequate from a structural safety viewpoint. Nevertheless, it is generally

assumed that the current probabilistic approach to the definition of design wind speeds

insures, at least in theory, a certain degree of consistency with regard to the effect of the

wind loads upon structural safety; i.e., all other relevant factors being equal, if appropriate recurrence intervals are used in design, the probabilities of failure of buildings in

'different wind climates will, on the average, be approximately the same.

I

In the practical application of the probabilistic definition of design wind speeds,

certain important questions arise. One such question pertains to the type of probability

distribution of the National Building Code of Canada [3] are based upon the assumption that

this behavior is best modeled by a Type I (Gumbel) distribution. The American National

Standard A58.1 [1], on the other hand, assumes that the appropriate models are Type II (Frechet) distributions with location parameters equal to zero and with tail length parameters dependent only upon type of storm.
A second important question is whether records of approximately 20-year length, i.e.,

of such length as has been used in developing wind intensity maps in the American National

Standard A58.1, are sufficient for making reliable predictions of extreme wind speeds.

The present work, which is part of an effort to evaluate and improve building code

provisions on design for wind, was undertaken with the intent of seeking an answer to these

two questions. A computer technique was developed for estimating the parameters of the

probability distribution function of the largest values which best fits any given set of

extreme wind speed data; using this program, an analysis was carried out of extreme winds

recorded at 21 U.S. weather stations and published by Court [5]. The data consisted of

5-minute averages of the largest yearly wind speeds recorded during 37 consecutive years.

All the data were obtained at stations where no change in the height and the exposure of the

wind recording instruments was noted throughout the period of record. The results of the

analyses are presented and, on their basis, answers to the two questions previously mentioned are suggested. These results point to the need for carefully taking into account the

possible errors inherent in the probabilistic approach to the definition of wind speeds for

purposes of structural design.

2. PROBABILITY DISTRIBUTIONS OF THE LARGEST YEARLY WIND SPEED

Probabilistic considerations [11, pp. 274-275], as well as available empirical evidence

[6, 17, 18] suggest that the asymptotic probability distributions of the largest values with

unlimited upper tail are an appropriate model for the behavior of the largest yearly wind

speed. There are two such distributions, known as the type I and type II distributions of

the largest values [11], whose cumulative distributions functions, F (v) and F (v) , re- spectively, are of the form

F (v) = exp {-exp[-(v-y)/a]

}

y < v < °°

-OO < |j < oo (1)

< a < °°

and

F (v) = exp {-[v-y)/a]~Y

} y < v < °° II

—00 < y < 00

-OO < (J < oo

Y >

(2)

in which y, a, and y are location, scale and tail length parameters, respectively. Actually,

the type I distribution may be shown to be a type II distribution with y = °°; however, it is

convenient to refer to it separately.

It is convenient in many applications to use the inverse function of the CDF, known as

the percent point function (PPF) . For Eqs. 1 and 2, the PPF's are, respectively,
< F < 1 (3)

< Fn< 1 (4)

v(F ) = y - a £n (-Jin F )

v(F

I]:

) = v + o (-In F

II

)"1/Y

It is customary to denote the CDF value F XXX or F as p and v(F) = G (p) . With these notaA

tions, Eq. 4 becomes

Gx

(p) = u + a (-£n p)"1/Y < p < 1 (5) Y The probability distribution of the largest value depends upon the form of the underlying (or initial) distribution, i.e., the distribution of the parent population of wind

speeds from which the largest values have been extracted. The underlying distribution is of

the exponential type if its CDF converges toward unity with increasing value of the variate

as fast as, or faster than, the CDF of the exponential distribution; otherwise, it is said

to be of the Cauchy type. Under the assumption of statistical independence, it can be shown

that, asymptotically, i.e., for increasingly larger sample sizes, the largest sample value

from an exponential and from a Cauchy type distributions have type I and type II distributions, respectively. Since there is empirical evidence to the effect that its parent population (say, the largest weekly wind speed) appears to follow a Rayleigh distribution, which

is of the exponential type, it has been argued that the largest yearly wind speed should

follow a type I distribution [6]. According to Thorn [16, 17, 18, 19], however, the largest

yearly wind speed follows type II distributions with location parameter \i = and with tail

length parameters y - 9.0 and y - 4.5 for winds associated with extratropical storms and

with tropical storms, respectively. At locations at which both types of storms occur, Thorn assumes that a mixed distribution holds,

F(v) = PE

F

E IZ

(v) + P T

FT

I]: (v) (6)

E T

in which F (v) , F (v) are two type II CDF's for extratropical and for tropical storm

winds, p = probability of largest yearly wind being produced by an extratropical storm (or

the proportion of extratropical storm extreme winds) and p = 1 - p . T E

3. DESCRIPTION OF PROCEDURE FOR ESTIMATING CDF, ITS

PARAMETERS AND THE N-YEAR WINDS

The purpose of this section is to describe the computer technique which was utilized

to estimate the CDF, the values of its parameters and the corresponding extreme wind speeds

I (i.e. wind speeds with given probabilities of being exceeded in any one year). The input to

I this procedure is the observed set of annual wind speeds from a given station. Based on any

given set of observed annual wind speeds, the principal output from this procedure is the

estimated wind speeds vrr for various mean recurrence intervals. In this study N = 50, 100,

500 and 1000 years were used.
The procedure consists of 3 distinct stages. In the first stage t^e value of y (Eqs.

2 and 4) is determined which yields the closest fit to the observed data set (recall that

Y = oo corresponds to an extreme value type I distribution) . The "closest fit" criterion

used in this stage is the so-called maximum probability plot correlation coefficient criterion [10]. The probability plot correlation coefficient is defined as

E(X - X) [M .(D) - M(D)]

r = Corr(X.M) — , (7)

[E(X - X)

Z E(M (D) - M(D))

Z

] ±/ ^

in which X = EX./n, M(D) = IM.(D)/n, n = sample size, D = probability distribution tested.

The quantities X. are obtained by ?. rearrangement of the data set: X is the smallest, X„ the

second smallest, X. the i-th smallest of the observations in the set. The quantities M.(D) are obtained as follows. Given a random variable X with probability distribution D and

given an integer sample size n, it is possible from probabilistic considerations, to derive

mathematically the distributions of the smallest, second smallest, and in general the i-th

smallest values of X in a sample of size n. There are various quantities that can be utilized to measure the location of the distribution of the i-th smallest value X. (e.g., the

mean, the median or the mode). As shown in Ref. 10, it is convenient to use the median as a measure of location in Eq. 7 - these medians of the distribution of the i-th smallest value

being denoted by M.(D).

If the d^ta set was generated by the distribution D, then aside from a location and

scale factor, X. will be approximately equal to M. (D) for all i ; and so the plot of X. versus M .(D) (referred to as probability plot) will be approximately linear. This linearity

will, in turn, result in a near unity value in r . Thus, the better the fit of the distribution D to the data, the closer r will be to unity [10,22].

The procedure just described makes use of 43 extreme value type T I distributions defined by various values of y from 1 to 25 in steps of 1, from 25 to 50 in steps of 5, from

50 to 100 in steps of 10, from 100 to 250 in steps of 50, y = 350, y = 500, y = 750, y = 1000 and y = °°. For any given data set, 43 probability plot correlation coefficients are computed corresponding to the&e distributions, and the distribution with the maximum probability plot correlation coefficient is chosen as the one which best fits the data. The

final result from this first stage is the value y of the y corresponding to the estimated

best fitting distributed.

The second stage in the procedure consists of estimating the location and scale parameters y and a, respectively, in Eqs. 1, 2, 3 and 4 for the observed data set and for the

determined optimal value y as determined in stage 1. Estimates of the location and scale

opt ° follow directly from the basic probability plot approach. If a least squares line is fit to
the probability plot corresponding to y , then the computed intercept and slope of the

fitted line serve as estimates for the unknown location and scale parameters y and a. In terms of the X. and M.(D), these estimated location and scale values u and O axe as follows:

1 1 E(X. - X)(M.(D) - M(D))

=

~

? (9)

Z[M.(D) - M(D)]

V = X. - M(D) (10)

The third and final stage in the procedure determines the predicted wind speed v— for

various intervals N of interest (say, N = 50, 100, 500, and 1000 years). The estimate for v— is given by

v^ = y + o Gx

(1 - 1/N) (11) Y

opt

where y is the optimal value of y (as determined in stage 1) , y and o are the estimates

of the location and scale parameters y and a in Eqs. 1, 2, 3, and 4 (as determined in stage

iY,

2), and G (p) is the percent point function of the best fitting extreme value distribute "opt

'opt

best fit) , then

»

tion. If Y 4 °° (that is, if a member of the extreme value type II family provides the

opt

l

Gx

(p) = (-Jin p)"1/Y (12a) Y

"opt

If y = oo (that is, if the extreme value type 1 distribution provides the best fit), then

Gx

(p) = -Jin (-An p) (12b) Y

'opt

In effect, the procedure described in this section is an automated equivalent of probability

paper plotting in which 45 types of probability paper, corresponding to 43 extreme value

distributions, would be used and in which fitting would be carried out on the basis of the

least squares method, rather than by eye.

It is noted that the procedure described is applicable without modification to the

extreme value analysis of any physical phenomenon, i.e., it is in no way restricted to the

analysis of extreme winds.

A listing of the computer program, and sample inputs and outputs, are given in an

Appendix.
4. SUMMARY AND INTERPRETATION OF RESULTS

The data analyzed were obtained from Ref . 5 and are listed in table 1. The main results

of the analysis are 3 isted in table 2. The quantity y «. is the value of y (see Eq. 2) for

which the best distributional fit of the largest values is obtained. The quantities vare extreme wind speeds corresponding to a N-year mean recurrence interval (N = 50, 100 and

1000) and were calculated using Eq. 2 with y

=

Y or, if y = °°, Eq. 1. The quantities

v, s , v are the mean, the standard deviation and the maximum value of the largest yearly

v max ' b j j

CO winds, respectively. The quantities in parentheses are the values vzz of the extreme winds

corresponding to a N-year recurrence interval calculated assuming y = °°. These quantities were omitted in those cases in which y = °° or y was sufficiently large for the

'opt 'opt

difference v- - v- to be insignificant. N N

Optimum Probabilistic Models

The results of the analysis may be conveniently divided into four categories, characterized by the following ranges of values y : (1) y > 40: (2) 10 < y < 39; (3) 5 < J a o i p t

> opt _ _ i opt; _ > _

Y _£ 9; (4) 2 <^ y .£ 4. Of the 21 37-year series of data considered, about 45% belong

to the first category, 25% to the second, 10% to the third, and 20% tn the fourth. Similar

percentages were obtained from the analysis of the 30-, 25- and 20-year sets of data listed

in table 2. It is seen that the assumption of a unique generally valid distribution, i.e., one characterized by a single value of y .. — whether it be a type I distribution, or a type II distribution with y = and a specified value of y as proposed in Ref. 17, 18 and

19 — is not confirmed by the results presented herein.

The results obtained also showed that a mixed probability distribution (Eq. 6) cannot

improve the empirical fit of the data in any significant way. In cddition, it is noted

that during a normal period of record (20-40 years) the frequency at any one station of winds

associated with hurricanes is small (of the order of one in 20 years or even less) and

therefore the sample size is insufficient for carrying out a meaningful statistical analysis.

That this is the case can be seen in table 1, which shows that in the period 1912-1948,

5-minute winds in excess of 70 mph were recorded only twice at both the Key West, Florida

and the Corpus Christi, Texas weather stations.

Length of Record and Reliability of Predictions

According to Shellard [16], for extreme wind predictions to be acceptable the length

of the record used should be of at least 15, or preferably 20 years. Thorn proposed isotach

maps for the United States using records of 15- to 21-year average length [17] . Implied in

Shellard' s statement and in Thorn's work is the assumption that, at any given station, the mean value, the standard deviation and the sample distribution of the largest annual winds

8
for a 20-year record are essentially the same as for any longer record at the station, i.e.,

that a 20-year record is a representative segment of a statistically stationary time series.

The extent to which this assumption is warranted was checked in each of the 21 cases

included in table 2. The 37-year records were broken up into 30-, 25- and 20-year overlapping records which were separately analyzed. The results of the analysis are included in

table 2 and are summarized in table 3, which shows the average values of the quantities

, max min. , min ,. ^, _ .. . . _ n - „ . „_ ,_

_, ^ _ (v- - vrr )/v— for the four ranges of y : v > 40; 10 < y < 39: 5 < y ._ < 9; N N N . opt opt — — opt — — opt —

, , max mm , . ,, .. — , ,. , , 2 < V < 4, where v— , v— are the maximum and the minimum N-year speed predicted on the — opt — N N basis of the 20-year records at one station and y is the tail length parameter of the

distribution estimated from the 37 years of data at that station.

_._-„_ . ... , , max mm. , min ,. TABLE 3 - Average Values of (v|r - v— )/vrr for

Various Ranges of Y and Various Mean

opt

Recurrence Intervals

N (years) 50 100 1000

Y > 40 .13 .17 .39

opt —

10 < y < 39 .23 .29 .55 — opt —

5 < y < 9 .46 .53 2.40 — opt —

2 < y < 4 .54 .96 7.20 — opt —

The lower bounds for the error in the estimation of the N-year winds may be calculated

in the case of a type I distribution on the basis of a mathematical statistical result,

viz., the Cramer-Rao relation, which states that for the EVI distribution [see Ref. 11, p. 282]

,\ . 1.10867a „.,. var(y) > (13) n

" . 0.60793a2 .... var(a) > (14) — n

where var(y), var(a) are the variances of the estimated values of \i, a', is the actual

value of the scale parameter and n is the sample size. Using Eqs. 13, 14, the quantity

3SD(V ), where SD(V ) is the standard deviation of the error in the estimation of the

100-year wind, was calculated for n = 20 years and n = 50 years and for typical values of q

obtained from the analysis of the data. The results of the calculations are shown in Fig. 1 and show that even for type I distributions, the estimation errors are not negligible. The

results presented in tables 2 and 3 suggest 66 that as Y decreases such errors become intol-

'opt

erably large and that extreme caution is thus in order in the interpretation and use in
o

If)

in

CM

Q.

E

-) O

2 o

e

o

o

If)

p

u

c

0)c

01ft

<u o

60

Pn

o o o

ro O

CVJ

qdui (00l A)aS£

10
structural design of probabilistically computed extreme winds. It is therefore the belief

of the writers that further research is needed into the question of the validity of current

probabilistic approaches to the definition of design wind speeds.

5. CONCLUSIONS

From the analysis of the sets of data reported herein, the following results were

obtained:

1. No single distribution was universally applicable to all the data sets. The type

I distribution was applicable in about 45% of the cases. In about 25% of the cases, the tail length parameter was 10 <_ y £ 39, in about 10% of the cases 5

< Y < 9 and in about 20% of the cases, 2 < Y < 4. — opt — — opt —

\ 2. No necessary correlation was noted between type of wind climate and the magnitude

of the tail length parameter, i.e., both type I distributions and type II distributions with small tail length parameters were found to fit series of data generated by tropical storms, as well as data generated by extratropical storms.

3. Predictions of extreme wind speeds based on records of 20-year length were found

to vary, on the average, by about 15%-100% for 50-years or 100-year recurrence

intervals and 40% to a few hundred percent for a 1000-year recurrence interval

between different 20-year sets of data taken at the same station.

These results suggest that while . in principle, probabilistic methods provide the most

rational approach to specifying design wind speeds, it is of the utmost importance that the errors inherent in this approach be carefully taken into account. In particular, predictions of wind speeds corresponding to 1000-year recurrence intervals appear in most cases

to be far too unreliable to be used with any reasonable degree of confidence for purposes

of structural design. Therefore, unless Very carefully substantiated .the use for such

purposes of 1000-year winds in conjunction with reduced values of the safety factor, which

has been suggested recently, would, in the writers' opinion, be unwise. In the light of the

results obtained, it also appears that the reliability of predictions based on records of

20-year length or so may in certain cases be quite unsatisfactory. This question is therefore believed to merit further investigation and is currently under study at the National

Bureau of Standards.

11
6. REFERENCES

1. "American National Standard Building Requirements for Minimum Design Loads in Buildings

and Other Structures A58.1", American National Standards Institute, New York, N.Y.,

1972.

2. Benjamin, J.R. , and Cornell, C.A., "Probability, Statistics, and Decision for Civil

Engineers" , McGraw-Hill Book Company, New York, 1970.

3. Canadian Structural Design Manual , Supplement No. 4 to the National Building Code of

Canada, National Research Council of Canada, 1970.

4. Caspar, W. , "Maximale Windgeschwindigkeinten in der Bundesrepublik Deutchland", Bautechnik, Vol. 47, No. 10, October 1970, pp. 335-340.

5. Court, A., "Wind Extremes as Design Factors", Journal of the Franklin Institute , Vol.

256, July 1953, pp. 39-55.

6. Davenport, A.G. , "The Dependence of Wind Loads upon Meteorological Parameters", Proceedings , Wind Effects on Buildings and Structures, Ottawa, Vol. 1, University of

Toronto Press, Toronto, Canada, 1978.

7. Duchene-Marullaz , P., "Etude des vitesses maximales du vent", Cahiers du Centre

Scientifique et Technique du Batiment , Cahier 1118, No. 131, July-Aug, 1972.

8. Durs,t, C.S., "Wind Speeds over Short Periods of Time", Meteorological Magazine , vol.

89, 1960.

9. Filliben, J.J., "Simple and Robust Linear Estimation of the Location Parameter of a Symmetric Distribution", Ph.D. Dissertation, Princeton University, Princeton, N. J.

1969.

10. Filliben, J.J. , "Techniques for Tail Length Analysis", Proceedings , 18th Conference on

the Design of Experiments in Army Research, Development and Testing, Aberdeen, Md.,

October 1972.

11. Johnson, N.L. , and Kotz, S., Continuous Univariate Distributions - 1 , John Wiley and

Sons, New York, 1970.

12. Lieblein, J., "Method of Analyzing Extreme Value Data", Technical Note 3053 , National

Advisory Committee for Aeronautics, Washington, D.C.

13. Local Climatological Data Annual and Monthly Summaries , U.S. Department of Commerce,

National Environmental Data Service, National Climatic Center, Ashville, N.C.

14. Sachs, P., Wind Forces in Engineering , Pergamon Press, Oxford, 1972.

15. Shellard, H.C., "Extreme Wind Speeds over Great Britain and Northern Ireland", Meteorological Magazine, vol. 37, 1958.

16. Shellard, H.C., "The Estimation of Design Wind Speeds", Proceedings , Wind Effects on Buildings and Structures, Vol. 1, Her Majesty's Stationery Office, London, England,

1965.

17. Thorn, H.C.S., "New Distributions of Extreme Winds in the United States", Journal of

The Structural Division , ASCE, Vol. 94, No. ST7, Proc. Paper 6038, July 1968, pp. 1787-

1801.

18. Thorn, H.C.S., "Engineering Climatology if Design Winds, with Special Reference to the

Chicago Area", Proceedings , Symposium on Wind Effects on High-Rise Buildings, Northwestern University, Evanston, IL., March 1970.

12
19. Thom, H.C.S., "Distributions of Extreme Winds Over Oceans", Journal of the Waterways,

Harbors and Coastal Engineering Division , ASCE, Vol. 99, No. WW1, Proc. Paper 9556,

Feb., 1973, pp. 1-17.

20. U.S. Weather Bureau, Manual of Surface Observations, Washington, D.C. , 1951, p. 92.

21. Vellozzi, J., and Cohen, E. , "Dynamic Response of Tall Flexible Structures to Wind

Loading", Proceedings , Technical Meeting Concerning Wind Loads on Buildings and Structures, Building Science Series 30, National Bureau of Standards, SD Catalog No. C 13.29/2:30, U.S. Government Printing Office, Washington, D.C, 1970.

22. Filliben, J.J., "The Probability Plot Correlation Coefficient Test for Normality",

Technometrics, Vol. 17, No. 1, Feb., 1975, pp. 111-117.

13
0) H Dl

2 H- Z

rr rr Q.3 rr 3* -

Ml CA fD

ON UJ On roUiUJUiUJUiUJUi

ON UJ Ul UJ UJ

00 O h- Ln v£>

.p- UJ On UJ uj O 00 CO 00 Ui O H UJ 00

On UJ ON

On UJ CO uj u> on UJ ON UJ UJ

ro on o i—CO 00 M

Ul M \D CO h- -C- CO CO o o

\0 CO

CO ^ H N) \D

00 On ON 00

CO UJ ro NJ vO ^O

ro uj uj uj

0o ui uj ^O \£> NJ 00

O CO 00 -^1 -*J **j Ul ON U> ON UJ ro

UJ ro UJ UJ UJ

00 *-J Ul ON \D ON UJ ON I—

^O UJ 00

J>"»JUJU-1UJUJUJUJ

v© UJ N3 ON

00 CO ro U*i Ui ro |—

ui uj uj uj ON ON W H*

ON UJ UI U> UJ UJ UJ

U*i -^J ^ UJ --J 0O O

Ul UJ UJ Ul

ro uj uj ui uj ui uj onujuiujioujujuj

M -O I—

--J Ul 00 M -J

ro ro oo oo

oo on on Ui

ui w ro w ONUJuiuJrorouJUJUJUJ

uj ui uj ro uj uj

vO ON CO UJ

ONUJUJUJUJUiUiUiUJUJ^- O'^JUiroOh-'ONI-'UivOON

-^jro^jo-ouJON-cuj-p-onujonujujuj >»oONcouirororooo

on oo ro \o

ro uj ro ro -v UlUJONUJOOUJUJUJ

UJ ui UJ UJ

ro uj UJ uj

so |_» ^j

vO UJ oo

£» 4> UJ UJ Ui UJ ON UJ 00 -J kD 00

ON <*D UJ I—

UJ UJ UJ On

00 UJ \D h-1

ON UJ ON Ui Ui

h-* ui ui ro -vj

rouJUJONUJUitouj CO CO UJ

c* *o. ui to

on ui ro vo

ui ro u> uj uj uj on

UJ Ui UJ uj ** **

00 Ui I—' CO on I—

(-•oooNi—UJroOvororo

oouirouJUJONUJuiuJUJUi uiONUiujrorouix^^j

ui-t-ON-c-ONUiujroujui

on «X> -o

Oo u> -j on ^o UJ ui ui ro ujH 00 U

UJ ON ON UJ UJ

Ul Ui UJ UJ UJ UJ

^J ON UJ UJ Ui 4>

UJ UJ Ul

CO ^O P

*- I— "-J

ro ro on vo M M H- O

14
TABLE 2 - Results of Analysis

Station

Cairo

(111.)

Alpena

(Mich.)

Tatoosh

Island

(Wash.)

Williston

(N.D.)

Richmond

(Virginia)

No. of

Years

37

30

25

20

37

30

25

20

37

30

25

20

3730

25

20

3730

25

20

Years of

Record

1912-48

1912-41

1919-48

1912-36

1924-48

1912-31

1917-36

1923-42

1929-48

1912-48

1912-41

1919-48

1912-36

1924-48

1912-31

1917-36

1923-42

1929-48

1912-48

1912-41

1919-48

1912-36

1924-48

1912-31

1917-36

1923-42

1929-48

1912-48

1912-41

1919-48

1912-36

1924-48

1912-31

1917-36

1923-42

1929-48

1912-48

1912-41

1919-48

1912-36

1924-48

1912-31

1917-36

1923-42

1929-48

'opt 50 100 1000

(miles per hour)

max

oo 53 56 65 39.9 4.8 51

oo 53 56 65 39.9 4.7 51

oo 53 56 64 40.6 4.6 51

oo 52 54 62 39.8 4.3 47

oo 53 56 65 40.2 4.6 51

oo 53 55 64 40.4 4.5 47

oo 52 55 62 40.8 4.2 47

CO 54 56 66 40.1 4.8 51

oo 52 55 64 39.6 4.6 51

CO 51 53 59 41.9 3.4 50

CO 50

N 52

52 57 41.4 3.1 47

oo 54 61 42.2 3.6 50

oo 50 51 57 41.4 3.0 47

00 52 54 61 42.2 3.5 50

CO 50 51 57 41.5 2.9 47

oo 51 52 59 41.7 3.2 47

oo 50 52 58 41.5 3.1 47

oo 53 55 62 42.5 3.7 50

CO 83 86 99 64.8 6.6 84

CO 84 88 101 64.8 7.0 84

18 84 88 104 65.3 6.2 84

CO 83 87 99 63.9 6.9 84

oo 79 83 93 64.8 5.2 78

oo 84 89 103 64.0 7.4 84

12 86(84) 91(88) 111(101) 64.2 7.0 84

oo 81 85 96 65.5 5.7 78

oo 81 84 95 65.0 5.6 78

00 53 55 64 39.8 4.8 50

oo 53 56 65 39.9 4.8 50

CO 53 56 65 39.7 4.8 50

18 54 57 69 39.9 4.8 50

CO 52 55 64 39.2 4.7 50

12 55(54) 59(57) 73(66) 39.8 5.1 50

13 54(53) 58(56) 61(65) 39.8 4.7 50

24 53(52) 56(55) 67(64) 39.0 4.7 50

CO 51 54 62 39.3 4.3 46

CO 52 54 63 39.9 4.7 48

OO 52 54 63 38.9 4.7 48

CO 49 51 58 37.7 4.0 46

CO 52 55 63 39.8 4.6 48

OO 50 52 60 38.3 4.1 46

oo 53 56 65 40.0 4.7 48

oo 51 54 62 39.0 4.4 47

CO 49 51 59 38.1 3.9 46

CO 51 53 62 38.1 4.5 46

15
TABLE 2 - Results of Analysis (continued)

Station

No. of

Years

Years of

Record

opt 50 100 1000

(miles per hour)

37 1912-48 OO 53 55 62 47.3 3.5 53

30 1912-41 CO 54 56 63 44.1 3.7 53

1919-48 OO 53 55 62 44.7 3.5 53

Burlington 25 1912-36 OO 54 57 63 44.5 3.6 53

(Vermont) 1924-48 CO 50 52 57 42.8 2.8 50

20 1912-31 OO 56 58 65 45.1 3.8 53

1917-36 OO 55 57 64 45.1 3.6 53

1923-42 8 54(54) 58(56) 71(63) 43.5 3.7 53

1929-48 OO 48 49 53 42.2 2.2 47

37 1912-48 OO 62 65 74 48.5 4.9 60

30 1912-41 OO 62 65 77 48.5 5.1 60

1919-48 OO 60 62 71 47.7 4.5 57 Eastport 25 1912-36 OO 64 67 77 48.9 5.5 60

(Maine) 1924-48 OO 57 59 66 47.4 3.8 55

20 1912-31 OO 65 69 79 49.5 5.8 60

1917-36 CO 62 65 74 48.2 5.0 57 1923-42 OO 55 57 63 46.9 3.3 52

1929-48 OO 57 59 66 47.7 3.4 55

37 1912-48 OO 61 64 77 42.5 6.6 62

30 1912-41 CO 62 65 77 43.9 6.4 62

1919-48 5 63(59) 70(63) 101(75) 41.0 6.5 62

Canton 25 1912-36 CO 62 66 78 45.2 6.2 62

(N.Y.) 1924-48 2 66(56) 81(59) 186(70) 39.7 6.0 62

20 1912-31 CO 64 67 79 46.6 6.1 62

1917-36 7 65(62) 71(66) 95(78) 44.3 6.4 62

1923-42 2 70(58) 87(62) 203(73) 41.3 6.3 62

1929-48 13 49 52 62 38.1 3.7 48

37 1912-48 80 40 42 48 32.5 2.9 41

30 1912-41 OO 39 41 45 32.2 2.5 38

1919-48 60 41 43 48 32.5 3.1 41

Yuma 25 1912-36 OO 39 41 45 32.4 2.5 38

(Ariz.) 1924-48 18 42 44 52 32.4 3.3 41

20 1912-31 OO 39 41 46 32.3 2.6 38

1917-36 OO 40 42 47 32.7 2.7 38

1923-42 OO 40 41 46 32.1 2.8 38

1929-48 150 42 44 51 32.9 3.4 41

37 1912-48 35 65 68 77 51.4 5.0 68

30 1912-41 90 65 68 78 51.1 5.1 68

1919-48 200 66 69 79 51.9 5.3 68

Duluth 25 1912-36 OO 62 64 72 50.2 4.3 59

(Minn.) 1924-48 22 68 72 84 52.0 5.5 68

20 1912-31 OO 61 64 72 49.9 4.4 57 1917-36 OO 63 66 74 50.4 4.7 59

1923-42 OO 68 71 82 52.2 5.7 68

1929-48 7 71(68) 76(71) 98(82) 52.1 5.9 68

16
TABLE 2 - Results of Analysis (continued)

No. of Years of

Station Years Record Y v r „ v v V s v

'opt 50 100 1000 V ma: (miles per hour)

37 1912-48 30 54 56 65 41.1 4.6 56

30 1912-41 70 55 58 67 41.4 4.8 56

1919-48 23 55 59 70 41.2 4.9 56

Valentine 25 1912-36 CO 56 59 68 41.9 5.0 56

(Neb.) 1924-48 6 57(54) 62(57) 82(67) 40.6 4.9 56

20 1912-31 00 53 55 62 42.1 3.8 49

1917-36 70 57 61 71 42.0 5.4 56

1923-42 7 59(56) 64(59) 84(69) 41.1 5.3 56

1929-48 4 61(56) 67(59) 101(69) 40.1 5.4 56

37 1912-48 24 66 72 91 42.8 8.2 66

30 1912-41 10 66(64) 72(68) 96(83) 42.3 7.8 66

1919-48 9 69(66) 76(71) 104(88) 42.3 8.7 66

Charleston 25 1912-36 00 60 64 76 42.3 6.5 55

(S.C.) 1924-48 13 71(68) 77(74) 102(91) 43.4 9.0 66

20 1912-31 00 59 63 75 41.3 6.3 55

1917-36 11 60(58) 65(62) 84(74) 40.9 6.3 55

1923-42 4 73(66) 84(71) 137(87) 42.4 8.4 66 — 1929-48 12 73(71) 81(76) 108(95) 43.5 6.9 66

37 1912-48 23 46 48 55 35.6 3.7 46

30 1912-41 11 47(46) 50(48) 61(56) 35.3 4.0 46

1919-48 00 43 44 50 34.8 2.9 41

Eureka 25 1912-36 12 48(47) 51(50) 64(58) 35.3 4.3 46

(Calif.) 1924-48 00 42 44 49 34.6 2.9 41

20 1912-31 11 49(48) 53(51) 65(59) 30.1 4.3 46

1917-36 00 43 45 51 34.3 3.2 40

1923-42 00 41 43 48 34.1 2.7 38

1929-48 oo 43 44 50 34.6 2.9 41

37 1912-48 15 60(59) 65(63) 85(77) 37.7 7.7 57

30 1912-41 19 62(60) 67(65) 86(80) 39.1 7.8 57 1919-48 00 48 51 60 35.0 4.8 46

Oklahoma 25 1912-36 00 62 67 82 40.4 7.9 57 City 1924-48 00 45 47 55 33.8 4.1 46

(Oklahoma) 20 1912-31 oo 64 68 83 42.2 7.7 57 1917-36 6 60(56) 66(60) 94(73) 38.2 6.5 57

1923-42 9 47(45) 50(48) 62(55) 34.7 3.9 46

1929-48 00 42 44 50 32.8 3.3 38

37 1912-48 11 42(41) 45(43) 56(50) 30.1 4.0 40

i 30 1912-41 9 43(41) 46(44) 59(51) 30.2 4.2 40

1 1919-48 7 42(40) 46(42) 59(49) 29.8 3.7 40

i Baker 25 1912-36 8 43(42) 47(44) 61(52) 30.2 4.2 40

(Oregon) 1924-48 10 42(41) 46(43) 58(51) 30.0 4.0 40

( 20 1912-31 6 45(42) 49(45) 67(53) 30.3 4.5 40

I 1917-36 7 44(42) 48(44) 63(52) 30.2 4.2 40

1923-42 10 44(42) 47(45) 61(53) 30.3 4.3 40

1929-48 7 40(39) 43(41) 56(47) 29.5 3.3 38

17
TABLE 2 - Results of Analysis (continued)

No. of Years of

Station Years Record Y v rri V v V s V

'opt 50 100 1000 V ma: (miles per hour)

37 1912-48 7 69(65) 77(70) 311(87) 39.8 9.4 66

30 1912-41 4 54(50) 61(53) 92(63) 36.0 5.2 53

1919-48 7 71(67) 80(73) 117(91) 39.8 10.3 66

Sheridan 25 1912-36 12 49(48) 52(50) 64(58) 35.7 4.3 46

(Wyoming) 1924-48 8 74(69) 83(76) 119(96) 40.4 11.0 66

20 1912-31 200 49 52 61 36.3 4.6 46

1917-36 2 51(44) 59(46) 124(52) 34.3 3.6 46

1923-42 2 62(50) 96(53) 182(63) 35.3 5.8 53

1929-48 35 75 82 107 42.2 11.7 66

37 1912-48 5 83(79) 91(83) 126(97) 58.4 7.6 82

30 1912-41 3 81(66) 91(78) 145(89) 56.9 6.5 82

1919-48 4 86(80) 96(85) 144(99) 58.2 8.2 82

Block 25 1912-36 oo 68 71 78 56.8 4.2 66

Island 1924-48 5 88(82) 97(87) 138(103) 59.1 8.5 82

(R.I.) 20 1912-31 00 68 71 79 56.3 4.4 66

1917-36 oo 65 67 74 56.0 3.5 62

1923-42 2 90(76) 107(80) 238(93) 56.45 7.4 82

1929-48 6 91(86) 100(91) 140(109) 59.6 9.5 82

37 1912-48 4 56(53) 74(56) 97(66) 37.3 5.7 58

30 1912-41 2 58(50) 70(52) 156(62) 36.2 5.2 58

1919-48 4 59(54) 67(58) 103(69) 37.6 6.1 58

Winnemucca 25 1912-36 oo 45 47 53 35.7 3.4 43

(Nevada) 1924-48 3 63(55) 73(59) 131(71) 37.4 6.6 58

20 1912-31 oo 45 47 54 36.0 3.3 43

1917-36 00 46 48 55 35.9 3.6 43

1923-42 2 68(55) 85(58) 211(71) 36.6 6.9 58

1929-48 4 64(58) 73(62) 118(75) 38.1 7.1 58

37 1912-48 4 94(89) 103(93) 146(106) 69.3 7.3 95

30 1912-41 4 94(91) 107(95) 153(110) 69.7 7.8 95

1919-48 4 94(91) 107(96) 154(110) 69.7 7.9 95

North 25 1912-36 3 100(91) 113(95) 183(110) 69.2 8.0 95

Head 1924-48 4 92(87) 100(90) 137(102) 69.7 6.2 87

(Wash.) 20 1912-31 2 102(87) 120(91) 255(104) 67.7 7.5 95

1917-36 3 105(94) 119(99) 198(115) 69.6 8.8 95

1923-42 4 94(88) 102(92) 142(104) 70.1 6.5 87

1929-48 4 95(89) 103(93) 145(105) 70.3 6.7 87

37 1912-48 3 83(71) 99(77) 188(97) 43.4 10.7 84

30 1912-41 2 84(67) 106(72) 271(89) 42.4 9.9 84

1919-48 3 88(75) 106(82) 206(103) 44.3 11.7 84

Key West 25 1912-36 2 89(70) 114(75) 300(94) 43.2 10.4 84

(Florida) 1924-48 5 77(70) 87(76) 135(95) 43.4 9.7 73 20 1912-31 2 97(73) 125(79) 335(99) 44.1 11.3 84

1917-36 2 92(74) 126(80) 336(100) 44.3 11.3 84

1923-42 10 64(61) 69(66) 92(79) 41.4 7.1 61

1929-48 3 82(73) 96(79) 162(99) 43.2 10.6 73

18
TABLE 2 - Results of Analysis (continued)

Station

No. of

Years

Years of

Record

'opt 50 100 1000

(miles per hour)

max

37 1912-48 2 97(78) 125(85) 327(107) 46.3 12.7 95

30 1912-41 2 101(78) 131(85) 351(107) 45.5 13.6 95

1919-48 2 93(74) 119(80) 310(100) 45.6 11.4 95

Corpus 25 1912-36 2 109(82) 142(89) 390(114) 46.2 14.7 95

Christi 1924-48 70 63 67 82 44.0 6.9 61

(Texas) 20 1912-31 2 119(86) 159(94) 448(122) 46.3 16.4 95

1917-36 1 132(76) 228(82) 1952(103) 44.9 12.9 95

1923-42 5 65(61) 72(64) 104(77) 42.5 6.5 61

1929-48 00 65 68 83 45.2 7.1 61

19
APPENDIX

COMPUTER PROGRAM LISTING

SAMPLE INPUT AND OUTPUT

20
JJF6*SIMIU.MAIN

1 C

2 C 3 C 4 C

5 C 6 c 7 c 8 c 9 c 10 c

11 c 12 c 13 c 14 c 15 c 16 c 17 c 18 c 19 c 20 c 21 c 22 c 23 c 24 c 25 c 26 c 27 c 28 c 29 c 30 c 31 c 32 c 33 c 34 c 35 c 36 c 37 c 38 c 39 c 40 c 41 c 42 c 43 c 44 c 45 c 46 c 47 c 48 c 49 c 50 c 51 c 52 c 53 c 54 c 55 c 56 c 57 c

PURPOSE—THIS MAIN PROGRAM READS IN DATA UPON WHICH AN EXTREME VALUE ANALYSIS IS TO BE PERFORMED.

AFTER READING IN THE DATA. THIS PROGRAM

CALLS IN A SUBROUTINE (EXTREM) WHICH

PERFORMS THE EXTREME VALUE ANALYSIS.

INPUT DATA—THE NUMBER OF SETS OF DATA TO BE ANALYZED

(FORMAT— 12) THE IDENTIFYING TITLE FOR DATA SET 1 (FORMAT—80A1) THE NUMBER OF OBSERVATIONS IN DATA SET 1 (FORMAT—12) THE DATA FOR SET 1 (EACH DATA CARD HAVING A 16F5.1 FORMAT)

THE IDENTIFYING TITLE FOR DATA SET 2 (FORMAT—80A1) THE NUMBER OF OBSERVATIONS IN DATA SET 2 (FORMAT—12) THE DATA FOR SET 2 (EACH DATA CARD HAVING A 16F5.1 FORMAT)

THE IDENTIFYING TITLE FOR THE LAST DATA SET

(FORMAT—80A1) THE NUMBER OF OBSERVATIONS IN THE LAST DATA SET

(FORMAT— 12) THE DATA FOR THE LAST SET

(EACH DATA CARD HAVING A 16F5.1 FORMAT)

OUTPUT—THIS MAIN PROGRAM WILL (FOR EACH DATA SET)

SKIP TO A NEW PAGE* PRINT OUT THE TITLE* PRINT OUT THE NUMBER OF OBSERVATIONS*

AND PRINT OUT THE INPUT DATA. THIS WILL THEN BE FOLLOWED (FOR EACH DATA SET) BY 4 OR 5 (DEPENDING ON THE DATA)

PAGES OF AUTOMATIC OUTPUT RESULTING

FROM THE EXTREME' VALUE ANALYSIS

SUBROUTINE EXTREM WHICH IS CALLED BY THIS MAIN PROGRAM.

SUBROUTINES NEEDED—EXTREM* SORT* UNIMED* EV1PLT* EV2PLT* AND PLOT

LANGUAGE—ANSI FORTRAN

COMMENT—THIS MAIN ROUTINE AND ALL SUBROUTINES

ASSUME THAT THE INPUT AND OUTPUT UNITS

HAVE A NUMERICAL DESIGNATION OF 5 AND 6* RESPECTIVELY.

THIS DESIGNATION IS MADE WITH THE FORTRAN

STATEMENTS

IRD=5

IPR=6

ONE OR BOTH OF WHICH ARE FOUND AT THE BEGINNING OF THE EXECUTABLE CODE

IN THIS MAIN ROUTINE AND ALL SUBROUTINES.

IF 5 AND 6 ARE NOT THE PROPER DESIGNATIONS

FOR YOUR COMPUTER* THEN SIMPLY CHANGE THE

21
58 C 5 AND 6 IN THE IRD=5 AND IPR=6 STATEMENTS

59 C TO THE APPROPRIATE VALUE FOR YOUR COMPUTER.

60 C COMMENT—THIS MAIN ROUTINE AND ALL SUBROUTINES

61 C WILL* AS THEY STAND* ACCEPT DATA SETS

62 C WITH UP TO 200 OBSERVATIONS.

63 C IE COMPUTER STORAGE IS LIMITED AND IF SMALLER DATA

64 C SETS ARE EXPECTED* THEN COMPUTER STORAGE MAY BE SAVED

65 C BY RESETTING THE DIMENSION LIMITS OF THE VECTORS

66 C X IN THIS MAIN ROUTINE*

67 C Wi Y» AND Z IN THE SUBROUTINE EXTREM,

68 C W AND Y IN THE SUBROUTINE EVIPLTr AND

69 C W AND Y IN THE SUBROUTINE EV2PLT

70 C FROM THEIR PRESENT VALUE OF 200

71 C TO WHATEVER THE EXPECTED MAXIMUM DATA SET SIZE IS. 72 C THE STATEMENT IUPPER=200

73 C (ABOUT STATEMENT 122 IN THIS SUBROUTINE*

74 C AND UP NEAP THE BEGINNING OF THE EXECUTABLE

75 C CODE IN THE SUBROUTINES EXTREM, EV1PLT. AND

76 C EV2PLT) SHOULD ALSO BE CHANGED FROM 200 77 C TO THE EXPECTED MAXIMUM DATA SET SIZE.

78 C COMMENT—ON THE UNIVAC 1108* EXEC 8 COMPUTER SYSTEM -AT NBS*

79 C THIS MAIN ROUTINE AND THE 6 NEEDED SUBROUTINES

80 C HAVE A TOTAL (CODE + DIMENSIONS + COMMON) STORAGE

81 C REQUIREMENT OF APPROXIMATELY 13000 WORDS (DECIMAL).

82 C 83 C 84 C 85 C 86 C 87 C 88 C 89 C 90 C NOTE THE RELATIVELY LARGE STORAGE ALLOCATION

91 C FOR THE LAST SUBROUTINE (PLOT). IF THE AMOUNT

92 C OF USABLE STORAGE IN YOUR COMPUTER IS LESS THAN

93 C 13000* THEN AN ALTERNATIVE PLOT ROUTINE IS 94 C AVAILABLE FROM THE AUTHOR WHICH WILL REDUCE

95 C THE TOTAL STORAGE ALLOCATION FROM 13000

96 C TO APPROXIMATELY 6000.

97 C WRITTEN BY—JAMES J. FlLLIBEN (205.03)

98 C EMIL SIMIU (461.01)

99 C NATIONAL BUREAU OF STANDARDS

100 C WASHINGTON* C. C. 20234

101 C UPDATED—DECEMBER 1974

102 C 103 DIMENSION X (200 ) * ITITLE(80

104 C 105 IRD=5

106 IPR=6

107 C 108 C READ IN THE NUMBER OF SETS OF DATA TO RE ANALYZED

109 C 110 REAO(IRD*105)NUMSET

111 C 112 C OPERATE ON EACH SET

113 C 114 DO100ISET=l»NUMSET

115 C

CODE DIMENSIONS COMMON

MAIN PROGRAM 130 320 EXTREM 770 1680

SORT 340 180

UNIMED 140 80

EV1PLT 210 550

EV2PLT 280 560

PLOT 380 190 7150

22
116 C READ IN THE TITLE AND THE NUMBER OF OBSERVATIONS

117 C FOR THIS SET

118 C 119 READdRD»205) ( ITITLE< I >

I=l»80)

120 READ(IRD»210)N

121 C 122 C ZERO-OUT THE X VECTOR. AND THEN READ THE DATA FOR THIS SET

123 C INTO THE X VECTOR

124 C 125 DO200I=l.N

126 X(I)=0.0

127 200 CONTINUE

128 READ(IRD»215)(X(I)#I=1.N)

129 C 130 C WRITE OUT THE TITLE* THE NUMBER OF OBSERVATIONS.

131 C AND THE DATA FOR THIS SET

132 C 133 WRITEdPR.998)

134 NSKIPrlO

135 D0300ISKIP=1#NSKIP

136 WRJTEdPR.999)

137 300 CONTINUE

138 WRITEdPR»305) dTITLEd >» 1=1.80)

139 WRITE(IPR»999)

140 WRITECIPR.999)

141 WRITE(IPR.999)

142 WRITE(IPR.310)N

143 WRITE(IPR»999)

144 WRITE(IPR.315)

145 WRITE(IPR»320) (X(I).I=1.N)

146 C 147 C DO AN EXTREME VALUE ANALYSIS OF THE DATA FOR THIS SET

148 C 149 CALL EXTREM(X.N)

150 C 151 100 CONTINUE

152 C 153 105 FORMAT (12)

154 205 FORMAT(80A1)

155 210 F0RMAT(I2)

156 215 FORMAT (16F5.1)

157 305 FORMATdH .20X.80A1)

158 310 FORMATdH .29HTHE NUMBER OF OBSERVATIONS = .15)

159 315 FORMATdH .10HINPUT DATA)

160 320 FORMATdH .13X.16F5.D

161 998 FORMAT (1H1)

162 999 FORMATdH > 163 STOP

164 END

I3PRT.S SIMIU.EXTREM

23
JJF6*SIMIU.EXTREM

1 SUBROUTINE EXTREM(X*N)

2 C

3 C

4 C

5 C 6 C 7 C 8 C 9 C 10 C 11 C 12 C 13 C 14 C 15 C 16 c 17 c 18 c 19 c 20 c 21 c 22 c 23 c 24 c 25 c 26 c 27 c 28 c 29 c 30 c 31 c 32 c 33 c 31 c 35 c 36 c 37 c 38 c 39 c 40 c 41 c 42 c 13 c 44 c 15 c

<+6 c 47 c 18 c 49 c 50 c 51 c 52 c 53 c 54 c 55 c 56 c 57 c

PURPOSE—THIS SUBROUTINE PERFOMS AN EXTREME VALUE ANALYSIS

ON THE DATA IN THE INPUT VECTOR X. THIS ANALYSIS CONSISTS OF DETERMINING THAT PARTICULAR

EXTREME VALUE TYPE 1 OR EXTREME VALUE TYPE 2 DISTRIBUTION

WHICH BEST FITS THE DATA SET. THE GOODNESS OF FIT CRITERION IS THE MAXIMUM PROBABILITY

PLOT CORRELATION COEFFICIENT CRITERION.

AFTER THE BEST-FIT DISTRIBUTION IS DETERMINED*

ESTIMATES ARE COMPUTED AMD PRINTED OUT FOR THE LOCATION AND SCALE PARAMETERS.

TWO PROBABILITY PLOTS ARE ALSO PRINTED OUT—

THE BEST-FIT TYPE 2 PROBABILITY PLOT

(IF THE BEST FIT WAS IN FACT A TYPE 2)» AND THE TYPE 1 PROBABILITY PLOT.

PREDICTED EXTREMES FOR VARIOUS RETURN PERIODS ARE ALSO COMPUTED AND PRINTED OUT. INPUT ARGUMENTS—X = THE SINGLE PRECISION VECTOR OF

(UNSORTED OR SORTED) OBSERVATIONS.

N =. THE INTEGER NUMBER OF OBSERVATIONS

IN THE VECTOR X. OUTPUT—6 PAGES OF AUTOMATIC PRINTOUT

PRINTING—YES

RESTRICTION—THE MAXIMUM ALLOWABLE VALUE OF N AS INPUT TO THIS SUBROUTINE IS 7500

OTHER DATAPAC SUBROUTINES NEEDED—SORT. UNIMED»EV1PLT»EV2PLT»PL0T

FORTRAN LIBRARY SUBROUTINES NEEDED—SORT AND ALOG

MODE OF INTERNAL OPERATIONS—SINGLE PRECISION

LANGUAGE—MACHINE-INDEPENDENT ANSI FORTRAN

COMMENT—THIS SUBROUTINE AS IT STANDS WILL ACCEPT DATA SETS

WITH UP TO 200 OBSERVATIONS.

IF COMPUTER STORAGE IS LIMITED AND IF SMALLER DATA

SETS ARE EXPECTED* THEN COMPUTER STORAGE MAY BE SAVED

BY RESETTING THE DIMENSION LIMITS OF THE VECTORS

W. Y» AND Z BELOW FROM THEIR PRESENT VALUE OF 200

TO WHATEVER THE EXPECTED MAXIMUM DATA SET SIZE IS. THE STATEMENT IUPPER=200

(ABOUT STATEMENT 122 IN THIS SUBROUTINE)

SHOULD ALSO BE CHANGED FROM 200 TO THE EXPECTED

MAXIMUM DATA SET SIZE.

' STORAGE SAVINGS CAN ALSO BE ACHIEVED BY JUDICIOUSLY

REDUCING THE DIMENSION LIMITS OF THE VECTORS W AND Y IN THE NEEDED SUPPORT SUBROUTINES EV1PLT

AND EV2PLT.

REFERENCE—FILLIBEN (1972)* 'TECHNIQUES FOR TAIL LENGTH

ANALYSIS'* PROCEEDINGS OF THE EIGHTEENTH

CONFERENCE ON THE DESIGN OF EXPERIMENTS IN ARMY RESEARCH AND TESTING* PAGES 125-150.

—FILLIBEN* 'THE PERCENT POINT FUNCTION'

UNPUBLISHED MANUSCRIPT.

—JOHNSON AND KOTZ (1970)* CONTINUOUS UNIVARIATE

DISTRIBUTIONS-1* PAGES 272-295.

WRITTEN BY—JAMES J FILLIBEN (JUNE 1972)

STATISTICAL ENGINEERING LABORATORY (205.03)

NATIONAL BUREAU OF STANDARDS

WASHINGTON* D. C. 20231

24
58 C PHONE—-301-921-2315

59 C UPDATED—DECEMBER* 1974

60 C 61 INTEGER BLANK » ALPHAM. ALPHAA . ALPHAY

62 INTEGER ALPHAI » ALPHAN. ALPHAF. ALPHAT. ALPHAY

63 INTEGER ALPHAG.EQUAL

61 DIMENSION X(l)

65 DIMENSION W (200 ) » Y ( 200

)

tl

200

)

66 DIMENSION GAMTAB < 50 ) .CORR (50

)

67 DIMENSION YK5C) »YS(50) »T(50)

68 DIMENSION IFLAG1 ( 50 ) » IFLAG2 ( 50 ) » IFLAG3 ( 50

)

69 DIMENSION C ( 10 ) . AM (50 ) .SCRAT (50

)

70 DIMENSION AINDEX(50)

71 DIMENSION P0(10)

72 DIMENSION H(6n»2)

73 DATA BLANK.ALPHAM.ALPHAA.ALPHAX/1H » lHMi IMA r 1HX/

74 DATA ALPHAI.ALPHAN.ALPHAF.ALPHAT.ALPHAY/1HI.1HN.1HF.1HT.1HY/

75 DATA ALPHAG»EQUAL/1HG.1H=/

76 DATA GAMTAB(l) »GAMTA8(2)»GAMTAP(3) »GAMTAB(4) .GAMTAB(5).

77 1GAMTAB(6) .GAMTAB(7) »GAMTAB(8) »GAMTAB(9) .GAMTAR(IO)

78 lGAMTAB(ll) .GAMTAB ( 12 ) . GAMTAR ( 1 3) .GAMTAR(14) .GAMTAB (15).

79 1GAMTAB(16) » GAMTAB ( 17) .GAMTAB ( 18) .GAMTAB (19) .GAMTAB (20).

80 1GAMTAB(2D >GAMTAB(22) »GAMTAB(?3) .GAMTAR(24) .GAMTAB(25)

81 1/1.»2.»3.»4.»5.»6.»7..8.»9..10.»11.»12..

82 113.. 14. .15. .16. .17. .18.. 19. .20. .21.. 22. .23.. 24. .25.'

83 DATA GAMTAB(26) .GAMTAB<27) .GAMTAB (28) .GAMTAB(29) »GAMTAB(30)

8*+ 1GAMTAB(3D »GAMTAB(32) »GAMTAB<33) .GAMTAB(34) »GAMTAB(35) . 85 1GAMTAB(36) »GAMTAB(37) »GAMTAB(3fl) .GAMTAB<39) .GAMTAB(40 )

86 1GAMTABU1) .GAMTAB(42)

87 1/30. . 35. . 40. . 45.. 50.. 60. » 70. » 80. »90. » 100. » 150. .200. .250..

88 1350. .500. .750. .1000./

89 DATA C(1)»C(2) »C (3) »C (4) »C (5) ,C (6) »C (7) iC (8) »C (9) »C ( 10)

90 1/60. .75. .100. . 150. » 250. » 500. .1 000. » 10000. » 100000. rlOOOOOO./

91 DATA P0(D »P0(2) »P0(3) .P0(4) »P0(5) .P0(6) »P0(7) »P0(8) »P0 (9) »P0 ( 10)

92 1/ . . . 5 . . 75 » . 9 . . 95 » . 975 . . 99 . . 999 . . 9999 . . 99999/

93 DATA AINDEX(l) .AINDEX(2) » AINDEX (3) » AINDEX (4 ) »AINPEX(5)»

94 1AINDEX(6) .AINDEX(7) »AINDEX(8) » AINPEX(9) .AINDEX(IO)

95 lAINDEX(ll) .AINDEXQ2) .AINDEXQ3) .AINDEX(14) .AINDEXQ5)

96 UINDEXU6) »AINDEX(17) .AINDEX(18) .AINDEX(IQ) rAIMDEX(20) . 97 1A INDEX (2D » A INDEX (22) » A INDEX (23) . A INDEX (24) .AINDEX (25)

98 1/1.»2.»3.»4.»5.»6.»7.»8..9..10..11..12.»

99 113. .14. .15. .16. .17. .18. .19. .20.. 21. .22.. 23. .24. .25./

100 DATA AINDEX(26) » AINDEX (27) » AINOEX (28) . AINDEX ( 29) . AINDEX (30)

101 1AINDEXOD »AINDEX(32> »AINPEX(33) .AINDEX(34) .AINDEX(35)

102 1AINDEX(36) »AINDEX(37) »AINOEX(38) »AINDEX(39) »AINDEX(40)

103 1AINDEX(41) .AINDEX(42)#AINnEX(43) .AINDEX (44) . AINDEX (45) . 104 1AINDEX(46) .AINDEX(47) »AINDEX(48) .AINDEX(49) »AINDEX(50)

105 l/26..27.»28.»29.»30..31.r32..33.r34.»35.»36.»37..38..

106 139. .40. »41. .42. .43. .44. .45. .46. .47. .48. .49. . 50 ./

107 DATA T(l).T(2)»T(3)»T<4)»T(5)»T(6).T(7)rT(8)»T(9)»T<10)»

108 1T(11).T(12)»T(13)»T(14).T(15).T(16)»T(17).T(18)»T(19).T(20)»

109 1T(21).T(22).T(23).T(24)»T(25)

110 1/1 0. 1801 If 3. 39672 . 2. 47043* 2. 14609. 1.98712.1. 89429 ,1. 83394

111 11. 79175 » 1.76069. 1.73691. 1.71814 .1.70297.1. 69045 rl. 67996

112 11. 67103 . 1. 66335.1. 65667 f 1.65002.1.64564 .1.64102 .1.63689.

113 11.63316.1.62979.1.62672.1.62391/

114 DATA T(26> »T(27) »T(28) »T(29) . TOO).

115 IT (31) .T(32) .T(33) .T(34) .T(35) .T(36) »T(37) .T(3B) . T( 39) »T(40)

25
116 1T(41)»T(42)»T(43)

117 1/1. 61287 . 1.60516.1. 59947 »1. 59510 .1.59164.1. 58651. 1.58289.

118 11.58019.1.57811.1.57645.1.57152.1.56908.1.56763.) .56666.

119 11 . 56546 . 1 . 56377 . 1 . 56.330 . 1 . 561 R7/

120 C 121 IPR=6

122 IUPPER=200

123 NUMDIS=43

124 AN=N

125 C 126 C CHECK THE INPUT ARGUMENTS FOR ERRORS

127 C 128 IF(N.LT.1.0R.N.GT.IUPPER)GOT050

129 IF(N.EQ.DG0T055

130 HOLD=X(l)

131 D060I=2»N

132 IF(X(I).NE.HOLD)GOTO90

133 60 CONTINUE

134 WRITEdPR. 9)H0LD

135 RETURN

136 50 WRITEdPR. 17) IUPPER

137 WRITE(IPR.47)N

138 RETURN

139 55 WRITEdPR. 18)

140 RETURN

141 90 CONTINUE

142 9 FORMAT (1H »109H***** NON-FATAL DIAGNOSTIC—THE FIRST INPUT ARGUME

143 INT (A VECTOR) TO THE EXTREM SUBROUTINE HAS ALL ELEMENTS = .E15.8.6

144 1H *****)

145 17 FORMAT (1H . 98H***** FATAL ERROR—THE SECOND INPUT ARGUMENT TO THE

146 1 EXTREM SUBROUTINE IS OUTSIDE THE ALLOWABLE (1..I6.16H) INTERVAL * 147 1****)

148 18 FORMAT (1H .100H***** NON-FATAL DIAGNOSTIC—THE SECOND INPUT ARGUME

149 INT TO THE EXTREM SUBROUTINE HAS THE VALUE 1 *****)

150 47 FORMATdH 35H***** THE VALUE OF THE ARGUMENT IS .18 .6H *****)

151 C 152 C COMPUTE THE SAMPLE MINIMUM AND SAMPLE MAXIMUM

153 C 154 XMIN=X(1)

155 XMAX=X(1)

156 DO140I=2.N

157 IF(X(I).LT.XMIN)XMIN=X(I>

158 IF(X(I).GT.XMAX)XMAX=X(D

159 140 CONTINUE

160 C 161 C COMPUTE THE PROB PLOT CORRELATION COEFFICIENTS FOR THE VARIOUS VALUES

162 C OF GAMMA

163 C 164 CALL SORT(X.N.Y)

165 CALL UNIMED(N.Z)

166 C 167 DO100IDIS=l»NUMDIS

168 IF(IDIS.EQ.NUMDIS)GOT0150

169 A=GAMTAB(IDIS)

170 DO110I=l»N

171 W(I) = (-ALOG(Zd)))**(-1.0/A)

172 110 CONTINUE

173 GOT0170

26
174 150 D0160I=1»W

175 W(I)=-ALOG(ALOG(1.0/Z(I)))

176 160 CONTINUE

177 C 178 170 SUM1=0.0

179 SUM2=0.0

180 D0200I=1»N

181 SUM1=SUM1+Y(I)

182 SUM2=SUM2+W(I)

183 200 CONTINUE

184 YBAR=SUM1/AN

185 WBAR=SUM2/AN

186 SUM1=0.0

187 SUM2=0.0

188 SUM3=0.0

189 DO300I=l»N

190 SUM2=SUM2+ ( Y ( I

) -YBAR ) * ( W ( T

) -WR AR

)

191 SUM1=SUM1+(Y(I)-YBAR)*(Y(1)-YRAR)

192 Sl'M3=SUM3+(W(I)-WRAR)*(W(T)-WRAR)

193 300 CONTINUE

194 SY=SQRT(SUM1/(AN-1.0)

)

^ 195 CC=SUM2/S0RT(SUM3*SUM1)

196 YSLOPF=SUM2/SUM3

197 YINT=YBAR-YSLOPE*WBAR

198 CORR(IDIS)=CC

199 YI(IOIS)=YINT

200 YS(lOTS)=YSLOPE

201 100 CONTINUE

202 C 203 C DETERMINE THAT DISTRIBUTION WITH THE MAX PROB PLOT CORR COEFFICIENT

204 C 205 IOISMX=l

206 C0RRMX=C0RR<1)

207 D0400IDIS=1»NUMDIS

208 IF(CORRdDIS) .GT.CORRMX ) IOISMX=IDIS

209 IF(CORRCIDIS) .GT.CORRMX)CORRMX=CORR ( IDTS)

210 400 CONTINUE

211 D0500IDIS=1>NUMDIS

212 IFLAGKIDIS)=BLANK

213 IFLAG2(IDIS)=BLANK

214 IFLAG3(IDIS)=BLANK

215 IF(IDI5.EO.IDISMX)G.OT0550

216 GOT0500

217 550 IFLAG1(IDIS)=ALPHAM

218 IFLAG2(IDIS)=ALPHAA

219 IFLAG3(IDIS)=ALPHAX

220 500 CONTINUE

221 C 222 C WRITE OUT THE TABLE OF PROB PLOT CORR COEFFICIENTS FOR VARIOUS GAMMA

223 C 224 WRITE(IPR»998)

225 WRITE(IPR»305)

226 WRITE(IPR»999)

227 WRITE(IPR»310)N

228 WRITE(IPR»311)YBAR

229 WRITE(IPR»312)SY

230 WRITE(IPR»313)XMIN

231 WRITE(IPR»314)XMAX

27
232 WRIT£(IPR»999)

233 WRITE(IPR»323)

234 WRITE(IPR»324>

235 WRITE(IPR»325)

236 WRITE(IPR»999)

237 C 238 NUMDM1=NUMDIS-1

239 IF(NUMDM1.LT.1)GOT0850

240 D0800I=1»NUMDM1

241 WR ITE ( IPR. 805 ) GAMTAR ( I

, CORR ( I ) . IFLAG1 ( I» IFLAG2 ( I ) » IFLAG3 < I ) » 242 lYKI).YS(I)

243 800 CONTINUE

244 850 I=NUMDIS

245 WR ITE ( IPR • 806 ) ALPHA

, ALPH AN » ALPHAF t ALPHA I » ALPH AN • ALPHAI

246 1ALPHAT » ALPHAY . CORR ( I ) IFLAG1 ( I ) » IFLAG2 ( I

t IFLAG3 ( I )

247 lYKI).YS(I)

248 C 249 C PLOT THE PROB PLOT CORR COEFFICIENT VERSUS GAMMA VALUE INDEX

250 C 251 CALL PLOT(CORR>AlNDEX»NUMDIS)

252 WRITE ( IPR » 810 ) ALPHAG. ALPHA A t ALPHAM# ALPMAM. ALPHAA t EQUAL

253 IGAMTAB(l) »GAMTAB(12) »GAMTAB(23> »GAMTAB<34) t 254 1ALPHAI»ALPHAN» ALPHAF » ALPHA I »ALPHAN» ALPHAI »ALPHAT» ALPHAY

255 WRITE(IPR»999)

256 WRITE(IPR#812)

257 WRITE(IPR»813)

258 C 259 C IF THE OPTIMAL GAMMA IS FINITF. PLOT OUT THE EXTREME VALUE

260 C TYPE 2 PROBABILITY PLOT FOR THE OPTIMAL VALUE

261 C OF GAMMA.

262 C 263 IF(IDISMX.LT.NUMDIS)CALL EV2PLT(X»N»GAMTAB( IDISMX)

264 C 265 C PLOT OUT AN EXTREME VALUE TYPF 1 PROBABILITY PLOT

266 C 267 CALL EV1PLT(X»N)

268 C 269 C FORM THE VARIOUS RETURN PERIOD VALUES

270 C 271 1650 K=0

272 D02100I=1»4

273 D02200J=1»9

274 K=K+1

275 AM(K)=J*(10**(I-1)>

276 2200 CONTINUE

277 2100 CONTINUE

278 K=K+1

279 AM(K)=10000.

280 K=K+1

281 AM(K)=50000.

282 K=K+1

283 AM(K)=100000.

284 K=K+1

285 AM(K)=500000.

286 K=K+1

287 AM(K)=1000000.

288 K=K+1

289 AM(K)=N
290 NUMAM=K

291 CALL SORT(AM.NUMAM»SCRAT>

292 0023001=1 »NUMAM

293 AM(I)=SCRAT(I)

294 2300 CONTINUE

295 C 296 C IF THE OPTIMAL GAMMA IS FlNITEr COMPUTE THE

297 C PREDICTED EXTREME <= F(l-(1/M)> FOR VARIOUS RETURN PERIODS M 298 C FOR THE OPTIMAL EXTREME VALUE TYPE 2 DISTRIBUTION.

299 C 300 IF(IDISMX.EQ.NUMDIS)GOT02450

301 A=GAMTAB(IDISMX)

302 YINT=YI(IDISMX)

303 YSLOPE=YS(IDISMX)

304 D02400I=2»NUMAM

305 R=1.0/AM(I)

306 P=1.0-R

307 ARG=-AL0G(P)

308 IF(ARG.LE.0.0)GOTO2400 \

309 H (

t 1 ) =YINT+YSLOPE* ( ARG** (-1 . 0/A)

)

310 2400 CONTINUE

311 C 312 C COMPUTE THE PREDICTED EXTREME (= F(l-(1/M)) FOR VARIOUS RETURN

313 C PERIODS M FOR THE EXTREME VALUE TYPE 1 DISTRIBUTION.

31* C 315 2450 YINT=YI(NUMDIS)

316 YSLOPE=YS(NUMDIS)

317 DO2500I=2»NUMAM

318 R=1.0/AM(I)

319 P=1.0-R

320 ARG=-ALOG«P)

321 IF(ARG.LE.O.O)GOT02500

322 H(I»2)=YINT+YSL0PE*(-AL0G(ARG)

323 2500 CONTINUE

324 C 325 C WRITE OUT THE PAGE WITH THE RETURN PERIODS AND THE PREDICTED EXTREMES

326 C FOR THE 2 DISTRIBUTIONS—OPTIMAL EXTREME VALUE TYPE 2» AND EXTREME

327 C VALUE TYPE 1. 328 C 329 WRITE(IPR»998)

330 IF(IDISMX.EQ.NUMDIS)GOT02750

331 WRITE* IPR»2602)

332 WRITE(IPR»2604>

333 WRITE(IPR»2606)

334 WRITE(IPR»2608)

335 WRITE ( IPR » 2610 ) GAMTAB < IDISMX

)

336 WRITE(IPR»999)

337 D02700I=2»NUMAM

338 WRITE(IPR.2705)AM(I)»H(I»1)»H(I»2)

339 J=I-1

340 JSKIP=J-5*(J/5)

341 IF(JSKIP.EQ.0)WRITE(IPR»999)

342 2700 CONTINUE

343 RETURN

344 C 345 2750 WRITE(IPRi2802)

346 WRITE(IPR»2804)

347 WRITE(IPR»2806)

29
348

349

350

351

352

353

354

355

356

357

358

359

360

361

362

363

364

365

366

367

368

369

370

371

372

373

374

375

376

377

378

379

380

381

382

383

384

385

386

387

388

389

390

391

392

393

WRITE(IPR»2808)

WRITE(IPR»999)

D02900I=2»NUMAM

WRITE <IPR '2705 >AM( I) »H(I»?)

J=I-1

JSKIP=J-5*(J/5)

IF(JSKIP.EQ.0)WRITE(IPR»999>

2900 CONTINUE

998

999

305

310

311

312

313

314

323

324

325

805

806

810

812

813

2602

2604

2606

2608

2610

2705

2802

2804

2806

2608

FORMAT

FORMAT

FORMAT

FORMAT

FORMAT

FORMAT

FORMAT

FORMAT

FORMAT

INFORMAT

IEFORMAT

FORMAT

FORMAT

FORMAT

115Xf8A

FORMAT

1LATION

FORMAT

FORMAT

1 27H

FORMAT

1 20H

FORMAT

1 27H

FORMAT

1 22H

FORMAT

FORMAT

FORMAT

FORMAT

FORMAT

FORMAT

RETURN

END

(1H1)

(1H ) (1H >40X»22HEXTREME VALUF ANALYSIS)

(1H »37X»20HTHE SAM°LE SIZE N = .17)

(1H »34X»18HTHE SAMPLE MEAN = »F14.7)

(1H »28X»32HTHE SAMPLE STANDARD DEVIATION = fF14.7)

<1H »32X»21HTHE SAMPLE MINIMUM = »F14.7)

(1H »32X»21HTHE SAMPLE MAXIMUM = .F14.7)

<1H »67H EXTREME VALUE

SCALE)

(1H r69H TYPE 2 TAIL LENGTH

ESTIMATE)

(1H »37H PARAMETER (GAMMA)

(1H »3X»F10.2»13X»FP.5f lXr3Al»2X»F14.7»2X»F14.7)

(1H »5X»8Al»13X»F8.5»lX.3Al»2X»F14.7»2XrF14.7)

(1H »12X»5Al»lX»AlrF14.7,liy,Fl(+.7»llX»F14.7»llX»F14.7.

1)(1H »96HTHE ABOVE IS A PLOT OF THE 43 PROBABILITY PLOT CORRE

COEFFICIENTS (FROM THE PREVIOUS PAGF)

(1H f 16X»41HVERSUS T HE q-3 EXTREME VALUE DISTRIBUTIONS)

PROBABILITY PLOT

CORRELATION

COEFFICIENT)

LOCATIO

ESTIMAT

(1H »43H RETURN PERIOD

PREDICTED EXTREME WIND)

(1H »43H (IN YEARS)

BASED ON)

(1H »42H EXTREME VALUE TYPE 1) (1H r43H

DISTRIBUTION)

(1H »30H

PREDICTED EXTREME WIND.

BASED ON OPTIMAL

EXTREME VALUE TYPE 2» DISTRIBUTION

(GAMMA = *F12.5»1H)>

(1H r2X.F9.1»13X»F10.2,17X»F10.2)

(1H »43H RETURN PERIOD PREDICTED EXTREME WIND)

36H (IN YEARS) RASED ON)

EXTREME VALUE TYPE 1) DISTRIBUTION)

(1H

(1H »42H

(1H »38H

OPRT.S SIMIU.SORT

30
JJF6*SIMIU.S0RT

12 C

3 C 1 C 5 C 6 C 7 C 8 C 9 C 10 C

11 C 12 C 13 C 11 C 15 C 16 C 17 C 18 c 19 c 20 c 21 c 22 c 23 c 21 c 25 c 26 c 27 c 28 c 29 c 30 c 31 c 32 c 33 c 31 c 35 c 36 c 37 c 38 c 39 c 10 c 11 c 12 c 13 c 11 c 15 c 16 c 17 c 18 c 19 c 50 c 51

52

53 c 51

55 c 56 c 57 c

subroutine sort(x.n.y)

this routine sorts the elements of the input vector x and puts the sorted

elements into the vector y. the input to this routine is the singlf precision vector x of

(unsorted) observations, the integer value n (= sample size),

and an empty single precision vector y into which the sorted observations

will be placed.

the output from this routine is the single precision vector y into which

the sorted observations have been placed.

restrictions on the maximum allowable value of n—the dimensions

of vectors iu and il (defined and used internally within this routine)

determine the maximum allowable value of n for this

routine. if iu and il each have dimension k» then n may not exceed

2**(k+1) - 1. for this routine as written* the dimensions of iu and il

have been set to 36. thus the maximum allowable value of n is

approximately 137 billion. since this exceeds the maximum allowable

value for an integer variable in many computers* and since a sort of 137

billion elements is presently impractical and unlikely* therefore no test for whether the input sample size n exceeds 137 billion has been

incorporated into this routine. it is thus assumed that there is no (practical) restriction on the maximum value of n for this routine,

printing—none unless an error condition exists

this routine is single precision in internal operation,

subroutines needed—none sorting method—binary sort reference—cacm march 1969* page 186 (binary sort algorithm by richard

c. singleton. —cacm january 1970. page 51. —cacm october 1970. page 621. —jacm january 1961* page 11.

the binary sort algorithm used herein is extremely fast as the

following time trials (performed by sorting random numbers)

on the univac 1108 exec 8 system indicate.

these time trials were carried out in august* 1971.

by way of comparison. the time trial values for the easy-to-program

but extremely inefficient bubble sort method have also been

included:

number of random binary sort • bubble sort

numbers sorted

.002 SEC .002 SEC

.011 SEC .015 SEC

.111 SEC 1.332 SEC

.176 SEC 37.683 SEC

1.887 SEC NOT COMPUTED

WRITTEN BY JAMES J. FILLIREN* STATISTICAL ENGINEERING LABORATORY (205.03)

NATIONAL BUREAU OF STANDARDS* WASHINGTON. D.C. 20231 JUNE 1972

DIMENSION X(l) *Y(1)

DIMENSION IU(36)»IL(36)

IPR=6

CHECK THE INPUT ARGUMENTS FOR ERRORS

N = 10 N = 100

N = 1000

N = 3000

N = 10000

31
58 IF(N.LT.DG0TO50

59 IF(N.EQ.DG0T055

60 HOLD=X(l)

61 D060I=2rN

62 IF(X(I) .NE.HOLD)GOT090

63 60 CONTINUE

61 WRITEdPR* 9)H0LD

65 D061I=1»N

66 Y(I)=X(I)

67 61 CONTINUE

68 RETURN

69 50 WRITE(IPR»15)

70 WRITE(IPR»17)N

71 RETURN

72 55 WRITE(IPR»18)

73 Y(1)=X<1)

71 RETURN

75 90 CONTINUE

76 9 FORMAT (1H . 10BH***** NON-FATAL DIAGNOSTIC —THE FIRST INPUT ARGUME

77 INT (A VECTOR) TO THE SORT SUBROUTINE HAS ALL ELEMENTS = »E15.8»6

78 1H *****)

79 15 FORMAT (1H » 91H***** FATAL ERROR —THE SECOND INPUT ARGUMENT TO THE

80 1 SORT SUBROUTINE IS NON-POSITIVE *****)

81 18 FORMAT (1H . 100H***** NON-FATAL DIAGNOSTIC—THE SECOND INPUT ARGUME

82 INT TO THE SORT SUBROUTINE HAS THE VALUE 1 *****)

83 17 FORMAT (1H » 35H***** THE VALUE OF THE ARGUMENT IS »I8 »6H *****)

81 C 85 C COPY THE VECTOR X INTO THE VECTOR Y 86 DO100I=l«N

87 Y(I)=X(I)

88 100 CONTINUE

89 C 90 C CHECK TO SEE IF THE INPUT VECTOR IS ALREADY SORTED

91 C 92 NM1=N-1

93 DO200I=lrNMl

91 IP1=I+1

95 IF(Y(I).LE.Y(IP1))GOT0200

96 GOT0250

97 200 CONTINUE

98 RETURN

99 250 M=l

100 1=1 101 J=N

102 305 IF(I.GE.J)GOT0370

103 310 K=I 104 MID=(I+J)/2

105 AMED=Y(MIO)

106 IF(Y(I) ,LE.AMED)GOT0320

107 Y(MID)=Y(I)

108 Y(I)=AMED

109 AMED=Y(MID)

110 320 L=J

111 IF(Y(J) ,GE.AMED)G0T0310

112 Y(MID)=Y(J)

113 Y(J)=AMED

111 AMED=Y(MID)

115 IF(Y(I) .LE.AMED)GOTO310

32
116 Y(MIO)=Y(I)

117 Y(I)=AMED

118 AMED=Y(MIO)

119 GOTO340

120 330 Y(L)=Y(K)

121 Y(K)=TT

122 340 L=L-1

123 IF(Y(L) .GT.AMFD)GOT0340

124 TT=Y(L)

125 350 K=K+1

126 IF(Y(K) .LT.AmED)GOTO350

127 IF(K.LE.L)GOT0330

128 LVI=L-I

129 JMK=J-K

130 IF ( LM I . LE . JMK ) GOTO360

131 IL(M)=I

132 IU(M)=L

133 I=K 134 M=M+1 y 135 GOTO3R0

136 360 IL(M)=K

137 IU(M)=J

138 J=L

139 M=M+1

140 GOT0380

mi 370 M=M-1

142 IF (M.EQ.O) RETURN

m3 I=IL(M)

144 J=IU(M)

145 380 JMI=J-I

146 IF(JMI.GE.ll)GOTO310

1*7 IF(I.EQ.DGOT0305

148 1 = 1-1 119 390 1 = 1+1 150 IF(I.EQ.J)GOT0370

151 AMED=Y(I + D

152 IF(Y(I).LE.AMED)GOTO390

153 K=I

154 395 Y(K+1)=Y(K)

155 K=K-1

156 IF(AMED.LT.YtK) )G0T0395

157 Y(K+1)=AMEP

158 GOTO390

159 END

BPRT»S SIMIU.UNIMED

33
JJF6*SIMiU.UNIMED

i subroutine: uuimfdcniX)

2 C

3 C THIS ROUTINE COMPUTES AN APP'jnx IMATION TO THE MEDIAN OF THE I-TH ORDER

4 C STATISTIC (For' I = 1 »?»... »N) FRO" t UMIFopM DISTRIBUTION (ON THE UNIT

5 C INTERVAL (0»D).

6 C THIS IS IDENTICAL TO THE "EDMN OF THE BETA DISTRIBUTION WITH PARAMETERS

7 C I AND N-I+l FOR I=1»2»...»N.

8 C THE INPUT TO THIS ROUTINE IS THE OFSIRFD TNTEGFR SAMPLE SIZE M

9 C AND AM EMPTY SINGLE PRECISION VECTOR X (OF DIMENSION «T LEAST N) INTO

10 C WHICH THE N GENERATFO UNIFORM ORDER STATISTIC MEDIANS WILL BE PLACED.

11 C THE OUTPUT FROM THIS ROUTINE TS THE SINGLE PRFCISION VECTOR X 12 C INTO WHICH THE N GENERATED UuIFOR" ORDF R STATISTIC MEDIANS

13 C HAVE BEEN PLACED.

It C ALL OF THE PROBABILITY PLOT ROIITINFS MAKE USE OF THIS ROUTINE.

15 C JUSTIFICATION AND ACCURACY OF THE ALGORITHM USED IS FO'UND IN AN

16 C UNPUBLISHED JJF MANUSCRIPT.

17 C THERE IS NO RESTRICTION ON T'lE MAXIMUM VALUE OF w FOR THIS ROUTINE.

18 C PRINTING—NONE UNLESS AN FRROd CONDITION EXISTS

19 C THIS ROUTINE IS SINGLE PRECISION IN INTERNAL OPERATION

20 C SUBROUTINES NEEDED—NONE

21 C RFFERENCE--UNPUBLISHED JJF MANUSCRIPT

22 C WRITTEN BY JAMES J. FILLlnENr STATISTICAL ENGINEERING LABORATORY (205.03)

23 C NATIONAL BUREAU OF STANDARDS* WASHINGTON D.C. 20234 JUNE 1972

24 C 25 DIMENSION X(l)

26 C 27 AN=N

28 IPR=6

29 C 30 C CHECK THE INPUT ARGUMENTS FOR ERRORS

31 C 32 IF(M.LT.l)GOTO50

33 IF(N.EQ.DG0T055

34 GOT090

35 50 WRITE(IPR» 5) 36 WRITE(IPR»47)N

37 RETURN

38 55 WRITE(IPR» 8) 39 90 CONTINUE

40 5 FORMAT (1H , vlH***** FATAL ERROR—THE FIRST INPUT ARGUMENT TO THE

41 1 IJNIMEO SUBROUTINE IS NON-POSITIVE *****)

42 fi FORl'ATdH , lOOh***** NON-FATAL DIAGNOSTIC —THF FIRST INPUT APGUME

43 INT TO THE UNIMEO SUBROUTINE HAS THE VALUE 1 *****)

44 47 FORMATdH . jSH***** THE VALUF OF THE ARGUMENT IS »I8 »6H *****)

45 C 46 X(N)=P.5**(1.0/AN)

47 X(l)=t.0-X(N)

48 NHALF=(N/?)+l

49 NEVOD0=2* (N/2) 50 I c ( N . NE . NEVOID ) X ( NHALF ) =0 .

51 IF(N.LE.3)RET')RN

52 GAM=U.3175

53 I«'AXrr'/2 54 D0100I=2»IVAa

55 AI = I 56 I^?EV=^|-I + l 57 X(I)=(AI-GAM)/ (AN-2.0*GAW+1.3)

34
58 X(IR£V)=1.0-X(I)

59 100 CONTINUE

60 RETURN

61 END

OPRT»S SIMIU.EV1PLT

35
JJF6*SIMIU.EV1PLT

1 SUBROUTINE EVlPLT(XrN)

2 C 3 C THIS ROUTINE GENERATES AN EXTREME VALUF TYPE 1 (EXPONENTIAL TYPE)

4 C PROBABILITY PLOT

5 C THE INPUT TO THIS ROUTINE IS THE SINGLE PRECISION VECTOR X OF

6 C (UNSORTED OR SORTED) OBSERVATIONS AND THE INTEGER VALUE N (= SAMPLE SIZE)

7 C THE OUTPUT FROM THIS ROUTINE IS A ONE-PAGE EXTREME VALUE TYPE 1 8 C PROBABILITY PLOT

9 C PRINTING—YES

10 C SUBROUTINES NFEDED—SORT» UNITED* AND PLOT

11 C REFERENCE—UNPUB. JJF MANUSCRIPT 'THE PERCENT POINT FUNCTIONS PAGE 26

12 C WRITTEN BY JAMES J. FILLIREN> STATISTICAL ENGINEERING LABORATORY (205.03)

13 C NATIONAL BUREAU OF STANDARDSr WASHINGTON* D.C. 20234 JUNE 1972

14 C 15 DIMENSION X(l)

16 DIMENSION Y (200

t W (200)

17 C 18 DATA TAU/1. 56166687/

19 C 20 AN=N

21 IPR=6

22 IUPPER=7500

23 C 24 C CHECK THE INPUT ARGUMENTS FOR ERRORS

25 C 26 IF(N.LT.l.OR.N.GT.IUPPER)GOTOS0

27 IF(N.EQ.DG0T055

26 HOLD=X(l)

29 D060I=2fN

30 IF(X(I).NE.HOLD)GOT090

31 60 CONTINUE

32 WRITE(IPR» 9)H0LD

33 GOT090

34 50 WRITE(IPR»17)IUPPER

35 WRITE(IPR»47)N

36 RETURN

37 55 WRITE(IPR»18)

38 RETURN

39 90 CONTINUE

40 9 FORMAT (1H »109H***** NON-FATAL DIAGNOSTIC —THE FIRST INPUT ARGUME

41 INT (A VECTOR) TO THE EV1PLT SUBROUTINE HAS ALL ELEMENTS = »E15.8»6

42 1H *****)

43 17 FORMAT (1H » 98H***** FATmL ERROR—THE SECOND INPUT ARGUMENT TO THE

44 1 EV1PLT SUBROUTINE IS OUTSIDE THE ALLOWABLE (1»»I6»16H) INTERVAL * 45 1****)

46 18 FORMAT! 1H »100H***** NON-FATAL DIAGNOSTIC —THE SECOND INPUT ARGUME

47 INT TO THE EV1PLT SUBROUTINE HAS THE VALUE 1 *****)

48 47 FORMATdH t 35H***** THE VALUF OF THE ARGUMENT IS »I8 »6H *****)

49 C 50 CALL SORT(X»N»Y)

51 CALL UNIMED(N»W)

52 DO100I=l»N

53 W(I)=-ALOG(ALOG(1.0/W(I) )

54 100 CONTINUE

55 CALL PLOT(Yri,VfN)

56 WRITE<IPR»105)TAU*N

57 SUM1=0.0

36
58 SUM2=0.0

59 DO200I=l»N

60 SUM1=SUM1+Y<1>

61 SUM2=SUM2+W(1)

62 200 CONTINUE

63 YBAR=SUM1/AN

6<* WBAR=SUM2/AN

65 SUM1=0.0

66 SUM2=0.0

67 SUM3=0.0

68 D0300I=1»N

69 SUM1=SUM1+ ( Y ( I

) -YBAR ) * { Y < I

> -YB AR

)

70 SUM2=SUM2+ ( Y ( I

) -YBAR ) * ( W ( I

> -WP AR

)

71 SUM3=SUM3+ ( W ( I

) -WBAR ) * ( W ( I

> -WB AR

)

72 300 CONTINUE > 73 CC=SUM2/SQRT(SUM3*SUM1>

7«* YSLOPE=SUM2/SUM3

75 YINT=YBAR-YSLOPE*WBAR

76 WRITE(IPRr305)COYINT»YSL0PE

77 105 FORMATdH »6<+HEXTREME VALUE TYPE 1 (EXPONENTIAL TYPE) PROBABILITY

78 1PLOT (TAU = .E15.8. 1H) #23X»20HTHE SAMPLE SIZE N = »I7)

79 305 FORMATdH »«*3HPROBABILITY PLOT CORRELATION COEFFICIENT = »F8.5»5X.

80 122HESTIMATED INTERCEPT = »E15.8»3X r 18HESTIMATED SLOPE = »E15.8>

81 RETURN • 82 END

0PRT»S SIMIU.EV2PLT

37
JJF6*SIMIU.EV2PLT

1 SUBROUTINE EV2PLT(X#N»GAMMA)

2 C 3 C THIS ROUTINE GENERATES AN EXTREME VALUE TYPE 2 (CAUCHY TYPE)

4 C PROBABILITY PLOT

5 C THE INPUT TO THIS ROUTINE IS THE SINGLE PRECISION VECTOR X OF

6 C (UNSORTED OR SORTED) OBSERVATIONS* THE INTEGER VALUE N (= SAMPLE SIZE)

7 C AND THE SINGLE PRECISION VALUE GAMMA (THE EXPONENT PARAMETER)

8 C THE OUTPUT FROM THIS ROUTINE IS A ONE-PAGE EXTREME VALUE TYPE 2 9 C PROBABILITY PLOT

10 C THE MAXIMUM ALLOWABLE VALUE OF N FOR THIS ROUTINE IS 7500

11 C PRINTING—YES

12 C SUBROUTINES NEEDED—SORT, UNIMED, AND PLOT

13 C REFERENCE—UNPUB. JJF MANUSCRIPT 'THE PERCENT POINT FUNCTION', PAGE 26

14 C WRITTEN BY JAMES J. FILLIBEN, STATISTICAL ENGINEERING LABORATORY (205.03)

15 C NATIONAL BUREAU OF STANDARDS, WASHINGTON* D.C. 20234 DECEMBER 1972

16 C 17 DIMENSION X(l)

18 DIMENSION Y (200 ) » W (200)

19 C 20 AN=N

21 IPR=6

22 IUPPER=7500

23 C 24 C CHECK THE INPUT ARGUMENTS FOR ERRORS

25 C 26 IF(N.LT.1.0R.N.GT.IUPPER)GOT050

27 IF(N.EQ.DG0T055

28 HOLD=X(l)

29 D060I=2»N

30 IF(X(I).NE.HOLD)GOT090

31 60 CONTINUE

32 WRITE(IPR» 9)H0LD

33 GOT090

34 50 WRITE(IPR»17)IUPPER

35 WRITE(IPR»47)N

36 RETURN

37 55 WRITEdPR.18)

38 RETURN

39 90 CONTINUE

40 9 FORMATUH »109H***** NON-FATAL DIAGNOSTIC—THE FIRST INPUT ARGUME

41 INT (A VECTOR) TO THE EV2PLT SUBROUTINE HAS ALL ELEMENTS = »E15.8»6

42 1H *****)

43 17 FORMAT (1H , 98H***** FATAL ERROR—THE SECOND INPUT ARGUMENT TO THE

44 1 EV2PLT SUBROUTINE IS OUTSIDE THE ALLOWABLE (1,»I6»16H) INTERVAL * 45 1****)

46 18 FORMAT (1H »100H***** NON-FATAL DIAGNOSTIC—THE SECOND INPUT ARGUME

47 INT TO THE EV2PLT SUBROUTINE HAS THE VALUE 1 *****)

48 47 FORMATUH » 35H***** THE VALUE OF THE ARGUMENT IS ,18 »6H *****)

49 C 50 CALL SORT(X,N,Y>

51 CALL UNIMED(N#W)

52 DO100I=l»N

53 W ( I ) = ( -ALOG ( W ( I )>)**( -1 . 0/GAMVA

54 100 CONTINUE

55 CALL PLOT(Y»W»N)

56 Q=.9975

57 PP9975=(-AL0G(Q) )**(-!. 0/GAWA)

38
58 Q=.0025

59 PP0025=(-ALOG(Q) )**(-1.0/GAMMA)

60 Q=.975

61 PP975 =(-ALOG(Q))**(-1.0/GAMMA)

62 Q=.025

63 PP025 =(-ALOG(Q))**(-1.0/GAMMA)

6H TAU= ( PP9975-PP0025 ) / ( PP975-PP025

)

65 WRITE(IPR»105)GAMMA»TAUrN . 66 SUM1=0.0 ' 67 SUM2=0.0

68 D0200I=1#N

69 SUM1=5UM1+Y(I)

70 SUM2=SUM2+W(I)

71 200 CONTINUE

72 YBAR=SUM1/AN

73 WBAR=SUM2/AN

71 SUM1=0.0

75 SUM2=0.0

76 SUM3=0.0

77 DO300I=l»N

78 SUM1=SUM1+ (Y(I) -YBAR ) * ( Y ( I

> -Y8AR

)

79 SUM2=SUM2+ < Y ( I

) -YBAR ) * ( W ( I

> -WBAR

)

80 SUM3=SUM3+ < W ( I ) -WBAR )*(W(D-WBAR)

81 300 CONTINUE

82 CC=SUM2/SQRT(SUM3*SUM1)

83 YSL0PE=SUM2/SUM3

84 YINT=YBAR-YSLOPE*WBAR

85 WRITE(IPR»305)CC»YINTrYSLOPE

86 105 FORMAT* 1H »63HEXTREME VALUE TYPE 2 (CAUCHY TYPE) PROB. PLOT WITH E 87 1XP. PAR. = »E17.10»1X»7H(TAU = »E15.8» 1H) »lXr 16HSAMPLE SIZE N = »I 88 17)

89 305 FORMATdH »«43HPROBABILITY PLOT CORRELATION COEFFICIENT = »F8.5»5X»

90 122HESTIMATED INTERCEPT = ,E15.8»3X»18HESTIMATED SLOPE = »E15.8)

91 RETURN

92 END

I3PRT.S SIMIU.PLOT

39
JJF6*SIMIU.PL0T

12 C 3 C

4 C

5 C 6 C 7 C 8 C 9 C 10 C 11 C 12 C 13 C 14 C 15 C 16

17

18

19

20

21 C 22

23

2i+ C 25 C 26 C 27

28 C 29

30

31

32

33

34

35

36

37

38

39

to

41

12

13

44

45

46

47

48

49

5051

52

53

54

55

56

57

SUBROUTINE PLOT(Y.X.N)

THIS ROUTINE YIELDS A ONE-PAGE PLOT OF Yd) VERSUS X(I). THE INPUT TO THIS ROUTINE IS THE SINGLE PRECISION VECTOR Y OF OBSERVATIONS* THE SINGLE PRECISION VECTOR X OF CORRESPONDING

OBSERVATIONSf AND THE INTEGER VALUF N (= SAMPLE SIZE).

MULTIPLE PLOT POINTS ARE NOT INDICATED AS SUCH. THERE IS NO RESTRICTION ON THE MAXIMUM VALUE OF N FOR THIS ROUTINE.

PRINTING—YES

SUBROUTINES NEEDED—NONE

WRITTEN BY JAMES J. FILLIBEN. STATISTICAL ENGINEERING LABORATORY (205.03)

NATIONAL BUREAU OF STANDARDS. WASHINGTON. D.C. 20234 JUNE 1972

UPDATED OCT 1974

UPDATED NOV 1974

INTEGER BLANK . HYPHEN . ALPHA 1 » ALPHAX INTEGER ALPHAM.ALPHAA.ALPHAD.ALPHAN. EQUAL

DIMENSION Y(1).X(1)

COMMON IGRAPH(55»130)

DIMENSION YLABLE(ll)

DATA BLANK. HYPHEN. ALPHAI . ALPHAX/lH .1H-.1HI.1HX/

DATA ALPHAM . ALPHAA . ALPHAD . ALPHAN . EQUAL/1HM. 1HA . 1HD . 1HN» 1H=/

CHECK THE INPUT ARGUMENTS FOR ERRORS

IPR=6

IF(N.LT.l)GOTO50

IF(N.EQ.DG0T055

HOLD=Y(l)

D060I=2.N

IF(Y(I) ,NE.H0LD)G0T065

60 CONTINUE

WRITEdPR. 9)H0LD

65 HOLD=X(l)

D070I=2.N

IF(X(I) .NE.HOLD)GOT090

70 CONTINUE

WRITEdPR. 19)H0LD

GOT090

50 WRITEdPR. 25)

WRITEdPR. 47)N

RETURN

55 WRITEdPR. 28)

RETURN

90 CONTINUE

9 F0RMAT(1H .10RH***** NON-FATAL DIAGNOSTIC—THE FIRST INPUT ARGUME

INT (A VECTOR) TO THE PLOT SUBROUTINE HAS ALL ELEMENTS = »E15.8»6

1H *****)

19 FORMAT (1H .108H***** NON-FATAL DIAGNOSTIC—THE SECOND INPUT ARGUME

INT (A VECTOR) TO THE PLOT SUBROUTINE HAS ALL ELEMENTS = »E15.8»6

1H *****)

25 FORMAT (1H . 91H***** FATAL ERROR—THE THIRD INPUT ARGUMENT TO THE

1 PLOT SUBROUTINE IS NON-POSITIVE *****)

28 FORMATdH . 100H***** NON-FATAL DIAGNOSTIC—THE THIRD INPUT ARGUME

INT TO THE PLOT SUBROUTINE HAS THE VALUE 1 *****)

40
58 47 FORMATUH t 35H***** THE VALUE OF THE ARGUMENT IS >I8 »6H *****)

59 C 60 C DETERMINE THE Y VALUES TO BE LISTED ON THE LEFT VERTICAL AXIS

61 C 62 YMIN=Y(1)

63 YMAX=Y(1)

64 D0105I = lfN

65 IF(Y(I).LT.YMIN)YMIN=Y(I)

66 IF(Y(I).GT.YMAX)YMAX=Y(I)

67 105 CONTINUE

68 DO110I=l»ll

69 AIM1=I-1

70 YLABLE(I)=YMAX-(AIM1/10.0)*(YMAX-YMIN)

71 110 CONTINUE

72 C 73 C DETERMINE XMIN» XMAX> XMIO» X25 (=THE 25* POINT) » AND 74 C X75 (=THE 75% POINT)

75 XMIN=X(1)

76 XMAX=X(1)

77 D0115I=1»N

78 IF(X(D .LT.XMIN)XMIN=X(I)

79 IF(X(I) .GT.XMAX)XMAX=X(I)

80 115 CONTINUE

81 XMID=(XMIN+XMAX)/2.0

82 X25=0.75*XMIN+0.25*XMAX

83 X75=0.25*XMIN+0.75*XMAX

84 C 85 C BLANK OUT THE GRAPH

86 C 87 D0100I=1»55

88 D0200J=1»129

89 IGRAPH(I»J)=BLANK

90 200 CONTINUE

91 100 CONTINUE

92 C 93 C PRODUCE THE Y AXIS

94 C 95 D0300I=3»43

96 IGRAPH(I»25)=ALPHAI

97 IGRAPH(I»129)=ALPHAI

98 300 CONTINUE

99 D0350I=3»43»4

100 IGRAPH(I»25)=HYPHEN

101 IGRAPH(I»129)=HYPHEN

102 350 CONTINUE

103 IGRAPH(3r21)=EQUAL

104 IGRAPH<3»22)=ALPHAM

105 IGRAPH(3»23)=ALPHAA

106 IGRAPH(3»24)=ALPHAX

107 IGRAPH(23»21)=EQUAL

108 IGRAPH(23»22>=ALPHAM

109 IGRAPH<23»23)=ALPHAI

110 IGRAPH<23»24)=ALPHAD

111 IGRAPH(43.21)=EQUAL

112 IGRAPH<43>22)=ALPHAM

113 IGRAPHU3t23)=ALPHAI

114 IGRAPH(43»2t*)=ALPHAN

115 C

41
116 C PRODUCE THE X AXIS

117 C 118 D0400J=27>127

119 IGRAPHd»J>=HYPHEN

120 IGRAPH<45»J)=HYPHEN

121 (+00 CONTINUE

122 D0450J=27d27»25

123 I6RAPH(1»J)=ALPHAI

124 IGRAPH(45»J)=ALPHAI

125 450 CONTINUE

126 D0460J=40»127#25

127 IGRAPH(1»J>=ALPHAI

128 IGRAPH(45»J>=ALPHAI

129 160 CONTINUE

130 C 131 C DETERMINE THE (X»Y) PLOT POSITIONS

132 C 133 RATIOY=40.0/(YMAX-YMIN)

134 RATIOX=100.0/(XMAX-XMIN)

135 D0600I=1»N

136 MX=RATIOX*(X(I)-XMIN)+0.5

137 MX=MX+27

138 MY=RATIOY*(Y(I)-YMIN)+0.5

139 MY=43-MY

140 IGRAPH(MY»MX)=ALPHAX

141 600 CONTINUE

112 C 143 C WRITE OUT THE GRAPH

144 C 145 WRITE(IPR»998)

146 D0700I=1»45

147 IP1=I+1

148 IFLAG=IPl-(IPl/4)*4

149 K=IPl/4

150 IF(IFLAG.NE.0)WRITE(IPR.705) (IGRAPH(I»J)»J=1»129)

151 IF(IFLAG.EQ.0)WRITE(IPR.706)YLABLE(K)f (IGRAPH(I»J)»J=21»129)

152 700 CONTINUE

153 WRITE (IPRr 707 ) XMIN»X25»XMID» X75»XMAX

154 705 FORMATdH rl29AD

155 706 FORMATdH »F20.7» 109A1

156 707 FORMATdH r 14X»F20 .7, 5XtF20.lt 5X»F20.7r 5X#F20.7» 1X»F20.7)

157 998 FORMAT* 1H1)

158 RETURN

159 END

I3PRT»S SIMIU.DATA4

42
MAXIMUM YEARLY WIND SPEEDS. CORPUS CHRISTI. TEXAS* 1912-1948

THE NUMBER OF OBSERVATIONS = 37

INPUT DATA

47.0 41.0 41.0 40.0 90.0 38.0 50.0 95.0 41.0 51.0 43.0 39.0 44.0 38.0 37.0 39.0

37.0 36.0 43.0 35.0 47.0 49.0 47.0 39.0 49.0 42.0 50.0 38.0 42.0 37.0 61.0 54.0

45.0 56.0 43.0 51.0 39.0

43

w —-^r
EXTREME VALUE ANALYSIS

THE SAMPLE SI2E N = 37 THE SAMPLE MEAN = 46.3243241

THE SAMPLE STANDARD DEVIATION = 12.7498453

THE SAMPLE MINIMUM = 35.0000000

THE SAMPLE MAXIMUM = 95.0000000

EXTREME VALUE PROBABILITY PLOT LOCATION SCALE

TYPE 2 TAIL LENGTH CORRELATION ESTIMATF ESTIMATE

PARAMETER (GAMMA) COEFFICIENT

1.00 .91022 40.9147968 1.2478256

2.00 .97191 MAX 31.0718093 9.3875747

3.00 .96594 20.9988256 19.2656157

4.00 .95b01 11.1366539 29.1725538

5.00 .94787 1.4036670 38.9815693

6.00 .94158 -8.2550049 48.7112250

7.00 .93668 -17.8684359 58.3850513

8.00 .93279 -27.4526310 68.0197020

9.00 .92965 -37.0169468 77.6262465

10.00 .92706 -46.5671854 87.2121058

11.00 .92489 -56.1070957 96.7823524

12.00 .92305 -65.6392174 106.3405638

13.00 .92147 -75.1653366 115.8893194

14.00 .92010 -84.6867085 125.4305038

15.00 .91891 -94.2042866 134.9655533

16.00 .91785 -103.7187805 144.4955559

17.00 .91691 -113.2307281 154.0213642

13.00 .91607 -122.7405472 163.5436363

19.00 .91531 -132.2485828 173.0629253

20.00 .91463 -141.7550869 182.5796452

21.00 .91400 -151.2603016 192.0941658

22.00 .91344 -160.7643795 201.6067638

23.00 .91292 -170.2674828 211.1176987

21.00 .91244 -179.7697258 220.6271515

25.00 .91200 -189.2712097 230.1353092

30.00 .91022 -236.7699699 277.6610641

35.00 .90894 -284.2588081 325.1694984

40.00 .90797 -331.7414551 372.6669960

45.00 .90721 -379.2199707 420.1571693

50.00 .90660 -426.6956291 467.6422195

60.00 .90569 -521.6410751 562.605 8524

70.00 .90503 -616.5816116 657.5525970

80.00 .90454 -711.5190506 752.4977341

90.00 .90415 -806.4543991 847.4390717

100.00 .90384 -901.3884811 942.3779907

150.00 .90291 -1376.0462189 1417.0503082

200.00 .90245 -1850.6961517 1891.7076111

250.00 .90217 -2325.3416443 2366.3575439

350.00 .90185 -3274.6292419 3315.6502075

500.00 .90161 -4698.5545654 4739.5795898

750.00 .90142 -7071.7604980 7112.7887573

1000.00 .90132 -9444.9562988 9485.9863281

INFINITY .90104 44. 0333295 9.4928209

44
>- _l

1- <

t- >

_I

11 111 CD I

<£ bj o m ir o o O 1- o o a. x o o a uj o o o o 10 m o o * a « • o

«H • U-l UJ 1 I I z 1- »-

•-» s. U. l/l II O 3

•o U)

(M IIK <r * O bj o < _l >

45
tj

m

3)X

O -i a J}

*> pi

CD 3 h rn rH <

-t >

-<

r OI c OI

T3m • r o O -1 o H -< o

T3 o o m o o o

33 fU o

33 II m — s

r- o t-l j» > z HC 1 i O OI' M X z -aj\o O —( o O -< 3> m tj IV) ti m CD T| — -

t-H r> -o h 33 m o

Z. 03 -H •

II T3

r-oH

-o s. -1 *+ - H

n£> I

h. I\> n • »H x ro

T3

O • IV) -• nmui

1/1 > 0* H *3 ->l

46
1 UJ hi

o_ ->

>- <_> t- >->

O u. * _l li- 00 < Id

Ch t-i

O

CM 1- (_) CO 2

ro Lui 2 X < • 2 O

•H

O t-i 1 1 0. 1- z X <

|-H hi -1

5"

w

u

II or o -< cc o o o UJ u o 0. o >- 1- o 1- o o _l • UJ 0. in z> 10 _l V < 1- > -1

_l

UJ >-< s: co Ul <

OH CD 1- O X CC

uj a.

47
RETURN PERIOD

(IN YEARS)

PREDICTED EXTREME WIND BASED ON OPTIMAL

EXTREME VALUE TYPE 2 DISTRIBUTION

(GAMMA = 2.00000)

PREDICTED EXTREME WIND

BASED ON EXTREME VALUE TYPE 1 DISTRIBUTION

2.0

3.0

1.0

5.0

6.0

42.35

45.81

48.57

50.94

53.06

44.51

49.60

52.86

55.27

57.19

7.0

8.0

9.0

10.0

20.0

54.98

56.76

58.43

59.99

72.52

58.78

60.15

61.34

62.40

69.23

30.0

37.0

40.0

50.0

60.0

82.06

87.79

90.07

97.12

103.48

73.16

75.18

75.93

78.07

79.82

70.0

80.0

90.0

100.0

200.0

109.33

114.77

119.88

124.71

163.67

81.30

82.57

83.70

84.70

91.31

300.0

400.0

500.0

600.0

700.0

193.53

218.71

240.68

260.92

279.36

95.16

97.90

100.02

101.75

103.21

800.0

900.0

1000.0

2000.0

3000.0

296.51

312.62

327.86

450.85

545.21

104.48

105.60

106.60

113.19

117.04

4000.0

5000.0

6000.0

7000.0

8000.0

624.76

694.86

758.23

816.50

870.73

119.77

121.89

123.62

125.08

126.35

9000.0

10000.0

50000.0

100000.0

500000.0

921.66

969.87

2130.33

2999.87

6674.48

127.47

128.47

143.74

150.32

165.62

1000000.0 9426.27 172.20

48

NBS-114A (REV. 7-73)

U .S. DEPT. OF COMM. BIBLIOGRAPHIC DATA

SHEET

1. PUBLICATION OK REPORT NO.

::iir tw-868

2. Gov't Accession

No.

3. Recipient's Accession No.

4. TITLE AND SUBTITLE

Statistical Analysis of Extreme Winds

5. Publication Date June 1975

6. Performing Organization Code

7. AUTHOR(S)

Emil Simiu and James J. Filliben

8. Performing Organ. Report No.

9. PERFORMING ORGANIZATION NAME AND ADDRESS

NATIONAL BUREAU OF STANDARDS

DEPARTMENT OF COMMERCE

WASHINGTON, D.C. 20234

10. Project/Task/Work Unit No.

11. Contract /Grant No.

12. Sponsoring Organization Name and Complete Address (Street, City, State, ZIP)

Same as 9.

13. Type of Report & Period Covered Final

14. Sponsoring Agency Code

15. SUPPLEMENTARY NOTES

Library of Congress Catalog Card Number: 75-600028

16. ABSTRACT (A 200-word or less factual summary of most significant information. If document includes a significant

bibliography or literature survey, mention it here.) With a view to assessing the validity of current probabilistic approaches to the

definition of design wind speeds, a study was undertaken of extreme wind speeds

based on records taken at 21 U.S. weather stations. For the purpose of analyzing

extreme value data, a computer program was developed, which is described herein.

The following results were obtained: (1) the assumption that a single probability

distribution is universally applicable to all extreme wind data sets in a given

type of climate was not confirmed, and (2) predictions of 100-year wind speeds

based on overlapping 20-year sets of data taken at the same station differed

between themselves by as much as 100%. Similar predictions for 1000-year winds

differed by as much as a few hundred percent. Since wind pressures are proportional

to the square of the wind speeds, errors of such magnitude are unacceptably high

for structural design purposes. It is therefore suggested that while, in principle,

probabilisitc methods provide the most rational approach to specifying design wind

speeds, it is of the utmost importance that the possible errors inherent in this

approach be carefully taken into account.

17. KEY WORDS (six to twelve entries; alphabetical order; capitalize only the first letter of the first key word unless a proper name; separated by semicolons) Building code s ; extreme value distr ibutions ; hurricanes;

probability distribution functions; reliability; risk; statistical analysis; storms;

structural engineering; wind loads; wind speeds

18. AVAILABILITY [)( Unlimited

I' or Official Distribution. Do Not Release to NTIS

[^Cj Order From Sup. of Doc, U.S. Government Printing Office Washington, D.C. 20402, SD Cat. No. Cl\ 46:868

_ Order From National Technical Information Service (NTIS)

Springfield, Virginia 22IM

IV. SECURITY CLASS

(THIS REPORT)

UNCLASSIFIED

20. SECURITY CLASS

(THIS PAGE)

UNCLASSIFIED

21. NO. OF PAGES

52

22. Price

$1.50

USCOMM-DC 29042-P74
NBS TECHNICAL PUBLICATIONS

PERIODICALS JOURNAL OF RESEARCH reports National Bureau

of Standards research and development in physics,

mathematics, and chemistry. It is published in two sections, available separately:

• Physics and Chemistry (Section A)

Papers of interest primarily to scientists working in these fields. This section covers a broad range of physical and chemical research, with major emphasis on

standards of physical measurement, fundamental constants, and properties of matter. Issued six times a

year. Annual subscription: Domestic, $17.00; Foreign,

$21.25.

• Mathematical Sciences (Section B)

Studies and compilations designed mainly for the mathematician and theoretical physicist. Topics in mathematical statistics, theory of experiment design, numerical analysis, theoretical physics and chemistry, logical design and programming of computers and computer

systems. Short numerical tables. Issued quarterly. Annual subscription: Domestic, $9.00; Foreign, $11.25. DIMENSIONS/NBS (formerly Technical News Bulletin) —This monthly magazine is published to inform

scientists, engineers, businessmen, industry, teachers,

students, and consumers of the latest advances in science and technology, with primary emphasis on the work at NBS. The magazine highlights and reviews such

issues as energy research, fire protection, building tech- nology, metric conversion, pollution abatement, health and safety, and consumer product performance. In addition, it reports the results of Bureau programs in measurement standards and techniques, properties of matter and materials, engineering standards and services, instrumentation, and automatic data processing. Annual subscription: Domestic, $9.45; Foreign, $11.85.

NONPERIODICALS

Monographs—Major contributions to the technical liter- ature on various subjects related to the Bureau's scientific and technical activities. Handbooks—Recommended codes of engineering and

industrial practice (including safety codes) developed

in cooperation with interested industries, professional

organizations, and regulatory bodies.

Special Publications—Include proceedings of confer- ences sponsored by NBS, NBS annual reports, and other

special publications appropriate to this grouping such

as wall charts, pocket cards, and bibliographies.

Applied Mathematics Series—Mathematical tables, manuals, and studies of special interest to physicists, engineers, chemists, biologists, mathematicians, computer programmers, and others engaged in scientific and technical work.

National Standard Reference Data Series—Provides

quantitative data on the physical and chemical properties of materials, compiled from the world's literature and critically evaluated. Developed under a world-wide

•program coordinated by NBS. Program under authority

of National Standard Data Act (Public Law 90-396). NOTE: At present the principal publication outlet for these data is the Journal of Physical and Chemical

Reference Data (JPCRD) published quarterly for NBS

by the American Chemical Society (ACS) and the American Institute of Physics (AIP). Subscriptions, reprints, and supplements available from ACS, 1155 Sixteenth

St. N. W., Wash. D. C. 20056.

Building Science Series—Disseminates technical infor- mation developed at the Bureau on building materials, components, systems, and whole structures. The series presents research results, test methods, and performance criteria related to the structural and environmental functions and the durability and safety characteristics of building elements and systems.

Technical Notes—Studies or reports which are complete

in themselves but restrictive in their treatment of a

subject. Analogous to monographs but not so comprehensive in scope or definitive in treatment of the subject area. Often serve as a vehicle for final reports of work performed at NBS under the sponsorship of other government agencies.

Voluntary Product Standards—Developed under pro- cedures published by the Department of Commerce in Part 10, Title 15, of the Code of Federal Regulations.

The purpose of the standards is to establish nationally

recognized requirements for products, and to provide

all concerned interests with a basis for common under- standing of the characteristics of the products. NBS

administers this program as a supplement to the activities of the private sector standardizing organizations. Federal Information Processing Standards Publications (FIPS PUBS)—Publications in this series collectively constitute the Federal Information Processing Standards Register. Register serves as the official source of information in the Federal Government regarding stand- ards issued by NBS pursuant to the Federal Property and Administrative Services Act of 1949 as amended,

Public Law 89-306 (79 Stat. 1127), and as implemented

by Executive Order 11717 (38 FR 12315, dated May 11, 1973) and Part 6 of Title 15 CFR (Code of Federal

Regulations). Consumer Information Series—Practical information, based on NBS research and experience, covering areas

of interest to the consumer. Easily understandable language and illustrations provide useful background

knowledge for shopping in today's technological

marketplace.

NBS Interagency Reports (NBSIR)—A special series of interim or final reports on work performed by NBS for

outside sponsors (both government and non-govern- ment). In general, initial distribution is handled by the sponsor; public distribution is by the National Technical Information Service (Springfield, Va. 22161) in paper

copy or microfiche form.

Order NBS publications (except NBSIR's and Biblio- graphic Subscription Services) from: Superintendent of Documents, Government Printing Office, Washington,

D.C. 20402. BIBLIOGRAPHIC SUBSCRIPTION SERVICES

The following current-awareness and literature-survey

bibliographies are issued periodically by the Bureau:

Cryogenic Data Center Current Awareness Service A literature survey issued biweekly. Annual subscription: Domestic, $20.00; foreign, $25.00.

Liquefied Natural Gas. A literature survey issued quarterly. Annual subscription: $20.00. Superconducting Devices and Materials. A literature

survey issued quarterly. Annual subscription: $20.00. Send subscription orders and remittances for the pre- ceding bibliographic services to National Technical Information Service, Springfield, Va. 22161.

Electromagnetic Metrology Current Awareness Service Issued monthly. Annual subscription: $100.00 (Special rates for multi-subscriptions). Send subscription order and remittance to Electromagnetics Division, National Bureau of Standards, Boulder, Colo. 80302.
U.S. DEPARTMENT OF COMMERCE

National Bureau of Standards

Washington. D.C. 20234

OFFICIAL BUSINESS

Penalty for Private Use, $300

POSTAGE AND FEES PAID

U.S. DEPARTMENT OF COMMERCE COM-215

SPECIAL FOURTH-CLASS RATE

BOOK