{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Script to explore the Mathematics of Marriage project.\n",
    "\n",
    "#Here is the data from Table 3.11 for the male 1940-44 cohort, percentage\n",
    "#data rescaled to fractions:\n",
    "\n",
    "men4044 = [[0, 21.1/100], [5, 66.1/100], [10, 83.1/100], [15, 88.8/100], [20, 91.2/100], [25, 92.7/100], [30, 94.0/100]];\n",
    "N = len(men4044);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#A plot\n",
    "plt1 = scatter_plot(men4044,facecolor='red');\n",
    "show(plt1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#The function P(t) that might fit this data, according to the model, is\n",
    "var('t A b');\n",
    "P0 = men4044[0][1];\n",
    "P = function('P')(t,A,b)\n",
    "P(t,A,b) = P0/(P0+(1-P0)*exp(-A*(b^t-1)/log(b)));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Form a sum of squares to fit the data:\n",
    "SS = function('SS')(A,b);\n",
    "SS(A,b) = add((P(men4044[j][0],A,b)-men4044[j][1])^2 for j in range(N));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#A contour plot of log(SS(A,b))\n",
    "contour_plot(log(SS(A,b)),(A,0,1),(b,0,1),cmap='winter',fill=False,contours=20).show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Something near A = 0.6, b = 0.9 looks promising.\n",
    "#We could minimize using Calculus 3 techniques, but Sage has a\n",
    "#built-in command \"find_fit\" to find the least-squares model fit to the data:\n",
    "model(t) = P(t,A,b); #Specify the model to be fit, with \"t\" as the independent variable\n",
    "sol = find_fit(men4044,model,parameters=[A,b],initial_guess=[0.6,0.9]) #Fit the model by adjusting A and b\n",
    "Abest = sol[0].rhs();\n",
    "bbest = sol[1].rhs();\n",
    "bestP(t) = model(A=Abest,b=bbest) #Define the best fit function of the form in \"model\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Plot the function bestP(t) with this best choice for A and b.\n",
    "plt2 = plot(bestP(t), t, 0, 30);\n",
    "pp = plt1 + plt2;\n",
    "show(plt1+plt2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#The data for the 1945-49 men cohort:\n",
    "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]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#The data for the 1940-44 women cohort:\n",
    "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]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#The data for the 1945-49 women cohort:\n",
    "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]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 9.2",
   "language": "sage",
   "name": "sagemath"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}