{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# edapy_12_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",
    "# function to make a numpy N-by-1 column vector\n",
    "# c=eda_cvec1(v) from entity v that is array-like,\n",
    "# including a number, a list of numbers, a tuple\n",
    "# of numbers, an N-by-0 np array and a N-by1 np array.\n",
    "# It promotes integers to floats, and integers\n",
    "# and floats to complex, by context. The function\n",
    "# also insures that, if v is an np array, that\n",
    "# c is a copy, as contrasted to a view, of it.\n",
    "# This version takes just one input argment.\n",
    "# whereas c=eda_cvec(v1,v2,...) concatenates\n",
    "# many argiments.\n",
    "def eda_cvec1(v):\n",
    "    if isinstance(v, int):\n",
    "        w = np.zeros((1,1),dtype=int);\n",
    "        w[0,0] = v;\n",
    "        return w;\n",
    "    elif isinstance(v, float):\n",
    "        w = np.zeros((1,1),dtype=float);\n",
    "        w[0,0] = v;\n",
    "        return w;\n",
    "    elif isinstance(v, complex):\n",
    "        w = np.zeros((1,1),dtype=complex);\n",
    "        w[0,0] = v;\n",
    "        return w;\n",
    "    elif isinstance(v, np.ndarray):\n",
    "        s = np.shape(v);\n",
    "        if len(s) == 1:\n",
    "            return np.copy(np.reshape(v,(s[0],1)));\n",
    "        else:\n",
    "            [r,c]=s;\n",
    "            if( c==1 ):\n",
    "                return(np.copy(v));\n",
    "            elif(r==1):\n",
    "                return(np.copy(v.T));\n",
    "            else:\n",
    "                raise TypeError(\"eda_cvec: %d by %d ndarray not allowed\" % (r, c));\n",
    "    elif isinstance(v, list):\n",
    "        r = len(v);\n",
    "        t = int;\n",
    "        for vi in v:\n",
    "            if isinstance(vi,int):\n",
    "                pass;\n",
    "            elif isinstance(vi,float):\n",
    "                t=float;\n",
    "            elif isinstance(vi,complex):\n",
    "                t=complex;\n",
    "                break;\n",
    "            else:\n",
    "                raise TypeError(\"eda_cvec: list contains unsupported type %s\" % type(vi));\n",
    "        w = np.zeros((r,1),dtype=t);\n",
    "        w[:,0] = v;\n",
    "        return w;\n",
    "    elif isinstance(v, tuple):\n",
    "        r = len(v);\n",
    "        t = int;\n",
    "        for vi in v:\n",
    "            if isinstance(vi,int):\n",
    "                pass;\n",
    "            elif isinstance(vi,float):\n",
    "                t=float;\n",
    "            elif isinstance(vi,complex):\n",
    "                t=complex;\n",
    "                break;\n",
    "            else:\n",
    "                raise TypeError(\"eda_cvec: tuple contains unsupported type %s\" % type(vi));\n",
    "        w = np.zeros((r,1),dtype=t);\n",
    "        w[:,0] = v;\n",
    "        return w;\n",
    "    else:\n",
    "        raise TypeError(\"eda_cvec: %s not supported\" % type(v));\n",
    "        \n",
    "# function to make a numpy N-by-1 column vector\n",
    "# c=eda_cvec(v1, v2, ...) from a list of several\n",
    "# array-like entities v1, v2, including a number\n",
    "# a list of numbers, a tuple of numbers, an N-by-0 np array\n",
    "# and a N-by-1 np array. The function\n",
    "# also insures that, if v is an np array, that\n",
    "# c is a copy, as contrasted to a view, of it\n",
    "# It promotes integers to floats, and integers\n",
    "# and floats to complex, by context.\n",
    "# This version concatenates many argments,\n",
    "# whereas c=eda_cvec1(v1) takes just one argiment.\n",
    "# I recommend always using eda_cvec(v1, v2, ...)\n",
    "def eda_cvec(*argv):\n",
    "    t = int;\n",
    "    Nt = 0;\n",
    "    for a in argv:\n",
    "        v = eda_cvec1(a);\n",
    "        N,M = np.shape(v);\n",
    "        Nt = Nt + N;\n",
    "        if( N==0 ):\n",
    "            continue; # skip vector of zero length\n",
    "        if (t==int) and isinstance(v[0,0],float):\n",
    "            t=float;\n",
    "        elif isinstance(v[0,0],complex):\n",
    "            t=complex;\n",
    "    w = np.zeros((Nt,1),dtype=t);\n",
    "    Nt = 0;\n",
    "    for a in argv:\n",
    "        v = eda_cvec1(a);\n",
    "        N,M = np.shape(v);\n",
    "        w[Nt:Nt+N,0] = v[0:N,0];\n",
    "        Nt = Nt + N;\n",
    "    return w;\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",
    "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": 4,
   "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": [
    "# edapy_12_01: p.d.f. of square of a random variable\n",
    "\n",
    "# standard set-up of e-axis\n",
    "L=101;       # length L\n",
    "emin=0.01;   # start value; avoid zero\n",
    "emax=3.00;   # end value\n",
    "De=(emax-emin)/(L-1); # sampling interva;\n",
    "e = eda_cvec( emin+De*np.linspace(0,L-1,L)); # vector e\n",
    "\n",
    "# Nornal p.d.f. for p(e), for e>=0\n",
    "pe = sqrt(2/pi)*np.exp(-0.5*np.power(e,2));\n",
    "\n",
    "Ae = De*np.sum(pe); # area under p.d.f.\n",
    "\n",
    "# stadard setup of E axis\n",
    "Emin=0.01;   # start value, avoid zero\n",
    "Emax=3;      # end value\n",
    "DE=(Emax-Emin)/(L-1);  # sampling interval\n",
    "E = eda_cvec(Emin+DE*np.linspace(0,L-1,L)); # vector E\n",
    "\n",
    "pE = np.sqrt(np.reciprocal(2*pi*E)) * np.exp(-0.5*E); # chi2 p.d.f.\n",
    "AE = DE*sum(pE); # area under p.d.f.\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": 5,
   "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": [
    "# edapy_12_02: chi-squared p.d.f.\n",
    "\n",
    "# standard set-up of X2-axis\n",
    "L=1001;                    # length\n",
    "X2min=0.0;                 # start value\n",
    "X2max=10.0;                # end value\n",
    "DX2=(X2max-X2min)/(L-1);   # sampling\n",
    "X2 = eda_cvec(np.linspace(X2min,X2max,L)); # vector X2\n",
    "\n",
    "# plot setup\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",
    "# plot p(X2,N) p.d.f for several values of N\n",
    "for N in range(1,6):  # loop over degrees of freedom N\n",
    "    pX2 = eda_cvec(stats.chi2.pdf(X2,N)); # chi-squared p.d.f.\n",
    "    plt.plot(X2,pX2,'k-');\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "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": [
    "# edapy_12_03: student t p.d.f.\n",
    "\n",
    "# standard set-up of t-axis\n",
    "L=1001;             # length\n",
    "tmin=-5.0;          # start value\n",
    "tmax=5.0;           # end value\n",
    "Dt=(tmax-tmin)/(L-1); # sampling interval\n",
    "t = eda_cvec(np.linspace(tmin,tmax,L));\n",
    "\n",
    "# plot setup\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",
    "# plot p(t,N) for several values of N\n",
    "pt = np.zeros((L,1))\n",
    "for N in range(1,6):\n",
    "    pt = eda_cvec(stats.t.pdf(t,N)); # t p.d.f.\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": [
    "# edapy_12_04: F p.d.f.\n",
    "\n",
    "# make a vector F\n",
    "L=1001;                # length\n",
    "Fmin=0.0;              # start value\n",
    "Fmax=5.0;              # end value\n",
    "DF=(Fmax-Fmin)/(L-1);  # sampling interval\n",
    "F=eda_cvec(np.linspace(Fmin,Fmax,L));\n",
    "\n",
    "# plot setup\n",
    "fig1 = plt.figure(1,figsize=(12,12));\n",
    "\n",
    "# plot of p(F,N,M) has several frames, N constant for each frame\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 = eda_cvec(stats.f.pdf(F,N,M));\n",
    "        plt.plot(F,pF,'k-');\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# edapy_12_05,  simulated particle size data 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 mytestA.txt and mytestB.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(\"../Data/mytestA.txt\",dobs1, delimiter=\"\\t\");\n",
    "\n",
    "# calibration test 2\n",
    "dobs2 = np.random.normal(d,sdtrue,(N,1));\n",
    "np.savetxt(\"../Data/mytestB.txt\",dobs2, delimiter=\"\\t\");\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "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": [
    "# edapy_12_06:, Z and chi-squared tests of particle size data\n",
    "\n",
    "# Z distribution, two-sided test\n",
    "Z=0.278;\n",
    "Pi = stats.norm.cdf(abs(Z),0.0,1.0)-stats.norm.cdf(-abs(Z),0.0,1.0);\n",
    "Po = 1.0-(stats.norm.cdf(abs(Z),0.0,1.0)-stats.norm.cdf(-abs(Z),0.0,1.0));\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\" % (Pi) );\n",
    "print(\"  The probability that Z is outside the interval %.4f to %.4f\" % (-Z, Z));\n",
    "print(\"    is %.4f\" % (Po) );\n",
    "\n",
    "# Chi-squared distribution, one-sided test\n",
    "X2=21.9;\n",
    "N=25;\n",
    "Pi=stats.chi2.cdf(X2,N);\n",
    "Po=1.0-stats.chi2.cdf(X2,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\" % (Pi) );\n",
    "print(\"   The Probability that X2 is outside the interval 0 to %.4f\" % (X2) );\n",
    "print(\"     is %.4f\" % (Po));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "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": [
    "# edapy_12_07: statistical tests of particle size data\n",
    "\n",
    "# Part A: first dataset\n",
    "\n",
    "# load data\n",
    "D=np.genfromtxt('../Data/testA.txt', delimiter='\\t');\n",
    "s = np.shape(D);\n",
    "NA = s[0];\n",
    "dobsA=eda_cvec(D);\n",
    "\n",
    "dtrueA = 100.0;     # known size of calibrations standard'\n",
    "dbartrueA = dtrueA; # known mean of measurement\n",
    "sd2trueA=1.0;       # known variance of measurement\n",
    "sdtrueA=1.0;        # 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 = (dbarestA-dbartrueA)/(sdtrueA/sqrt(NA));\n",
    "print(\"Z is %.4f\" % (ZA) );\n",
    "\n",
    "PZA = 1 - (stats.norm.cdf(abs(ZA),0.0,1.0)-stats.norm.cdf(-abs(ZA),0.0,1.0));\n",
    "print(\"Probability of |Z| exceeding %.4f is %.4f\" % (abs(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 = (dbarestA-dbartrueA)/(sdestA/sqrt(NA));\n",
    "print(\"t is %.4f\" % (tA) );\n",
    "\n",
    "PtA = 1.0 - (stats.t.cdf(abs(tA),NA)-\n",
    "             stats.t.cdf(-abs(tA),NA));\n",
    "print(\"Probability of |t| exceeding %.4f is %.4f\" % (abs(tA), PtA) );\n",
    "\n",
    "print(' ');\n",
    "\n",
    "# Part B: second dataset\n",
    "\n",
    "# load data\n",
    "D=np.genfromtxt('../Data/testB.txt', delimiter='\\t');\n",
    "s = np.shape(D);\n",
    "NB = s[0];\n",
    "dobsB=eda_cvec(D);\n",
    "\n",
    "dtrueB = 100.0;           # 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 = (dbarestB-dbartrueB)/(sdtrueB/sqrt(NB));\n",
    "print(\"Z is %.4f\" % (ZB) );\n",
    "\n",
    "PZB = 1 - (stats.norm.cdf(abs(ZB),0.0,1.0)-stats.norm.cdf(-abs(ZB),0.0,1.0));\n",
    "print(\"Probability of |Z| exceeding %.4f is %.4f\" % (abs(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 = (dbarestB-dbartrueB)/(sdestB/sqrt(NB));\n",
    "print(\"t is %.4f\" % (tB) );\n",
    "\n",
    "PtB = 1.0 - ( stats.t.cdf(abs(tB),NB) - stats.t.cdf(-abs(tB),NB) );\n",
    "print(\"Probability of |t| exceeding %.4f is %.4f\" % (abs(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 = 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(abs(Z),0,1)-stats.norm.cdf(-abs(Z),0,1));\n",
    "print(\"Probability of |Z| exceeding %.4f is %.4f\" % (abs(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 = 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(abs(t),M)-stats.t.cdf(-abs(t),M));\n",
    "print(\"Probability of |t| exceeding %.4f is %.4f\" % (abs(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",
    "PF = 1 - (stats.f.cdf(F,NA,NB)-stats.f.cdf(1/F,NA,NB));\n",
    "print(\"1/F is %.4f and F is %.4f\" % (1/F, F) );\n",
    "print(\"Probability of F outside interval (%.4f,%.4f) is %.4f\" % (1/F, F, PF) );\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 20072 characters]\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 20288 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 cannot be excluded to 95% confidence\n"
     ]
    }
   ],
   "source": [
    "# edapy_12_08: F-test to assess difference in fit\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",
    "# load data from file\n",
    "D=np.genfromtxt('../Data/brf_fragment.txt', delimiter='\\t');\n",
    "s = np.shape(D);\n",
    "N = s[0];\n",
    "t = eda_cvec(10.0+np.linspace(0,N-1,N));  # hours from start of orignal record\n",
    "dobs = eda_cvec(D);\n",
    "\n",
    "# plot setup\n",
    "fig1 = plt.figure(1,figsize=(10,6));\n",
    "plt.plot(t,dobs,'ko');\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;   # number of mode parameters\n",
    "G=np.zeros((N,M));  # set up data kernel\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) ); # least squares fit\n",
    "dpreA = np.matmul( G, mestA ); # predicted data\n",
    "plt.plot(t,dpreA,'k--');\n",
    "EA = np.matmul( (dobs-dpreA).T, (dobs-dpreA) );  # total error\n",
    "vA = N-M; # degrees of freedom\n",
    "plt.show();\n",
    "\n",
    "# plot setup\n",
    "fig2 = plt.figure(2,figsize=(10,6));\n",
    "plt.plot(t,dobs,'ko');\n",
    "plt.title('cubic fit');\n",
    "plt.xlabel('time t, hours');\n",
    "plt.ylabel('d(i)');\n",
    "\n",
    "# fit 2, cubic fit\n",
    "M=4; # number of model parameters\n",
    "G=np.zeros((N,M));  # set up data kernel\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) ); # least squares solution\n",
    "dpreB = np.matmul( G, mestB ); # predicted data\n",
    "plt.plot(t,dpreB,'k--');\n",
    "EB = np.matmul( (dobs-dpreB).T, (dobs-dpreB) ); # total error\n",
    "vB = N-M; # degrees of freedom\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",
    "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",
    "if( P>0.05):\n",
    "    print('Null hypothesis cannot be excluded to 95% confidence')\n",
    "else:\n",
    "    print('Null hypothesis can be excluded to 95% confidence')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "true      data variance 64.0000\n",
      "estimated data variance 69.7200\n",
      "variance of data estimated from PSD 69.83\n",
      " \n",
      "power spectral density\n",
      "true mean 3.1938 variance 10.2001\n",
      "estimated mean 3.4844 variance 13.4551\n",
      " \n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 48984 characters]\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 53340 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",
      "94.1520 percent of points are below 0.9500 confidence\n",
      " \n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 28004 characters]\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_12_09: variance of power spectral density (untapered)\n",
    "# case of untapered random timeseries\n",
    "\n",
    "# generic time set up of t-axis\n",
    "N=floor(2**10);  # length of time vector\n",
    "Dt=0.025;        # sampling inteval\n",
    "T = N*Dt;        # duration\n",
    "t = eda_cvec(Dt*np.linspace(0,N-1,N)); # time vector\n",
    "\n",
    "# create an uncorrelated random timeseries\n",
    "sd2true=64.0;          # variance\n",
    "sdtrue=sqrt(sd2true);  # standard deviation\n",
    "draw = np.random.normal(0,sdtrue,(N,1)); # time series of random data\n",
    "\n",
    "# generic frequency set up\n",
    "No2 = floor(N/2);   # integer version of N.2\n",
    "fmax=1/(2.0*Dt);    # nyquist frequency\n",
    "Df=fmax/No2;        # frequency sampling\n",
    "Nf=No2+1;           # number of non-negative frequencies\n",
    "# full (pos and neg) frequency vector\n",
    "f=eda_cvec(Df*np.concatenate((np.linspace(0,No2,Nf),np.linspace(-No2+1,-1,No2-1)),axis=0));\n",
    "Dw=2.0*pi*Df;       # angular feqweuncy sampling\n",
    "w=2.0*pi*f;         # full (pos and neg) angular fedquency vector\n",
    "Nw=Nf;              # number of non-negative angular frequencies\n",
    "fpos=eda_cvec(Df * np.linspace(0,No2,Nf)); # vector of non-negative frequencies\n",
    "wpos=2.0*pi*fpos;   # vector of non-negative angular frequencies\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;  # estimated variance\n",
    "print(\"true      data variance %.4f\" % (ff*sd2true) );\n",
    "print(\"estimated data variance %.4f\" % (sd2est));\n",
    "\n",
    "# compute power spectral density\n",
    "fftraw = np.fft.fft(d,axis=0);                # Fourier transform\n",
    "dbar = eda_cvec(Dt*fftraw[0:Nf,0]);           # keep non-negative frequencies 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);        # mean of p.s.d.\n",
    "s2varest=np.std(s2)**2;         # variance of p.s.d.\n",
    "s2sigmaest=sqrt(s2varest);      # standard deviation\n",
    "sd2psd=Nf*Df*s2meanest;         # variance of data estimated from PSD \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;            # number of bins in histogram\n",
    "s2min = np.min(s2);   # minumun bin value\n",
    "s2max = np.max(s2);   # maximum bin value\n",
    "ct, ed = np.histogram(s2,Nhist,(s2min,s2max)); # histogram\n",
    "Nc = len(ct);         # length of count vector ct\n",
    "Ne = len(ed);         # length of edges vector ed\n",
    "s2hist = eda_cvec(ct);                          # counts as a vector   \n",
    "centers = eda_cvec(0.5*(ed[1:Ne]+ed[0:Ne-1]));  # bins centers as a vector\n",
    "\n",
    "# chi-squared pdf\n",
    "Db=(centers[1,0]-centers[0,0])/c;              # sampling divided by c\n",
    "s2histtrue=Nf*Db*stats.chi2.pdf(centers/c,p);  # normalied chi-squared p.d.f.\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",
    "# plot histograms\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": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "true      data variance 25.4092\n",
      "estimated data variance 24.5415\n",
      "variance of data estimated from PSD 24.58\n",
      " \n",
      "power spectral density\n",
      "true mean 1.2680 variance 1.6078\n",
      "estimated mean 1.2265 variance 1.3261\n",
      " \n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 42708 characters]\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 61232 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.9064 percent of points are below 0.9500 confidence\n",
      " \n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 38700 characters]\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_12_10: variance of power spectral density (tapered)\n",
    "# case of random timeseries tapered with Hamming window\n",
    "\n",
    "# generic time set up\n",
    "N=floor(2**10);  # number of points in time series\n",
    "Dt=0.025;        # sampling interval\n",
    "T = N*Dt;        # duration\n",
    "t = eda_cvec(Dt*np.linspace(0,N-1,N)); # time vector\n",
    "\n",
    "# create an uncorrelated random timeseries\n",
    "sd2true=64.0;          # variance\n",
    "sdtrue=sqrt(sd2true);  # std dev\n",
    "draw = np.random.normal(0,sdtrue,(N,1)); # time series of random data\n",
    "\n",
    "# generic frequency set up\n",
    "No2 = floor(N/2);   # integer version of N.2\n",
    "fmax=1/(2.0*Dt);    # nyquist frequency\n",
    "Df=fmax/No2;        # frequency sampling\n",
    "Nf=No2+1;           # number of non-negative frequencies\n",
    "# full (pos and neg) frequency vector\n",
    "f=eda_cvec(Df*np.concatenate((np.linspace(0,No2,Nf),np.linspace(-No2+1,-1,No2-1)),axis=0));\n",
    "Dw=2.0*pi*Df;       # angular feqweuncy sampling\n",
    "w=2.0*pi*f;         # full (pos and neg) angular fedquency vector\n",
    "Nw=Nf;              # number of non-negative angular frequencies\n",
    "fpos=eda_cvec(Df * np.linspace(0,No2,Nf)); # vector of non-negative frequencies\n",
    "wpos=2.0*pi*fpos;   # vector of non-negative angular frequencies                              \n",
    "\n",
    "#  Hamming window function\n",
    "W = eda_cvec( 0.54 - 0.46*np.cos(2*pi*np.linspace(0,N-1,N)/(N-1)) );\n",
    "d = np.multiply( W, draw );\n",
    "# window changes variance of signal by ff\n",
    "\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);                # Fourier transform\n",
    "dbar = eda_cvec(Dt*fftraw[0:Nf,0]);           # keep non-negative frequencies, 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);                      # mean\n",
    "s2varest=np.std(s2)**2;                       # vaiance\n",
    "s2sigmaest=sqrt(s2varest);                    # standard deviation\n",
    "sd2psd=Nf*Df*s2meanest;                       # variance of data estimated from PSD\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;                # mean of chi-squaed p.d.f., with normalization\n",
    "s2vartrue = 2*p*(c**2);          # variance of chi-squaed p.d.f., with normalization\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 p.s.d.\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               # number of bins in histogram\n",
    "s2min = np.min(s2);     # start bin\n",
    "s2max = np.max(s2);     # end bin\n",
    "ct, ed = np.histogram(s2,Nhist,(s2min,s2max)); # histogram\n",
    "Nc = len(ct);           # number of counts ct\n",
    "Ne = len(ed);           # number of edges ed\n",
    "s2hist = eda_cvec(ct);  # counts as a vector\n",
    "centers = eda_cvec(0.5*(ed[1:Ne]+ed[0:Ne-1])); # edges as a vector\n",
    "\n",
    "# chi-squared pdf\n",
    "Db=(centers[1,0]-centers[0,0])/c;               # normalization factor\n",
    "s2histtrue=Nf*Db*stats.chi2.pdf(centers/c,p);   # true hisogram from chi-squared p.d.f.\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": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "true      data variance 25.4092\n",
      "estimated data variance 24.8075\n",
      "variance of data estimated from PSD 24.81\n",
      " \n",
      "power spectral density\n",
      "true mean 1.2680 variance 1.6078\n",
      "estimated mean 1.2380 variance 1.9533\n",
      " \n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 47364 characters]\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 47536 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.5166 percent of points are below 0.9500 confidence\n",
      " \n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 25980 characters]\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "largest spectral peak: 1.668021e+01\n",
      "probability of laregst peak 0.99999806\n",
      "probability of PSD exceeding maximum observed value 0.00000194\n",
      "Nf probability of PSD exceeding maximum observed value 0.00099263\n"
     ]
    }
   ],
   "source": [
    "# edapy_12_11, confidence of spectral peak\n",
    "# for random noise plus a small ampliude sinusoid\n",
    "\n",
    "# generic time set up\n",
    "N=floor(2**10);  # number of points in time series\n",
    "Dt=0.025;        # sampling interval\n",
    "T = N*Dt;        # duration\n",
    "t = eda_cvec(Dt*np.linspace(0,N-1,N)); # time vector\n",
    "\n",
    "# generic frequency set up\n",
    "No2 = floor(N/2);   # integer version of N.2\n",
    "fmax=1/(2.0*Dt);    # nyquist frequency\n",
    "Df=fmax/No2;        # frequency sampling\n",
    "Nf=No2+1;           # number of non-negative frequencies\n",
    "# full (pos and neg) frequency vector\n",
    "f=eda_cvec(Df*np.concatenate((np.linspace(0,No2,Nf),np.linspace(-No2+1,-1,No2-1)),axis=0));\n",
    "Dw=2.0*pi*Df;       # angular feqweuncy sampling\n",
    "w=2.0*pi*f;         # full (pos and neg) angular fedquency vector\n",
    "Nw=Nf;              # number of non-negative angular frequencies\n",
    "fpos=eda_cvec(Df * np.linspace(0,No2,Nf)); # vector of non-negative frequencies\n",
    "wpos=2.0*pi*fpos;   # vector of non-negative angular frequencies \n",
    "\n",
    "# low amplitude periodic signal buried in the noise\n",
    "sd2true=64.0;          # variance\n",
    "sdtrue=sqrt(sd2true);  # standard deviation\n",
    "# data = signal + noise\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 = eda_cvec(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.fft.fft(d,axis=0);         # Fourier transform\n",
    "dbar = eda_cvec(Dt*fftraw[0:Nf,0]);    # keep non-negative frequencies, 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);     # mean\n",
    "s2varest=np.std(s2)**2;      # variance\n",
    "s2sigmaest=sqrt(s2varest);   # standard deviation\n",
    "sd2psd=Nf*Df*s2meanest;      # variance of data estimated from PSD\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;                # true mean of normalizenchi-suqared distribution\n",
    "s2vartrue = 2*p*(c**2);          # true varoance of normalized chi-suqared distribution\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;            # number of bins in histogram\n",
    "s2min = np.min(s2);   # start bin\n",
    "s2max = np.max(s2);   # end bin\n",
    "ct, ed = np.histogram(s2,Nhist,(s2min,s2max)); # histogram\n",
    "Nc = len(ct);         # number of counts ct\n",
    "Ne = len(ed);         # number of edges ed\n",
    "s2hist = eda_cvec(ct);                         # counts as a vector\n",
    "centers = eda_cvec(0.5*(ed[1:Ne]+ed[0:Ne-1])); # bin centers as a vector\n",
    "\n",
    "# chi-squared pdf\n",
    "Db=(centers[1,0]-centers[0,0])/c;\n",
    "# true histogram derived from normalized chi-squared distribution\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([s2meantrue, s2meantrue], [0, 0.1*max(s2hist)],'k:');\n",
    "plt.plot([cv*c, cv*c], [0, 0.1*max(s2hist)],'k:');\n",
    "plt.show();\n",
    "\n",
    "# probability that PSD exceeds maximum observed value\n",
    "speak = np.max(s2);\n",
    "print(\"largest spectral peak: %e\" % (speak) );\n",
    "ppeak = stats.chi2.cdf(speak/c,p);\n",
    "print(\"probability of laregst peak %.8f\" % (ppeak) );\n",
    "print(\"probability of PSD exceeding maximum observed value %.8f\" % (1-ppeak) );\n",
    "print(\"Nf probability of PSD exceeding maximum observed value %.8f\" % (1-(ppeak**Nf)) );"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p.s.d. of AR1 process with these betas\n",
      "[[0.25 0.5  0.75]]\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 30760 characters]\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_12_12, p.s.d. of AR1 process\n",
    "\n",
    "# tunable constants\n",
    "beta = 0.60;     # beta or AR1 process\n",
    "sigma = 100;     # sigma of x\n",
    "N=1024;          # length of time series\n",
    "\n",
    "# standard set up a time axis, t\n",
    "Dt=0.6;\n",
    "T = N*Dt;\n",
    "t= eda_cvec(Dt*np.linspace(0,N-1,N)); \n",
    "\n",
    "# standard setup of frequency axis, f\n",
    "No2 = floor(N/2);   # integer version of N.2\n",
    "fmax=1/(2.0*Dt);    # nyquist frequency\n",
    "Df=fmax/No2;        # frequency sampling\n",
    "Nf=No2+1;           # number of non-negative frequencies\n",
    "# full (pos and neg) frequency vector\n",
    "f=eda_cvec(Df*np.concatenate((np.linspace(0,No2,Nf),np.linspace(-No2+1,-1,No2-1)),axis=0));\n",
    "Dw=2.0*pi*Df;       # angular feqweuncy sampling\n",
    "w=2.0*pi*f;         # full (pos and neg) angular fedquency vector\n",
    "Nw=Nf;              # number of non-negative angular frequencies\n",
    "fpos=eda_cvec(Df * np.linspace(0,No2,Nf)); # vector of non-negative frequencies\n",
    "wpos=2.0*pi*fpos;   # vector of non-negative angular frequencies \n",
    "\n",
    "# delete nonnegative frequencies from axes\n",
    "f=f[0:Nf,0:1];\n",
    "w=w[0:Nf,0:1];\n",
    "\n",
    "# plot of spectra with different beta's\n",
    "betalist = eda_cvec( 0.25, 0.5, 0.75);\n",
    "NB,i = np.shape(betalist);\n",
    "print('p.s.d. of AR1 process with these betas');\n",
    "print(betalist.T);\n",
    "\n",
    "fig1 = plt.figure(1,figsize=(10,6));\n",
    "plt.axis( [0, fmax, 0, 1/(1+(betalist[NB-1,0]-2*betalist[NB-1,0])) ]  );\n",
    "plt.xlabel('f');\n",
    "plt.ylabel('1/|v(f)|^2');\n",
    "\n",
    "for beta in betalist:\n",
    "    denom = (1 + (beta**2) - 2*beta*np.cos(2*pi*0.5*f/fmax) );\n",
    "    plt.plot( f, np.divide(1,denom), 'k-', lw=2 );\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "plot time series x(t) and y(t)\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 75872 characters]\n",
      "text/plain": [
       "<Figure size 720x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "plot enlarged portions of time series x(t) and y(t)\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 64776 characters]\n",
      "text/plain": [
       "<Figure size 720x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "plot power spretral density, |x(f)|^2 and |y(f)|^2\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 59892 characters]\n",
      "text/plain": [
       "<Figure size 720x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "plot power spretral density, |x(f)|**2 and |y(f)|**2\n",
      "with other stuff overlain\n",
      "scaling constant c: 5.988304e+03\n",
      "pc: 1.197661e+04\n",
      "2pc**2: 1.434391e+08\n",
      "Approximate 95% confidence level used\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 70652 characters]\n",
      "text/plain": [
       "<Figure size 720x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_12_13, 95% confidence limit of AR1 process\n",
    "\n",
    "# tunable constants\n",
    "beta =    0.6;   # beta of AR1 process\n",
    "sigma = 100.0;   # sigma of x\n",
    "N=1024;          # length of time series\n",
    "tmaxplot=120.0;  # max time on enlarged plot\n",
    "ADDSINUSOID=0;   # add sinusoid to random timeseries on 1\n",
    "EXACT95=0;       # Exact confidence on 1, (mean+2sigma) approx on 0\n",
    "\n",
    "# standard set up a time axis, t\n",
    "Dt=0.6;\n",
    "T = N*Dt;\n",
    "t=eda_cvec(Dt*np.linspace(0,N-1,N)); \n",
    "\n",
    "# standard setup of frequency axis, f\n",
    "No2 = floor(N/2);   # integer version of N.2\n",
    "fmax=1/(2.0*Dt);    # nyquist frequency\n",
    "Df=fmax/No2;        # frequency sampling\n",
    "Nf=No2+1;           # number of non-negative frequencies\n",
    "# full (pos and neg) frequency vector\n",
    "f=eda_cvec(Df*np.concatenate((np.linspace(0,No2,Nf),np.linspace(-No2+1,-1,No2-1)),axis=0));\n",
    "Dw=2.0*pi*Df;       # angular feqweuncy sampling\n",
    "w=2.0*pi*f;         # full (pos and neg) angular fedquency vector\n",
    "Nw=Nf;              # number of non-negative angular frequencies\n",
    "fpos=eda_cvec(Df * np.linspace(0,No2,Nf)); # vector of non-negative frequencies\n",
    "wpos=2.0*pi*fpos;   # vector of non-negative angular frequencies \n",
    "\n",
    "# uncorrelated random time series x(t)\n",
    "x = np.random.normal(loc=0,scale=sigma,size=(N,1));\n",
    "\n",
    "# AR1 process y(t)\n",
    "y = np.zeros((N,1));\n",
    "y[0,0] = x[0,0];\n",
    "for i in range(1,N):\n",
    "    y[i,0] = beta*y[i-1,0] + x[i,0];\n",
    "\n",
    "if( ADDSINUSOID ):\n",
    "    sigmay2 = np.var(y);\n",
    "    sigmay = sqrt(sigmay2);\n",
    "    Ac=0.25*sigmay;\n",
    "    fc=fmax/4.0;\n",
    "    y = y + Ac*np.cos(2.0*pi*fc*t);\n",
    "    print(\"A small amplitude sinusoid of amplitude %.2f and frequency %.2f has been added\" % (Ac,fc));\n",
    "\n",
    "# plot entire time series\n",
    "print('plot time series x(t) and y(t)');\n",
    "fig1 = plt.figure(1,figsize=(10,6));\n",
    "\n",
    "ax1 = plt.subplot(2,1,1);\n",
    "plt.axis( [0, Dt*(N-1), min(x), max(x) ] );\n",
    "plt.plot( t, x, 'k-', lw=2 );\n",
    "plt.xlabel('t');\n",
    "plt.ylabel('x(t)');\n",
    "ax1 = plt.subplot(2,1,2);\n",
    "plt.axis( [0, Dt*(N-1), min(y), max(y) ] );\n",
    "plt.plot( t, y, 'k-', 'LineWidth', 2 );\n",
    "plt.xlabel('t');\n",
    "plt.ylabel('y(t)');\n",
    "plt.show();\n",
    "\n",
    "# plot enlarged portion of time series\n",
    "print('plot enlarged portions of time series x(t) and y(t)');\n",
    "fig2 = plt.figure(2,figsize=(10,6));\n",
    "ax1 = plt.subplot(2,1,1);\n",
    "plt.axis( [0, tmaxplot, min(x), max(x) ] );\n",
    "plt.plot( t, x, 'k-', lw=2 );\n",
    "plt.xlabel('t');\n",
    "plt.ylabel('x(t)');\n",
    "ax1 = plt.subplot(2,1,2);\n",
    "plt.axis( [0, tmaxplot, min(y), max(y) ] );\n",
    "plt.plot( t, y, 'k-', lw=2 );\n",
    "plt.xlabel('t');\n",
    "plt.ylabel('y(t)');\n",
    "plt.show();\n",
    "\n",
    "# compute power spectral density\n",
    "xfftraw = np.fft.fft(x,axis=0);         # raw DFT output\n",
    "xbar = Dt*xfftraw[0:Nf,0:1];            # Fourier transform, + freqs only\n",
    "xs2 = (2/T)*np.power(np.abs(xbar),2);   # power spectral density\n",
    "xs2=xs2[0:Nf,0:1];                      # nonnegative f's only\n",
    "yfftraw = np.fft.fft(y,axis=0);         # raw DFT output\n",
    "ybar = Dt*yfftraw[0:Nf,0:1];            # Fourier transform, + freqs only\n",
    "ys2 = (2/T) * np.power(abs(ybar),2);    # power spectral density\n",
    "\n",
    "# delete nonnegative frequencies from axes\n",
    "f=f[0:Nf];\n",
    "w=w[0:Nf];\n",
    "\n",
    "# plot of spectra\n",
    "print('plot power spretral density, |x(f)|^2 and |y(f)|^2');\n",
    "fig3 = plt.figure(3,figsize=(10,6));\n",
    "ax1 = plt.subplot(2,1,1);\n",
    "plt.axis( [0, np.max(f), 0, np.max(xs2) ] );\n",
    "plt.plot( f, xs2, 'k-', lw=2 );\n",
    "plt.xlabel('f');\n",
    "plt.ylabel('|x(f)|**2');\n",
    "ax1 = plt.subplot(2,1,2);\n",
    "plt.axis( [0, np.max(f), 0, np.max(ys2) ] );\n",
    "plt.plot( f, ys2, 'k-', lw=2 );\n",
    "plt.xlabel('f');\n",
    "plt.ylabel('|y(f)|**2');\n",
    "plt.show();\n",
    "\n",
    "# plot of spectra, onto which other stiff will be overlain\n",
    "print('plot power spretral density, |x(f)|**2 and |y(f)|**2');\n",
    "print('with other stuff overlain');\n",
    "fig4 = plt.figure(4,figsize=(10,6));\n",
    "ax1 = plt.subplot(2,1,1);\n",
    "plt.axis( [0, np.max(f), 0, np.max(xs2) ] );\n",
    "plt.plot( f, xs2, 'k-', lw=2 );\n",
    "plt.xlabel('f');\n",
    "plt.ylabel('|x(f)|**2');\n",
    "ax1 = plt.subplot(2,1,2);\n",
    "plt.axis( [0, np.max(f), 0, np.max(ys2) ] );\n",
    "plt.plot( f, ys2, 'k-', lw=2 );\n",
    "plt.xlabel('f');\n",
    "plt.ylabel('|y(f)|**2');\n",
    "\n",
    "Nr = 1000;\n",
    "Nr95 = 950;\n",
    "Nr50 = 500;\n",
    "S = np.zeros((Nf,Nr));\n",
    "\n",
    "# repeat for Nr realizations\n",
    "for ir in range(Nr):\n",
    "    # uncorrelated time series\n",
    "    x = np.random.normal(loc=0.0,scale=sigma,size=(N,1));\n",
    "    # build AR1 time sereis\n",
    "    y = np.zeros((N,1));\n",
    "    y[0,0] = x[0,0];\n",
    "    for i in range(1,N):\n",
    "        y[i,0] = beta*y[i-1,0] + x[i,0];\n",
    "    yfftraw = np.fft.fft(y,axis=0);  # FFT of y\n",
    "    ybar = Dt*yfftraw[0:Nf,0:1]; # FT of nonnegative f's\n",
    "    ys2 = (2/T) * np.power(np.abs(ybar),2); # p.s.d.\n",
    "    S[0:Nf,ir:ir+1]=ys2[0:Nf,0:1]; # save p.s.d.\n",
    "\n",
    "# statistics of the ensemble of p.s.d.'s\n",
    "S95 = np.zeros((Nf,1));   # 95 percent level\n",
    "Smean = np.zeros((Nf,1)); # mean\n",
    "for ifr in range(Nf):   # loop over frequencies\n",
    "    v = S[ifr:ifr+1,0:Nr].T; # the p.s.d. at one frequency\n",
    "    v = eda_cvec(np.sort(v,axis=0)); # sort the values\n",
    "    S95[ifr,0]=v[Nr95,0];  # 95% less than this value\n",
    "    Smean[ifr,0]=np.mean(v); # mean\n",
    "\n",
    "plt.plot( f, S95,   '-', lw=5, color=[0.8,0.8,0.8] );\n",
    "plt.plot( f, Smean, '-', lw=5, color=[0.8,0.8,0.8] );\n",
    "\n",
    "# theoretical spectrum\n",
    "p=2;    # degrees of freedom\n",
    "# in this code, the window function is just w(i)=1\n",
    "ff = 1.0; # variance of window function\n",
    "c = ( ff*(sigma**2) ) / (2*Nf*Df);  # scaling constant\n",
    "xs2meantrue = p*c;        # mean of p.s.d. of x\n",
    "xs2vartrue = 2*p*(c**2);  # variance of p.s.d. of x\n",
    "print(\"scaling constant c: %e\" % (c) );\n",
    "print(\"pc: %e\" % (xs2meantrue) );\n",
    "print(\"2pc**2: %e\" % (xs2vartrue) );\n",
    "# the AR1 filter v(f)\n",
    "denom = (1.0 + (beta**2) - 2.0*beta*np.cos(pi*f/fmax) );\n",
    "# mean of p.s.d. of y\n",
    "ys2mean =  xs2meantrue*np.reciprocal(denom);\n",
    "# variance of p.s.d. of y\n",
    "if( EXACT95 ):\n",
    "    C95 = stats.chi2.ppf(0.95,p);\n",
    "    ys295 = C95*c*np.reciprocal(denom);\n",
    "    print('Exact 95%% confidence level used');\n",
    "else:\n",
    "    ys295 = ys2mean + 2.0*sqrt(xs2vartrue)*np.reciprocal(denom);\n",
    "    print('Approximate 95% confidence level used');\n",
    "plt.plot( f, ys2mean, 'k--', lw=2 );\n",
    "plt.plot( f, ys295, 'k-', lw=2 );\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 23864 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.5250 +/- 0.0176 (2 sigma)\n",
      "bootstrap method:  slope between 0.5081 and 0.5429 (95 percent confidence)\n"
     ]
    }
   ],
   "source": [
    "# edapy_12_14: bootstrap statistics of slope\n",
    "# in least squares fit of a straight line to a small\n",
    "# fragment of Black Rock Forest temperature data\n",
    "\n",
    "# load data and set up time axis\n",
    "D=np.genfromtxt('../Data/brf_fragment.txt', delimiter='\\t');\n",
    "s = np.shape(D);\n",
    "N = s[0];  # number of data\n",
    "torig = eda_cvec(10+np.linspace(0,N-1,N));  # time vector, hours from start of orignal record\n",
    "dorig = eda_cvec(D);                        # data vector\n",
    "\n",
    "# least-squares slope and its variance of original data\n",
    "M=2;                    # number of model parameters\n",
    "G=np.zeros((N,M));      # data kernel setup\n",
    "G[:,0]=np.ones((N,));\n",
    "G[:,1]=torig.ravel();\n",
    "GTG = np.matmul(G.T,G); # G_transpose time G\n",
    "mest = la.solve( GTG, np.matmul(G.T,dorig) ); # least squares solution\n",
    "slopeorig=mest[1,0];    # break out estimated slope\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);  # individual errors\n",
    "E = np.matmul(e.T,e);           # total error\n",
    "sigma2dorig = E/N;              # posterior variance of data\n",
    "Cmorig = sigma2dorig * la.inv(GTG);  # posterior covariance of model parametres\n",
    "sigma2morig = Cmorig[1,1];      # breal out variance of slope\n",
    "\n",
    "# number of realizations;\n",
    "Nr = 10000;\n",
    "slope = np.zeros((Nr,1));    # vecor of bootstraped slopes\n",
    "\n",
    "# loop over realizations\n",
    "for p in range(Nr):\n",
    "    # resample\n",
    "    # N random integers, each bewteen 1 and N-1\n",
    "    rowindex=np.random.randint(N,size=N); \n",
    "    # random resampling of time with duplication\n",
    "    t = eda_cvec(torig[rowindex,0]);\n",
    "    # random resampling of data with duplication\n",
    "    d = eda_cvec(dorig[rowindex,0]);\n",
    "    # straight line fit\n",
    "    M=2;                     # number of model parameters\n",
    "    G=np.zeros((N,M));       # data akernel setup\n",
    "    G[:,0]=np.ones((N,));\n",
    "    G[:,1]=t.ravel();\n",
    "    GTG = np.matmul(G.T,G);  # G-transpose times G\n",
    "    mest=la.solve(GTG,np.matmul(G.T,d)); # least sq soln\n",
    "    slope[p,0]=mest[1,0];    # save slope in vector\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;      # number of bins\n",
    "slopemin = 0.2;  # staring bin\n",
    "slopemax = 0.8;  # ending bin\n",
    "ct, ed = np.histogram(slope,Nhist,(slopemin,slopemax)); # histogram\n",
    "Nc = len(ct);    # length of counts ct\n",
    "Ne = len(ed);    # length of edges ed\n",
    "slopehist = eda_cvec(ct); # counts as a vector\n",
    "# bin centers as a vector\n",
    "centers = eda_cvec(0.5*(ed[1:Ne]+ed[0:Ne-1]));\n",
    "# approximate pdf via normalized histogram\n",
    "Ds = centers[1,0] - centers[0,0];  # bin sampling\n",
    "norm = Ds*np.sum(slopehist);       # normalization factor\n",
    "pbootstrap =slopehist/norm;        # p.d.f. from histogram\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",
    "# approximate integral for mean as sum\n",
    "mb = Ds*np.sum( np.multiply( centers, pbootstrap) );\n",
    "# approximate integral for variance as sum\n",
    "vb = Ds*np.sum( np.multiply( np.power(centers-mb,2) , pbootstrap) ); \n",
    "\n",
    "# 5% and 95% confidence of\n",
    "# cumulative probability function\n",
    "Pbootstrap = eda_cvec(Ds*np.cumsum(pbootstrap));\n",
    "# find 0.025 cumulative probability\n",
    "ilovec = np.where(Pbootstrap>=0.025); \n",
    "# np.where() returns a bunch of stuff,\n",
    "# but what we want, ilo, is the first thing\n",
    "ilox = ilovec[0];\n",
    "ilo = ilox[0];\n",
    "blo = centers[ilo];\n",
    "# find 0.975 cumulative probability\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": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CaO/Na2O ratio of Factor 2\n",
      "mean 0.485601 variance 0.000192\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 20656 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.4856 +/- 0.0277 (2 sigma)\n",
      "bootstrap method:  ratio between 0.4574 and 0.5126 (95 percent confidence)\n"
     ]
    }
   ],
   "source": [
    "# edapy_12_15: bootstrap statistics using Atlantic Rock data set\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('../Data/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",
    "    # randomly resample\n",
    "    rowindex=np.random.randint(N,size=N);  # N random integers, each between 0 and N-1\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",
    "# make histogram of ratios\n",
    "Nhist=500;       # number of bins\n",
    "ratiomin = 0.2;  # start bin\n",
    "ratiomax = 0.8;  # end bin\n",
    "ct, ed = np.histogram(ratio,Nhist,(ratiomin,ratiomax)); # histogram\n",
    "Nc = len(ct);    # number of counts ct\n",
    "Ne = len(ed);    # number rof edges ed\n",
    "ratiohist = eda_cvec(ct);  # counts as a vector\n",
    "centers = eda_cvec(0.5*(ed[1:Ne]+ed[0:Ne-1])); # bn centers as a vector\n",
    "\n",
    "# approximate pdf via normalized histogram\n",
    "Dr = centers[1,0] - centers[0,0];   # bin spacing\n",
    "norm = Dr*np.sum(ratiohist);        # normalization of p.d.f.\n",
    "pbootstrap = ratiohist/norm;        # p.d.f. estimated from histogram\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 = eda_cvec(Dr*np.cumsum(pbootstrap));\n",
    "\n",
    "# 2.5% connfidence\n",
    "ilovec = np.where( (Pbootstrap>=0.025).ravel() );\n",
    "# np.where() returns a bunch of stuff, we want the first, ilo\n",
    "ilox = ilovec[0];\n",
    "ilo = ilox[0];\n",
    "blo = centers[ilo];\n",
    "\n",
    "# 97.5% connfidence\n",
    "ihivec = np.where(Pbootstrap>=0.975);\n",
    "# np.where() returns a bunch of stuff, we want the first, ihi\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"
   ]
  },
  {
   "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
}
