sim_effect_size.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.stats import pearsonr\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "\n",
    "def simulate_sample(true_effect, n, seed=None):\n",
    "    rng = np.random.default_rng(seed=seed)\n",
    "    cov = np.array([[1, true_effect],\n",
    "                    [true_effect, 1]])\n",
    "    #discovery\n",
    "    y, yhat = zip(*rng.multivariate_normal(mean=[0,0], cov=cov, size=n))\n",
    "    stat1 = pearsonr(x=y, y=yhat)\n",
    "\n",
    "    # replication\n",
    "    y, yhat = zip(*rng.multivariate_normal(mean=[0,0], cov=cov, size=n))\n",
    "    stat2 = pearsonr(x=y, y=yhat)\n",
    "\n",
    "    return stat1[1]/2, stat2[1]/2   # one-sided"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.01\n",
      "0.016681005372000592\n",
      "0.027825594022071243\n",
      "0.046415888336127774\n",
      "0.0774263682681127\n",
      "0.1291549665014884\n",
      "0.21544346900318834\n",
      "0.3593813663804626\n",
      "0.5994842503189409\n",
      "1.0\n"
     ]
    },
    {
     "data": {
      "text/plain": "      p_rep      n   true_r2\n0  0.803995  359.0  0.016681\n1  0.809403  129.0  0.046416\n2  0.805627   46.0  0.129155\n3  0.858525   16.0  0.359381\n4  1.000000   10.0  1.000000",
      "text/html": "<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>p_rep</th>\n      <th>n</th>\n      <th>true_r2</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>0</th>\n      <td>0.803995</td>\n      <td>359.0</td>\n      <td>0.016681</td>\n    </tr>\n    <tr>\n      <th>1</th>\n      <td>0.809403</td>\n      <td>129.0</td>\n      <td>0.046416</td>\n    </tr>\n    <tr>\n      <th>2</th>\n      <td>0.805627</td>\n      <td>46.0</td>\n      <td>0.129155</td>\n    </tr>\n    <tr>\n      <th>3</th>\n      <td>0.858525</td>\n      <td>16.0</td>\n      <td>0.359381</td>\n    </tr>\n    <tr>\n      <th>4</th>\n      <td>1.000000</td>\n      <td>10.0</td>\n      <td>1.000000</td>\n    </tr>\n  </tbody>\n</table>\n</div>"
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "alpha=0.05\n",
    "n_iter = 1000\n",
    "true_effects = np.logspace(-2,0,10)\n",
    "ns = np.logspace(1,3,10).astype(int)\n",
    "\n",
    "df = pd.DataFrame(columns=['p_rep', 'n', 'true_r2'])\n",
    "\n",
    "for true_r in true_effects:\n",
    "    print(true_r)\n",
    "    for n in ns:\n",
    "        p1s = np.zeros(n_iter)\n",
    "        p2s = np.zeros(n_iter)\n",
    "        for iter in range(n_iter):\n",
    "            p1, p2 = simulate_sample(true_r, n=n)\n",
    "            p1s[iter] = p1\n",
    "            p2s[iter] = p2\n",
    "        p_rep = np.sum(p2s[p1s<alpha]<alpha)/np.sum(p1s<alpha)\n",
    "        if p_rep>0.8:\n",
    "            df.loc[len(df.index)] = [p_rep, n, true_r*true_r]\n",
    "            break\n",
    "\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "[<matplotlib.lines.Line2D at 0x7f984ac49b20>]"
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEKCAYAAADjDHn2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAi1klEQVR4nO3deXRU9d0/8Pdksm9kn+zLsIawyRLWJEhFkID8RFArClXwQdlFW+vPB4rU0v6UAqYgi9SNKq32Z5USClYwCfsaiBDWZJJMEpKQfd9m7vOHTxlRycxNZubOzH2/zvEcMjB3Pud43ry58733exWCIAggIiLqgpPUAxARke1jWRARkVEsCyIiMoplQURERrEsiIjIKJYFEREZ5Sz1AJai1+uh03XvqmClUtHt9xI5GuZBPlxclPf8PYctC51OQG1tc7fe6+fn2e33Ejka5kE+goN97vl7/BqKiIiMYlkQEZFRLAsiIjKKZUFEREbZxQJ3c3MzXn/9dbi4uCAxMREPP/yw1CMREcmKZGcWr776KsaOHYvp06ff9XpWVhamTJmCyZMnY+fOnQCAr776ClOmTMEbb7yBw4cPSzEuEZGsSVYWs2bNwq5du+56TafTYd26ddi1axfS09Oxb98+3Lx5E+Xl5QgLCwMAKJX3vg7YHErqWtDc3mnRzyAisjeSfQ01atQoFBcX3/VaTk4OYmJiEBUVBQBITU3FoUOHoFKpUFZWhvj4eOj1epOOr1Qq4OfnKXqueR9nw1nphE8WJMK9ixtUiORCqXTqVpbIsdjUmkV5eTlCQ0Pv/KxSqZCTk4Onn34av/3tb5GRkYH777/fpGN196a8BaOj8fKXl/HK3y9i7dT+UCgUoo9B5Eh4U558dHVTnk2Vxb14enri97//vVU+K6VPIFZM6oO3D99E/xBvPDki0iqfS0Rky2zq0tn/fN30H+Xl5VCpVFafY3FKb9zfNwhvZ+bjVGGN1T+fiMjW2FRZDB48GAUFBdBqtWhvb0d6ejomTZpk9TmcnBRYO7U/4gI98dq+KyiubbH6DEREtkSysli1ahWeeOIJaDQaJCcn47PPPoOzszPWrFmDhQsXYtq0aXjooYfQt29fSebzdFViw8wECABe/vIymtt1ksxBRGQLFIIgOOTewx0dOrPsOnuqsAbL//+3mNgnCH+YEc8Fb5IdLnDLB3ed7YHRMf5YnqzG4RuVeO9UkdTjEBFJgmVhgidHROCh+BBsP1aIrLwqqcchIrI6loUJFAoF/u/kvohXeWPN/qvQVPGUnIjkhWVhIncXJd58eCDcnJ3w8peX0dDKLUGISD5YFiKE+rrj/80YiNK6Vvz3/ivQ6R3y2gAioh9hWYg0LLIXfjmpN45rarDtWIHU4xARWYVdbPdha2YNDce1iiZ8eFqLfsFeeHBAiNQjERFZFM8suunlSb0xLMIX6w5ex7WKRqnHISKyKJZFN7konfCHGQPRy90Zv/zyMmqbO6QeiYjIYlgWPRDo5Yq3Ziagqqkdr+7LRafOtGdtEBHZG5ZFDw0M9cFrD/bDWW0dNmfmSz0OEZFFcIHbDKYNVOFaRSM+OVeCfiHeeHhQqPE3ERHZEZ5ZmMmyZDUSo/3wh69v4NKteqnHISIyK5aFmTg7KfC76fEI9nbDr/bmorKxTeqRiIjMhmVhRn4eLvjjzAQ0tnXiV3tz0d7JBW8icgwsCzPrE+yFtVP749tbDXjz0E046ONCiEhmWBYWMKlfMJ4dE40vL5Xhswu3pB6HiKjHWBYWsmhcDJLUAdiYkYdz2lqpxyEi6hGWhYU4KRRYN20Aovzc8et/XsGt+lapRyIi6jaWhQV5uzljw8wEdOr1+OWXuWjt0Ek9EhFRt7AsLCwmwBNvTIvH9YpG/PbgdS54E5FdYllYwXh1ABZPiMVX125j95liqcchIhKNZWEl8xOjMLl/MLYc0eC4plrqcYiIRGFZWIlCocDqKf3QJ9gLr6VfQVFNi9QjERGZjGVhRR4uSmyYmQBnJye8/MVlNLZ1Sj0SEZFJWBZWFt7LHX+YEY+immb85l/XoOeCNxHZAZaFBEZE+eHFib2RlVeFd48XSj0OEZFRLAuJPHZfOGYkqLDrZBEO36iUehwioi6xLCSiUCjw6wf6YlCYD9b+6ypuVjZJPRIR0T2xLCTk6uyENx8eCE9XZ7z8xWXUtXRIPRIR0U9iWUgs2NsNbz48EBWNbXgt/Qo69VzwJiLbw7KwAUPCffHKz/rgVGEttmRppB6HiOhHnKUegL4zc3AYrlU04eNzxeiv8sJD8SqpRyIiuoNnFjZk1UQ1hkf2wu++uoEr5Q1Sj0NEdAfLwoY4K53whxnx8PdwwS+/zEV1c7vUIxERAWBZ2Bx/T1dsmJmA2pYO/HpvLjp0eqlHIiJiWdii/ipvrH6wH7JL6vHHb/KkHoeIiAvctmpKfAiu327ER2eK0T/EG48MCZN6JCKSMZ5Z2LDFE+IwNtYfbx66iYsldVKPQ0QyxrKwYUonBd5IHYAwXzf8am8uyhvapB6JiGSKZWHjfN1dsOH/JKC1Q49f7c1FWycXvInI+lgWdkAd6IV10/ojt6wBv//3dQh8BgYRWRnLwk6k9AnCf42LQXpuBd746joqm3gPBhFZD6+GsiMLxkSjuV2HPedL8O9rt/H0qCg8NTISHi5KqUcjIgenEBz0O42ODh1qa5u79V4/P89uv9caimpasPWIBodvVCLIyxWLxsVgxqBQKJ0UUo9GDsjW80DmExzsc8/fY1n8BHsJx8WSOrydqcG3t+qhDvTE8mQ1xsX5Q6FgaZD52EseqOdYFiLZUzgEQcA3Nyqx5YgG2tpWjIz2w4rkOAxQ3ft/OpEY9pQH6hmWhUj2GI4OnR6fX7yFd08Uoq61Ew/Fh+CFCbEI83WXejSyc/aYB+oeloVI9hyOxrZOfHBai7+eL4EgCHj8vgg8MzoaPu68loG6x57zQOKwLERyhHCU1bdi+7EC7M+tgK+7M54dE405w8LhouTV0iSOI+SBTMOyEMmRwnGtohFpmfk4XVSLiF7uWJIUhwf6BXERnEzmSHmgrrEsRHK0cAiCgBMFNfhTlgY3K5swKMwHK5LVGBbZS+rRyA44Wh7o3lgWIjlqOHR6Aem55dh+rAC3G9sxsU8gliTFITbAU+rRyIY5ah7ox1gWIjl6OFo7dPjkXAk+PK1FW6cOjwwJw3PjYhDg6Sr1aGSDHD0PZMCyEEku4ahubse7xwvxj5xbcHdRYt6oKDw5IgLu3D6EvkcueSCWhWhyC0dBdTO2HtEg42YVQrxdsWh8LFIHqrh9CAGQXx7kjGUhklzDkV1ch7SsfFy61YA+QV5YnhKHsbEBUo9FEpNrHuSIZSGSnMMhCAK+vl6JrUc0KKlrxegYPyxLVqN/iLfUo5FE5JwHuWFZiMRwAO2devz9YineO1mE+tZOTBsYgufHxyKU24fIDvMgHywLkRgOg4bWTrx/qgh/yy6BQqHAE8Mj8IvEKHi7cfsQuWAe5INlIRLD8WO36lux7WgB/nWlAn4eLlg4JhqzhoZx+xAZYB7kg2UhEsNxb1fLG/B2lgZni2oR5eeOpUlxuL8vtw9xZMyDfLAsRGI4uiYIAo5rapCWlY/8qmYMCffF8uQ4DI3g9iGOiHmQD5aFSAyHaTr1AvZdKsOO44WobGrHpL5BWJIUh2h/D6lHIzNiHuSDZSESwyFOS4cOH58txkdntGjXCZg9NAwLx8TAz9NF6tHIDJgH+WBZiMRwdE9VUzvePVGIL/53+5BfJEbhieHcPsTeMQ/ywbIQieHoGU1VM7Yc0SAr77vtQ16YEItpA1Vw4iK4XWIe5INlIRLDYR7ntLVIy9Igt6wBfYO9sCJZjdGx/lKPRSIxD/LBshCJ4TAfvSDg62u3sfWIBqX1bRgT64/lyXHoG8ztQ+wF8yAfLAuRGA7za+/U47MLpXjvVBEaWjsxPUGF58fHIsTHTerRyAjmQT5YFiIxHJZT39qB905q8emFEjgpFJg7IgJPj+L2IbaMeZAPloVIDIflldS1YNvRAhy8ehv+Hi54blwMHhkcCmduH2JzmAf5YFmIxHBYT25ZA97OzMf54jpE+3tgWVIcUvoEcvsQG8I8yAfLQiSGw7oEQcDR/Gr8KUsDTXUzhkX4YkWKGoPCfKUejcA8yAnLQiSGQxqdegF7L5Vhx7ECVDd34IF+wViSFItIP24fIiXmQT5YFiIxHNJqbtfhL2e12H2mGJ16AbOHhWPBmGj4eXD7ECkwD/LBshCJ4bANlY1t2HG8EHsvlcHTVYlnEqPx+PAIuDlzEdyamAf5YFmIxHDYlrzKJmw5osHR/GqE+rjhhQmxmBofwu1DrIR5kI+uyqLLf6Jdu3YNjz32GFJSUrB69WrU1dXd+b3Zs2ebb0KiLvQO8sKmRwZh25wh8PNwwW/+dQ3z/5KNM0U1Uo9GJBtdlsXatWuxdOlS/POf/0RsbCyefPJJFBUVAQA6OzutMiDRf4yM9sOHT92HddP6o7alA4s/+xYrP7+EvMomqUcjcnhd3jbb1NSE5ORkAMCCBQuQkJCAhQsX4s033+R18CQJJ4UCD8WrMKlvMD7NLsF7p4rw5EfnMGNQKJ4fF4Mgb24fQmQJRvdYaGhogI/Pd99jjRkzBmlpaVi+fPldX0kRWZubsxOeHhWFGYNC8d7JInx2oRQHr1TgqZGReHpUFDxd+QwNInPq8muo5557Dnl5eXe9NmDAAHzwwQeYPHmyRQcjMoWfhwtW3d8bnz0zEkm9A7HrZBEe+fNpfH6xFJ16h7x2g0gSRq+G0ul02LBhA1555RVrzWQWvBpKni7dqkdaZj6yS+oRF+CJpclxSFIH8GvTHmAe5KPbV0MBgFKpxLlz58w6EJGlDArzxY7Hh2LDzIHQCwJe+uIynv80B5fLGqQejciumXSfxW9+8xuUl5dj6tSp8PT0vPP6gw8+aNHheoJnFtSp0+OLb8vw7olCVDd34MH+wVicFIuIXtw+RAzmQT66OrMw6SEC7e3t8Pf3x6lTp+563ZbLgshZ6YTZw8Lx0MAQfHSmGB+fLcY3NysxZ1g4nh0djV7cPoTIZGa5g3vHjh1YtGiROeYxG55Z0A9VNLRhx/EC/PNSOXzcnfHM6Gg8Niwcrtw+pEvMg3z0aM3CFAcOHDDHYYgsKsTHDaun9Mcn80YgIdQHb2fmY877Z3DwSgX0jrnrDZHZmKUsHHR7KXJQfYK9kPboYGyZPRjebs747/1X8YuPs3FOWyv1aEQ2yyxlwcsSyR6NjvHH7qeHY+3U/qhu7sDzn+Zg1T8uQVPFr1yIfohnFiRrTgoFUhNU+PszI7E0KQ7ni+vw8w/P4vf/voGqpnapxyOyGWYpi6lTp5rjMESScXdRYn5iFL5YkIjZw8Lx5aUyzPrzGbx/qghtnXqpxyOSnElXQ2k0GqxduxZVVVXYt28frl69isOHD2Px4sXWmLFbeDUU9URRTQv+lJWPjJtVCPN1w9KkOEzuHyzLr1yZB/no8dVQq1evxksvvQRn5+9uyxgwYAD2799vnumIbFC0vwfempmA7Y8Nga+7C15Lv4oFey7g29J6qUcjkoRJZdHS0oIhQ4bc9ZpSyV09yfGNiPLDh3Pvw5op/XCrvg3P7rmA1/Zdwa36VqlHI7Iqk+7g9vf3R1FR0Z1T8AMHDiA4ONiigxHZCqWTAjMGheJn/YKx+4wWu88WI+NmJZ4cEYn5iVHwdjMpRkR2zaQ1C61Wi9WrVyM7Oxu+vr6IjIzEW2+9hcjISGvM2C1csyBLKW9owztHNdifW4EATxcsGh+LmYNCoXRyzPUM5kE+ulqzELXdR3NzM/R6Pby9vc0ymCWxLMjScssasDkjD9kl9egT5IWVKWqMjvWXeiyzYx7ko8dlsWXLlp98fenSpd2fysJYFmQNgiDgmxuVSMvSoKSuFePjArAiRY24QE/jb7YTzIN89PhqKE9Pzzv/KZVKHDlyBCUlJWYbkMheKRQKTOoXjE9/MRIrUtS4WPrdTX1vHrqJ2uYOqccjMptu7Trb3t6OBQsWYPfu3ZaY6Ue0Wi22bduGxsZGpKWlmfQenlmQFGqa2/HuiSJ8frEUHq5KLBgTY/c72zIP8mH2XWdbWlpQVlZm0p999dVXMXbsWEyfPv2u17OysjBlyhRMnjwZO3fu7PIYUVFRWL9+fXdGJbIqf09X/OpnfbBn/kgMDe+FtzPz8fiHZ3H4RiW3xSG7ZtI1fzNmzLjza71ej+rqaixZssSkD5g1axaeeuqpu57hrdPpsG7dOrz//vtQqVSYPXs2Jk2aBJ1Oh40bN971/vXr1yMwMNCkzyKyFXGBntg8axBOFlRjc2Y+Xtmbi/sifLFyYm8MDL33v96IbJVJZbF9+3bDG5ydERgYeOdubmNGjRqF4uLiu17LyclBTEwMoqKiAACpqak4dOgQFi1ahB07dpg6O5HNGxMbgL9E+2PvpTLsOFaA+R9nI3VgCF6YEAeVj5vU4xGZzOjf+DqdDgsWLDDrA47Ky8sRGhp652eVSoWcnJx7/vmamhps2rQJubm5Jj+VT6lUwM+ve1ekKJVO3X4v0U95Nrk35iTGYOeRfLx3vACHblRi4fg4PJcUB09X276pj3kgwISyUCqViIuLQ2lpKcLDw60x04/4+/tj3bp1ot6j0wlc4Cabs2BUJB7qF4StRzTYkpGHv57RYvGEWKQmqOBko5sUMg/y0dUCt0n/pKmvr0dqaiqGDBkCDw+PO69//+spMVQq1V0L5OXl5VCpVN06FpG9Ce/ljt9Nj8fjwyOwKSMP6w5ex9+yS/HiRDVGRPlJPR7RTzKpLNra2u5aSxAEARs2bOj2hw4ePBgFBQXQarVQqVRIT0/HH//4x24fj8geDQn3xXs/H4Z/X7uNP2Vp8PynOZjYJxDLktWI9vcwfgAiKzKpLHQ6HRITE+96rbXVtF03V61ahdOnT6OmpgbJyclYtmwZ5syZgzVr1mDhwoXQ6XR49NFH0bdvX/HTE9k5hUKBBweEILl3IPacL8EHp7R47IOzeGxYOBaOjYavu4vUIxIBMHJT3ieffII9e/ZAq9XeuXIJAJqamjB8+PAenV1YGm/KI3tU1dSO7ccKsPdSGXzcnLFwbAxmDw2Ds1K6m/qYB/no9t5QDQ0NqKurw8aNG/HSSy/ded3Lywt+fn5mHdLcWBZkz27cbsTmjHycLqpFtL8HVqSokaQOkORJfcyDfJht11l7wrIgeycIAo5pqvF2Zj4KqlswMtoPK1PU6B9i3V2fmQf5YFmIxHCQLenU6fF5Thl2Hi9AfWsnHh4UiufHxyDI2zo39TEP8sGyEInhIFtU39qB905q8bfsErgoFZifGIW5IyLh7mLZRxwzD/LBshCJ4SBbVlzbgrQsDb65UQmVjxuWJMViyoAQi93UxzzIB8tCJIaD7MH54lpszsjHlfJGJIT64MWJagyN6GX2z2Ee5INlIRLDQfZCLwj4V24Fth7V4HZjOx7oF4SlyXGI6GW+m/qYB/lgWYjEcJC9aenQ4S9nivHRGS10goCfD4/AM6Oj4e3W800KmQf5YFmIxHCQvapoaMO2YwVIv1wOPw8XLBofg5mDw+Ds1P31DOZBPlgWIjEcZO+uljdgU0Y+zhfXIS7QEytT1BgXF9CtYzEP8sGyEInhIEcgCAIyb1YhLSsf2tpWjI31x4oUNXoHeYk6DvMgHywLkRgOciQdOj0+u1CKXSeK0NTeiUeGhOG/xsUgwNPVpPczD/LBshCJ4SBHVNvSgV0nCvH3C6Vwd1Hi2dHReHx4BNycu96kkHmQD5aFSAwHObKCqmakZeXjSH41wnu5Y1lSHH7WL+iemxQyD/LBshCJ4SA5OFVYg7cz83HjdhOGhPti1UQ1EsJ8f/TnmAf5YFmIxHCQXOj0AvZdLsM7RwtQ3dyBqfEhWDIhFqG+7nf+DPMgHywLkRgOkpum9k58dFqLj8+VAADmjojAvMQoeLk6Mw8ywrIQieEguSqrb8XWowU4cKUCgV6ueGF8DJ4ar0ZDfYvUo5EVsCxEYlmQ3F26VY9NGfnIKa3HAJUPliXFIjHGX+qxyMJYFiKxLIi+u6nv0PVKbD1agOLaFiSpA7A8RY3YAE+pRyMLYVmIxLIgMvDwcsOOb27ivVNFaO3UY/bQMCwcGwM/DxepRyMzY1mIxLIgMvhPHqqb27HzeCH+kXMLXq7OWDg2GnOGhcNF2fVNfWQ/WBYisSyIDH6Yh7zKJmzOzMfJghpE+bljebIaKX0C73lTH9kPloVILAsig3vl4bimGpsz86GpasaIqF5YmaLGANW9/7Ih28eyEIllQWTQVR469QK+yLmFHccLUdfSgdQEFRZPiEWwt5uVpyRzYFmIxLIgMjAlD41tnXj/VBH2nC+BUqHAvMQoPDUyEh4uSitNSebAshCJZUFkICYPxbUt2HpEg6+vVyLE2xWLJ8ThoYEhcOJ6hl1gWYjEsiAy6E4eLpbUYWNGPnLLGhCv8sbKiWoMj/SzzIBkNiwLkVgWRAbdzYNeEHDwagW2ZGlQ0diO+/sGYXlyHCL9PCwwJZkDy0IklgWRQU/z0Nqhw8fnivHhaS06dAIevy8CC8ZEw8fd2YxTkjmwLERiWRAZmCsPlY1t2H6sEHsvlcHX3RmLxsfikSFhcHbieoatYFmIxLIgMjB3Hq5VNGJzRh7OausQF+iJlSlqjIsLMNvxqftYFiKxLIgMLJEHQRCQlVeFzZn5KK5txbg4f6xM6Y24QG5SKCWWhUgsCyIDS+ahvVOPTy+U4s8nC9HSrsOjQ8Px3DhuUigVloVILAsiA2vkoeZ7mxR6cpNCybAsRGJZEBlYMw95lU3YnJGPk4U1iPb3wIoUNZLUAdyk0EpYFiKxLIgMrJ0HQRBwXFODzZl5KKhuwahoP7w4UY2+wd5Wm0GuWBYisSyIDKTKQ6dOj89zbmHn8UI0tHXi4UGheH58LAK9XK0+i1ywLERiWRAZSJ2H+tYO7DpRhE8vlMLd2QnPjI7GE8Mj4ObM9QxzY1mIJHU4iGyJreShoLoZaZn5OJJfjfBe7lieHIdJfYO4nmFGLAuRbCUcRLbA1vJwqrAGmzPycbOyCfdF+OLF+3sjng9dMguWhUi2Fg4iKdliHnR6AV9eKsP2owWobenAtAQVlvChSz3GshDJFsNBJBVbzsMPH7o0/38fuuTOhy51C8tCJFsOB5G12UMeimtb8KcsDQ7fqITKxw1Lk+IwZUAw1zNEYlmIZA/hILIWe8rD+eJabPwmH9cqGjEozAerJvbG4HBfqceyGywLkewpHESWZm950AsC0i+XY+vRAlQ1tWPKgGAsTYpDqK+71KPZPJaFSPYWDiJLstc8NLfr8OEZLT4+WwwAmDsyEvNHRcHTlesZ98KyEMlew0FkCfaeh7L6Vmw5osHBq7cR5OWKxRNikZqgghPXM36EZSGSvYeDyJwcJQ85pfXYlJGHS7caMCDEGy/er8bwSD+px7IpLAuRHCUcRObgSHnQCwK+unobW45oUN7Qhkl9g7AsOQ6Rfh5Sj2YTWBYiOVI4iHrKEfPQ2qHDX84W48PTWugEAT8fHoFnRkfD281Z6tEkxbIQyRHDQdRdjpyHioY2vHOsAOmXy+Hv4YLnJ8Ri5qBQKJ3kuZ7BshDJkcNBJJYc8pBb1oBNGXm4UFKPPkFeeHGiGokx/lKPZXUsC5HkEA4iU8klD4Ig4PCNSqRl5qO0vg1J6gCsSFEjJsBT6tGshmUhklzCQWQKueWhrVOPv54vwfunitDaqcdjw8KxcGw0fN1dpB7N4lgWIsktHERdkWseqprasf1YAfZeKoOPmzP+a1wMZg0Jg7PScR+6xLIQSa7hIPopcs/D9YpGbMrMx9miWsQGeGDlxN4YHxcg9VgWwbIQSe7hIPo+5uG79YysvGqkZeWjqKYFY2L9sTJFjd5BXlKPZlYsC5EYDiID5sGgQ6fHZxdKsetEEZrbO/HIkDAsGhcLP0/HWM9gWYjEcBAZMA8/VtvcgZ0nCvH5xVJ4uCqxcEwMHrsvHC52vp7BshCJ4SAyYB7uLb+qCZsz8nGioAZRfu5YnqxGSp9Au33oEstCJIaDyIB5MO64phqbM/KhqW7GyKheeHFib/QL8ZZ6LNFYFiIxHEQGzINpOvUCPr94CzuPF6C+tRMPDw7F8+NjEeTlKvVoJmNZiMRwEBkwD+LUt3bgzyeL8LfsUrgpnfCL0VF4ckQk3Jxtfz2DZSESw0FkwDx0T1FNC9Iy85GZV4VwXzcsS1bjZ/2CbHo9g2UhEsNBZMA89MyZohpsysjHjdtNGBbhixcn9sbA0Hv/pSwlloVIDAeRAfPQczq9gL2XyrD9WAGqmzuQOjAEiyfEIcTHTerR7sKyEInhIDJgHsynsa0T75/SYs/5YigVCsxLjMLTIyPh7qKUejQALAvRGA4iA+bB/ErqWrAlS4Ovr1cixNsVS5LiMDU+BE4Sr2ewLERiOIgMmAfLyS6uw6aMPFwpb0RCqA9enKjG0Iheks3DshCJ4SAyYB4sSy8I2J9bjneOFuB2Yzsm9w/GsuQ4hPm6W30WloVIDAeRAfNgHS0dOnx0WovdZ4shCALmjozE/MQoeLk6W20GloVIDAeRAfNgXWX1rdh6tAAHrlQg0MsVi8fHIjVBBaWT5dczWBYiMRxEBsyDNC7dqsfGb/Lw7a0G9Av2wqr7e2NElJ9FP5NlIRLDQWTAPEhHEAT8+9ptpGVpUN7Qhol9ArEiRY1IPw+LfB7LQiSGg8iAeZBea4cOH58rxoentejUC3j8vggsGBMNbzfzrmewLERiOIgMmAfbcbuxDduOFmDf5XL4ebhg0fgYzBwcBmczrWewLERiOIgMmAfbc7W8ARsz8pFdXAd1oCdWTeyN0bH+PT4uy0IkhoPIgHmwTYIg4JsblXg7S4PSulZMUAdgRbIasYGe3T4my0IkhoPIgHmwbe2devwtuwR/PlmE1k49Xhgfi/mJUd06VldlYb27PYiIyOxcnZ3w9KgopCaosPN4IW7Vt1rkc1gWREQOIMDTFb9+oK/Fjm/7z/kjIiLJsSyIiMgolgURERnFsiAiIqNYFkREZBTLgoiIjGJZEBGRUSwLIiIyymG3+yAiIvPhmQURERnFsiAiIqNYFkREZBTLgoiIjGJZEBGRUSwLIiIyimVBRERGsSyIiMgoPinPBFqtFtu2bUNjYyPS0tKkHodIMl9//TUyMjLQ2NiI2bNnY8KECVKPRFYi2zOLV199FWPHjsX06dPvej0rKwtTpkzB5MmTsXPnTgBAVFQU1q9fL8WYRBYnJgsPPPAA3njjDbz++uvYv3+/FOOSRGRbFrNmzcKuXbvuek2n02HdunXYtWsX0tPTsW/fPty8eVOiCYmsoztZ2LZtG+bOnWvtUUlCsi2LUaNGoVevXne9lpOTg5iYGERFRcHV1RWpqak4dOiQRBMSWYeYLAiCgLfeegvJyclISEiQaGKSgmzL4qeUl5cjNDT0zs8qlQrl5eWoqanBmjVrkJubix07dkg4IZF13CsLu3fvxokTJ3DgwAHs2bNHwgnJ2rjAbQJ/f3+sW7dO6jGIJDdv3jzMmzdP6jFIAjyz+B6VSoWysrI7P5eXl0OlUkk4EZE0mAX6IZbF9wwePBgFBQXQarVob29Heno6Jk2aJPVYRFbHLNAPyfbhR6tWrcLp06dRU1ODwMBALFu2DHPmzEFmZibWr18PnU6HRx99FC+88ILUoxJZFLNAppBtWRARken4NRQRERnFsiAiIqNYFkREZBTLgoiIjGJZEBGRUSwLIiIyimVBRERGsSyIiMgobiRIZCXFxcV47rnnMGLECGRnZ0OlUuGdd96Bu7u71KMRGcUzCyIrKiwsxNy5c5Geng4fHx8cPHhQ6pGITMKyILKiyMhIxMfHAwASEhJQUlIi8UREpmFZEFmRq6vrnV8rlUrodDoJpyEyHcuCiIiMYlkQEZFR3KKciIiM4pkFEREZxbIgIiKjWBZERGQUy4KIiIxiWRARkVEsCyIiMoplQURERrEsiIjIqP8BAd63x14Y+egAAAAASUVORK5CYII=\n"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ax=sns.lineplot(data=df, x='n', y='true_r2')\n",
    "ax.loglog(True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "(array([0.35030195]), array([0.74831034]))"
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.stats import norm\n",
    "r = 0.5\n",
    "n=100\n",
    "\n",
    "z = np.arctanh(r)\n",
    "norm.interval([0.95], loc=z, scale=1/np.sqrt(n-3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "[84.92761044616469, 29.01223670641077, 782.6482175199623]"
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import statsmodels.stats.power as smp\n",
    "transf = lambda r, n: np.arctanh(r)+r/(2.*(n-1))\n",
    "[smp.NormalIndPower().solve_power(np.arctanh(r), alpha=0.05, ratio=0,\n",
    "            power= 0.8, alternative='two-sided') + 3 for r in [0.3, 0.5, 0.1]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "[10152.001766984604,\n 7236.050547314269,\n 5157.250542505102,\n 3675.2607405548974,\n 2618.7405019017388,\n 1865.5399624034374,\n 1328.5776423103227,\n 945.7724955166505,\n 672.8663939867668,\n 478.30720568404064,\n 339.60087555224635,\n 240.71074773510136,\n 170.20343692606727,\n 119.92685648625017,\n 84.06749879095432,\n 58.477803809399134,\n 40.19542075785188,\n 27.097329627814133,\n 17.64315068994928,\n 10.64366898990942]"
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import statsmodels.stats.power as smp\n",
    "\n",
    "true_effects = np.logspace(np.log(0.23),0,20, endpoint=False)\n",
    "\n",
    "n_rep = [smp.NormalIndPower().solve_power(np.arctanh(r), alpha=0.05*0.05*2, ratio=0,\n",
    "            power= 0.8, alternative='larger') + 3 for r in true_effects]\n",
    "n_rep"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "outputs": [
    {
     "data": {
      "text/plain": "<Figure size 432x144 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAACOCAYAAAA7IRMWAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAi4ElEQVR4nO3deXSUdZ7v8XctqUpCKvtWBUkIgewkQMIuNoK22m0L2nbrxYZp79zuuXqPc3sc5ujMMDLnCuOh2+nT03faduxFL+O02m4gARVRESGyFZB9hRQJqUpCJSErqSRVdf+oJBBIIImVpKryfZ3DMb/61fL75nn88PDU8/x+CqfT6UQIIYTPUU73AIQQQkwOCXghhPBREvBCCOGjJOCFEMJHScALIYSP8oiAdzqd2Gw25IIeIYRwH48I+N7eXoqLi+nt7R3X60pKSiZpRFPD28c/Gm+vy9vHPxpfqMsXaphKHhHwE9XT0zPdQ/hGvH38o/H2urx9/KPxhbp8oYap5NUBL4QQYnQS8EII4aMk4IUQwkdJwAshhI+SgBdCCB8lAS+EED7K6wO+r99BZW0r3T190z0UIYTwKOrpHsA3daq0gRf/3ykUCkiIDSYlIYzUhDBSEsKZHRWEUqmY7iEKIcS08PqAX7lQz//56UrKTS2U17ZytMDMJ8cvAjArwI+U+LCB0A8nOT6UoEDNNI9YCCGmxoQCft26dWg0GrRaLQBbt25lzZo1nDt3jueffx6bzcbs2bP5xS9+QUREhFsHfCOFQsHilGgWp0QD4HA4qb/cScXFFsovtlJxsZW3Pq1gcJqbuJggUuLDSZ3rOsqPi9GhkqN8IYQPmvAR/K9//WuSk5OH2g6Hg7/7u7/jxRdfJDc3l5dffpmXXnqJF1980S0DHSulUkFcjI64GB13L0sAoLunj6raK5TXtlBuauVESQOHTtUCEKBVkxwfSkpCOKkJYSTHhxESpJ3SMQshxGRw2yma4uJitFotubm5ADz22GOsX79+ygN+JIH+fmQnR5GdHAW4Zq+0WLsGjvBdR/rvfl6Fw+E6zNdHzho6j5+ZFEF8jA6FQo7yhRDeZcIBv3XrVpxOJzk5OTzzzDNYLBYMBsNQf3h4OA6HgytXrhAaGuqOsbqNQqHAEBWEISqIdblxAPTY+qm6dIWKi62Um1o4W3GZL4yXAAjVacmaH0nW/CiyF0QSGzFrOocvhBBjonBOYBJ2i8WCXq+nt7eXnTt30tXVxT333MN7773Hq6++OvS87Oxsvvzyy9sGvM1mo7i4eNyDn0xOp5PWLjumRhs1DT3UNNro7HEAEDpLRWKslsQYfxJjtOgCVNM8WiHETJWTkzNq34SO4PV6PQAajYZNmzbx5JNPsmXLFsxm89BzWlpaUCqV4zp6z8zMHPridiyMRuMti3Mnp9NJXWMHhdVWCqouU3S+mbPnuwGIi9GRvcB1hL9wfiRBAX5jes+pHP9U8va6vH38o/GFunyhhqk07oDv7u7Gbrej0+lwOp0cOHCAtLQ0MjMz6enp4fTp0+Tm5vLWW29x3333TcaYp4VCoSA+Npj42GAeuGMedoeTC/VXKKyyUlht5dOTteQdrUGpgHlzQsmeH0nWgijSE8Px13j91ahCCC807uRpbm7m6aefxm6343A4SEpKYvv27SiVSn7+85+zffv2YZdJ+iqVUsGCuDAWxIXx/XULhu6oLai6TGG1lb1HzvPeF9WoVQpSEsLJXhDF4uQoFsSHyWWZQogpMe6Aj4uLY8+ePSP2LVmyhH379n3TMXklP7WSjHkRZMyLYNO9ri9tS2taKKy+TEHVZd48WM6fPiknKMCPxSnR5KRGo7LZp3vYQggfJucOJom/Vs2S1GiWpLpuwOro7uVc5WXOlDdxpqKRr87VA/DBicMsSXUFfurccNQqr58eSAjhISTgp4guUMOaRbNZs2g2TqcTk6WdDz87Q0OHmg8OV/Pu51UEaNUsSo4iJzWaJSkxRIUFTNt4Dxvr2P1RGdbWq0SGBbDl/jTW5sRN23iEEOMnAT8NFAoFiYYQ7kgPJicnh+6ePgqqrBjLGzlT0cTXRRbAdXVOzsDRfca8CPzUU3M55mFjHb966yz2gRu/Lrde5VdvnQWQkBfCi0jAe4BAfz9WLtSzcqEep9PJpaZOjOWNGMubyDtaw54vz6PVqFiYFEluajRL02OJDg+ctPG8uqdoKNwH2R1OXt1TJAEvhBeRgPcwCsW1uXQ2fms+PbZ+is5bOVPehLGiidNljbzyQREJsTqWpseyND2GlIRwt16Z09E98tz6oz0uhPBMEvAezl+rHgjyWADqL3dyqrSRU6UNQ+fudYEactKiWZYWy+LU6DHfaCWE8G0S8F5mdlQQs78VxMZvJdF1tY+zlU2cLGngdFkTh42XUCoVZCRGsDQ9hqXpMcyJ1o37M3SBfiMeresC5S8OIbyJBLwXmxXgxx3Zs7kjezZ2h5PKi62cKmvgVGkjf9xXwh/3lWCInDV0Kic9MQI/9e0vw/zpxoX829tn6bdfOw+vVin46caFk1mOEMLNJOB9hEqpIC0xnLTEcLZ8J52mlm5OlblO5RzIr2HvkfME+qtZnBzNsowYctNiCZ418upWg1+kymWSQng3CXgfFR0eyHdXJ/Ld1Yn02PopqLo8FPjHCs0oFZCWGMGKzFiWZ+jRRw6fAnltTpwEuhBeTgJ+BvDXqlmeqWd5ph6Hw0n1pSucLGngREkDf/iwhD98WEJcjI4VmbEsy4glOS5MFisXwgdIwM8wSqWC5HjX0oQ/uj+NhuauobB/74tq3vmsijCdlmUZrrDPXhCF1k/muxfCG0nAz3CxEbN48M4kHrwzic7uXk6XNXK8pIEjZ+v55PhFtBoVi5OjWJ6hZ2l6jKxXK4QXkYAXQ4ICNUPn3vv67RRVN3OixMLJkgaOFzegVEDq3HCWZ+hZsTAWQ2TQdA9ZCHELEvBiRH5q1dBsmP/z4SzO17dxoriBEyUWXssr4bW8EhJidazI1LMiU0/SnBBZmFwIDzPugN+1axeffPIJ9fX17Nu3j+TkZABqamp47rnnhhbZ3rVrF3PnznX3eMU0UCgUzJ8Tyvw5oTx+XyqNLd2cKLbwdbGFdz6r5O1DlUSGBrAiM5YVmfqb5rERQkyPcQf8+vXr2bJlC48//viwx7dv386mTZvYsGEDe/fu5fnnn2f37t1uG6jwHDHhgUPn7ds6bZwqbeR4sYWDxy+Sd7QGf42ClVVGVi7Uszg5Gn+t/ENRiOkw7v/zcnNzb3qsubmZ0tJSXnvtNQAeeOABXnjhBVpaWggPD//moxQeKyRIy93L4rl7WTw9tn7OVjax/8sSTpU28oXxEhq1ksUp0azIlC9phZhqbjm0slgsxMTEoFK5LqdTqVRER0djsVjGFfDFxcXj/myj0Tju13gSbx//jTTAQyvDsTuc1F62UVbXQ1nNZU6UNKBQQHyUhpBAFTUNNjp6HIQEqlifHUxW4qzbvvdU8rXtMsgX6vKFGtwpJydn1D6P+rdzZmYmWu3Yj/CMRuMti/N03j7+0RiNRpYtzWXZQNvpdHK+vo3jxRY+O1nLxaarQ89t67az71QbiYmJHnPnrC9vF2+vyxdqmEpuWQBUr9fT2NiI3e5aRNput9PU1IRer3fH2wsvN/gl7Y/uS0Mxwh2yff0OfvXWWXYfKKWythWnU76kFcId3HIEHxERQVpaGnl5eWzYsIG8vDzS0tLk/Lu4ibX16oiP2x3OoTtpB6/IWbXQQHpiOCpZiFyICRl3wO/YsYODBw9itVp54oknCA0NZf/+/fzzP/8zzz33HC+//DLBwcHs2rVrMsYrvFxkWACXRwj5qLAA/u2ZtZwqbSC/8NoVObpAjevyy4V6Fi2IQiPTJggxZuMO+G3btrFt27abHk9KSuKdd95xy6CE79pyfxr//k4Btj770GNaPxVb7k9DF6hhXW4863JdV+QYK5o4XmThWKGZT0/WEqBVkZMaw6qFBnLSogn0dy1ActhYJ1MbCzECj/qSVfi+sc41769VszrLwOosA339DoqqrXxdbOF4kYWjBWb81EoWJ0cTFqzli9N19PY7ALjcepV/f6dg2GcJMVO5PeDXrVuHRqMZuhpm69atrFmzxt0fI7zYeOea91Mrh02bUG5qIb/IzNdFFk6W3ny6x9ZnZ/dHZRLwYsablCP4X//610NTGAjhTiqlgox5EWTMi+B/PJjJg1s/HPF5l1uv0tDcRWyEZ11fL8RU8qzLE0pLXf+12+H116Gw0NXu63O1B2+E6umB11/Hv6bG1e7udvVXVLjanZ2udnW1q93W5mpfuOBqt7a62iaTq221utp1da52U5OrXV/vajc0uNoNDa52fb2r3dTkatfVudpWq6ttMrnara2u9oULrnZbm6tdXQ2vv46yu9vVrqhw9Q+2y8pc7Z4eV7u42NXuG1gIu7DQ1R64LJVz51ztQUYjXD9NxKlT8MYb19rHj8Obb15r5+fD229fax89Cu++e6395Zfw/vvX2l98AXv2XGsfOgT79g01dcePw/791/o//tj1Z9D+/XDw4LX2vn2u9xi0Z4/rMwa9/75rDIPefReOHkWhUBAVFsAdFUdJNZcPdd9ZfoQUcwU/+ZdD/O9fHub4s7+g8ePD116/e7frdzTo9dddv0MAu52IffugsJCO4iPU/uonNH0/F8vf/4CO4iND+x5lZa7ne+m+R2enq+1l+17Q6dO33Pc4eHBK9r0hb7/tqmHQm2+6ahz0xhuu38Gg2+x7Y8m9m/a9W5iUI/itW7fidDrJycnhmWeeITg4eEyvq6qqwuZwuP4nM5no9vfnal8f9PcTYTLRFRBAj82GwmYj3GSCzEyMRiOKnh7CTSY6dTpsnZ0ou7sJM5noDA3F1taGsrOTMJOJjrAweltbUbW3E2oy0VFURG9zM6orVwg1mWgvLKSvqQl1Swshg+2GBtRWKyEmE20FBfTX1+PX1ETwYDs8HL+GBoJNJq4UFGAPDUVjNqMzmbhy7hz24GA0ly6hM5loPXcOR1AQ2ro6gkwmSErCaDSiNZkIMploOXsWp78//jU1zDKZaDlzBqdWi391NbNMJprPnAG1moCqKgJNJpqNRlCpCKiouNYGAsvK8K+poeX69sWLQ+1Z5eVo6utpvb7d2HitXVGBn9XKlYF2UEUF6ra2a+3KSlSdnbQNtHWVlShttmttoKq6mvaBdnBVFQDtUVGudnU1TrWajogIAEKqq3FotXSEhbna589jDwqic2C/CT1/nn6rlc6goKF2X1sbXQEBrEn1p+Xk8P1IpVSwPGUWCfNDKKu7yrFCM3++dJorx1rIVlbz7cIDBJ7R0WZ+hD5Dpmtf02i4are79j2g8vCHqLqKUPT1AtDf1UrTvt9gSq5AZ6qnKyiInu5ur9r3CAqiuLiYIJOJ1rNncQQGet++B1SfPz/6vldVhaK/f0r2PYCwCxfo7eqia+CUdNiFC/TabHT5uS4ACK+pocfhoFupvNZWKhn46/TmfW8MuXfjvjeX0Smcbr6rxGKxoNfr6e3tZefOnXR1dfHSSy/d8jU2m43i4mK5k9VHTHVdt7uKprntKl8e/IpjRhPVfVE4UBKh7CBbW8+d65ayeO1dw5YoNBqNROW/Sn+7dcTPUwdHEnbX4+gy75z02tzJF/Y3X6hhKrn9CH7w7lWNRsOmTZt48skn3f0RQgxzuy9tI0ICyG14j0U6K50OLcV9cRT0xvPl1QV8vr+D8K8OsnKhnlVZejISXUd2/e3No75ff7sV6/5XALwu5MXM4taA7+7uxm63o9PpcDqdHDhwgLS0NHd+hBATMhjYQUobK7TVrNBWc9XhR0lfHNUJP+DTk7XsP1ZD8CwN82P9WKFOIamvErXCMeL7OftttH7xXxLwwqO5NeCbm5t5+umnsdvtOBwOkpKS2L59uzs/QogJUQdH3HTKJUDZx4qodn7442VDN1blF5o5XmTmTP8yAhSLyPSrI0tTS6qfGY3CPuz1tzrKF8ITuDXg4+Li2HP9N9xCeIiwux7Huv8VnP22occUai1hd7kWrrn+xqrjJ0+jDJrDl1+e5fSFOE71JqGhj3RNPYv8LpKuqUer6EcdHDFd5QgxJnInq5gRBk+ltH7xX/S3N6MOjiDsrscx2hLZvePgsC9odSoFOemxLEu/n367g5OHvuDIkXMU2GZzrncufl39pGoauCMjmfCrfQQF+E1zdUKMTAJezBi6zDuHnTM/bKwbNi/O4DQH380NZvBCDbVKyap717Nwth/Nn/8XlS0qikim0D6P337Vze/yPyJ7QRSrsgwsz4iVFauER5GAFzPW7o/Khk16Bq5pDj4raOeJR4Y/d/Avh7nAtwGHw0llXSv5hRbyC8383z+f4zdKBZnzIlidbWBlpp6wYP+pKkWIEUnAixlrtLnp27rtIz5+PaVSQWpCOKkJ4TzxQDoX6ts4Vmgmv9DCb98r5JX3C0mbG87qLAMrFxqICguY8Dg7io/cdGpJrt4RYyEBL2as0eamDwm89Zzzo91YlTQnlM33p1Hb2DF0ZP+7vcX8bm8xyfGhrM4ysCrLMK75cTqKj3A57zdg7wdc1+BfzvsNINfgi9uTgBcz1mhz06/PHn1qjdHO24PrhiuFQkFCbDAJscH8t2+nYL7c6TqyL7LwWl4pr+WVMm92CKuy9KzOMjAnWnfLMTYf/ONQuA+x99N88I8S8OK2JODFjDXa3PQ6mkZ9zWjn7UebntgQFcQP1ifzg/XJNLZ0k19oJr/QzBsflfPGR+XEx+qGLs+Mj9WhUAxfs9ZxtWPEcYz2uBDXk4AXM9pI0xwYjaMH/Gjn7Ud7/Hox4YE8tHY+D62dj/XKVb4eWK3qrU8rePNgBbOjZrFq4DRO0uyQm8JeiPGSgBdiHEY7bx85zi9RI0MD+N6aeXxvzTxa23s4Xmwhv9AytPB4THggq7IMJCrjiLfXcWPWK/yDvkkZYoaQgBdiHG61puxEhQX7c/+qRO5flUhbp40TJQ3kF5rZ99V5+u13EarsItvvItmaWhLVTSgVCiLv/Ut3lCN8nAS8EOMw1jVlJyokSMu3lyfw7eUJdF7t48gnh/nqxCWO2VL40pZOsKKbbO0l7rT0sTTNgUrlWWv2CM8iAS/EOI13TdmJCgrwI7PuXVKDrPQ41ZT2zuFcbwLHe+bx1cedhBz9hBWZelZlGciaH4lawl7cQAJeCA82OGOlv6KfJVoTS7QmbE41ZX2zOT9/E0fOXuKT4xcJCvBjRaae1dkGshdE4qe+9bX8YmaQgBfCg400zbFW0U9uZBcPb87F1mfnbEXTwLX2Zg6dqiXQX82yjFhWZxlYnBKN1k/CfqZya8Dv2rWLTz75hPr6evbt20dycrI7316IGed20xxr/VSsyNSzIlNPX7+dgiorxwrMHC+2cNh4iQCtiqVpsazKNqDsH3nxEuG73Brw69evZ8uWLTz++OPufFshZqzRpjke6S5WP7WK3LQYctNi+F/2bIqqrRwrNPN1kYUj5+rxUylYVn6K1VkGctKiCfSXaY59nVsDPjc3151vJ4Tg5mmOx0KtUrI4JZrFKdE8+XAWJTXN7DlUQGlNM8cKzWjUSpakRrM6y8DS9FhmyZz2PsmzzsE/9hi0tl5r//CH8NRT0N0N3/nOTU+PWLsWcnLAaoVHHrmpnyefhEcfhbo62Lz55v6//Vv43vegogL+6q9u7t+2De6+G86dg5/97Ob+f/kXWLUK8vPhH/7h5v5f/QoWLYJDh2DHjpu6tU8/7Rr/vn3wr/968+v/8z8hLg7efht++9ub+999FyIj4fXXXX9udOAABAbCyy/Dn/98c//hw67/vvQS5OUN7wsIgI8+cv38wgvw2WfD+yMi4L33XD///d/D118PdSV3dEBaGrzxhuuBn/3M9Tu8XnIyvPqq6+ef/hQqK4f3L1rk+v0B/OhHcOnS8P6VK+HFF10/f//70HzD8nnr18M//ZPr5/vvh6s33Jz0wAOwdavr57Vrhw+towP+8i9vue/x4x+7/njBvqfasYMsILGjgyCdju6ePvI2P8tHdVew7/2QyDMfEhTgR0iQFt0sDWqlwmP3veSODpg7d9R9D4A5c7x23wNum3s37XuDv8sReFTAd3V10dtxbY6N1tpaLhuNKHp6WNAx8twbRqMR1ZUrJI3Qf/nCBVqNRvwaGkgcob+xupo2oxGtyUTCCP2Wyko6wsIIqKggboT++vJyurRaZpWXM3uE/rrSUq7a7egqK9HfYvwh1dXEjNBfU1hIX1MTYRcuEDVC//mCAuyhoUSYTESM0F919ixOf3+iamsJG6G/0mgEIObSJUJu6Hf09VE90B9rNhN8Q3+/UsmFgX5DQwNBN/Q3t7RgGuif09RE4A39PVYrtQP98VYr/jf0dzc1cWmgf25LC5ob+jsbGjAP9M9rbUV9Q3+72UzDQP/89naUPT3D+tsuXaJxoD95hN9N7W32vWaTiWYv3Pc6B37OiO1j3nIDdkcIujIl3T19dHT3AhCgUXLiwxPELbIwxwP3vdbW1lvue71evu/dLvdu3PduNV2dwul0Om/RPyHr1q3jlVdeGfOXrDabjeLiYjIzM9Fqx74ijtFoJGdw6R0v5O3jH4231+Xt4x/NrepyOp1U1V0hv9DMsUIzDc3dKJUKFiZFsDrLwIqFesJ007+Aia9um8niUUfwQojpoVAoSI4PIzk+jL/47vULmJh5+b1Cfvt+IRnzIgYWMNETETLxBUzE1HFrwO/YsYODBw9itVp54oknCA0NZf/+/e78CCHEJFMoFCTNCb22gElDB8cKzRwtMPMfHxTx6p4iUhPCWZ1tYNU3XK1KTC63Bvy2bdvYtm2bO99SCDGNFAoFCfpgEvTBbLo3ldqGdvKLLBwrMPP7vcX8fm8xKfFhA9Mc68e1WpWYfHKKRggxZvGxwcTHBvPYPSnUX+4cOmf/Wl4Jr+WVMH9OCKuyDKzONmCIlCmNp5sEvBBiQmZft1pVQ3PXUNjvPlDG7gNlJBqCh9ahjYu59dKEYnJIwAshvrHYiFk8fNcCHr5rAU0t3eQXmckvtPDGx+W88bFracI7sgysyjYQH3Pz0oRickjACyHcKjo8kI3fms/Gb82nue0q+YWupQnf/LSCPx2sYE50kGsd2mwDc/XBEvaTSAJeCDFpIkKGL034dbHrC9p3Pqvk7UOV6CNnDS06njRH1qF1Nwl4IcSUCAv25zurEvnOwNKEx4stHC0w8/7hat79vIro8MCBsNeTHB8mYe8GEvBCiCkXEqTl3hVzuXfFXNq7ejlR7DqNs++r83xwuJrI0ABWZem5I2s2KQlhKJUS9hMhAS+EmFbBszTcszyBe5Yn0NndO7DouIUDx0x8eOQC4cH+rMrSszrLgMPh9plVfJoEvBDCYwQFali/NJ71S+Pp7unjZGkj+YVmDh6/SN7RGmb5K7mztoDVWQYy50XIouO3IQEvhPBIgf5+rF0yh7VL5nDV1s/pskbyDpfw+ek6Pso3ERKkca1Dm2VgoSw6PiIJeCGExwvQqlmzaDaB9gYyMrMxVjSRX2AeWnRcF+hadHxVloHsBVH4qSXsQQJeCOFl/LXqoUsrr190/FihmU9P1jIrwI/lGbGszjawODkKP/XMXXRcAl4I4bVuXHT8bOVljhWYOVHSwOen6wj0V7MsPZZVWQaWpEaj9ZtZYS8BL4TwCX5qFcvSY1mWHktfv4PCalfYHy+2cPjMJfw1Kpamu47sc1Kj8df4fvz5foVCiBnHT60kJzWGnNQYnnokm6JqK8cKzXxdZOGrc/VoNSpyU2NYnWUgNz2GAK1vRqFvViWEEAPUKiWLU6JZnBLNkw9nUVLTzLECs2te+0IzGrWSnLQYVmUZWJYeQ6C/33QP2W0k4IUQM4ZKpSRrfhRZ86P46UNZlNU0DyxNaOHrIgtqlZIlKdGsztazLENPUIB3h70EvBBiRlIpFWQmRZKZFMlPNiyk/GLL0MyXJ0sbUKvOkb0gijuyDSzP1KML1Ez3kMdNAl4IMeMplQrSEyNIT4zgv38vg6q6Vo4NhP2/vX0O1TsFZM2PZHW2gRWZekKCtNM95DGRgBdCiOsolQpSEsJJSQjniQfSOX+pjaMF9eQXWvj3dwp4+d0CMpMiuSPbwIqFesJ0/tM95FFJwAshxCgUCgXz40KZHxfKX3w3nRpzu+umqoJ6Xn6vkN++X0jGvAhWZxlYuVBPREjAdA95GAl4IYQYA4VCwbzZIcybHcKP7kultqGDowWuO2j/44MiXt1TRNrc8IGwNxAVNv1hLwEvhBDjpFAoSNAHk6AP5vH7UqltaHdddllg5nd7i/nd3mJSEsK4I9vAqoUGosMDp2WcEvBCCPENxccGEx8bzGP3pHCpqWPoapw/fFjCHz4sYUFc6NA6tLERs6ZsXBLwQgjhRnOidfzwbh0/vDsZi7VraCK01/eX8vr+UpLmhAxNlmaICprUsYwp4Gtqanjuuee4cuUKoaGh7Nq1i7lz5w57ztGjR/nlL39JZWUlmzdv5tlnn52M8QohhNfQR87ikXULeGTdAhpbuskfCPvdB8rYfaCMufpg7lhk4OG18ydl1ssxBfz27dvZtGkTGzZsYO/evTz//PPs3r172HPi4uLYuXMnH3/8Mb29vW4fqBBCeLOY8EAeWjufh9bOp6m1m68Hztn/6eNyUuPDyU6Ocvtn3jbgm5ubKS0t5bXXXgPggQce4IUXXqClpYXw8PCh5yUkJABw6NAhCXghhLiF6LBANtyZxIY7k+i3OyZtNarbBrzFYiEmJgaVyvXPB5VKRXR0NBaLZVjAu0NxcfG4X2M0Gt06hqnm7eMfjbfX5e3jH40v1OULNbhTTk7OqH0e9SVrZmYmWu3YbwE2Go23LM7Tefv4R+PtdXn7+EfjC3X5Qg1T6bYBr9fraWxsxG63o1KpsNvtNDU1odfr3TYIp9MJMKFTOzabzW3jmA7ePv7ReHtd3j7+0fhCXb5Qg7tpNBoUCsVNj9824CMiIkhLSyMvL48NGzaQl5dHWlqaW0/P9PX1AVBZWTnu107ktI4n8fbxj8bb6/L28Y/GF+ryhRrcbbSzHwrn4OHzLZw/f57nnnuO9vZ2goOD2bVrF/PmzeMnP/kJf/3Xf83ChQs5ffo0zzzzDJ2dnTidTnQ6HTt37mTNmjW3HZzD4aCrqws/P78R/xYSQggxutGO4McU8EIIIbzP5FybI4QQYtpJwAshhI+SgBdCCB8lAS+EED5KAl4IIXyUBLwQQvgoCXghhPBRXhHwu3btYt26daSkpAy727WmpoZHH32Ue++9l0cffRSTyTR9g7yNp556igcffJCNGzeyadMmysrKAO+qYSLbwRPrm8i28LQ63L0tpqM+d28HT9tGHsHpBU6dOuU0m83Ou+66y1lRUTH0+ObNm5179uxxOp1O5549e5ybN2+eriHeVnt7+9DPn376qXPjxo1Op9O7apjIdvDE+iayLTytDndvi+moz93bwdO2kSfwioAfdP3ObLVanTk5Oc7+/n6n0+l09vf3O3NycpzNzc3TOcQx+eCDD5wPPfSQ19Yw1u3gDfWNZVt4ch3u2BaeUN833Q6eUIMn8qjpgsdjKuepd5d//Md/5NixYzidTn7/+997ZQ03ulUNTqfTY+sbz7bw5DquN9EaprM+d20Hb9lGU80rzsH7ip07d3L48GH+5m/+hp///OfTPZwZTbaFZ5DtMLm8NuCvn6cemJR56ifLxo0bOXHiBLGxsV5bw6BbbQdv2EZj2RbeUAdMfFt4Qn3fdDt4Qg2eyGsD/vp56oFJmafeXbq6urBYLEPtzz//nJCQEK+qYTS3qsET65vItvDEOkYy0Rqmoz53bwdv2UZTzSumC96xYwcHDx7EarUSFhZGaGgo+/fvH3Week9jtVp56qmnuHr1KkqlkpCQEJ599lkyMjK8pgaY2HbwtPomui08rQ53b4uprm8ytoOnbSNP4BUBL4QQYvy89hSNEEKIW5OAF0IIHyUBL4QQPkoCXgghfJQEvBBC+CgJeCGE8FES8EII4aMk4IUQwkf9fwdFwY4AV+xMAAAAAElFTkSuQmCC\n"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "sns.set(rc={\"figure.figsize\":(6, 2)})\n",
    "sns.set_style(\"whitegrid\")\n",
    "ax = sns.lineplot(x=n_rep, y=np.square(true_effects))\n",
    "ax.loglog(True)\n",
    "ax.axhline(0.05, color='red', linestyle=':', alpha=0.5)\n",
    "ax.axhline(0.01, color='red', linestyle='--')\n",
    "ax.set_xticks([10, 30, 100,300,1000,3000])\n",
    "ax.set_xticklabels([10, 30, 100,300,1000,3000])\n",
    "ax.set_yticks([0.001, 0.01, 0.05, 0.1, 0.5])\n",
    "ax.set_yticklabels([0.1, 1, 5,10,50])\n",
    "ax.scatter(x=[75, 75, 125, 200, 275, 375], y=[0.26, 0.25, 0.081, 0.058, 0.028, 0.021])\n",
    "ax.scatter(x=[250, 350, 325, 700, 700, 700], y=[0.042, 0.039, 0.044, 0.021, 0.022, 0.015])\n",
    "sns.despine()\n",
    "plt.savefig('fig/analitical_solution.pdf')"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "outputs": [
    {
     "data": {
      "text/plain": "0.0025000000000000005"
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.stats import binom\n",
    "binom.pmf(2, n=2,p=0.05)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "0.05*0.05"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "np.log(0.1)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "true_effects**2"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-2.3025850929940455"
      ]
     },
     "execution_count": 92,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.log(0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.01      , 0.01584893, 0.02511886, 0.03981072, 0.06309573,\n",
       "       0.1       , 0.15848932, 0.25118864, 0.39810717, 0.63095734])"
      ]
     },
     "execution_count": 95,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "true_effects**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}