The wolves and moose of Isle Royale, Michigan, have been well studied:
http://www.isleroyalewolf.org/http%3A//www.cpc.ncep.noaa.gov/products/precip/CWlink/pna/nao.shtml

We begin with the InsightMaker model which is available at
https://insightmaker.com/insight/2068/Isle-Royale-Predator-Prey-Interactions, and then tune the parameters so as to get a best fit to the historical data. 

Andy Long
Spring, 2018

Thanks to:
https://mathematica.stackexchange.com/questions/74214/the-lotka-volterra-predator-prey-model

Clear everything possibly in our way:

Clear["Global`*"]

Now let’s get the actual wolf/moose data, extracted from a much more extensive data set, available from 
http://www.isleroyalewolf.org/data/data/home.html

IsleRoyale=Import["http://ceadserv1.nku.edu/longa//classes/mat375/data/isleroyale_graph_data_28Dec2011.csv"];

A look at the actual data:

year=Rest[IsleRoyale[[All,1]]];
minYear=Min[year]
year=year-minYear;
Max[year]
wolfdata=Rest[IsleRoyale[[All,2]]];
moosedata=Rest[IsleRoyale[[All,3]]];
dataratio=wolfdata/moosedata;
ListPlot[Transpose[{year+minYear,moosedata}]];
ListPlot[Transpose[{year+minYear,wolfdata}]];
(* Do you see a ring-around-the-rosey? *)
ListPlot[Transpose[{moosedata,wolfdata}]]

Now let’s make some moose and wolves! 
Model and starting parameters from the InsightMaker model:

WolfDeathRate=0.25;
WolfBirthRateFactor=0.000203;
MooseBirthRate=0.25;
MooseDeathRateFactor=0.0081;
mo=421.5;
wo=20;

Here’s the solution for the ratio of wolves to moose, with starting parameters from the InsightMaker model:

lv=
NDSolve[
{(*  moose: *)
m'[t]==MooseBirthRate m[t]-MooseDeathRateFactor m[t]w[t] ,
(* wolves: *)
w'[t]==WolfBirthRateFactor m[t] w[t]-WolfDeathRate w[t],
m[0]==mo,
w[0]==wo
},
{m[t],w[t]},
{t,0,52}];
ratio[t_]=Evaluate[{w[t]/m[t]}/.lv];
p1=ListPlot[Transpose[{year,dataratio}]];
p2=Plot[ratio[t],{t,0,52}];
Show[p2,p1]

I got these parameters from optimization:

0.0086539

params={0.00020621277957375955`,0.23996157218610736`,0.26993905426786`,0.010532995023990771`,660.0000000006974`,20.00000026558482`};
{WolfBirthRateFactor,WolfDeathRate,MooseBirthRate,MooseDeathRateFactor,mo,wo}=params;
plotRatio[WolfBirthRateFactor_,WolfDeathRate_,MooseBirthRate_,MooseDeathRateFactor_,mo_,wo_]:=
Module[{lv,m,w,ratio,t},
lv=
NDSolve[
{(*  moose: *)
m'[t]==MooseBirthRate m[t]-MooseDeathRateFactor m[t]w[t] ,
(* wolves: *)
w'[t]==WolfBirthRateFactor m[t] w[t]-WolfDeathRate w[t],
m[0]==mo,
w[0]==wo
},
{m[t],w[t]},
{t,0,52}];
Plot[Evaluate[{w[t]/m[t]}/.lv],{t,0,52}]]
p3=Apply[plotRatio,params];
Show[p2,p3,p1]
Show[p3,p1]

0.00821259

params={0.000222876,0.239896,0.269851,0.0110315,660.,20.};
{WolfBirthRateFactor,WolfDeathRate,MooseBirthRate,MooseDeathRateFactor,mo,wo}=params
p4=Apply[plotRatio,params];
Show[p2,p3,p4,p1]
Show[p4,p1]

0.00806343

params={0.00023195,0.239829,0.269769,0.0113139,660.,20.};
{WolfBirthRateFactor,WolfDeathRate,MooseBirthRate,MooseDeathRateFactor,mo,wo}=params
p5=Apply[plotRatio,params];
Show[p2,p3,p4,p5,p1]
Show[p1,p5]

0.00775043

