laplace_calculations.ipynb

{
 "metadata": {
  "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.6.3-final"
  },
  "orig_nbformat": 2,
  "kernelspec": {
   "name": "python36364bitpython3condac51aacc59df947699a98ee01d7b768f0",
   "display_name": "Python 3.6.3 64-bit ('python3': conda)",
   "language": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2,
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [],
   "source": [
    "class function:\n",
    "    def __init__(self,function_as_real_map):\n",
    "        self.real_map = function_as_real_map\n",
    "    def function_as_discrete_map(self, start, end, points_num):\n",
    "        return [self.real_map(start+i*(end-start)/(points_num-1)) for i in range(points_num)]\n",
    "    def function_as_discrete_map_extra_points(self, start, end, points_num):\n",
    "        return self.function_as_discrete_map(start, end+5, points_num+5)\n",
    "\n",
    "def laplace_integration(list):\n",
    "    k_7_0 = 41393.0/60480.0\n",
    "    k_7_1 = -23719.0/60480.0\n",
    "    k_7_2 = 22742.0/60480.0\n",
    "    k_7_3 = -14762.0/60480.0\n",
    "    k_7_4 = 5449.0/60480.0\n",
    "    k_7_5 = -863.0/60480.0\n",
    "\n",
    "    k_7 = [k_7_0, k_7_1, k_7_2, k_7_3, k_7_4, k_7_5]\n",
    "    \n",
    "    b_coor = len(list) - 6\n",
    "    result = sum([list[i] for i in range(0,b_coor)])\n",
    "    result += sum([k_7[i] * list[b_coor + i] for i in range(0,6)])\n",
    "    result -= sum([k_7[i] * list[i] for i in range(0,6)])\n",
    "    \n",
    "    return result\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_squared = function(lambda x: x*x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": [
       "[0.0,\n",
       " 1.0,\n",
       " 4.0,\n",
       " 9.0,\n",
       " 16.0,\n",
       " 25.0,\n",
       " 36.0,\n",
       " 49.0,\n",
       " 64.0,\n",
       " 81.0,\n",
       " 100.0,\n",
       " 121.0,\n",
       " 144.0,\n",
       " 169.0,\n",
       " 196.0,\n",
       " 225.0]"
      ]
     },
     "metadata": {},
     "execution_count": 85
    }
   ],
   "source": [
    "grid_values = x_squared.function_as_discrete_map_extra_points(0.0,10.0,11)\n",
    "grid_values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": [
       "333.33333333333337"
      ]
     },
     "metadata": {},
     "execution_count": 86
    }
   ],
   "source": [
    "laplace_integration(grid_values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": [
       "(333.33333333333326, 3.700743415417188e-12)"
      ]
     },
     "metadata": {},
     "execution_count": 87
    }
   ],
   "source": [
    "import scipy.integrate as integrate\n",
    "integrate.quad(x_squared.real_map, 0, 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.integrate as integrate\n",
    "\n",
    "def compare_integration_methods(integrand_as_lambda,start,end,num_point_evals):\n",
    "    f = function(integrand_as_lambda)\n",
    "    grid_values = f.function_as_discrete_map_extra_points(start,end,num_point_evals)\n",
    "    laplace_int = laplace_integration(grid_values)\n",
    "    theoretical_int, theoretical_dev = integrate.quad(f.real_map, start, end)\n",
    "    abs_dev = np.abs(theoretical_int-laplace_int)\n",
    "    rel_dev = np.abs(theoretical_int-laplace_int)/theoretical_int\n",
    "    print(\"Percentage deviation:\", rel_dev)\n",
    "    print(\"Absolute deviation:\", abs_dev)\n",
    "    print(\"Uncertainty in theoretical value:\", theoretical_dev)\n",
    "    print(\"--------------------------------------------------\")\n",
    "\n",
    "def compare_integration_methods_verbose(integrand_as_lambda,start,end,num_point_evals):\n",
    "    f = function(integrand_as_lambda)\n",
    "    grid_values = f.function_as_discrete_map_extra_points(start,end,num_point_evals)\n",
    "    laplace_int = laplace_integration(grid_values)\n",
    "    theoretical_int, theoretical_dev = integrate.quad(f.real_map, start, end)\n",
    "    abs_dev = np.abs(theoretical_int-laplace_int)\n",
    "    rel_dev = np.abs(theoretical_int-laplace_int)/theoretical_int\n",
    "    print(\"Theoretical value:\", theoretical_int)\n",
    "    print(f\"Interval from {start} to {end}\")\n",
    "    print(\"Percentage deviation:\", rel_dev)\n",
    "    print(\"Absolute deviation:\", abs_dev)\n",
    "    print(\"Uncertainty in theoretical value:\", theoretical_dev)\n",
    "    print(\"--------------------------------------------------\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "Percentage deviation: 3.41060513165e-16\nAbsolute deviation: 1.13686837722e-13\nUncertainty in theoretical value: 3.700743415417188e-12\n--------------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "compare_integration_methods(lambda x: x*x,0.0,10.0,11)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": [
       "[1.0,\n",
       " 0.7984992210323577,\n",
       " 0.6361919882852219,\n",
       " 0.5039718605924449,\n",
       " 0.3933165453538113,\n",
       " 0.29528485238517843,\n",
       " 0.20045054639867776,\n",
       " 0.10475631057767379,\n",
       " 0.027752547646899984,\n",
       " 0.001259078477171793,\n",
       " 5.83717348811404e-07,\n",
       " 1.8277053679093564e-15]"
      ]
     },
     "metadata": {},
     "execution_count": 90
    }
   ],
   "source": [
    "# Example with annuity of a zero year old - intentisty constant values from Norberg eq. (3.25)\n",
    "interest_rate = 0.02 # arbitrary choice \n",
    "alpha = 0.0005\n",
    "beta = 7.5858e-05\n",
    "gamma = np.math.log(1.09144)\n",
    "p_survival_x_year_old = lambda t,x: np.math.exp(-alpha*t-beta*(np.math.exp(gamma*(t+x))-1)/gamma)\n",
    "p_survival_zero_year_old = lambda t: p_survival_x_year_old(t,0)\n",
    "annuity_x_year_old_integrand = lambda t,x: np.math.exp(-1*interest_rate*t)*p_survival_x_year_old(t,x)\n",
    "annuity_zero_year_old_integrand = lambda t: np.math.exp(-1*interest_rate*t)*p_survival_zero_year_old(t)\n",
    "\n",
    "function(annuity_zero_year_old_integrand).function_as_discrete_map(0,120,12)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "Percentage deviation: 3.78374019353e-15\nAbsolute deviation: 1.42108547152e-13\nUncertainty in theoretical value: 3.395439631872089e-10\n--------------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "compare_integration_methods(annuity_zero_year_old_integrand,0.0,120.0,121)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "x cubed\nPercentage deviation: 1.81898940355e-16\nAbsolute deviation: 4.54747350886e-13\nUncertainty in theoretical value: 2.775557561562892e-11\n--------------------------------------------------\nx to the power 4\nPercentage deviation: 0.0\nAbsolute deviation: 0.0\nUncertainty in theoretical value: 2.220446049250313e-10\n--------------------------------------------------\nx to the power 5\nPercentage deviation: 0.0\nAbsolute deviation: 0.0\nUncertainty in theoretical value: 1.8503717077085944e-09\n--------------------------------------------------\nx to the power 6\nPercentage deviation: 1.62981450558e-16\nAbsolute deviation: 2.32830643654e-10\nUncertainty in theoretical value: 1.5860328923216523e-08\n--------------------------------------------------\nx to the power 7\nPercentage deviation: 4.58333333334e-05\nAbsolute deviation: 572.916666668\nUncertainty in theoretical value: 1.387778780781446e-07\n--------------------------------------------------\nx to the power 8\nPercentage deviation: 0.000296047\nAbsolute deviation: 32894.1111111\nUncertainty in theoretical value: 1.2335811384723963e-06\n--------------------------------------------------\nexponential\nPercentage deviation: 0.122354406265\nAbsolute deviation: 2694.91279003\nUncertainty in theoretical value: 6.239389118119916e-10\n--------------------------------------------------\nexponential\nPercentage deviation: 0.0061622947631\nAbsolute deviation: 0.0113329008526\nUncertainty in theoretical value: 9.92116330910929e-12\n--------------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "# Examples with basic functions\n",
    "print(\"x cubed\")\n",
    "compare_integration_methods(lambda x: x**3,0.0,10.0,11)\n",
    "print(\"x to the power 4\")\n",
    "compare_integration_methods(lambda x: x**4,0.0,10.0,11)\n",
    "print(\"x to the power 5\")\n",
    "compare_integration_methods(lambda x: x**5,0.0,10.0,11)\n",
    "print(\"x to the power 6\")\n",
    "compare_integration_methods(lambda x: x**6,0.0,10.0,11)\n",
    "print(\"x to the power 7\")\n",
    "compare_integration_methods(lambda x: x**7,0.0,10.0,11)\n",
    "print(\"x to the power 8\")\n",
    "compare_integration_methods(lambda x: x**8,0.0,10.0,11)\n",
    "print(\"exponential\")\n",
    "compare_integration_methods(lambda x: np.exp(x),0.0,10.0,11)\n",
    "print(\"sine\")\n",
    "compare_integration_methods(lambda x: np.sin(x),0.0,10.0,11)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "For a 80 year old doing the integral from 0 to 120\nTheoretical value: 2.458672942142802\nInterval from 0.0 to 120.0\nPercentage deviation: 3.37622794691e-09\nAbsolute deviation: 8.30104029959e-09\nUncertainty in theoretical value: 2.025232637842267e-10\n--------------------------------------------------\nFor a 80 year old doing the integral from 0 to 30\nTheoretical value: 2.458672075352603\nInterval from 0.0 to 30.0\nPercentage deviation: 4.16876335106e-09\nAbsolute deviation: 1.024962204e-08\nUncertainty in theoretical value: 1.0727505925202197e-12\n--------------------------------------------------\nFor a 100 year old doing the integral from 0 to 120\nTheoretical value: 0.007348811477706047\nInterval from 0.0 to 120.0\nPercentage deviation: 5.11575113186e-06\nAbsolute deviation: 3.75946906349e-08\nUncertainty in theoretical value: 4.145520167647774e-09\n--------------------------------------------------\nFor a 100 year old doing the integral from 0 to 30\nTheoretical value: 0.007348811477706047\nInterval from 0.0 to 30.0\nPercentage deviation: 5.11575113186e-06\nAbsolute deviation: 3.75946906349e-08\nUncertainty in theoretical value: 4.145520166947887e-09\n--------------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "# Examples with different x for annuities\n",
    "print(\"For a 80 year old doing the integral from 0 to 120\")\n",
    "compare_integration_methods_verbose(lambda t: annuity_x_year_old_integrand(t,80),0.0,120.0,121)\n",
    "print(\"For a 80 year old doing the integral from 0 to 30\")\n",
    "compare_integration_methods_verbose(lambda t: annuity_x_year_old_integrand(t,80),0.0,30.0,31)\n",
    "print(\"For a 100 year old doing the integral from 0 to 120\")\n",
    "compare_integration_methods_verbose(lambda t: annuity_x_year_old_integrand(t,100),0.0,120.0,121)\n",
    "print(\"For a 100 year old doing the integral from 0 to 30\")\n",
    "compare_integration_methods_verbose(lambda t: annuity_x_year_old_integrand(t,100),0.0,30.0,31)"
   ]
  }
 ]
}