These are the Department of Mathematical and Statistical Sciences Legacy Pages

This is an informal discussion on the Wishart distribution and law of large numbers for Wishart samples with application to the convergence of the sample covariance. Other topics to be discussed include laws of large numbers in Hilbert and Banach spaces.

NOTE: The seminar starts at 10AM. 1. Discussion of "The Geometry of Statistics" presentation from the Statistics Seminar

2. Central limit theorem in Banach spaces

This is a follow-up to the previous seminar "A roadmap to the proof of Kalman filter formulas in a Hilbert space". We'll try to solve some of the problems posed there.

I will introduce the basics of morphing and image registration as we use them in data assimilation.

The purpose of this talk is to give enough background and information to use Loren's code (under development).

NOTE UNUSUAL DAY.

This is mostly a hands-on coding session to question, discuss, and test combinations of our codes: 1. enkf (matlab) 2. morphing (fortran) 3. fft data assimilation (matlab) 3. epidemic simulation (java) 4. weather-wildfire (fortran) If time allows, further discussion will deal with: 1. requirements for the cross-covariance 2. covariance models in the frequency domain 3. covariogram and its computation by the FFT (Marcotte, 1996)

We'll compare our notes and codes. PLEASE NOTE UNUSUAL DAY.

This talk will begin a survey of EnKF literature. Arter the talk, we'll continue working on the ICCS10 and MS10 conference papers, connect and test codes, etc.

This talk will continue a survey of EnKF literature. Arter the talk, we'll continue working on the ICCS10 conference paper, connect and test the epidemic simulation and FFT EnKF codes, etc.

Possible topics include: 1. FFT EnKF (Jan), 2. Continue the review (Ashok), 3. Epidemic model (Loren), 4. EnKF in Hilbert spaces (Jan), a graphical EnKF demo (Volodymyr)

Oscar will give a brief presentation about the basics of 3D VAR. Further topics, time and interest permitting: background for FFTEnKF and 1D EnKF demo (Volodymyr), FFTEnKF and covariance fitting (Jan), continuation of comparison of deterministic EnKFs (Ashok), discussion of epidemic models (Loren),...

After the talk, ongoing discussion of FFT EnKF and of the epidemic and wildfire applications.

This talk will present the basics of spatial statistics: stationary random fields and covariance estimation by the variogram and the covariogram. A classical application of the variogram is Kriging.

This talk will explain the Marcotte (1996) algorithm.

In this talk, we introduce a state-space tracking algorithm, based on combined particle learning (PL) and sequential Bayesian inference. The proposed algorithm is particularly well-suited to on-line learning and surveillance of infectious diseases -- it is capable of assessing the probability of an pandemic, while simultaneously accounting for uncertainty in disease parameters and producing predictions in real-time. As compared to the now widely used MCMC-based methods, this PL method, which is based on efficient use of an essential state vector, is easier to implement, computationally faster, as well as more readily generalizable to problems with complex non-linear dynamics. We illustrate this algorithm for tracking influenza with Google Flu Trends data, taking a closer look at the spread of flu in the US during 2003-2009, and in New Zealand during 2006-2009.

Informal discussion.

I'll look into what people do about benchmarking data assimilation in the literature and possibly propose some benchmark problems. Discussion.

Discussion.

I will review the main ideas from the paper Castronovo, E. and Harlim, J. and Majda, A. J., Mathematical test criteria for filtering complex systems: plentiful observations, Journal of Computational Physics, 227, 2008, 3678--3714 and possibly some related papers by the authors on the topic.

1. Live walk through the data assimilation code: assimilation cycle, netcdf files, data structures, fftenkf, connecting a model 2. How to tell if the code works correctly: statistical validation on a simple gaussian linear model 3. Linear test models from FDKF papers 4. Nonlinear benchmarks 5. Only the last, the problems we actually want to deal with.

Work sessions follow: Jonathan Beezley - Matlab classes for FFT EnKF, Jan Mandel - Stochastic validation of EnKF

Work presentations & dicusssion follow: Jonathan Beezley - classes in data assimilation code Jan Mandel - stochastic validation Note: the seminar will move to room 626 if it is available

This is a run-through Loren's ICCS2010 presentation next week in Amsterdam. Further presentation and discussion: Jonathan Beezley, Matlab classes for ensemble filters

