PAGE 1 FIRST AN ADVERTISEMENT Neil Grossbard(x6356),proficient in - PowerPoint PPT Presentation

PAGE 1 FIRST AN ADVERTISEMENT Neil Grossbard(x6356),proficient in numerical analysis and computer programming. Thanks to the efforts of Jack Jasperse , Bamandas Basu , Paul Rothwel and many other scientists, my work has involved finding solutions to differential equations eigenvalue, eigenvector problems, integrals sometimes multi-dimensional, Pic codes,wavelets, and performed symbolic manipulation using Fortran to create Fortran code to quickly solve linear equations and take derivatives of the solution. I am an (expert?) at estimating PSD’s. (POWER SPECTRAL DENSITIES) and I am available to help in most numerical or simulation problems. SEE "Power spectral artifacts in published balloon data and implications regarding saturated gravity wave theories" J.G.R. VOL 105, NO. D4,PAGES 4667-4683 Feb 27,2000 BY Dewan and Grossbard

PAGE 2 MONTE CARLO METHODS,JACKKNIFE AND BOOTSTRAP TO FIND ERROR ESTIMATES. IN PARTICULAR FOR VARIOUS METHODS OF ESTIMATING PSD’S. (POWER SPECTRAL DENSITIES) BY NEIL GROSSBARD

PAGE 3 I AM NEIL GROSSBARD. IN THE PREVIOUS CENTURY I RECEIVED THE FOLLOWING SCIENCE DEGREES: BACHELORS IN PHYSICS (M.I.T.) AND 3 MASTERS (NORTHEASTERN) PHYSICS,MATHEMATICS AND ELECTRICAL ENGINERING. RECENTLY MY HOBBY IS STUDYING NUMERICAL METHODS, ESPECIALLY SIGNAL PROCESSING AND MACRO-ECONOMICS. DOES ANYONE WISH TO ARGUE,TALK,E-MAIL ABOUT ANY OF THESE TOPICS? I THINK THE MOST USFUL PART OF THIS PRESENTATION IS AN EXPLAINATION OF MONTE CARLO,JACKKNIFE AND BOOTSTRAP METHODS FOR ERROR ANALYSIS.

PAGE 4 Consider any data which is used to estimate some parameter(s) how can you use Monte Carlo to assign errors to the results? THE PLAN First consider a simple Monte Carlo problem with no correlation. Second how to estimate PSD (autocorrelation). Then how to do Monte Carlo with correlated data. Finally I will suggest the use of autoreggessive methods to estimate the PSD (BURG).

PAGE 5 What is the Monte-Carlo Method Using Google from "A Practical Guide to Monte Carlo Simulation" by Jon Wittwer "A Monte Carlo method is a technique that involves using random numbers and probability to solve problems" I expect that usually you would use Monte Carlo when you use a possibly complicated procedure to calculate values from experimental data and wish to assign errors to the final results. In Monte Carlo it is hoped that you have an estimate for the errors in the experimental data. You then can add an error to the actual data using psudo-random values and process the modified data to get a new result then repeat this proceedure many times and keep track of each result. You then can calculate estimates of the error statistics. (Histograms,standard deviation, confidence regions etc...) I will now show a simple example. The problem is to find the formula for a straight line which is estimated from a set of data points where the relative accuracy of each data points is known. This is a weighted least squares problem.

PAGE 6 NUMBER OF MONTE CARLO SAMPLES=10000 4 CONSTANT= 1.8336 VARIANCE= 0.0017 3 ACTUAL 1.8333 ACTUAL 0.0017 2 1 0 1 2 3 4 STRAIGHT LINE -3 -2 -1 0 1 2 3 LINE 1 STAN DEV STANDARD DEVIATIONS 4 SLOPE= 0.9994 VARIANCE= 0.0020 ACTUAL 1.0000 ACTUAL 0.0020 3 2 1 0 1 2 3 4 -3 -2 -1 0 1 2 3 ESTIMATED VALUES STANDARD DEVIATIONS

