Efficient Table Generation & NDSolve Error

I'm currently working on solving an SEIR model with data from the CDC. I've created a function to calculate the Square difference for the model and want to create a table of values for the 4 variable combinations, along with their Goodness of fit characteristic from the square difference.

I'm getting this alert when running the function:
"NDSolve has computed initial values that give a zero residual for the differential-algebraic system, but some components are different from those specified. If you need them to be satisfied, giving initial conditions for all dependent variables and their derivatives is recommended."

It sounds like something I could ignore, but I don't want to mess anything up by doing so.

Also... the big question. Is there any recommendations to make creating this table of values more efficient rather than waiting a half hour for this to work itself out?

ClearAll["Global*"]

CDCdata = {{0, 161}, {1, 197}, {2, 298}, {3, 339}, {4, 409}, {5,
614}, {6, 848}, {7, 1359}, {8, 1985}, {9, 3346}, {10, 4590}, {11,
5938}, {12, 6511}, {13, 6514}, {14, 6952}, {15, 5982}, {16,
5167}, {17, 4077}, {18, 3126}, {19, 2457}, {20, 2136}, {21,
1848}, {22, 1680}, {23, 1497}, {24, 1269}, {25, 1007}, {26,
898}, {27, 644}, {28, 459}, {29, 380}};

ListPlot[CDCdata]

p = Interpolation[CDCdata];
Plot[p[t], {t, 0, 29}]

Population = 20328;
Squarediff[Mu_, Beta_, Alpha_,
Nu_] := (sol =
First@NDSolve[{s'[t] ==
Mu*Population - Mu*s[t] - Beta*(i[t]/Population)*s[t],
e'[t] == Beta*(i[t]/Population)*s[t] - (Mu + Alpha) e[t],
i'[t] == Alpha*(e[t]) - (Nu + Mu)*(i[t]),
r'[t] == Nu*i[t] - Mu*r[t],
i[0] == 161,
s[0] == 20328,
r[0] == 0, e[0] == 0},
{s, e, i, r}, {t, 0, 29}];
Sum[(((CDCdata[[j, 2]] - i[CDCdata[[j, 1]]]) /. sol)^2), {j, 1,
Length[CDCdata]}])

values = Flatten[
Table[{{a, b, c, d}, Squarediff[a, b, c, d]}, {a, 0, 0.0005,
0.00001}, {b, 0, 0.2, 0.01}, {c, 0, 0.1, 0.01}, {d, 0, 0.1,
0.01}], 1];

-
What about Beta? – Artes Nov 10 '13 at 4:21
Don't create symbols that start with a capital letter as those might shadow/conflict with Mathematica defined symbols. In your case Beta is a Mathematica defined function! – Matariki Nov 10 '13 at 6:19
@Artes are you trying to say the same thing as the comment below? – Kevin Murphy Nov 10 '13 at 14:26

EDIT

ParametricNDSolveValue may be a useful way to explore the parameter space and then you can quantify goodness of fit. Using the same definitions as your original post (with change of variables):

fun = ParametricNDSolveValue[{s'[t] ==
mu*Population - mu*s[t] - beta*(i[t]/Population)*s[t],
e'[t] == beta*(i[t]/Population)*s[t] - (mu + alpha) e[t],
i'[t] == alpha*(e[t]) - (nu + mu)*(i[t]),
r'[t] == nu*i[t] - mu*r[t], i[0] == 161, s[0] == 20328, r[0] == 0,
e[0] == 0}, {s, e, i, r}, {t, 0, 29}, {mu, beta, alpha, nu}]


This creates set of interpolating functions for: susceptible, infected etc. As you are interested in i[t] this is the third function. You can dynamically explore the parameter space ( I have put the funcion p[t] as the reference):

Manipulate[Module[{foi},
foi = fun[a, b, c, d][[3]];
Plot[{foi[t], p[t]}, {t, 0, 29}]],
{{a, 0.001, \[Mu]}, 0.001, 1, Appearance -> "Labeled"},
{{b, 0.01, \[Beta]}, 0.01, 1, Appearance -> "Labeled"},
{{c, 0.01, \[Alpha]}, 0.01, 1, Appearance -> "Labeled"},
{{d, 0.01, \[Nu]}, 0.01, 1, Appearance -> "Labeled"}
]


Not rigorous but a starting point. This may narrow the parameter space and you can calculate or use other methods to find best fit.

-
Thanks for your answer! I've changed the variable names. The issue I have now is, I want to get all 4 variable names a bit more accurate. By changing the step size in the values` table, I increase the calculation time by a ton. It ran over night and still isn't done. Any recommendations? – Kevin Murphy Nov 10 '13 at 14:25
@KevinMurphy see edit – ubpdqn Nov 11 '13 at 1:52