{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Circular statistics\n", "\n", "This notebook reproduces the examples of Berens (2009) implemented in the Matlab [CircStats](https://www.mathworks.com/matlabcentral/fileexchange/10676-circular-statistics-toolbox-directional-statistics) toolbox.\n", "\n", "## Part I | Artificial data" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import pandas as pd\n", "\n", "alpha_deg = np.array([13,15,21,26,28,30,35,36,41,60,92,103,\n", " 165,199,210, 250,301,320,343,359])\n", "\n", "beta_deg = np.array([1,13,41,56,67,71,81,85,99,110,119,131,\n", " 145,177,199,220,291,320,340,355])\n", "\n", "# Convert to radians\n", "alpha_rad = np.deg2rad(alpha_deg)\n", "beta_rad = np.deg2rad(beta_deg)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Descriptive statistics" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Alpha mean: 0.41, beta mean: 1.27\n", "Alpha std: 1.26, beta std: 1.44\n" ] } ], "source": [ "from scipy.stats import circmean, circstd\n", "print('Alpha mean: %.2f, beta mean: %.2f' % (circmean(alpha_rad), circmean(beta_rad)))\n", "print('Alpha std: %.2f, beta std: %.2f' % (circstd(alpha_rad), circstd(beta_rad)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plot" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAARQAAAEUCAYAAADqcMl5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztnXd4VNX6tu+dCgkljRZapNgIHIqUn6DSFNAAIgcE8VBExYBdEOsnnOPBQ1QUFRAEUVRQbCBBQCkqUgKYCAQQpBhII9SQIKnzfn9sZpiEJCRkZvaeybqva19DdtbseWbIfuZda73rXZqIoFAoFI7Ay2gBCoXCc1CGolAoHIYyFIVC4TCUoSgUCoehDEWhUDgMZSgKhcJhKENRKBQOQxmKQqFwGD5GC1B4BpqmvQ9cADKAe4DrgFxgK/C8iCQaKE/hIlSEUsXRNO0jTdOmVPIaGtAfWA50B2YDNwM9gQJgraZpIZVTqnAHlKEoSkXTtL2apkkpxxS7ph2BasCvItJHRBaKSKKI7Ab+BdQBuhrwFhQuRhmKoizuvvh4J9AACAf+BsYC04u1WykiBSVcoyb639kZJ+pUmARlKIqyqAcIsFFE0oFAIAA9Erlg124gsKyUa8wEfge2OFOowhyoQdkqhqZpLwAv2J3yB0TTtIl25/qJyEbgH8BhEcm+eL4teoRy0O56LYBmwJoSXmsG0A3oJiKFDn0jClOiDKXq8T6w1O7n6UAK8I7duZSLj22AnXbn2wKJImKxO3c3sE5Eztu/iKZpbwHDgB4icthB2hUmRxlKFUNETgOnrT9rmpYFnBaRgyU0bwOssvu5LUUNBvTuzsf2JzRNm4luJt1F5A9H6Fa4B2oMRVEimqZ5AZHALrvTzYEkuzZ1gC7ACrtzs4AxwHDgjKZp9S8eNVwiXGEoylAUpdEcfRDW3lB2A5M1Tet38ef+wHYROW7XZjz6zM46IM3usB+jUXgomioBqbhaNE1bDmwSkRijtSjMgYpQFJVhE7DEaBEK86AiFIVC4TBUhKJQKByGMhSFQuEwlKEoFAqHoQxFoVA4DGUoCoXCYShDUSgUDsPVa3nUHLUbs/qnrfTt3sVoGYqrR3P2C6gIRVFuWl17jdESFCZHGYqi3Pj7+xktQWFylKEoys38xd8ZLUFhcpShKIogAnFxMGoUhISAt7f+OHo09O48CrVSQ1EWylAUNvLzYexYGDYMIiMhMRFyc/XHVq1gwMAcxo7V2ykUJeHqxYHq+82kiOhmkpoKX38NgYGXt1m2ahtz3u5Ew4awYAFoTp8zUDgYNcujKJuyuijbtlHuLsq2bbBhQ+lmAnB3v058843ebvt2h70FhQehDMWNuVIX5d57KXcXZc4cGD++ZDMREQoKCvjvux8RECBER+vtFYriqC6Pm1KeLsr583DPPRTpouTm5pKenk5aWprtSE9PJyZmMt26RXPmzB7S09PJzMykoKCAgoICCgv1HTCCQutw9tQJNK0RIr8TGNiUwMBA6tevT4MGDWyP9of1XGBpYY/ClTi9y6MMxU2Ji9Mjk8TE0qOK1NRUNm/eycMP/x9t2kzjyJEvOHbsWClXLEDfoqfk7XO8vb2pU78R6SlJ6AnWFwDfcuutU6cO7dq1o3379rbHZs2a4eWlgmQXogxFUTKjR+vdmkmTdPM4cuQI8fHxJCQkEB8fT3x8PBkZGRdbTwJuAB7Ax8eHevXqXRZBvPHG87z33s/ceGMQDRo0ICgoCF9fX3x8fPDx8UHTNN6ct4SnHxrG0aOFtG/vTVLSec6dO1ck4ikp+klLSyM3N/ey91CrVi3atWtnM5j27dtz/fXX4+3t7cJPskqhDEVRMiEhwjvvbGDr1m+IjY0lKSnpsjZBQUG0b9+eli2789lnz7J9+xFatmxZ4g1rb1BXIiYG9u2DhQvLp1VEOHr0aBGzi4+PJy0t7bK2oaGh3HnnnURFRXHHHXcQFBRUvhdRlAdlKIpLpKamsnLlSmJjY/nuu2+w76KEhobSqVOnIl2KiIgINE0jPx+qV4eCkrYyv8iVulAAy9f8Qq+utxIZCUuXQqdOlXs/aWlpJCQk2Ixm+/btRbpkPj4+3HLLLdx1111ERUVx7bXXoqm56sqgDKUqIyIkJCSwfPlyYmNjiY+Pt/vtKSIjhzFoUBeioqK46aabSh2PSEmBNm3g1KmyXksf5E1JgW++KdlUvl/7OzNfb+u0PBQRYf/+/cTGxhIbG8uvv/5qGxAGaNGiBVFRUfTv35/u3bur8ZeK43w3FhFXHopycO7cOZk7d660a9dO0E1YAKlevboMGDBA5s2bJ0OGZEtMTPmuN326yOjRV26XlycyZoxIRIT+nORk/Vxysv5z06b67/PyKvX2ys3p06dlyZIlMmLECAkJCSnyWVxzzTUybdo0SU9Pd40Yz8Dp97gyFBMRHx8v48aNkxo1athunJCQEImOjpbvv/9e/v77b1vbrVv1Gz87u+xrZmXpRhAXVz4NFovedtQokZAQEW9v/XH0aJGxj3971e+tshQUFMivv/4qkydPloiICNvn4+PjI0OGDJG1a9dKYWGhYfrcBGUonk52drYsWLBAOnXqVOQb+JZbbpFPP/1ULly4UOLzLBY9WrjjjtJNJStL//2YMXr7ynL2XFblL+IACgoKZNWqVTJw4EDx8vKyfWYtW7aU119/XU6cOHHZcywW3YRHjhQJDhbx8tIfR43SDdQRn48boAzFUzl27Jg8/vjjUrt2bdsNUbt2bXn88cclMTGxXNdwdRflzyPHHHMhB3Ls2DGZMmWKNGrUyPY5+vn5yX333Sc7d+4UkaKfU0yMSEqKSH6+/hgTo593ZVfOQJSheBrHjx+XJ598Uvz9/W03QOfOnWXhwoVy/vz5Cl+vrC7Ktm2O1f7uwi8de0EHkp+fL8uXL5c777xTNE0TQDRNk2HDhsvgwWelT5/SI7nsbMdGciZGGYqncPr0aXnhhRckMDDQZiRDhw6VhIQEo6V5HH/99Zc89thj4ufnJ9BJ4LCMHBktSUlJpT4nO1uPVMo71uSmKENxd7KysuTVV1+VoKAgm5FERUW5pZF8/t1aoyVUiKSkJGnZ8lfRtGdtXaHHHntM0tLSSmxf3tkwN0YZirty4cIFeeutt6ROnTo2I+nRo4ds3rzZaGlXzaYdu4yWUGGCg0U2bjwkw4cPt3WFAgICZPLkyXLq1KkibZOT9e6iB6MMxR1Zs2ZNkanNLl26yLp164yWVSXx8tIHYEVEdu3aJQMHDrT9vwQFBcn8+fPFcnHgJC9PH4PyYJShuBOnTp2S0aNH2/5gW7duLStWrLD9wbo7k6fNMlpChQkO1mdz7ImLi5NevXrZ/p969eolhw4dUhGKMhTz8NVXX0m9evUEEH9/f/nf//4n+davRg/hQk6u0RIqzKhRUmJGscVikc8++0xCQ0Nt3aCoqF9k1CiPTo5ThmJ20tLSZPDgwUUS0v744w+jZTmF3X8cMlpChblSRnFGRoYMHz5cIFDgiLRqNabceUBuiDIUs2KxWGThwoUSHBwsgNSoUUNmzZrl0enf8xYvN1pChSlvRnG7dsclIGCJAOLr6yv//ve/JTfX/SKyK6AMxYykp6dLnz59bFFJ3759y8xxUBhLeTOKT5w4Kw8//HCRMbDdu3cbLd+RKEMxG1u3bpWGDRvaFu4tWrTIYwZdr8Sir743WsJVU5GM4vXr10vz5s0FkMDAQFm6dKkhmp2AMhQzMX/+/IvZl0i3bt0kNTXVaEkuZceufUZLcBnnz5+XESNG2KKV5557TgoKCoyWVVmUoZiB3NxciY6Otv1xTZgwwRP714piWCwWeeutt8Tb21sA6dOnz2XJcG6GMhSjSUtLk65du9qmgz/88EOjJRnGk1PeNlqCIaxfv17CwsIEkGbNmtlWMbshylCMZMuWLRIeHi6ANGzYULY5evmuwm1ISkqS9u3b23JWPv/8c6MlXQ3KUIxi/vz54uvrK4DceuutqtSgVK0xlJL4+++/5V//+pet6ztp0iR3SxNQhuJqLBaL/Pe//7X90Tz22GOSVwUq75QHd57lcRQWi0VmzpxpG1cZMWKEO/19KENxJRaLRZ599llbcZ7333/faEkKk7J27Vpb7d+BAweWWqrTZChDcRUFBQUybtw4W+HjJUuWGC3JdLhjpqwz2bp1qy1TulevXpKVZY6au2WgDMUVFBQU2PrG/v7+Ehsba7QkU+KOa3mcza5du2yLQm+++WbJzDxn5mLYTr/Hq/xOSYWFhTzwwAN88sknBAYGsmrVKu666y6jZZmSFhGNjJZgOlq3bs3GjRtp3Lgxmzdvo2XLX7j3XguRkfoujLm5+mOrVnDvvfpmavn5Rqt2Iq5wLbvDVBQWFtrqlwQEBMjPP/9stCRT4471UFzFwYOHJCDgc4FVcvPNt5fY/TFBMWzV5XEWFotFxo4dazOTn376yWhJCjdm61aRRo3ypEGDFrZUg5J2MTC4GLbq8jiLqVOnsmDBAqpXr05sbCy33Xab0ZJMz+bfdhstwbTMmQOPP+7LL7+sIjw8nF9++YWRI0disViKtAsMhOhovb1H4grXsjtMweLFiwUQLy8vNQBbAdyt6r0rsS81uXfvXtsGbi+88MJlbQ0sNam6PI5my5Yttk223n67aq5NUTge+2LYInqhcmvy26JFi4q0NbAYtjIUR5KUlGSb4hs3blyVqWPiKMy8c6DRlFQMe9asWRcrwPnJ++8n2KaSNU0/DJhKdvo9romIS3tYrnwxe7KysujWrRu7du2iV69erFq1Cl9fX6PkuCUH/0pWU8elMHq0PjU8aVLR8+PHP8GcOf/Ay6sX998fxN9/1yY2FnJyoHp18PaG7t3h66/Bz8/pMjVnv0CVGJQtLCzkvvvuY9euXVx77bV8+eWXykyugjqhQUZLMC3R0TB7Npw/f+mcCPz991uEhv4Di+VnlizJ4rrrLhAWBps2wcGD8NxzsHYtREZCXp5x+h2GK8Igu8MQnnnmGQEkODhYDhw4YJQMt+f/vfGB0RJMS0nFsK0V9++9N1cCAzcK1JHQ0B0yalRhkW7O8eMi1aqJREU5vfvj9Hvc4yOUTz/9lDfffBMfHx+++eYbWrZsabQkt2XqMw8aLcG0aBrMnQsNG+rRRkwMvPEG9OsHmzb58fDDN+LltYNTp3YSGvoCml3no25dePFF+Okn2L7dsLfgGFzhWnaHS0lKSpJatWoJoFYOO4D1m34zWoLpsS+GbR18DQjQi2HPm5cg3t7eomnaZYmUycmX2jkRp9/jHmsohYWFtu0mBw4cqGZ0HMCy1WppQkXw8hIJCio6+/Pyyy8LIBEREXLu3DnbeetUspPzU5ShXC3vvvuuABIWFibHjx935UsrFCJyaYrYPj8lNzdX2rVrJ4A89NBDtvPWZDcn56c4/R73yDGUAwcO8OyzzwIwd+5c6tata7Aiz+DNeUuMluBWDBgA1apBRsalc35+fixatAg/Pz8++OADZszYxKhR0LIlnD4NFos+Bb1tmz5L5G54nKEUFBQwcuRILly4wP33388999xjtCSPYWhUT6MluBXR0XqeycKFRc9HRkYydeo0YAHPPtuYhg3PExoK48fDkCFuXurAFWGQ3eF0rPVgGzZsKGfOnHHFS1YZjp88bbQEt8Ji0aeCq1XTp4btz48eXShBQVsF6kj9+r/L/feLNGlyaRWyk0odqC5PRfj999+ZMmUKAAsXLiQoSCViOZL5i78zWoJboWl6BmzjxtC0Kbz6KqSkwObNsHatF2PGtETTtpOevpMtW07i6wt9+uhRTePGEBYGq1fr3R+3wRWuZXc4jfz8fGnTpo0AMn78eGe+lEJRIXJz9UglMFCfGtY0kerVRYYOFbnttt0C5wSy5MUXsyQlRR/ETUkRiYnRB3ZbtNBngRyAilDKy4cffsiuXbuIiIggJibGaDkeyeqfthotwS3x84PvvoP16/UxEtDT7NeuhZSUVtSufQyoR2bm84SHg48PhIfr64Li4iApCcaNc49BWo8wlOzsbF555RUA/ve//xEYGGiwIoWiKJoGnTrBRx/p/87Jge+/h4ICjR9+EDQthzlz4hk0KJOQEL3bExIC//kPFBbChg3ukUXrEYYyY8YM0tPT6dixI0OHDjVajsfSt3sXoyV4BLVr61PJc+boMzvt2rWiRYufKSz8lIMHlxcpbt2kif6csDCYNctY3eXB7Q3l+PHjti7O66+/jqY5fYV2lWXaux8bLcEjGDAAPvtM7wbdd5/enQkPv4lq1TqSmDiKpKQttm5PrVp6m4AA+OILN+j2uGKgxu5wOI888ogA0r9/f2dcXmGHmjZ2DNZVyJom8uuv+r+zs0VefPFFAaRr165isVgkK0ukaVN9KvnMGRGodHFrp9/jbm0o+/btE29vb/Hy8pI9e/Y4+vKKYhxNURvGOwJrqQMfH5FBg/TZHBGRzMxMCQsLE0AWL/6uSB6KgxYPqlmesnj++ecpLCxk7Nix3HjjjUbL8XiWxq43WoJHYC11EBEBy5ZBZqaen1K9ei2eeCIGmMTIke1o0MDC3Ll6+88+g6govZtkalzhWnaHw9i4caNtT53U1FRHXlqhcAlbtuh9hOHDLy0MDAmxSM2aXwncJHPmzBERsXV9Nm2q9OJBFaGUhjUjduLEiTRo0MBYMVWE5Wt+MVqCR9G5s56jkpoKR49CQQGcOqVdXPuzg//85z+cOZPP4MHQs6c+41O7ttGqy8Yti1Tv3buXVq1aERAQQHJyMsHBwY64rOIKbNgcT4+b2xstw6MYORL27YOTJ/XFhCNGQJ06wo039uDQocH4+0cj4kN+vr5yOTwcFi+Gjh3hKiY0VZHqknjvvfcA+Ne//qXMxIUoM3E8EybAiRN6wtvevdCmDQQEaCQlrQVGExKyjCNH9NIGoaF6SUkzr0R2uwglMzOThg0bcv78eXbv3k1kZKQjdCnKwStvzld1ZR2MiG4OKSnwzTd6vsnYsXDsWAFxcY3Jykpn06adTJ3ahoYNYcEC+PtvuOcebD9XIFJxfpKWKwZq7I5K8/bbbwsgPXr0cMTlFBXg7LksoyV4JHl5+vRwRITIhAkijRrpeSdjx/4/gUlSo8YJGTOm6ALBq9x03en3uFtFKBaLheuuu46DBw/yzTffMGjQIEfpUpQDtdGX8xDR1+rcdx+kpemp9zVrFnL27CL8/BaQmrqc0NDQIs+JidHHX4oXcCoDNYZiz5o1azh48CBNmjShf//+RsupcqjVxs7Dunjw9Gn48099xufMGW/69fuSvLxNzJ8//7LnjBhhvrwUtzKUd999F4Do6Gh8fHwMVlP1eHT0P42W4PFkZur79Fh59NFHAZg9ezYFBQVF2tatq7c3E25jKAcPHmTVqlX4+/vz4INqYNAIvlixzmgJHo91JbKVvn370qJFC44ePUpsbGyRthkZ5stLcRtD+eijjwAYPnw4YWFhxoqpojQOV7sHOBvrSmQrXl5eTJgwAYAFCxYUafvZZ3p7M+E2g7KRkZHs2bOHH3/8kd69eztSk0JhGuLiYNgwvRaKtU5Yeno64eHh+Pn5cfLkSWrUqEF2tr7l6dKl+thLOamag7Ii+gc7ahQXq1cJe/b8gq/vZwQEdMe1Hqiw8txrs42W4PF06gQ9euh5JufP6+fq1atPq1YPkJs7lwYN/PDyguBg8PXV7xUz3Q+mM5T8fD2xZ9gw3YETE2HatBlAJK1awYgRPqbNEvR0pjytxq6cTfFN1197Tb8XkpPfAI4SEbGKJk1g8OBL94mp7gdXJLvYHWVirRPRp4+euGOla9euAsjSpUudtV+Johzs/uOQ0RKqDBaLXoipRQsRX18Rb2+LwEnx8/tMNm++tLdpBe+HqrXaeNs2vRjv119f6j8eP36czZs34+fnR9++fQkM1FOU3aVoryexJT7RaAlVBms6vZ6Pohezvv76buTljSAnZ6OtndnuB1MZirVor33R+hUrViAi9O7dm5o1awL676Oj9fYK1/HQcJNNKXg49veDCHTq9BjwEf36dbJVxR89GvbsgUceMcf9YCpD+e47PfvPnmXLlgFw9913FzlvxixBT+eTr1cZLaFKYb0frOOKP/74AJBIaGh3cnKExMRL+yDHx5vjfjDVtLG3t76GwZoEm5WVRZ06dcjLyyM1NZX69evb2ubnQ/XqekiocA2/7f6DDq2vN1pGlcHbW9+/Z9w4vQjTl19auO66RqSlpREfH0+7du0AfTbo7rth3Tp9D58yVh9XrWnj4lmC69atIzc3ly5duhQxEzBnlqCno8zEtdSuDWvWXBpXrFnTi4EDBwL6UICVwECYPVs3EqPHUUxlKMWzBOPi4gDo2bPnZW3NmCXo6Tw1dabREqoUAwbAlClFxxV79OgBwLZiO6h/+y20a2f8OIqpujzFswR79+7NunXr+Pbbb4uMoVxllqDCk1m9GgYN0vvBAQH6H1DNmvrXfFCQPoIZEqJnhNWsqRdpveEGo1WXSVwc3Hwz7N8PLVro5w4fPkzz5s2pW7cu6enpaJpmux/ee09PBj11qtRLVq0CS9Y8lDvuEMnKskjt2rUFkOTkZFubrCyVh2IUO3btM1pC6ezfL+LnZ00cLf3w9tZ32Jo40WjFV8Ri0SXffvulvCyLxSIhISECSFJSUpH7ITf3ilXxq1Yein2W4A03FJCZ+TB16rSlTp1wUlL0gjKRkfrvrfuVKFzH3gNHjJZQOtdeC926Xbmdlxe0bQuvvup8TZVE0y4FV5GR+t9/aqpG+/adgXBeeimzyP1w4oQJxhVd4Vp2R7md+d//Xi2wUHx9My/uV6LvmrZtW3mvoqhy/PSTSGBg6dGJpomEh4ucPGm00nIzapTI9Ol6qcdRo/T7QNMKBE5K69bbi9wP06dfcWdBp9/jpjQUEZGnnnpKAJk6dWpFnqZwIvMWLzdaQun8/bfIBx/oMX9phlKrlsiBA0YrrRDWfZDtl6J8++23AkivXr1s5+z3QS6DqtXlsWfHjh0AdOzY0WAlCiv/196EOwz89Rc89RSEhcFDD+mJGDVqXN6uenWIjYWWLV0usTKUtPrYek/s2LEDESE7G9tmYEbfLqY0lMLCQuLj4wG46aabDFajsGKaAtUisHYt9O4N11wDb7+t7y1xww36Bjd+fkXbBwTAvHlwyy2GyK0MxVcfx8SASDj16jUiMzOQyZNPmWtc0RVhkN1RLhITEwWQpk2blvcpChcwedosYwWcOyfyzjv6PhPWboyXl8g//1k01v/vf0WqV9d/HxAg8vLLxml2EBZLyeMot956qCLjilVzDOXLL78UQPr371/uT0rhwfzxh8hDD4n4+18ykpAQkalTRY4fv7z96dO6oVSrJjJkiEfmF7z00ksCyEsvvVSRpzn9Hjdl6fjU1FQAGjdubLAShT2bf9vNzR1au+bFCgvh++9h2jTYard9R4cO8PzzMHDgpUVfxQkOhkcf1Z/36acm6Ac4nkaN9O5nWlqawUqKYmpDCQ8PN1iJwp5jqRnQwckvcvo0fPABzJhxaWGXry8MHw4TJ0LrchpaTIzzNJoA671hvVfMgjIURbm5t38v5118507dBJYuvbSEvH59eOYZePBBPcNLYcOshmLKWR7rh9SwYUODlSjsee+jrxx7wfx8+OILPXO1bVtYvFg3k27dYOVKfQfxiROVmZSAWQ1FRSiKctO3exfHXCg9XV8W+847cPasfq5aNRgzRs8pcbNcESOoW7cuXl5enDhxgry8PPyKT5UbhDIURbmpE1qJSEFEHySdPh1WrACLRT/ftCk8+yyMHFlyQpqiRLy9valfvz6pqamkp6fTpEkToyUBJuzynD9/nszMTPz9/QkODjZajsKOGfM+r/iTcnJg4UK4/np9Lf7y5bqZ3HGHXjnoyBG94Icykwpjxm6P6SIU6zRYeHg4mgdO97kzU5+pwL48R4/CzJl6+qY1ZzwwUK+m/PjjYJJvVHfGjIZiugjF+uE0aNDAYCWK4mzYHF92AxFYv16PPpo21ad+z5+H667TU+JPnoQ33lBm4iCshpKSkuK019A0bbymaUc0TcvRNO03TdPKXL9gugglKysLgCA1sm86zmVll/yL7Gz4+GN92vfoUf2cl5eefDZ5sr7CTUWbDqf2xeIn2dml/L9UEk3T7gVmAuOBXy8+rtI07UYROVrik8qTTgvsRS/fWNIxpQKpuVdk2bJlAsiAAQMqklKsMIL9+0UeeaRoSnxQkMgrr4ikpxutzmOx7irYunW8wCnRtEIJDtbX+cTFlbnSoGLrciAO+KDYuT+B10p7Tnm7PNaCrncCDYBw4G9gLDC9AqZ3RQouJjX5lJZWrTCMN+ct0QdUV67Uc0Wuuw7ef1/f+6RdOz0pLSNDr6xcr57Rcj0S+72/69c/CUTy4ov/KbJHjyP2OtY0zQ89L/qHYr/6Abi5tOeV11DqoUcjG0UkHQgEAoBfReRCxeWWTv7FT8LX19eRl1VUlnPnGHomBcLDISoKNm3S19Lcfz/8/ru+09SQIXqavMIpiFzaoycxEW67bRuQRmFhLuHhMGmSfj4lRW8nlas/HwZ4A8eLnT8O1L+8uU55w4B/AIdFxNpZa4seoRysoMgrYo1QvL29HX1pRWV47z38Z7wNp05A3bqXUuJDQoxWVmWw7v1t3RXC+qWbbxeOWPc6jozU9+hxwK4QxW1JK+GcjfJGKG2AnXY/twUSRcRSEWWrf9rK6p/0laPT3v2YjFNnOJZ6XA+lgeVrfuH4WT3gCW3WlsysbA7+lWxL+f5ixTo2/7YbgOdem01Obh6J+w/zwRJ9D8ZPvl7Fb7v/AC7tIfPb7j9sW2h+sOQ7EvcfJic3j+demw3oK2i/WLEO0FPLD/6VTGZWNq+8OR/QZzaWr/kF0EP+Y6nHyTh1hmnvflzu92SdHXnlzfnu+55u78vrE56HFStY/vFSNnTrDSEh7v2e3Oz/6Y03LzBy1N/M/FB/T+fyvOjSvS9eXl5F3tP7ny0hOhqefym9yHuqICeBQi6PRupyedRyiXIOzmwFXrH7ORaYV9FBnvIMOH3xxRcCyJAhQ8rTXKGoMgQHi6SkXPr5P//WztdnAAAYE0lEQVT5jwDy4osvXtY2OVkvGVOMqxmUnVfs3AEqMyiraZoXEAnssjvdHEgqh8tVGOtgbIHatNh0WL/hFcaQman3Nq1YuzolTWDUrau3ryQzgNGapj2oadoNmqbNRJ+Qeb+0J5RnDKU5+iCsvaHsBiZrmhYvIqsqo7g4JfULFQrFpb2/rUvcyprAcMTe3yLyhaZpocBL6LO7icCdIlJqMHHFCEVE/hQRTUQO2Z0bKiK1HG0mcOnDycvLc/SlFZXEYauNFVdF8b2/y4pQHLX3t4jMFpEIEfEXkQ4i8ktZ7U2Xeh8aGgrAiRMnDFaiKI51gFNhDNHRMHv2paVR1nskLCysSLvsbL1ddLSrFZrQUMy44Emh8+B9DvjKU1w1xffoKanMh9F79JjOUOrVq4emaWRkZKhxFJORm6u6oUZSfI+eXbv6AOGm2vvbdIbi4+NDvXr1EBGOHy99ulvhepbGrjdaQpXH1xcWLNArZ5450wDYRZcubWjTBvbtgy+/hA8/NC5hWZNK5udWkHK9WIcOHYiPjycuLo5ODkj1Uyg8jQsXLhAQEICvry85OTl4eZUrNnB6zGK6CAUuFadW4yjmwpq1qTAe+/GTcpqJSzCPEjvUwKw5qVVTlWk0C2atu6wMRVFuetzc3mgJiosoQ6kArihtp6g4V7HATOEklKFUAOu+rYcPHzZYicKepx8eZrQExUWOHDkCXLpXzIIpDaVt27YAJCQkYLFUqEKCwomcOHXWaAmKi+zYsQO4dK+YBVMaSv369WnUqBFZWVkcOHDAaDmKi6jVxuagoKCAhIQEAG666SaD1RTFlIYClz6o7du3G6xEYeXR0f80WoIC2LNnDzk5OTRv3pwQk1XMM62hdLy4EMEa2imMx1pdTGEs1nvCbNEJmNhQVIRiPhqH171yI4XTsd4THY1Y/XcFTG8oCQkJqnqbSbi5Q2ujJShQEcpVERISQrNmzcjJyWHPnj1Gy1GArbiywjhyc3PZtWsXmqbRvr35Eg1NayhwKaRT3R5zMOXpCmyWrnAKO3fuJD8/nxtuuIGaNWsaLecyTG0o1pBu27ZtBitRABz8K9loCVUOEYiLg1Gj9C2QunS5CThFbu77bNtW6c28HI6pDeXWW28FYNWqVbi4zIKiBLbEJxotoUphv+1oZKS+wVefPgOBSNq183fYtqOOxJT1UKxYLBYaN25MamoqO3bsoEOHDs7SpVCYChHdLFJT4euv9R0Bs7KyCAsLIz8/n/T0dAID63LPPXqFtgULylWhrWrWQ7GiaV506fIE8BHdut2At7ce9o0ejSnDPU/HuguewvlYtx21mgnA6tWrycvLo2vXrtStW9e27eiGDfq2o2bAtIZiDfc2bXoMSKRp0yhyc3H4LvOK8nPjtdcYLaHKMGcOjB9/yUwAli1bBsDdd99tOxcYqFe3nzPH1QpLxpRdHvtwb8mSPCIi6nDu3Dn+/PNPWrRoAehVvysY7ikUbkNIiP7laa1OkJeXR926dcnMzCxyHwCkpECbNnDq1BUvWzW7PPbhXnCwH3fddRcAy5cvt7UxY7jn6Vg3AVc4n+Lbjv78889kZmbSqlWrImYCDtt21CGY0lCKh3vWEM8a8lkxW7jn6bz1yhNGS6gyWLcdtVJSd8eKI7YddRSm7PIUD/fOnTtHnTp1KCgoIC0tjbp21l2BcE9RSX7b/QcdWl9vtIwqwejR+ljhpEkgIjRu3JiUlBS2b99+Wcp9TIy+hcbChVe8bNXs8hQP92rVqkWvXr2wWCysWLGiSFszhXuezt4DR4yWUGWw33Z0x44dpKSk0LBhw8tSJ4zcdrQkTGkoxcM9gMGDBwMwf37RuqZmCvc8nX8N7me0hCqD/bajc+YsAvR7QLObfTB629GSMKWhFN9lHmDYsGEEBQWxdevWIjVSHLXLvOLKfLDkO6MlVBms246GheXy0UfPAJMYNOhR8vMxzbajJWFKQym+yzxAYGAgDzzwAADvvfceYL5wz9P5v/aRRkuoUvj6Qvv2sxAZSoMGvRk8uCXVq2OabUdLwpSDstY8lJQUfWrYOttz+PBhWrRogZ+fH/v2HeORR+qoPBQXkpObRzV/P6NlVBkKCwu59tprOXz4MMuXL2dA5UPxqjkoW3yX+ZgY3VwaN25Gz57/Ijf3cdq39zVduOfpTJmh9uVxJatWreLw4cNERETYcrHMjikjFFtj0ZPWZs+GFSv02ZyAgDyysj6jbt2vSUlZho+Pj7O0KhSG0qdPH3744Qdef/11Jk6c6IhLOv2r19SGUhIWi4UbbriBAwcO8OWXX/LPf6pK7K5i82+7VRlIF7F//36uv/56qlevTnJysqOq21fNLk9ZeHl58eijjwKXBmcVruFYasaVGykcwqxZswC4//77TbdVRlm4XYQCeuZsw4YNyc7OJiEhwXS7pykUleHs2bM0adKErKwsdu7cSZs2bRx1aRWhlEStWrV48EG9vumLL75osJqqw3sffWW0hCrB9OnTycrKomfPno40E5fglhEKQEZGBs2btyA7+wZuv30ZO3Y0IDNTz5odMEBfXNixo5oBciQH/0qmRYS5Nuf2NI4dO8a1115LTk4OcXFxdOrUyZGXVxFKaQQH1+XaazcCn7N37xfs2mVRBZicTJ3QIKMleDwvv/wyOTk53HvvvY42E5fgloYiAuPGQUhIJPXr305KylP88svn+PjoK5QnTdKNJSVFb6dKRTqGGfM+N1qCR7Nz504WLVqEr68v06ZNM1rO1SEirjwcwtatIhERItnZIgsWLBBAmjZtKjk5OUXaZWfr7eLiHPXKCoXz6NOnjwDy5JNPOuslnH6Pu2WEYl+AadSoUbRq1YqkpCTbVJsVVYDJsWzYHG+0BI/lxx9/ZM2aNdSuXZuXXnrJaDlXjytcy+5wCMHBIikpl35euXKlABIcHCynT58u0jY5WSQkxFGvXLVZtvpnoyV4JIWFhdK2bVsB5H//+58zX8rp97hbzvJ4e0NuLliz7kWE3r17s379eiZOnMjrr79ua5ufD9Wrg9pvXWFWPvnkE0aOHEmjRo04cOAA1atXd9ZLqVmekihegEnTNGJiYgCYOXMmu3btsv1OFWByHG/OW2K0BI/j1KlTPPvsswC8+uqrzjQTl+CWhlJSAaYOHTrwyCOPkJ+fz8iRI8nLywNUASZHMjSqp9ESPI4JEyaQnp5Ot27duP/++42WU3lc0a+yOxyC/SyPPVlZWdKsWTMB5IUXXpCsLJGmTdUsj6M4fvL0lRspys2SJUsEkMDAQDl06JArXlLN8pSEfb1N+6puNWrUYNGiRWiaxmuvvUuvXmdNVW/T3Zm/WJWAdBSpqamMHz8egBkzZtCsWTODFTkGtxyUBX2wddw4faOv6GgYMUKvgJ+RASNGfM/PP99IrVo7SErqR1BQ4JUvqFC4CBHhrrvuYtWqVfTr14+VK1cWKT7tRNSgbGn4+uqlH7/4Avbu1etsWuttNmnSh2bNnuPcuSG8/PJzRkv1GFb/tNVoCR7B/PnzWbVqFcHBwcyfP99VZuIS3DZCuRIJCQl06tSJgoICfvzxR3r37u2ql/ZYVv+0lb7duxgtw605fPgw//jHP8jOzmbx4sUMHz7clS+vIpSrpV27dkyZMgWAMWPGcPbsWWMFeQDKTMqHCMTFwahR+i6Y3t7646hRFgYPnk52djZDhw5l2LBhRkt1OB5rKACTJ0+mc+fOJCcnM3LkSAoLC42W5NZMe/djoyWYnvx8fZX7sGF6gfXERGyr4I8dW8Pvvz9H9eqLmTlztkd1dWy4YirJ7nA5f/75pwQHBwsgEydONEKCx6CmjcvGYhEZM0akT5/LUxo+/vhjAcTLq6Z06HBSxozR27sYp9/jHm8oIiIbNmwQHx8fAWT+/PlGyXB7jqakGy3B1JSWH7Vx40bx8/MTQGbPnm3kKniVh+IIunfvzpyLS44feeQRfvrpJ2MFuSlLY9cbLcHU2K+Ct3L48GEGDRpEXl4ejz32GNHR0Z69Ct4VrmV3GMrTTz8tgISEhMiff/5ptByFh1F8FfzZs2flxhtvFED69u0r+fn5tt8ZtApedXkcSUFBgURFRQkg11133WWlDhRlo8oXlI2Xl4jVM/Lz86Vv374CyI033ihnz54t0jYvT8Tb2+USVZfHkXh7e7N48WLatGnD/v37GTJkCPmq6Gy5qVWzhtESTI39Kvinn36a1atXExYWRmxsLLWLLXn31FXwHpvYVhZJSUl06tSJjIwMHnzwQebNm+eZU3gKlzJ6tF4g3d//HZ544gn8/PxYv349Xbt2vaxtTAzs2wcLF7pUovP/yF0RBtkdpmHLli1SrVo1AWT8+PFiMWAOz934f298YLQEU7N1q0hoaKZAgADy8ccfl9jOwFXwagzFmaxevVr8/f0FkEcffVSZyhU4ey7LaAmmZs6c9wUWCKyW11+fXWKbrCyRO+4QlYfioMN0fP/997Ycgccee0yZShn8eeSY0RJMy7x58wQQ8JFOnRIlIkJk+nR9NicvT3+cPl2PTMaM0c8ZgDIUV7By5UqbqTzyyCNSWFhotCRT8u7CL42WYEreeeedi2aCzJgxQywWvTszapQ+NeztrT+OHi2ybZuhUpWhuIqVK1fauj/33Xef5Bn0FaJwHywWi7z66qs2M3nrrbeMlnQllKG4kvXr10uNGjUEkAEDBsiFCxeMlmQqPv9urdESTIPFYpFJkyYJIJqmyQcfuMWAtTIUVxMXF2dbTNi9e3c5ceKE0ZJMw6Ydu4yWYApycnJkzJgxAoiPj498/vnnRksqL06/x6tkHsqVSExM5Pbbbyc9PZ2mTZvy7bff0q5dO6NlKUxASkoKgwcPJi4ujmrVqvHVV19x1113GS2rvKgCS0YQGRnJ9u3b6dSpE0lJSXTt2pXFixcbLctwnnttttESDOXXX3+lQ4cOxMXF0aRJEzZt2uROZuIaXBEG2R1uxYULF+SBBx6wDbo99dRTRRZ4VTUu5OQaLcEQLBaLzJo1y1YCo0ePHpKRkWG0rKtBjaEYjcVikdmzZ3vCH1Ol2f2HS/aOMRUe9qWiDMUsbNy4UerVqyeANGnSRH777TejJbmceYuXGy3BpRw7dkw6duwogFSvXl0+/fRToyVVFmUoZiI5OVk6d+4sgFSrVk3mzp2rMms9lNWrV0vdunUFkKZNm0p8fLzRkhyBMhSzkZOTIw8++KAtBO7Zs6ertpE0nEVffW+0hEphsegL+EaO1IsheXnpj6NG6ZmtFovIqVOnZOTIkbb/3169enlS6oAyFLOyZMkSCQsLs4XDM2bMkIKCAqNlOZUdu/YZLeGqycvT19BERIjExOiV1fLz9ceYGP18jx6HpU6dcFsEOn36dHceLykJZShm5sSJE3LffffZvs06d+4siYmJRstSFKOsavQiIqmpqTJgwHCB1QILpFu3W2T//v2uF+p8lKG4AytWrJCGDRsKIL6+vjJ16lTJzfW8KdYnp7xttISrorRq9BaLRT788EMJCgoSQAID60poaKZs2eKxi0OVobgLZ8+elXHjxtmilcjISNm6davRshSij5HExBQ9d+jQIbn99ttt/1/9+vWTpKQkmT5dXxXsoShDcTc2bNggzZs3t/2hDho0yGO6Qe46hmJfjT4lJUWio6PF19dXAAkNDZVPPvnENltnUDV6V6EMxR05f/68PP/887YSk5qmyf333+/2W3e46yyPl5dIWtoJeeaZZ4r8n4wcOVKOHz9epK1B1ehdhTIUdyYlJUUmTJhg+zb09vaWhx56SI4ePWq0tCrD2bNnpVq18xIQ0MIWNQ4ePFj27NlTYnsVoShDMT1HjhyR0aNHi5eXlwDi7+8vTz755GXfjiLly5UwCnfKlM3OzpbXXnvtYimKhQITpV+/flfMcFZjKMpQ3IZ9+/bJ0KFDbd+UgYGB8sQTT8jevXtFpHy5EgbWI3WLtTwpKSny73//27ZMApC2bcdJgwYXSpwytsfAavSuQhmKJ5KQkGDbwdB63HLLrXLbbQfl9tsLS/3Dz842tGK6aVcbFxYWyurVq2XQoEHi7e1t+0w7duwoa9askcJCi4wZo392pX22BlejdxXKUDyZ+Ph4efjhhyUwMFCgk8BhCQ1tIhMnTpQDBw6U+JzsbD1SMeJbdPK0Wa5/0TJIT0+X1157Ta655hqbifj4+MjgwYPlhx9+KLLOyj76M2E1elehDKUqkJmZKV26/CENGswoErX06tVLli5delmSnIf388uksLBQ1q1bJ0OGDLENdlsX8P33v/+VtLS0Up9r4mr0rsLp97gqAWkSQkJg924hOXkbc+fO5fPPP+fChQsA1KpViz59+hAVFUW/fv3Iy6tDmzZw6pRrNW7+bTc3d2jt2hcFsrOzWbt2LbGxsaxcuZL09HQAvLy86N+/P+PGjeOOO+7A29vb5drcDLUVaVXBy0sfgLVy5swZeeedd6RNmzZFohZN06Rz527i5VUov//+u0vLJ7iy6v2hQ4fknXfekTvuuMO2Z5L1aNKkiUyZMkWOHVMbj1UQFaFUFUJCIDERwsMv/92RI0dYuXIlsbGxbNiwgby8MGAXEEajRo1skUvnzp2pV6+eq6U7hLNnzxIfH8+aNWuIjY1l7969tt9pmkaXLl2IiooiKiqK1q1bq83trw6nf2jKUEzC6NHQqhVMmlR2u+zsbCZM+IstWzLJyvqnLfy3Eh4eTrt27Wjfvr3taNy4sUNuwPc++opHR/+z0tfJyMggISGB+Ph44uPjSUhI4NChQ0Xa1KpVi759+3LXXXfRr18/6tSpU+nXVShDqTLExcGwYXqUEhhYervsbIiMhKVL4aabLCQkJBAbG8v69etJSEggKyvrsueEhobaTKZVq1Y0aNDAdoSEhJTbbA7+lUyLiEZXbCcinDt3jvT0dNLS0khLS2P//v02E0lOTr7sOf7+/rRp04Zu3brRv39/unXrhq+vb7l0KcqNMpSqggiMHQspKfDNNyWbSnY2DB4MDRvCggVQ3AcsFguHDh0q8u0fHx/PqTJGb/38/Khfvz7169cvYjTBwcH4+PgUOfILBS8sFBQU2I7ixmE9rAPKJVGjRg3atm1bJIq6/vrrlYE4H2UoVYn8fBg3DjZsgOhoGDEC6taFjAz47DOYPRt69oS5c6G8956IkJycTHx8PL/99huHDh0qcuNnZmaWW9+Dz0xl/puvlKttYGAgDRo0sBlVRESEzTxatGiBl5faEsoAlKFUNURg+3bdPFasgMxMqF0bBgyA8eOhY0fHvt6FCxdKjDDOnTtXJBIpKCjAYrHg6+tri1h8fX0JCAiwRTX2UU7NmjUdK1ThCJShKMzDhs3x9Li5vdEyFFeP2opUYR7OZWUbLUFhclSEolBUHVSEojAPb85bYrQEhclREYqi3BxLPU7jcPfMxFUAKkJRmAl/fz+jJShMjjIURbmZv/g7oyUoTI7q8igUVQfV5VGYh9U/bTVagsLkuDpCUSgUHoyKUBQKhcNQhqJQKByGMhSFQuEwlKEoFAqHoQxFoVA4DGUoCoXCYShDUSgUDkMZikKhcBjKUBQKhcPwMVqAwjPQNO194AKQAdwDXAfkAluB50Uk0UB5ChehIhRFpdH0jX36A8uB7sBs4GagJ1AArNU0LcQwgQqXodbyKEpF07S9wA2l/HqqiEy52K4TsAqoJyIFxa5RA8gE7haRFU6UqzABKkJRlMXdFx/vBBoA4cDfwFhgerF2K4ubyUVqov+dnXGiToVJUIaiKIt66DVsNopIOhAIBAC/ioj91oADgWWlXGMm8DuwxZlCFeZAGUoVQ9O0VzVNkysc3S82/wdwWESs+2e0RY9QDtpdrwXQDFhTwmvNALoBg0Wk0JnvS2EO1CxP1eNt4NMrtDl68bENsNPufFsgUUQsdufuBtaJyHn7C2ia9hYwDOghIocrJ1nhLihDqWKIyEngZDmbt0EfbLXSlqIGA3p352P7E5qmzUQ3k+4i8sdVSlW4IarLoygRTdO8gEhgl93p5kCSXZs6QBdghd25WcAYYDhwRtO0+hePGi4RrjAUZSiK0miOPghrbyi7gcmapvW7+HN/YLuIHLdrMx59ZmcdkGZ3THS6YoXhqDwUxVWjadpyYJOIxBitRWEOVISiqAybALU/qcKGilAUCoXDUBGKQqFwGMpQFAqFw1CGolAoHIYyFIVC4TCUoSgUCoehDEWhUDgMZSgKhcJh/H8UoYRL2wfbcAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt\n", "from matplotlib.patches import Circle\n", "from pingouin import circ_r\n", "\n", "%matplotlib inline\n", "\n", "# Convert angles to unit vector\n", "z = np.exp(1j * alpha_rad)\n", "\n", "# Plot circle\n", "fig, ax = plt.subplots(figsize=(4, 4))\n", "circle = Circle((0, 0), 1, edgecolor='k', facecolor='none', linewidth=2)\n", "ax.add_patch(circle)\n", "ax.axvline(0, lw=1, ls=':', color='slategrey')\n", "ax.axhline(0, lw=1, ls=':', color='slategrey')\n", "ax.plot(np.real(z), np.imag(z), 'bo', mfc='none', ms=12)\n", "\n", "# Plot mean resultant vector\n", "r = circ_r(alpha_rad)\n", "phi = circmean(alpha_rad)\n", "zm = r * np.exp(1j * phi);\n", "ax.arrow(0, 0, np.real(zm), np.imag(zm), width=0.01, head_width=0.1, \n", " head_length=0.1, fc='r', ec='r')\n", "\n", "# X and Y ticks in radians\n", "ax.set_xticks([])\n", "ax.set_yticks([])\n", "ax.spines['top'].set_visible(False)\n", "ax.spines['right'].set_visible(False)\n", "ax.spines['left'].set_visible(False)\n", "ax.spines['bottom'].set_visible(False)\n", "ax.text(1.2, 0, '0', fontsize=14, verticalalignment='center')\n", "ax.text(-1.3, 0, '$\\pi$', fontsize=14, verticalalignment='center')\n", "ax.text(0, 1.2, '$+\\pi/2$', fontsize=14, horizontalalignment='center')\n", "_ = ax.text(0, -1.3, '$-\\pi/2$', fontsize=14, horizontalalignment='center')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Inferential statistics" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Alpha:\tz = 4.06, p = 0.015\n", "Beta:\tz = 2.53, p = 0.078\n" ] } ], "source": [ "from pingouin.circular import circ_rayleigh\n", "\n", "# Rayleigh test for uniformity\n", "print('Alpha:\\tz = %.2f, p = %.3f' % circ_rayleigh(alpha_rad))\n", "print('Beta:\\tz = %.2f, p = %.3f' % circ_rayleigh(beta_rad))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Correlations" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Circular-circular correlation:\tr = 0.67, p = 0.007\n", "Circular-linear correlation:\tr = 0.64, p = 0.017\n" ] } ], "source": [ "from pingouin import circ_corrcc, circ_corrcl\n", "\n", "# Circular-circular correlations of alpha and beta\n", "print('Circular-circular correlation:\\tr = %.2f, p = %.3f' % circ_corrcc(alpha_rad, beta_rad))\n", "\n", "# Circular-linear correlation of alpha with range(20)\n", "print('Circular-linear correlation:\\tr = %.2f, p = %.3f' % circ_corrcl(alpha_rad, np.arange(alpha_rad.size)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Part II | Neuron firing\n", "\n", "The dataset provides the orientation tuning properties of three neurons recorded from the primary visual cortex of awake macaques. The number of action potentials is modulated by the orientation of a visual stimulus. The main variables are (1) the stimulus orientations (spaced 22.5 degrees apart) and (2) the number of spikes fired in response to each orientation of the stimulus." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
OrientationN1SpikesN2SpikesN3Spikes
00.0632510
122.566155
245.079125
367.517120
490.0101122
\n", "
" ], "text/plain": [ " Orientation N1Spikes N2Spikes N3Spikes\n", "0 0.0 63 25 10\n", "1 22.5 66 15 5\n", "2 45.0 79 12 5\n", "3 67.5 171 2 0\n", "4 90.0 101 12 2" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from pingouin.datasets import read_dataset\n", "\n", "# Load Berens (2009) neuron dataset\n", "df = read_dataset('circular')\n", "\n", "df.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Descriptive statistics\n", "\n", "Warning: the scipy.stats circular functions do not accept binned angle data." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Circular mean:\t2.36\n" ] } ], "source": [ "from pingouin import circ_axial, circ_mean\n", "\n", "# Convert the orientation to radians.\n", "ori = circ_axial(np.deg2rad(df['Orientation'].values), 2)\n", "spacing = np.diff(ori)[0]\n", "\n", "# We will only focus on the first neuron.\n", "spk = df['N1Spikes'].values\n", "\n", "# Circular mean angle\n", "print('Circular mean:\\t%.2f' % circ_mean(ori, spk))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Inferential statistics" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "z = 42.83, p = 0.000\n" ] } ], "source": [ "z, pval = circ_rayleigh(ori, spk, spacing)\n", "print('z = %.2f, p = %.3f' % (z, pval))" ] } ], "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.6.5" } }, "nbformat": 4, "nbformat_minor": 2 }