%% The Mathematics of Marriage
%% Kurt Bryan and SIMIODE
%
% Here is the data from Table 3.11 for the male 1940-44 cohort, percentage data rescaled to fractions:

men4044 = [0.211, 0.661, 0.831, 0.888, 0.912, 0.927, 0.940];
times = [0 5 10 15 20 25 30];
N = length(times);

% A plot of the data

plot(times,men4044)

% The function P(t) that might fit this data, according to the model, is

syms t;
syms A;
syms b;
P0 = men4044(1);
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(A,b)=sum((P(times,A,b)-men4044).^2);

% Now minimize in A and b. A contour plot of log(SS) may be helpful.

fcontour(log(SS(A,b)), [0 1 0 1])

% Something near A = 0.6, b = 0.9 looks promising. We can set dSS/dA = 0 and dSS/db = 0 to find this point.

dSSdA = diff(SS(A,b),A);
dSSdb = diff(SS(A,b),b);
vpasolve([dSSdA==0,dSSdb==0],[A,b],[0.6,0.9]);
fprintf("A = %f\n",Absol.A)
fprintf("b = %f\n",Absol.b)

% Plot P(t) with these values, compare to the data

plot(times,men4044)
hold on
fplot(P(t,Absol.A,Absol.b),[0 30],'-r')

% Here is the data for the 1945-49 men

men4549 = [0.223, 0.655, 0.801, 0.861, 0.893, 0.913, 0.925]

% Here is the data for the 1940-44 women:

women4044 = [0.481, 0.782, 0.868, 0.897, 0.914, 0.925, 0.932]

% Here is the data for the 1945-49 women:

women4549=[0.431, 0.769, 0.850, 0.884, 0.902, 0.915, 0.922]