#The function P(t) that might fit this data, according to the model, is
var('t A b');
P0 = men4044[0][1];
P = function('P')(t,A,b)
P(t,A,b) = P0/(P0+(1-P0)*exp(-A*(b^t-1)/log(b)));

#Form a sum of squares to fit the data:
SS = function('SS')(A,b);
SS(A,b) = add((P(men4044[j][0],A,b)-men4044[j][1])^2 for j in range(N));

#A contour plot of log(SS(A,b))
contour_plot(log(SS(A,b)),(A,0,1),(b,0,1),cmap='winter',fill=False,contours=20).show()

#Something near A = 0.6, b = 0.9 looks promising.
#We could minimize using Calculus 3 techniques, but Sage has a
#built-in command "find_fit" to find the least-squares model fit to the data:
model(t) = P(t,A,b); #Specify the model to be fit, with "t" as the independent variable
sol = find_fit(men4044,model,parameters=[A,b],initial_guess=[0.6,0.9]) #Fit the model by adjusting A and b
Abest = sol[0].rhs();
bbest = sol[1].rhs();
bestP(t) = model(A=Abest,b=bbest) #Define the best fit function of the form in "model"

#Plot the function bestP(t) with this best choice for A and b.
plt2 = plot(bestP(t), t, 0, 30);
pp = plt1 + plt2;
show(plt1+plt2)

#The data for the 1945-49 men cohort:
men4549 = [[0, 22.3/100], [5, 65.5/100], [10, 80.1/100], [15, 86.1/100], [20, 89.3/100], [25, 91.3/100], [30, 92.5/100]]

#The data for the 1940-44 women cohort:
women4044 = [[0, 48.1/100], [5, 78.2/100], [10, 86.8/100], [15, 89.7/100], [20, 91.4/100], [25, 92.5/100], [30,
93.2/100]]

#The data for the 1945-49 women cohort:
women4549 = [[0, 43.1/100], [5, 76.9/100], [10, 85.0/100], [15, 88.4/100], [20, 90.2/100], [25, 91.5/100], [30,
92.2/100]]