What if we use seekPeriods on known data, with noise?

Andy Long
Spring, 2018

We will see that seekPeriods can detect some of our periods, but doesn’t work as well on some others.

Simulate data, with yearly oscillation:

Numyears=55
years=Table[t,{t,1/24,Numyears,1/12}];
data=Table[.012 t+5*Sin[2Pi t]+15*Cos[2 Pi t],{t,1/24,Numyears,1/12}];
data=data+RandomVariate[NormalDistribution[0,3],Length[data]];
datamat=Transpose[{years,data}];
ListPlot[datamat]

A function we need:

seekPeriods[p_,data_]:=
Module[
{lm=LinearModelFit[data,{x,Sin[2π x/p],Cos[2π x/p]},x]},
lm["AdjustedRSquared"]
]

seekPeriodsRSq[p_,data_]:=
Module[
{lm=LinearModelFit[data,{x,Sin[2π x/p],Cos[2π x/p]},x]},
lm["RSquared"]
]

Let’s look for our periods!

Plot[seekPeriods[p,datamat],{p,1/6,1.1},PlotRange->All]
Plot[seekPeriods[p,datamat],{p,1.1,25},PlotRange->All]
Plot[seekPeriods[p,datamat],{p,25.1,50},PlotRange->All]

lreg=LinearModelFit[Transpose[{years,data}],{x,Sin[2Pi x],Cos[2 Pi x]},x]
resids=Transpose[{years,lreg["FitResiduals"]}];

Plot[seekPeriods[p,resids],{p,1/6,1.1},PlotRange->All]
Plot[seekPeriods[p,resids],{p,1.1,25},PlotRange->All]
Plot[seekPeriods[p,resids],{p,25.1,50},PlotRange->All]

Now let’s simulate data with a period 12 oscillation, as well:

data=Table[.012 t+5*Sin[2Pi t]+15*Cos[2 Pi t]
+1*Sin[2Pi t/12]-.4*Cos[2 Pi t/12]
,{t,1/24,Numyears,1/12}];
data=data+RandomVariate[NormalDistribution[0,3],Length[data]];
datamat=Transpose[{years,data}];
ListPlot[Transpose[{years,data}]]

Let’s look for our periods!

Plot[seekPeriods[p,datamat],{p,1/6,1.1},PlotRange->All]
Plot[seekPeriods[p,datamat],{p,1.1,25},PlotRange->All]
Plot[seekPeriods[p,datamat],{p,25.1,50},PlotRange->All]

seekPeriodsRSq[8,datamat]
seekPeriods[8,datamat]

Let’s do an incomplete removal, taking the dominant period -- then what happens to the period 12 piece?

lreg=LinearModelFit[Transpose[{years,data}],{x,Sin[2Pi x],Cos[2 Pi x]},x]
resids=Transpose[{years,lreg["FitResiduals"]}];

Plot[seekPeriods[p,resids],{p,1/6,1.1},PlotRange->All]
Plot[seekPeriods[p,resids],{p,1.1,25},PlotRange->All]
Plot[seekPeriods[p,resids],{p,25.1,50},PlotRange->All]

lreg=LinearModelFit[Transpose[{years,data}],{x,Sin[2Pi x],Cos[2 Pi x],Sin[2Pi x/12],Cos[2 Pi x/12]},x]
resids=Transpose[{years,lreg["FitResiduals"]}];

Plot[seekPeriods[p,resids],{p,1/6,1.1},PlotRange->All]
Plot[seekPeriods[p,resids],{p,1.1,25},PlotRange->All]
Plot[seekPeriods[p,resids],{p,25.1,50},PlotRange->All]

Now let’s simulate data with a period 11 oscillation, as well:

data=Table[.012 t+5*Sin[2Pi t]+15*Cos[2 Pi t]
+1*Sin[2Pi t/11]-.4*Cos[2 Pi t/11]
,{t,1/24,Numyears,1/12}];
data=data+RandomVariate[NormalDistribution[0,3],Length[data]];
datamat=Transpose[{years,data}];
ListPlot[Transpose[{years,data}]]

Let’s look for our periods!

Plot[seekPeriods[p,datamat],{p,1/6,1.1},PlotRange->All]
Plot[seekPeriods[p,datamat],{p,1.1,25},PlotRange->All]
Plot[seekPeriods[p,datamat],{p,25.1,50},PlotRange->All]

lreg=LinearModelFit[Transpose[{years,data}],{x,Sin[2Pi x],Cos[2 Pi x],Sin[2Pi x/11],Cos[2 Pi x/11]},x]
resids=Transpose[{years,lreg["FitResiduals"]}];

Plot[seekPeriods[p,resids],{p,1/6,1.1},PlotRange->All]
Plot[seekPeriods[p,resids],{p,1.1,25},PlotRange->All]
Plot[seekPeriods[p,resids],{p,25.1,50},PlotRange->All]

Now let’s simulate data with a period 15 oscillation, as well:

data=Table[.012 t+5*Sin[2Pi t]+15*Cos[2 Pi t]
+1*Sin[2Pi t/15]-.4*Cos[2 Pi t/15]
,{t,1/24,Numyears,1/12}];
data=data+RandomVariate[NormalDistribution[0,3],Length[data]];
datamat=Transpose[{years,data}];
ListPlot[Transpose[{years,data}]]

Let’s look for our periods!

Plot[seekPeriods[p,datamat],{p,1/6,1.1},PlotRange->All]
Plot[seekPeriods[p,datamat],{p,1.1,25},PlotRange->All]
Plot[seekPeriods[p,datamat],{p,25.1,50},PlotRange->All]
lreg=LinearModelFit[Transpose[{years,data}],{x,Sin[2Pi x],Cos[2 Pi x],Sin[2Pi x/15],Cos[2 Pi x/15]},x]
resids=Transpose[{years,lreg["FitResiduals"]}];
Plot[seekPeriods[p,resids],{p,1/6,1.1},PlotRange->All]
Plot[seekPeriods[p,resids],{p,1.1,25},PlotRange->All]
Plot[seekPeriods[p,resids],{p,25.1,50},PlotRange->All]

Now let’s simulate data with JUST a period 15 oscillation, as well:

data=Table[.012 t
+1*Sin[2Pi t/15]-.4*Cos[2 Pi t/15]
,{t,1/24,Numyears,1/12}];
data=data+RandomVariate[NormalDistribution[0,3],Length[data]];
datamat=Transpose[{years,data}];
ListPlot[Transpose[{years,data}]]

Let’s look for our periods!

Plot[seekPeriods[p,datamat],{p,1/6,1.1},PlotRange->All]
Plot[seekPeriods[p,datamat],{p,1.1,25},PlotRange->All]
Plot[seekPeriods[p,datamat],{p,25.1,50},PlotRange->All]
lreg=LinearModelFit[Transpose[{years,data}],{x,Sin[2Pi x/15],Cos[2 Pi x/15]},x]
resids=Transpose[{years,lreg["FitResiduals"]}];
Plot[seekPeriods[p,resids],{p,1/6,1.1},PlotRange->All]
Plot[seekPeriods[p,resids],{p,1.1,25},PlotRange->All]
Plot[seekPeriods[p,resids],{p,25.1,50},PlotRange->All]