Calculate proportion of values for dates

Question

At work, we were discussing when is it the best time to change to winter tires for bikes and/or cars.

Using WeatherData[] and DateListPlot[], it was fairly straightforward for me to create the diagram below:

Fig 1 Mean temperature per day in Stockholm, red is negative and a risk without winter tires.

The code for this is

cityTemp = WeatherData["Stockholm", "MeanTemperature", {{1977, 1, 1}, {2011, 12, 31}, "Day"}];
iceRiskDays = Select[cityTemp, Last[#] < 0 &];

yearStrip[ dataItem_] := {{0, Part[First[dataItem], 2], Last[First[dataItem]]}, Last[dataItem]}

DateListPlot[{yearStrip[#] & /@ cityTemp, yearStrip[#] & /@ iceRiskDays} ]

My question is: How do I calculate for each day the proportion of values for that day that are below 0 Celsius? (e.g. for a date at the end of November, the proportion is likely to be bigger than 0.5)?

My attempts to do this ended with trying to create separate lists for each date, but I felt that this was a less elegant way and also creates less "fit" with DateListPlot.

Answer 1

Nice question. First I wrote your code in this way. Where we don't need the function yearStrip

{date,year,month,day,temp}={1,1,2,3,2}
cityTemp=WeatherData["Stockholm","MeanTemperature",{{1977,1,1},{2011,12,31},"Day"}];
cityTemp[[All,date,year]]=0
iceRiskDays=Select[cityTemp,#[[temp]]<0&];
DateListPlot[{cityTemp,iceRiskDays}]

Where we get the same result:

Now we can group the day-month information using GatherBy in this way:

dataGather=GatherBy[cityTemp,#[[date,{month,day}]]&];
getNegativeProportion[data_]:=N@Total[UnitStep[-data]]/Length[data]
negativeProportionList=Sort[{#[[1,1]],getNegativeProportion[#[[All,temp]]]}&/@dataGather];
movAvg=MovingAverage[{N@AbsoluteTime@#[[date]],#[[temp]]}&/@negativeProportionList,14];
DateListPlot[{negativeProportionList,movAvg},Joined-> True,PlotStyle-> {Blue,Red}]

Where we get:

With the red line as the moving average for 14 days!

Answer 2

Following @J.M. if you are after the likelihood of icy days as a function of month

ff[x_] := Select[cityTemp, #[[1, 2]] == x &] // (Count[#, {d_List, t_?Negative}]/Length[#]) &

{Range[12] - 1/2, ff /@ Range[12]} // Transpose // 
ListLinePlot[#, InterpolationOrder -> 0,AxesLabel -> {"Month", "likelihood"}] &

which suggests in January its about 15% (from where I live) and 60 % in Stockholm (Brrr...)

Comparing Paris to Stockholm (using the nice package autoLegend)

ff[x_, temp_] := Select[temp, #[[1, 2]] == x &] // (Count[#, {d_List, t_?Negative}]/
 Length[#]) &

{{Range[12] - 1/2, ff[#, cityTempS] & /@ Range[12]} // Transpose,
{Range[12] - 1/2, ff[#, cityTempP] & /@ Range[12]} // Transpose} // 
ListLinePlot[#, InterpolationOrder -> 0,
AxesLabel -> {"Month", "likelihood"}] &//autoLegend[#, {"Stockholm","Paris"}] &

Stealing yet again another suggestion of @J.M. we can look at global warming ;-) This time we look at the faction of days with ice as a function of year:

 dat = Map[Count[#, {d_List, t_?Negative}]/Length[#] &, 
  SplitBy[cityTempP, {#[[1, 1]] &, #[[1, 2]] &}], {2}] // Rest;
 fP = LinearModelFit[Mean /@ dat, {1, x}, x];

 dat = Map[Count[#, {d_List, t_?Negative}]/Length[#] &, 
 SplitBy[cityTempS, {#[[1, 1]] &, #[[1, 2]] &}], {2}] // Rest;
 fS = LinearModelFit[Mean /@ dat, {1, x}, x];

 Show[Plot[{fS[x], fP[x]}, {x, 1, 31}, AxesLabel -> {"year", "Likelihood"}],
 ListPlot[{fS["Data"], fP["Data"]}]]// autoLegend[#, {"Stockholm", "Paris"}] &

And finally for the sake of contributing something original (the code for FilledErrorPlot), let us launch the climate-skeptic versus climato-believer debate while illustrating the famous Lies, damned lies, and statistics methodology on two different ways of plotting the same data; guess which is which ;-)

Show[dat // Transpose // FilledErrorPlot[#] &,
Plot[{fS[x]}, {x, 1, 33}, 
PlotStyle -> Directive[AbsoluteThickness[2], ColorData[10][1]], 
AxesLabel -> {"year", "Likelihood"}], AspectRatio -> 1.5, 
PlotRange -> {0.035, 0.45}, AxesOrigin -> {0., 0.035}, 
AxesLabel -> {"year", "Fraction of \n icy days"}]

versus

Show[dat // Transpose // FilledErrorPlot[#, ErrorOnMean -> False] &,
Plot[{fS[x]}, {x, 1, 33}, 
PlotStyle -> Directive[AbsoluteThickness[2], ColorData[10][1]], 
AxesLabel -> {"year", "Likelihood"}], AspectRatio -> 0.5, PlotRange -> {-0.5, 1}, 
AxesLabel -> {"year", "Fraction of \n icy days"}]

Answer 3

Here is a slightly different approach using GatherBy.

Reformat dates to be just {month, day}:

shortDates = {DateString[First@#, { "Month","Day",}], Last@#} & /@  cityTemp;

Group the values by day and month, count the days on which the temperature is not positive:

probs = Count[#[[All, 2]], t_?NonPositive]/Length@# & /@    
        GatherBy[Sort@shortDates, #[[1]] &];

Plot the likelihoods:

   ListLinePlot[probs, AxesLabel -> { "Day Number","Likelihood"}]

You can calculate different likelihoods by adjusting the short date conversion string, for a month value just switch the date reformatting to months:

shortDates = {DateString[First@#, {"Month"}], Last@#} & /@ cityTemp;

Answer 4

If my understanding of the question is correct, this does the job:

cityTemp = WeatherData["Stockholm", "MeanTemperature",
                       {{1977, 7, 1}, {2011, 12, 31}, "Day"}];

(* proportions; I ignore February 29 for this *)
props = Array[(Count[#, {_List, _?Negative}]/Length[#]) &[Cases[cityTemp,
              {Prepend[Rest[DatePlus[{2011, 1, 1}, #]], yr_], temp_}, Infinity]] &, 365, 0];

Take[N[props], 5] (* sample *)
   {0.705882, 0.588235, 0.705882, 0.676471, 0.705882}

To plot these proportions, you can do this:

DateListPlot[props, {{2011, 1, 1}, {2011, 12, 31}, "Day"},
             FrameLabel -> {None, "Proportion"}, Joined -> True]

Pick the day of the year with the most number of "cold" instances:

Pick[Array[Rest[DatePlus[{2011, 1, 1}, #]] &, 365, 0], props, Max[props]]
   {{2, 18}}

So, at least for Stockholm, February 18 is often a cold day.