Followed by a work session and planning of work over Summer.

1. The EnKF EpiSim paper

2. Laplace operator and isotropic smooth random fields

3. miscellaneous

Preceded by a brief discussion of the the EnKF EpiSim paper and presentation for JSM (Ashok Krishnamurthy). NOTE: The seminar is Tuesday this week ONLY. We will return to Wednesdays afterwards.

This is part 2 of presentation and work discussion based on

J. Mandel, Optimal statistical interpolation and morphing, working notes maintained in the git repository osimorph/doc/osimorph.tex.

J. Mandel, J. D. Beezley, K. Eben, P. Jurus, V. Y. Kondratenko, and J. Resler, Data assimilation by morphing fast Fourier transform ensemble Kalman filter for precipitation forecasts using radar images, CCM Report 289, April 2010

We will also discuss presentation in preparation: Ashok Krishnamurthy, Bayesian Tracking of Emerging Epidemics Using Ensemble Kalman Filters Presentation to be given at 2010 Joint Statistical Meetings, Vancouver, Canada, July 31- August 5, 2010

Discussion topics: A scalable domain domain decomposition solver for stochastic PDEs (i.e. a two level preconditioner with a coarse grid).

Data assimilation (nonlinear filtering and some full Bayesian approaches)

Nonlinear (aeroelastic) oscillators (ODEs exhibiting Hopf bifurcation) describing the flutter of aircraft wings. We looked at parameter estimation and combined state and parameter estimation problems from Wind Tunnel Test data.

Nonlinear Shallow water equations (depth average Navier Stokes equation) for forecasting floods when a dam breaks and Advection diffusion equation for contaminant tracking. State estimation problem using Ensemble Kalman filter.

Long term interest is on how to couple the domain decomposition solver with the sequential data assimilation method

Links:

http://www.abhijitsarkar.net

http://hpcu.dyndns.org/temp/resume.pdf

http://hpcu.dyndns.org/PAPERS/domain_decomposition_preconditioner_2level_coarse _grid.pdf

http://hpcu.dyndns.org/PAPERS/domain_decomposition_preconditioner_1level.pdf

http://hpcu.dyndns.org/PAPERS/LCO_parameter_estimate.pdf

http://hpcu.dyndns.org/PAPERS/paper_Hydro_Montreal08.pdf

http://hpcu.dyndns.org/thesis/thesis_HM.pdf

http://hpcu.dyndns.org/PAPERS/Flutter_margin.pdf

Modeling extreme events in hydrology such as heavy rainfalls or extreme temperatures implies the use of heavy tailed distributions. In other cases, hydrologists look for bounded support distributions (e.g for humidity rates). None of those situations can be addressed with Gaussian distributions. Kalman filtering allows one to estimate the state of a first order auto regressive, hidden time process which residuals are Gaussian. Observations are linearly dependent on the hidden state, with Gaussian perturbations. We present a more general framework in which the hypothesis of normality can be relaxed : noises are assumed to be elliptically distributed, without any extra hypothesis about their generator's class. The need for stability of the model implies limiting ourselves to finite time series.

For hydrologic purposes, generators can be chosen among Generalized Pareto distributions (GPD). Stability properties of this elliptical family allow for explicit expressions of marginal and conditional laws (ie the estimates). An appropriate choice of parameters induces heavy tailed processes, possibly without second moments, or bounded processes.

Adapting the Kalman filter algorithm for elliptical noises does not change evolution equations for estimates and covariance matrices. Nevertheless, updated generators must be computed. Confidence regions differ from those of a Gaussian hidden state space model.

Finally, a possible application to rainfall modeling is discussed. Distributions of rainfall have positive, unbounded support, which cannot be the case for elliptical distributions. As a solution, we use an elliptical copula model instead of the original elliptical model, which allows for any margins specification.

* I ll bring the team up to date on the JSM proceedings paper * I ll give a brief literature review of Optimal Statistical Interpolation

* I ll be looking at the following articles/book chapters:

1) Book Chapter - Bayesian estimation. Optimal interpolation. Statistical linear estimation - Talagrand - 2003

2) Optimal interpolation and the Kalman filter - Cohn et al., - 1981 (Proceedings of the Fifth Conference on Numerical Weather Prediction)

