Checking out the BG Climate Normals

Andy Long
Spring, 2020

Let’s load in the data, and see what we’ve got. These were obtained from NOAA, and I stashed a copy so that I would have some control over it:

BGNormals=Import["http://ceadserv1.nku.edu/longa//classes/mat375/days/docs/Fletcher/data/BowlingGreenDailyClimateNormals.csv","CSV"];
BGheader=BGNormals[[1]]
BGNormalData=BGNormals[[2;;]];

Now I want to pick off the important variables: max, min, and dtr.

(* Dumb counting method...:) *)
navg=13*4+6+1;
ndtr=navg+2;
nmax=navg+4;
nmin=navg+6;
navgstddev=navg+8;
ndtrstddev=navg+10;
nmaxstddev=navg+12;
nminstddev=navg+14;
indices={navg,ndtr,nmax,nmin,navgstddev,ndtrstddev,nmaxstddev,nminstddev};
vnames=BGheader[[indices]]

BGNormalTemps=BGNormalData[[All,indices]];
BGNormalTemps[[1,All]]
{avg,dtr,max,min,avgstddev,dtrstddev,maxstddev,minstddev}=Transpose[BGNormalTemps];

Averages and DTR, presented as data in the data set, are actually just derived from max and min normals, as these commands show:

days=Table[i,{i,1,365}];
p1=ListPlot[Transpose[{days,min}],PlotLabel->"BG Normals"];
p2=ListPlot[Transpose[{days,avg}],PlotLabel->"BG Normal Mean"];
p3=ListPlot[Transpose[{days,max}],PlotLabel->"BG Normal Mean"];
p4=ListPlot[Transpose[{days,dtr}],PlotLabel->"BG Normals"];
p5=ListPlot[Transpose[{days,max-min}],PlotLabel->"BG Normals",PlotStyle->Green];
p6=ListPlot[Transpose[{days,(max+min)/2}],PlotLabel->"BG Normals",PlotStyle->Orange];
Show[p1,p2,p3,p4,p5,p6,PlotRange->All]

ListPlot[Transpose[{days,dtr-(max-min)}],PlotLabel->"BG Normals DTR calcs"]
ListPlot[Transpose[{days,avg-(max+min)/2}],PlotLabel->"BG Normals AVG calcs"]

Now let’s do some regression:

Regressions for Minimum Normals

We can do the regressions either as linear regressions, or as non-linear regressions -- and we’ll get exactly the same results! It’s because of the magic of the trig identity, which allows us to break the Sin[2π(t-t0)] term into Sine and Cosine terms. 

If we do them as linear regressions, the R-squared terms are more useful.

This first regression just fits 
a + b Cos[2 Pi x] + c Sin[2Pi x]

pointOfYear=days/365;
data=Transpose[{pointOfYear,min}];
minplot=ListPlot[data,PlotLabel->"Minimum -- BG Normals"];
minlm=LinearModelFit[Cases[data,{_?NumericQ..}],{Cos[2 Pi x],Sin[2Pi x]},x]
minfits=minlm["FitResiduals"];
minlm["ParameterTable"]
minlm["RSquared"]
minlm["AdjustedRSquared"]
minlm["ParameterConfidenceIntervals"]
Show[minplot,Plot[minlm[x],{x,0,1}]]
Histogram[minfits]

This second regression adds a linear term (d x): 
a + d x +  b Cos[2 Pi x] + c Sin[2Pi x]
and we discover that the d is significantly different from 0, and positive -- 1.9, which, on an annual basis would represent quite an amazing increase in minimum temperature!

minlm=LinearModelFit[Cases[data,{_?NumericQ..}],{x, Cos[2 Pi x],Sin[2Pi x]},x]
minfits=minlm["FitResiduals"];
minlm["ParameterTable"]
minlm["RSquared"]
minlm["AdjustedRSquared"]
minlm["ParameterConfidenceIntervals"]
Show[minplot,Plot[minlm[x],{x,0,1}]]
Histogram[minfits]

