Easy Full-Factorial Ensemble Experiments

This is another installment of the Lorenz63 series. In this post I wanted to demonstrate an ’experiment’ on the Lorenz ’63 system. The basics of these sorts of computer experiments as they relate to climate modeling are covered at climateprediction.net.

The main factors underlying our experimental design are initial conditions, parameters and forcing. Since we followed Roache’s advice about making our code Verifiable with the method of manufactured solutions (MMS) the ability to specify forcings is already included. Python has some nice built-in functions in the itertools module that make doing a computational experiment really easy. First all of the arguments for the time integrator are defined as lists (the ones which aren’t part of the experiment only have one element, so they remain constant):

# factors in our experimental design (all of the inputs to the system 
# integrator, some only have one level): 
nt = 2000 
savevery = [10] 
nsaves = nt / savevery[0] + 1 
T = 2.0 
# the first element of the array is the initial condition, the rest 
# are over-written with the solution: 
x = [1.05*sp.ones(nsaves,dtype=float,order=’F’), 
     0.95*sp.ones(nsaves,dtype=float,order=’F’)] 
y = [-1.05*sp.ones(nsaves,dtype=float,order=’F’), 
      -0.95*sp.ones(nsaves,dtype=float,order=’F’)] 
z = [-10.5*sp.ones(nsaves,dtype=float,order=’F’), 
         -9.5*sp.ones(nsaves,dtype=float,order=’F’)] 
t = [sp.zeros(nsaves, dtype=float, order=’F’)] 
Ra = [27.0, 28.0, 29.0] 
Pr = [9.0, 10.0, 11.0] 
b = [8.0/3.0 - 0.1, 8.0/3.0 + 0.1] 
h = [T / float(nt)] 
tol = [1e-14] 
maxits = [20] 
forcing = [0, 1] 
bx = [0.0] 
by = [-1.0] 
bz = [-10.0]

Then a call to itertools.product() gives us the Cartesian Product of all of those lists, which is a full-factorial set of the inputs (2 ⋅ 2 ⋅ 2 ⋅ 3 ⋅ 3 ⋅ 2 ⋅ 2 = 288 treatments in this example).

# this gives us an iterator of every combination of our factor levels: 
full_factorial = itertools.product( 
        x, y, z, t, Ra, Pr, b, h, tol, maxits, forcing, bx, by, bz, savevery)

Once we have the iterator defined, we can calculate all of the trajectories associated with the input treatments in the experiment:

# iterate over all of the ’treatments’ in our design, the iterator 
# yields the arguments directly: 
for args in full_factorial: 
    L63.time_loop(*args) 
    xresult.append(args[0].copy()) 
    yresult.append(args[1].copy()) 
    zresult.append(args[2].copy()) 
    tresult.append(args[3].copy())

The resulting trajectories are shown in Figure 1. It’s pretty clear that realistic uncertainties in parameters and initial conditions lead quickly to pretty diffuse predictive distributions. Part of the challenge of attribution studies is determining the effect of each change in parameter or forcing. In realistic models, the number of ’factors’ to vary is quite large.

---

(a) X component (b) Y component (c) Z component Figure 1: Trajectories for the Full-Factorial Experiment

---

With such a low-dimension system, and few parameters it is easy enough to run the whole full-factorial, smarter experimental designs are required for experiments based on large simulation runs. A good reference on this is Design and Analysis of Computer Experiments.

Converging and Diverging Views

I was brushing up on my maximum entropy and probability theory the other day and came across a great passage in Jaynes' book about convergence and divergence of views. He applies basic Bayesian probability theory to the concept of changing public opinion in the face of new data, especially the effect prior states of knowledge (prior probabilities) can have on the dynamics. The initial portion of section 5.3 is reproduced below.

5.3 Converging and diverging views (pp. 126 – 129)

Suppose that two people, Mr A and Mr B have differing views (due to their differing prior information) about some issue, say the truth or falsity of some controversial proposition S. Now we give them both a number of new pieces of information or ’data’, D₁,D₂,…,D_n, some favorable to S, some unfavorable. As n increases, the totality of their information comes to be more nearly the same, therefore we might expect that their opinions about S will converge toward a common agreement. Indeed, some authors consider this so obvious that they see no need to demonstrate it explicitly, while Howson and Urbach (1989, p. 290) claim to have demonstrated it.

