Innovation, Entropy and Exoplanets

I enjoy Shipulski on Design for the short articles on innovation. They are generally not technical at all. I like to think of most of the posts as innovation poetry to put your thoughts along the right lines of effort. This recent post has a huge, interesting technical iceberg riding under the surface though.

If you run an experiment where you are 100% sure of the outcome, your learning is zero. You already knew how it would go, so there was no need to run the experiment. The least costly experiment is the one you didn’t have to run, so don’t run experiments when you know how they’ll turn out. If you run an experiment where you are 0% sure of the outcome, your learning is zero. These experiments are like buying a lottery ticket – you learn the number you chose didn’t win, but you learned nothing about how to choose next week’s number. You’re down a dollar, but no smarter.

The learning ratio is maximized when energy is minimized (the simplest experiment is run) and probability the experimental results match your hypothesis (expectation) is 50%. In that way, half of the experiments confirm your hypothesis and the other half tell you why your hypothesis was off track.
Maximize The Learning Ratio

A One-Equation Local Correlation-Based Transition Model

This article is available for free download until 15 Oct 2015, h/t ANSYS blog.

Here's the Abstract:

A model for the prediction of laminar-turbulent transition processes was formulated. It is based on the LCTM (‘Local Correlation-based Transition Modelling’) concept, where experimental correlations are being integrated into standard convection-diffusion transport equations using local variables. The starting point for the model was the γ-Re θ model already widely used in aerodynamics and turbomachinery CFD applications. Some of the deficiencies of the γ-Re θ model, like the lack of Galilean invariance were removed. Furthermore, the Re θ equation was avoided and the correlations for transition onset prediction have been significantly simplified. The model has been calibrated against a wide range of Falkner-Skan flows and has been applied to a variety of test cases.

Keywords: Laminar-turbulent transition, Correlation, Local variables
Authors: Florian R. Menter, Pavel E. Smirnov , Tao Liu, Ravikanth Avancha

Transition location, and subsequent turbulence modeling remain the largest source of uncertainty for most engineering flows. Even for chemically reacting flows the source of uncertainty is often less the parameters and reactions for the chemistry, and more the uncertainty in the fluid state driven by shortcomings in turbulence and transition modeling.

Guidelines for Planning and Evidence for Assessing a Well-Designed Experiment

This paper is full of great guidance for planning a campaign of experimentation, or assessing the sufficiency of a plan that already exists. The authors break up the effort into four phases:

1. Plan a Series of Experiments to Accelerate Discovery
2. 1. Design Alternatives to Span the Factor Space
   2. Decide on a Design Strategy to Control the Risk of Wrong Conclusions
3. Execute the Test
4. Analyze the Experimental Design

They give a handy checklist for each phase (reproduced below). The checklists are comprehensive (significantly more than my little list of questions) and I think they stand-alone, but the whole paper is well worth a read. Design of experiments is more than just math, as this paper stresses it is a strategy for discovery.

Gaussian Processes for Machine Learning

What a great resource for learning about Gaussian Processes: The Gaussian Processes Web Site.

• The Gaussian Processes for Machine Learning book.
  
• Software around the web, and to go with the book
  
• Data Sets
• Tutorials
  

Validation & Verification in Physics of Plasmas

Physics of Plasmas is making a collection of 20 papers on verification and validation available for free download for a limited time.

Theoretical models, both analytical and numerical, are playing an increasingly important role in predicting complex plasma behavior, and providing a scientific understanding of the underlying physical processes.

Since the ability of a theoretical model to predict plasma behavior is a key measure of the model’s accuracy and its ability to advance scientific understanding, it is Physics of Plasmas’ Editorial Policy to encourage the submission of manuscripts whose primary focus is the verification and/or validation of codes and analytical models aimed at predicting plasma behavior.

Phil Roe: Colorful Fluid Dynamics

Echos of Tufte in one of his introductory statements: "It's full of noise, it's full of color, it's spectacular, it's intended to blow your mind away, it's intended to disarm criticism." And further on the dangers of "colorful fluid dynamics":

These days it is common to see a complicated flow field, predicted with all the right general features and displayed in glorious detail that looks like the real thing. Results viewed in this way take on an air of authority out of proportion to their accuracy.
--Doug McLean

This lecture is sponsored by MConneX.

Roe wraps up the lecture by referencing a NASA sponsored study, CFD Vision 2030, that addresses whether CFD will be able to reliably predict turbulent separated flows by 2030. The conclusion is that advances in hardware capability alone will not be enough, but that significant improvements in numerical algorithms are required.

Starcraft, Jaynes, and Bayesian UQ

Interesting content covered on Nuite Blanche of a recent Paris Machine Learning meetup. The work on applying Bayesian Programming and Learning for Multi-Player Video Games was really neat. It's about developing a bot for playing Starcraft. Some additional related presentations:

• A Bayesian Model for Plan Recognition in RTS Games applied to StarCraft
  
• A Bayesian Model for RTS Units Control applied to StarCraft
  

