xml version 1.0 encoding UTF8 standalone no
record xmlns http:www.loc.govMARC21slim xmlns:xsi http:www.w3.org2001XMLSchemainstance xsi:schemaLocation http:www.loc.govstandardsmarcxmlschemaMARC21slim.xsd
leader nam Ka
controlfield tag 001 001441459
003 fts
006 med
007 cr mnuuuuuu
008 031203s2003 flua sbm s0000 eng d
datafield ind1 8 ind2 024
subfield code a E14SFE0000132
035
(OCoLC)53961818
9
AJM5899
b SE
SFE0000132
040
FHM
c FHM
090
QC21.2
1 100
Looper, Jason K.
0 245
Semiparametric estimation of unimodal distributions
h [electronic resource] /
by Jason K. Looper.
260
[Tampa, Fla.] :
University of South Florida,
2003.
502
Thesis (M.S.)University of South Florida, 2003.
504
Includes bibliographical references.
516
Text (Electronic thesis) in PDF format.
538
System requirements: World Wide Web browser and PDF reader.
Mode of access: World Wide Web.
500
Title from PDF of title page.
Document formatted into pages; contains 93 pages.
520
ABSTRACT: One often wishes to understand the probability distribution of stochastic data from experiment or computer simulations. However, where no model is given, practitioners must resort to parametric or nonparametric methods in order to gain information about the underlying distribution. Others have used initially a nonparametric estimator in order to understand the underlying shape of a set of data, and then later returned with a parametric method to locate the peaks. However they are interested in estimating spectra, which may have multiple peaks, where in this work we are interested in approximating the peak position of a singlepeak probability distribution. One method of analyzing a distribution of data is by fitting a curve to, or smoothing them. Polynomial regression and leastsquares fit are examples of smoothing methods. Initial understanding of the underlying distribution can be obscured depending on the degree of smoothing. Problems such as under and oversmoothing must be addressed in order to determine the shape of the underlying distribution.Furthermore, smoothing of skewed data can give a biased estimation of the peak position. We propose two new approaches for statistical mode estimation based on the assumption that the underlying distribution has only one peak. The first method imposes the global constraint of unimodality locally, by requiring negative curvature over some domain. The second method performs a search that assumes a position of the distribution's peak and requires positive slope to the left, and negative slope to the right. Each approach entails a constrained leastsquares fit to the raw cumulative probability distribution.We compare the relative efficiencies [12] of finding the peak location of these two estimators for artificially generated data from known families of distributions Weibull, beta, and gamma. Within each family a parameter controls the skewness or kurtosis, quantifying the shapes of the distributions for comparison. We also compare our methods with other estimators such as the kerneldensity estimator, adaptive histogram, and polynomial regression. By comparing the effectiveness of the estimators, we can determine which estimator best locates the peak position. We find that our estimators do not perform better than other known estimators. We also find that our estimators are biased. Overall, an adaptation of kernel estimation proved to be the most efficient.The results for the work done in this thesis will be submitted, in a different form, for publication by D.A. Rabson and J.K. Looper.
590
Adviser: Ph.D, David A. Rabson
653
probability.
single peak constrained least squares.
bumphunting.
690
Dissertations, Academic
z USF
x Physics
Masters.
773
t USF Electronic Theses and Dissertations.
4 856
u http://digital.lib.usf.edu/?e14.132
PAGE 1
Semiparametric Estimation of Unimo dal Distributions b y Jason K. Lo op er A thesis submitted in partial fulllmen t of the requiremen ts for the degree of Master of Science in Ph ysics Departmen t of Ph ysics College of Arts and Science Univ ersit y of South Florida Ma jor Professor: Da vid A. Rabson, Ph.D. W ei Chen, Ph.D. Ch unMin Lo, Ph.D. Date of Appro v al: August 20, 2003 Keyw ords: Probabilit y Bumph un ting, Single P eak, Constrained Least Squares c r Cop yrigh t 2003, Jason K. Lo op er
PAGE 2
DEDICA TION T o m y father who nev er stopp ed encouraging me, to m y other father (Stephen D. Smith) who nev er stopp ed funding me, and to m y grann y Lo op er who nev er lost her sense of w onder.
PAGE 3
A CKNO WLEDGMENTS Without the instruction and sup ervision of Dr. Da vid Rabson this thesis could not ha v e b een accomplished, and man y results could not ha v e b een concluded.
PAGE 4
T ABLE OF CONTENTS LIST OF T ABLES iii LIST OF FIGURES iv CHAPTER 1 INTR ODUCTION AND MOTIV A TION 1 1.1 In tro duction 1 1.2 Motiv ation 5 1.3 Researc h Plan 6 CHAPTER 2 P ARAMETRIC AND NONP ARAMETRIC ESTIMA TORS 7 2.1 Kernel Densit y Estimator 7 2.2 Optimal Kernel Densit y Estimator 10 2.3 Adaptiv e Histogram 10 2.4 Optimal Cheb yshev P olynomial Fit 11 2.5 Un biased P arametric Estimator 13 CHAPTER 3 SEMIP ARAMETRIC ESTIMA TORS 22 3.1 Unimodal3 22 3.2 Unimodal2 24 CHAPTER 4 MODEESTIMA TION EFFICIENCY 30 4.1 Setting a Scale 30 4.2 Mirror Distributions 31 4.3 STD of the STD as an Error assessmen t 32 4.4 Choice of 34 CHAPTER 5 RESUL TS F OR DISTRIBUTIONS HA VING UNIT MEAN 35 5.1 W eibull and Gamma Estimators 38 5.2 OKS and OCPF 38 5.3 FKS and AHIST 39 5.4 Unimo dal3 and Unimo dal2 39 CHAPTER 6 RESUL TS F OR BET A DISTRIBUTIONS HA VING NEGA TIVE AND POSITIVE KUR TOSIS 50 6.1 Beta Estimator 50 6.2 OKS and OCPF 51 6.3 FKS and AHIST 53 6.4 Unimodal3 and Unimodal2 54 i
PAGE 5
CHAPTER 7 CONCLUSION AND FUTURE DIRECTIONS 66 REFERENCES 67 APPENDICES 68 App endix A Momen ts, Prop erties and Cum ulan ts 69 App endix B T est Probabilit y Distributions 73 App endix C Cheb yshev P oly omial 77 App endix D Finite Dierences 79 ii
PAGE 6
LIST OF T ABLES T able 1. W eibull P arameters for = 1 36 T able 2. Gamma P arameters for = 1 36 T able 3. Beta P arameters for Flat P eaks 52 T able 4. Beta P arameters for Sharp P eaks 52 iii
PAGE 7
LIST OF FIGURES Figure 1. Probabilit y Distribution 2 Figure 2. Cum ulativ e Distribution 3 Figure 3. Ra w Data Set and Ra w Cum ulativ e Distribution 4 Figure 4. Ov ersmo othed and Undersmo othed Distribution b y Kernel 8 Figure 5. Double P eak Distribution Smo othed b y Kernel 9 Figure 6. Ov er and Undersmo othing with AHIST 11 Figure 7. AHIST with T en P ercen t Bin P arameter 12 Figure 8. OCPF with Negativ e Curv ature at Start of Probabilit y Distribution 18 Figure 9. OCPF to a Probabilit y Distribution 19 Figure 10. Absolute Error on a Gamma Distribution With Lo w Sk ewness 20 Figure 11. Absolute Error on a Gamma Distribution With High Sk ewness 21 Figure 12. P olynomial t to cum ulativ e with 1st and 3rd deriv ativ es 23 Figure 13. Ra w Cum ulativ e Distribution After Unimodal3 's FirstStage Estimator 24 Figure 14. Constrained LeastSquares Fit to Filtered Ra w Cum ulativ e Distribution and Deriv ativ e 25 Figure 15. P olynomial Fit to Cum ulativ e With 1st and 2nd Deriv ativ es 27 Figure 16. Before and After Unimodal2 's FirstStage Estimation 28 Figure 17. Exhaustiv e Searc h on Decimated and Binary on Undecimated on 200 Deviate Samples Dra wn from W eibull and Gamma 28 Figure 18. Exhaustiv e Searc h on Decimated and Binary on Undecimated on 200 Deviate Samples Dra wn from Beta 29 Figure 19. Exhaustiv e and Binary Searc h on Decimated for 200 and 1000Deviate Samples Dra wn from Gamma 29 iv
PAGE 8
Figure 20. Sharp and Flat P eaks 31 Figure 21. Mirror Images of Sk ew ed Distribution 32 Figure 22. Gamma Results for 1000 Deviates 37 Figure 23. Gamma Results for 200 Deviates 41 Figure 24. W eibull Results for 200 Deviates 42 Figure 25. W eibull Results for 1000 Deviates 43 Figure 26. Gamma Estimator on a W eibull Distribution 44 Figure 27. Optimal Kernel on 1000 and 200 Deviate Sets 45 Figure 28. Kernel and Adaptiv e Histogram Estimator 46 Figure 29. Bias of Unimodal3 46 Figure 30. Histogram of Bias of Unimodal3 on a Gamma Distribution 47 Figure 31. RMS error VS Num b er of Deviates for Unimodal3 47 Figure 32. Bias of Unimodal3 on W eibull Distribution 48 Figure 33. Bias of Unimodal2 on W eibull and Gamma Distributions 48 Figure 34. Histogram of Bias of Unimodal2 on Gamma Distributions 49 Figure 35. Histogram of Biasness of Unimodal2 on W eibull Distributions 49 Figure 36. Histogram of Biasness of Beta estimator on Beta Distributions p = 1 : 1 53 Figure 37. Histogram of Bias of Beta Estimator on Beta Distributions p = 1 : 8 54 Figure 38. Beta p = 1 : 1 Results for Sharp P eaks With 1000 Deviates 56 Figure 39. Beta p = 1 : 1 Results for Flat P eaks With 200 Deviates 57 Figure 40. Beta p = 1 : 8 Results for Sharp P eaks With 1000 Deviates 58 Figure 41. Beta p = 1 : 8 Results for 200 Flat P eaks With 200 Deviates 59 Figure 42. OCPF on a Beta Distribution p = 1 : 1 and q = 1 : 1 60 Figure 43. OKS with Large Negativ e Kurtosis 60 Figure 44. Ra w Deviates for Beta p = 1 : 1 and q = 1 : 1 61 Figure 45. Unimodal3 's FirstStage Estimation on The Beta p = 1 : 1 and q = 1 : 1 62 v
PAGE 9
Figure 46. Unimodal2 's Bias on Beta p = 1 : 8 62 Figure 47. Unimodal2 's Bias on Beta p = 1 : 1 63 Figure 48. Unimodal2 's Bias Ov er All Sk ewness on Beta F amilies 63 Figure 49. Unimodal3 's Bias on Beta p = 1 : 1 64 Figure 50. Unimodal3 's Bias on Beta p = 1 : 8 64 Figure 51. Unimodal3 's Bias Ov er All Sk ewness of Beta F amilies 65 Figure 52. W eibull Distribution Sk ewness Examples 71 vi
PAGE 10
SEMIP ARAMETRIC ESTIMA TION OF UNIMOD AL DISTRIBUTIONS Jason Keith Lo op er ABSTRA CT One often wishes to understand the probabilit y distribution of sto c hastic data from exp erimen t or computer sim ulations. Ho w ev er, where no mo del is giv en, practitioners m ust resort to parametric or nonparametric metho ds in order to gain information ab out the underlying distribution. Others ha v e used initially a nonparametric estimator in order to understand the underlying shap e of a set of data, and then later returned with a parametric metho d to lo cate the p eaks. Ho w ev er they are in terested in estimating sp ectra, whic h ma y ha v e m ultiple p eaks, where in this w ork w e are in terested in appro ximating the p eak p osition of a singlep eak probabilit y distribution. One metho d of analyzing a distribution of data is b y tting a curv e to, or smo othing them. P olynomial regression and leastsquares t are examples of smo othing metho ds. Initial understanding of the underlying distribution can b e obscured dep ending on the degree of smo othing. Problems suc h as under and o v ersmo othing m ust b e addressed in order to determine the shap e of the underlying distribution. F urthermore, smo othing of sk ew ed data can giv e a biased estimation of the p eak p osition. W e prop ose t w o new approac hes for statistical mo de estimation based on the assumption that the underlying distribution has only one p eak. The rst metho d imp oses the global constrain t of unimo dalit y lo cally b y requiring negativ e curv ature o v er some domain. The second metho d p erforms a searc h that assumes a p osition of the distribution's p eak and requires p ositiv e slop e to the left, and negativ e slop e to the righ t. Eac h approac h en tails a constrained leastsquares t to the ra w cum ulativ e probabilit y distribution. 1 1 Refer to c hapter one for the ra w cum ulativ e probabilit y distribution. vii
PAGE 11
W e compare the relativ e eciencies [12 ] of nding the p eak lo cation of these t w o estimators for articially generated data from kno wn families of distributions W eibull, b eta, and gamma. Within eac h family a parameter con trols the sk ewness or kurtosis, quantifying the shap es of the distributions for comparison. W e also compare our metho ds with other estimators suc h as the k erneldensit y estimator, adaptiv e histogram, and p olynomial regression. By comparing the eectiv eness of the estimators, w e can determine whic h estimator b est lo cates the p eak p osition. W e nd that our estimators do not p erform b etter than other kno wn estimators. W e also nd that our estimators are biased. Ov erall, an adaptation of k ernel estimation pro v ed to b e the most ecien t. The results for the w ork done in this thesis will b e submitted, in a dieren t form, for publication b y D.A. Rabson and J.K. Lo op er. viii
PAGE 12
CHAPTER 1 INTR ODUCTION AND MOTIV A TION 1.1 In tro duction W e are in terested in the p eak p osition of a singlep eak, con tin uous, univ ariate probabilit y distribution. A probabilit y distribution, p ( x ), on some domain [ a; b ], denes the probabilit y p ( x ) dx that a random v ariable will b e measured in the range [ x; x + dx ]. There are t w o constrain ts: p ( x ) 0 (1) Z b a p ( x ) dx = 1 (2) Fig. 1 sho ws an example of a probabilit y distribution. If w e ha v e a con tin uous distribution, p ( x ), w e can describ e the con tin uous cum ulativ e distribution, C ( x ), b y C ( x ) = Z x 1 p ( x 0 ) dx 0 : (3) An example of a con tin uous cum ulativ e distribution is sho wn in Fig. 2 If w e ha v e a set of random measuremen ts, x i tak en from a distribution, p ( x ), w e can use the discrete form of the cum ulativ e distribution C ( x ) = X x i
PAGE 13
Figure 1. Probabilit y distribution. in the discrete cum ulativ e distribution. The discrete cum ulativ e distribution is another w a y of describing ho w m uc h probabilit y is at, or to the left of, eac h datum. 1 It is not dicult to see that when the data are group ed close together, the ra w cum ulativ e distribution will ha v e a large slop e. Where the slop e is largest corresp onds to the p eak p osition in the probabilit y distribution. W e m ust rst smo oth the ra w cum ulativ e distribution in order to estimate the p eak in the probabilit y distribution. In analyzing data w e are often in terested in determining the underlying distribution asso ciated with observ ed or articially pro duced measuremen ts. Giv en a nite set of measuremen ts, w e can emplo y sev eral dieren t statistical metho ds in order to estimate the p eak of the distribution. These metho ds can b e used individually or together. The relativ e eciency of these metho ds dep ends on our kno wledge of the system b eing measured, or of the articially pro duced data. The t w o statistical metho ds that will b e discussed are parametric and nonparametric estimators. A com bination of these metho ds to lo cate a p eak p osition of a unimo dal distribution is the main topic of this w ork. 1 Random measuremen ts from a distribution will b e referred to as ra w data sets, and the discrete cum ulativ e distribution will b e referred to as the ra w cum ulativ e distribution. 2
PAGE 14
Figure 2. Cum ulativ e distribution. Often w e ha v e partial or no kno wledge of the underlying functional form of the distribution asso ciated with a data set. In this situation, nonparametric metho ds are generally used b ecause they imp ose no restrictions on the data. Nonparametric estimators can b e eectiv e on an y family of distributions, ev en when the distribution is not kno wn. The coun terparts are the parametric metho ds, in whic h w e assume a functional form with a small n um b er of parameters. The estimation then consists of nding the b est t for these parameters. If w e kno w something ab out the distribution, w e can use a h ybrid tec hnique, whic h consists of b oth parametric and nonparametric metho ds. Riedel prop osed a piecewise con v ex metho d of estimating the shap e of an unkno wn function, whic h uses a t w o stage smo othing tec hnique [9 ]. The rst stage uses a nonparametric smo othing metho d for determining the lo cation of inrection p oin ts and then a secondstage metho d to smo oth the function b et w een these p oin ts. Riedel's metho d parallels our prop osal in that w e also use a rststage estimation to determine the domain of a p eak and then a nalstage estimation to smo oth. Riedel is also using an initial smo othing metho d to gain information ab out a sp ectrum of m ultiple p eaks, where w e already understand that our distributions ha v e only one p eak. Unlik e Reidel's metho d, 3
PAGE 15
Figure 3. The gure on the left is a set of random measuremen ts tak en form a b eta distribution ha ving p = 1 : 8 and q = 5. The set to the righ t is its discrete (ra w) cum ulativ e distribution.our rststage estimator is used for lo cation of a p eak and not to determine the n um b er of p eaks. Riedel's prop osal also diers from ours in the nal metho d of smo othing [10 ]. W e prop ose t w o dieren t metho ds for determining the p eak p osition of an unkno wn probabilit y distribution. The distributions w e are in terested in estimating are unimo dal and therefore con tain no more than t w o inrection p oin ts. The rst metho d, named Unimodal3 uses an estimation tec hnique to determine the domain of the p eak p osition. Once this domain is found, a constrained leastsquares t is formed in that domain on the ra w cum ulativ e distribution, where the constrain t is on the third deriv ativ e of the t [3 ]. Rasm us, Nic holaos and Sidirop oulos use w eigh ted least squares, ho w ev er they are also in terested in estimating sp ectra [8]. The second prop osed metho d, named Unimodal2 p erforms an exhaustiv e or binary searc h for the p eak p osition, sub ject to the condition of p ositiv e slop e to the left, and negativ e slop e to the righ t of the p eak lo cation. Once the t has b een estimated, the residuals from the tted cum ulativ e distribution can b e compared for ev ery prop osed p eak lo cation. 2 2 Residuals are the square dierence of the tted cum ulativ e distribution and the ra w cum ulativ e distribution. 4
PAGE 16
1.2 Motiv ation One motiv ation for this pro ject came from mo deling nite nanometerscale quan tum wires. It is found that wires with a nite n um b er of hopping sites act unlik e systems with an innite n um b er of sites. In innite c hains, electronic transp ort is b eliev ed to undergo an immediate transition from ballistic conduction to diusiv e when innitesimal secondneigh b or coupling is turned on. In small c hains, Rabson, Narozhn y and Millis no w ha v e n umerical evidence instead for a crosso v er region of the lev elspacing probabilit y distribution, P ( E ), where is the mean lev el spacing [7 ]. The Hamiltonian with only rstneigh b or in teractions is in tegrable. This is due the large n um b er of conserv ed quan tities whic h determine the b eha vior of the system. Here energy lev els can generically cross. 3 Therefore the lev elspacing distributions for the rstneigh b or in teraction follo w P oisson statistics, P ( E ) = 1 e ( E = ) : (5) When w e turn on secondneigh b or in teraction, w e break the in tegrabilit y of the system b y reducing the n um b er of conserv ed quan tities. Here the lev el spacings generically do not cross, so that small and large lev el spacings b ecome unlik ely The lev el spacing with the secondneigh b or in teraction follo ws WignerDyson statistics, P ( E ) = E 2 e E 2 = 4 2 : (6) The crosso v er describ es the region of c hange from P oisson to WignerDyson statistics. Our understanding of this region is limited, but it is lik ely that the lev elspacing distributions p ossesses only one most probable lev el. Therefore the unkno wn distributions within this crosso v er region eac h p ossess only one p eak. W e hop e that the p eak of the lev elspacing distributions in the crosso v er region will giv e us insigh t in to the c hange of conductances 3 F or example as a function of magnetic rux. 5
PAGE 17
and its nitesize scaling. Since the functional forms of these distributions are not kno wn, w e m ust resort to statistical metho ds suc h as parametric and nonparametric estimators. 1.3 Researc h Plan In order to test the eciencies [12 ] of Unimodal3 and Unimodal2 in determining the p eak p osition of an unkno wn distribution, w e generate sets of deviates from unimo dal distributions on [0 ; 1 ). Here w e will emplo y a randomn um b er generator to pro duce deviates from kno wn probabilit y distributions W eibull, gamma, and b eta. 4 These articial data sets will range from a small n um b er of deviates p er data set to a large n um b er, all ha ving m ultiple realizations so that w e ma y analyze the statistics of the estimators. Since scien tists are unable to gather an innite n um b er of measuremen ts in an actual ph ysical system, w e are in terested in the eciencies of our estimators on data sets that ha v e a small n um b er of data p oin ts. W e will compare the relativ e eciencies of our estimators to those of other estimators suc h as k ernel smo othing, p olynomial smo othing, histograms, and the un biased estimators of the distributions. W e will v ary the parameters of eac h generated distribution to ha v e a range of sk ewness and kurtosis. 5 By c hanging the parameters, w e will pro duce underlying distributions ranging from those with p eaks near the y axis to others that are ratter. 4 Deviate is a sim ulated measuremen t or datum. 5 Refer to App endix A for examples of distributions ha ving dieren t sk ewness. 6
PAGE 18
CHAPTER 2 P ARAMETRIC AND NONP ARAMETRIC ESTIMA TORS Tw o t yp es of estimators w ere p oin ted out in the in tro duction, parametric and nonparametric. This c hapter fo cuses on describing traditional parametric and nonparametric estimators that are used in analyzing data. Some of these estimators ha v e b een mo died to impro v e their abilit y to nd the p eak of a distribution. 2.1 Kernel Densit y Estimator The k ernel densit y estimator is a nonparametric approac h to estimating a probabilit y distribution. Because there is no imp osed parametric mo del, k ernel smo othing allo ws the data to sp eak for themselv es. The basic principles of the k ernel estimator w ere indep enden tly in tro duced b y Fix and Ho dges (1951) and Ak aik e (1954) [13 ]. If w e ha v e a set of data, x i from a distribution P ( x ), then the k ernel estimator of that distribution is giv en b y F ( x ; h ) = ( nh ) 1 n X i =1 K x x i h (7) where h is the parameter kno wn as the bandwith and K is the k ernel function satisfying R K ( x ) dx = 1. W e can condense this form ulation b y in tro ducing a rescaling notation K h ( u ) = h 1 K ( u=h ) ; (8) whic h allo ws us to write F ( x ; h ) = ( n ) 1 n X i =1 K h ( x x i ) : (9) 7
PAGE 19
Because the k ernel function K is c hosen to b e a symmetric unimo dal probabilit y densit y function, then F ( x ; h ) is also a densit y function. W e ha v e c hosen the k ernel function, K to b e a parab ola. The cen tral parameter in k ernel smo othing is the bandwidth h of the estimator. The bandwidth represen ts the base size of the parab ola, or k ernel, that is cen tered on a datum. The v alue of the function F ( x ; h ) at eac h p oin t is simply the a v erage of the k ernel functions at eac h p oin t. The bandwidth determines the degree of smo othing. F or example, t w o smo othing extremes are a consequence of v arying the bandwidth size. If the bandwidth is set to a minimal v alue, undersmo othing of the underlying data is the result. As seen in Fig. 4, ob vious problems of undersmo othing include the presence of m ultiple p eaks and spurious structure in the b o dy of the distribution. T aking this problem to the extreme results in delta functions cen tered at eac h datum. Con trary to undersmo othing is the problem of o v ersmo othing. Ov ersmo othing results when the bandwidth is set to a large v alue. As seen in Fig. 4, problems p ertaining to o v ersmo othed data sets include the shifting of the estimated mo de p osition and the rattening of features. Figure 4. Example of o v er and undersmo othing b y k ernel estimation on a 1000deviate sample dra wn from a b eta distribution ha ving p = 1 : 8 and q = 5. Ra w data from the b eta distribution are plotted as crosses on the x axis. The dotted line represen ts the underlying b eta distribution. T o the left is an example of a distribution that has b een o v ersmo othed. It is noticeable that the smo othed p eak p osition has shifted relativ e the the true p eak p osition. T o the righ t is an example of undersmo othing. 8
PAGE 20
Heuristically it app ears that the most ecien t bandwidth for the purp ose of lo cating the p osition of a single p eak is one that lies just o v er the threshold for undersmo othing. F or example, Fig. 4 sho ws that o v ersmo othing can result in the shifting of the smo othed p eak p osition. Therefore w e wish to nd an h that is close to, but do es not result in undersmo othing. Ho w ev er undersmo othing can result in m ultiple p eaks, but if these p eaks are small compared to a large dominating p eak, then these small p eaks are not lik ely to b e candidates for the p eak of the distribution. Therefore w e can ignore small p eaks that app ear in undersmo othing. No w w e can nd an h that pro duces a distribution that is not o v ersmo othed and not quite undersmo othed. Fig. 5 sho ws an example. Figure 5. Example of a k ernelsmo othed b eta distribution p = 1 : 8 and q = 5. The ra w data are plotted on the x axis. P 1 is a mo de candidate and P 2 is not. If P 2 is less than half the size of P 1, then P 2 is not a mo de candidate. All p eaks that are less than half the size of the most dominan t p eak are ignored. Here the h parameter is set to b e as small as p ossible, consisten t with there b eing only a single p eak candidate. 9
PAGE 21
2.2 Optimal Kernel Densit y Estimator The Optimal Kernel Smo othing (OKS) metho d is an alternativ e to traditional k ernel smo othing. Unlik e traditional k ernel smo othing, where the user m ust test for the b est bandwidth, OKS determines h that b est satises the criterion for a unimo dal distribution. Sp ecically t w o bandwidths m ust b e determined in order to brac k et the most ecien t bandwidth. The rst brac k eting bandwidth, h 1 is set so that the resulting smo othed distribution will b e o v ersmo othed. This can b e ac hiev ed b y setting h 1 to the standard deviation of the distribution or some m ultiple. The second brac k eting bandwidth, h 2 is set so that the resulting smo othed distribution is undersmo othed. Here h 2 is a small fraction of the standard deviation. Once these t w o bandwidths ha v e b een found, the optimization routine then searc hes for a bandwidth h o b et w een h 1 and h 2 that is as small as p ossible consisten t with unimo dalit y The optimal bandwidth h o is the one that lies just under the threshold for undersmo othing. Again undersmo othing can result in m ultiple p eaks; therefore OKS also compares the p eaks using a heigh t restriction for mo de candidates to help determine the optimal h o Here the heigh t restriction is 1 = 2 the heigh t of the highest p eak. An y p eak that falls b ello w this heigh t is not a mo de candidate. This metho d determines h o for ev ery data set. This diers from the k ernel estimation in section 2.1 where w e hold h constan t o v er all 500 realizations, and run sev eral trials v arying h to determine h o W e will from no w on refer to the k ernel estimator in 2.1 as Fixed Kernel Smo othing (FKS). 2.3 Adaptiv e Histogram Another nonparametric estimator is the adaptiv e histogram (AHIST) estimator. This metho d also has one parameter and do es not rely on a parametric mo del. The cen tral parameter in AHIST estimation is the n um b er of data p oin ts in a bin. Unlik e traditional histogram bins, eac h adaptiv e bin holds the same n um b er of data p oin ts. The width of the bin is increased to encompass the giv en n um b er of p oin ts. The area under the bin 10
PAGE 22
equals the n um b er of p oin ts and is held constan t, therefore the heigh t v aries in v ersely as the width. As sho wn in Fig. 6, data p oin ts that are close pro duce a taller bin than those that are spread farther apart. One do es not wish to set the bin size to b e to o large b ecause this w ould decrease the n um b er of bins b eing used. Similar to the problem of o v ersmo othing, a small n um b er of bins rattens out an y features that migh t b e presen t in the underlying distribution. Also one do es not wish to use a bin size that is to o small. This w ould result in the undersmo othing of the distribution. The optimal bin size is one that puts the smallest n um b er of data in eac h bin but do es not undersmo oth the distribution. Fig. 7 sho ws an example. Figure 6. The graph on the left represen ts a b eta distribution that has b een o v ersmo othed b y AHIST. The graph on the righ t represen ts the same b eta distribution undersmo othed. The underlying distribution is plotted as the dotdashed line. Both graphs are dra wn from a sample of 1000deviates. 2.4 Optimal Cheb yshev P olynomial Fit Optimal Cheb yshev P olynomial Fitting (OCPF) is the next nonparametric estimator that w e will in tro duce. In this section w e will explain ho w w e use OCPF to smo oth the ra w cum ulativ e distribution. W e are calling this metho d optimal b ecause w e ha v e included a t w otest pro cedure for tting a curv e to the ra w cum ulativ e distribution. These tests insure that the t results in a singlemo de probabilit y distribution. 11
PAGE 23
Figure 7. The histogram ab o v e is of a b eta distribution with an optimal bin parameter of .1, whic h is 10 p ercen t of the data set. The underlying distribution is plotted as the dotdashed line. If w e ha v e a ra w cum ulativ e distribution, C ( x ), formed from a set of random deviates from a unkno wn distribution, OCPF rst determines the domain of the data p oin ts, 1 then calculates the Cheb yshev co ecien ts using ether a cubic or linear in terp olation. Here w e ha v e decided to use a linear in terp olation, due to smo othing anomalies pro duced when using a cubic in terp olation. 2 The resulting n um b er of co ecien ts equals n + 1, where n is the order of the t determined b y OCPF. OCPF starts with a n = 4 order t and then increases n b y 1, dep ending on the t w o tests. Once all of the co ecien ts ha v e b een calculated, OCPF can appro ximate C ( x ) and then b y taking the deriv ativ e, the underlying probabilit y distribution. In order to tell whether the nal t to the cum ulativ e distribution has b een optimized, OCPF incorp orates t w o tests on its deriv ativ e. As in the previous estimation metho ds there are t w o ma jor problems, o v er and undersmo othing. A large n can result in undersmo othing while a small n can result in o v ersmo othing. What w e w ould lik e is the highestorder n that do es not result in under1 App endix C describ es a c hange of v ariables necessary to use Cheb yshev p olynomials. 2 This smo othing anomaly in v olv es a h ump at the b eginning of the cum ulativ e distribution t, that is most lik ely not part of the underlying distribution. 12
PAGE 24
smo othing. Therefore the rst test after the appro ximation of C ( x ) is one that determines if the 1 st deriv ativ e of the t is unimo dal. As sho wn in the k ernel estimation, undersmo othing can result in m ultiple p eaks, therefore OCPF incorp orates the same p eakcandidate tests as OKS. The second test that OCPF p erforms is a test on the b eginning of the 1 st deriv ativ e of the t. Because w e understand that the distributions w e are trying to estimate are unimo dal, the 1 st deriv ativ e of the t close to the y axis should not ha v e a negativ e slop e. 3 If an n th order t pro duces a curv e that has a negativ e slop e at the start of a estimated probabilit y distribution, then OPCF thro ws a w a y the t and pro cceeds to a higher order estimation. What is considered to b e an optimal t is one that pro duces the highest order estimation, do es not ha v e a negativ e slop e at the b eginning of its 1 st deriv ativ e, and do es not undersmo oth the distribution. Fig. 8 and 9 sho w examples of a go o d and a bad t. 2.5 Un biased P arametric Estimator Eac h of our three test distributions p ossesses t w o parameters that determine its shap e. 4 The follo wing estimators are parametric un biased estimators of the distributions. These estimators determine the distributions' parameters in order to lo cate the p eak p osition for W eibull, gamma, and b eta distributions. Since an y distribution or p opulation sample on [0 ; 1 ), can b e scaled to ha v e unit mean, w e require that t w o of the test distributions ha v e a mean equal to one. Therefore when w e pro duce a random data set from the W eibull and gamma distributions, w e ha v e built in the constrain t that the mean ha v e a v alue of unit y Because w e ha v e a nite n um b er of data, the mean of the sample is not exactly 1, but is close. W e will use the fact that the mean is equal to one in our W eibull and gamma parametric estimators. The reason w e are calling these parametric estimators un biased is that the mean error in estimating the mo de p osition has a v alue v ery close to zero. Lunneb org giv es the bias 3 This will fail on a P oisson distribution. 4 Refer to App endix B for distribution parameters. 13
PAGE 25
of an estimator as B ( t j ) = 1 N N X i =1 ( t i ) (10) where in our case t is the estimator's predicted p eak p osition and is the true p eak p osition [4 ]. W e do not nd the p eak estimations to b e more hea vily w eigh ted to either side of the true p eak p osition. Ho w ev er this is the case only for an estimator that acts on its o wn distribution. F or example, if w e ha v e an W eibull estimator that can appro ximate the parameters from a unkno wn W eibull distribution, then the W eibull estimator will b e un biased in estimating the p eak p osition. Ho w ev er this ma y not b e the case for the W eibull estimator on a b eta distribution. The gamma distribution is dened as G ( x; a; b ) = 1 b a ( a ) x ( a 1) exp[ ( x b )] (11) where a and b are the parameters to b e estimated and are greater than zero, and ( a ) is the Gamma function. 5 The mean of a gamma distribution as a function of its parameters is = ab: (12) Since w e are dealing with a = 1, then a = 1 =b and therefore w e only ha v e one parameter left to estimate. Johnson and Kotz [1 ] giv e an un biased estimator for the gamma distribution log ( =G ) = log ( a ) ( a ) (13) where G is the geometric mean dened as G = N Y i =1 x 1 = N i (14) 5 Refer to App endix B.1 for denition of gamma function 14
PAGE 26
and ( a ) is the digamma function dened as ( a ) = 0 ( a ) ( a ) : (15) Since the geometric mean, G is liable to blo w up on an y signican t data set, (13) can b e written as log ( ) 1 N N X i =1 log ( x i ) = log ( a ) ( a ) : (16) In order to estimate a w e m ust rst simplify (16) : log ( a ) + ( a ) = 0 (17) where is a constan t n um b er for eac h data set and is the com bination of the mean and the geometric mean. F rom (17) the problem is reduced to a ro otnding routine. In order to nd the ro ot of (17) w e m ust rst nd t w o a v alues. The rst a 1 m ust result in a negativ e v alue for the lefthand side of (17) and the second a 2 m ust result in a p ositiv e v alue. Once a 1 and a 2 ha v e b een found, the gamma estimator incorp orates a ro otnding routine that con v erges on the b est a v alue b et w een a 1 and a 2 After the parameter a has b een determined, w e set b = 1 =a F rom (11) w e estimate the p eak p osition X M = b ( a 1) = 1 b: (18) Fig. 10 sho ws a histogram of the absolute errors for the gamma estimator in estimating the p eak of an almostzerosk ewness gamma distribution with a = 100 and b = 0 : 01. Fig. 11 sho ws a histogram of the absolute errors for the gamma estimator on estimating the p eak of a gamma distribution that has a higher sk ewness with a = 2 and b = 0 : 5. The W eibull distribution is dened as W ( x; a; c ) = ca 1 ( x a ) ( c 1) exp[ ( x a ) 2 ] (19) 15
PAGE 27
where a and c are the parameters to b e estimated. The mean of the W eibull distribution is = a (1 + c 1 ) (20) where ( c ) is the Gamma function. Johnson and Kotz [1 ] giv e un biased estimators for the W eibull distribution a = exp 1 N N X i =1 log ( x i ) + r c # (21) and 1 c = p 6 vuut 1 N 1 N X i =1 (log ( x i ) ) 2 (22) where r in (21) is Euler's constan t dened b y r = lim n !1 n X i =1 1 i ln( n ) 0 : 577215664901 532 86 06 06 51 209 00 8 (23) and N is the n um b er of data. As can b e seen from (21) in order to estimate a w e need again the log ( G ), where G is the geometric mean. If w e ha v e the sp ecial case when = 1, then w e ma y use (22) to solv e for c then substitute this in to the new equation, a = c= (1 =c ), to nd a With these t w o parameters w e can estimate the p eak p osition as X M = ( c 1 c ) (1 =c ) a: (24) The b eta distribution is B ( x; p; q ) = x p 1 (1 x ) q 1 ( p; q ) (25) dened on the in terv al [0 ; 1] ; where p and q are the parameters to b e estimated. Since it is dened on a xed nite in terv al, w e cannot rescale it to unit mean as w e did with W eibull 16
PAGE 28
and gamma. Johnson and Kotz [1] giv e un biased estimators for b eta [1 ] ( p ) ( p + q ) = 1 N X j log ( x j ) (26) ( q ) ( p + q ) = 1 N X j log (1 x j ) : (27) In order to calculate p or q w e can either incorp orate a t w odimensional ro otnding routine or replace one of the parameters. Here w e decide to replace the q parameter b y using the mean for a b eta distribution as a function of its parameters giv en as = p q + p : (28) No w q = p (1 ) = and (27) can b e written as ( p ) ( p ) = 1 N X j log ( x j ) : (29) In order to estimate p w e m ust simplify (29) ( p ) ( p ) = 0 (30) where is the log ( G ). The problem is reduced again to a ro otnding routine for p This routine is the same pro cedure as the gamma ro ot nding for p When the p v alue has b een estimated w e can use (28) to calculate q and then estimate the p eak p osition as X M = p 1 q + p 2 : (31) 17
PAGE 29
0.5 1 1.5 2 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Figure 8. Cheb yshev p olynomial t to a W eibull distribution with a = 1 : 12838 and c = 2 ha ving 1000 deviates. The order of the t, n = 12, is not optimal b ecause of the negativ e slop e at the b eginning of the distribution. The true distribution is sho wn in the top righ t. 18
PAGE 30
0.5 1 1.5 2 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Figure 9. OCPF to a W eibull distribution with a = 1 : 12838 and c = 2 ha ving 1000 deviates. Here the order of the t is n = 4, and there are no negativ eslop e problems at the start of the distribution and no m ultiple p eaks. The true distribution is sho wn in the top righ t. 19
PAGE 31
80 100 120 140 0.01 0.02 0.03 0.04 Figure 10. The histogram bin sho ws the absolute error for the gamma estimator on a sample of 1000 deviates from a gamma distribution with p = 100 and b = 0 : 01. The underlying distribution is represen ted in the top left corner. As can b e seen, the absolute errors are cen tered ab out the 0 p oin t on the x axis, indicating a lac k of bias. The small spread ab out zero sho ws furthermore that the estimation is ecien t. 20
PAGE 32
1 2 3 4 5 6 0.05 0.1 0.15 0.2 0.25 0.3 0.35 Figure 11. The histogram sho ws the absolute error for the gamma estimator on a sample of 1000 deviates from a gamma distribution with p = 2 and b = 0 : 5 While the estimator again app ears (n umerically) to b e un biased, it is also less ecien t. The underlying distribution is represen ted in the top left corner. Here the underlying distribution has a higher sk ewness than in Fig. 10. 21
PAGE 33
CHAPTER 3 SEMIP ARAMETRIC ESTIMA TORS In the previous c hapter, w e in tro duced parametric and nonparametric estimators. W e will in tro duce t w o estimation algorithms that use partial kno wledge of the geometry of the underlying distribution. The t w o estimators, Unimodal3 and Unimodal2 use parametric tting and exhaustiv e searc h routines to lo cate the domain around the p eak of a unimo dal probabilit y distribution and then estimate the p eak lo cation. W e w ould lik e to note here that all of the deriv ativ es are estimated using nite dierences as explained in app endix D.3.1 Unimodal3 Unimodal3 is the rst estimator w e will discuss. As men tioned b efore, this estimator uses a t w ostage smo othing tec hnique in order to lo cate the p eak p osition of a probabilit y distribution. The fundamen tal idea in this estimator is imp osing unimo dalit y Unimodal3 ac hiev es this with a lo cal constrain t in its secondstage estimation. The second stageestimation is made b y tting a curv e to a domain in the ra w cum ulativ e distribution, using constrained leastsquares, the constrain t b eing that the third deriv ativ e to the t ha v e only negativ e v alues. 1 The rst stage is only a necessary adjunct. Since tting all of the data will most lik ely not result in negativ e v alues in the third deriv ativ e, w e w ould lik e rst to estimate the domain of the data with only negativ e curv ature and th us con taining the p eak. After w e ha v e collected the data in to a ra w cum ulativ e distribution, w e can apply an OCPF routine to t the data. By observing where the third deriv ativ e to the t of the ra w 1 Refer to app endix D for nite dierences related to the 3 r d deriv ativ e. 22
PAGE 34
cum ulativ e distribution drops in to the negativ e, w e lo cate the domain in the probabilit y distribution that has only negativ e curv ature. Fig. 12 sho ws an example of a cum ulativ e distribution and its 1 st and 3 r d deriv ativ es. Figure 12. The cum ulativ e distribution is represen ted b y the solid line, the 1 st deriv ativ e b y the dashed line, and the 3 r d deriv ativ e b y the dotdashed line. W e need to nd the domain along the x axis in whic h the 3 r d deriv ativ e drops in to the negativ e. The plots are an ideal case, in order to illustrate the cum ulativ e distribution and its 1 st and 3 r d deriv ativ es. The 1 r d and 3 st deriv ativ es ha v e b een plotted in arbitrary units. W e are lo oking for the turning p oin ts whic h brac k et the mo de of the probabilit y distribution. A t the turning p oin ts to the probabilit y distribution, the third deriv ativ e of the t to the cum ulativ e crosses the x axis. Fig. 13 sho ws the output of the rststage estimation on the ra w cum ulativ e distribution. W e will call the output of the rststage estimator for Unimodal3 the ltered ra w cum ulativ e distribution. The second stageestimator of Unimodal3 mak es use of the w ell kno wn metho d of constrained leastsquares tting. Here Unimodal3 p erforms a constrained leastsquares t to the ltered ra w cum ulativ e distribution, where the 3 r d deriv ativ e is constrianed to p ossess only negativ e v alues. Fig. 14 sho ws the leastsquares t to the ltered ra w cum ulativ e distribution and the deriv ativ e of that t. As can b e noticed in Fig. 14, Unimodal3 is not useful in determining the shap e of the probabilit y distribution ab out 23
PAGE 35
Figure 13. The graph to the left is the ra w cum ulativ e distribution tak en from a b eta distribution, p = 1 : 8 and q = 1 : 8, ha ving 1000 deviates. The graph to the righ t is the ra w cum ulativ e distribution after the rststage estimator has estimated the range of negativ e curv ature in the probabilit y distribution. The true mo de p osition for this distribution is at x = : 5. The mo de p osition falls within the ltered domain. the p eak. This is b ecause the deriv ativ e of the constrained leastsquares t empirically is piecewiselinear. Increasing the size of the data set b eing used do es little to impro v e the smo othness of the output of Unimodal3 3.2 Unimodal2 Similar to Unimodal3 Unimodal2 can also use a t w ostage smo othing tec hnique to determine the p osition of the p eak for an unkno wn distribution. The second stage of the estimator p erforms an exhaustiv e or a binary searc h for the p eak p osition. The secondstage estimation assumes a p eak p osition, then ts the ra w cum ulativ e data sub ject to the constrain t that the second deriv ativ e b e p ositiv e to the left and negativ e to the righ t of the assumed p eak p osition. The p eak p osition resulting in the smallest residual is the estimated mo de. Fig. 15 sho ws an example of a cum ulativ e distribution and its 1 st and 2 nd deriv ativ es. If w e c ho ose to use a rststage estimation, it is to decrease the run time of the second stage. There are three dieren t metho ds w e can use to searc h for the p eak p osition, the c hoice en tailing a tradeo b et w een accuracy and run time. W e can p erform an exhaustiv e searc h 24
PAGE 36
Figure 14. The graph to the left is the constrained leastsquares t to the ltered ra w cum ulativ e distribution of Fig.13. The graph to the righ t is the deriv ativ e of the leastsquares t. The mo de p osition for this distribution is at x = : 5, and Unimodal3 estimated the p eak to b e at x = 0 : 56034. for the p eak p osition, whic h tests ev ery datum as a p eak estimation, or w e can p erform a binary searc h, whic h starts at b oth ends of the domain and con v erges on the p eak estimation. The binary is more time ecien t than the exhaustiv e searc h; ho w ev er it could con v erge on a lo cal (rather than global) minim um in the ob jectiv e function. If w e c ho ose not to decimate the data (rststage estimation), exhaustiv e searc h tak es prohibitiv ely long. The binary searc h on undecimated data for a range of n um b er of deviates from 400{1000 a v erages 14 seconds p er realization. F or 500 realizations, it w ould tak e 14 hours to run a full set of samples. 2 The exhaustiv e searc h on decimated data a v eraging 300 deviates a v erages 4 seconds p er realization. F or 500 realizations this w ould only tak e 4 hours to run on a full set of samples. Figs. 17 and 18 sho w that for 200deviate samples, the RMS errors are indistinguishable (within error bars) b et w een the binary searc h on undecimated data and the exhaustiv e searc h on decimated data. Therefore using binary searc h on decimated data pro v es to b e the most ecien t metho d for run time with no tradeo in accuracy Fig. 19 sho ws that w e can impro v e the run time eciency without a tradeo in accuracy b y using a binary searc h on the decimated data. Using a binary searc h on decimated data pro v es to b e more time ecien t than exhaustiv e on decimated. 2 This run time is on a P en tiumIV pro cessor 2.4GHz. 25
PAGE 37
In theory binary can nev er b e more ecien t than exhaustiv e. Ho w ev er Fig. 19 sho ws sligh tly smaller errors for the binary searc h. This ma y b e due to the bias of Unimodal2 None the less,to reduce the c hance of estimations falling in to lo cal minima, w e c ho ose to run an exhaustiv e searc h on decimated data. The rststage estimation for Unimodal2 decimates the n um b er of data p oin ts, using a distance scale to compare eac h datum to its neigh b ors. F or example, when the measuremen ts are group ed closely together, it is lik ely that this range of grouping can hold a p ossible p eak. Because w e are dealing with unimo dal distributions, there should b e only one ma jor domain of group ed measuremen ts. Here w e set a length scale that is a fraction of the a v erage dierence b et w een measuremen ts. Unimodal2 compares ev ery p oin t with its neigh b ors using this length scale as a limit. If a datum falls within the limit of its neigh b ors, then that datum is sa v ed. Ho w ev er if that datum is not within the limit of its neigh b ors, then that datum is remo v ed. This pro cess remo v es data that fall out on the tails of a distribution, where a p eak is less lik ely Fig. 16 sho ws examples b efore and after the decimation. As can b e seen in Fig. 16, some of the data in the middle of the distribution ha v e b een ltered out. A t rst, this app ears to b e a problem. W e decided to compare Unimodal2 's rststage estimator, whic h remo v es data from the edges and the middle, against Unimodal3 's rststage estimator, whic h remo v es data from only the edges, and see ho w the exhaustiv e searc h for Unimodal2 compares. W e pro duced t w o sets of ninet y realizations ha ving n um b ers of deviates ranging from 400 to 1000 and ran the rststage estimations for Unimodal2 and Unimodal3 on eac h set. The size of the output les from Unimodal3 's rst stage ranged from 50 to 800 deviates, while the a v erage size of Unimodal2 's rststage output w as 250 deviates. Therefore it generally tak es longer to run Unimodal2 's exhaustiv e searc h on the output of Unimodal3 's rst stage than on Unimodal2 's rst stage. W e nd that the time it tak es to run Unimodal2 's exhaustiv e searc h on a data set with X n um b er of deviates scales as X 3 : 3 So for a data set with 1000 deviates, the run tak es 20 times as long as a data set 26
PAGE 38
Figure 15. The smo othed cum ulativ e distribution is represen ted b y the solid line, the 1 st deriv ativ e estimate b y the dotdashed line, and the 2 nd deriv ativ e b y the dashed line. The plots are an ideal case in order to illustrate the p oin t where the 2 nd deriv ativ e crosses the x axis. This p oin t represen ts the p eak p osition in the 1 st deriv ativ e. The 1 st and 2 nd deriv ativ es are plotted in arbitrary units. with 400 deviates. Ho w ev er, the o v erall eciency in calculating the p eak p osition do es not dier b y m uc h when w e use Unimodal2 's or Unimodal3 's rststage estimation. Therefore to run in a reasonable time, w e use Unimo dal2's rststage estimator. After w e ha v e assumed a mo de p osition and estimated the cum ulativ e probabilit y distribution sub ject to the secondderiv ativ e constrain t, w e compare this estimation to the ra w cum ulativ e distribution and nd the sum of square dierences. When w e ha v e a p eak candidate that pro duces the smallest sum of square dierences, w e ha v e found the b est estimate for the p eak p osition. 27
PAGE 39
Figure 16. The graph to the left sho ws the ra w cum ulativ e distribution for a gamma distribution with p = 100 and b = 0 : 01 tak en from a sample of 1000 deviates. The graph to the righ t is after Unimo dal2's rststage estimation ha ving 311 deviates. As can b e seen, the tails of the distribution ha v e b een thinned due to the unlik eliho o d of their con taining a p eak candidate. Figure 17. The dashed lines are Unimodal2 's exhaustiv e searc h on decimated data, and the solid lines are the binary searc h on undecimated data. The graph to the left is of 200deviate samples dra wn from a family of gamma distributions. There are 500 realizations p er p oin t. The graph to the righ t is of 200deviate samples dra wn from a family of W eibull distributions. There are 500 realizations p er p oin t. 28
PAGE 40
Figure 18. The dashed lines are Unimodal2 's exhaustiv e searc h on decimated data, and the solid lines are the binary searc h on undecimated data. The graph to the left is of 200deviate samples dra wn from a family of sharpp eak ed b eta distributions. There are 500 realizations p er p oin t. The graph to the righ t is of 200deviate samples dra wn from a family of ratp eak ed b eta distributions. There are 500 realizations p er p oin t. Figure 19. The dashed lines are Unimodal2 's exhaustiv e searc h on decimated data, and the solid lines are the binary searc h on decimated data. The graph to the left is of 200deviate samples dra wn from a family of gamma distributions. There are 500 realizations p er p oin t. The graph to the righ t is of 1000deviate samples dra wn from a family of gamma distributions. There are 500 realizations p er p oin t. 29
PAGE 41
CHAPTER 4 MODEESTIMA TION EFFICIENCY In this c hapter w e will discuss ho w w e determine the relativ e eciency of our estimators in lo cating the true p eak p osition of an underlying probabilit y distribution. 4.1 Setting a Scale In this section w e will explain wh y w e decide to use the standard deviation as a comparison scale. It will b e sho wn that a set of realizations that ha v e broad rat p eaks can pro duce higher errors than realizations with sharp er, p oin tier p eaks. Therefore, w e can use the standard deviation of the p opulation as a scale for comparing the estimators' eciencies. F or example, if w e ha v e t w o dieren tly shap ed distributions, as in g.20, without using the standard deviations of the distributions as a scale, all estimators will p erform b etter on the sharp p eak than on the broad rat p eak. Ho w ev er, if w e divide the error b y the standard deviation of the p opulation, w e can b etter compare the dierences in the estimations for eac h p eak. Therefore w e set the comparison function of ho w w ell the estimator estimates the true p eak p osition for sev eral realizations to b e ( ) = R M S S T D ; (32) where c haracterizes the distribution, R M S is the ro ot mean square error, and S T D is the standard deviation of W e w ould exp ect the R M S and S T D to decrease with a sharp er skinner p eak and increase with a rat broader one. 30
PAGE 42
Figure 20. Without using the standard deviation of the p opulation as a comparison scale, all estimators w ould p erform w ell on this sharp p eak and p o orly on the rat p eak. 4.2 Mirror Distributions In this section w e will explain wh y w e c ho ose to divide the RMS b y STD. Ev en though w e are in terested in ho w w ell our h ybrid estimators estimate the true p eak p osition of a highly sk ew ed distribution, w e m ust consider data sets that ha v e negativ e and p ositiv e sk ewnesses on the same comparison scale. 1 That is, w e cannot just compare the estimators' eciency in lo cating p eak p ositions that are only close to the y axis. F or example, a distribution with a sk ewness v alue of S can ha v e a mirror coun terpart with sk ewness S Fig. 21 sho ws a b eta distribution with sk ewness v alues that are absolutely equal but ha v e dieren t signs. These distributions also ha v e iden tical kurtosis and v ariance v alues. The only prop erties that will dier for these distributions are the mean and the signs of the sk ewness and higher o ddorder cum ulan ts. 2 Our estimators should pro duce similar errors for a distribution that has S or S Therefore w e cannot consider the p ercen tage error in estimating a p eak p osition as a w a y of comparing the eciencies of our estimators. F or example if w e ha v e t w o distributions as in Fig. 21, and w e ha v e the same absolute error estimation on b oth, 1 P ositiv e sk ewness is describ ed b y the p eak p osition b eing shifted asymmetrically to w ard the y axis. 2 Refer to app endix A for denition of cum ulan ts. 31
PAGE 43
0.2 0.4 0.6 0.8 1 1 2 3 4 0.2 0.4 0.6 0.8 1 1 2 3 4 Figure 21. This gure sho ws a mirror image of a b eta distribution cen tered ab out x = : 5. The distribution to the left has parameters p = 1 : 8 and q = 10. The distribution to the righ t has p = 10 and q = 1 : 8. w e cannot divide the errors b y the p eak p ositions and still ha v e the same fractional error. Because of this, w e m ust nd another w a y of determining the estimators' eciencies. Since the t w o distributions in the example ha v e the same standard deviation, w e can divide the absolute error b y the standard deviation and pro duce the same fractional result. Therefore our estimators will pro duce the same errors for mirror distributions. 4.3 STD of the STD as an Error assessmen t Because w e are going to b e plotting the RMS errors o v er the STD, w e need to kno w ho w signican t our comparisons are. If w e assume that the distribution of mo de estimations o v er man y realizations is Gaussian ab out the true mo de p osition, w e m ust assess the standard deviation of the standard deviation of the mo de estimate. Keeping giv es the 32
PAGE 44
v ariance of the v ariance estimate as 2 k 2 = 2 N + 1 k 2 2 (33) where k 2 is the v ariance of an arbitrary distribution [2 ]. F or our purp ose w e are going to call k 2 the mean square error giv en as k 2 = 1 N X i ( x i ) 2 ; (34) where x i is the mo de estimation for the i th realization and is the true mo de p osition. If w e wish to nd the standard deviation k 2 w e ha v e k 2 = k 2 r 2 N + 1 (35) whic h giv es us the standard deviation as a function k 2 Ho w ev er w e are in terested in p k 2 whic h is the standard deviation of the standard deviation of the mo de estimates. Therefore w e can use the relation giv en b y T a ylor for uncertain t y in p o w er, q j q j = j n j x j x j ; (36) where x is measured with uncertain t y x and is used to calculate the p o w er q = x n [11 ]. The fractional uncertain t y in q is j n j times that in x If w e mak e the relations to (36) q is p k 2 q is p k 2 n is 1 2 x is k 2 and x is k 2 Solving (36) for q and substituting the previous relations giv es p k 2 = s k 2 2( N + 1) : (37) No w w e ha v e an expression for the standard deviation of the standard deviation of the mo de estimate where p k 2 is the RMS error of the mo de estimate for N realizations. 33
PAGE 45
4.4 Choice of The nal decision to b e made on ho w to compare the relativ e eciencies of our estimators is what should the set b e in (32) ? What prop ert y of the distribution should w e compare against? In the quan tumwire lev elspacing distributions, w e are in terested in probabilit y distributions that lie b et w een the P oisson and WignerDyson distributions. Some ma y ha v e high sk ewness. Therefore one of the sets w e will use is the sk ewness of the underlying distribution. Since w e kno w nothing outside of unimo dalit y in these distributions, w e can also use the kurtosis of the underlying distributions as a second set. 34
PAGE 46
CHAPTER 5 RESUL TS F OR DISTRIBUTIONS HA VING UNIT MEAN In this c hapter w e will sho w the results of p eak lo cation for eac h estimator on three dieren t family of distributions. Our data sets are pro duced from the W eibull, b eta, and gamma distributions. W e v aried eac h of the distributions' parameters in order to ac hiev e a range of sk ewness and kurtosis v alues, thereb y c hanging the shap es of the distributions. In order to observ e the statistics of these estimators, w e ran a series of tests on data sets eac h con taining 500 realizations. Eac h realization within a set con tained the same n um b er of deviates, and for dieren t sets, the n um b er of deviates ranged from 100 to 1000. W e will sho w ho w eac h estimator compared in relativ e eciency in lo cating the p eak p ositions of data sets with 1000 and 200 deviates b y plotting the mean of the fractional error, (32) o v er all of the realizations, v ersus the sk ewness of the distributions. Because an y data set on the semiinnite in terv al [0 ; 1 ) can b e scaled so that the mean is equal to one, w e ha v e prepared data sets from W eibull and gamma distributions where the mean is equal to one. W e compare the relativ e eciencies of eigh t dieren t estimators on lo cating the p eak p osition of these distributions. These eigh t estimators are k ernel estimation (FKS), optimal k ernel smo othing (OKS), adaptiv e histogram (AHIST), optimal Cheb yshev p olynomial (OCPF), W eibull, gamma, Unimodal2 and the Unimodal3 estimator. W e compare the fractional error of (32) against the sk ewness of the distributions for a set of gamma and W eibull distributions listed in T ables 1 and 2. A ratter p eak has a negativ e kurtosis and a sharp er p eak a p ositiv e kurtosis. As to b e exp ected, if w e are giv en more random deviates in a data set, the mo de estimations ha v e a lo w er fractional error. The more information w e ha v e, the b etter the estimation. Figs. 22, 23, 24, and 35
PAGE 47
a c S k ew ness K ur tosis W eibull 1 : 11142 3 : 5 0 : 0249969 0 : 286522 1 : 11985 3 0 : 168098 0 : 270617 1 : 12706 2 : 5 0 : 358619 0 : 143214 1 : 12838 2 0 : 631109 0 : 245085 1 : 12777 1 : 95 0 : 665277 0 : 311061 1 : 1245 1 : 8 0 : 778745 0 : 556757 1 : 11536 1 : 6 0 : 961966 1 : 04394 1 : 09719 1 : 4 1 : 19847 1 : 83907 1 : 09053 1 : 35 1 : 26919 2 : 11427 1 : 08275 1 : 3 1 : 34593 2 : 4322 1 : 07747 1 : 27 1 : 39518 2 : 6475 1 : 07367 1 : 25 1 : 42955 2 : 80217 1 : 06309 1 : 2 1 : 5211 3 : 2356 1 : 05075 1 : 15 1 : 62209 3 : 74807 1 : 03636 1 : 1 1 : 73397 4 : 35981 T able 1. Fifteen dieren t W eibull distribution, with unit mean. W ( x; a; c ) = ca 1 ( x a ) ( c 1) exp[ ( x a ) 2 ]. p b S k ew ness K ur tosis Gamma 100 0 : 01 0 : 2 0 : 06 20 0 : 05 0 : 447214 0 : 3 10 0 : 1 0 : 632456 0 : 6 4 0 : 25 1 1 : 5 2 0 : 5 1 : 41421 3 1 : 1 0 : 909091 1 : 90693 5 : 45455 T able 2. Six dieren t gamma distributions with unit mean. G ( x; a; b ) = 1 b a ( a ) x ( a 1) exp[ ( x b )]. 36
PAGE 48
25 sho w the a v eraged results o v er 500 realizations of all of the estimators for gamma and W eibull data sets con taining samples of 1000 and 200 deviates. Figure 22. Results from a family of gamma distributions for realizations dra wing 1000 deviates eac h. See also T able 2. F or UNI2 w e ha v e used the exhaustiv e searc h metho d for decimated data. 37
PAGE 49
5.1 W eibull and Gamma Estimators In b oth the 1000and 200deviate samples, the gamma and W eibull estimators pro duced the smallest fractional errors o v er all ranges of sk ewness on their resp ectiv e distributions. This is to b e exp ected since these are parametric un biased estimators of their distributions. Ho w ev er, for a small range of sk ewness v alues in the W eibull distributions, the gamma estimator pro duced lo w er fractional errors than W eibull. This w as not exp ected. If w e observ e the W eibull distributions at these sk ewness v alues and then observ e the gamma distributions with the parameters estimated from the gamma estimator on these W eibull distributions, w e can see that the gamma and W eibull distributions are similar. This w ould explain wh y the gamma estimator do es so w ell for the W eibull distributions in this region of sk ewness. Fig. 26 sho ws an example. 5.2 OKS and OCPF OKS and OCPF had the nextlo w est fractional errors o v er all sk ewness v alues. Ho w ev er in the 1000deviate case, OKS's fractional error diered o v er all the sk ewness range b y only four times the a v erage 1 p k 2 while OCPF's fractional error diered b y eigh teen times its a v erage p k 2 As will b e sho wn in the next c hapter, OCPF suers when estimating the p eaks of distributions that ha v e negativ e kurtosis. 2 Ho w ev er, for the p ositiv e kurtosis v alues of the gamma and W eibull distributions, listed in T able 1, OCPF is rat within : 02 p k 2 In the 200deviate case, OKS suers at high sk ewness v alues. The explanation is simple. When OKS tries to estimate a data set that do es not ha v e enough information, the optimal routine m ust k eep increasing the size of the k ernel parameter to o v ercome undersmo othing. This is often the case with a small n um b er of deviates. If the k ernel parameter is set to o high, o v ersmo othing is the result, and therefore there is a shift of the smo othed p eak p osition relativ e to the true p eak p osition. This shift do es not aect 1 p k 2 is represen ted b y the error bars in gs. 22, 23, 24, 25; refer to (37) for denition. 2 Negativ e kurtosis is represen ted b y a ratter p eak. 38
PAGE 50
distributions that are near zero sk ewness as m uc h as those that are highly sk ew ed. Fig. 27 sho ws an example. 5.3 FKS and AHIST Kernel estimation (FKS) 3 and AHIST estimation follo w similar fractional error patterns. Both of these are nonparametric estimators, and from sev eral trials w e can c ho ose the b est parameter v alue for eac h. With FKS, 10 p ercen t of the standard deviation of the sample app ears to b e a go o d o v erall parameter. This p ercen tage of the standard deviation of the sample yielded the lo w est fractional error o v er all ranges of sk ewness. With the AHIST metho d, 10 p ercen t of the p opulation of deviates w as a go o d o v erall parameter v alue. This is a little more in tuitiv e considering that the parameter v alue represen ts the p ercen tage of data in eac h bin. T o o great or to o small of a p ercen tage can result in under or o v ersmo othing. Fig. 28 sho ws these estimates. As to b e exp ected, when there is a smaller n um b er of deviates in eac h realization, the fractional error increases for eac h estimator. W e see the same error pattern in the 1000deviate case as w e do in the 200deviate case, except the error has increased when there are only 200 deviates in eac h realization. 5.4 Unimo dal3 and Unimo dal2 Our t w o estimators, Unimodal2 and Unimodal3 sho w ed similar fractionalerror problems in b oth the W eibull and gamma distributions. W e encoun tered sev eral problems with the Unimodal3 estimator. The rst problem is that Unimo dal3 sho ws biased results for the mo de estimates. Fig. 29 sho ws an example of fourteen Unimodal3 outputs on a gamma distribution with p = 100 and b = 0 : 01. As can b e seen from this gure, Unimodal3 sho ws a p ositiv e bias, meaning that the ma jorit y of the p eak estimates are to the righ t of the true mo de p osition. W e can observ e this p ositiv e bias o v er all 500 realizations through a histogram plot of the absolute error in Fig. 30 3 Fixedwidth k ernel smo othing (FKS), from c hapter 2.2 39
PAGE 51
The second problem with Unimodal3 is sho wn in the estimations on a gamma distribution with lo w sk ewness and large n um b er of deviates. Fig. 31 sho ws us a RMS error plot v ersus n um b er of deviates. W e w ould exp ect the RMS error to b e high for a small n um b er of deviates in eac h realization and decrease as the n um b er of deviates increase. Ho w ev er w e do not observ e this eect. Here w e notice a dip in the RMS error at 200 deviates and then a rise in the RMS error as the n um b er of deviates increase. This is not what w e w ould exp ect, since more information should lead to a b etter appro ximation. If w e pro duce another 500 realizations with the same distribution parameters sho wn also in Fig. 31, w e see the same systematic results. This problem could b e the result of the p ositiv e bias of Unimodal3 on the gamma distribution. W e do not notice the same systematic eect, in Fig. 31, when Unimodal3 is estimating a W eibull distribution. Ho w ev er w e do notice a bias in estimating the mo de for distributions with small sk ewness. Fig. 32 sho ws a histogram plot of the absolute error o v er 500 realizations for a W eibull distribution ha ving a = 1 : 11142 and c = 3 : 5 and ho w this bias v aries with sk ewness. F or greater sk ewness, the bias is m uc h lo w er than for the nearzerosk ewness distributions. W e can see from Fig. 32 that the magnitude of this bias is less for the W eibull distribution. Ho w ev er, the fractional error for sk ewness v alues ranging b et w een 0 : 5 r 1 2 are similar on b oth distributions in the 1000and 200deviate cases. Again this is most lik ely due to the similarities in the shap es of the distributions within this sk ewness range. In Fig. 33, Unimodal2 sho ws small c hanges in bias o v er a range of sk ewness for b oth gamma and W eibull families. Ho w ev er these c hanges are small. Unimodal2 also sho ws little to no bias for distributions ha ving nearzero sk ewness. This is dieren t from Unimodal3 in Figs. 32 and 30. Figs. 34 and 35 sho w histograms of the bias of Unimodal2 on extremes in sk ewness for the gamma and W eibull families. 40
PAGE 52
Figure 23. Results from a family of gamma distributions for realizations dra wing 200 deviates. See also T able 2. F or UNI2 w e ha v e used the exhaustiv e searc h metho d for decimated data. 41
PAGE 53
Figure 24. Results from a family of W eibull distributions for realizations dra wing 1000 deviates eac h. See also T able 1. F or UNI2 w e ha v e used the exhaustiv e searc h metho d for decimated data. 42
PAGE 54
Figure 25. Results from a family of W eibull distributions for realizations dra wing 200 deviates eac h. See also T able 1. F or UNI2 w e ha v e used the exhaustiv e searc h metho d for decimated data. 43
PAGE 55
Figure 26. The dotted line represen ts a W eibull distribution with a = 1 : 06309 and c = 1 : 2. The solid line is a gamma distribution with p = 1 : 33 and b = : 75. Both are tak en from 1000deviates sample. These parameters for the gamma distribution w ere calculated b y the gamma estimator on the W eibull distribution sho wn b y the solid line. As can b e seen these distributions are close to the same shap e. 44
PAGE 56
Figure 27. In the graph to the left, the dotted line represen ts optimal k ernel estimation on a gamma distribution with p = 20 and b = 0 : 05, and the solid line represen ts the true gamma distribution with the same parameters. This distribution has a lo w sk ewness, and therefore o v ersmo othing rattens the p eak p osition but do es not shift it m uc h. The graph to the righ t is another gamma distribution ha ving, p = 1 : 1, b = 0 : 909091. This distribution has a high sk ewness, and therefore o v ersmo othing rattens the p eak and shifts its p osition. Both are tak en from 200deviate samples. 45
PAGE 57
Figure 28. In the graph to the left, the solid line represen ts a gamma distribution with p = 20 and c = 0 : 05. The dotted line is the k ernel appro ximation with the k ernel parameter equal to 10 p ercen t of the standard deviation of the p opulation. The graph to the righ t is the adaptiv e histogram estimator with the same gamma distribution. Both are tak en from 1000deviate samples. Figure 29. The solid v ertical line represen ts the true mo de p osition for a gamma distribution with p = 100 and b = 0 : 01. The dotted lines are the output of the Unimodal3 distribution on fourteen realization. As can b e seen, the ma jorit y of the p eak estimations fall to the righ t of the true mo de p osition, indicating a p ositiv e bias. All are tak en from 1000deviate samples. 46
PAGE 58
Figure 30. The histogram ab o v e is of the absolute error o v er all 500 realizations of the Unimodal3 estimator on a gamma distribution tak en from 1000deviates samples with p = 100 and b = 0 : 01 and sk ewness of r 1 = 0 : 2. W e can see in the graph to the righ t ho w the bias c hanges with sk ewness. Figure 31. The graph to the left sho ws the Unimodal3 results for the a v erage of 500 realizations on a gamma distribution p = 100 and b = 0 : 01. W e see the same systematic results for the a v erage of another 500 realizations in the graph to the righ t. All realizations tak en from 1000deviate samples. 47
PAGE 59
Figure 32. In the Wiebull distribution w e can see ho w Unimodal3 has a negativ e bias for the mo de estimations on distributions that ha v e near zero sk ewness. The graph to the left is a histogram of the biases for all 500 realizations tak en from a sample of 1000 deviates for a = 1 : 11142 and c = 3 : 5 ha ving a sk ewness of r 1 = 0 : 025. The graph to the righ t demonstrates ho w this bias c hanges with sk ewness. Figure 33. The graph to the left sho ws a plot of the bias of Unimodal2 in estimating the p eak p osition v ersus the sk ewness for the W eibull distribution tak en from 1000deviate samples. The graph to the righ t sho ws the same for a gamma distribution. As can see seen from the graphs, b oth distributions sho w a bias o v er a range of sk ewness v alues; ho w ev er the bias is small. F or Unimodal2 w e ha v e used the exhaustiv e searc h metho d on decimated data. 48
PAGE 60
Figure 34. The histogram to the left sho ws the absolute errors from Unimodal2 of all 500 realizations for a gamma distribution tak en from 1000deviate samples ha ving a zero sk ewness v alue. T o the righ t is the same for a gamma distribution ha ving a large sk ewness of r 1 = 1 : 90693. F or Unimodal2 w e ha v e used the exhaustiv e searc h metho d on decimated data. Figure 35. The histogram to the left sho ws the absolute errors from Unimodal2 of all 500 realizations for a W eibull distribution tak en from 1000deviate samples ha ving a zero sk ewness v alue. T o the righ t is the same for a W eibull distribution ha ving a high sk ewness v alue of r 1 = 1 : 73397. F or Unimodal2 w e ha v e used the exhaustiv e searc h metho d on decimated data. 49
PAGE 61
CHAPTER 6 RESUL TS F OR BET A DISTRIBUTIONS HA VING NEGA TIVE AND POSITIVE KUR TOSIS In this c hapter w e will concen trate on the estimators' relativ e eciency on t w o families of b eta distributions. These distributions are dieren t from the W eibull and gamma distributions in that the means of these distributions are not equal to one, and the domain is constrained to 0 x 1 W e will lo ok at a family of b eta distributions where the p eaks are rat, T able 3, and a family of distributions where the p eaks are sharp, T able 4. 6.1 Beta Estimator The b eta estimator pro duced the lo w est fractional error o v er most ranges of sk ewness. Ho w ev er for the family of rat distributions, in Fig. 38, the b eta estimator p erforms p o orly compared to OKS, and in Fig. 39, the b eta estimator p erforms p o orly compared to AHIST, UNI2, and FKS in estimating the p eaks of nearzero sk ewness. W e found that (30) has m ultiple ro ots, only one of whic h results in a p eak in the in terv al 0 x 1. W e found sev eral realizations resulted in p eak estimates that where outside this in terv al therefore resulting in large fractional errors. Because of these errors w e mo died (30) instead to con v erge on the mo de p osition and not the p v alue. W e used the fact that the domain of the b eta distribution is limited to 0 x 1 in order to set the upp er and lo w er mo de p ositions limits, m 1 and m 2 These limits are similar to the p 1 and p 2 limits in c hapter 2.5. W e can determine these m 's b y calculating the mean, of eac h realization. If for a realization is less than 1 = 2, then m 1 = 0 and m 2 = where is some tolerance. 1 If for a realization is greater than 1 = 2, then m 1 = + and m 2 = 1. Once w e ha v e 1 W e used = 10 8 50
PAGE 62
determined these limits w e can use the relation p = 2 m 1 m (38) substituted in to (30) to giv e log ( 2 m 1 m ) + ( 2 m 1 m ) = 0 (39) again where is a constan t n um b er for eac h data set and is the com bination of the mean and the geometric mean. By using the same ro otnding routine w e can con v erge on the estimated mo de for eac h realization. Ho w ev er, on some realizations w e ran in to problems when nding the upp er and lo w er limits of m If our limits did not pro duce opp osing signs in (39) and < 1 = 2, then w e set the estimated mo de for that realization equal to zero. If our limits did not pro duce opp osing signs in (39) and > 1 = 2, then w e set the estimated mo de for that realization equal to one. By observing the spik es on the edges of the histogram of absolute errors in Fig. 36, w e can see sev eral realizations of the b eta distribution ha ving p = 1 : 1 and q = 1 : 1 pro duced this problem. 6.2 OKS and OCPF OKS and OCPF resulted in the next smallest fractional errors o v er all sk ewnesses. Ho wev er, in the case of large negativ e kurtosis, OCPF pro duces a larger fractional error than Unimodal3 FKS, and AHIST. Fig. 42 sho ws a ratp eak ed distribution, ha ving steep sides, with a sk ewness of zero. OCPF m ust use a highorder p olynomial in order to t the steep sides of this distribution. When the order is increased, the p eak of the distribution is undersmo othed. Therefore OCPF cannot t the sharp sides without undersmo othing the distribution. OCPF do es not suer from this problem when the sk ewness increases b ecause the ratness of the p eak decreases. Fig. 38 and Fig. 40 sho w examples. OKS results in the lo w est fractional error for all of the estimation metho ds at the largest negativ e kurtosis. In Fig. 43, w e can see that OKS do es not suer from the 51
PAGE 63
p q S k ew ness K ur tosis Beta 1 : 1 1 : 1 0 1 : 15385 1 : 1 1 : 2 0 : 0735413 1 : 12549 1 : 1 1 : 3 0 : 140178 1 : 08709 1 : 1 1 : 4 0 : 201008 1 : 04132 1 : 1 1 : 5 0 : 256887 0 : 990119 1 : 1 1 : 6 0 : 308494 0 : 934923 1 : 1 1 : 7 0 : 356378 0 : 876821 1 : 1 1 : 8 0 : 400988 0 : 816641 1 : 1 1 : 9 0 : 442697 0 : 755024 1 : 1 2 0 : 481818 0 : 69247 T able 3. T en dieren t b eta distribution families. B ( x; p; q ) = x p 1 (1 x ) q 1 ( p; q ) p q S k ew ness K ur tosis Beta 1 : 8 1 : 8 0 0 : 909091 1 : 8 1 : 9 0 : 0531526 0 : 893363 1 : 8 2 : 3 0 : 181951 0 : 802405 1 : 8 3 0 : 36578 0 : 594268 1 : 8 6 0 : 773718 0 : 259259 1 : 8 15 1 : 14018 1 : 54851 1 : 8 50 1 : 37243 2 : 6643 T able 4. Sev en dieren t b eta distribution families. B ( x; p; q ) = x p 1 (1 x ) q 1 ( p; q ) 52
PAGE 64
Figure 36. The histogram on the left is of absolute errors in b eta mo de estimations for a b eta distribution with p = 1 : 1 and q = 1 : 1 ha ving 1000 deviates, 500 realizations and sknew ess r 1 = 0. The histogram on the righ t is of the b eta mo de estimations for a b eta distribution with p = 1 : 1 and q = 2 ha ving 1000 deviates, 500 realizations and sk ewness r 1 = : 481818. The spik es on either side of the histograms are due to mo de estimation errors found when estimating ratp eak ed families. undersmo othing problems of OCPF. Ho w ev er, from Fig. 43, OKS do es a p o or job in estimating the steep sides of the ratp eak ed b eta distribution. This is not a problem that aects the p eak estimation. If w e reduce the n um b er of deviates to 200 in eac h realization, OKS's and OCPS's fractional errors increase. Ho w ev er OKS's fractional error in the ratp eak ed families, Fig. 38, increases m uc h more within the range of sk ewness than in the sharp erp eak ed families, Fig. 39. This is b ecause OKS suers from o v ersmo othing when trying to estimate a realization ha ving a small n um b er of deviates. 6.3 FKS and AHIST FKS and AHIST pro duced nearly the same fractional errors for the t w o families of b eta distributions. As to b e exp ected, the fractional errors increased as the n um b er deviates decreased. 53
PAGE 65
Figure 37. The histogram on the left is of absolute errors in b eta mo de estimations for a b eta distribution with p = 1 : 8 and q = 1 : 8 ha ving 1000 deviates, 500 realizations and sknew ess r 1 = 0. The histogram on the righ t is of the b eta mo de estimations for a b eta distribution with p = 1 : 8 and q = 50 ha ving 1000 deviates, 500 realizations and sk ewness r 1 = 1 : 37243. 6.4 Unimodal3 and Unimodal2 Unimodal2 2 and Unimodal3 pro duced the largest fractional errors on b oth families of b eta distributions o v er all ranges of sk ewness. F or the family of distributions in T able 3, Unimodal3 suers from problems in the rststage estimations. Unimodal3 pro duces the largest fractional error for all ranges of sk ewness b ecause the rststage estimation p erforms p o orly on the negativ ekurtosis distributions. Since Unimodal3 uses OCPF as its rststage estimator, problems in lo cating the p eak domain can parallel problems in OCPF estimations. Fig. 45 sho ws an OCPF estimation of the rst deriv ativ e of the ra w cum ulativ e distribution. The output for the rststage estimation do es not include the true mo de p osition. Therefore the secondstage will fault on estimating the mo de. When the sk ewness increases, Unimodal3 's rststage estimator no longer has the problem of excluding the mo de, and therefore the fractional error decreases. F or the sharpp eak ed families, Unimodal2 sho ws no bias in estimating the mo de. F or this family of distributions, Fig. 48 sho ws the bias for all sk ewness, and Fig. 46 sho ws 2 F or Unimodal2 w e ha v e used the exhaustiv e searc h metho d on decimated data. 54
PAGE 66
the histograms of these biases for the extremes in sk ewness. Ho w ev er, for the ratp eak ed families, Unimodal2 do es sho w a bias as the sk ewness increases. Fig. 48 sho ws the bias for all sk ewnesses, and Fig. 47 sho ws the histograms of these biases for the extremes in sk ewness. F or the sharpp eak ed families, Unimodal3 sho ws no bias in estimating the mo de. Fig. 51 sho ws this bias for all sk ewnesses, and Fig. 50 sho w the histograms of these biases for the extremes in sk ewness. Ho w ev er, for the ratp eak ed families, Unimodal3 do es sho w a bias as the sk ewness increases. Fig. 51 sho ws this bias for all sk ewnesses, and Fig. 49 sho ws the histograms of these biases for the extremes in sk ewness. 55
PAGE 67
Figure 38. Results for a family of b eta distributions p = 1 : 1 dra wing 1000 deviates. Eac h p oin t has 500 realizations. See also T able 3. F or Unimodal2 w e ha v e used the exhaustiv e searc h metho d on decimated data. 56
PAGE 68
Figure 39. Results for a family of b eta distributions with p = 1 : 1 dra wing 200 deviates. Eac h p oin t has 500 realizations. See also T able 3. F or Unimodal2 w e ha v e used the exhaustiv e searc h metho d on decimated data. 57
PAGE 69
Figure 40. Results for a family of b eta distributions with p = 1 : 8 dra wing 1000 deviates. Eac h p oin t has 500 realizations. See also T able 4. F or Unimodal2 w e ha v e used the exhaustiv e searc h metho d on decimated data. 58
PAGE 70
Figure 41. Results for a family of b eta distributions with p = 1 : 8 dra wing 200 deviates. Eac h p oin t has 500 realizations. See also T able 4. F or Unimodal2 w e ha v e used the exhaustiv e searc h metho d on decimated data. 59
PAGE 71
Figure 42. In the graph to the left is an OCPF estimation on a b eta distribution ha ving p = 1 : 1 and q = 1 : 1 tak en from a 1000deviate sample. T o the righ t is an example of ho w the p olynomial t has to undersmo oth a rat p eak distribution in order to t the sharp sides. Figure 43. The solid line is the actual b eta distribution with p = 1 : 1 and q = 1 : 1 ha ving 1000 deviates. The dotted line is the OKS estimation. 60
PAGE 72
Figure 44. Plotted as x's on the x axis is the 500deviate sample from a b eta distribution ha ving p = 1 : 1 and q = 1 : 1. The actual distribution is sho wn in the top left hand corner. 61
PAGE 73
Figure 45. The graph on the righ t sho ws an OCPF to a 1000deviate sample dra wn from a b eta distribution with p = 1 : 1 and q = 1 : 1. The solid line is the underlying distribution. The graph on the left is the output from Unimodal3 's rststage estimation. The true mo de p osition at .5 is not in the data domain. Figure 46. The histogram to the left sho ws the absolute errors of mo de estimation divided b y the standard deviation for Unimodal2 on a b eta distribution with p = 1 : 8 and q = 1 : 8 and a sk ewness r 1 = 0. Eac h of the histograms has 500 realizations of 1000deviate samples. The histogram to the righ t sho ws the absolute errors of mo de estimation for Unimodal8 on a b eta distribution with p = 1 : 8 and q = 50 and a sk ewness r 1 = 1 : 37243. Again there are 500 realizations of 1000deviate samples. F or Unimodal2 w e ha v e used the exhaustiv e searc h metho d on decimated data. 62
PAGE 74
Figure 47. The histogram to the left sho ws the absolute errors of mo de estimation divided b y the standard deviation for Unimodal2 on a b eta distribution with p = 1 : 1 and q = 1 : 1 and a sk ewness r 1 = 0. The histogram to the righ t sho ws the absolute errors of mo de estimation for Unimodal2 on a b eta distribution with p = 1 : 1 and q = 2 and a sk ewness r 1 = 0 : 481818. Both are from 500 realizations of 1000deviate samples. F or Unimodal2 w e ha v e used the exhaustiv e searc h metho d on decimated data. Figure 48. The graph to the left sho ws the absolute errors divided b y standard deviation for Unimodal2 o v er all sk ewnesses, for the family of b eta distributions with p = 1 : 1 and q v arying from 1.1 to 2. The graph to the righ t is the absolute errors for Unimodal2 o v er all sk ewnesses for the family of b eta distributions with p = 1 : 8 and q v arying from 1.8 to 50. Eac h are from 500 realizations of 1000deviates samples. F or Unimodal2 w e ha v e used the exhaustiv e searc h metho d on decimated data. 63
PAGE 75
Figure 49. The histogram to the left sho ws the absolute errors of mo de estimation divided b y the standard deviation for Unimodal3 on a b eta distribution with p = 1 : 1 and q = 1 : 1 and a sk ewness r 1 = 0. The histogram to the righ t is the absolute errors of mo de estimation divided b y the standard deviation for Unimodal3 on a b eta distribution with p = 1 : 1 and q = 2 and a sk ewness r 1 = 0 : 481818. Eac h are tak en from 500 realizations of 1000deviate samples. Figure 50. The histogram to the left sho ws the absolute errors of mo de estimation divided b y the standard deviation for Unimodal3 on a b eta distribution with p = 1 : 8 and q = 1 : 8 and a sk ewness r 1 = 0. The histogram to the righ t is the absolute errors of mo de estimation divided b y the standard deviation for Unimodal3 on a b eta distribution with p = 1 : 8 and q = 50 and a sk ewness r 1 = 1 : 37243. Eac h are tak en from 500 realizations of 1000deviate samples. 64
PAGE 76
Figure 51. The graph to the left sho ws the absolute errors for Unimodal3 o v er all sk ewnesses for the family of b eta distributions with p = 1 : 1 and q v arying from 1.1 to 2. The graph to the righ t is the absolute errors for Unimodal3 o v er all sk ewnesses for the family of b eta distributions with p = 1 : 8 and q v arying from 1.8 to 50. Eac h is tak en from 500 realizations of 1000deviate samples. 65
PAGE 77
CHAPTER 7 CONCLUSION AND FUTURE DIRECTIONS Unimodal3 and Unimodal2 sho w a bias in their estimates on distributions ha ving rat p eaks and zero sk ewness, whic h con tributed to their large fractional errors. These errors w ere larger than estimates on more highly sk ew ed data. This w as con trary to our initial assumption that a higher sk ewness w ould lead to larger fractional errors. Most estimators, excluding OKS, ga v e large fractional errors for zero sk ewness and smaller fractional errors as the sk ewness increased. OKS and OCPF pro v ed to b e the most ecien t estimators. P ossible future w ork could include a c hange in the rststage estimator for Unimodal3 It w as sho wn that OCPF fails in extracting the true p eak domain for families of distributions ha ving a rat p eak and zero sk ewness. This con tributed to the magnitude of Unimodal3 's secondstage estimation errors. Ho w ev er, OKS did not suer in lo cating the p eak of rat distributions with zero sk ewness, and therefore using this metho d in Unimodal3 's rststage estimator ma y help in correctly extracting the true p eak domain for ratp eak ed distributions. Problems to o v ercome, if OKS is used for the Unimodal3 's rststage estimation, w ould include extracting the domain of negativ e curv ature from OKS. 66
PAGE 78
REFERENCES [1] N. Johnson and S. Kotz, Continuous univariate distributions Wiley 1970. [2] E.S. Keeping, Intr o duction to statistic al infer enc e V an Nostrand, 1962. [3] C. La wson and R. Hanson, Solving le ast squar es pr oblems Pren ticeHall, 1974. [4] C. Lunneb org, Data analysis by r esampling: Conc epts and applic ations Duxbury 2000. [5] P Pitman, Pr ob ability Springer, 1992. [6] William H. Press, Saul A. T euk olsky William T. V etterling, and Brian P Flannery Numeric al r e cip es in c Cam bridge, 1992. [7] Da vid A. Rabson, B.N. Narozhn y and A.J. Millis, Cr ossover fr om p oisson to wignerdyson level statistics in spin chains with inte gr ability br e aking Submitted for publication. [8] Bro Rasm us, Nic holaos D. Sidirop oulos, and Age K. Smilde, Maximum likeliho o d tting using or dinary le ast squar es algorithms J. Chemometrics 16 (2002), 387{400. [9] Kurt S. Riedel, Pie c ewise c onvex function estimation and mo del sele ction Appro ximation theory VI I I, V ol. 1 (College Station, TX, 1995), Ser. Appro x. Decomp os., v ol. 6, W orld Sci. Publishing, Riv er Edge, NJ, 1995, pp. 467{475. [10] Alexander Sidorenk o and Kurt S. Riedel, A daptive kernel estimation of the sp e ctr al density with b oundary kernel analysis Appro ximation theory VI I I, V ol. 1 (College Station, TX, 1995), Ser. Appro x. Decomp os., v ol. 6, W orld Sci. Publishing, Riv er Edge, NJ, 1995, pp. 519{527. [11] J. T a ylor, A n intr o duction to err or analysis Univ ersit y Science Bo oks, 1982. [12] D. W ac k erly W. Mendenhall I I I, and R. Sc heaer, Mathematic al statistics with applic ations Duxbury 1996. [13] M. W and and M. Jones, Kernel smo othing Chapman & Hall, 1995. 67
PAGE 79
APPENDICES 68
PAGE 80
App endix A Momen ts, Prop erties and Cum ulan ts Man y distributions ha v e the same mean and standard deviation. Because of this fact w e m ust consider other n umerical descriptions that c haracterize data set and its distribution. The ra w momen ts 0n of a data set dene a probabilit y distribution. The n th momen t of a random v ariable Y tak en ab out the origin is dened b y E ( Y n ) = Z y n p ( Y = y ) dy = 0n (40) where E ( Y ) is the exp ectation v alue of Y = y The momen tgenerating function m ( t ), pac ks all the momen ts of a random v ariable Y in to a nice simple expression and is dened b y m ( t ) = E ( e tY ) : (41) This function exists if there is a p ositiv e constan t b suc h that m ( t ) is nite for j t j b The imp ortance of m ( t ) is that if it exists, w e can calculate an y of the momen ts for Y F or example, if m ( t ) exists, then for an y k th deriv ativ e with resp ect to t d k m ( t ) dt k = m ( k ) (0) = 0k (42) when k is p ositiv e and t is set to zero. Ho w ev er when the range in the distribution is large, the higherorder momen ts increase in size rapidly Therefore instead of the momen ts b eing expressed ab out zero, they ma y also b e expressed ab out the mean of the distribution E ( Y n ) = Z ( y ) n P ( Y = y ) dy = n : (43) If one can nd the momen tgenerating function of a random v ariable Y then one ma y use the ra w momen ts, or the momen ts ab out the mean, to determine the prop erties of the probabilit y distribution suc h as the mean, v ariance, sk ewness, and kurtosis. These 69
PAGE 81
App endix A (Con tin ued) prop erties expressed as functions of the momen ts ab out the mean are as follo ws: = 1 2 = 2 r 1 = 3 3 = 2 2 r 2 = 4 3 22 22 : (44) W e describ e a probabilit y distribution b y its mean, v ariance, sk ewness, and kurtosis. The mean of a sample, x i is giv en b y = 1 N N X i =1 x i (45) where N is the total n um b er of sample measuremen ts. This is simply the a v erage o v er all of the sample measuremen ts. The v ariance of a sample x i is giv en b y 2 = 1 ( N 1) N X i =1 ( x i ) 2 : (46) The v ariance of a sample of measuremen ts is the sum of the square dierences b et w een a measuremen ts and their mean, divided b y N 1. One ma y also use sk ewness and kurtosis to describ e a probabilit y distribution. The sk ewness is a measure of asymmetry in a sample distribution. Fig. 52 sho w a family of W eibull distributions with dieren t sk ewnesses. The kurtosis is a measure of ratness or sharpness of the p eak of a probabilit y distribution. 70
PAGE 82
App endix A (Con tin ued) Figure 52. W eibull distributions with a = 1 and c ranging from 1.4 to 2.6 As the c parameter increases the p eak p osition mo v es a w a y from the y axis. Therefore when c = 1 : 4, the distribution has a larger sk ewness, and the p eak is closer to the y axis. Another n umerical description of the uniqueness of a distribution, is b y the use of the cum ulan ts n The cum ulan ts are deriv ed from the ra w momen ts, 0n and are dened b y ln[ ( t )] = 1 X n =0 n ( it ) n n (47) where ( t ) is the c haracteristic function dened as the F ourier transform of the probabilit y densit y function p ( x ) ( t ) = F x [ p ( x )]( t ) = Z 1 1 e itx p ( x ) dx: (48) 71
PAGE 83
App endix A (Con tin ued) Here the cum ulan ts are the real co ecien ts of a Maclaurin series of ln[ ( t )]. The rst four cum ulan ts expressed as functions of ra w momen ts n are as follo ws: 1 = 1 2 = 2 21 3 = 2 31 3 1 2 + 3 4 = 6 41 + 12 21 2 3 22 4 1 3 + 4 : (49) The unique prop erties of a distribution suc h as mean, v ariance, sk ewness, and kurtosis can also b e expressed as a com bination of cum ulan ts = 1 2 = 2 r 1 = 3 (3 = 2) 2 r 2 = 4 22 : (50) 72
PAGE 84
App endix B T est Probabilit y Distributions The W eibull distribution is named after W alo ddi W eibull, who oered it as an analytical to ol for mo deling the breaking strengths of materials. Other uses include reliabilit y and lifetime mo deling [5 ]. The W eibull distribution is dened as W ( x; a; c ) = ca 1 ( x a ) ( c 1) exp[ ( x a ) 2 ] ; (51) where a; c > 0 and 0 x < 1 The a parameter alters the v ariance and the scale of the distribution, while the c parameter alters the sk ewness as w ell as the kurtosis. The momen ts of the W eibull distribution ab out the mean are 1 = a (1 + c 1 ) 2 = a 2 (1 + 2 c 1 ) 3 = a 3 (1 + 3 c 1 ) 4 = a 4 (1 + 4 c 1 ) ; (52) and the mean, v ariance, sk ewness and kurtosis are = a (1 + c 1 ) (53) 2 = a 2 [(1 + 2 c 1 ) (1 + c 1 ) 2 ] (54) r 1 = 2(1 + c 1 ) 3(1 + c 1 )(1 + 2 c 1 ) [(1 + 2 c 1 ) (1 + c 1 ) 2 ] (3 = 2) + (1 + 3 c 1 ) [(1 + 2 c 1 ) (1 + c 1 )] (3 = 2) : (55) r 2 = 6(1 + c 1 ) 4 + 12(1 + c 1 ) 2 (1 + 2 c 1 ) 3(1 + 2 c 1 ) 2 4(1 + c 1 )(1 + 3 c 1 ) + (1 + 4 c 1 ) [(1 + 2 c 1 ) (1 + c 1 ) 2 ] 2 : (56) 73
PAGE 85
App endix B (Con tin ued) Here ( p ) is the Gamma function dened as ( p ) = Z 1 0 t p 1 exp[ t ] dt: (57) If w e in tegrate ( p ) from 0 t < 1 w e can see that (1) = 1 ; and for an y p > 1, ( n ) = ( n 1)!, as long n is an in teger. The Gamma distribution is based on t w o parameters and is used to mo del sk ew ed frequency distributions. Suc h examples of sk ew ed frequency p opulations include the lengths of time b et w een malfunctions for aircraft engines, and arriv als at a sup ermark et c hec k out coun ter [12 ]. The Gamma distribution is dened as G ( x; a; b ) = 1 b a ( a ) x ( a 1) exp[ ( x b )] : (58) The sk ewness and kurtosis of the distribution are go v erned b y the v alue of a while the mean and v ariance dep end on b oth a and b The momen ts of the distribution ab out the mean are 1 = b (1 + a ) ( a ) 2 = b 2 (2 + a ) ( a ) 3 = b 3 (3 + a ) ( a ) 4 = b 4 (4 + a ) ( a ) ; (59) and the mean, v ariance, sk ewness, and kurtosis are = ab (60) 2 = ab 2 (61) 74
PAGE 86
App endix B (Con tin ued) r 1 = 2 p ( a ) (62) r 2 = 6 a : (63) Unlik e the gamma and W eibull distributions, the b eta distribution is dened o v er a closed in terv al 0 x 1. The b eta distribution function is often used to mo del prop ortions. Suc h examples include the prop ortion of c hemical impurities in a sample, and the prop ortion of time a mac hine is under repair [12 ]. The b eta distribution is dened as B ( x; p; q ) = x p 1 (1 x ) q 1 ( p; q ) ; (64) where ( p; q ) is the Beta function, ( p; q ) = Z 1 0 t p 1 (1 t ) q 1 dt; (65) and p; q > 0. Compared to the gamma and W eibull distributions, the b eta distribution can assume a wide range of shap es. The ra w momen ts ab out the mean for the b eta distribution are 1 = (1 + p )( p ) ( p; q )(1 + p + q ) 2 = (2 + p )( p ) ( p; q )(2 + p + q ) 3 = (3 + p )( p ) ( p; q )(3 + p + q ) 4 = (4 + p )( p ) ( p; q )(4 + p + q ) ; (66) and the mean, v ariance, sk ewness, and kurtosis are = p ( p + q ) (67) 75
PAGE 87
App endix B (Con tin ued) 2 = pq ( p + q ) 2 ( p + q + 1) (68) r 1 = 2( q p ) ( p + q + 2) r p + q + 1 pq (69) r 2 = 6( a 3 + a 2 (1 2 c ) + c 2 (1 + c ) 2 ac (2 + c )) ac (2 + a + c )(3 + a + c ) : (70) 76
PAGE 88
App endix C Cheb yshev P oly omial The Cheb yshev p olynomial is dened as T n ( x ) = cos ( n arccos x ) ; (71) where n is the degree of the p olynomial. T n ( x ) can ha v e explicit expressions; Press, T euk olsky V etterling and Flannery sho w the rst four for n = 0 to n = 3 [6 ]: T 0 ( x ) = 1 T 1 ( x ) = x T 2 ( x ) = 2 x 2 1 T 3 ( x ) = 4 x 3 3 x T n + 1( x ) = 2 xT n ( x ) T n 1( x ) n 1 : (72) These p olynomials are orthogonal on the the in terv al [1,1] o v er a w eigh t (1 x 2 ) 1 = 2 Z 1 1 T i ( x ) T j ( x ) p 1 x 2 dx = 8>>>><>>>>: 0 : i 6 = 2 : i = j 6 = 0 : i = j = 0 (73) and ha v e n zeros on the in terv al [1,1] lo cated at x = cos ( k 1 = 2) n k = 1 ; 2 ; ::::; n: (74) The Cheb yshev p olynomials not only satisfy a con tin uous relation, Equ. 73, but also a discrete orthogonalit y relation. If x k ( k = 1 ; :::; m ) are the m zeros of T m giv en b y equations 77
PAGE 89
App endix C (Con tin ued) 5, and if i; j < m then m X k =1 T i ( x k ) T j ( x k ) = = 8>>>><>>>>: 0 : i 6 = m= 2 : i = j 6 = 0 m : i = j = 0 (75) Press, T euk olsky V etterling, and Flannery sho w us that for an y arbitrary function on the in terv al [1,1], and if N co ecien ts, c j where j = 0 ; ::::; N 1, are dened b y c j = 2 N N X k =1 f ( x k ) T j ( x k ) (76) then the appro ximation form ula is f ( x ) [ N 1 X k =0 c k T k ( x )] 1 2 c 0 (77) and is exact for x equal to all of the N zeros of T N ( x ). Ho w ev er not all of our distributions are dened on the in terv al [1,1]. Therefore w e m ust use an eectiv e c hange of v ariable giv en b y x 0 x 1 2 ( b + a ) 1 2 ( b a ) (78) where a and b are the arbitrary limits instead of [1,1]. 78
PAGE 90
App endix D Finite Dierences In Unimodal2 w e imp osed constrain ts on the 2 nd deriv ativ e, and in Unimodal3 w e imp osed constrain ts on the 3 r d deriv ativ e. W e use the nite dierences of a data set, x i where i = 0 :::N 1, in order to estimate the deriv ativ es. W e express the 1 st 2 nd and 3 r d nite dierences of a data set x i as a i = f i +1 f i i (79) b i = a i +1 a i i (80) c i = b i +1 b i i (81) where i = x i +1 x i and f i is the discrete cum ulativ e distribution f i = 1 N 1 i : (82) W e wish to express the function f i in terms of the nite dierences. W e can in v ert (79) (80) and (81) to get f i = f i 1 + a i 1 i 1 (83) a i = a i 1 + b i 1 i 1 (84) b i = b i 1 + c i 1 i 1 : (85) If w e are using Unimodal2 w e will express f i in terms of a i and b i (second and third dierences) and get f i = f i 1 + i (86) 79
PAGE 91
App endix D (Con tin ued) where i = a i 1 i 1 By recursiv ely plugging in the expression for a i 1 tak en from (84) in to (83) w e get i = h a 0 + i 2 X k =0 b k k i i 1 : (87) The rst three expressions for f i are f 0 = f 0 f 1 = f 0 + a 0 0 f 2 = f 0 + a 0 0 + a 0 1 + b 0 0 1 : :: (88) W e can write (88) in matrix form. The nal problem is reduced to a design matrix for leastsquares minimization 1 expressed as D w = f (89) where D a matrix dened b y (88) and the v ector w is the v ector with comp onen ts w 0 = f 0 w 1 = a 0 w j = b j 2 w her e 2 j p 0 + 2 w j = b j 2 w her e j > p 0 + 2 (90) and p 0 is the index of the assumed mo de p osition x p 0 1 Refer to Numerical Recip es in C for denition of design matrix. 80
PAGE 92
App endix D (Con tin ued) The leastsquares minimization of R 2 = X i ([ D w ] i f i ) 2 (91) is iden tical to the problem (89) [6 ]. If w e replace w with e w w e can write a new iden tit y D e w = e f : (92) If w e enforce no constrain ts on e w then e w = w e f = f and e f w ould reconstruct the original discrete cum ulativ e distribution f Requiring e w to b e nonnegativ e leads to a constrained leastsquares problem, and e f is the closest cum ulativ e probabilit y distribution ( in the leastsquares sense) to the original data. If w e are using Unimodal3 then w e ha v e to extend (87) to include the third nite dierences, c i : i = h a 0 + i 2 X k =0 b 0 + k 1 X l =0 c l l k ] i 1 (93) where the problem is again reduce to a design matrix for leastsquares minimization, and w comp onen ts are no w w 0 = f 0 w 1 = a 0 w 2 = b 0 w 3 = c 0 w 4 = c 1 : : :: (94) 81
PAGE 93
App endix D (Con tin ued) Again w e can replace w with e w and ha v e a new iden tit y as in (92) where e w is constrained to b e nonnegativ e. 82
