1 / 55

# Environmental Data Analysis with MatLab

Environmental Data Analysis with MatLab. Lecture 8: Solving Generalized Least Squares Problems. SYLLABUS.

Download Presentation

## Environmental Data Analysis with MatLab

An Image/Link below is provided (as is) to download presentation Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

### Presentation Transcript

1. Environmental Data Analysis with MatLab Lecture 8: • Solving Generalized Least Squares Problems

2. SYLLABUS Lecture 01 Using MatLabLecture 02 Looking At DataLecture 03Probability and Measurement ErrorLecture 04 Multivariate DistributionsLecture 05Linear ModelsLecture 06 The Principle of Least SquaresLecture 07 Prior InformationLecture 08 Solving Generalized Least Squares ProblemsLecture 09 Fourier SeriesLecture 10 Complex Fourier SeriesLecture 11 Lessons Learned from the Fourier Transform Lecture 12 Power SpectraLecture 13 Filter Theory Lecture 14 Applications of Filters Lecture 15 Factor Analysis Lecture 16 Orthogonal functions Lecture 17 Covariance and AutocorrelationLecture 18 Cross-correlationLecture 19 Smoothing, Correlation and SpectraLecture 20 Coherence; Tapering and Spectral Analysis Lecture 21 InterpolationLecture 22 Hypothesis testing Lecture 23 Hypothesis Testing continued; F-TestsLecture 24 Confidence Limits of Spectra, Bootstraps

3. purpose of the lecture use prior information to solve exemplary problems

4. review of last lecture

5. failure-proof least-squares add information to the problem that guarantees that matrices like [GTG] are never singular such information is called prior information

6. examples of prior information soil has density will be around 1500 kg/m3 give or take 500 or so chemical components sum to 100% pollutant transport is subject to the diffusion equation water in rivers always flows downhill

7. linear prior information with covariance Ch

8. simplest examplemodel parameters near known values Hm= h with H=I h = [10, 20]T Ch= m1 = 10 ± 5 m2 = 20 ± 5 m1 and m2 uncorrelated

9. use Normal p.d.f. to representprior information

10. Normal p.d.f.defines an“error in prior information” individualerrors weighted by their certainty

11. now suppose that we observe some data: d= dobs with covariance Cd

12. represent the observations with a Normal p.d.f. p(d) = observations mean of data predicted by the model

13. this Normal p.d.f.defines an“error in data” weighted by its certainty prediction error

14. Generalized Principle of Least Squaresthe best mest is the one thatminimizes the total error with respect to mjustified by Bayes Theoremin the last lecture

15. generalized least squaressolution pattern same as ordinary least squares … … but with more complicated matrices

16. Generalized least squares is equivalent to solvingF m = fby ordinary least squares = m

17. top partdata equation weighted by its certainty = m σd-1{Gm=d } data equation certainty of measurement

18. bottom partprior information equation weighted by its certainty = m σh-1{Hm=h } prior information equation certainty of prior information

19. exampleno prior information butdata equation weighted by its certainty = m called “weighted least squares”

20. straight line fitno prior information butdata equation weighted by its certainty fit data with low variance data with high variance

21. straight line fitno prior information butdata equation weighted by its certainty fit data with low variance data with high variance

22. another exampleprior information that the model parameters are small m≈0H=Ih=0assume uncorrelated with uniform variancesCd = σd2I Ch = σh2I

23. Fm=h m = [FTF]-1FTm=f • m=[GTG + ε2I]-1GTd with ε= σd/σm

24. called “damped least squares” • m=[GTG + ε2I]-1GTd with ε= σd/σm • ε=0: minimize the prediction error • ε→∞: minimize the size of the model parameters • 0<ε<∞: minimize a combination of the two

25. advantages:really easy to codemest = (G’*G+(e^2)*eye(M))\(G’*d);always works • m=[GTG + ε2I]-1GTd with ε= σd/σm • disadvantages: • often need to determine ε empirically • prior information that the model parameters are small not always sensible

26. smoothness as prior information

27. model parameters represent the values of a functionm(x)at equally spaced increments along the x-axis

28. function approximated by its values at a sequence of x’s mi mi+1 m(x) x xi xi+1 Δx • m(x) → m=[m1, m2, m3, …, mM]T

29. rough function has large second derivative a smooth functionis one that is not rough a smooth function has a small second derivative

30. approximate expressions for second derivative

31. m(x) x xi • i-th row of H: (Δx)-2 [ 0, 0, 0, … 0, 1, -2, 1, 0, …. 0, 0, 0] 2nd derivative at xi column i

32. what to do about m1 and mM?not enough points for 2nd derivativetwo possibilitiesno prior information for m1 and mMorprior information about flatness(first derivative)

33. m(x) x x1 • first row of H: (Δx)-1 [ -1, 1, 0, … 0] 1st derivative at x1

34. example problem: • to fill in the missing model parameters so that • the resulting curve is smooth m = d x

35. the model parameters, man ordered list of all model parameters m=

36. data equationGm=d = data kernel “associates” a measured model parameter with an unknown model parameter data are just model parameters that have been observed

37. The prior information equation, Hm=h • “smooth interior” / “flat ends” • h=0

38. put them together into theGeneralized Least Squares equation F = f = • choose σd/σm to be << 1 • data takes precedence over prior information

39. the solution using MatLab

40. graph of the solution solution passes close to data solution is smooth m = d x

41. Two MatLab issues Issue 1: matrices like G and F can be quite big, but contain mostly zeros. Solution 1: Use “sparse matrices” which don’t store the zeros Issue 2: matrices like GTG and FTF are not as sparse as G and F Solution 2: Solve equation by a method, such as “biconjugate gradients” that doesn’t require the calculation of GTG and FTF

42. Using “sparse matrices” which don’t store the zeros: N=200000; M=100000; F=spalloc(N,M,3*M); creates a 200000× 100000 matrix that can hold up to 300000 non-zero elements. “sparse allocate” note that an ordinary matrix would have 20,000,000,000 elements

43. Once allocated, sparse matrices are used just like ordinary matrices … … they just consume less memory.

44. Issue 2: Use biconjugate gradient solver to avoid calculating GTG and FTF Suppose that we want to solve FTF m = FTf The standard way would be: mest = (F’F)\(F’f); but that requires that we compute F’F

45. a “biconjugate gradient” solver requires only that we be able to multiply a vector, v, by GTG, where the solver supplies the vector, v. so we have to calculate y=GTG v the trick is to calculate t=Gv first, and then calculate y=G’t this is done in a Matlab function,afun()

More Related