Now I’ll do the regression using non-linear regression. In a way, it’s cleaner: I don’t need to deduce the amplitude and phase shift -- I’m staring at them! The results are identical, however -- and that’s one of the cool things. 

Again, the first model is -- well, there it is in the command for NonlinearModelFit! I like the way we write the regression model with NonlinearModelFit.

minlm=NonlinearModelFit[data,{a + amp Sin[2Pi (x - x0)]},{a, amp,{x0,.3}},x]
minfits=minlm["FitResiduals"];
minlm["ParameterTable"]
minlm["RSquared"]
minlm["AdjustedRSquared"]
minlm["ParameterConfidenceIntervals"]
Show[minplot,Plot[minlm[x],{x,0,1}]]
Histogram[minfits]

Our model has an amplitude of 22.24 (a total variation of twice that), and the curve hits its mean at .3063 of the year; that means it bottoms out at .3063-.25 of the year, or .0563*365.25 -- about the 21st of January. 

As before, we now add the linear term, only doing so with NonlinearModelFit.

minlm=NonlinearModelFit[data,{a + b x + amp Sin[2Pi (x - x0)]},{a,b, amp,{x0,.3}},x]
minfits=minlm["FitResiduals"];
minlm["ParameterTable"]
minlm["RSquared"]
minlm["AdjustedRSquared"]
minlm["ParameterConfidenceIntervals"]
Show[minplot,Plot[minlm[x],{x,0,1}]]
Histogram[minfits]

Regressions for Maximum Normals

I’m now going to reproduce the same type of results for the maxima, rather than the minima.

data=Transpose[{pointOfYear,max}];
maxplot=ListPlot[data,PlotLabel->"Maximum -- BG Normals"];
maxlm=LinearModelFit[data,{Cos[2 Pi x],Sin[2Pi x]},x]
maxfits=maxlm["FitResiduals"];
maxlm["ParameterTable"]
maxlm["RSquared"]
maxlm["AdjustedRSquared"]
maxlm["ParameterConfidenceIntervals"]
Show[maxplot,Plot[maxlm[x],{x,0,1}]]
Histogram[maxfits]

maxlm=LinearModelFit[data,{x,Cos[2 Pi x],Sin[2Pi x]},x]
maxfits=maxlm["FitResiduals"];
maxlm["ParameterTable"]
maxlm["RSquared"]
maxlm["AdjustedRSquared"]
maxlm["ParameterConfidenceIntervals"]
Show[maxplot,Plot[maxlm[x],{x,0,1}]]
Histogram[maxfits]

So we notice that the linear term is 1.28 or so; it’s significantly different from zero. Maxima are increasing over the course of the year.  

But not nearly as much as the minima, for which it was almost 2: minima were increasing twice as fast as maxima, according to these models.

maxlm=NonlinearModelFit[data,{a + amp Sin[2Pi (x - x0)]},{a,amp,{x0,.3}},x]
maxfits=maxlm["FitResiduals"];
maxlm["ParameterTable"]
maxlm["RSquared"]
maxlm["AdjustedRSquared"]
maxlm["ParameterConfidenceIntervals"]
Show[maxplot,Plot[maxlm[x],{x,0,1}]]
Histogram[maxfits]

Our model has an amplitude of 26.40 (a total variation of twice that), and the curve hits its mean at .3030 of the year; that means it bottoms out at .3030-.25 of the year, or .0530*365.25 -- about the 19th of January. 

As before, we now add the linear term, only doing so with NonlinearModelFit.

maxlm=NonlinearModelFit[data,{a + b x + amp Sin[2Pi (x - x0)]},{a,b, amp,{x0,.3}},x]
maxfits=maxlm["FitResiduals"];
maxlm["ParameterTable"]
maxlm["RSquared"]
maxlm["AdjustedRSquared"]
maxlm["ParameterConfidenceIntervals"]
Show[maxplot,Plot[maxlm[x],{x,0,1}]]
Histogram[maxfits]