PAGE 7 NUMBER OF MONTE CARLO SAMPLES=10000 4 AVERAGE CONSTANT= 1.8335 AVERAGE VARIANCE= 0.0017 3 ACTUAL 1.8333 ACTUAL 0.0017 2 1 0 1 2 3 4 STRAIGHT LINE -3 -2 -1 0 1 2 3 LINE 1 STAN DEV STANDARD DEVIATIONS 4 AVERAGE SLOPE= 0.9998 AVERAGE VARIANCE= 0.0021 ACTUAL 1.0000 ACTUAL 0.0020 3 2 1 0 1 2 3 4 -3 -2 -1 0 1 2 3 ESTIMATED VALUES STANDARD DEVIATIONS

PAGE 8 The Jackknife, the Bootstrap and Other Resampling Plans BRADLEY EFRON Department of Statistics Stanford University CBMS-NSF REGIONAL CONFERENCE SERIES IN APPLIED MATHEMATICS

PAGE 9 The following is a fancy way of showing how to estimate the standard deviation using Jackknife (obvious?). ∧ n n-1 σ JACK = ( 2 ) - (i) -x - (.) ) 1/2 ∑ i=1 (x n The following is a fancy way of showing how to estimate the standard deviation using Bootstrap (obvious?). ∧ _ { } - . 2 1/2 ] B ∧ ∧ 1 ∗ [θ θ ∗ ∑ b SD= B-1 b=1 DID NOT KNOW WHERE TO PUT THIS Also note you can get the Discrete Fourier Transform for any number of values (NO 2^N) using the Chirp z- Transform Algorithm. Found in Digital Signal Processing by Openheim and Schafer p321.

NUMBER OF SAMPLES FOR DRAW=20000 EACH DRAW=100 PAGE 10 NUMBER OF MONTE CARLO SAMPLES=10000 4 AVERAGE CONSTANT= 1.8341 AVERAGE VARIANCE= 0.0017 3 ACTUAL 1.8333 ACTUAL 0.0017 2 1 0 1 2 3 4 STRAIGHT LINE -3 -2 -1 0 1 2 3 LINE 1 STAN DEV STANDARD DEVIATIONS 4 AVERAGE SLOPE= 1.0014 AVERAGE VARIANCE= 0.0021 ACTUAL 1.0000 ACTUAL 0.0020 3 2 1 0 1 2 3 4 -3 -2 -1 0 1 2 3 ESTIMATED VALUES STANDARD DEVIATIONS

PAGE 11 RESULTS OF LINE FIT CONSTANT SLOPE SIMULATED 1.8333E+00 1.0000E+00 USING MONTE CARLO 1.8336E+00 9.9940E-01 USING SEPARATE MONTE CARLO 1.8335E+00 9.9982E-01 USING BOOT STRAP MONTE CARLO 1.8341E+00 1.0014E+00 To understand classical PSD estimation we need the convolution theorums. CONTINUOUS CONVOLUTION INTEGRAL ∞ ∞ ∫ ∫ τ τ τ τ τ τ g(w)= h( )f(w- )d = h(w- )f( )d ∞ ∞ - - DISCRETE CONVOLUTION SUMMATION ∞ ∑ [ ] [ ] [ ] [ ] ∗ [ ] g n = h n-k f k =h n f n ∞ - CONVOLUTION IN FREQUENCY DOMAIN IS MULTIPLICATION IN TIME DOMAIN. [ ] Consider some h n functions (tapers).

PAGE 12 TAPER TAPER EXPANDED SCALE FFT EXPANDED SCALE FFT EXPANDED SCALE FFT

PAGE 13 TIME-BANDWIDTH PRODUCT=2 THETA= 5.72E-05 LAMBDA= 1.00E+00 TAPER EXPANDED SCALE FFT THETA= 2.44E-03 LAMBDA= 9.98E-01 TAPER EXPANDED SCALE FFT FROM MULTITAPER SUBROUTINE

PAGE 14 Now a fundimental problem when we use the Periodogram to estimate the PSD. From Digital Signal Proceesing by Rabiner and Gold { } [ ω ] 2 1+ ( sin N α ) The variance of the Periodogram= 4 ω N sin ω Here N is the number of points and is the frequency. α is the standard deviation of the data. All classical methods use some averaging technique (explicitly or implicitly) to cut this variation. All the methods I will mention assume the data is collored Gaussian. This is generally a good assumption due to the central limit theorum. When this is not true the simplest technique is to approximate the distribution as a sum of gaussians. I had heard of this before but I was reminded by Richard Hegbloom (thank you).

PAGE 15 Some methods of classical PSD estimation. All of the estimates attempt to cut the variance of the PSD estimates at the expense of increasing the resolution frequency. Original Blackman-Tukey Calculate the biased estimate of the auto-correlation function. ∑ τ n is the number of data points. X(i)*X(i+ ) τ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ τ A( )= varies between -AND+ N (a number <n). n Use all possible i values Apply a Bartlett (Triangular) Window. The Triangle goes to 0 long before the end of the possible auto-correlation estimates. Tukey suggests something like 1/10 of the data length. Apply a cosine transform to the windowed auto-correlation. Welch’s Method Divide the data into usually overlapping segments. For each segment apply your favorite window. Find the magnitude squared of the discrete Fourier Transform for each windowed segment. Average the separate PSD estimates. Welch is recently the favored method, even though the next method blows it out of the water.

PAGE 16 The Multitaper method of PSD estimation. SEE THE FOLLOWING SPECTRAL ANALYSIS FOR PHYSICAL APPLICATIONS Multitaper and Conventional Univariate Techniques By Donald B. Percival and Andrew T. Walden The general idea for Multitaper is to get PSD estimates from a set of windows each separately applied to the initial data. In particular Thomson developed a set of multiple orthogonal tapers. Discrete prolate spheroidal sequence data tapers. A nice feature of this method and Burg’s method is that in principle they can be used when data values are missing. Interpolating the missing data values for Burg’s method is usually not necessary.

THIS IS FILE RUNSIM RUN ON DATE May 14 12 KZ*= 2.09E-03 RAD/M N= 2.00E-02 RAD/SEC�˝ DX= 5.00E+01 METERS USE NUMBER OF POINTS USED IN SIM=4096 NUMBER OF POINTS USED IN ANALYSIS=256 R.N. SEED=-12345 AVERAGE OF 1000 SCENES SLOPE= 3.00E+00 FIT LINE BETWEEN FREQUENCIES 4.00E-03 AND 4.00E-02 106 105 THE PSD TO BE ESTIMATED 2 (RAD/M)-1 104 THEORETICAL PSD(M /S ) 2 103 102 101 100 10-1 10-4 10-3 10-2 10-1 WAVENUMBER(RAD/METER)

THIS IS FILE RUNSIM RUN ON DATE May 14 12 KZ*= 2.09E-03 RAD/M N= 2.00E-02 RAD/SEC�˝ DX= 5.00E+01 METERS USE NUMBER OF POINTS USED IN SIM=4096 NUMBER OF POINTS USED IN ANALYSIS=256 R.N. SEED=-12345 AVERAGE OF 1000 SCENES SLOPE= 3.00E+00 FIT LINE BETWEEN FREQUENCIES 4.00E-03 AND 4.00E-02 AVERAGE SLOPE= -2.76E+00 ESTIMATED STANDARD DEVIATION AVERAGE= 1.05E-02 EACH= 3.31E-01 106 105 2 (RAD/M)-1 THE MONTE CARLO RESULT� PERIODOGRAM ESTIMATE OF PSD(M /S ) 104 NO SPECIAL TAPER 2 NO PREWHITTENING 103 102 101 100 10-1 10-4 10-3 10-2 10-1 WAVENUMBER(RAD/METER)