This problem is real enough to be interesting, and simple enough to tackle. I really liked how this work dealt with the separation of strategic/tactical level decision making.

11th World Congress on Computational Mechanics (WCCM XI)

The 11th World Congress on Computational Mechanics (WCCM XI) is coming up in July in Barcelona. Dates: abstracts due 29 Nov 2013, and full papers due 28 Feb 2014. The deadline for mini-symposium proposals has passed. There are some interesting looking ones (ECCM, ECFD),

• Uncertainty Quantification Techniques for Fluid-flow Problems
  
• Verification & Validation
  
• Frontiers of Numerical Uncertainty Quantification in Computational Mechanics
  
• New Trends in Topology Optimization
  

Environmental Decisions in the Face of Uncertainty

New report from the National Academies Press on decision making under uncertainty.

Description: The U.S. Environmental Protection Agency (EPA) is one of several federal agencies responsible for protecting Americans against significant risks to human health and the environment. As part of that mission, EPA estimates the nature, magnitude, and likelihood of risks to human health and the environment; identifies the potential regulatory actions that will mitigate those risks and protect public health1 and the environment; and uses that information to decide on appropriate regulatory action. Uncertainties, both qualitative and quantitative, in the data and analyses on which these decisions are based enter into the process at each step. As a result, the informed identification and use of the uncertainties inherent in the process is an essential feature of environmental decision making.

National Strategy for Advancing Climate Modeling

NAP has a new report out on a national strategy for advancing climate modeling. I've just done a quick skim so far. Some tidbits below that touch on V&V and uncertainty.

Validating the Prediction of Unobserved Quantities

I was going to just post a bit of the abstract and a link to this recent tech report as a comment on my V&V workshop notes, but Validating the Prediction of Unobserved Quantities was interesting enough to deserve its own post. The basic approach is well-founded in probability theory, and there are some novel concepts as well.

Here's the abstract:

In predictive science, computational models are used to make predictions regarding the response of complex systems. Generally, there is no observational data for the predicted quantities (the quantities of interest or QoIs) prior to the computation, since otherwise predictions would not be necessary. Further, to maximize the utility of the predictions it is necessary to assess their reliability|i.e., to provide a quantitative characterization of the discrepancies between the prediction and the real world. Two aspects of this reliability assessment are judging the credibility of the prediction process and characterizing the uncertainty in the predicted quantities. These processes are commonly referred to as validation and uncertainty quantification (VUQ), and they are intimately linked. In typical VUQ approaches, model outputs for observed quantities are compared to experimental observations to test for consistency. While this consistency is necessary, it is not sufficient for extrapolative predictions because, by itself, it only ensures that the model can predict the observed quantities in the observed scenarios. Indeed, the fundamental challenge of predictive science is to make credible predictions with quantified uncertainties, despite the fact that the predictions are extrapolative. At the PECOS Center, a broadly applicable approach to VUQ for prediction of unobserved quantities has evolved. The approach incorporates stochastic modeling, calibration, validation, and predictive assessment phases where uncertainty representations are built, informed, and tested. This process is the subject of the current report, as well as several research issues that need to be addressed to make it applicable in practical problems.

VV&UQ for Historic Masonry Structures

What a neat application of verification, validation and uncertainty quantification (VV&UQ) methods! The paper is Uncertainty quantification in model verification and validation as applied to large scale historic masonry monuments.

Abstract: This publication focuses on the Verification and Validation (V&V) of numerical models for establishing confidence in model predictions, and demonstrates the complete process through a case study application completed on the Washington National Cathedral masonry vaults. The goal herein is to understand where modeling errors and uncertainty originate from, and obtain model predictions that are statistically consistent with their respective measurements. The approach presented in this manuscript is comprehensive, as it considers all major sources of errors and uncertainty that originate from numerical solutions of differential equations (numerical uncertainty), imprecise model input parameter values (parameter uncertainty), incomplete definitions of underlying physics due to assumptions and idealizations (bias error) and variability in measurements (experimental uncertainty). The experimental evidence necessary for reducing the uncertainty in model predictions is obtained through in situ vibration measurements conducted on the masonry vaults of Washington National Cathedral. By deploying the prescribed method, uncertainty in model predictions is reduced by approximately two thirds.

Highlights:

• Developed a finite element model of Washington National Cathedral masonry vaults.
  
• Carried out code and solution verification to address numerical uncertainties.
  
• Conducted in situ vibration experiments to identify modal parameters of the vaults.
  
• Calibrated and validated model to mitigate parameter uncertainty and systematic bias.
  
• Demonstrated a two thirds reduction in the prediction uncertainty through V&V.
  

Keywords: Gothic Cathedral; Modal analysis; Finite element modeling; Model updating; Bayesian inference; Uncertainty quantification

I haven't read the full-text yet, but it looks like a coherent (Bayesian) and pragmatic approach to the problem.

Notre Dame V&V Workshop

