INFR.collocation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One benefit of modelling dependent errors is that the waiting time for physical calibration and validation of a new satellite can be reduced considerably when two instead of three collocations are employed (Su et al. 2014).  Here, we demonstrate the retrieval of regression slope using two collocation methods: standard triple collocation (Stoffelen 1998) and a collocation of only two datasets, but where one dataset offers five samples (the other four samples being \"nearby\" in time or space).  We show that both methods (triple collocation and a so-called INFERS method) provide good estimates of regression slope.\n",
    "\n",
    "Triple collocation employs this error model:\n",
    "\n",
    "\\begin{eqnarray}\n",
    "  \\begin{array}{r} \\mathrm{in\\ situ}\\ \\\\\n",
    "                   \\mathrm{nowcast} \\\\\n",
    "                   \\mathrm{satellite} \\end{array}\n",
    "  \\begin{array}{r}   I   \\\\   N   \\\\   S   \\end{array}\n",
    "  \\begin{array}{c} \\ = \\ \\\\ \\ = \\ \\\\ \\ = \\ \\end{array}\n",
    "  \\begin{array}{l} \\color{white}{\\alpha_N + \\beta_N} \\color{black} t + \\epsilon_I \\\\\n",
    "                   \\alpha_N + \\beta_N t + \\epsilon_N \\\\\n",
    "                   \\alpha_S + \\beta_S t + \\epsilon_S \\end{array} \\\\\n",
    "  \\nonumber\n",
    "\\end{eqnarray}\n",
    "\n",
    "where $t$ is truth and $\\epsilon_I \\epsilon_N \\epsilon_S$ are independent errors.  The regression slope of interest is $\\beta_N$ (or $\\beta_S$).  Following McColl et al. (2014), we can retrieve $\\beta_N$ from the $INS$ collocations.  In [julia](https://julialang.org) this is"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "estimated bn  0.398 should be close to  0.400\n",
      "estimated bn  0.599 should be close to  0.600\n",
      "estimated bn  0.801 should be close to  0.800\n",
      "estimated bn  0.997 should be close to  1.000\n",
      "estimated bn  1.193 should be close to  1.200\n",
      "estimated bn  1.396 should be close to  1.400\n",
      "estimated bn  1.595 should be close to  1.600\n",
      "estimated bn  1.792 should be close to  1.800\n",
      "estimated bn  1.998 should be close to  2.000\n",
      "estimated bn  2.202 should be close to  2.200\n",
      "estimated bn  2.403 should be close to  2.400\n",
      "estimated bn  2.594 should be close to  2.600\n",
      "estimated bn  2.815 should be close to  2.800\n",
      "estimated bn  2.993 should be close to  3.000\n",
      "estimated bn  3.192 should be close to  3.200\n"
     ]
    }
   ],
   "source": [
    "ct = 9.0 ;            ci = 2.0;                         # rand() is random between 0 and 1\n",
    "an = 1.0 ;            cn = 1.0;                         # randn() is Gaussian with 0 mean and SD of 1\n",
    "as = 2.0 ; bs = 1.5 ; cs = 1.5;\n",
    "numb = 10^5;                                            # number of collocations\n",
    "\n",
    "for bn = 0.4:0.2:3.2\n",
    "  TT = ct *  rand(numb);                                # truth\n",
    "  DI = ci * randn(numb);                                # in situ error\n",
    "  DN = cn * randn(numb);                                # nowcast error\n",
    "  DS = cs * randn(numb);                                # satellite error\n",
    "  II =           TT + DI;\n",
    "  NN = an + bn * TT + DN;\n",
    "  SS = as + bs * TT + DS;\n",
    "  @printf(\"estimated bn %6.3f should be close to %6.3f\\n\", cov(NN,SS)/cov(II,SS), bn)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The INFERS method employs a more complicated model, but with the simplest (AR-1) type of shared or propagated error.  Note that with application of AR-1, the minimum number of equations (and samples of the second dataset) to match the number of unknowns is four (i.e., identification is possible with NFER or NFRS), but the symmetry of five samples (NFERS) facilitates a retrieval of all parameters, which is done numerically:\n",
    "\n",
    "\\begin{eqnarray} \\\\\n",
    "  \\begin{array}{r} \\mathrm{in\\ situ}\\ \\\\\n",
    "                   \\mathrm{nowcast} \\\\\n",
    "                   \\mathrm{forecast} \\\\\n",
    "                   \\mathrm{extended\\ forecast} \\\\\n",
    "                   \\mathrm{revcast} \\\\\n",
    "                   \\mathrm{extended\\ revcast} \\end{array}\n",
    "  \\begin{array}{r} I \\\\ N \\\\ F \\\\ E \\\\ R \\\\ S \\end{array}\n",
    "  \\begin{array}{c} = \\\\  =  \\\\ = \\\\ = \\\\ = \\\\ = \\end{array}\n",
    "  \\begin{array}{l} \\color{white}{\\alpha_N + \\beta_N} \\color{black} t + \\color{white}{\\lambda_E (  \\lambda_F ( \\lambda_N} \\color{black} \\epsilon_I \\\\\n",
    "                                 \\alpha_N + \\beta_N  t + \\color{white}{\\lambda_E (                \\lambda_F (} \\color{black} \\lambda_N \\epsilon_I + \\epsilon_N \\\\\n",
    "                                 \\alpha_F + \\beta_F  t + \\color{white}{\\lambda_E (} \\color{black} \\lambda_F ( \\lambda_N \\epsilon_I + \\epsilon_N ) + \\epsilon_F \\\\\n",
    "                                 \\alpha_E + \\beta_E  t +               \\lambda_E (                \\lambda_F ( \\lambda_N \\epsilon_I + \\epsilon_N ) + \\epsilon_F ) + \\epsilon_E \\\\\n",
    "                                 \\alpha_R + \\beta_R  t + \\color{white}{\\lambda_E (} \\color{black} \\lambda_R ( \\lambda_N \\epsilon_I + \\epsilon_N ) + \\epsilon_R \\\\\n",
    "                                 \\alpha_S + \\beta_S  t +               \\lambda_S (                \\lambda_R ( \\lambda_N \\epsilon_I + \\epsilon_N ) + \\epsilon_R ) + \\epsilon_S \\end{array} \\\\\n",
    "  \\nonumber\n",
    "\\end{eqnarray}\n",
    "\n",
    "We thus allow for correlated errors, possibly because the in situ data is assimilated by the analysis ($\\lambda_N \\epsilon_I$) or because there is common error in adjacent analysis windows (i.e., taking \"persistence\" as a reasonable forecast or revcast, shared error is accommodated by nonzero $\\lambda_F$ and $\\lambda_R$).  Retrieval of $\\beta_N$ employs a cost function defined by the covariance of FR and ES.  In particular, the method requires error correlation between E and S (cf. discussion of instrumental variable regression in Su et al. 2014).  This is in stark contrast to the triple collocation method, where error correlation is to be avoided; here, error correlation is not only helpful, it is required!  In practice, analyses admit observations in windows that are typically as large or larger than the interval between analyses (note that fractional ES error covariance is also quantified by the product $\\lambda_E \\lambda_F \\lambda_R \\lambda_S$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "estimated bn  0.401 should be close to  0.400 and st  6.762 should be close to  6.750\n",
      "estimated bn  0.598 should be close to  0.600 and st  6.781 should be close to  6.750\n",
      "estimated bn  0.751 should be close to  0.800 and st  6.651 should be close to  6.750\n",
      "estimated bn  1.001 should be close to  1.000 and st  6.786 should be close to  6.750\n",
      "estimated bn  1.190 should be close to  1.200 and st  6.706 should be close to  6.750\n",
      "estimated bn  1.403 should be close to  1.400 and st  6.778 should be close to  6.750\n",
      "estimated bn  1.603 should be close to  1.600 and st  6.728 should be close to  6.750\n",
      "estimated bn  1.801 should be close to  1.800 and st  6.744 should be close to  6.750\n",
      "estimated bn  2.005 should be close to  2.000 and st  6.753 should be close to  6.750\n",
      "estimated bn  2.198 should be close to  2.200 and st  6.706 should be close to  6.750\n",
      "estimated bn  2.388 should be close to  2.400 and st  6.735 should be close to  6.750\n",
      "estimated bn  2.596 should be close to  2.600 and st  6.763 should be close to  6.750\n",
      "estimated bn  2.817 should be close to  2.800 and st  6.814 should be close to  6.750\n",
      "estimated bn  3.015 should be close to  3.000 and st  6.748 should be close to  6.750\n",
      "estimated bn  3.205 should be close to  3.200 and st  6.703 should be close to  6.750\n"
     ]
    }
   ],
   "source": [
    "TI = 9.0 ;                      ; CI = 1.0              # rand() is random between 0 and 1\n",
    "AN = 1.0 ;            LN = 0.20 ; CN = 1.5              # randn() is Gaussian with 0 mean and SD of 1\n",
    "AF = 2.0 ; BF = 1.5 ; LF = 0.85 ; CF = 0.4\n",
    "AE = 2.4 ; BE = 1.0 ; LE = 1.20 ; CE = 0.5\n",
    "AR = 3.0 ; BR = 0.8 ; LR = 0.95 ; CR = 0.6\n",
    "AS = 0.0 ; BS = 0.5 ; LS = 1.06 ; CS = 0.8\n",
    "numb = 10^5                                             # number of collocations\n",
    "\n",
    "for BN = 0.4:0.2:3.2\n",
    "  TT = TI *  rand(numb)                                 # truth\n",
    "  EI = CI * randn(numb)                                 # in situ error\n",
    "  EN = CN * randn(numb)                                 # nowcast error\n",
    "  EF = CF * randn(numb)                                 # forecast error\n",
    "  EE = CE * randn(numb)                                 # extended forecast error\n",
    "  ER = CR * randn(numb)                                 # revcast error\n",
    "  ES = CS * randn(numb)                                 # extended  revcast error\n",
    "  ii =           TT +                  EI\n",
    "  nn = AN + BN * TT +             LN * EI + EN\n",
    "  ff = AF + BF * TT +       LF * (LN * EI + EN) + EF\n",
    "  ee = AE + BE * TT + LE * (LF * (LN * EI + EN) + EF) + EE\n",
    "  rr = AR + BR * TT +       LR * (LN * EI + EN) + ER\n",
    "  ss = AS + BS * TT + LS * (LR * (LN * EI + EN) + ER) + ES\n",
    "\n",
    "  vari = cov(ii, ii)\n",
    "  varn = cov(nn, nn)\n",
    "  varf = cov(ff, ff)\n",
    "  vare = cov(ee, ee)\n",
    "  varr = cov(rr, rr)\n",
    "  vars = cov(ss, ss)\n",
    "  cvin = cov(ii, nn)\n",
    "  cvif = cov(ii, ff)\n",
    "  cvie = cov(ii, ee)\n",
    "  cvir = cov(ii, rr)\n",
    "  cvis = cov(ii, ss)\n",
    "  cvnf = cov(nn, ff)\n",
    "  cvne = cov(nn, ee)\n",
    "  cvnr = cov(nn, rr)\n",
    "  cvns = cov(nn, ss)\n",
    "  cvfe = cov(ff, ee)\n",
    "  cvfr = cov(ff, rr)\n",
    "  cvfs = cov(ff, ss)\n",
    "  cver = cov(ee, rr)\n",
    "  cves = cov(ee, ss)\n",
    "  cvrs = cov(rr, ss)\n",
    "\n",
    "  function cost(st::Float64, bn::Float64)\n",
    "    ibt = vari -        st\n",
    "    nbt = varn - bn^2 * st ; inbt = cvin - bn * st\n",
    "    tmp = st * (nbt - bn * inbt)\n",
    "    bf = (cvif * nbt - cvnf * inbt) / tmp ; lf = (cvnf - bf * bn * st) /       nbt\n",
    "    be = (cvie * nbt - cvne * inbt) / tmp ; le = (cvne - be * bn * st) / (lf * nbt)\n",
    "    br = (cvir * nbt - cvnr * inbt) / tmp ; lr = (cvnr - br * bn * st) /       nbt\n",
    "    bs = (cvis * nbt - cvns * inbt) / tmp ; ls = (cvns - bs * bn * st) / (lr * nbt)\n",
    "    fbt = varf - bf^2 * st ; ifbt = cvif - bf * st\n",
    "    ebt = vare - be^2 * st ; iebt = cvie - be * st\n",
    "    rbt = varr - br^2 * st ; irbt = cvir - br * st\n",
    "    sbt = vars - bs^2 * st ; isbt = cvis - bs * st\n",
    "    mask = 0.0 ; (ibt <= 0 || nbt <= 0 || inbt <= 0 || fbt <= 0 || ifbt <= 0 || ebt <= 0 || iebt <= 0 || rbt <= 0 || irbt <= 0 || sbt <= 0 || isbt <= 0) && (mask = 1.0)\n",
    "    fecov = cvfe - bf * be * st - le *                fbt\n",
    "    frcov = cvfr - bf * br * st -      lf *      lr * nbt\n",
    "    fscov = cvfs - bf * bs * st -      lf * ls * lr * nbt\n",
    "    ercov = cver - be * br * st - le * lf *      lr * nbt\n",
    "    escov = cves - be * bs * st - le * lf * ls * lr * nbt\n",
    "    rscov = cvrs - br * bs * st -           ls *      rbt\n",
    "    log(abs(fecov)) + log(abs(frcov)) + log(abs(fscov)) + log(abs(ercov)) + log(abs(escov)) + log(abs(rscov)) + 99 * mask\n",
    "  end\n",
    "\n",
    "  minst = minimum([cvin^2 / varn, cvif^2 / varf, cvir^2 / varr])\n",
    "  maxst = vari\n",
    "  intst = collect(linspace(minst, maxst, 1000))\n",
    "  minbn = cvin / vari                                   # check cost for values between OLS and RLS\n",
    "  maxbn = varn / cvin\n",
    "  intbn = collect(linspace(minbn, maxbn, 1000))\n",
    "  costs = Array(Float64, 1000, 1000)\n",
    "  for (a, vala) in enumerate(intst), (b, valb) in enumerate(intbn)\n",
    "    costs[a,b] = cost(vala, valb) ; (isinf(costs[a,b]) && costs[a,b] < 0) && (costs[a,b] = Inf)\n",
    "  end\n",
    "  indst, indbn = ind2sub(costs, findmin(costs)[2])\n",
    "  st = intst[indst] ; bn = intbn[indbn]\n",
    "  @printf(\"estimated bn %6.3f should be close to %6.3f\", bn, BN)\n",
    "  @printf(     \" and st %6.3f should be close to %6.3f\\n\", st, TI^2/12)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Raw Cell Format",
  "kernelspec": {
   "display_name": "Julia 0.4.6-pre",
   "language": "julia",
   "name": "julia-0.4"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.4.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}