Nevertheless, let us see for ourselves whether probability theory can reproduce such phenomena. Denote the prior information by I_A, I_B, respectively, and let Mr A be initially a believer, Mr B be a doubter:

(5.16)

after receiving data D, their posterior probabilities are changed to

(5.17)

(5.17)

If D supports S, then since Mr A already considers S almost certainly true, we have P(D|S I_A), and so

(5.18)

Data D have no appreciable effect on Mr A’s opinion. But now one would think that if Mr B reasons soundly, he must recognize that P(D|S I_B) > P(D|I_B), and thus

(5.19)

Mr B’s opinion should be changed in the direction of Mr A’s. Likewise, if D had tended to refute
S, one would expect that Mr B’s opinions are little changed by it, whereas Mr A’s will move in the direction of Mr B’s. From this we might conjecture that, whatever the new information D, it should tend to bring different people into closer agreement with each other, in the sense that

(5.20)

Although this can be verified in special cases, it is not true in general.

Is there some other measure of ‘closeness of agreement’ such as log[P(S|D I_a)∕P(S|D I_B], for which this converging of opinions can be proved as a general theorem? Not even this is possible; the failure of probability theory to give this expected result tells us that convergence of views is not a general phenomenon. For robots and humans who reason according to the consistency desiderata of Chapter 1, something more subtle and sophisticated is at work.

Indeed, in practice we find that this convergence of opinions usually happens for small children; for adults it happens sometimes but not always. For example, new experimental evidence does cause scientists to come into closer agreement with each other about the explanation of a phenomenon.

Then it might be thought (and for some it is an article of faith in democracy) that open discussion of public issues would tend to bring about a general consensus on them. On the contrary, we observe repeatedly that when some controversial issue has been discussed vigorously for a few years, society becomes polarized into opposite extreme camps; it is almost impossible to find anyone who retains a moderate view. The Dreyfus affair in France which tore the nation apart for 20 years, is one of the most thoroughly documented examples of this (Bredin, 1986). Today, such issues as nuclear power, abortion, criminal justice, etc., are following the same course. New information given simultaneously to different people may cause a convergence of views; but it may equally well cause a divergence.

This divergence phenomenon is observed also in relatively well-controlled psychological experiments. Some have concluded that people reason in a basically irrational way; prejudices seem to be strengthened by new information which ought to have the opposite effect. Kahneman and Tversky (1972) draw the opposite conclusion from such psychological tests, and consider them an argument against Bayesian methods.

But now in view of the above ESP example, we wonder whether probability theory might also account for this divergence and indicate that people may be, after all, thinking in a reasonably rational, Bayesian way (i.e. in a way consistent with their prior information and prior beliefs). The key to the ESP example is that our new information was not

S ≡ fully adequate precautions against error or deception were taken, and Mrs Stewart did in fact deliver that phenomenal performance.

It was that some ESP researcher has claimed that S is true. But if our prior probability for S is lower than our prior probability that we are being deceived, hearing this claim has the opposite effect on our state of belief from what the claimant intended.

The same is true in science and politics; the new information a scientist gets is not that an experiment did in fact yield this result, with adequate protection against error. It is that some colleague has claimed that it did. The information we get from TV evening news is not that a certain event actually happened in a certain way; it is that some news reporter claimed that it did.

Scientists can reach agreement quickly because we trust our experimental colleagues to have high standards of intellectual honesty and sharp perception to detect possible sources of error. And this belief is justified because, after all, hundreds of new experiments are reported every month, but only about once in a decade is an experiment reported that turns out later to have been wrong. So our prior probability for deception is very low; like trusting children, we believe what experimentalists tell us.

In politics, we have a very different situation. Not only do we doubt a politician’s promises, few people believe that news reporters deal truthfully and objectively with economic, social, or political topics. We are convinced that virtually all news reporting is selective and distorted, designed not to report the facts, but to indoctrinate us in the reporter’s socio-political views. And this belief is justified abundantly by the internal evidence in the reporter’s own product – every choice of words and inflection of voice shifting the bias invariably in the same direction.

Not only in political speeches and news reporting, but wherever we seek for information on political matters, we run up against this same obstacle; we cannot trust anyone to tell us the truth, because we perceive that everyone who wants to talk about it is motivated either by self-interest or by ideology. In political matters, whatever the source of information, our prior probability for deception is always very high. However, it is not obvious whether this alone can prevent us from coming to agreement.

Jaynes, E.T., Probability Theory: The Logic of Science (Vol 1), Cambridge University Press, 2003.

Antarctic Ice Melt Lowest Ever Measured

That's the sensational headline anyway. Is it part of a significant downward trend though? Here's a graph of ice melt anomaly from the paper: An updated Antarctic melt record through 2009 and its linkages to high-latitude and tropical climate variability

You can use the g3data software to pull data points (that's what I did) if you want to run your own analysis with R. Does a simple statistical analysis support the claim that "there seems to be a downward trend"?

The R to generate the above graph is shown below.

melt = read.table("melt.dat")
attach(melt)

m.1 = lm(V2 ~ V1)

# Make confidence and prediction intervals
m.1.cinterval = predict(m.1, level=0.95, interval="confidence")
m.1.pinterval = predict(m.1, level=0.95, interval="prediction")

# plot the data and the fits/intervals to a png file
png("melt.png", width=640, height=480)

plot(V1, V2, ylab="Melting Anomaly", xlab="Year")
lines(V1, m.1.cinterval[,1], lty=1)
lines(V1, m.1.cinterval[,2], lty=2)
lines(V1, m.1.cinterval[,3], lty=2)
lines(V1, m.1.pinterval[,2], lty=2)
lines(V1, m.1.pinterval[,3], lty=2)
title("Antarctic Summer Melt Anomaly")
dev.off()

And here's the summary table for the linear model.

Call:
lm(formula = V2 ~ V1)

Residuals:
Min 1Q Median 3Q Max
-1.7507 -0.7252 -0.1028 0.7953 2.2894

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 60.54322 42.88214 1.412 0.169
V1 -0.03035 0.02150 -1.412 0.169

Residual standard error: 0.983 on 28 degrees of freedom
Multiple R-squared: 0.06646, Adjusted R-squared: 0.03312
F-statistic: 1.993 on 1 and 28 DF, p-value: 0.1690

What does that model summary output mean? There really is NOT a significant downward trend of melt anomaly with years (there's no first order trend in fact). It's not measurably different from noise. How does this square with the recent reports of runaway melting though?

Whether the ice is melting to much or too little,

The public is not well served by this constant drumbeat of false alarms fed by computer models manipulated by advocates.
-- DR. DAVID WOJICK, UN IPCC expert reviewer, PhD in Philosophy of Science, co-founded Department of Engineering and Public Policy at Carnegie-Mellon University

These recent posts about climate change stuff were inspired by a post I read about climate change skeptics, which I found because of my Google alerts on things related to 'computational fluid dynamics'. In the post she mentions Freeman Dyson, he's a pretty smart guy.

Concerning the climate models, I know enough of the details to be sure that they are unreliable. They are full of fudge factors that are fitted to the existing climate, so the models more or less agree with the observed data.

It is much easier for a scientist to sit in an air-conditioned building and run computer models, than to put on winter clothes and measure what is really happening outside in the swamps and the clouds. That is why the climate model experts end up believing their own models.
-- Dr. Freeman Dyson, Professor Emeritus of Physics at the Institute for Advanced Study at Princeton, fellow of the American Physical Society, member of the US National Academy of Sciences, and a fellow of the Royal Society of London

My favourite quote on the whole mess, a level-headed engineer from MIT quoted in a short article:

Mort Webster, assistant professor of engineering systems, warned that public discussion over climate change policies gets framed as a debate between the most extreme views on each side. "The world is ending tomorrow, versus it's all a myth," he said. "Neither of those is scientifically correct or socially useful."

Global Cooling?

While I don't appreciate most of the pseudo-religious global warming hysterics, I can't say the recent BBC piece on 'global cooling' is anything but sensationalism. It's curious that they don't just show their readers a simple graph of the data.

I don't know about you, but the recent years that they are claiming show a cooling trend look well within the expected variability. The solid line on the graph is a simple linear regression, the inner dashed lines are a 0.95-level confidence interval and the outer dashed lines are a 0.95-level prediction interval. The prediction interval gives a range of where you would expect temperature measurements to fall based on the fitted trend and residuals. You can download the NASA temperature data yourself and have a look.

The R code to make this graph is shown below

temps = read.csv("temps.csv", header=TRUE)
attach(temps)

# fit two models, one for January to December annual averages, and one
# for December to November annual averages
m.1 = lm(JD ~ Year)
m.2 = lm(DN ~ Year)

# Make confidence and prediction intervals
m.1.cinterval = predict(m.1, level=0.95, interval="confidence")
m.2.cinterval = predict(m.2, level=0.95, interval="confidence")
m.1.pinterval = predict(m.1, level=0.95, interval="prediction")
m.2.pinterval = predict(m.2, level=0.95, interval="prediction")

# plot the data and the fits/intervals to a png file
png("temps.png", width=640, height=480)

plot(Year, JD, ylab="Degrees C")
lines(Year, m.1.cinterval[,1], lty=1)
lines(Year, m.1.cinterval[,2], lty=2)
lines(Year, m.1.cinterval[,3], lty=2)
lines(Year, m.1.pinterval[,2], lty=2)
lines(Year, m.1.pinterval[,3], lty=2)
title("NASA GISS Global Average Temperature")

dev.off()

In case you're interested here's the summary of the fitted model shown in the graph. Just looking at the graph shows there's probably something more complicated than just a linear trend going on though.

Call:
lm(formula = JD ~ Year)

Residuals:
Min 1Q Median 3Q Max
-0.317448 -0.093632 0.001046 0.087362 0.298468

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.9770910 0.5788968 5.143 1.01e-06 ***
Year 0.0056581 0.0002977 19.009 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1244 on 126 degrees of freedom
Multiple R-squared: 0.7414, Adjusted R-squared: 0.7394
F-statistic: 361.3 on 1 and 126 DF, p-value: < 2.2e-16

To be fair, the BBC piece does go into greater detail that it isn't really global cooling, that the Pacific Decadal Oscillation might be starting to kick in, so we'll probably see a little 20 year dip like that between 1940 and 1960. I'd say you need a few more years of data to actually be able to call it though (kind of like calling an economic recession).

My Own Little Pareto Distribution

I signed up for Google Analytics, which collects bunches of neat data. One of the interesting reports is on the number of views each page on your website gets. I plotted them (in Octave of course) with pages across the x-axis, number of views on the y-axis, it looks something like a power law distribution.

x = [11,6,4,2,1,1,1,1,1,1];
bar(x);
print("-dpng","pageviews.png");

Judging by the numbers of pageviews I fear my website is in the trivial many, but that's ok a few people found the octave posts useful.

The Shuffle

The shuffle, also known as the Fisher's Exact Test, is a permutation test that can be used to estimate the sampling distribution of a statistic without relying on parametric assumptions. This is especially important when sample sizes are small. The other neat thing about permutation tests is that you don't have to know what the distribution of your statistic is. So if you have a really odd function of your data that you want to use rather than one of the classical statistics you can.

Octave has a great built-in function called nchoosek which makes shuffling a breeze. Called with scalar arguments it returns the value of the binomial coefficient, which is the number of ways you can choose k things from n things (k<=n). For fixed n, nchoosek is maximum when k=n/2. That is plotted below for n=2:24.

Notice that's a log-log scale, as n increases it quickly becomes intractable to perform a complete permutation test.

Suppose we have two data vectors, and we want to know if they are from populations with different means. We can use the normrnd() function to draw 8 samples from a normal distribution with 0 mean and 8 samples from a normal distribution with a mean of 1.

n=8;
x1 = normrnd(0,1,n,1);
x2 = normrnd(1,1,n,1);
perms = nchoosek( [x1;x2], n );

We stored all 12870 permutations of the vectors (16 choose 8) in the matrix perms. Now we use the built-in function complement to find the difference of the means under each of those labellings.

for(i=1:length(m1) )
comps(i,:) = complement( perms(i,:), [x1;x2] );
endfor
m1 = mean( perms, 2 ); % take row-wise averages, DIM=2
m2 = mean( comps, 2 ); % take row-wise averages, DIM=2
m_shuff = m1 - m2;

Now we can use the distribution of m_shuff to decide if the difference in the means of our two data vectors is significant.


Octave script with all of these calculations: shuffle.m.