ex2p1_ex2p2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "### Inverse Problems and Parameter Estimation, GEOS 627/427, University of Alaska Fairbanks\n",
    "\n",
    "- script ex2p1_ex2p2.ipynb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%run lib_header.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "import scipy.io\n",
    "import scipy.special\n",
    "from numpy.linalg import inv\n",
    "\n",
    "from lib_geos import corrcov\n",
    "from lib_inverse import plot_ellipse\n",
    "from lib_peip import * # chi2cdf, chi2inv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load data\n",
    "dat    = scipy.io.loadmat('data/data1.mat')\n",
    "data1  = dat['data1']\n",
    "print(f\"data1: {data1.shape}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ndata  = len(data1)\n",
    "ta     = data1[:,0]\n",
    "ya     = data1[:,1]\n",
    "sigmaa = data1[:,2]\n",
    "t      = ta.reshape(ndata,1)\n",
    "y      = ya.reshape(ndata,1)\n",
    "sigma  = sigmaa.reshape(ndata,1)\n",
    "ones   = np.ones((ndata,1))\n",
    "ndata  = len(t)\n",
    "\n",
    "# t-value for plotting continuous curves\n",
    "dt = 0.05\n",
    "tplot = np.arange(np.min(t)-1, np.max(t)+1+dt, dt)\n",
    "\n",
    "print('displaying t, y, sigma:')\n",
    "showmat(np.hstack((t,y,sigma)),2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the data and error bars\n",
    "plt.figure(figsize=(9,7))\n",
    "plt.grid()\n",
    "plt.errorbar(ta,ya,yerr=sigmaa, ls='none',marker='o',mfc='none',capsize=3,mec='k',ms=3)\n",
    "plt.xlabel('Time (s)')\n",
    "plt.ylabel('Elevation (m)');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Build the parabolic system matrix\n",
    "G = np.hstack(((ones, t, -0.5*t**2)))\n",
    "showmat(G)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Apply the weighting\n",
    "yw = y / sigma\n",
    "Gw = G / np.hstack((sigma,sigma,sigma))\n",
    "showmat(Gw,2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solve for the least-squares solution\n",
    "print('Least-squares solution, m =')\n",
    "m = inv(Gw.T@Gw) @ Gw.T @ yw\n",
    "showmat(m,1)\n",
    "m1 = m[0][0]\n",
    "m2 = m[1][0]\n",
    "m3 = m[2][0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# alternative functions: fit with polynomial\n",
    "# note: this does NOT take into account the sigma associated with each y-value\n",
    "P = np.polyfit(ta,ya,2)\n",
    "mpoly = np.flipud(P)\n",
    "showmat(mpoly,1)\n",
    "# check: the quadratic coefficient should match -0.5*m3\n",
    "-0.5*m3 - mpoly[2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the data and model predicted data\n",
    "\n",
    "yplot     = m1 + m2*tplot - 0.5*m3*tplot**2\n",
    "yplotpoly = mpoly[0] + mpoly[1]*tplot + mpoly[2]*tplot**2\n",
    "\n",
    "plt.figure(figsize=(9,7))\n",
    "plt.plot(tplot,yplot,'k')\n",
    "plt.plot(tplot,yplotpoly,'r--')\n",
    "plt.grid()\n",
    "#plt.plot(t,y,'o-')\n",
    "plt.errorbar(ta,ya,yerr=sigmaa, ls='none',marker='o',mfc='none',capsize=3,mec='k',ms=3)\n",
    "# vectors would have to be flattened, unfortunately:\n",
    "#plt.errorbar(t.flatten(),y.flatten(),yerr=sigma.flatten(), ls='none',marker='o',mfc='none',capsize=3,mec='k',ms=3)\n",
    "plt.xlabel('Time (s)')\n",
    "plt.ylabel('Elevation (m)');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ginv = inv(Gw.T@Gw) @ Gw.T\n",
    "showmat(ginv,2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Covariance matrix')\n",
    "covm = ginv @ ginv.T\n",
    "showmat(covm,2)\n",
    "C = covm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# see notes_tarantola.pdf\n",
    "PCONF = 0.95\n",
    "DELTA = scipy.special.erfinv(PCONF)*np.sqrt(2)\n",
    "print(DELTA)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the 1.96-sigma (95%) conf intervals\n",
    "print('%.1f%% parameter confidence intervals (m-delm, mest, m+delm)' % (PCONF*100))\n",
    "delm = DELTA*np.sqrt(np.diag(covm)).reshape((3,1))\n",
    "showmat(delm,2)\n",
    "showmat(np.hstack((m-delm, m, m+delm)),2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Because there are 3 parameters to estimate, we have ndata-3 degrees of freedom.\n",
    "dof = ndata-3\n",
    "chi2 = np.linalg.norm((y - G@m)/sigma)**2\n",
    "print('chi-square misfit for %i dof = %.2f' % (dof,chi2))\n",
    "\n",
    "# Find the p-value for this data set\n",
    "p = 1-chi2cdf(chi2,dof)\n",
    "print('chi-square p-value = %.2f' % p)\n",
    "\n",
    "print('correlation matrix =')\n",
    "r = corrcov(covm)\n",
    "showmat(r,2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Displaying Data and Model Fit (fig 1)')\n",
    "\n",
    "# Output covm and the eigenvalues/eigenvectors of covm.\n",
    "print('Covariance matrix for fitted parameters.')\n",
    "\n",
    "lam0,u0 = np.linalg.eig(inv(covm))\n",
    "# we need to sort them to match Aster\n",
    "# note: the sign of u1 is flipped\n",
    "indices = np.argsort(lam0)\n",
    "lam, u = lam0[indices], u0[:, indices]\n",
    "print('Eigenvalues of the covariance matrix:')\n",
    "showmat([lam],4)\n",
    "print('Eigenvectors of the covariance matrix:')\n",
    "showmat(u,2)\n",
    "\n",
    "print('%.1f%% confidence ellipsoid semiaxis lengths:' % (PCONF*100))\n",
    "semi_axes = np.sqrt(chi2inv(PCONF,3)/lam)\n",
    "showmat([semi_axes],2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Monte Carlo Section\n",
    "y0 = G@m       # unweighted data \n",
    "nreal = 1000\n",
    "mmc = np.zeros((3,nreal))\n",
    "chimc = np.zeros(nreal)\n",
    "Ginv = np.linalg.pinv(Gw)    # note: pseudoinverse of weighted G matrix\n",
    "\n",
    "for i in range(nreal):\n",
    "    # Generate a trial data set of perturbed, weighted data\n",
    "    noise = np.random.randn(ndata, 1) * sigma\n",
    "    ytrial = y0 + noise              # unweighted data + noise\n",
    "    ywtrial = ytrial / sigma         # weighted data + noise\n",
    "    # KEY COMMAND: SOLVE FOR m FOR A SET OF DATA WITH ERRORS ADDED\n",
    "    mmx = Ginv @ ywtrial\n",
    "    mmc[:,i] = mmx.flatten()\n",
    "    chimc[i] = np.linalg.norm((G @ mmx - ytrial) / sigma) ** 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the histogram of chi squared values\n",
    "plt.figure()\n",
    "plt.hist(chimc,30)\n",
    "plt.ylabel('N')\n",
    "plt.xlabel(r'$\\chi^2$')\n",
    "print('Displaying 1000 Monte-Carlo Chi-square Values (fig 2)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the histograms of the model parameters\n",
    "plt.figure(figsize=(10,5))\n",
    "\n",
    "plt.subplot(1,3,1)\n",
    "plt.hist(mmc[0,:])\n",
    "plt.title(r'$m_1$ (m)')\n",
    "\n",
    "plt.subplot(1,3,2)\n",
    "plt.hist(mmc[1,:])\n",
    "plt.title(r'$m_2$ (m/s)')\n",
    "\n",
    "plt.subplot(1,3,3)\n",
    "plt.hist(mmc[2,:])\n",
    "plt.title(r'$m_3$ (m/s$^2$)')\n",
    "print('Displaying Monte-Carlo Model Histograms (fig 3)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the realizations of each pair of model parameters with the other\n",
    "plt.figure(figsize=(10,5))\n",
    "\n",
    "plt.subplot(1,3,1)\n",
    "plt.plot(mmc[0,:],mmc[1,:],'k*')\n",
    "plt.xlabel(r'$m_1$ (m)')\n",
    "plt.ylabel(r'$m_2$ (m/s)')\n",
    "\n",
    "plt.subplot(1,3,2)\n",
    "plt.plot(mmc[0,:],mmc[2,:],'k*')\n",
    "plt.xlabel(r'$m_1$ (m)')\n",
    "plt.ylabel(r'$m_3$ (m/s$^2$)')\n",
    "\n",
    "plt.subplot(1,3,3)\n",
    "plt.plot(mmc[1,:],mmc[2,:],'k*')\n",
    "plt.xlabel(r'$m_2$ (m/s)')\n",
    "plt.ylabel(r'$m_3$ (m/s$^2$)')\n",
    "\n",
    "plt.subplots_adjust(wspace=0.5)\n",
    "plt.tight_layout()\n",
    "\n",
    "print('Displaying Projections of 1000 Monte-Carlo models (fig 4)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the 95% error ellipses for each pair of parameters\n",
    "# Note that because we're doing pairs of parameters there are 2 degrees of freedom in the Chi-square here, rather than 3.  \n",
    "\n",
    "deltachi2 = chi2inv(PCONF,2)\n",
    "delta = np.sqrt(deltachi2)\n",
    "print('PCONF = %.2f, deltachisq = %.2f, delta = %.2f' % (PCONF,deltachi2,delta))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# THESE NUMBERS DO NOT QUITE MATCH ASTER EQUATION 2.56\n",
    "print('%.1f%% parameter confidence intervals (m, delm)' % (PCONF*100))\n",
    "delm = delta*np.sqrt(np.diag(covm)).reshape((3,1))\n",
    "showmat(np.hstack((m, delm)),2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot_ellipse() can be found in lib_inverse.py\n",
    "def plot_ellipse_custom(inds,DELTA2,C,m):\n",
    "    indc = np.ix_(inds,inds)\n",
    "    C = covm[indc]\n",
    "    [x,y] = plot_ellipse(DELTA2,C,np.array([m[inds[0]],m[inds[1]]]))\n",
    "    plt.plot(x,y,'r')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10,5))\n",
    "\n",
    "plt.subplot(1,3,1)\n",
    "plot_ellipse_custom([0,1],deltachi2,C,m)\n",
    "plt.xlabel(r'$m_1$ (m)')\n",
    "plt.ylabel(r'$m_2$ (m/s)')\n",
    "\n",
    "plt.subplot(1,3,2)\n",
    "plot_ellipse_custom([0,2],deltachi2,C,m)\n",
    "plt.xlabel(r'$m_1$ (m)')\n",
    "plt.ylabel(r'$m_3$ (m/s$^2$)')\n",
    "\n",
    "plt.subplot(1,3,3)\n",
    "plot_ellipse_custom([1,2],deltachi2,C,m)\n",
    "plt.xlabel(r'$m_2$ (m/s)')\n",
    "plt.ylabel(r'$m_3$ (m/s$^2$)')\n",
    "\n",
    "plt.subplots_adjust(wspace=0.5)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "a4e78903152c5c8c0a1611fb4b22aa78baa308a7734308b75819bd49ff5f4ec4"
  },
  "kernelspec": {
   "display_name": "inverse [conda env:inverse]",
   "language": "python",
   "name": "conda-env-inverse-py"
  },
  "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.10.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}