Probability distribution function of the gestation period

Question

The problem is : The "gestation period" for humans averages 40 weeks, with a standard deviation of 10 days. Using CERN's root package, calculate and plot the probability that a woman will give birth tomorrow, given that she has made it all the way to today without giving birth yet, as a function of the day. You can start about two weeks before the due date, and go until two weeks afterward.

I am trying to solve this problem mathematically before solving it numerically, the approach I am taking is: Let the Guassian distribution give ideally the normalized number of women in a population N0 who give birth during each unit (day) with mean 40 weeks and standard deviation ten days ( probability distribution for large numbers of women of giving birth at 40 weeks with standard deviation 10 days). This can be represented by a normalized distribution N0 * ( 1/ 10 sqrt(2pi))* exp{ -((x- 280)/10) (sqaure)} [ 280 days is basically 40 weeks: the mean ]

The number who give birth during a period dx at time T is dN = N0 * (1/100 * Sqrt(2pi)) * exp ((T-280)/10) square dT

The number of women who have given birth until a time T is then: N(T) = N0 * Integral of the exponential distribution until T which is the error function of T with some modifications for the coefficients.

therefore the chance (probability)that a women will be giving birth at T having NOT done so up to that point is : dN(T)/ [ N0- N(T)] giving at T a function F(T)dT where F(T) is the probability asked for. Since the denominator is the likelihood that no birth has taken place as it is the number who have not given birth up to that point. N0 will cancel out.

I don’t know if my approach works because the distribution I should get must go to 1 after the mean because logically the probability of giving birth should not decrease. Is there a flaw in the logic above ? Any help would be appreciated.

Answer 1

One issue with assuming a normal distribution is that it is a continuous distribution while days are discrete. You need to decide early there is a integer time T where the probability of birth before $T$ is $\frac12$ and after $T$ is $\frac12$, or whether there is a due day indicated by $T$ which is more likely than any other day (in which case you may have the probability of birth before $T+\frac12$ being $\frac12$ and after $T+\frac12$ being $\frac12$). By choose a mean of $280$, you seem implicitly to have chosen the former approach.

N0 * (1/100 * Sqrt(2pi)) * exp ((T-280)/10) square dT looks odd and difficult to read. Perhaps you meant $N_0 (\frac1{100} \sqrt{2\pi}) \exp ((T-280)/10)^2$, which is not quite correct. I suspect you may have intended $N_0 \frac1{10 \sqrt{2\pi}} \exp\left(-\frac12(\frac{T-280}{10})^2\right)$

What you want to find is $$\mathbb P(t \lt T \le t+1 \mid t \lt T) = \dfrac{\mathbb P(t \lt T \le t+1)}{ \mathbb P(t \lt T)} \\ =\dfrac{F(t+1)-F(t)}{1-F(t)} =\dfrac{\int\limits_t^{t+1} f(x) \,dx}{\int\limits_t^{\infty} f(x) \,dx}\\ =\dfrac{\int\limits_t^{t+1} \frac1{10 \sqrt{2\pi}} \exp\left(-\frac12(\frac{x-280}{10})^2\right) \,dx}{\int\limits_t^{\infty} \frac1{10 \sqrt{2\pi}} \exp\left(-\frac12(\frac{x-280}{10})^2\right) \,dx}$$

This integration cannot be done by hand, but most computers can do it using preprogrammed functions. I do not know about CERN's ROOT package but in R you use something like

conditionalprobnextday <- function(t){
   (pnorm(t+1, 280, 10) - pnorm(t, 280, 10)) / (1 - pnorm(t, 280, 10))
   }

and perhaps a chart of this might look like