params={0.000240036,0.239826,0.269764,0.0113423,611.,20.};
{WolfBirthRateFactor,WolfDeathRate,MooseBirthRate,MooseDeathRateFactor,mo,wo}=params;
p6=Apply[plotRatio,params];
Show[p2,p3,p4,p5,p6,p1]
Show[p1,p6]

0.00690574

params={0.000267702,0.239817,0.269739,0.0113796,591.,23.};
{WolfBirthRateFactor,WolfDeathRate,MooseBirthRate,MooseDeathRateFactor,mo,wo}=params;
p7=Apply[plotRatio,params];
Show[p2,p3,p4,p5,p6,p7,p1]
Show[p1,p7]

0.00685691 (from non-linear fit):
params=
{0.0004566167858993027`,
0.4720627720562842`,
0.13518440641048893`,
0.004767210896693055`,
805.2511688586087`,
29.917958639847548`}

({{"", "Estimate", "Standard Error", "t-Statistic", "P-Value"}, {x, 1.9712433301742296`*^-8, 0.00001907406921239978`, 0.001033467640398805`, 0.9991793668829085`}})
({{"", "Estimate", "Standard Error", "t-Statistic", "P-Value"}, {x, 0.00004189231838418318`, 0.014307391809529545`, 0.0029280192324278483`, 0.9976749862348391`}})
({{"", "Estimate", "Standard Error", "t-Statistic", "P-Value"}, {x, -0.000013195043173165362`, 0.003031159304945539`, -0.004353134179268231`, 0.9965433700242906`}})
({{"", "Estimate", "Standard Error", "t-Statistic", "P-Value"}, {x, -7.202592525457392`*^-7, 0.00019896877989490763`, -0.003619961146297271`, 0.9971255472222801`}})
({{"", "Estimate", "Standard Error", "t-Statistic", "P-Value"}, {x, 0.055704739317598166`, 30.885654360631246`, 0.0018035797029640676`, 0.9985678538649706`}})
({{"", "Estimate", "Standard Error", "t-Statistic", "P-Value"}, {x, 0.002045818037591286`, 1.7944273126907373`, 0.0011400952399256501`, 0.9990946984330663`}})

params={0.0004566167858993027`,0.4720627720562842`,0.13518440641048893`,0.004767210896693055`,805.2511688586087`,29.917958639847548`}
{WolfBirthRateFactor,WolfDeathRate,MooseBirthRate,MooseDeathRateFactor,mo,wo}=params;
p8=Apply[plotRatio,params];
Show[p2,p3,p4,p5,p6,p7,p8,p1]
Show[p1,p8]

Actually I’m starting with some parameter values from an “eyeball optimization” of the wolf/moose ratio.

WolfDeathRate=0.24;
WolfBirthRateFactor=0.0002;
MooseBirthRate=0.27;
MooseDeathRateFactor=0.01;
mo=660;
wo=20;

Thanks to 
https://mathematica.stackexchange.com/questions/30781/can-graph-interpolating-function-but-not-use-it-any-thoughts