3) Multivariate Error Covariance Estimates by Monte Carlo Simulation for Assimilation Studies in the Pacific Ocean - Borovikov et al., - 2005

4) Ensemble Optimal Interpolation - multivariate properties in the Gulf of Mexico - Counillon and Bertino - 2009

5) Relationships between statistical and deterministic methods of data assimilation - Thacker - 1986 (Proceedings of the International Symposium on Variational Methods in Geosciences)

6) Spectral Characteristics of Kalman Filter Systems for Atmospheric Data Assimilation - Daley and Menard - 1993

Other topics: BibTeX, repository.

The presentation will be followed by a demonstration and discussion of strategies for transition to our new CVS repository.

The presentation will take about one hour, and will be followed by a work session on other topics.

This presentation will introduce the EAKF and LEKF and their implementation in progress here.

This work aims at improved computation of covariances and of multiscale structures such as clouds, in the Weather Research and Forecasting (WRF) data-assimilation (WRFDA) system, in particular the horizontal factor of the control-variable transform. Better representation can be achieved in the horizontal transform by wavelet-compression techniques that have been proven in many other applications. In this work, two past obstacles to effective incorporation of wavelets in limited-area models such as WRF are resolved: isometric-injective (i.e., energy preserving, left-invertible) wavelets avoid boundary-condition assumptions at any scale; and these wavelets can be applied to non-dyadic data lengths. A summary technical description of these improved wavelets and their implementation into WRFDA is presented. By retaining only a diagonal background-covariance matrix in wavelet space, appropriate heterogeneity is obtained for the model-space covariances. A second wavelet application is to partition observation error into a part due to poor representation (e.g., too-coarse resolution), and a residual, using a novel criterion in wavelet space. Other methods to construct inhomogeneous anisotropic covariance models are cited, and other potential technical improvements are discussed.

I will report on several presentations and discussion related to exascale computing at the IMA workshop "Numerical solution of PDE: fast solutions techniques" last week. Exascale means 10^18 floating point operations per second (Flops), and a computer capable of doing that is expected to be built by 2018. This is about 10^9 times more than a processor used in a typical laptop or PC. Current processors are limited to about 10Gflops (10^10 Flops) per core, and a handful of cores per node.

The need for higher parallelism and design constraints (communication, power consumption, failure tolerance) demand new methodologies and demise of existing ones, which changes will trickle down to everyday systems as well. Such emerging computer methodologies include the GPU and the cloud programming models, while the importance of the now prevalent MPI will diminish. Quantification of uncertainty and stochastic computing in general are examples of mathematical methodologies that will gain increased importance. The presentation will contain excerpts from the videos of IMA talks and discussions, as well as from various other materials.

See http://ccm.ucdenver.edu/wiki/Jan_Mandel/Blog/2010_Oct_Nov#Exascale_computing for links and other material.

I will give a dry run my two 15 min presentations the next day at the American Meteorological Society 91st Annual Meeting in Seattle:

Jan Mandel, Jonathan D. Beezley, and Loren Cobb, Spectral and morphing ensemble Kalman filters

Jan Mandel, Jonathan D. Beezley, and Adam K. Kochanski, An overview of the coupled atmosphere-wildland fire model WRF-Fire

The presentations will be followed by a discussion.

Abstract for the next Joint Statistical Meeting will be also discussed.

1. Wavelet and EnKF code demo and walkthrough (Jon) 2. Abstract for the upcoming AMS meeting (Ashok) 3. Notes from the AMS meeting in Seattle last week (Jan)

This will be a shorter meeting just to touch base on ongoing projects.

Slutsky's theorem states that if random elements u_n converge weakly to random element u, and v_n converge weakly to constant v (i.e., not random), then u_n+v_n converge to u+v, and similarly for product and quotient. The usual proof assumes that the random elements have values in finite dimensional vector space. However, there are proofs available that work in a metric space. Based on Van der Vaart's book, Wikipedia discussion and article, and Dudley's paper. Portmanteau theorem from Billingsley's book.

Please see http://ccm.ucdenver.edu/wiki/Jan_Mandel/Blog for links to the material.

Introduction following Da Prato's book, An Introduction to Infinite-Dimensional Analysis

