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}]]
moosedata[[1]]
wolfdata[[1]]

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]

params={0.0004566167858993027`,0.4720627720562842`,0.13518440641048893`,0.004767210896693055`,805.2511688586087`,29.917958639847548`}; (* nonlinear *)
params={0.0004582330863784787`,0.4751441941027219`,0.13430521341134866`,0.004720635218159999`,808.7832188526006`,30.048310913892962`}; (* from below *)
{WolfBirthRateFactor,WolfDeathRate,MooseBirthRate,MooseDeathRateFactor,mo,wo}=params;
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];
p1a=ListPlot[Transpose[{year,dataratio}]];
p2a=Plot[ratio[t],{t,0,52}];
Show[p1a,p2a]
p1a=ListPlot[Transpose[{year,moosedata}]];
p2a=Plot[Evaluate[{m[t]}/.lv],{t,0,52}];
Show[p1a,p2a]
p1a=ListPlot[Transpose[{year,wolfdata}]];
p2a=Plot[Evaluate[{w[t]}/.lv],{t,0,52}];
Show[p1a,p2a]

I got these parameters from optimization (Moose-n-Wolf-Insightmaker-BestFit.nb):

params={0.000240036,0.239826,0.269764,0.0113423,611.,20.};
params={0.0004566167858993027`,0.4720627720562842`,0.13518440641048893`,0.004767210896693055`,805.2511688586087`,29.917958639847548`}; (* nonlinear *)
{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}]]
p6=Apply[plotRatio,params];
Show[p2,p6,p1]
Show[p1,p6]

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

Clear[f]
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}}]]];
Table[ratio[i],{i,0,52}]
]
f[WolfBirthRateFactor,WolfDeathRate,MooseBirthRate,MooseDeathRateFactor,mo,wo]
estimates=Apply[f,params];
estimates

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

dt=10^(-3);
initialFactor=10^0;
sensitiveFactor=10^(-0);
f1[WolfBirthRateFactor_,WolfDeathRate_,MooseBirthRate_,MooseDeathRateFactor_,mo_,wo_]:=
(f[(1+sensitiveFactor dt)WolfBirthRateFactor,WolfDeathRate,MooseBirthRate,MooseDeathRateFactor,mo,wo]-f[WolfBirthRateFactor,WolfDeathRate,MooseBirthRate,MooseDeathRateFactor,mo,wo])/(sensitiveFactor dt WolfBirthRateFactor)
f2[WolfBirthRateFactor_,WolfDeathRate_,MooseBirthRate_,MooseDeathRateFactor_,mo_,wo_]:=
(f[WolfBirthRateFactor,(1+dt)WolfDeathRate,MooseBirthRate,MooseDeathRateFactor,mo,wo]-f[WolfBirthRateFactor,WolfDeathRate,MooseBirthRate,MooseDeathRateFactor,mo,wo])/(dt WolfDeathRate)
f3[WolfBirthRateFactor_,WolfDeathRate_,MooseBirthRate_,MooseDeathRateFactor_,mo_,wo_]:=
(f[WolfBirthRateFactor,WolfDeathRate,(1+dt)MooseBirthRate,MooseDeathRateFactor,mo,wo]-f[WolfBirthRateFactor,WolfDeathRate,MooseBirthRate,MooseDeathRateFactor,mo,wo])/(dt MooseBirthRate)
f4[WolfBirthRateFactor_,WolfDeathRate_,MooseBirthRate_,MooseDeathRateFactor_,mo_,wo_]:=
(f[WolfBirthRateFactor,WolfDeathRate,MooseBirthRate,(1+dt)MooseDeathRateFactor,mo,wo]-f[WolfBirthRateFactor,WolfDeathRate,MooseBirthRate,MooseDeathRateFactor,mo,wo])/(dt MooseDeathRateFactor)
f5[WolfBirthRateFactor_,WolfDeathRate_,MooseBirthRate_,MooseDeathRateFactor_,mo_,wo_]:=
(f[WolfBirthRateFactor,WolfDeathRate,MooseBirthRate,MooseDeathRateFactor,(1+initialFactor dt)mo,wo]-f[WolfBirthRateFactor,WolfDeathRate,MooseBirthRate,MooseDeathRateFactor,mo,wo])/(initialFactor dt mo)
f6[WolfBirthRateFactor_,WolfDeathRate_,MooseBirthRate_,MooseDeathRateFactor_,mo_,wo_]:=
(f[WolfBirthRateFactor,WolfDeathRate,MooseBirthRate,MooseDeathRateFactor,mo,(1+initialFactor  dt)wo]-f[WolfBirthRateFactor,WolfDeathRate,MooseBirthRate,MooseDeathRateFactor,mo,wo])/(initialFactor dt wo)

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.00020621277957375955`,0.19466437259371702`,0.2975792202921484`,0.012347035565225824`,660.0000000006974`,20.00000026558482`};
params={WolfBirthRateFactor,WolfDeathRate,MooseBirthRate,MooseDeathRateFactor,mo,wo}
Map[Norm,{Apply[f1,params],Apply[f2,params],Apply[f3,params],Apply[f4,params],Apply[f5,params],Apply[f6,params]}]

Okay, now we iterate:

params={0.0004234318101923372`,0.4714892395462701`,0.13535046346199842`,0.00443178985381275`,867.0801972243067`,32.21605332763454`}
Norm[dataratio-Apply[f,params]]^2
runs=1;
Do[
yhat=Apply[f,params];
w=dataratio-yhat;
diffLinMod=LinearModelFit[Transpose[
{Apply[f1,params],Apply[f2,params],Apply[f3,params],Apply[f4,params],Apply[f5,params],Apply[f6,params],w}],{a,b,c,x,y,z},{a,b,c,x,y,z},IncludeConstantBasis->False];
newparams=diffLinMod["BestFitParameters"];
params=params+newparams;
Print[diffLinMod["ParameterTable"]],
{i,runs}
]
params
Norm[dataratio-Apply[f,params]]

params={0.000267702,0.239817,0.269739,0.0113796,591.,23.}; (* 0.082952946598904` *)
params={0.00023195,0.239829,0.269769,0.0113139,660.,20.} (* 0.08294549472838839` *)

params={0.0004234318101923372`,0.4714892395462701`,0.13535046346199842`,0.00443178985381275`,moosedata[[1]],wolfdata[[1]]}
params={0.000240036,0.239826,0.269764,0.0113423,moosedata[[1]],wolfdata[[1]]}
params={0.000456813497789625`,0.4724329986468058`,0.13507930222204007`,0.004761597649195472`,805.6719414177632`,29.93405848015483`}
objective=1000;
newobjective=Norm[dataratio-Apply[f,params]];
newobjective^2
While[newobjective < objective,
objective=newobjective;
Do[
tmp=objective;
yhat=Apply[f,params];
w=dataratio-yhat;
fs={f1,f2,f3,f4,f5,f6};
diffLinMod=LinearModelFit[Transpose[
{Apply[fs[[j]],params],w}],{x},{x},IncludeConstantBasis->False];
Print[diffLinMod["ParameterTable"]];
fitparams=diffLinMod["BestFitParameters"];
params[[j]]=params[[j]]+fitparams[[1]],
{j,6}
];
newobjective=Norm[dataratio-Apply[f,params]];
Print[newobjective," ",objective," ",params]
]
params
Norm[dataratio-Apply[f,params]]^2

params
Norm[dataratio-Apply[f,params]]