Thanks Claude and Luca. That clears things up. So the actual (continuous) time of an egg hatching is only recorded by day. That means, for example, the 8 eggs recorded on day 9 hatched sometime between day 8 and day 9. To account for the fact that you have "censored" data one needs to construct the likelihood (or the log likelihood) a bit different than what one would do if there were continuous measurements.

Also, Mathematica has a built-in Erlang distribution as Claude has shown. However, Mathematica uses a slightly different parameterization that you use. But the conversion is simple: to use Mathematica's built-in Erlang distribution you would use

ErlangDistribution[k,k/d]

To construct the likelihood function so that maximum likelihood estimates of $k$ and $d$ could be obtained, you would use the following:

where $G(x)$ is the cumulative distribution function of the Erlang distribution and $f_i$ is the associated number of eggs counted on day $i$. The Mathematica code is

data = {{8, 9}, {9, 12}, {10, 16}, {11, 20}, {12, 25}, {13, 30}, {14, 35}, {15, 39}, {16, 43}, {17, 47},
{18, 49}, {19, 51}, {20, 52}, {21, 51}, {22, 50}, {23, 48}, {24, 46}, {25, 43}, {26, 40}, {27, 36}, {28, 33},
{29, 29}, {30, 26}, {31, 22}, {32, 19}, {33, 17}, {34, 14}, {35, 12}, {36, 10}, {37, 8}};
n = 1000;
count = data[[All, 2]];

logLikelihood = count[[1]] Log[CDF[ErlangDistribution[k, k/d], 8]] + 
  Sum[count[[i]] Log[(CDF[ErlangDistribution[k, k/d], i + 7] - CDF[ErlangDistribution[k, k/d], i + 6])], {i, 2, Length[data]}] +
  (n - Total[count]) Log[(1 - CDF[ErlangDistribution[k, k/d], 37])];

sol = NMaximize[{logLikelihood, d > 10 && k \[Element] Integers && k > 1}, {k, d}, MaxIterations -> 1000]
(* {-3308.429332485035`,{k -> 7,d -> 22.53458409510396}} *)

The maximum likelihood estimates do not match what you expect but let's plot the data and the fits with maximum likelihood estimates and the "expected" parameters.

Show[Plot[{PDF[ErlangDistribution[k, k/d] /. sol[[2]], x],
   PDF[ErlangDistribution[k, k/d] /. {k -> 20, d -> 35}, x]}, {x, 8, 37},
   PlotStyle -> {Blue, Red}, AxesOrigin -> {7, 0},
   PlotLegends -> {"k -> 7, d -> 22.5346", "k -> 20, d -> 35"}],
 ListPlot[Transpose[{data[[All, 1]], data[[All, 2]]/1000}]]]

So I think we're closer but not quite there yet as this doesn't match with your expected values for $k$ and $d$.

Thank You Claude, Each pair of data is {time; born}: if I consider a number of 1000 eggs laid in day 0, everyday I control if there are some hatching eggs, and I count it. For "born" I mean, according this experimental situation how many hatching eggs I register everyday, and this representation, should give an Erlang distribution. Parameter k is the number of substage not observable that every individual has to overtake to became a larva (hatching eggs), and d is the mean develop time for population. At this point, n, the constant considered by Jim, is the total number of individual who composed the initial cohort.

I hope to cleared the experimental situation now. I have built a model in excel who reproduces this real situation, not considering the mortality, and it seems to provides an Erlang, but I need to fit it to confront the parameter k and d (they must be equal if the numerical model run well).

According to my data, I have to obtain k=20 and d=35

What you say is partially true: my Erlang have an integer N=1000 who represent the initial cohort entering in the process. At time 0 a cohort composed from N=1000 individuals enter in the process, and then exit distributed as an Erlang function. N could be another parameter, but according to this question, I can divide for 1000 the fitting function, solving (partially) the problem.

But there are still some questions open: if I "normalize" the fitting function dividing for 1000, I should obtain my original parameter k=20 and d=35, but it seems to not happen, why?

Thank you really for your help!

Thank you Claude Mante and Jim Baldwind for your interest. Before i firget to insert the .txt with data that I'm using for the fit. It represents data from a cohort who develop itself in time (insects/day), and theoretically they have to represent an Erlang distribution (or Gamma, the problem is similar). When I tried to analyze experimental data, I noticed this error in Mathematica, and first to proceed I investigated the problem. To understand if my experimental data had something wrong, I built with excel two columns x, f(x) who represents the coordinates of an Erlang distribution, of which I know parameters. If mathematica works well, with NonlinearModelFit I have to obtain back initial parameters.

According this theory, I have a curve that resembles a probability density function, and not a random samples of a probability distribution. I hope that .txt can help you to understand the problem. Thank you everybody for the interest and to help me

data = {{8, 9}, {9, 12}, {10, 16}, {11, 20}, {12, 25}, {13, 30}, {14, 35}, \
{15, 39}, {16, 43}, {17, 47}, {18, 49}, {19, 51}, {20, 52}, {21, 51}, \
{22, 50}, {23, 48}, {24, 46}, {25, 43}, {26, 40}, {27, 36}, {28, 33}, \
{29, 29}, {30, 26}, {31, 22}, {32, 19}, {33, 17}, {34, 14}, {35, 12}, \
{36, 10}, {37, 8}}

Attachments:

It doesn't work, this is the answer:

FindDistributionParameters::ntsprt: One or more data points are not in support of the process or distribution GammaDistribution[k,a]. >>

In my NonlinearModelFit code, is there any error? I don't know what is happening