{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# eda12_00 clear all variables and import vatious modules\n",
    "%reset -f\n",
    "import os\n",
    "from datetime import date\n",
    "from math import exp, pi, sin, cos, tan, sqrt, floor, ceil, log, atan2\n",
    "import numpy as np\n",
    "import scipy.sparse.linalg as las\n",
    "import scipy.interpolate as ip\n",
    "import scipy.spatial as sp\n",
    "from scipy import sparse\n",
    "import scipy.linalg as la\n",
    "import scipy.signal as sg\n",
    "import scipy.stats as stats\n",
    "from matplotlib import pyplot as plt\n",
    "import matplotlib.cm as cm\n",
    "from matplotlib.colors import ListedColormap\n",
    "\n",
    "# eda_draw function makes a \"pictorial matrix equation\"\n",
    "# arguments are vectors, matrices and strings\n",
    "# which are plotted in the order that the appear\n",
    "# except that strings starting with 'title ' are plotted\n",
    "# under the subseqeunt matrix or vector\n",
    "# always returns a status of 1\n",
    "\n",
    "def eda_draw(*argv):\n",
    "    bw = np.zeros((256,4));\n",
    "    v = 0.9*(256 - np.linspace( 0, 255, 256 ))/255;\n",
    "    bw[:,0] = v;\n",
    "    bw[:,1] = v;\n",
    "    bw[:,2] = v;\n",
    "    bw[:,3] = np.ones(256);\n",
    "    bwcmap = ListedColormap(bw);\n",
    "    # size of plot\n",
    "    W = 16;\n",
    "    H = 4;\n",
    "    fig1 = plt.figure(1);\n",
    "    # figsize width and height in inches\n",
    "    fig1.set_size_inches(W,H);\n",
    "    ax1 = plt.subplot(1,1,1);\n",
    "    plt.axis([0, W, -H/2, H/2]);\n",
    "    plt.axis('off');\n",
    "    LM = W/6;    # matrix width and heoght\n",
    "    LV = W/40;   # vector width\n",
    "    FS = 0.12;    # character width\n",
    "    TO = 0.4;    # title vertical offset\n",
    "    SP = 0.2;    # space between objects\n",
    "    LS = 0.2;    # leading space\n",
    "    p = LS; # starting x-position\n",
    "    istitle=0; # flags presence of a title\n",
    "    for a in argv:\n",
    "        if isinstance(a,np.ndarray):\n",
    "            sh = np.shape(a);\n",
    "            if len(sh) == 1:  # conversion to nx1 array\n",
    "                n = sh[0];\n",
    "                m = 1;\n",
    "                ap = a;\n",
    "                a = np.zeros((n,1));\n",
    "                a[:,0] = ap;\n",
    "            else:\n",
    "                n = sh[0];\n",
    "                m = sh[1];\n",
    "            if m==1:\n",
    "                pold=p;\n",
    "                left=p;\n",
    "                right=p+LV;\n",
    "                bottom=-LM/2;\n",
    "                top=LM/2;\n",
    "                plt.imshow( a, cmap=bwcmap, vmin=np.min(a), vmax=np.max(a), extent=(left,right,bottom,top) );\n",
    "                p = p+LV;\n",
    "                pm = (p+pold)/2;\n",
    "                if istitle:\n",
    "                    plt.text(pm,-(LM/2)-TO,titlestr,horizontalalignment='center');\n",
    "                    istitle=0;\n",
    "                p = p+SP;\n",
    "            else:\n",
    "                pold=p;\n",
    "                left=p;\n",
    "                right=p+LM;\n",
    "                bottom=-LM/2;\n",
    "                top=LM/2;\n",
    "                plt.imshow( a, cmap=bwcmap, vmin=np.min(a), vmax=np.max(a), extent=(left,right,bottom,top) );\n",
    "                p = p+LM;\n",
    "                pm = (p+pold)/2;\n",
    "                if istitle:\n",
    "                    plt.text(pm,-(LM/2)-TO,titlestr,horizontalalignment='center');\n",
    "                    istitle=0;\n",
    "                p = p+SP;\n",
    "        elif isinstance(a,str):\n",
    "            ns = len(a);\n",
    "            istitle=0;\n",
    "            if( ns>=6 ):\n",
    "                if 'title ' in a[0:6]:\n",
    "                    istitle=1;\n",
    "                    titlestr=a[6:];\n",
    "            if( istitle != 1):\n",
    "                plt.text(p,0,a);\n",
    "                p = p + ns*FS + SP;\n",
    "    plt.show();\n",
    "    return 1;\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 4200 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "area under p(e): 1.0014\n",
      "area under p(E): 0.9090\n"
     ]
    }
   ],
   "source": [
    "# eda12_01\n",
    "# distribution of square of a random variable\n",
    "\n",
    "L=101;\n",
    "emin=0.01;\n",
    "emax=3;\n",
    "De=(emax-emin)/(L-1);\n",
    "e = np.zeros((L,1));\n",
    "e[:,0] = emin+De*np.linspace(0,L-1,L);\n",
    "\n",
    "pe = sqrt(2/pi)*np.exp(-0.5*np.power(e,2));\n",
    "Ae = De*np.sum(pe);\n",
    "\n",
    "Emin=0.01;\n",
    "Emax=3;\n",
    "DE=(Emax-Emin)/(L-1);\n",
    "E = np.zeros((L,1));\n",
    "E[:,0]=Emin+DE*np.linspace(0,L-1,L);\n",
    "\n",
    "pE = np.sqrt(np.reciprocal(2*pi*E)) * np.exp(-0.5*E);\n",
    "AE = DE*sum(pE);\n",
    "\n",
    "eda_draw(' ',pe,'caption p(e)',' ',pE,'caption p(E)');\n",
    "\n",
    "print(\"area under p(e): %.4f\" % (Ae) );\n",
    "print(\"area under p(E): %.4f\" % (AE) );"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 35068 characters]\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda12_02\n",
    "# chisquared distribution\n",
    "\n",
    "L=1001;\n",
    "X2min=0;\n",
    "X2max=10;\n",
    "DX2=(X2max-X2min)/(L-1);\n",
    "X2 = np.zeros((L,1));\n",
    "X2[:,0]=X2min+DX2*np.linspace(0,L-1,L);\n",
    "\n",
    "fig1 = plt.figure(1,figsize=(10,6));\n",
    "plt.axis([X2min, X2max, 0, 1]);\n",
    "plt.xlabel('X2');\n",
    "plt.ylabel('p(xX2)');\n",
    "plt.title('chi-squared pdf');\n",
    "\n",
    "pX2 = np.zeros((L,1))\n",
    "for N in range(1,6):\n",
    "    pX2[:,0] = stats.chi2.pdf(X2.ravel(),N);\n",
    "    plt.plot(X2,pX2,'k-');\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 47072 characters]\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda12_03\n",
    "# student t distribution\n",
    "\n",
    "L=1001;\n",
    "tmin=-5;\n",
    "tmax=5;\n",
    "Dt=(tmax-tmin)/(L-1);\n",
    "t = np.zeros((L,1));\n",
    "t[:,0]=tmin+Dt*np.linspace(0,L-1,L);\n",
    "\n",
    "fig1 = plt.figure(1,figsize=(10,6));\n",
    "plt.axis([tmin, tmax, 0, 0.5]);\n",
    "plt.xlabel('x');\n",
    "plt.ylabel('p(x)');\n",
    "plt.title('students t pdf');\n",
    "\n",
    "pt = np.zeros((L,1))\n",
    "for N in range(1,6):\n",
    "    pt[:,0] = stats.t.pdf(t.ravel(),N);\n",
    "    plt.plot(t,pt,'k-');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 79756 characters]\n",
      "text/plain": [
       "<Figure size 864x864 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda12_04\n",
    "# F distribution\n",
    "\n",
    "L=1001;\n",
    "Fmin=0;\n",
    "Fmax=5;\n",
    "DF=(Fmax-Fmin)/(L-1);\n",
    "F=Fmin+DF*np.linspace(0,L-1,L);\n",
    "\n",
    "fig1 = plt.figure(1,figsize=(12,12));\n",
    "\n",
    "pF = np.zeros((L,1))\n",
    "i = 1;\n",
    "for M in  [2, 5, 25, 50]:\n",
    "    ax1 = plt.subplot(4,1,i);\n",
    "    plt.axis([Fmin, Fmax, 0, 2]);\n",
    "    plt.xlabel('F');\n",
    "    plt.ylabel(\"p(F,%d,.)\" % (M) );\n",
    "    if(i==1):\n",
    "        plt.title('F pdf for M, N in [2, 5, 25, 50]' );\n",
    "    i=i+1;\n",
    "    for N in  [2, 5, 25, 50]:\n",
    "        pF[:,0] = stats.f.pdf(F.ravel(),N,M);\n",
    "        plt.plot(F,pF,'k-');\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# eda12_05\n",
    "# generates simulated data for particle size example\n",
    "\n",
    "# note: to prevent this code from overwriting the existing\n",
    "# files testA.txt and testB.txt, I have changed the filenames\n",
    "# to testA_new.txt and testB_new.txt\n",
    "\n",
    "N=25;\n",
    "d=100;\n",
    "\n",
    "s2dtrue=1;\n",
    "sdtrue=sqrt(s2dtrue);\n",
    "\n",
    "# calibration test 1\n",
    "dobs1 = np.random.normal(d,sdtrue,(N,1));\n",
    "np.savetxt(\"testA_new.txt\",dobs1, delimiter=\"\\t\");\n",
    "\n",
    "# calibration test 2\n",
    "dobs2 = np.random.normal(d,sdtrue,(N,1));\n",
    "np.savetxt(\"testB_new.txt\",dobs2, delimiter=\"\\t\");\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Z distribution, with Z=0.2780, two sided test:\n",
      "  The probability that Z is inside  the interval -0.2780 to 0.2780\n",
      "    is 0.2190\n",
      "  The probability that Z is outside the interval -0.2780 to 0.2780\n",
      "    is 0.7810\n",
      "Chi-squared distribution with 25 degrees of freedom, one sized test:\n",
      "   The Probability that X2 is inside the interval 0 to 21.9000\n",
      "     is 0.3585\n",
      "   The Probability that X2 is outside the interval 0 to 21.9000\n",
      "     is 0.6415\n"
     ]
    }
   ],
   "source": [
    "# eda12_06\n",
    "# A few exemplary calulations of probability\n",
    "\n",
    "# Z distribution, two-sided test\n",
    "Z=0.278;\n",
    "P = stats.norm.cdf(Z,0.0,1.0)-stats.norm.cdf(-Z,0.0,1.0);\n",
    "\n",
    "print(\"Z distribution, with Z=%.4f, two sided test:\" % (Z) );\n",
    "print(\"  The probability that Z is inside  the interval %.4f to %.4f\" % (-Z, Z) );\n",
    "print(\"    is %.4f\" % (P) );\n",
    "print(\"  The probability that Z is outside the interval %.4f to %.4f\" % (-Z, Z));\n",
    "print(\"    is %.4f\" % (1-P) );\n",
    "\n",
    "# Chi-squared distribution, one-sided test\n",
    "X2=21.9;\n",
    "N=25;\n",
    "P=stats.chi2.cdf(X2,N);\n",
    "\n",
    "print(\"Chi-squared distribution with %d degrees of freedom, one sized test:\" % (N) );\n",
    "print(\"   The Probability that X2 is inside the interval 0 to %.4f\" % (X2) );\n",
    "print(\"     is %.4f\" % (P) );\n",
    "print(\"   The Probability that X2 is outside the interval 0 to %.4f\" % (X2) );\n",
    "print(\"     is %.4f\" % (1-P));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "file testA.txt\n",
      "N is 25\n",
      "true mean 100.0000\n",
      "estimated mean 100.1620\n",
      "true variance 1.0000\n",
      "true root-variance 1.0000\n",
      "estimated variance (given true mean) 0.8108\n",
      "estimated root-variance (given true mean) 0.9004\n",
      "estimated variance (given estimated mean) 0.8172\n",
      "estimated root-variance (given estimated mean) 0.9040\n",
      "Z is 0.8098\n",
      "Probability of |Z| exceeding 0.8098 is 0.4181\n",
      "chi-squared is 20.2689\n",
      "Probability of chi2 exceeding 20.2689 is 0.7326\n",
      "t is 0.8994\n",
      "Probability of |t| exceeding 0.8994 is 0.3770\n",
      " \n",
      "file testB.txt\n",
      "N is 25\n",
      "true mean 100.0000\n",
      "estimated mean 99.7455\n",
      "true variance 1.0000\n",
      "true root-variance 1.0000\n",
      "estimated variance (given true mean) 1.1273\n",
      "estimated root-variance (given true mean) 1.0618\n",
      "estimated variance (given estimated mean) 1.1068\n",
      "estimated root-variance (given estimated mean) 1.0521\n",
      "Z is 1.2727\n",
      "Probability of |Z| exceeding 1.2727 is 0.2031\n",
      "chi-squared is 28.1837\n",
      "Probability of chi2 exceeding 28.1837 is 0.2995\n",
      "t is 1.1987\n",
      "Probability of |t| exceeding 1.1987 is 0.2419\n",
      " \n",
      "significance of difference between two means, using Z-test\n",
      "Z is 1.4348\n",
      "Probability of |Z| exceeding 1.4348 is 0.1514\n",
      " \n",
      "difference between two means with t-test\n",
      "t is 1.5014 and has 47 degrees of freedom\n",
      "Probability of |t| exceeding 1.5014 is 0.1400\n",
      " \n",
      "1/F is 0.7192 and F is 1.3905\n",
      "Probability of F outside interval (0.7192,1.3905) is 0.4156\n"
     ]
    }
   ],
   "source": [
    "# eda12_07\n",
    "# generates simulated data for particle size example\n",
    "\n",
    "# Part A: first dataset\n",
    "\n",
    "# load data\n",
    "D=np.genfromtxt('testA.txt', delimiter='\\t');\n",
    "s = np.shape(D);\n",
    "NA = s[0];\n",
    "dobsA=np.zeros((NA,1));\n",
    "dobsA[:,0] = D;\n",
    "\n",
    "dtrueA = 100.0;     # known size of calibrations standard'\n",
    "dbartrueA = dtrueA; # known mean of measurement\n",
    "sd2trueA=1;         # known variance of measurement\n",
    "sdtrueA=1;          # known root-variance of measurement\n",
    "\n",
    "# calculate useful statistics\n",
    "print('file testA.txt');\n",
    "\n",
    "# number of observations\n",
    "print(\"N is %d\" % (NA) );\n",
    "\n",
    "dbartrueA = dtrueA;\n",
    "print(\"true mean %.4f\" % (dbartrueA) );\n",
    "\n",
    "dbarestA = np.sum(dobsA)/NA;\n",
    "print(\"estimated mean %.4f\" % (dbarestA) );\n",
    "\n",
    "print(\"true variance %.4f\" % (sd2trueA) );\n",
    "print(\"true root-variance %.4f\" % (sdtrueA) );\n",
    "\n",
    "sd2estA = np.matmul( (dobsA-dbartrueA).T, (dobsA-dbartrueA)) / (NA);\n",
    "print(\"estimated variance (given true mean) %.4f\" % sd2estA );\n",
    "\n",
    "sdestA = sqrt(sd2estA);\n",
    "print(\"estimated root-variance (given true mean) %.4f\" % (sdestA) );\n",
    "\n",
    "sd2estpA= np.matmul( (dobsA-dbarestA).T, (dobsA-dbarestA)) / (NA-1);\n",
    "print(\"estimated variance (given estimated mean) %.4f\" % (sd2estpA) );\n",
    "\n",
    "sdestpA=sqrt(sd2estpA);\n",
    "print(\"estimated root-variance (given estimated mean) %.4f\" % (sdestpA));\n",
    "\n",
    "ZA = abs(dbarestA-dbartrueA)/(sdtrueA/sqrt(NA));\n",
    "print(\"Z is %.4f\" % (ZA) );\n",
    "\n",
    "PZA = 1 - (stats.norm.cdf(ZA,0.0,1.0)-stats.norm.cdf(-ZA,0.0,1.0));\n",
    "print(\"Probability of |Z| exceeding %.4f is %.4f\" % (ZA, PZA) );\n",
    "\n",
    "chi2A = np.matmul( ((dobsA-dtrueA).T/sdtrueA), ((dobsA-dtrueA)/sdtrueA) );\n",
    "print(\"chi-squared is %.4f\" % (chi2A) );\n",
    "\n",
    "Pchi2A = 1 - stats.chi2.cdf(chi2A,NA);\n",
    "print(\"Probability of chi2 exceeding %.4f is %.4f\" % (chi2A, Pchi2A) );\n",
    "\n",
    "tA = abs(dbarestA-dbartrueA)/(sdestA/sqrt(NA));\n",
    "print(\"t is %.4f\" % (tA) );\n",
    "\n",
    "PtA = 1.0 - ( stats.t.cdf(tA,NA) - stats.t.cdf(-tA,NA) );\n",
    "print(\"Probability of |t| exceeding %.4f is %.4f\" % (tA, PtA) );\n",
    "\n",
    "print(' ');\n",
    "\n",
    "# Part B: second dataset\n",
    "\n",
    "# load data\n",
    "D=np.genfromtxt('testB.txt', delimiter='\\t');\n",
    "s = np.shape(D);\n",
    "NB = s[0];\n",
    "dobsB=np.zeros((NB,1));\n",
    "dobsB[:,0] = D;\n",
    "\n",
    "dtrueB = 100;             # known size of calibrations standard'\n",
    "dbartrueB = dtrueB;       # known mean\n",
    "sd2trueB = 1.0;           # known variance of measurement\n",
    "sdtrueB = sqrt(sd2trueB); # known root-variance of measurement\n",
    "\n",
    "# calculate useful statistics\n",
    "print('file testB.txt');\n",
    "\n",
    "# number of observations\n",
    "print(\"N is %d\" % (NB) );\n",
    "\n",
    "print(\"true mean %.4f\" % (dbartrueB) );\n",
    "\n",
    "dbarestB = np.sum(dobsB)/NB;\n",
    "print(\"estimated mean %.4f\" % (dbarestB) );\n",
    "\n",
    "print(\"true variance %.4f\" % (sd2trueB) );\n",
    "print(\"true root-variance %.4f\" % (sdtrueB) );\n",
    "\n",
    "sd2estB = np.matmul( (dobsB-dbartrueB).T, (dobsB-dbartrueB)) / (NB);\n",
    "print(\"estimated variance (given true mean) %.4f\" % sd2estB );\n",
    "\n",
    "sdestB = sqrt(sd2estB);\n",
    "print(\"estimated root-variance (given true mean) %.4f\" % (sdestB) );\n",
    "\n",
    "sd2estpB= np.matmul( (dobsB-dbarestB).T, (dobsB-dbarestB)) / (NB-1);\n",
    "print(\"estimated variance (given estimated mean) %.4f\" % (sd2estpB) );\n",
    "\n",
    "sdestpB=sqrt(sd2estpB);\n",
    "print(\"estimated root-variance (given estimated mean) %.4f\" % (sdestpB));\n",
    "\n",
    "ZB = abs(dbarestB-dbartrueB)/(sdtrueB/sqrt(NB));\n",
    "print(\"Z is %.4f\" % (ZB) );\n",
    "\n",
    "PZB = 1 - (stats.norm.cdf(ZB,0.0,1.0)-stats.norm.cdf(-ZB,0.0,1.0));\n",
    "print(\"Probability of |Z| exceeding %.4f is %.4f\" % (ZB, PZB) );\n",
    "\n",
    "chi2B = np.matmul( ((dobsB-dtrueB).T/sdtrueB), ((dobsB-dtrueB)/sdtrueB) );\n",
    "print(\"chi-squared is %.4f\" % (chi2B) );\n",
    "\n",
    "Pchi2B = 1 - stats.chi2.cdf(chi2B,NB);\n",
    "print(\"Probability of chi2 exceeding %.4f is %.4f\" % (chi2B, Pchi2B) );\n",
    "\n",
    "tB = abs(dbarestB-dbartrueB)/(sdestB/sqrt(NB));\n",
    "print(\"t is %.4f\" % (tB) );\n",
    "\n",
    "PtB = 1.0 - ( stats.t.cdf(tB,NB) - stats.t.cdf(-tB,NB) );\n",
    "print(\"Probability of |t| exceeding %.4f is %.4f\" % (tB, PtB) );\n",
    "\n",
    "print(' ');\n",
    "\n",
    "# PART 3 comparison of the tw datsets\n",
    "\n",
    "# are the means of the two tests significantly different ?\n",
    "\n",
    "print(\"significance of difference between two means, using Z-test\");\n",
    "A = dbarestA - dbarestB;\n",
    "B = sqrt(sd2trueA/NA+sd2estpB/NB);\n",
    "Z = abs(A/B);\n",
    "print(\"Z is %.4f\" % (Z) );\n",
    "\n",
    "# are the means of the two tests significantly different ?\n",
    "PZ = 1 - (stats.norm.cdf(Z,0,1)-stats.norm.cdf(-Z,0,1));\n",
    "print(\"Probability of |Z| exceeding %.4f is %.4f\" % (Z, PZ) );\n",
    "print(' ');\n",
    "\n",
    "print(\"difference between two means with t-test\");\n",
    "A = dbarestA - dbarestB;\n",
    "B = sqrt(sd2estpA/NA+sd2estpB/NB);\n",
    "t = abs(A/B);\n",
    "C = (sd2estpA/NA+sd2estpB/NB)**2;\n",
    "D = ((sd2estpA/NA)**2)/(NA-1) + ((sd2estpB/NB)**2)/(NB-1);\n",
    "M  =floor(C/D+0.5);\n",
    "print(\"t is %.4f and has %d degrees of freedom\" % (t, M) );\n",
    "\n",
    "# are the means of the two tests significantly different ?\n",
    "Pt = 1 - (stats.t.cdf(t,M)-stats.t.cdf(-t,M));\n",
    "print(\"Probability of |t| exceeding %.4f is %.4f\" % (t, Pt) );\n",
    "print(' ');\n",
    "\n",
    "# are the variances of the two tests significantly different ?\n",
    "F = (chi2A/NA) / (chi2B/NB);\n",
    "if( F<1 ):\n",
    "    F=1/F;\n",
    "print(\"1/F is %.4f and F is %.4f\" % (1/F, F) );\n",
    "\n",
    "PF = 1 - (stats.f.cdf(F,NA,NB)-stats.f.cdf(1/F,NA,NB));\n",
    "print(\"Probability of F outside interval (%.4f,%.4f) is %.4f\" % (1/F, F, PF) );\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 26020 characters]\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 26580 characters]\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "linear error 31.52, degrees of freedom 49\n",
      "cubic  error 27.18, degrees of freedom 47\n",
      "improvement in error 0.14\n",
      "1/F 0.8991 F 1.1122\n",
      "P(F<0.8991 or F>1.1122) = 0.714073\n",
      "Null hypothesis: fits the same up to stastical flucuation\n",
      "Null hypothesi cannot be excluded to 95% confidence\n"
     ]
    }
   ],
   "source": [
    "# eda12_08\n",
    "# illustration of F-test to assess difference in fit of two models\n",
    "# using a small fragment of Black Rock Forest temperature data\n",
    "\n",
    "D=np.genfromtxt('brf_fragment.txt', delimiter='\\t');\n",
    "s = np.shape(D);\n",
    "N = s[0];\n",
    "t = np.zeros((N,1));\n",
    "t[:,0] = 10+np.linspace(0,N-1,N);  # hours from start of orignal record\n",
    "dobs = np.zeros((N,1));\n",
    "dobs[:,0] = D;\n",
    "\n",
    "fig1 = plt.figure(1,figsize=(10,6));\n",
    "plt.plot(t,dobs,'k-');\n",
    "plt.title('linear fit');\n",
    "plt.xlabel('time t, hours');\n",
    "plt.ylabel('d(i)');\n",
    "\n",
    "# fit 1, straight line\n",
    "M=2;\n",
    "G=np.zeros((N,M));\n",
    "G[:,0]=np.ones((N,1)).ravel();\n",
    "G[:,1]=t.ravel();\n",
    "mestA=la.solve( np.matmul(G.T,G), np.matmul(G.T, dobs) );\n",
    "dpreA = np.matmul( G, mestA );\n",
    "plt.plot(t,dpreA,'k--');\n",
    "EA = np.matmul( (dobs-dpreA).T, (dobs-dpreA) );\n",
    "vA = N-M;\n",
    "plt.show();\n",
    "\n",
    "fig2 = plt.figure(2,figsize=(10,6));\n",
    "plt.plot(t,dobs,'k-');\n",
    "plt.title('cubic fit');\n",
    "plt.xlabel('time t, hours');\n",
    "plt.ylabel('d(i)');\n",
    "\n",
    "# fit 2, cubic\n",
    "M=4;\n",
    "G=np.zeros((N,M));\n",
    "G[:,0]=np.ones((N,1)).ravel();\n",
    "G[:,1]=t.ravel();\n",
    "G[:,2]=np.power(t,2).ravel();\n",
    "G[:,3]=np.power(t,3).ravel();\n",
    "mestB=la.solve( np.matmul(G.T,G), np.matmul(G.T, dobs) );\n",
    "dpreB = np.matmul( G, mestB );\n",
    "plt.plot(t,dpreB,'k--');\n",
    "EB = np.matmul( (dobs-dpreB).T, (dobs-dpreB) );\n",
    "vB = N-M;\n",
    "plt.show();\n",
    "\n",
    "print(\"linear error %.2f, degrees of freedom %d\" % (EA, vA) );\n",
    "print(\"cubic  error %.2f, degrees of freedom %d\" % (EB, vB) );\n",
    "\n",
    "\n",
    "print(\"improvement in error %.2f\" % ((EA-EB)/EA) );\n",
    "F = (EA/vA) / (EB/vB);\n",
    "print(\"1/F %.4f F %.4f\" % (1/F,F) );\n",
    "\n",
    "if( F<1 ):\n",
    "    F=1/F\n",
    "\n",
    "P = 1 - (stats.f.cdf(F,vA,vB)-stats.f.cdf(1/F,vA,vB));\n",
    "print(\"P(F<%.4f or F>%.4f) = %f\" % (1/F, F, P) );\n",
    "\n",
    "print('Null hypothesis: fits the same up to stastical flucuation')\n",
    "if( P>0.05):\n",
    "    print('Null hypothesi cannot be excluded to 95% confidence')\n",
    "else:\n",
    "    print('Null hypothesi can be excluded to 95% confidence')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "true      data variance 64.0000\n",
      "estimated data variance 61.4060\n",
      "variance of data estimated from PSD 61.45\n",
      " \n",
      "power spectral density\n",
      "true mean 3.1938 variance 10.2001\n",
      "estimated mean 3.0664 variance 8.6467\n",
      " \n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 54292 characters]\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 65144 characters]\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9500 confidence level is at s2 = 9.5677\n",
      "95.9064 percent of points are below 0.9500 confidence\n",
      " \n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 40332 characters]\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda12_09\n",
    "# illustration of variance of power spectral density\n",
    "# case of untapered random timeseries\n",
    "\n",
    "# generic time set up\n",
    "N=floor(2**10);\n",
    "Dt=0.025;\n",
    "T = N*Dt;\n",
    "t = np.zeros((N,1));\n",
    "t[:,0] = Dt*np.linspace(0,N-1,N); \n",
    "\n",
    "# create an uncorrelated random timeseries\n",
    "sd2true=64.0;\n",
    "sdtrue=sqrt(sd2true);\n",
    "draw = np.random.normal(0,sdtrue,(N,1));\n",
    "\n",
    "# generic frequency set up\n",
    "No2 = floor(N/2);\n",
    "fmax=1/(2.0*Dt);\n",
    "Df=fmax/No2;\n",
    "Nf=No2+1;\n",
    "f=np.zeros((N,1));\n",
    "f[:,0] = Df * np.concatenate( ( np.linspace(0,No2,Nf), np.linspace(-No2+1,-1,No2-1) ), axis=0 );\n",
    "Dw=2*pi*Df;\n",
    "w=np.zeros((N,1));\n",
    "w[:,0] = Dw * np.concatenate( ( np.linspace(0,No2,Nf), np.linspace(-No2+1,-1,No2-1) ), axis=0 );\n",
    "Nw=Nf;\n",
    "fpos=np.zeros((Nf,1));\n",
    "fpos[:,0] = Df * np.linspace(0,No2,Nf);\n",
    "wpos=np.zeros((Nf,1));\n",
    "wpos[:,0] = Dw * np.linspace(0,No2,Nf);              \n",
    "                    \n",
    "# constant window function\n",
    "W = np.ones((N,1));\n",
    "d = np.multiply( W, draw );\n",
    "\n",
    "# window changes variance of signal by ff\n",
    "ff = np.sum(np.multiply(W,W))/N;\n",
    "\n",
    "# true and estimated variance of d\n",
    "sd2est=np.std(d)**2;\n",
    "print(\"true      data variance %.4f\" % (ff*sd2true) );\n",
    "print(\"estimated data variance %.4f\" % (sd2est));\n",
    "\n",
    "# compute power spectral density\n",
    "fftraw = np.zeros((N,1), dtype=np.complex);\n",
    "fftraw = np.fft.fft(d,axis=0);                # raw DFT output\n",
    "dbar = np.zeros((Nf,1), dtype=np.complex);\n",
    "dbar[:,0] = Dt*fftraw[0:Nf,0];                # Fourier transform, + freqs only\n",
    "s2 = (2/T) * np.power(np.abs(dbar),2);      # power spectral density\n",
    "\n",
    "# statistics of power spectral density\n",
    "s2meanest = np.mean(s2);\n",
    "s2varest=np.std(s2)**2;\n",
    "s2sigmaest=sqrt(s2varest);\n",
    "sd2psd=Nf*Df*s2meanest;\n",
    "print(\"variance of data estimated from PSD %.2f\" % (sd2psd));\n",
    "print(' ');\n",
    "\n",
    "# true mean and variance\n",
    "p=2; # degrees of freedom\n",
    "c = ( ff*sd2true ) / (2*Nf*Df);  # normalization factor\n",
    "s2meantrue = p*c;\n",
    "s2vartrue = 2*p*(c**2);\n",
    "print(\"power spectral density\");\n",
    "print(\"true mean %.4f variance %.4f\" % (s2meantrue,s2vartrue));\n",
    "print(\"estimated mean %.4f variance %.4f\" % (s2meanest,s2varest));\n",
    "print(' ');\n",
    "\n",
    "# plot data\n",
    "fig1 = plt.figure(1,figsize=(10,6));\n",
    "damax = np.amax(d);\n",
    "plt.axis([0, T, -1.5*damax, 1.5*damax])\n",
    "plt.plot(t,d,'k-');\n",
    "plt.xlabel('time t (s)');\n",
    "plt.ylabel('d(t)');\n",
    "plt.title('random time series');\n",
    "plt.plot([0, T], [-2*sdtrue*sqrt(ff), -2*sdtrue*sqrt(ff)], 'k-' );\n",
    "plt.plot([0, T], [0, 0], 'k-', 'LineWidth', 1 );\n",
    "plt.plot([0, T], [2*sdtrue*sqrt(ff), 2*sdtrue*sqrt(ff)], 'k-' );\n",
    "\n",
    "# plot PSD\n",
    "fig2 = plt.figure(2,figsize=(10,6));\n",
    "plt.xlabel('frequency f (Hz)');\n",
    "plt.ylabel('ss(f)');\n",
    "plt.title('power spectral density');\n",
    "plt.plot(fpos,s2,'k-');\n",
    "\n",
    "# plot mean on spectra\n",
    "plt.plot([0, fmax], [s2meantrue, s2meantrue],'k-');\n",
    "\n",
    "# histogram of values\n",
    "Nhist=100\n",
    "s2min = np.min(s2);\n",
    "s2max = np.max(s2);\n",
    "ct, ed = np.histogram(s2,Nhist,(s2min,s2max));\n",
    "Nc = len(ct);\n",
    "Ne = len(ed);\n",
    "s2hist = np.zeros((Nc,1));\n",
    "s2hist[:,0]=ct;\n",
    "centers = np.zeros((Ne-1,1));\n",
    "centers[:,0]=0.5*(ed[1:Ne]+ed[0:Ne-1]);\n",
    "\n",
    "# chi-squared pdf\n",
    "Db=(centers[1,0]-centers[0,0])/c;\n",
    "s2histtrue=Nf*Db*stats.chi2.pdf(centers/c,p);\n",
    "\n",
    "# plot 95% confidence level on spectra\n",
    "pv=0.95;\n",
    "cv = stats.chi2.ppf(pv,p);\n",
    "plt.plot([fpos[0,0], fpos[Nf-1,0]], [cv*c, cv*c],'k-');\n",
    "plt.show();\n",
    "print(\"%.4f confidence level is at s2 = %.4f\" % (pv,cv*c) );\n",
    "percent = 100 * np.sum( s2hist[centers<=cv*c]) / np.sum(s2hist);\n",
    "print(\"%.4f percent of points are below %.4f confidence\" % (percent,pv));\n",
    "print(' ');\n",
    "\n",
    "fig3 = plt.figure(3,figsize=(10,6));\n",
    "plt.axis([0, np.max(centers), 0, max(s2hist) ])\n",
    "plt.xlabel('s2)');\n",
    "plt.ylabel('p(s2)');\n",
    "plt.title('histogram of power spectral density');\n",
    "plt.plot(centers,s2hist,'k-');\n",
    "plt.plot(centers,s2histtrue,'k-');\n",
    "plt.plot([cv*c, cv*c], [0, 0.1*max(s2hist)],'k:');\n",
    "plt.show();\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "true      data variance 25.4092\n",
      "estimated data variance 26.6966\n",
      "variance of data estimated from PSD 26.72\n",
      " \n",
      "power spectral density\n",
      "true mean 1.2680 variance 1.6078\n",
      "estimated mean 1.3335 variance 1.5545\n",
      " \n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 45968 characters]\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 61088 characters]\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9500 confidence level is at s2 = 3.7985\n",
      "93.9571 percent of points are below 0.9500 confidence\n",
      " \n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 40776 characters]\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda12_10\n",
    "# illustration of variance of power spectral density\n",
    "# case of random timeseries tapered with Hamming window\n",
    "\n",
    "# generic time set up\n",
    "N=floor(2**10);\n",
    "Dt=0.025;\n",
    "T = N*Dt;\n",
    "t = np.zeros((N,1));\n",
    "t[:,0] = Dt*np.linspace(0,N-1,N); \n",
    "\n",
    "# create an uncorrelated random timeseries\n",
    "sd2true=64.0;\n",
    "sdtrue=sqrt(sd2true);\n",
    "draw = np.random.normal(0,sdtrue,(N,1));\n",
    "\n",
    "# generic frequency set up\n",
    "No2 = floor(N/2);\n",
    "fmax=1/(2.0*Dt);\n",
    "Df=fmax/No2;\n",
    "Nf=No2+1;\n",
    "f=np.zeros((N,1));\n",
    "f[:,0] = Df * np.concatenate( ( np.linspace(0,No2,Nf), np.linspace(-No2+1,-1,No2-1) ), axis=0 );\n",
    "Dw=2*pi*Df;\n",
    "w=np.zeros((N,1));\n",
    "w[:,0] = Dw * np.concatenate( ( np.linspace(0,No2,Nf), np.linspace(-No2+1,-1,No2-1) ), axis=0 );\n",
    "Nw=Nf;\n",
    "fpos=np.zeros((Nf,1));\n",
    "fpos[:,0] = Df * np.linspace(0,No2,Nf);\n",
    "wpos=np.zeros((Nf,1));\n",
    "wpos[:,0] = Dw * np.linspace(0,No2,Nf);                              \n",
    "\n",
    "#  Hamming window function\n",
    "W = np.zeros((N,1));\n",
    "W[:,0] = 0.54 - 0.46*np.cos(2*pi*np.linspace(0,N-1,N)/(N-1));\n",
    "d = np.multiply( W, draw );\n",
    "\n",
    "# window changes variance of signal by ff\n",
    "ff = np.sum(np.multiply(W,W))/N;\n",
    "\n",
    "# true and estimated variance of d\n",
    "sd2est=np.std(d)**2;\n",
    "print(\"true      data variance %.4f\" % (ff*sd2true) );\n",
    "print(\"estimated data variance %.4f\" % (sd2est));\n",
    "\n",
    "# compute power spectral density\n",
    "fftraw = np.zeros((N,1), dtype=np.complex);\n",
    "fftraw = np.fft.fft(d,axis=0);                # raw DFT output\n",
    "dbar = np.zeros((Nf,1), dtype=np.complex);\n",
    "dbar[:,0] = Dt*fftraw[0:Nf,0];                # Fourier transform, + freqs only\n",
    "s2 = (2/T) * np.power(np.abs(dbar),2);      # power spectral density\n",
    "\n",
    "# statistics of power spectral density\n",
    "s2meanest = np.mean(s2);\n",
    "s2varest=np.std(s2)**2;\n",
    "s2sigmaest=sqrt(s2varest);\n",
    "sd2psd=Nf*Df*s2meanest;\n",
    "print(\"variance of data estimated from PSD %.2f\" % (sd2psd));\n",
    "print(' ');\n",
    "\n",
    "# true mean and variance\n",
    "p=2; # degrees of freedom\n",
    "c = ( ff*sd2true ) / (2*Nf*Df);  # normalization factor\n",
    "s2meantrue = p*c;\n",
    "s2vartrue = 2*p*(c**2);\n",
    "print(\"power spectral density\");\n",
    "print(\"true mean %.4f variance %.4f\" % (s2meantrue,s2vartrue));\n",
    "print(\"estimated mean %.4f variance %.4f\" % (s2meanest,s2varest));\n",
    "print(' ');\n",
    "\n",
    "# plot data\n",
    "fig1 = plt.figure(1,figsize=(10,6));\n",
    "damax = np.amax(d);\n",
    "plt.axis([0, T, -1.5*damax, 1.5*damax])\n",
    "plt.plot(t,d,'k-');\n",
    "plt.xlabel('time t (s)');\n",
    "plt.ylabel('d(t)');\n",
    "plt.title('random time series');\n",
    "plt.plot([0, T], [-2*sdtrue*sqrt(ff), -2*sdtrue*sqrt(ff)], 'k-' );\n",
    "plt.plot([0, T], [0, 0], 'k-', 'LineWidth', 1 );\n",
    "plt.plot([0, T], [2*sdtrue*sqrt(ff), 2*sdtrue*sqrt(ff)], 'k-' );\n",
    "\n",
    "# plot PSD\n",
    "fig2 = plt.figure(2,figsize=(10,6));\n",
    "plt.xlabel('frequency f (Hz)');\n",
    "plt.ylabel('ss(f)');\n",
    "plt.title('power spectral density');\n",
    "plt.plot(fpos,s2,'k-');\n",
    "\n",
    "# plot mean on spectra\n",
    "plt.plot([0, fmax], [s2meantrue, s2meantrue],'k-');\n",
    "\n",
    "# histogram of values\n",
    "Nhist=100\n",
    "s2min = np.min(s2);\n",
    "s2max = np.max(s2);\n",
    "ct, ed = np.histogram(s2,Nhist,(s2min,s2max));\n",
    "Nc = len(ct);\n",
    "Ne = len(ed);\n",
    "s2hist = np.zeros((Nc,1));\n",
    "s2hist[:,0]=ct;\n",
    "centers = np.zeros((Ne-1,1));\n",
    "centers[:,0]=0.5*(ed[1:Ne]+ed[0:Ne-1]);\n",
    "\n",
    "# chi-squared pdf\n",
    "Db=(centers[1,0]-centers[0,0])/c;\n",
    "s2histtrue=Nf*Db*stats.chi2.pdf(centers/c,p);\n",
    "\n",
    "# plot 95% confidence level on spectra\n",
    "pv=0.95;\n",
    "cv = stats.chi2.ppf(pv,p);\n",
    "plt.plot([fpos[0,0], fpos[Nf-1,0]], [cv*c, cv*c],'k-');\n",
    "plt.show();\n",
    "print(\"%.4f confidence level is at s2 = %.4f\" % (pv,cv*c) );\n",
    "percent = 100 * np.sum( s2hist[centers<=cv*c]) / np.sum(s2hist);\n",
    "print(\"%.4f percent of points are below %.4f confidence\" % (percent,pv));\n",
    "print(' ');\n",
    "\n",
    "fig3 = plt.figure(3,figsize=(10,6));\n",
    "plt.axis([0, np.max(centers), 0, max(s2hist) ])\n",
    "plt.xlabel('s2)');\n",
    "plt.ylabel('p(s2)');\n",
    "plt.title('histogram of power spectral density');\n",
    "plt.plot(centers,s2hist,'k-');\n",
    "plt.plot(centers,s2histtrue,'k-');\n",
    "plt.plot([cv*c, cv*c], [0, 0.1*max(s2hist)],'k:');\n",
    "plt.show();\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "true      data variance 25.4092\n",
      "estimated data variance 25.8134\n",
      "variance of data estimated from PSD 25.92\n",
      " \n",
      "power spectral density\n",
      "true mean 1.2680 variance 1.6078\n",
      "estimated mean 1.2934 variance 1.6670\n",
      " \n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 47252 characters]\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 51276 characters]\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9500 confidence level is at s2 = 3.7985\n",
      "95.7115 percent of points are below 0.9500 confidence\n",
      " \n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 29132 characters]\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda12_11\n",
    "# illustration of variance of power spectral density\n",
    "# case of low amplitude cosine signal buried in random noise, tapered with Hamming window\n",
    "\n",
    "# generic time setup\n",
    "N=floor(2**10);\n",
    "Dt=0.025;\n",
    "T = N*Dt;\n",
    "t = np.zeros((N,1));\n",
    "t[:,0]=Dt*np.linspace(0,N-1,N); \n",
    "\n",
    "# generic time/frequency set up\n",
    "No2 = floor(N/2);\n",
    "fmax=1/(2.0*Dt);\n",
    "Df=fmax/No2;\n",
    "Nf=No2+1;\n",
    "f=np.zeros((N,1));\n",
    "f[:,0] = Df * np.concatenate( ( np.linspace(0,No2,Nf), np.linspace(-No2+1,-1,No2-1) ), axis=0 );\n",
    "Dw=2*pi*Df;\n",
    "w=np.zeros((N,1));\n",
    "w[:,0] = Dw * np.concatenate( ( np.linspace(0,No2,Nf), np.linspace(-No2+1,-1,No2-1) ), axis=0 );\n",
    "Nw=Nf;\n",
    "fpos=np.zeros((Nf,1));\n",
    "fpos[:,0] = Df * np.linspace(0,No2,Nf);\n",
    "wpos=np.zeros((Nf,1));\n",
    "wpos[:,0] = Dw * np.linspace(0,No2,Nf);     \n",
    "\n",
    "# low amplitude periodic signal buried in the noise\n",
    "sd2true=64.0;\n",
    "sdtrue=sqrt(sd2true);\n",
    "draw = 0.25*sdtrue*np.cos(2*pi*(fmax/4)*t) + np.random.normal(0,sdtrue,(N,1));                \n",
    "\n",
    "#  Hamming window function\n",
    "W = np.zeros((N,1));\n",
    "W[:,0] = 0.54 - 0.46*np.cos(2*pi*np.linspace(0,N-1,N)/(N-1));\n",
    "d = np.multiply( W, draw );\n",
    "\n",
    "# window changes variance of signal by ff\n",
    "ff = np.sum(np.multiply(W,W))/N;\n",
    "\n",
    "# true and estimated variance of d\n",
    "sd2est=np.std(d)**2;\n",
    "print(\"true      data variance %.4f\" % (ff*sd2true) );\n",
    "print(\"estimated data variance %.4f\" % (sd2est));\n",
    "\n",
    "# compute power spectral density\n",
    "fftraw = np.zeros((N,1), dtype=np.complex);\n",
    "fftraw = np.fft.fft(d,axis=0);                # raw DFT output\n",
    "dbar = np.zeros((Nf,1), dtype=np.complex);\n",
    "dbar[:,0] = Dt*fftraw[0:Nf,0];                # Fourier transform, + freqs only\n",
    "s2 = (2/T) * np.power(np.abs(dbar),2);      # power spectral density\n",
    "\n",
    "# statistics of power spectral density\n",
    "s2meanest = np.mean(s2);\n",
    "s2varest=np.std(s2)**2;\n",
    "s2sigmaest=sqrt(s2varest);\n",
    "sd2psd=Nf*Df*s2meanest;\n",
    "print(\"variance of data estimated from PSD %.2f\" % (sd2psd));\n",
    "print(' ');\n",
    "\n",
    "# true mean and variance\n",
    "p=2; # degrees of freedom\n",
    "c = ( ff*sd2true ) / (2*Nf*Df);  # normalization factor\n",
    "s2meantrue = p*c;\n",
    "s2vartrue = 2*p*(c**2);\n",
    "print(\"power spectral density\");\n",
    "print(\"true mean %.4f variance %.4f\" % (s2meantrue,s2vartrue));\n",
    "print(\"estimated mean %.4f variance %.4f\" % (s2meanest,s2varest));\n",
    "print(' ');\n",
    "\n",
    "# plot data\n",
    "fig1 = plt.figure(1,figsize=(10,6));\n",
    "damax = np.amax(d);\n",
    "plt.axis([0, T, -1.5*damax, 1.5*damax])\n",
    "plt.plot(t,d,'k-');\n",
    "plt.xlabel('time t (s)');\n",
    "plt.ylabel('d(t)');\n",
    "plt.title('random time series');\n",
    "plt.plot([0, T], [-2*sdtrue*sqrt(ff), -2*sdtrue*sqrt(ff)], 'k-' );\n",
    "plt.plot([0, T], [0, 0], 'k-', 'LineWidth', 1 );\n",
    "plt.plot([0, T], [2*sdtrue*sqrt(ff), 2*sdtrue*sqrt(ff)], 'k-' );\n",
    "\n",
    "# plot PSD\n",
    "fig2 = plt.figure(2,figsize=(10,6));\n",
    "plt.xlabel('frequency f (Hz)');\n",
    "plt.ylabel('ss(f)');\n",
    "plt.title('power spectral density');\n",
    "plt.plot(fpos,s2,'k-');\n",
    "\n",
    "# plot mean on spectra\n",
    "plt.plot([0, fmax], [s2meantrue, s2meantrue],'k-');\n",
    "\n",
    "# histogram of values\n",
    "Nhist=100\n",
    "s2min = np.min(s2);\n",
    "s2max = np.max(s2);\n",
    "ct, ed = np.histogram(s2,Nhist,(s2min,s2max));\n",
    "Nc = len(ct);\n",
    "Ne = len(ed);\n",
    "s2hist = np.zeros((Nc,1));\n",
    "s2hist[:,0]=ct;\n",
    "centers = np.zeros((Ne-1,1));\n",
    "centers[:,0]=0.5*(ed[1:Ne]+ed[0:Ne-1]);\n",
    "\n",
    "# chi-squared pdf\n",
    "Db=(centers[1,0]-centers[0,0])/c;\n",
    "s2histtrue=Nf*Db*stats.chi2.pdf(centers/c,p);\n",
    "\n",
    "# plot 95% confidence level on spectra\n",
    "pv=0.95;\n",
    "cv = stats.chi2.ppf(pv,p);\n",
    "plt.plot([fpos[0,0], fpos[Nf-1,0]], [cv*c, cv*c],'k-');\n",
    "plt.show();\n",
    "print(\"%.4f confidence level is at s2 = %.4f\" % (pv,cv*c) );\n",
    "percent = 100 * np.sum( s2hist[centers<=cv*c]) / np.sum(s2hist);\n",
    "print(\"%.4f percent of points are below %.4f confidence\" % (percent,pv));\n",
    "print(' ');\n",
    "\n",
    "fig3 = plt.figure(3,figsize=(10,6));\n",
    "plt.axis([0, np.max(centers), 0, max(s2hist) ])\n",
    "plt.xlabel('s2)');\n",
    "plt.ylabel('p(s2)');\n",
    "plt.title('histogram of power spectral density');\n",
    "plt.plot(centers,s2hist,'k-');\n",
    "plt.plot(centers,s2histtrue,'k-');\n",
    "plt.plot([cv*c, cv*c], [0, 0.1*max(s2hist)],'k:');\n",
    "plt.show();\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 22932 characters]\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "standard method:   slope 0.5252 +/- 0.0150 (2 sigma)\n",
      "bootstrap method:  slope 0.5251 +/- 0.0173 (2 sigma)\n",
      "bootstrap method:  slope between 0.5087 and 0.5423 (95 percent confidence)\n"
     ]
    }
   ],
   "source": [
    "# eda12_12\n",
    "# illustration of Bootstrap statistics of slope in least\n",
    "# squares fit of a straight line to a small fragment of\n",
    "# Black Rock Forest temperature data\n",
    "\n",
    "# load data and set up time axis\n",
    "D=np.genfromtxt('brf_fragment.txt', delimiter='\\t');\n",
    "s = np.shape(D);\n",
    "N = s[0];\n",
    "torig = np.zeros((N,1));\n",
    "torig[:,0] = 10+np.linspace(0,N-1,N);  # hours from start of orignal record\n",
    "dorig = np.zeros((N,1));\n",
    "dorig[:,0] = D;\n",
    "\n",
    "# least-squares slope and its variance of original data\n",
    "M=2;\n",
    "G=np.zeros((N,M));\n",
    "G[:,0]=np.ones((N,));\n",
    "G[:,1]=torig.ravel();\n",
    "GTG = np.matmul(G.T,G);\n",
    "mest=  la.solve( GTG, np.matmul(G.T,dorig) );\n",
    "slopeorig=mest[1,0];\n",
    "    \n",
    "# usual error estimate based on propagagtion\n",
    "# of posterior estimate of sigma2 of data into\n",
    "# a sigma2 of slope\n",
    "e = dorig - np.matmul(G,mest);\n",
    "E = np.matmul(e.T,e);\n",
    "sigma2dorig = E/N;\n",
    "Cmorig = sigma2dorig * la.inv(GTG);\n",
    "sigma2morig = Cmorig[1,1];\n",
    "\n",
    "# number of realizations;\n",
    "Nr = 10000;\n",
    "slope = np.zeros((Nr,1));\n",
    "\n",
    "# loop over realizations\n",
    "for p in range(Nr):\n",
    "    \n",
    "    # resample\n",
    "    rowindex=np.random.randint(N,size=N);\n",
    "    t = np.zeros((N,1));\n",
    "    t[:,0] = torig[rowindex,0];\n",
    "    d = np.zeros((N,1));\n",
    "    d[:,0] = dorig[rowindex,0];\n",
    "    # straight line fit\n",
    "    M=2;\n",
    "    G=np.zeros((N,M));\n",
    "    G[:,0]=np.ones((N,));\n",
    "    G[:,1]=t.ravel();\n",
    "    GTG = np.matmul(G.T,G);\n",
    "    mest=  la.solve( GTG, np.matmul(G.T,d) );\n",
    "    slope[p,0]=mest[1,0];\n",
    "\n",
    "# compute histogram of slopes and normalize\n",
    "# normalize to a probability density distribution\n",
    "# with Dt*sum(pbootstrap)=1;\n",
    "\n",
    "Nhist=1000\n",
    "slopemin = 0.2;\n",
    "slopemax = 0.8;\n",
    "ct, ed = np.histogram(slope,Nhist,(slopemin,slopemax));\n",
    "Nc = len(ct);\n",
    "Ne = len(ed);\n",
    "slopehist = np.zeros((Nc,1));\n",
    "slopehist[:,0]=ct;\n",
    "centers = np.zeros((Ne-1,1));\n",
    "centers[:,0]=0.5*(ed[1:Ne]+ed[0:Ne-1]);\n",
    "\n",
    "# approximate pdf via normalized histogram\n",
    "Ds = centers[1,0] - centers[0,0];\n",
    "norm = Ds*np.sum(slopehist);\n",
    "pbootstrap =slopehist/norm;\n",
    "\n",
    "# compute Normal distribution of slopes\n",
    "# using parameters estimated from original data\n",
    "norm = (1/sqrt(2*pi*sigma2morig));\n",
    "pnormal = norm * np.exp(-np.power((centers-slopeorig),2)/(2*sigma2morig));\n",
    "\n",
    "# plot pdf\n",
    "fig1 = plt.figure(1,figsize=(10,6));\n",
    "pdfmax = np.max(pbootstrap);\n",
    "plt.axis([0.4, 0.6, 0, 1.5*pdfmax]);\n",
    "plt.plot(centers,pbootstrap,'k-');\n",
    "plt.plot(centers,pnormal,'k:');\n",
    "plt.xlabel('slope b');\n",
    "plt.ylabel('p(b)');\n",
    "\n",
    "# mean (expectation) and variance of bootstrap distibution\n",
    "mb = Ds*np.sum( np.multiply( centers, pbootstrap) );\n",
    "vb = Ds*np.sum( np.multiply( np.power(centers-mb,2) , pbootstrap) );\n",
    "\n",
    "# 5% and 95% confidence of bootstrap distibution\n",
    "Pbootstrap = np.zeros((Nc,1));\n",
    "Pbootstrap[:,0] = Ds*np.cumsum(pbootstrap);\n",
    "\n",
    "ilovec = np.where( (Pbootstrap>=0.025).ravel() );\n",
    "ilox = ilovec[0];\n",
    "ilo = ilox[0];\n",
    "blo = centers[ilo];\n",
    "\n",
    "ihivec = np.where(Pbootstrap>=0.975);\n",
    "ihix = ihivec[0];\n",
    "ihi = ihix[0]\n",
    "bhi = centers[ihi];\n",
    "\n",
    "plt.plot( [blo, blo], [0, 0.2*pdfmax], 'k-'  );\n",
    "plt.plot( [bhi, bhi], [0, 0.2*pdfmax], 'k-'  );\n",
    "plt.plot( [mb, mb], [0, 0.3*pdfmax], 'k-'  );\n",
    "plt.show();\n",
    "\n",
    "print(\"standard method:   slope %.4f +/- %.4f (2 sigma)\" % (slopeorig, 2*sqrt(sigma2morig)));\n",
    "print(\"bootstrap method:  slope %.4f +/- %.4f (2 sigma)\" % (mb, 2*sqrt(vb)))\n",
    "print(\"bootstrap method:  slope between %.4f and %.4f (95 percent confidence)\" % (blo, bhi));\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CaO/Na2O ratio of Factor 2\n",
      "mean 0.485729 variance 0.000193\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 19284 characters]\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "bootstrap method:  ratio 0.4857 +/- 0.0278 (2 sigma)\n",
      "bootstrap method:  ratio between 0.4586 and 0.5126 (95 percent confidence)\n"
     ]
    }
   ],
   "source": [
    "# eda12_13\n",
    "# Atlantic Rocks dataset\n",
    "# bootstrap error analysis for\n",
    "# CaO/Na2O ratio of varimax factor 2\n",
    "# note CaO/Na2O are elements 6/7\n",
    "\n",
    "# load data\n",
    "Draw = np.genfromtxt('rocks.txt', delimiter='\\t')\n",
    "N, M = np.shape(Draw);\n",
    "\n",
    "# number of realizations\n",
    "Nr=10000;\n",
    "ratio = np.zeros((Nr,1));\n",
    "\n",
    "for ir in range(Nr):\n",
    "\n",
    "    # resample\n",
    "    rowindex=np.random.randint(N,size=N);\n",
    "    D = np.zeros((N,M));\n",
    "    D[:,:] = Draw[rowindex,:];\n",
    "\n",
    "    # compute factors and factor loadings using singular value decompostion\n",
    "    [U, s, VT] = la.svd(D,full_matrices=False);\n",
    "    sh = np.shape(s);\n",
    "    Ns = sh[0];\n",
    "\n",
    "    # keep only first few singular values\n",
    "    P=5;\n",
    "    F = np.copy(VT[0:P,:]);\n",
    "    C = np.matmul( U[:,0:P], np.diag(s[0:P]));\n",
    "    \n",
    "    # spike these factors using the varimax procedure\n",
    "    k = [1, 2, 3, 4];\n",
    "    Nk = len(k);\n",
    "\n",
    "    # initial rotated factor matrix, FP, abd rotation matrix, MR\n",
    "    MR=np.identity(P);\n",
    "    FP=np.copy(F);\n",
    "    \n",
    "    fA = np.zeros((M,1));\n",
    "    fB = np.zeros((M,1));\n",
    "    fAp = np.zeros((M,1));\n",
    "    fBp = np.zeros((M,1));\n",
    "    mA = np.zeros((P,1));\n",
    "    mB = np.zeros((P,1));\n",
    "    mAp = np.zeros((P,1));\n",
    "    mBp = np.zeros((P,1));\n",
    "    for iter in range(3):\n",
    "        for ii in range(Nk):\n",
    "            for jj in range(ii+1,Nk):\n",
    "                # spike factors i and j\n",
    "                i=k[ii];\n",
    "                j=k[jj];\n",
    "                # copy factors from matrix to vectors\n",
    "                fA[:,0] = np.copy(FP[i,:]);\n",
    "                fB[:,0] = np.copy(FP[j,:]);\n",
    "                # standard varimax procedure to determine rotation angle q\n",
    "                u = np.power(fA,2) - np.power(fB,2);\n",
    "                v = 2.0 * np.multiply(fA, fB);\n",
    "                AA = 2*M*np.matmul(u.T,v);\n",
    "                BB = np.sum(u)*np.sum(v);\n",
    "                top = AA - BB;\n",
    "                CC = M*(np.matmul(u.T,u)-np.matmul(v.T,v));\n",
    "                DD = (np.sum(u)**2)-(np.sum(v)**2);\n",
    "                bot =  CC - DD;\n",
    "                q = 0.25 * atan2(top,bot);\n",
    "                # rotate factors\n",
    "                cq = cos(q);\n",
    "                sq = sin(q);\n",
    "                fAp =  cq*fA + sq*fB;\n",
    "                fBp = -sq*fA + cq*fB;\n",
    "                # put back into factor matrix, FP\n",
    "                FP[i,:] = fAp.ravel();\n",
    "                FP[j,:] = fBp.ravel();\n",
    "                # accumulate rotation\n",
    "                mA[:,0] = np.copy(MR[i,:]);\n",
    "                mB[:,0] = np.copy(MR[j,:]);\n",
    "                mAp =  cq*mA + sq*mB;\n",
    "                mBp = -sq*mA + cq*mB;\n",
    "                MR[i,:] = mAp.ravel();\n",
    "                MR[j,:] = mBp.ravel();\n",
    "\n",
    "    # new factor loadings\n",
    "    CP = np.matmul(C,MR.T);\n",
    "    \n",
    "    # MATLAB: ratio of elemnts 6 and 7 of factor 2 \n",
    "    # Python: ratio of elemnts 5 and 6 of factor 1\n",
    "    ratio[ir,0]=FP[1,5]/FP[1,6];\n",
    "\n",
    "Nhist=500;\n",
    "ratiomin = 0.2;\n",
    "ratiomax = 0.8;\n",
    "ct, ed = np.histogram(ratio,Nhist,(ratiomin,ratiomax));\n",
    "Nc = len(ct);\n",
    "Ne = len(ed);\n",
    "ratiohist = np.zeros((Nc,1));\n",
    "ratiohist[:,0]=ct;\n",
    "centers = np.zeros((Ne-1,1));\n",
    "centers[:,0]=0.5*(ed[1:Ne]+ed[0:Ne-1]);\n",
    "\n",
    "# approximate pdf via normalized histogram\n",
    "Dr = centers[1,0] - centers[0,0];\n",
    "norm = Dr*np.sum(ratiohist);\n",
    "pbootstrap = ratiohist/norm;\n",
    "\n",
    "# mean (expectation) and variance of bootstrap distibution\n",
    "mr = Dr*np.sum( np.multiply( centers, pbootstrap) );\n",
    "vr = Dr*np.sum( np.multiply( np.power(centers-mr,2) , pbootstrap) );\n",
    "print(\"CaO/Na2O ratio of Factor 2\");\n",
    "print(\"mean %.6f variance %.6f\" % (mr,vr) );\n",
    "\n",
    "# 5% and 95% confidence of bootstrap distibution\n",
    "Pbootstrap = np.zeros((Nc,1));\n",
    "Pbootstrap[:,0] = Dr*np.cumsum(pbootstrap);\n",
    "\n",
    "ilovec = np.where( (Pbootstrap>=0.025).ravel() );\n",
    "ilox = ilovec[0];\n",
    "ilo = ilox[0];\n",
    "blo = centers[ilo];\n",
    "\n",
    "ihivec = np.where(Pbootstrap>=0.975);\n",
    "ihix = ihivec[0];\n",
    "ihi = ihix[0]\n",
    "bhi = centers[ihi];\n",
    "\n",
    "# plot pdf\n",
    "fig1 = plt.figure(1,figsize=(10,6));\n",
    "pdfmax = np.max(pbootstrap);\n",
    "plt.axis([0.4, 0.6, 0, 1.5*pdfmax]);\n",
    "plt.plot(centers,pbootstrap,'k-');\n",
    "plt.xlabel('CaO/Na2O ratio, r');\n",
    "plt.ylabel('p(r)');\n",
    "\n",
    "plt.plot( [blo, blo], [0, 0.2*pdfmax], 'k-'  );\n",
    "plt.plot( [bhi, bhi], [0, 0.2*pdfmax], 'k-'  );\n",
    "plt.plot( [mr, mr], [0, 0.3*pdfmax], 'k-'  );\n",
    "plt.show();\n",
    "\n",
    "print(\"bootstrap method:  ratio %.4f +/- %.4f (2 sigma)\" % (mr, 2*sqrt(vr)))\n",
    "print(\"bootstrap method:  ratio between %.4f and %.4f (95 percent confidence)\" % (blo, bhi));\n",
    "\n",
    "\n"
   ]
  },
  {
   "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
