The tale of intel reports and tail events

I remember posting something about bomb blasts on one occasion expressing utter dismay at how people have become so insensitive to such incidents. This is my second post on the same topic, this time with a completely disinterested and technical disposition.

Time and again we come across a news article claiming that the Center had warned the local or the state police about a possible terrorist attack. The recent Hyderabad blasts are no exception (see [1] for instance). The local government can take appropriate action based on such information provided by central intelligence teams. This time around the Chief minister decided not to take the information too seriously (``Andhra Pradesh Chief Minister Kiran Kumar Reddy had said those were general alerts which often keep coming from the Centre'' [1]).

I am not sure of the format in which such information is transmitted to the local ministers but let me hazard a few guesses: perhaps it is a lengthy report in teletype font typewritten on legal size eggshell finish papers, or maybe a frantic last minute phone call from PM Manmohan Singh saying something that is roughly equivalent to ``Take extra care dude.'' In any case, the minister (Mr. NKK Reddy, in the present case) probably makes up a numerical measure of the seriousness of the intel report.

Let us suppose that a report comes in every day and Mr. Reddy decides to assign a random positive number as the ``alert level of the day.'' This procedure may be as simple as just counting the number of pages in the daily intel report. Now some notation: let $X_n$ denote the alert level for the day $n$. Assume $X_n$ is exponentially distributed with a (normalized) mean of $\lambda=1$. Mr. Reddy likes being conservative---he sees no point in putting the state on high alert just because one particular $X_n$ came out really big. After all, when one observes long sequences of random values, it is quite natural to expect some really huge readings, purely out of luck.

Let $M_n = \max \{X_1, X_2, \cdots, X_n\}$ denote the largest of this sequence of i.i.d. exponentials. It is easy to derive the distribution $F_n(\cdot)$ of $M_n$ in terms of the common distribution $G(\cdot)$ of the i.i.d. exponential random variables:
\begin{array}{rcl}
F_n(t) &:=& Pr(M_n \leq t) \\
           &=& Pr(\max\{X_1,\cdots, X_n\} \leq t) \\
           &=& Pr(\bigcap_{i=1}^n \{X_i \leq t\}) \\
           &=& \prod_{i=1}^n Pr(X_i \leq t) \\
           &=& (G(t))^n \\
           &=& (1-\exp(-t))^n,
\end{array}
for $t\geq 0$.

In order to be conservative and stay on budget, Mr. Reddy can raise an alarm only if the alert number he sees is in the $\alpha \%$ tail of the distribution of the maximum. Suppose $n=822$, which is how many days Mr. Reddy has been the CM [2]. $(100-\alpha) \%$ area under the density of $M_{822}$ can be calculated: $t_\alpha = F_{822}^{-1}(1-\alpha/100)$. If $\alpha = 5\%$, for example, then Mr. Reddy should pronounce high alert when he sees a number at least as big as 9.68197. See Figure below and also see [3].

Figure: CDF of $M_n$ with location of 5% tail (click image for larger view)

I agree, this is an oversimplified (and probably useless) approach. But embellishments are possible. Foremost, the i.i.d. assumption may not be true. Some contextual information may be incorporated to create a mild correlation structure between the $X_n$'s. The other weakness in this model is that rare events like bomb blast may occur even when the alert level is quite low. (Remember what the great logodaedalist Mr. Chidambaram said in July 2011---the absence of intelligence is not failure of intelligence [4].)

References

[1] ``Specific alert was sent to Hyderabad Thursday morning: Centre'' in The Hindu, as reported by PTI available here.
[2] Wikipedia page on a list of Andhra Pradesh chief ministers is here.
[3] A calculation using Wolfram Alpha here.
[4] ``No intel about blasts: Chidambaram,'' July 14, 2011. Available here. 
(Note: All links last accessed on Feb 23, 2013.)

Lake Mendota Freeze Durations

Q. Is the freeze duration of Lake Mendota decreasing over the years?

Freeze duration data for Lake Mendota and other Madison lakes is available online. $N=156$ data points for Lake Mendota are available starting in the year 1855 [1]. Let us run a simple statistical test to check for a negative slope trend.

Consider a linear model, i.e.,
$$ y_n = \beta_1 n + \beta_0 + w_n, $$
where $y_n$ denotes the freeze duration at any data point $0 \leq n \leq N-1$ and $w_n$ is i.i.d. Gaussian noise with zero mean and unknown variance. The regression coefficients $\hat{\beta_0}$ and $\hat{\beta_1}$ can be calculated from data using linear regression (see Figure below).

Now, setting up the hypotheses:
$$
H_0: \beta_1 = 0
$$
and
$$
H_1: \beta_1 < 0.
$$

Figure: Linear fit to Lake Mendota freeze duration data

It is known that $$t = \frac{\hat{\beta_1}}{SE}$$ follows a t-distribution with $N-2$ degrees of freedom. Here $$SE = \sqrt{\frac{\frac{1}{N-2}\sum_{n=0}^{N-1} (y_n - \hat{\beta_1} n- \hat{\beta_0} )^2}{\;\;\;\;\;\;\;\sum_{n=0}^{N-1} \left(n-\frac{(N-1)(N-2)}{\;\;\;\;\;2N}\right)^2}}$$ is called the standard error [2].

A. Unless there is an algebra error in my calculations, the t-test shows that there is statistically significant negative slope trend in this data (at both 1% and 5% significance levels).

References:
[1] UW-Madison Department of Atmospheric and Oceanic Sciences. URL: http://www.aos.wisc.edu/~sco/lakes/Mendota-ice.html (Accessed: Jan 20, 2012)
[2] T. Tarpey, Course Notes for Stat-630 at Wright State University. URL: http://www.wright.edu/~thaddeus.tarpey/stt630chap8.pdf (Accessed: Jan 20, 2012).


Appendix:
I have copied data from [1] in a .csv file for easy access. You can download it here. And here's some Octave code. Make sure the .csv file is in the same directory.

data = csvread('MendotaFreezeDurations.csv');
N = numel(data(:,1));
ply = polyfit(data(:,1), data(:,2), 1);
% ply = [-0.18829   117.82053]
SSE = sum((data(:,2)-polyval(ply, data(:,1))).^2);
% SSE =  4.6469e+04
MSE = SSE/(N-2);
SE = sqrt(MSE/sum( (data(:,1)-mean(data(:,1))).^2 ));
% SE =  0.030491
tstatistic = ply(1)/SE;
% tstatistic = -6.1754
tthresh1 = -tinv(0.99,N-2)
tthresh5 = -tinv(0.95,N-2)
if tstatistic < tthresh1, 
  disp('H0 rejected at 1%'); 
else, 
  disp('H0 accepted at 1%'); 
end;
if tstatistic < tthresh5, 
  disp('H0 rejected at 5%'); 
else, 
  disp('H0 accepted at 5%'); 
end;
% H0 is rejected in both cases. 
% There is a statistically significant 
% downward sloping trend.
 
plot(data(:,1), data(:,2),... 
  data(:,1), polyval(ply, data(:,1)));
axis tight
ylabel('freeze duration (days)')
xlabel('data point')
title('lake Mendota freeze duration linear regression')
print -depsc freezedurplot.eps