The folks at Notre Dame are putting on a Verification and Validation Workshop. The preliminary agenda is up on the workshop site.

The purpose of the workshop is to bring together a diverse group of computational scientists working in fields in which reliability of predictive computational models is important. Via formal presentations, structured discussions, and informal conversations, we seek to heighten awareness of the importance of reliable computations, which are becoming ever more critical in our world.

The intended audience is computational scientists and decision makers in fields as diverse as earth/atmospheric sciences, computational biology, engineering science, applied mechanics, applied mathematics, astrophysics, and computational chemistry.

It looks very interesting. Unfortunately, I don't think I'll be able to attend. Hopefully they will post posters and papers. Update: I was able to attend. They will be posting the presenters' slides. I will put up some of my notes when they get the slides up.

Empiricism and Simulation

There are two orthogonal ideas that seem to get conflated in discussions about climate modeling. One is the idea that you’re not doing science if you can’t do a controlled experiment, but of course we have observational sciences like astronomy. The other is that all this new-fangled computer-based simulation is untrustworthy, usually because “it ain’t the way my grandaddy did science.” Both are rather silly ideas. We can still weigh the evidence for competing models based on observation, and we can still find protection from fooling ourselves even when those models are complex.

What does it mean to be an experimental as opposed to an observational science? Do sensitivity studies, and observational diagnostics using sophisticated simulations count as experiments? Easterbrook claims that because climate scientists do these two things with their models that climate science is an experimental science [1]. It seems like there is a motivation to claim the mantle of experimental, because it may carry more rhetorical credibility than the merely observational (the critic Easterbrook is addressing certainly thinks so). This is probably because the statements we can make about causality and the strength of the inferences we can draw are usually greater when we can run controlled experiments than when we are stuck with whatever natural experiments fortune provisions for us (and there are sound mathematical reasons for this, having to do with optimality in experimental design rather than any label we may place on the source of the data). This seeming motivation demonstrated by Easterbrook to embrace the label of empirical is in sharp contrast to the denigration of the empirical by Tobis in his three part series [2, 3, 4]. As I noted on his site, the narrative Tobis is trying to create with those posts has already been pre-messed with by Easterbrook, his readers just pointed out the obvious weaknesses too. One good thing about blogging is the critical and timely feedback.

The confusions of these two climate warriors are an interesting point of departure. I think they are both saying more than blah blah blah, so it’s worth trying to clarify this issue. The figure below is based on a technical report from Sandia [5], which is a good overview and description of the concepts and definitions for model verification and validation as it has developed in the computational physics community over the past decade or so. I think this emerging body of work on model V&V places the relative parts, experiment and simulation, in a sound framework for decision making and reasoning about what models mean.

---

Figure 1: Verification and Validation Process (based largely on [5])

---

The process starts at the top of the flowchart with a “Reality of Interest”, from which a conceptual model is developed. At this point the path splits into two main branches. One based on “Physical Modeling” and the other based on “Mathematical Modeling”. Something I don’t think many people realize is that there is a significant tradition of modeling in science that isn’t based on equations. It is no coincidence that an aeronautical engineer might talk of testing ideas with a wind-tunnel model or a CFD model. Both models are simplifications of the reality of interest, which, for that engineer, is usually a full-scale vehicle in free flight.

Figure 2 is just a look at the V&V process through my Design of Experiments (DoE) colored glasses.

---

Figure 2: Distorted Verification and Validation Process

---

My distorted view of the V&V process is shown to emphasize that there’s plenty of room for experimentalists to have fun (maybe even a job [3]) in this, admittedly model-centric, sandbox. However, the transferability of the basic experimental design skills between “Validation Experiments” and “Computational Experiments” says nothing about what category of science one is practicing. The method of developing models may very well be empirical (and I think Professor Easterbrook and I would agree it is, and maybe even should be), but that changes nothing about the source of the data which is used for “Model Validation.”

The computational experiments highlighted in Figure 2 are for correctness checking, but those aren’t the sorts of computational experiments Easterbrook claimed made climate science an experimental science. Where do sensitivity studies and model-based diagnostics fit on the flowchart? I think sensitivity studies fit well in the activity called “Pre-test Calculations”, which, one would hope, inform the design of experimental campaigns. Diagnostics are more complicated.

Heald and Wharton have a good explanation for the use of the term “diagnostic” in their book on microwave-based plasma diagnostics: “The term ‘diagnostics,’ of course, comes from the medical profession. The word was first borrowed by scientists engaged in testing nuclear explosions about 15 years ago [c. 1950] to describe measurements in which they deduced the progress of various physical processes from the observable external symptoms” [6]. With a diagnostic we are using the model to help us generate our “Experimental Data”, so that would happen within the activity of “Experimentation” on this flowchart. This use of models as diagnostic tools is applied to data obtained from either experiment (e.g. laboratory plasma diagnostics) or observations (e.g. astronomy, climate science), so it says nothing about whether a particular science is observational or experimental. Classifying scientific activities as experimental or observational is of passing interest, but I think far too much emphasis is placed on this question for the purpose of winning rhetorical “points.”

