homework_1_RMorehead.ipynb

{
 "metadata": {
  "language": "Julia",
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "##Robert Morehead\n",
      "##Homework 1"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "##1.\n",
      "I got julia, julia studio and ijulia + notebook working on my system. I don't plan on using the lab computer.  "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "##2a."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#Code from exercise\n",
      "srand(42)\n",
      "N = 10000\n",
      "true_mean = 10000.\n",
      "y = true_mean+randn(N)\n",
      "\n",
      "sample_mean = mean(y)\n",
      "sample_var = var(y)\n",
      "println(\"Mean: \",sample_mean,\"  Varience: \",sample_var)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Mean: 9999"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        ".998011659798  Varience: 1.0199227833780644\n"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "##2b."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "function mean_and_var(y,sample_mean,sample_var) #made this a function for part 2d.\n",
      "y32bit = convert(Array{Float32,1},y);\n",
      "y16bit = convert(Array{Float16,1},y);\n",
      "mean32,var32 = mean(y32bit),var(y32bit)\n",
      "mean16,var16 = mean(y16bit),var(y16bit)\n",
      "\n",
      "println(\"64 bit mean: \",sample_mean,\" var: \",sample_var)\n",
      "println()\n",
      "println(\"32 bit mean: \",mean32,\" var: \",var32)\n",
      "println(\"dif mean: \",sample_mean-mean32,\" dif var: \",sample_var-var32)\n",
      "println(\"error mean: \",(sample_mean-mean32)/sample_mean,\" error var: \",(sample_var-var32)/sample_var)\n",
      "println()\n",
      "println(\"16 bit mean: \",mean16,\" var: \",var16)\n",
      "println(\"dif mean: \",sample_mean-mean16,\" dif var: \",sample_var-var16)\n",
      "println(\"error mean: \",(sample_mean-mean16)/sample_mean,\" error var: \",(sample_var-var16)/sample_var)\n",
      "end \n",
      "\n",
      "mean_and_var(y,sample_mean,sample_var)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "64 bit mean: "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "9999.998011659798 var: 1.0199227833780644\n",
        "\n",
        "32 bit mean: 9999.998 var: 1.0199206\n",
        "dif mean: -3.521520193316974e-5 dif var: 2.1958383915610113e-6\n",
        "error mean: -3.52152089351513e-9 error var: 2.152945720349752e-6\n",
        "\n",
        "16 bit mean: 9995.983 var: 16.001316\n",
        "dif mean: 4.014613222298067 dif var: -14.981393287178577\n",
        "error mean: 0.0004014614020539912 error var: -14.688752453944625\n"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The 32 bit array has a \"small\" difference compared to the 64 bit array (e-5, e-6), while the 16 bit array results in a very large difference. The error is increasing in the varience calculation faster because the rounding errors are being squared, also there may be catastrophic cancellation in the the varience calculation (it contains a subtraction) whin the random error is small.  "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "##2c.\n",
      "Increasing N to $10^5$ should increase the errors in both 32 and 16 bits cases as there are more evaluations for rounding errors to accululate.\n",
      "\n",
      "The relative error will be smaller in the $\\mu = 10^5$ case, because $\\mu$ is larger. The absolute difference should stay the same. "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "##2d."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#CASE 1\n",
      "#Code from exercise\n",
      "srand(42)\n",
      "N = 100000\n",
      "true_mean = 10000.\n",
      "y = true_mean+randn(N)\n",
      "\n",
      "sample_mean = mean(y)\n",
      "sample_var = var(y)\n",
      "\n",
      "\n",
      "mean_and_var(y,sample_mean,sample_var)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "64 bit mean: 10000.0031598515 var: 1.0088883985077686\n",
        "\n",
        "32 bit mean: 10000.003 var: 1.0088867\n",
        "dif mean: 0.0002301639997313032 dif var: 1.7035996264791464e-6\n",
        "error mean: 2.301639270029202e-8 error var: 1.6885907588975298e-6\n",
        "\n",
        "16 bit mean: 9996.01 var: 16.000065\n",
        "dif mean: 3.9933942264997313 dif var: -14.991176451345748\n",
        "error mean: 0.00039933929646468567 error var: -14.859102823978318\n"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#CASE 2\n",
      "#Code from exercise\n",
      "srand(42)\n",
      "N = 10000\n",
      "true_mean = 100000.\n",
      "y = true_mean+randn(N)\n",
      "\n",
      "sample_mean = mean(y)\n",
      "sample_var = var(y)\n",
      "\n",
      "\n",
      "mean_and_var(y,sample_mean,sample_var)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "64 bit mean: 99999.9980116598 var: 1.0199227833780122\n",
        "\n",
        "32 bit mean: 100000.0 var: 1.0200287\n",
        "dif mean: -0.001988340198295191 dif var: -0.00010592698728317806\n",
        "error mean: -1.9883402378301593e-8 error var: -0.00010385784983873483\n",
        "\n",
        "16 bit mean: Inf var: NaN\n",
        "dif mean: -Inf dif var: NaN\n",
        "error mean: -Inf error var: NaN\n"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Ah, nice one! In case two it looks like we hit an overflow in the $\\mu$ calculation. I will blame always working in double precision for not seeing that one comming. \n",
      "\n",
      "In case one, the errors slightly improved, rather than getting worse. Maybe the rounding is calleling out, and I got lucky? A quick change in random seed doesn't see to support this. The law of large numbers at play? Maybe, the floating point errors are just another source of random error. So yeah I a going to go with that. \n",
      "\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "##2e.\n",
      "1. Be aware of the range of values allowed by your precision level.\n",
      "2. Be mindful of floating point arithmitic when designing algorithms. \n",
      "3. When working with random samples, bigger is better (most of the time)  "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "##3a."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "srand(42)\n",
      "N = 1000000\n",
      "true_mean = 1e6\n",
      "y = true_mean+randn(N);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 41
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$m = 1/N \\times \\sum_{i=1}^{N} y_i$ \n",
      "\n",
      "$s^2 = 1/(N-1) \\times \\sum_{i=1}^N (y_i-m)^2 = 1/(N-1)  \\times \\left[ \\left( \\sum_{i=1}^N y_i^2 \\right) - N m^2 \\right]$\n",
      "\n",
      "Note: A quick glance at http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance cleared up my confusion on what was meant by \"one-pass\" vs. \"two-pass\". "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "function onepass(y::Array)\n",
      "\n",
      "  N = length(y)\n",
      "  sUm = zero(y[1]) #oopse this killed the built in sum function\n",
      "  sUm_squared = zero(y[1])\n",
      "  \n",
      "  for i in 1:N\n",
      "    sUm = sUm + y[i]\n",
      "    sUm_squared = sUm_squared + y[i]^2\n",
      "  end \n",
      "  \n",
      "  return (sUm_squared - sUm^2/N)*(1/(N-1))\n",
      "\n",
      "end"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 5,
       "text": [
        "onepass (generic function with 1 method)"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "onepass(y)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 6,
       "text": [
        "0.9094409094409094"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "##3b."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "function twopass(y::Array)\n",
      "  N = length(y)\n",
      "  sUm = zero(y[1]) #oopse this killed the built in sum function\n",
      "  sUm_squared = zero(y[1])  \n",
      "  \n",
      "  for i in 1:N\n",
      "    sUm = sUm + y[i]\n",
      "  end\n",
      "  \n",
      "  m = sUm/N\n",
      "    \n",
      "  for i in 1:N\n",
      "    sUm_squared = sUm_squared + (y[i]-m)^2\n",
      "  end\n",
      "  return (sUm_squared)*(1/(N-1))\n",
      "\n",
      "end"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 7,
       "text": [
        "twopass (generic function with 1 method)"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "##3c."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "println(\"One pass: \",onepass(y),\", Two pass:\", twopass(y),\", Julia var(): \",var(y))\n",
      "println(\"One pass error: \",(var(y)-onepass(y))/var(y),\", Two pass error:\", (var(y)-twopass(y))/var(y))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "One pass: "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "0.9094409094409094, Two pass:1.001568444572786, Julia var(): 1.0015684445727868\n",
        "One pass error: 0.0919832644799167, Two pass error:8.867875425917323e-16\n"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The one pass has the larger errors under the given test conditions."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "##3e.\n",
      "\n",
      "Well, first off you might say never use the one pass. But what about under different $N$ or $\\mu$?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#Smaller N\n",
      "srand(42)\n",
      "N = 10000\n",
      "true_mean = 1e6\n",
      "y = true_mean+randn(N);\n",
      "println(\"One pass: \",onepass(y),\", Two pass:\", twopass(y),\", Julia var(): \",var(y))\n",
      "println(\"One pass error: \",(var(y)-onepass(y))/var(y),\", Two pass error:\", (var(y)-twopass(y))/var(y))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "One pass: 1"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        ".0147014701470147, Two pass:1.0199227833772877, Julia var(): 1.0199227833772826\n",
        "One pass error: 0.005119322085323428, Two pass error:-5.0072672132735025e-15\n"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#Smaller mean\n",
      "srand(42)\n",
      "N = 1000000\n",
      "true_mean = 1e3\n",
      "y = true_mean+randn(N);\n",
      "println(\"One pass: \",onepass(y),\", Two pass:\", twopass(y),\", Julia var(): \",var(y))\n",
      "println(\"One pass error: \",(var(y)-onepass(y))/var(y),\", Two pass error:\", (var(y)-twopass(y))/var(y))\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "One pass: 1"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        ".001568381207053, Two pass:1.0015684445727955, Julia var(): 1.0015684445728124\n",
        "One pass error: 6.326652926186049e-8, Two pass error:1.6848963309242483e-14\n"
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "If you knew for sure (or made sure) were dealing with a big array of small values, then you might prefer the faster one-pass algorithm. But also the two pass algorithm has a danger of catastrophic cancellation, so the one pass may be more reliable when $y \\sim \\mu$.   "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "##3d.\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "function var_online(y::Array)\n",
      "  n = length(y);\n",
      "  sum1 = zero(y[1]);\n",
      "  mean = zero(y[1]);\n",
      "  M2 = zero(y[1]);\n",
      "  for i in 1:n\n",
      "\t  diff_by_i = (y[i]-mean)/i;\n",
      "\t  mean = mean +diff_by_i;\n",
      "\t  M2 = M2 + (i-1)*diff_by_i*diff_by_i+(y[i]-mean)*(y[i]-mean); \n",
      "  end;  \n",
      "  variance = M2/(n-1);\n",
      "end"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 14,
       "text": [
        "var_online (generic function with 1 method)"
       ]
      }
     ],
     "prompt_number": 14
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "println(\"One pass: \",onepass(y),\", Two pass:\", twopass(y),\", online(): \",var_online(y))\n",
      "println(\"One pass error: \",(var(y)-onepass(y))/var(y),\", Two pass error:\", (var(y)-twopass(y))/var(y),\", Online error:\",(var(y)-var_online(y))/var(y)  )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "One pass: 0."
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "9094409094409094, Two pass:1.001568444572786, online(): 1.001568444614351\n",
        "One pass error: 0.0919832644799167, Two pass error:8.867875425917323e-16, Online error:-4.149899663066529e-11\n"
       ]
      }
     ],
     "prompt_number": 21
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Seems to work better than the regular onepass, but with all the subtractions there is a larger risk of catastrophic cancellation, so this may not be a good choice for narrow distributions. "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "##4a."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "function make_set(N,mu)\n",
      "    srand(42) #Random Seed\n",
      "    true_vel = zeros(N)\n",
      "    m = ones(N)*2\n",
      "    vel_err = randn(N).*m  #Did somebody say z-score? This takes me back to stats 101! :) \n",
      "    #0 velocity with no planet, right?\n",
      "    vel_obs = true_vel+vel_err\n",
      "  return vel_obs,true_vel,vel_err,m\n",
      "end\n",
      "\n",
      "y,z,err,mu = make_set(100,2)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 109,
       "text": [
        "([-1.1120537522876766,-0.8887667142145911,0.054310676018293196,-0.5989681807170408,3.55572219611071,-2.289803063449925,-0.9372117643337359,0.31228692527926727,-5.2839820161513655,2.0066198029163016  \u2026  0.2750796038556071,6.299151848253794,-1.4428243969728147,-1.154185393973309,0.9187675505910903,0.4944154121093404,-0.05132498760782233,-1.2641619360530398,-2.0408387097357954,-2.623015600406212],[0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0  \u2026  0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0],[-1.1120537522876766,-0.8887667142145911,0.054310676018293196,-0.5989681807170408,3.55572219611071,-2.289803063449925,-0.9372117643337359,0.31228692527926727,-5.2839820161513655,2.0066198029163016  \u2026  0.2750796038556071,6.299151848253794,-1.4428243969728147,-1.154185393973309,0.9187675505910903,0.4944154121093404,-0.05132498760782233,-1.2641619360530398,-2.0408387097357954,-2.623015600406212],[2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0  \u2026  2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0])"
       ]
      }
     ],
     "prompt_number": 109
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "##4b."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# I don't like one line functions as rule.\n",
      "Pnorm(y,sigma,z) = 1/sqrt(2*pi*sigma^2)*exp(-(y-z)^2/(2*sigma^2))\n",
      "\n",
      "#does it work?\n",
      "[Pnorm(y[i],2,z[i]) for i in 1:length(y)]\n",
      "#probably "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 108,
       "text": [
        "100-element Array{Any,1}:\n",
        " 0.170901  \n",
        " 0.180717  \n",
        " 0.199398  \n",
        " 0.190723  \n",
        " 0.0410699 \n",
        " 0.103572  \n",
        " 0.17873   \n",
        " 0.197054  \n",
        " 0.00608369\n",
        " 0.120585  \n",
        " 0.111041  \n",
        " 0.196013  \n",
        " 0.174414  \n",
        " \u22ee         \n",
        " 0.129381  \n",
        " 0.0275091 \n",
        " 0.197593  \n",
        " 0.001399  \n",
        " 0.153769  \n",
        " 0.168874  \n",
        " 0.179496  \n",
        " 0.193468  \n",
        " 0.199405  \n",
        " 0.163352  \n",
        " 0.118515  \n",
        " 0.0844064 "
       ]
      }
     ],
     "prompt_number": 108
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "##4c.\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "function Pnorm_set(ys::Array,sigmas::Array,zs::Array)\n",
      "  n = length(ys)\n",
      "  ps = [Pnorm(ys[i],sigmas[i],zs[i]) for i in 1:n]\n",
      "  return prod(ps)\n",
      "end"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 95,
       "text": [
        "Pnorm_set (generic function with 1 method)"
       ]
      }
     ],
     "prompt_number": 95
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "##4d."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "y,z,err,mu = make_set(100,2)\n",
      "println(Pnorm_set(y,err,z))\n",
      "y,z,err,mu = make_set(600,2)\n",
      "println(Pnorm_set(y,err,z))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "3."
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "776467160177524e-61\n",
        "0.0\n"
       ]
      }
     ],
     "prompt_number": 110
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "##4e."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "function Pnorm_set_log(ys::Array,sigmas::Array,zs::Array)\n",
      "  n = length(ys)\n",
      "  ps = [log(Pnorm(ys[i],sigmas[i],zs[i])) for i in 1:n] # I LOVE comprehensions! \n",
      "  return sum(ps)\n",
      "end"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 114,
       "text": [
        "Pnorm_set_log (generic function with 1 method)"
       ]
      }
     ],
     "prompt_number": 114
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Well, it can't handel big N at least.  "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "y,z,err,mu = make_set(100,2)\n",
      "println(Pnorm_set_log(y,err,z))\n",
      "y,z,err,mu = make_set(600,2)\n",
      "println(Pnorm_set_log(y,err,z))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "-139"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        ".1289017137387\n",
        "-895.3145308517971\n"
       ]
      }
     ],
     "prompt_number": 117
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Ah that is better! "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "##4f.\n",
      "When working with small-valued floats it can be useful to work in log-space."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "##5. \n",
      "Write it in Python! (just kidding)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}