{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Taking means of arrays with missing values\n",
    "\n",
    "Sometimes you do a thing that leaves you scratching your head for a few minutes until you realize \"Yes, of course that's how it works.\" \n",
    "\n",
    "Today at the (virtual) SciPy 2020 xarray tutorial, one of the exercises involved taking the mean of a 2D array that had some missing values. For some reason, my initial inclination was to take the mean of the 2D array over the first dimension and then taking the mean of the resulting 1D array. But that actually gave a different result from taking the mean over the entire array at once! \n",
    "\n",
    "The reason in the end is a very simple answer: because the data set has missing values that get skipped, calculating the means successively effectively weights measurements differently. \n",
    "\n",
    "But I thought it was a fun example of a seemingly innocuous mistake that could lead to large errors, particularly when analyzing large datasets. So let's break it down!\n",
    "\n",
    "To start, let's import `numpy` and build a random 4 by 4 array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np \n",
    "b=np.random.rand(4,4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To find the mean, we can just use `b.mean()`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5961189081171596"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b.mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "but it's also true that if take the means of just the column and then take a mean of the result, you'll get the same answer. When you take the mean along an axis, it will reduce the dimensionality by 1, so we'll get a 1D array here:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.67246434, 0.68927864, 0.7035238 , 0.31920885])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b.mean(axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And if we take the mean of that result, we'll get the same answer as before (with a tiiiiiiny bit of error in the final digit):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5961189081171595"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b.mean(axis=0).mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we get the same answer if we take the mean along the other axis first: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.5961189081171596\n"
     ]
    }
   ],
   "source": [
    "print(b.mean(axis=1).mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**But what happens when we have some missing values?** Let's put some in!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "b[2,3]=np.nan \n",
    "b[2,1]=np.nan \n",
    "b[3,3]=np.nan "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we have `nan` values, we can use `np.nanmean`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6060431635413454"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.nanmean(b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ok, and when we take the mean of the mean?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5516577241118834"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.nanmean(b,axis=0).mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6079093503218946"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.nanmean(b,axis=1).mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Whoa! Very different results! What's going on here?\n",
    "\n",
    "Let's think a minute about what a mean is. Given a vector, $B$, with $N$ elements, we add up the elements and divide by $N$:\n",
    "\n",
    "$\\frac{1}{N}\\sum_i^N(B_i)$\n",
    "\n",
    "But when we have `nan` values, $N$ changes when we take the mean over axes successively!\n",
    "\n",
    "Let's take a look at `b`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.80532163, 0.36377322, 0.72169393, 0.13140395],\n",
       "       [0.94655188, 0.91627271, 0.64348503, 0.10323544],\n",
       "       [0.39911748,        nan, 0.74945358,        nan],\n",
       "       [0.53886636, 0.85992325, 0.69946265,        nan]])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So when we take the mean over the rows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.67246434, 0.71332306, 0.7035238 , 0.1173197 ])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.nanmean(b,axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first column uses N=4, the second column uses N=3, which we can confirm:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6724643365974177"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b[:,0].sum()/4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.7133230623983705"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "col2=b[:,1]\n",
    "col2[~np.isnan(col2)].sum()/3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**This means the values in the second column have a larger weight than the values in the frist column!** Which is not what we want. There may be situations in which you *want* to do that, but most cases you just want the overall mean.\n",
    "\n",
    "What we *could* do is sum over the axes successively but divide by the number of total not-`nan` values in the original array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "13"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "NnotNan=(1-np.isnan(b)).sum()\n",
    "NnotNan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6060431635413454"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.nansum(b) / NnotNan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6060431635413455"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.nansum(b,axis=0).sum() / NnotNan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6060431635413455"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.nansum(b,axis=1).sum() / NnotNan"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we see these answers agree once again! \n",
    "\n",
    "But, just use "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6060431635413454"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.nanmean(b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**and remember**: *means of means won't give you the mean when you have missing values.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}