The more interesting issue from a V&V perspective is introducing a new connection in the flowchart that shows how a dependency between model and experimental data could exist (Figure 3). Most of the time the diagnostic model, and the model being validated are different. However, this case where they are the same is an interesting and practically relevant one that is not addressed in the current V&V literature that I know of (please share links if you “know of”).

---

Figure 3: V&V process including model-based diagnostic

---

It should be noted that even though the same model may be used to make predictions and perform diagnostics, it will usually be run in a different way for those two uses. The significant changes between Figure 1 and Figure 3 are the addition of a “Experimental Diagnostic” box and the change to the mathematical cartoon in the “Validation Experiment” box. The change to the cartoon is to indicate that we can’t measure what we want directly (u), so we have to use a diagnostic model to estimate it based on the things we can measure (b). An example of when the model-based diagnostic is relatively independent of the model being validated might be using laser-based diagnostic for fluid flow. The equations describing propagation of the laser through the fluid are not the same as those describing the flow. An example of when the two codes might be connected would be if you were trying to use ultrasound to diagnose a flow. The diagnostic model and the predictive model could both be Navier-Stokes with turbulence closures. Establishing the validity of which is the aim of the investigation. I’d be interested in criticisms of how I explained this / charted this out.

Afterward

Attempt at Answering Model Questions

I’m not in the target population that professor Easterbrook is studying, but here’s my attempt at answering his questions about model validation[7].

1. “If I understand correctly–a model is ’valid’ (is that a formal term?) if the code is written to correctly represent the best theoretical science at the time...”

   I think you are using an STS flavored definition for “valid.” The IEEE/AIAA/ASME/US-DoE/US-DoD definition differs. “Valid” means observables you get out of your simulations are “close enough” to observables in the wild (experimental results). The folks from DoE tend to argue for a broader definition of valid than the DoD folks. They’d like to include as “validation” activities of a scientist comparing simulation results and experimental results without reference to an intended use.

2. “– so then what do the results tell you? What are you modeling for–or what are the possible results or output of the model?”

   Doing a simulation (running the implementation of a model) makes explicit the knowledge implicit in your modeling choices. The model is just the governing equations, you have to run a simulation to find solutions to those governing equations.

3. “If the model tells you something you weren’t expecting, does that mean it’s invalid? When would you get a result or output that conflicts with theory and then assess whether the theory needs to be reconsidered?”

   This question doesn’t make sense to me. How could you get a model output that conflicted with theory? The model is based on theory. Maybe this question is about how simplifying assumptions could lead to spurious results? For example, if a simulation result shows failure to conserve mass/momentum/energy in a specific calculation possibly due to a modeling assumption (more likely due to a more mundane error), I don’t think anyone but a perpetual-motion machine nutter would seriously reconsider the conservation laws.

4. “Then is it the theory and not the model that is the best tool for understanding what will happen in the future? Is the best we can say about what will happen that we have a theory that adheres to what we know about the field and that makes sense based on that knowledge?”

   This one doesn’t make sense to me either. You have a “theory,” but you can’t formulate a “model” of it and run a simulation, or just a pencil and paper calculation? I don’t think I’m understanding how you are using those words.

5. “What then is the protection or assurance that the theory is accurate? How can one ‘check’ predictions without simply waiting to see if they come true or not come true?”

   There’s no magic; the protection from fooling ourselves is the same as it has always been, only the names of the problems change.

Attempt at Understanding Blah Blah Blah

• “The trouble comes when empiricism is combined with a hypothesis that the climate is stationary, which is implicit in how many of their analyses work.” [8]

  The irony of this statement is extraordinary in light of all the criticisms by the auditors and others of statistical methods in climate science. It would be a valid criticism, if it were supported.

• “The empiricist view has never entirely faded from climatology, as, I think, we see from Curry. But it’s essentially useless in examining climate change. Under its precepts, the only thing that is predictable is stasis. Once things start changing, empirical science closes the books and goes home. At that point you need to bring some physics into your reasoning.” [2]

  So we’ve gone from what could be reasonable criticism of unfounded assumptions of stationarity to empiricism being unable to explain or understand dynamics. I guess the guys working on embedding dimension stuff, or analogy based predictions would be interested to know that.

• “See, empiricism lacks consilience. When the science moves in a particular direction, they have nothing to offer. They can only read their tea leaves. Empiricists live in a world which is all correlation, and no causation.” [3]

  Lets try some definitions.

  empiricism

  knowledge through observation

  consilience

  unity of knowledge, non-contradiction

  How can the observations contradict each other? Maybe a particular explanation for a set of observations is not consilient with another explanation for a different set of observations. This seems to be something that would get straightened out in short order though: it’s on this frontier that scientific work proceeds. I’m not sure how empiricism is “all correlation.” This is just a bald assertion with no support.

