A look at the Gamma distribution, possibly for rainfall modeling

Andy Long
Spring, 2018

Issues:

Rainfall data has a lot of zeros, which are troublesome for gamma distributions. On the other hand, zeros aren’t really zeros. That’s a curious notion, but actually rainfall data at zero just means “below the detection limit”. So it seems to me that we could find the detection limit dl, and then replace zeros with random uniform values in the interval (0,dl). Then do a fit, perhaps to the cumulative (cdf, which is better behaved than the pdf) to estimate  α and β, the shape and scale parameters.

Thanks to 
http://community.wolfram.com/groups/-/m/t/86994
for this function, that turns manipulate commands into animated gifs:

ManToGif[man_,name_String,step_Integer: 1]:=Module[{moviename},moviename=name<>If[$OperatingSystem==="MacOSX",".mov",".avi"];
Export[name<>".gif",Import[Export[moviename,man],"ImageList"][[1;;-1;;step]]];
DeleteFile[moviename];]

Now let’s get some data, and we’ll see if detect any significant changes over time in a particular city and its chosen month’s distribution.

CityList={"Lome","Tabligbo","Kouma-Konda","Atakpame","Sotouboua","Sokode","Kara","Niamtougou","Mango","Dapaong"};
MonthList={"January", "February", "March", "April", "May", "June", "July", "August", "September", "October", "November", "December"};

Choose your city:

city=1
cityChoice=CityList[[city]]
RainData=Import[
StringJoin[{"http://ceadserv1.nku.edu/longa/classes/mat375/data/BinaoData/Pluviometry/",cityChoice,".csv"}]
];

Choose your month:

month=8
monthChoice=MonthList[[month]]
month=Rest[Transpose[RainData][[1+month]]];
year=Rest[Transpose[RainData][[1]]];
data=Transpose[{year,month}];
ListPlot[data]

Now we’ll create a cumulative, and then use non-linear regression to fit a gamma model to that cumulative, plotting it against the data:

Clear[α,β];
cumulative=Table[i,{i,Length[year]}]/Length[year]*1.0;
janraindata=Transpose[{Sort[month],cumulative}];
meanMonth=Mean[month];
variMonth=Variance[month];
betastart=variMonth/meanMonth;
alphastart=meanMonth/betastart;
nlm=NonlinearModelFit[janraindata,CDF[GammaDistribution[α,β],x],{{α,alphastart},{β,betastart}},x];
params=nlm["BestFitParameters"]
α=α/.params
β=β/.params
gammaMean=α*β;
gammaVariance=α*β^2;
gammaMax=(α-1)β;
lp=ListPlot[janraindata];
bigplot=Plot[CDF[GammaDistribution[α,β],x],{x,0,Max[month]},Filling->Axis,PlotRange->{{0,Max[month]},{0,1}},PlotLabel->"α shape, β scale"];
plot=Show[bigplot,lp]

Export[StringJoin[{cityChoice,monthChoice,".png"}],plot,ImageSize->500,ImageResolution->1000]

In the following animation, I subset the data 20 years at a time, and then compute a model. One can move through time, and get a succession of pictures of the fits:

I start in the middle, and move backwards to 1 -- then reverse the plots, and append the final plots, to get a set that runs from year 1 to year 36 or so. I do this because I figure that the middle 20 will look the most like the average, and I start with the parameters at the global average. 

Then, as I move forward or backwards one year, I use the previous α and β, assuming that the optimal next pair of parameters shouldn’t change too much (hopefully swapping one point for another won’t change things too dramatically, but it sometimes does!).

alphakeep=α;
betakeep=β;
animationruns={};
alphastart=alphakeep;
betastart=betakeep;
For[yr=19,yr>0,yr--,
Module[
{cumulative=Table[i,{i,20}]/20.0,
data,
rain,
params,
α,
β,
lp,
plot,
nlm,
gammaMean,
gammaVariance,
gammaMax,
run
},
rain=Sort[month[[yr;;yr+19]]];
data=Transpose[{rain,cumulative}];
nlm=NonlinearModelFit[data,CDF[GammaDistribution[Abs[a],Abs[b]],x],{{a,alphastart},{b,betastart}},x];
params=nlm["BestFitParameters"];
α=a/.params;
β=b/.params;

(* Reset alphastart and betastart: *)
alphastart=α;
betastart=β;

gammaMean=α*β;
gammaVariance=α*β^2;
gammaMax=(α-1)β;
run={1960+yr,α,β,gammaMean,gammaVariance};

lp=ListPlot[data,PlotLabel->run,PlotRange->{{0,Max[cumulative]},{0,1}}];
plot=Plot[nlm[x],{x,0,Max[cumulative]},Filling->Axis,PlotRange->{{0,Max[cumulative]},{0,1}},PlotLabel->"α shape, β scale"];
animationruns=Append[animationruns,Show[lp,plot,bigplot]];
]
]
animationruns=Reverse[animationruns];
alphastart=alphakeep;
betastart=betakeep;
For[yr=20,yr<=Length[year]-19,yr++,
Module[
{cumulative=Table[i,{i,20}]/20.0,
data,
rain,
params,
α,
β,
lp,
plot,
nlm,
gammaMean,
gammaVariance,
gammaMax,
run
},
rain=Sort[month[[yr;;yr+19]]];
data=Transpose[{rain,cumulative}];
nlm=NonlinearModelFit[data,CDF[GammaDistribution[Abs[a],Abs[b]],x],{{a,alphastart},{b,betastart}},x];
params=nlm["BestFitParameters"];
α=a/.params;
β=b/.params;

(* Reset alphastart and betastart: *)
alphastart=α;
betastart=β;

gammaMean=α*β;
gammaVariance=α*β^2;
gammaMax=(α-1)β;
run={1960+yr,α,β,gammaMean,gammaVariance};

lp=ListPlot[data,PlotLabel->run,PlotRange->{{0,Max[cumulative]},{0,1}}];
plot=Plot[nlm[x],{x,0,Max[cumulative]},Filling->Axis,PlotRange->{{0,Max[cumulative]},{0,1}},PlotLabel->"α shape, β scale"];
animationruns=Append[animationruns,Show[lp,plot,bigplot]];
]
]

man=Manipulate[
Show[animationruns[[yr]]]
,{{yr,1},1,Length[animationruns],1}
,Deployed->True,FrameMargins->0
];

ManToGif[man,StringJoin[{cityChoice,monthChoice}],1]