Correlation of Powerspectrum.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['seaborn-dark', 'seaborn', 'seaborn-white', 'seaborn-ticks', 'dark_background', 'Solarize_Light2', 'grayscale', 'seaborn-pastel', 'seaborn-poster', '_classic_test', 'seaborn-paper', 'seaborn-muted', 'seaborn-deep', 'bmh', 'ggplot', 'seaborn-talk', 'seaborn-dark-palette', 'seaborn-colorblind', 'classic', 'seaborn-darkgrid', 'fast', 'seaborn-notebook', 'seaborn-bright', 'fivethirtyeight', 'seaborn-whitegrid']\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import langevin\n",
    "import lmfit as lm\n",
    "from lmfit.models import ExponentialModel\n",
    "from scipy.stats import pearsonr\n",
    "\n",
    "SMALL_SIZE = 16\n",
    "MEDIUM_SIZE = 18\n",
    "BIGGER_SIZE = 20\n",
    "\n",
    "plt.rc('font', size=SMALL_SIZE)          # controls default text sizes\n",
    "plt.rc('axes', titlesize=SMALL_SIZE)     # fontsize of the axes title\n",
    "plt.rc('axes', labelsize=MEDIUM_SIZE)    # fontsize of the x and y labels\n",
    "plt.rc('xtick', labelsize=SMALL_SIZE)    # fontsize of the tick labels\n",
    "plt.rc('ytick', labelsize=SMALL_SIZE)    # fontsize of the tick labels\n",
    "plt.rc('legend', fontsize=SMALL_SIZE)    # legend fontsize\n",
    "plt.rc('figure', titlesize=BIGGER_SIZE)  # fontsize of the figure title\n",
    "\n",
    "SEED = 350107321 # from random.org\n",
    "np.random.seed(SEED)\n",
    "\n",
    "print(plt.style.available)\n",
    "plt.style.use('seaborn-white')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "# function to calculate A and B from the dataset\n",
    "def OUanalytic(data):\n",
    "    N = data.size\n",
    "    data1sq = data[0]**2\n",
    "    dataNsq = data[-1]**2\n",
    "    datasq = np.sum(data[1:-1]**2)\n",
    "    datacorr = np.sum(data[0:-1]*data[1:])\n",
    "    coef = [(N-1)*datasq,\n",
    "       (2.0-N)*datacorr,\n",
    "       -data1sq-(N+1)*datasq-dataNsq,\n",
    "       N*datacorr]\n",
    "    B=np.roots(coef)[-1]\n",
    "    Q=(data1sq+dataNsq)/(1-B**2)\n",
    "    Q=Q+datasq*(1+B**2)/(1-B**2)\n",
    "    Q=Q-datacorr*2*B/(1-B**2)\n",
    "    A = Q/N\n",
    "    P2A = -N/A**2/2\n",
    "    Btmp = B**2*(1+2*N)\n",
    "    tmp = (1+Btmp)*(data1sq+dataNsq) + (2*Btmp + N + 1 -B**4*(N-1))*datasq - 2*B*(1+B**2+2*N)*datacorr\n",
    "    P2B = -tmp/((1-B**2)**2*(data1sq+dataNsq + (1+B**2)*datasq - 2*B*datacorr))\n",
    "    PAB = (N-1)*B/A/(1-B**2)\n",
    "    dA = np.sqrt(-P2B/(P2A*P2B-PAB**2))\n",
    "    dB = np.sqrt(-P2A/(P2A*P2B-PAB**2))\n",
    "    return A,dA,B,dB\n",
    "\n",
    "def OUresult(data,deltat):\n",
    "    A, dA, B ,dB = OUanalytic(data)\n",
    "    tau = -deltat/np.log(B)\n",
    "    dtau = deltat*dB/B/np.log(B)**2\n",
    "    return A,dA,tau,dtau"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "A,D = 1.0,1.0\n",
    "delta_t=0.01\n",
    "M=1000 # number of data sets\n",
    "N=10000 # length of data set\n",
    "P=500 # range to fit acf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "# let us first calculate the standard deviation for each point of the acf\n",
    "# all acfs are stored in acf_all for later analysis\n",
    "\n",
    "ps_all = None\n",
    "for i in range(M):\n",
    "    data = langevin.time_series(A=A, D=D, delta_t=delta_t, N=N)\n",
    "    \n",
    "    ps = np.abs(np.fft.fft(data))**2\n",
    "\n",
    "    if ps_all is None:\n",
    "        ps_all = ps\n",
    "    else:\n",
    "        ps_all = np.vstack((ps_all, ps))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "freqs = np.fft.fftfreq(data.size, delta_t)\n",
    "idx = np.argsort(freqs)\n",
    "\n",
    "ps_std = np.std(ps_all, axis=0)\n",
    "ps_mean = np.mean(ps_all, axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/hstrey/anaconda/envs/python3/lib/python3.5/site-packages/ipykernel/__main__.py:1: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future\n",
      "  if __name__ == '__main__':\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD6CAYAAACvZ4z8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3X9UVHX+P/Dn5YcJCOj2wyz5CLWbkkLASYu2FvxZmhLbD90+EB9Ww4+nVdYtIdet0KzWPkuasP7CSDDYLdQW03bXZE4gatvuWbBdO665nRBkS/wmovyIGWbu9w+4w/yEGebOnZk7z8c5dPJ9h7nvK/i873nd931fQRRFEUREpFoBnu4AERG5F4OeiEjlGPRERCrHoCciUjkGPRGRygV5ugOWvvvuO5w+fRo33ngjAgMDPd0dIiKvp9frcenSJUybNg2jR4+22u51QX/69GlkZGR4uhtERD6nsrISd999t1W71wX9jTfeCKC/wzfffLOHe0NE5P2++eYbZGRkGPPTktcFvVSuufnmmzFx4kQP94aIyHfYK3fzYiwRkcox6ImIVI5BT0Skcgx6IiKVY9ATEakcg56ISOUY9EREKqeuoE9N7f8iIiIjdQU9ERFZYdATEamc1y2B4JJTpzzdAyIir6OuoO/o8HQPiIi8Dks3REQqx6AnIlI5Bj0Rkcox6ImIVI5BT0Skcgx6IiKVU2fQcxkEIiIj1cyjn7u5FufyDxn//IPNtTj6bKrnOkRE5CVUMaKfu7kW59q6AEEwfp1r68LczbWe7hoRkcepIujPtXU51U5E5E9UEfRERGQfg56ISOUUDfqamhr8+Mc/xvz58/Hyyy9DFEUld09E5JcUC/rm5masX78e27Ztw6FDh3D27Fl8/PHHSu2eiMhvKTa98ujRo1iwYAFuueUWAMDmzZsxatQopXZPROS3FB3R6/V6PPPMM0hLS0NlZSUiIyOV2j0Rkd9SLOj1ej1OnDiBjRs3oqqqCv/85z9x4MABt++3urHV7fsgIvJmigX9DTfcgOTkZFx//fUYPXo05syZg3/84x+yvLcwxLb1H3wuyz6IiHyVYkE/c+ZMnDhxAh0dHdDr9Th+/DimTp0qy3sPNXfnSo9Oln0QEfkqp4Neo9EgMTHRqr2qqgrz5s1DfHw8lixZgsbGRrPtd911F5YvX46MjAwsWLAAkyZNwpIlS0becxO3jg2R5X2IiNTIqaBvaGhAXl6eVXt1dTUKCgqQlpaG4uJihIeHY9myZWhpaTF73eOPP47Dhw/jyJEjWLt2LQRhqKKL4/IenCzL+xARqZFDQa/VarF7925kZWUhKMh8RqYoiigqKsLixYuxcuVKpKSkYMeOHRg3bhzKy8vd0mlL6Ym3KrIfIiJf5FDQHzt2DCUlJcjPz0dmZqbZtvPnz6O1tRWzZs0ytgUHByM1NRX19fXy9paIiJzmUNDHxcVBo9EgKyvLqtzS1NQEAJg0aZJZe1RUlHHuPBEReY5Dd8aOHz/e7rbOzk4AQFhYmFl7WFgYDAYDenp6MGbMGBe6SERErnB5eqW0MJnlSN9eu7vY24syeyci8l4uB314eDgAoKvL/CEf3d3dCAgIQGhoqKu7cIi9ufRcH5OI/J3LQS/V5i2nUra0tCAmJkaxET0REdnmctBHR0djwoQJqKmpMbbpdDrU1tYiOTnZ1bcnIiIXubxMsSAIyMnJwcaNGxEZGYmkpCRUVFSgvb0d2dnZMnSRiIhcIct69BkZGejt7cXevXtRVlaG2NhYlJaWIioqSo63JyIiFzgd9KtWrcKqVaus2pcuXYqlS5fK0ikiIpIPHw5ORKRyfhH0L1T/09NdICLyGNUEfeAQ0zgr/9KsYE+IiLyLaoL+yXuiANH27VG8aYqI/Jlqgv6V9DhPd4GIyCupJuiJiMg2dQW9ndINEZE/U1fQc1kdIiIr6gp6Jj0RkRWVBT0REVli0BMRqRyDnohI5Rj0REQqx6AnIlI5vwl6LmxGRP7Kb4KeC5sRkb/ym6DnPbNE5K9UFfRje656ugtERF5HVUG/vqaE690QEVlQVdCnn6nzdBeIiLyO0w8Hd0Vubi7Onj2L0aNHAwBWrlyJuXPnKtkFIiK/o2jQnz59GgcPHkR4eLiSuyUi8muKlW7a2trQ2dmJn//851i0aBGKi4thMBiU2j0Rkd9SLOi//fZbJCcno7CwEO+99x7+9re/oaqqSqndA+BNU0TknxQL+tjYWGzduhXf+973EBoaiqysLNTW1iq1ewBABW+aIiI/pFjQnzp1Ch9//LHxz319fQgKUvQSARGRX1Is6HU6HV577TV0dnZCq9Xi97//vVtm3ASKrPsTEZlyOug1Gg0SExOt2quqqjBv3jzEx8djyZIlaGxsNNs+ffp0LF68GE888QQWLVqE+Ph4pKWljbzndjzZ+EfeNEVEZMKpoG9oaEBeXp5Ve3V1NQoKCpCWlobi4mKEh4dj2bJlaGlpMXtdTk4O/vSnP+HIkSN47rnnIAjyP+P1lZpdsr8nEZEvcyjotVotdu/ejaysLKu6uiiKKCoqwuLFi7Fy5UqkpKRgx44dGDduHMrLy93SaSIicpxDQX/s2DGUlJQgPz8fmZmZZtvOnz+P1tZWzJo1y9gWHByM1NRU1NfXy9tbIiJymkNBHxcXB41Gg6ysLKtyS1NTEwBg0qRJZu1RUVFobm6GXq+Xp6cyqW5s9XQXiIgU5VDQjx8/HhERETa3dXZ2AgDCwsLM2sPCwmAwGNDT0+NiF+WVt++Up7tARKQol6dXigMzXCxH+vbalRAwxKwbHWdfEpGfcTnopQXKurq6zNq7u7sREBCA0NBQV3fhtP9u/JBTLImIBrgc9FJt3nIqZUtLC2JiYjwyoucUSyKiQS4HfXR0NCZMmICamhpjm06nQ21tLZKTk119eyIicpHLi80IgoCcnBxs3LgRkZGRSEpKQkVFBdrb25GdnS1DF+WXsfsTVObwJERE/kGWVcUyMjLQ29uLvXv3oqysDLGxsSgtLUVUVJQcby+7E19e9nQXiIgU43TQr1q1CqtWrbJqX7p0KZYuXSpLp4iISD6qejg4ERFZU23Q//CrBvtTLDn1koj8iGqDvnJfgae7QETkFVQb9MPJ2P2Jp7tARKQI/wx6QcCJf3/r6V4QESnCP4OeiMiPqDroMxsO88IrEfk9VQf9cGvevPA/LyvUEyIiz1F10A9JEFBxs/VDzomI1MZ/gx4AoPzKmkRESvPzoOejBYlI/VQf9LdevWR/oyBg3TsnlesMEZEHqD7o8+rKh5x50x0comBviIiUp/qgTz9T5+kuEBF5lOqD3hHVd6Z6ugtERG7jF0E//mqb/fKNIGD9nOXKdoiISEF+EfSf7lw25PYrIREK9YSISHl+EfQOETinnojUiUE/YG52kae7QETkFh4J+k2bNmHNmjWK7nPIBc4EAeduiuGonohUSfGgr6+vR3V1tdK7HXaBM6PUVLf2g4hIaYoG/bfffoutW7dixYoVSu7WYRlPbACOH/d0N4iIZKVY0IuiiHXr1mHt2rWIiPDMLJfhyjcnYpIAvV7ZThERuZliQV9WVoYpU6bg7rvvVmqXVhwu3wQFubcjREQKUizRPvzwQ/T29qK2thYdHR3o7u7Ghg0bUFBQoFQXHBK/qhL/KM7wdDeIiGSjWNDv37/f+P/vv/8+Tp486ZGQz2w4jIqkhbZn2AgCrko3T40dC1y5omzniIjcwOnSjUajQWKi9ZOZqqqqMG/ePMTHx2PJkiVobGyUpYNyc7h809HBGThEpApOBX1DQwPy8vKs2qurq1FQUIC0tDQUFxcjPDwcy5YtQ0tLi833efTRR1FYWDiyHitgyuqq/v+p48qXROT7HAp6rVaL3bt3IysrC0EWFypFUURRUREWL16MlStXIiUlBTt27MC4ceNQXl7ulk67arjZN9+NCkF1bIqynSIichOHgv7YsWMoKSlBfn4+MjMzzbadP38era2tmDVrlrEtODgYqampqK+vl7e3Mhm2fCMIyJufa/x/IiJf5lDQx8XFQaPRICsrC4JF8DU1NQEAJk2aZNYeFRWF5uZm6L10Xvp1fdoht+uCRnFUT0Sq4FDQjx8/3u5NTp2dnQCAsLAws/awsDAYDAb09PS42EX3eP1PRUM+YtBqVM8Ls0Tko1y+YUocCEvLkb69dm/hyCMGzUb1dXUMeyLySS4HfXh4OACgq6vLrL27uxsBAQEIDQ11dRdu88OvGoYd1f/ywZ8p1yEiIjdwOeil2rzlVMqWlhbExMR47YgeACr3DX/DVs+okME/cFRPRD7I5aCPjo7GhAkTUFNTY2zT6XSora1FcnKyq2/vdkGGvmFfk/HEhsE/cG49EfkYl4NeEATk5OTg3XffxZYtW1BXV4dnnnkG7e3tyM7OlqGL7lX44ZvDlm9OxCRZtRER+QpZVq/MyMhAfn4+Dh48iNzcXFy7dg2lpaWIioqS4+3dKv1MHcZfbRs67NG/2JmZsWPd2CsiIvkIojhMwinswoULmD17NjQaDSZOnOjcN7sw0o7OPzT094simv5vkVUbEZGnDZebfDj4gIieqxzVE5EqMegHDLsG/cASxmYXZjs6WK8nIq/HoDfhyLz6EzFJ5mEPcMolEXk1Br0JR+bVS2Fvtg4O59cTkRdj0Ft481Dh8BdZTdfBkTDsichLMegtOLIGDtC/Do5VCYc3UxGRF2LQ2zDkg0kktm6kGmgnIvImDHobHH6uLIB7VpRaNzLsiciLMOjtcLRWfzHiJrww539tbmPNnoi8AYPejvQzdQjVOvDQFEFARdJC20+jYs2eiLwAg34Irx3Z5tgyB4KA1Qufs7uNiMiTGPRDSD9T59iFWQAQBHz/2QN2tyE1laUcIvIIBv0wXqnZ5fAsnL6gUbbr9UB/GefUKYY9ESmOQe+AV2p2AaJh+BcO1Ovthn1HB3D8uLydIyIaBoPeQW8e3uxwCcfuxVkA0Os5I4eIFMWgd1D6mTr8oO0r1y/OSurquMwxESmCQe+Eo2W5ToV9dP4h+yN7oL+UQ0TkZgx6Jx0tyx1+OWOgvzwjCFi9aI39mr30uqAgeTtJRGSCQT8ClfsKENSndapmP2TYS3V7IiI3YNCP0L83P+b4i4e7QGvyOl6kJSK5KRb0oihi06ZNWLBgARYsWIDKysrhv8nL9a+H48C0S8BYxklcVTl04NfVcXRPRLJSLOg1Gg2++OILHD58GAcOHEBpaSnOnTun1O7dIv1MHd489AZgMDhcxmkPjcTqRWus17K38VrOyiEiOSgW9HPmzMGuXbsQEBCAy5cvo6+vD6GhoUrt3m3Sz9Sh6Tdpzn3TwFr2Q9btgcGHjzPwicgFitbog4ODUVhYiAULFiA5ORm33HKLkrt3K4eWNTblaN0e6A98hj0RjZDiF2PXrFmDTz75BF9//TV+97vfKb17t3FqATTJQN1+2DIOMDi6lxZH40VbInKQYkF/9uxZfPHFFwCA0NBQzJs3D2fPnlVq94p4pWZX/8je0Zo94HgZR1JX179ezqlTI+8oEfkVxYL+3Llz2LBhA3Q6HbRaLY4ePYrp06crtXvFjLRmX5G0EHOzixx7vV7Pcg4ROczpoNdoNEhMTLRqr6qqwrx58xAfH48lS5agsbHRbPvChQuRkJCARx55BI899himT5+ORYsWjbznXm4kZZxzN8UMv2yCKV6sJSIHOBX0DQ0NyMvLs2qvrq5GQUEB0tLSUFxcjPDwcCxbtgwtLS1mr8vLy8Mf//hHHDp0CCtXrnSt517O4XXsTZksm+Dw6B4wr98TEVlwKOi1Wi12796NrKwsBFmsyyKKIoqKirB48WKsXLkSKSkp2LFjB8aNG4fy8nK3dNpXmNXsnTEwup+yusq575NutgoKYugTkZFDQX/s2DGUlJQgPz8fmZmZZtvOnz+P1tZWzJo1y9gWHByM1NRU1NfXy9tbH5R+pg5vHn4D0Pc5Pbr/blQIovMPOTYrx5RezweTE5GRQ0EfFxcHjUaDrKwsCBa35zc1NQEAJk2aZNYeFRWF5uZm6PV6eXrqw9LP1KGpMH3EpZwTMUmIzv/A8dq9xfdzdE/k3xwK+vHjxyMiIsLmts7OTgBAWFiYWXtYWBgMBgN6enpc7KJ6jKhuDwwEdoDztXsJSzpEfs3l6ZXiQGhZjvTttfu7Ec21lwzU7mOe+8PIdi6VdMaOZeAT+RGXgz48PBwA0NXVZdbe3d2NgIAAVaxnIzdprr3DT6syJQgQA4MQnX8IMfkfOH6jlamOjsFRPqdmEqmey0Ev1eYtp1K2tLQgJiaGI/ohHC3LHdnofqD2LgoBqEhaiHtWlI68Ex0d/SUdaZTPkT6R6rgc9NHR0ZgwYQJqamqMbTqdDrW1tUhOTnb17VVPGt2Pv9rm/OgeAAQBFyNuGtnsHIl0p21dXf8Xw55IVVwOekEQkJOTg3fffRdbtmxBXV0dnnnmGbS3tyM7O1uGLvqHT3cuc6l2Pzg7x4k7a+3hw0+IVEWWtW4yMjKQn5+PgwcPIjc3F9euXUNpaSmioqLkeHu/IcfoXrqzNjr/0PBPs3Lw/YyzdTjSJ/JJgiiOJFHc58KFC5g9ezY0Gg0mTpzo3DeraBT6wpz/RUXSQtePSRSR2XAYr9TskqdjkpQUoLZW3vckohEZLjf5cHAvNeI595YGVsZ0ejmF4dTVDV7EJSKvFjT8S8hTpFF4ReLDg2WUkTBZTgEAwrQ9ePXINqSfcXGZBOkiriD0j/AlHOkTeRWO6L3cKzW70PSbNLx5qBDBut6Rj/ClE4UgoOu6UGMdf0Tz8G2RZuxw1g6R12HQ+4j0M3U4t/mx/nLOSGbmWBoI/Yqkhbgt76DrM3VMSbN2pC/O0SfyKAa9jzEd4Tu9IqYtggBDQKBxhD/52QPyhj5gPkc/iNVCIqUx6H2U2YqY0gjfldAfGH33Bl+H1YvW4PvPHpCvs6b0evMpm2PH8oIukZtxeOXjXqnZZbxoG7+qEldDIlyfkikI6Asahej8Q/hB21c4WpYrQ08tWK6ZHxQE3H//4J95QZdINgx6FflHcUb//HtXZ+kAxu+VnmMrufXqJeTVlbs+Y8eSZfBb1vMZ/EQjxhumVKo6NgWrFz7neuDbMvAr88OvGlC5r0De97YnMLB/xM/AJ7LCG6b8lLScwg+/anD9gq0li7V1XF5qwRHSiN90Ng8ROYRBr3KV+wrw5qFC3NrRBogyXLQ1NRC47aGRxlk7I3oClov7t5rCyWmcRGZYo/cD6WfqzGrqGU9swImYpP4/yDUytqjpC6KILYffkL+Wb480hVMSFNT/KSAwEOjrU6YPRF6KQe+HpLq6bAunmRp4L3FgFc3Vi9aYbR7XcxUFNSXuPwFID6WXpnNKIiMH///KFff2gchLsHTjx6Tn14b0dg+WdOSs55uWVizKPIqWeEx1dAx+cRlm8hMc0fs5y7KOW0b5lgYecm46bRNQeBaPKdOpnZYLtEk424d8GIOezEg3X/0+cQH0gskHPrmD38b7nYhJwtzsIvfcoOWMOhtlJcvyD8s+5EMY9GTF9G5bwE0Xb22xMdJ32w1arpDKPpZ4AiAvxaCnYUnllOrYFKyfsxxXQiIGN7p5pN8aeZPVRV2vDH/A9gmAT+IiL8A7Y2nEqmNTsPrhXwABgYONSv0MTH5tFZvJI6fISCAhgScBksVwualo0O/Zswf79++HIAiYNm0aXn75ZYwaNcqpDg+JQe9R96woxcWImzzzcxj4NZbt6VmeEBnZ/6kAGLwgXFs7OBuIJwWyY7jcVKx089lnn+H999/Hvn37EBISgvz8fOzduxdPP/20Ul0gN/t05zJUx6Ygb34udEGjlA38gX1JT88yLfUIEJHR8KH8D0iXmxTygPksIInpcs78NEBOUCzoIyIi8OKLLyI0NBQAMHnyZHz99ddK7Z4UIk3XtFnPl7j7BGDx/iL6n6RVkbSwfzNEiBC8t9Zvj+WJwNbfo7T42/HjwJgxvDhMADxUo29ra8Pjjz+OTZs24b777jPbxtKNus3NLsK5m2IGGzz9M7P49fe58B+pwMDBu4e5MqjP85rSjaS1tRU5OTl48sknrUKe1M90jrxZmUeidPAPM8snUDTgycY/en/Zx1lSyEv/b+8TAjA4bdTyWgGvHfgMRYP+zJkzWLFiBZ5++mk89dRTSu6avJDNu3Klh6YAnhvtm+xXLwSalX1MZTYcVt8JwBbLaaOWPxdbPydOK/UqipVuLl++jEWLFuHFF1/EQw89ZPd1LN0QYGfOvsRbfs42/ukEGfpQ+OGb6i/9yEEqGQHAqVOD7QkJ5m3SnwHrkwc/VQBww/RKjUaDNWvWoLGx0ay9qqoKb731Fr755hvExsZi7dq1SExMNG7fvHkzysvLER0dbWxLTU3FL37xC6c6PCRvCQByG8Xu0nWFnX9SflP/V4LlNYYxY8wvVtt7jelsJVsnCTlPHAqehGQN+oaGBjz99NMQRdEs6Kurq/HLX/4SP/vZzxAXF4d33nkHDQ0NOHjwIKKiomTt8JC89R8+yU6Ru3Tdwc4/N5+86ctfBQ7cICg970A6mUg3wVl+OrG1dpLppxnJqVMjnjYrS9BrtVqUl5dj69atCA0NhU6nMwa9KIqYPXs2HnjgAWzYsAEAoNPp8NBDD2HmzJl44YUXZO3wkHzhHzq5hc0Lu4Bv/U7Y+KfIE4AfGkE1XZZZN8eOHUNJSQny8/Nx5coV7Nmzx7jt/PnzaG1txaxZs4xtwcHBSE1NRX19vdMdJhoJywu7gMXFXdPAF0XvPAHY6JO0fr/lA1yu69Pi9T8V8QRADnEo6OPi4qDRaBAREYHi4mKzbU1NTQCASZMmmbVHRUWhubkZer0egYGBIFKatApndWwKfpPyP/hPxA245er/w8x/f4qKhPn9a/R4Y+BbstHH3uDrbJ4AbPHYOv/kNRwK+vHjx9vd1tnZCQAICwszaw8LC4PBYEBPTw/GjBnjQheJXGNrtC+dAH754M/QMyqkv1EUAQEY+I/3nwQc7N+JmCSrh7zY4tPrBNGQXJ5HL5X4Bcvbzu20E3kLWycAifQpoDXiBhiD35Kv/G472E9b6wSZ8pv7BlTI5aAPDw8HAHR1deGGG24wtnd3dyMgIMC4tg2RLxnqJABYLNHsK4HviCGOxd6NY8PhtFLPcznopdp8S0uLWZ2+paUFMTExHNGTKpku3iaN/AUREG39vqvl38AIj8PWw2OGZXLB3N7MI8trLzyZ2Ody0EdHR2PChAmoqanB/QPzQnU6HWpra5Eq3TBApFKOjPxt3uGrlvB3lLPHa/J6ezOPTF/n6MkkQBSR3NSIhltjjddmAkQR/934IYDBZyWrbY0jl4NeEATk5ORg48aNiIyMRFJSEioqKtDe3o7s7GwZukjku2ydCGyGv+WUT387EQzHkb8PB15jEIT+O6tNXmsQhMGS1ED7UGscuVvYS3/Gqz+OQ3rirbK9pyyLmmVkZKC3txd79+5FWVkZYmNjUVpa6vRdsUT+YMSfAuzhScE5jpbXPPT32qXVY/V7/XfXyhX2fGYskQ+rjk3B6oXPWd8URqrQtOlhh143XG4GyN0xIlJO+pk6NP0mDW8eKsStHW0QRAPGdncgtLe7vxxk+WXKVhupkuIPHiEi+Q1XDgJsz1IBYH7TmKP46cGnMOiJ/IS9k4GzUxIduplMwhOCV2DQE5FTHPn0APQvKidNV5QIogiIBogBNta/4knBbRj0ROQW0qJyjjAtKwX16cyWm76uTwttUDAie66hK3i09VLUQ/Hhk0fmvf8l23sx6InI4xz9lOAM05NHZM81XB0VCkOgdeSFaHvQEzzae04KgoDMe/8Lr6THyfaWDHoiUiV3nDwU4YaZUJxeSUSkcgx6IiKVY9ATEakcg56ISOW87mKsXq8HAHzzzTfOf3OQ1x0OEZFzLlxw+lukvJTy05LXJeOlS5cA9K+I6bTbbpO5N0RECps9e8TfeunSJbMHQEm8bvXK7777DqdPn8aNN96IwEAbd88REZEZvV6PS5cuYdq0aRg9erTVdq8LeiIikhcvxhIRqZxqgr6qqgrz5s1DfHw8lixZgsbGRk93yS00Gg0SExPN2kRRxI4dO5Camoq77roLP/3pT/Hll196qIfy0Ov12LNnD+bPn4+EhAQsWLAAFRUVkD6AqvGYtVottmzZgpkzZyIhIQFZWVn4/PPPjdvVeMymtFot5s+fj7Vr1xrb1HjM7e3tmDx5stVXbm4uADcds6gCf/jDH8QpU6aIxcXFYm1trbhs2TIxMTFRbG5u9nTXZPX3v/9dTExMFBMSEszai4uLxbi4OLG8vFysqakRH3vsMfH+++8Xr1696qGeuq6oqEicNm2auH37dvHkyZNiUVGRGBsbK5aUlIiiqM5jXr9+vZiYmChWVlaK9fX14vLly8WkpCTxwoULoiiq85hNvfHGG+Idd9whPv/888Y2NR7zyZMnxTvuuEOsr68XGxsbjV9fffWVKIruOWafD3qDwSDOnDlTfOmll4xtWq1WnDVrlrhx40YP9kw+vb29YklJiTh16lRx+vTpZkF/7do1MSEhQdy1a5ex7cqVK2JiYqL49ttve6K7LtPr9WJiYqK4ZcsWs/b169eL9957ryqP+erVq+LUqVPN+t/T0yPGx8eL27ZtU+Uxm/r888/FhIQE8Z577jEGvVqPec+ePeJ9991nc5u7jtnnSzfnz59Ha2srZs2aZWwLDg5Gamoq6uvrPdgz+Rw7dgwlJSXIz89HZmam2bbPPvsM3d3dmG0yJSsyMhIzZszw2eO/du0a0tPTMW/ePLP2mJgYXL58GX/5y19Ud8whISGoqqrCo48+amwLCgqCIAjQarWq/DlL+vr6sG7dOixbtgzjx483tqv1mM+ePYvJkyfb3OauY/b5oG9qagIAq7mjUVFRaG5utnsDgS+Ji4uDRqNBVlYWBIulVKXjj4qKMmufOHGicZuviYyMxEsvvYQ777zTrP3jjz/GzTffjIsXLwJQ1zEHBQXhzjvvRGRkJAwGA1paWrBu3ToIgoC0tDRV/pwlu3fvhk6nw/Lly83a1XrMZ8+eRU9PD37yk59rUpxJAAAESUlEQVQgLi4OP/rRj7B7926Ioui2Y/a6G6ac1dnZCQAICwszaw8LC4PBYEBPTw/GjBnjia7JxnSUY6mzsxOjRo3CqFHmD2MICwsz/t2owb59+3Dy5Em88MILqj/m7du3o7i4GACQm5uL2267DUePHlXlMX/55ZfYuXMnysrKrI5NjT9ng8GAL7/8EiEhIXj++ecxYcIE1NXVYfPmzejt7UVwcLBbjtnng14cmIVhOdK11642oijaPUa1HPsHH3yAgoICPPjgg8jMzMSuXbtUfcxz5szBjBkz8Omnn2L79u3Q6XQYPXq06o7ZYDDgV7/6FR5//HGrmWSAOn+3RVHEzp07ccsttxirEPfeey+6u7vx1ltvYcWKFW45Zp8v3YSHhwMAurq6zNq7u7sREBCA0NBQT3RLMeHh4dBqtdDpdGbtXV1dxr8bX1ZWVob8/HykpqaisLAQgiCo/pinTJmCGTNmYNWqVXjqqadQWlqKkJAQ1R3zO++8g//85z/Izc1FX18f+vr6APSHYV9fnyp/zoGBgUhOTrYqNT/wwAPo6elx28/Z54Ne+gtraWkxa29paUFMTIzPnvkdNWnSJIiiiAsWCyFduHABMTExHuqVPDZv3oxf//rXeOSRR1BUVGT8OKvGY7506RIOHDhg9fE8NjYWWq0WkZGRqjvmmpoaXLx4ETNmzMDUqVMxdepU/Otf/0J1dTWmTp2KoKAg1R3zxYsX8d577+Hy5ctm7b29vQDgtp+zzwd9dHQ0JkyYgJqaGmObTqdDbW0tkpOTPdgzZSQmJuK6664zO/6Ojg789a9/9enjLy8vx65du5CVlYVNmzYhyGRlUjUe89WrV7Fu3TocOXLErP3EiRO4/vrrMWfOHNUd84YNG7B//36zr+joaMycORP79+/Hww8/rLpj1mq1eOmll/DBBx+YtR85cgTR0dGYO3euW47Z52v0giAgJycHGzduRGRkJJKSklBRUYH29nZkZ2d7untuFxYWhszMTGzduhUBAQGIjo7Gzp07MWbMGDzxxBOe7t6ItLW1obCwEHfccQcefvhhfPbZZ2bbp02bprpjvv322/Hggw/i9ddfh06nQ1RUFD766CMcPHgQr732GsaMGaO6Y77Nxmqzo0ePxtixYxEX1/9gbLUdc1RUFBYuXIitW7dCEATcfvvt+POf/4yPPvoI27Ztc9u/Z58PeqB/SePe3l7s3bsXZWVliI2NRWlpqdUUJbV69tlnERAQgLfffhvd3d1ITEzEpk2bfLaOefz4cWi1WnzxxRdYsmSJ1fZPPvlEdccMAK+//jp++9vfoqSkBG1tbfj+97+PrVu34qGHHgKgvp+zI9R4zK+++iq2b9+O8vJyXLp0CbfffjuKi4uNc+fdccxcvZKISOV8vkZPRERDY9ATEakcg56ISOUY9EREKsegJyJSOQY9EZHKMeiJiFSOQU9EpHIMeiIilfv/HyuAs0SUAI8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1a15cecda0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.errorbar(freqs[idx][N/2:],ps_mean[idx][N/2:],yerr=ps_std[idx][N/2:],fmt=\"o\",ecolor='red')\n",
    "plt.yscale('log')\n",
    "#plt.xscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(-0.01890757228690829, 0.55036351028500929)\n"
     ]
    }
   ],
   "source": [
    "print(pearsonr(ps_all.T[9000],ps_all.T[9001]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [default]",
   "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.5.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}