• “While empiricism is an insufficient model for science, while not everything reduces to statistics, empiricism offers cover for a certain kind of pseudo-scientific denialism. [...] This is Watts Up technique asea; the measurements are uncertain; therefore they might as well not exist; therefore there is no cause for concern!” [4]

  Tobis: Empiricism is an insufficient model for science. Feynman: The test of all knowledge is experiment. Tobis: Not everything reduces to statistics. Jaynes: Probability theory is the logic of science. To be fair, Feynman does go on to say that you need imagination to think up things to test in your experiments, but I’m not sure that isn’t included in empiricism. Maybe it isn’t included in the empiricism Tobis is talking about.

  So that’s what all this is about? You’re upset at Watts making a fallacious argument about uncertainty? What does empiricism have to do with this? It would be simple enough to just point out that uncertainty doesn’t mean ignorance.

Not quite blah blah blah, but the argument is still hardly thought out and poorly supported.

References

[1]   Easterbrook, S., “Climate Science is an Experimental Science,”http://www.easterbrook.ca/steve/?p=1322, February 2010.

[2]   Tobis, M., “The Empiricist Fallacy,” http://initforthegold.blogspot.com/2010/11/empiricist-fallacy.html, November 2010.

[3]   Tobis, M., “Empiricism as a Job,”http://initforthegold.blogspot.com/2010/11/empiricism-as-job.html, November 2010.

[4]   Tobis, M., “Pseudo-Empiricism and Denialism,”http://initforthegold.blogspot.com/2010/11/pseudo-empiricism-and-denialism.html, November 2010.

[5]   Thacker, B. H., Doebling, S. W., Hemez, F. M., Anderson, M. C., Pepin, J. E., and Rodriguez, E. A., “Concepts of Model Verification and Validation,” Tech. Rep. LA-14167-MS, Los Alamos National Laboratory, Oct 2004.

[6]   Heald, M. and Wharton, C., Plasma Diagnostics with Microwaves, Wiley series in plasma physics, Wiley, New York, 1965.

[7]   Easterbrook, S., “Validating Climate Models,”http://www.easterbrook.ca/steve/?p=2032, November 2010.

[8]   Tobis, M., “Empiricism,”http://initforthegold.blogspot.com/2010/11/empiricism.html, November 2010.

Thanks to George Crews and Dan Hughes for their critical feedback on portions of this.

[Update: George left a comment with suggestions on changing the flowchart. Here's my take on his suggested changes.

A slightly modified version of George's chart. I think it makes more sense to have the "No" branch of the validation decision point back at "Abstraction", which parallels the "No" branch of the verification decision pointing at "Implementation". Also switched around "Experimental Data" and "Experimental Diagnostic." Notably absent is any loop for "Calibration"; this would properly be a separate loop with output feeding in to "Computer Model."

]

Mathematical Science Foundations of Validation, Verification, and Uncertainty Quantification

There is an interesting study being conducted by the National Academies.  Here's their "Statement of Task":

• A committee of the NRC will examine practices for verification and validation (V&V) and uncertainty quantification (UQ) of large-scale computational simulations in several research communities.
• Identify common concepts, terms, approaches, tools, and best practices of V&V and UQ.
• Identify mathematical sciences research needed to establish a foundation for building a science of V&V and for improving practice of V&V and UQ.
• Recommend educational changes needed in the mathematical sciences community and mathematical sciences education needed by other scientific communities to most effectively use V&V and UQ.

Here's a list of the folks on the committee.  It should be interesting to see the study results (it's an 18 month long effort).

Zen Uncertainty

Zen Uncertainty: Attempts to understand uncertainty are mere illusions; there is only suffering.
-- WARNING: Physics Envy May Be Hazardous To Your Wealth!

Should we give up? No, there's plenty we can do to make the suffering more bearable. Lo and Mueller give an uncertainty taxonomy of five levels in their 'Physics Envy' paper:

1. Complete Certainty: the idealized deterministic world
2. Risk without Uncertainty: an honest casino
3. Fully Reducible Uncertainty: the odds in the honest casino are not posted, we have to learn them from limited experience
4. Partially Reducible Uncertainty: we're not quite sure which game at the casino we're playing so we have to learn that as well as the odds based on limited experience
5. Irreducible Uncertainty: we're not even sure if we're in the casino, we might be outside splashing around in the fountain...

At the bottom of the decent we find level infinity, Zen Uncertainty.

Section 2 of the paper provides a nice historical overview of the early work of Paul A. Samuelson, who single-handedly brought statistical mechanics to the economists, and they have never been the same since. Samuelson acknowledged the deep connection between his work and physics:

Perhaps most relevant of all for the genesis of Foundations, Edwin Bidwell Wil- son (1879–1964) was at Harvard. Wilson was the great Willard Gibbs’s last (and, essentially only) protege at Yale. He was a mathematician, a mathematical physicist, a mathematical statistician, a mathematical economist, a polymath who had done first-class work in many fields of the natural and social sciences. I was perhaps his only disciple . . . I was vaccinated early to understand that economics and physics could share the same formal mathematical theorems (Euler’s theorem on homogeneous functions, Weierstrass’s theorems on constrained maxima, Jacobi determinant identities underlying Le Chatelier reactions, etc.), while still not resting on the same empirical foundations and certainties.