f[WolfBirthRateFactor_,WolfDeathRate_,MooseBirthRate_,MooseDeathRateFactor_,mo_,wo_]:=
Module[{lv,ratio,m,w,wolf,moose,t,residual,ratios,p1,p2},
lv=
NDSolve[
{(*  moose: *)
m'[t]==MooseBirthRate m[t]-MooseDeathRateFactor m[t]w[t] ,
(* wolves: *)
w'[t]==WolfBirthRateFactor m[t] w[t]-WolfDeathRate w[t],
m[0]==mo,
w[0]==wo
},
{m[t],w[t]},
{t,0,52}];
ratio[x_]:=First[First[MapThread[#1/.lv/.t->#2&,{({w[t]/m[t]}),{x}}]]];
p1=ListPlot[Transpose[{year,dataratio}]];
p2=Plot[ratio[t],{t,0,52}];
ratios=Table[ratio[i],{i,0,52}];
residual=Table[ratio[i]-dataratio[[i+1]],{i,0,52}];
Show[p2,p1];
residual.residual
]
f[WolfBirthRateFactor,WolfDeathRate,MooseBirthRate,MooseDeathRateFactor,mo,wo]

Cheap and dirty numerical derivatives: I’m in a hurry! This could be improved....

dt=.01;
f1[WolfBirthRateFactor_,WolfDeathRate_,MooseBirthRate_,MooseDeathRateFactor_,mo_,wo_]:=
(f[1.01WolfBirthRateFactor,WolfDeathRate,MooseBirthRate,MooseDeathRateFactor,mo,wo]-f[WolfBirthRateFactor,WolfDeathRate,MooseBirthRate,MooseDeathRateFactor,mo,wo])/(.01WolfBirthRateFactor)
f2[WolfBirthRateFactor_,WolfDeathRate_,MooseBirthRate_,MooseDeathRateFactor_,mo_,wo_]:=
(f[WolfBirthRateFactor,1.01WolfDeathRate,MooseBirthRate,MooseDeathRateFactor,mo,wo]-f[WolfBirthRateFactor,WolfDeathRate,MooseBirthRate,MooseDeathRateFactor,mo,wo])/(.01WolfDeathRate)
f3[WolfBirthRateFactor_,WolfDeathRate_,MooseBirthRate_,MooseDeathRateFactor_,mo_,wo_]:=
(f[WolfBirthRateFactor,WolfDeathRate,1.01MooseBirthRate,MooseDeathRateFactor,mo,wo]-f[WolfBirthRateFactor,WolfDeathRate,MooseBirthRate,MooseDeathRateFactor,mo,wo])/(.01MooseBirthRate)
f4[WolfBirthRateFactor_,WolfDeathRate_,MooseBirthRate_,MooseDeathRateFactor_,mo_,wo_]:=
(f[WolfBirthRateFactor,WolfDeathRate,MooseBirthRate,1.01MooseDeathRateFactor,mo,wo]-f[WolfBirthRateFactor,WolfDeathRate,MooseBirthRate,MooseDeathRateFactor,mo,wo])/(.01MooseDeathRateFactor)
f5[WolfBirthRateFactor_,WolfDeathRate_,MooseBirthRate_,MooseDeathRateFactor_,mo_,wo_]:=
(f[WolfBirthRateFactor,WolfDeathRate,MooseBirthRate,MooseDeathRateFactor,1.01mo,wo]-f[WolfBirthRateFactor,WolfDeathRate,MooseBirthRate,MooseDeathRateFactor,mo,wo])/(.01mo)
f6[WolfBirthRateFactor_,WolfDeathRate_,MooseBirthRate_,MooseDeathRateFactor_,mo_,wo_]:=
(f[WolfBirthRateFactor,WolfDeathRate,MooseBirthRate,MooseDeathRateFactor,mo,1.01wo]-f[WolfBirthRateFactor,WolfDeathRate,MooseBirthRate,MooseDeathRateFactor,mo,wo])/(.01wo)

Now for some starting values, define a gradient vector (roughly), and we’re ready to do gradient descent:

{WolfBirthRateFactor,WolfDeathRate,MooseBirthRate,MooseDeathRateFactor,mo,wo}={0.0004566167858993027`,0.4720627720562842`,0.13518440641048893`,0.004767210896693055`,805.2511688586087`,29.917958639847548`};
params={WolfBirthRateFactor,WolfDeathRate,MooseBirthRate,MooseDeathRateFactor,mo,wo}
objective=Apply[f,params]
gradrough[WolfBirthRateFactor_,WolfDeathRate_,MooseBirthRate_,MooseDeathRateFactor_,mo_,wo_]:=
{f1[WolfBirthRateFactor,WolfDeathRate,MooseBirthRate,MooseDeathRateFactor,mo,wo],f2[WolfBirthRateFactor,WolfDeathRate,MooseBirthRate,MooseDeathRateFactor,mo,wo],f3[WolfBirthRateFactor,WolfDeathRate,MooseBirthRate,MooseDeathRateFactor,mo,wo],f4[WolfBirthRateFactor,WolfDeathRate,MooseBirthRate,MooseDeathRateFactor,mo,wo],f5[WolfBirthRateFactor,WolfDeathRate,MooseBirthRate,MooseDeathRateFactor,mo,wo],f6[WolfBirthRateFactor,WolfDeathRate,MooseBirthRate,MooseDeathRateFactor,mo,wo]}
newparams=params-dt^3*Apply[gradrough,params]
newobjective=Apply[f,newparams]

Here’s where we do the steepest descent:

While[newobjective < objective,
objective=newobjective;
params=newparams;
newparams=params-dt^3*Apply[gradrough,params];
newobjective=Apply[f,newparams];
Print[newobjective, params]
]

Apply[f,{0.000203,0.25,0.25,0.0081,421.5,20.}]
Apply[f,params]
Apply[f,newparams]
params