I will discuss preliminary research on doing statistics, tracking, and filtering of data sets which are manifold valued or are functions of a manifold. We will work with Riemannian manifolds, which will give us many tools of analysis and calculus in order to extend our notions of statistics. I will define key terms, discuss some preliminary results, some ideas for research, and some of the work of other researchers. I will also discuss particular applications of statistics on Riemannian manifolds in forecasting and image analysis.

Weak convergence is a generalization of convergence in distribution of random variables to random elements with values in a metric space. Portmanteau theorem gives several equivalent conditions to weak convergence. Weak convergence is useful because it is preserved by continuous mappings, unlike some other types of convergence, such as convergence in L^p norm.. However, the standard calculus theorem on the limit of the sum does not completely carry over; only a somewhat restricted version (Slutsky's theorem) holds.

This is a introductory presentation with proofs. Please see http://ccm.ucdenver.edu/wiki/Jan_Mandel/Blog/2011_Apr_May for links.

This is the third part of the series of convergence theorems in probability to use in ongoing work on data assimilation in infinite dimensional spaces. This theory is generally presented in finite dimension. We are going over the proofs in detail and seek suitable variants to make sure that they still hold. See http://ccm.ucdenver.edu/wiki/Jan_Mandel/Blog

Vitalli convergence theorem assumes only uniform integrability, which is weaker than the inequality assumed in the dominated convergence theorem. We show some interesting relations between uniform integrability, convergence in probability, and in Lp. This is the fourth part of the series of convergence theorems in probability to use in ongoing work on data assimilation in infinite dimensional spaces. This theory is generally presented in finite dimension in the literature. We are going over the proofs in detail and seek suitable variants to make sure that they still hold for random elements with values in Banach spaces. See http://ccm.ucdenver.edu/wiki/Jan_Mandel/Blog for more information.

We'll revisit the proofs in our earlier paper (to appear) from the point of view of the previous seminars on probability measures. Most arguments carry over to infinitely dimensional state but some require additional assumptions.

This presentation will walk through the basic Matlab classes for data assimilation in the Osimorph project. The flexible object oriented framework simplifies changing methodologies and applications, and it is the cornerstone of our current effort towards wavelet-based ensemble methods with position correction. While the objects hold metadata, the actual data are stored in NetCDF files, and brought into memory only as needed. With recent advances in computer hardware and software, production-size problems can be handled in Matlab.

The presentation will be followed by an update by Ashok Krishnamurthy on the epidemic model paper.

Please note the change of day and time!

Jan: Building wrf_model class in Jonathan's object oriented framework

Volodymyr: Update on perimeter ignition

I will discuss some work done during my visit in Prague over the last month.

1. Presentation "Spectral and morphing ensemble Kalman filters and applications" at the 31st International Symposium on Forecasting

2. Morphing and assimilation of Opera radar data for precipitation in WRF

3. A new quick and dirty method - Wavelet optimal statistical interpolation. Inspired by Laplace operator-based FFT OSI.

4. Progress in Jon's Matlab classes for data assimilation.

For some of this, see also http://ccm.ucdenver.edu/wiki/Jan_Mandel/Blog

Followed by discussion and catch-up on projects.

This informal presentation is an outcome of my attempt to understand the paper "Fire propagation in a 2-D random medium" by Abinet, Searby, and Stauffer, J. Physique, 1986. This paper describes an early cell-based (a.k.a. automata or game) method, which is one of the basic approaches used for the modeling of both fires and epidemics. A matlab demo will be shown. I will try to compare the results of the demo with the heuristic analysis in the paper, using critical exponents.

Similar approaches are used in epidemics simulation. These models build on the percolation theory in physics, and I will review its basic concepts.

See http://ccm.ucdenver.edu/wiki/Jan_Mandel/Blog for the papers.

The presentation will be followed by a discussion of our ongoing projects.

Long-term criticality, particle filter schemes for cell models.

Topics: 1. Self-organizing polynomial networks for data mining (Jan Mandel) 2. Epidemics spread visualization in Google Earth (Brian Smith)

1. Jon's Prague visit next week 2. Current and upcoming proposals

Sometimes in statistics one comes across data sets which are manifold valued. For example, in measuring the orientation of an object in 3D, the set off all possible orientations lies on a sphere. The concepts of mean and variance of distributions on a manifold are easy to define, but can be costly to compute. Given a discrete set of points on a manifold, a Karcher mean is defined as a point which minimizes the sum of square distances to the data points (called a minimum of the variance function). Distances on a manifold are taken to be the length of the shortest path on the manifold between two points. The Karcher mean is not necessarily unique. For example, on a sphere, given a point at the north pole and a second point at the south pole, the Karcher means are the entire equator. However, under certain bounding conditions on the distribution of the data points, the Karcher mean does become unique and the variance function becomes convex. For example, for the sphere, if all of our data points lie north of the equator, then there exists just one Karcher mean point, also north of the equator. Computing the Karcher mean becomes a convex constrained nonlinear optimization problem in these cases. I will discuss formulations of the problem of finding a Karcher mean and discuss computational complexity, as well as theorems related to the uniqueness of the Karcher mean.

From Jon's presentation at the IDAACS'11 conference and visit at the Institute of Computer Science, Czech Academy of Sciences in Prague.

Sometimes in statistics one comes across data sets which are manifold valued. For example, in measuring the orientation of an object in 3D, the set off all possible orientations lies on a sphere. The concepts of mean and variance of distributions on a manifold are easy to define, but can be costly to compute. Given a discrete set of points on a manifold, a Karcher mean is defined as a point which minimizes the sum of square distances to the data points (called a minimum of the variance function). Distances on a manifold are taken to be the length of the shortest path on the manifold between two points. The Karcher mean is not necessarily unique. For example, on a sphere, given a point at the north pole and a second point at the south pole, the Karcher means are the entire equator. However, under certain bounding conditions on the distribution of the data points, the Karcher mean does become unique and the variance function becomes convex. For example, for the sphere, if all of our data points lie north of the equator, then there exists just one Karcher mean point, also north of the equator. Computing the Karcher mean becomes a convex constrained nonlinear optimization problem in these cases. I will discuss formulations of the problem of finding a Karcher mean and discuss computational complexity, as well as theorems related to the uniqueness of the Karcher mean.

I will give a summary and analysis of the preliminary results of the Berkeley Earth Surface Temperature Project, which has recently reported a complete reanalysis of global land surface temperatures, using kriging and other modern methods of interpolation. I will also suggest how their methods could be extended using Bayesian data assimilation.

This will be an informal discussion of the Bayes interpretation of Tikhonov regularization (known in statistics as "ridge regression") and its extension to infinite-dimensional spaces.

This is a work session covering the basics to reconnect with this project.

This will be a slow-paced tutorial covering the basics and with focus on getting the math right, following Livings et al. Physica D 2007 http://dx.doi.org/10.1016/j.physd.2008.01.005

I will also outline some ideas for the a probabilistic analysis and a convergence proof.

The ENKF attempts to predict the trajectory of a particle by combining uncertain measurements of the particle's state vector (x(t_2)) with uncertain information from a previous state of the state vector (x(t_1)). The particle is assumed to be a random variable undergoing an affine transformation from t_1 to t_2. We wish to extend this method of prediction to particles travelling in a Riemannian or even a more general manifold. It is not immediately clear how we can extend the notion of an affine transformation on a manifold, but it can be done. We consider the simplest case in which the particle undergoes a parallel translation. In R^n, this corresponds to x(t_2) = x(t_1) + w, where w is some vector in R^n. One can visualize this as a vector field over R^n, with all vectors being parallel and equal to w, and each point x(t_1) being translated along this vector field. In a Riemannian manifold, we can extend this notion of parallel translation. This is actually one version of the Fundamental Theorem of Riemannian Geometry. We will explain how this theorem relates to parallel translation. We will also derive transformation laws for the mean and covariance matrix of a random vector on a Riemannian Manifold. We will also explain the concept of tensors, which has many applications in analysis, linear algebra, and pde's.

This talk will present an updated and simplified version of the proof of EnKF convergence in the gaussian case from Mandel, Cobb, and Beezley http://dx.doi.org/10.1007/s10492-011-0031-2 . The proof relies on the fact that the ensemble is a finite exchangeable sequence of random elements. Generalizations and obstacles to a generalization will be discussed. While all seems to be good for an ensemble in a Hilbert space with finite dimensional data, extension to the case of infinitely dimensional data runs into a difficulty when the data perturbation is white noise. Square root ensemble filters do not involve a data perturbation, but the generated ensemble may not be exchangeable, except for the ensemble adjustment Kalman filter (EAKF), which does make exchangeable ensembles.