Related to this theme, there's an interesting recent article over on Mobjectivist site about using ideas from physics to model income distributions.

Lo and Mueller propose to operationalize their uncertainty taxonomy with a 2-D checklist (table). The levels provide the columns across the top, and there is a row for each business component of the activity being evaluated, here's their description:

The idea of an uncertainty checklist is straightforward: it is organized as a table whose columns correspond to the five levels of uncertainty of Section 3, and whose rows correspond to all the business components of the activity under consideration. Each entry consists of all aspects of that business component falling into the particular level of uncertainty, and ideally, the individuals and policies responsible for addressing their proper execution and potential failings.

This seems like an idea that could be adapted and combined with best practices for model validation (and checklist sorts of approaches) in helping to define what sorts of uncertainties we are operating under when we make decisions using science-based decision support products.

Their final paragraph echos Lindzen's sentiments about climate science:

While physicists have historically been inspired by mathematical elegance and driven by pure logic, they also rely on the ongoing dialogue between theoretical ideals and experimental evidence. This rational, incremental, and sometimes painstaking debate between idealized quantitative models and harsh empirical realities has led to many breakthroughs in physics, and provides a clear guide for the role and limitations of quantitative methods in financial markets, and the future of finance.
-- WARNING: Physics Envy May Be Hazardous To Your Wealth!

Uncertain Rate in FFT-based Oil Extraction Model

This is an extension to the FFT-based oil extraction model (see the Mobjectivist blog for more more details). The basic approach remains the same, each phase of the process is modeled by a convolution of the rate in the previous phase with an exponential decay. Now we are going to apply some of the ideas I presented in Uncertainty Quantification with Stochastic Collocation.

Here is the original Python function which applied the convolutions using an FFT,

 
def exponential_convolution(x, y, r): 
    """Convolve␣an␣exponential␣function␣e^(-r*x)␣with␣y,␣uses␣FFT.""" 
    expo = sp.exp(-r*x) 
    expo = expo / sp.sum(expo) # normalize 
    return(ifft(fft(expo) * fft(y)).real)

and here is the function modified to use the complex step derivative calculation method:

 
def expo_conv_csd(x, y, r): 
    """Convolve␣an␣exponential␣function␣e^(-r*x)␣with␣y,␣uses␣FFT.␣Use 
␣␣␣␣the␣complex␣step␣derivative␣method␣to␣estimate␣the␣slope␣with 
␣␣␣␣respect␣to␣the␣rate. 
 
␣␣␣␣See␣Cervino␣and␣Bewley,␣’On␣the␣extension␣of␣the␣complex-step 
␣␣␣␣derivative␣technique␣to␣pseudospectral␣algorithms’, 
␣␣␣␣J.Comp.Phys.␣187␣(2003). 
␣␣␣␣""" 
    expo = sp.exp(-r*x) 
    # normalize: 
    expo = expo / sp.sum(expo) 
    # need to transform the real and imag parts seperately to avoid 
    # problems due to subtractive cancellation: 
    expo.real = ifft(fft(expo.real) * fft(y)).real 
    expo.imag = ifft(fft(expo.imag) * fft(y)).real 
    return(expo)

Notice that we need to transform the real and imaginary parts of our solution separately so that the nice subtractive-cancellation-avoidance properties of the method are retained (see [1] for a more detailed discussion).

Now we have the result of the convolution in the real part of the return value and the sensitivity to changes in rate in the imaginary part (see Figure 1).

---

Figure 1: Rate sensitivity with a normally distributed input and normally distributed rate

---

Now we proceed as shown before in the UQ post, except this time we have a separate probability density for the result at each time. Figure 2 shows the mean prediction as well as shading based on confidence intervals (at the 0.3, 0.6, and 0.9 level). The Python to generate this figure is shown below, it uses the alpha parameter (transparency) available in the matplotlib plotting package in a similar way to the vizualization approach shown in this post.

 
p.figure() 
p.plot(x, y, label="input") 
p.plot(x, y_mean, ’g’, label="mean␣output") 
p.fill_between(x, y_ci30[0], y_ci30[1], where=None, color=’g’, alpha=0.2) 
p.fill_between(x, y_ci60[0], y_ci60[1], where=None, color=’g’, alpha=0.2) 
p.fill_between(x, y_ci90[0], y_ci90[1], where=None, color=’g’, alpha=0.2) 
p.legend(loc=0) 
p.xlabel("time") 
p.ylabel("rate") 
p.savefig("uncertain_rate.png")

---

Figure 2: Results with an uncertain rate

---

Compare the result in Figure 2 with this Monte Carlo analysis. It’s not quite apples-to-apples, because that Monte Carlo includes uncertain discovery (the input) as well. Also, this result is just a first-order collocation, so it assumes that the rate sensitivity is constant with variations in the result. This is not actually the case, and you can see if you look closely that the 0.9-interval actually dips into negative territory, which is unphysical. This is a good example of the sort of common-sense checking Hamming suggested in his “N+1” essay as a requirement for any successful computation: Are the known conservation laws obeyed by the result? [2]. Clearly we are not going to “overshoot” on extraction and begin pumping oil back into the ground.

On a somewhat related note, Hamming had another insightful thing to say about the importance of checking the correctness of a calculation: It is the experience of the author that a good theoretician can account for almost anything produced, right or wrong, or at least he can waste a lot of time worrying about whether it is right or wrong [2]. Don’t worry,verify!

References

[1]   Cervino, L.I., Bewley, T.R., “On the extension of the complex-step derivative technique to pseudospectral algorithms,” Journal of Computational Physics, 187, 544-549, 2003.

[2]   Hamming, R.W., “Numerical Methods for Scientists and Engineers,” 2nd ed., Dover Publications, 1986.

(If you are serious about the art of scientific computing, that Dover edition of Hamming’s book is the best investment you could make with a very few bucks.)

Uncertainty Quantification with Stochastic Collocation

I’ve been posting lots of links to uncertainty quantification (UQ) references lately in comments (eg here, here, here), but I don’t really have a post dedicated only to UQ. So I figured I’d remedy that situaton. This post will be a simple example problem based on methods presented in some useful recent papers on UQ (along with a little twist of my own):

• A good overview of several of the related UQ methods [1]
• Discussion of some of the sampling approaches for high dimensional random spaces, along with specifics about the method for CFD [2]
• Comparison of speed-up over naive Monte Carlo sampling approaches [3]
• A worked example of stochastic collocation for a simple ocean model [4]

The example will be based on a nonlinear function of a single variable

(1)

The first step is to define a probability distribution for our input (Gaussian in this case, shown in Figure 1).

(2)

---

Figure 1: Uncertain Input Distribution

---

Then, since our example is a simple function, we can analytically calculate the resulting probability density for the output. To do that we need the derivative of the model with respect to the input (in the multi-variate case we’d need the Jacobian).

(3)

Then we divide the input density by the slope to get the new output density.

(4)

We’ll use equation 4 to measure the convergence of our collocation method.

The idea of stochastic collocation methods is that the points which the model (equation 1 in our example) is evaluated at are chosen so that they are orthogonal with the probability distributions on the inputs as a weighting function. For normally distributed inputs this results in choosing points at the roots of the Hermite polynomials. Luckily, Scipy has a nice set of special functions, which includes he_roots to give us the roots of the Hermite polynomials. What we really want is not the function evaluation at the collocation points, but the slope at the collocation points. One way to get this information in a nearly non-intrusive way is to use a complex step method (this is that little twist I mentioned). So, once we have values for this slope at the collocation points we fit a polynomial to those points and use this polynomial as a surrogate for equation 3. Figure 2 shows the results for several orders of surrogate slopes.

---

Figure 2: Stochastic Collocation Estimates of Output Uncertainty

---

The convergence of the method is shown in Figure 3. For low dimensional random spaces (only a few random inputs) this method will be much more efficient than random sampling. Eventually the curse of dimensionality takes its toll though, and random sampling becomes the faster approach.

---

Figure 3: Convergence of the Collocation Method

---

The cool thing about this method is that the complex step method which I used to estimate the Jacobian is also useful for sensitivity quantification which can be used for design optimization as well. On a related note, Dan Hughes has a post describing an interesting approach to bounding uncertainties in the output of a calculation given bounds on the inputs by using interval arithmetic.

References

[1]   Loeven, G.J.A., Witteveen, J.A.S., Bijl, H., Probabilistic Collocation: An Efficient Non-Intrusive Approach For Arbitrarily Distributed Parametric Uncertainties, 45th Aerospace Sciences Meeting, Reno, NV, 2007.

[2]   Hosder, S., Walters, R.W., Non-Intrusive Polynomial Chaos Methods for Uncertainty Quantification in Fluid Dynamics, 48th Aerospace Sciences Meeting, Orlando, FL, 2010.

[3]   Xiu, D., Fast Numerical Methods for Stochastic Computations, Comm. Comp. Phys., 2009.

[4]   Webster, M., Tatang, M.A., McRae, G.J., Application of the Probabilistic Collocation Method for an Uncertainty Analysis of a Simple Ocean Model, MIT JPSPGC Report 4, 2000.

Visualizing Confidence Intervals

This post is about visualizing confidence intervals. Matplotlib has some really neat capabilities that are useful in that regard, but before we get into the pictures, here’s a short digression on statistics.

I like resampling-based statistics a lot, as an engineer their practicality and intuitiveness appeals to me. It has been shown that students learn better and enjoy statistics more when they are taught using resampling methods. Here’s a nice description of the motivation for resampling (ok, maybe it's a bit over the top):

For more than a century the inherent difficulty of formula-based inferential statistics has baffled scientists, induced errors in research, and caused million of students to hate the subject.

Complexity is the disease. Resampling (drawing repeated samples from the given data, or population suggested by the data) is a proven cure. Bootstrap, permutation, and other computer-intensive procedures have revolutionized statistics. Resampling is now the method of choice for confidence limits, hypothesis tests, and other everyday inferential problems.

In place of the formidable formulas and mysterious tables of parametric and non-parametric tests based on complicated mathematics and arcane approximations, the basic resampling tools are simulations, created especially for the task at hand by practitioners who completely understand what they are doing and why they are doing it. Resampling lets you analyze most sorts of data, even those that cannot be analyzed with formulas.

– Resampling Stats

One of the useful things to do is resampling of residuals. The assumptions underlying most models are that the magnitude and sign of the residuals are not a function of the independent variable. This is the hypothesis which we’ll base our resampling on. First we fit a model (a line in these examples) to some data, then calculate all the residuals (difference between the data and the model). Then we can apply a couple of different resampling approaches towards understanding the confidence intervals.

A permutation-based way of establishing confidence intervals is easily accomplished in Python using the itertools module’s permutation function. This is an exact method, but the number of permutations grows as the factorial of the sample size. The six-sample example shown in figure 1 has 6! = 720 possible permutations of the residuals. With only ten samples the number of permutations grows to 3628800 (over three million). The point of this post is visualization, so plotting three million lines may not be worth our while.

---

Figure 1: Residual Permutations

---

Figure 1 shows the least-squares-fit line in solid black, with the lines fit using the permuted residuals slightly transparent. Here’s the Python to accomplish that:

 
nx = 6 
x = sp.linspace(0.0, 1.0, nx) 
data = x + norm.rvs(loc=0.0, scale=0.1, size=nx) 
yp = sp.polyfit(x, data, 1) 
y = sp.polyval(yp,x) 
r = data - y 
 
nperm = factorial(nx) 
for perm in permutations(r,nx): # loop over all the permutations of the resids 
    pc = sp.polyfit(x, y + perm, 1) 
    p.plot(x, sp.polyval(pc,x), ’k-’, linewidth=2, alpha=2.0/float(nperm)) 
p.plot(x, y, ’k-’) 
p.plot(x, data, ’ko’)

We’ve used the alpha argument to set the level of transparency in each of the permutations. Where more of the lines overlap, the color is darker. This gives a somewhat intuitive display of the ’weight’ of the variation in the model fit (confidence).

As mentioned above, doing exact permutation-based statistics becomes computationally prohibitive rather quickly, so we have to go to the random-sampling based approximations. The bootstrap is one such approach.

Rather than generating all possible permutations of the residuals, the bootstrap method involves drawing random samples (with replacement) from the residuals. This is done in Python easily enough by using the random_integer function to index into our array of residuals.

 
import scipy as sp

 
bootindex = sp.random.random_integers

 
nboot = 400 
for i in xrange(nboot): # loop over n bootstrap samples from the resids 
    pc = sp.polyfit(x, y + r[bootindex(0, len(r)-1, len(r))], 1) 
    p.plot(x, sp.polyval(pc,x), ’k-’, linewidth=2, alpha=3.0/float(nboot)) 
p.plot(x, y, ’k-’) 
p.plot(x, data, ’ko’)

Figure 2 shows the method applied to the same six data points as with the permutation.

---

Figure 2: Residual Bootstraps

---

Since the sampling in the bootstrap approach is with replacement we can get sets of residuals that have repeated values. This tends to introduce a bit of bias in the direction of that repeated residual. You can see that behaviour exhibited in Figure 2 by the spreading of the lines around the middle of the graph, in contrast to the tight intersection in the middle of Figure 1.

Of course the reason the bootstrap method is useful is we often have large samples that are impractical to do with permutations. The Python to generate the 12-sample example in Figure 3 is shown below.

 
nx = 12 
x = sp.linspace(0.0, 1.0, nx) 
data = x + norm.rvs(loc=0.0, scale=0.1, size=nx) 
yp = sp.polyfit(x, data, 1) 
y = sp.polyval(yp,x) 
r = data - y 
p.figure() 
for i in xrange(nboot): # loop over n bootstrap samples from the resids 
    pc = sp.polyfit(x, y + r[bootindex(0, len(r)-1, len(r))], 1) 
    p.plot(x, sp.polyval(pc,x), ’k-’, linewidth=2, alpha=3.0/float(nboot)) 
p.plot(x, y, ’k-’) 
p.plot(x, data, ’ko’)

---

Figure 3: Residual Bootstraps, large sample

---

These sorts of graphs probably won’t replace the standard sorts of confidence intervals (see the graph in this post for example), but it’s a kind of neat way of looking at things, and a good demo of some of the cool stuff you can do really easily in Python with Matplotlib and Scipy.