{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# 7. Maximum Likelihood fit\n\nSince version ``0.6.12`` SciKit-GStat implements an utility function factory\nwhich takes a Variogram instance and builds up a (negative) maximum likelihood\nfunction for the associated sample, distance matrix and model type.\nThe used function is defined in eq. 14 from Lark (2000). Eq. 16 from same\npublication was adapted to all available theoretical models available in\nSciKit-GStat, with the exception of the harmonized model, which\ndoes not require a fitting.\n\nThis tutorial helps you to build your own maximum likelihood approximation\nfor use in your geostatistical application. Note that SciKit-GStat is build\naround the method of moments and does not support maximum likelihood beyond\nthe approach presented here.\n\n\n**References**\n\nLark, R. M. \"Estimating variograms of soil properties by the method\u2010of\u2010moments and maximum\nlikelihood.\" European Journal of Soil Science 51.4 (2000): 717-728.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import skgstat as skg\nfrom skgstat.util.likelihood import get_likelihood\nimport numpy as np\nimport matplotlib.pyplot as plt\nfrom scipy.optimize import minimize\nimport warnings\nfrom time import time\nwarnings.filterwarnings('ignore')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We use the pancake dataset, sampled at 300 random locations to produce a quite dense sample.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "c, v = skg.data.pancake(N=300, seed=42).get('sample')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "First of, the variogram is calculated. We use Scott's rule to determine\nthe number of lag classes, explicitly set Trust-Region Reflective as\nfitting method (although its default) and limit the distance matrix to\n70% of the maximum separating distance.\nAdditionally, we capture the processing time for the whole variogram\nestimation. Note, that this also includes the calculation of the\ndistance matrix, which is a mututal step.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "t1 = time()\nV = skg.Variogram(c,v, bin_func='scott', maxlag=0.7, fit_func='trf')\nt2 = time() # get time for full analysis, including fit\nprint(f\"Processing time: {round((t2 - t1) * 1000)} ms\")\nprint(V)\nfig = V.plot()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Maximum likelihood using SciKit-GStat\nFirst step to perform the fitting is to make initial guesses for the parameters.\nHere, we take the mean separating distance for the effective range, the sample\nvariance for the sill and 10% of the sample variance for the nugget.\nTo improve performance and runtime, we also define a boundary to restrict\nthe parameter space.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# base initial guess on separating distance and sample variance\nsep_mean = V.distance.mean()\nsam_var = V.values.var()\nprint(f\"Mean sep. distance:  {sep_mean.round(1)}    sample variance: {sam_var.round(1)}\")\n\n# create initial guess\n#    mean dist.  variance    5% of variance\np0 = np.array([sep_mean, sam_var, 0.1 * sam_var])\nprint('initial guess: ', p0.round(1))\n\n# create the bounds to restrict optimization\nbounds = [[0, V.bins[-1]], [0, 3*sam_var], [0, 2.9*sam_var]]\nprint('bounds:        ', bounds)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Next step is to pass the Variogram instance to the function factory.\nWe find optimal parameters by minimizing the returned negative\nlog-likelihood function. Please refer to SciPy's minimize function to learn\nabout attributes. The returned function from the utility suite is built with\nSciPy in mind, as the function signature complies to SciPy's interface and,\nthus can just be passed to the minimize function.\nHere, we pass the initial guess, the bounds and set the solver method to\nSLSQP, a suitable solver for bounded optimization.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# load the likelihood function for this variogram\nlikelihood = get_likelihood(V)\n\n# minimize the likelihood function \nt3 = time()\nres = minimize(likelihood, p0, bounds=bounds, method='SLSQP')\nt4 = time()\n\n# some priting\nprint(f\"Processing time {np.round(t4 - t3, 2)} seconds\")\nprint('initial guess:     ', p0.round(1))\nprint('optimal parameters:', res.x.round(1))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Here, you can see one of the main limitations for ML approaches: runtime.\nA sample size of 300 is rather small and the ML is running \nconsiderably slower that MoM.\n\nApply the optimized parameters. For comparison, the three method-of-moment methods\nfrom SciKit-GStat are applied as well. Note that the used sample is quite dense.\nThus we do not expect a different between the MoM based procedures.\nThey should all find the same paramters.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# use 100 steps\nx = np.linspace(0, V.bins[-1], 100)\n\n# apply the maximum likelihood fit parameters\ny_ml = V.model(x, *res.x)\n\n# apply the trf fit\ny_trf = V.fitted_model(x)\n\n# apply Levelberg marquard\nV.fit_method = 'lm'\ny_lm = V.fitted_model(x)\n\n# apply parameter ml\nV.fit_method = 'ml'\ny_pml = V.fitted_model(x)\n\n# check if the method-of-moment fits are different\nprint('Trf and Levenberg-Marquardt identical: ', all(y_lm - y_trf < 0.1))\nprint('Trf and parameter ML identical:        ', all(y_pml - y_trf < 0.1))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Make the result plot\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.plot(V.bins, V.experimental, '.b', label='experimental')\nplt.plot(x, y_ml, '-g', label='ML fit (Lark, 2000)')\nplt.plot(x, y_trf, '-b', label='SciKit-GStat TRF')\nplt.legend(loc='lower right')\n#plt.gcf().savefig('compare.pdf', dpi=300)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Build from scratch\nSciKit-GStat's utility suite does only implement the maximum likelihood\napproach as published by Lark (2000). There are no settings to adjust\nthe returned function, nor use other implementations. If you need to use\nanother approach, the idea behind the implementation is demonstrated below\nfor the spherical variogram model. This solution is only build on SciPy and\ndoes not need SciKit-GStat, in case the distance matrix is build externally.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from scipy.spatial.distance import squareform\nfrom scipy.linalg import inv, det\n\n# define the spherical model only dependent on the range\ndef f(h, a):\n    if h >= a:\n        return 1.\n    elif h == 0:\n        return 0.\n    return (3*h) / (2*a) - 0.5 * (h / a)**3\n\n# create the autocovariance matrix \ndef get_A(r, s, b, dists):\n    a = np.array([f(d, r) for d in dists])\n    A = squareform((s / (s + b)) * (1 - a))\n    np.fill_diagonal(A, 1)\n\n    return A\n\n# likelihood function\ndef like(r, s, b, z, dists):\n    A = get_A(r, s, b, dists)\n    n = len(A)\n    A_inv = inv(A)\n    ones = np.ones((n, 1))\n    z = z.reshape(n, -1)\n    m = inv(ones.T @ A_inv @ ones) @ (ones.T @ A_inv @ z)\n    b = np.log((z - m).T @ A_inv @ (z - m))\n    d = np.log(det(A))\n    if d == -np.inf:\n        print('invalid det(A)')\n        return np.inf\n    loglike = (n / 2)*np.log(2*np.pi) + (n / 2) - (n / 2)* np.log(n) + 0.5* d + (n / 2) * b\n    return loglike.flatten()[0]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "You can adjust the autocorrelation function above to any other model and\nimplement other approaches by adjusting the ``like`` function.\nFinally, minimizing this function is the same like before. You also have to \ntake care of the SciPy interface, as you need to wrap the custom likelihood\nfunction to provide the parameters in the right format.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from scipy.optimize import minimize\nfrom scipy.spatial.distance import pdist\n\n# c and v are coordinate and values array from the data source\nz = np.array(v)\n\n# in case you use 2D coordinates, without caching and euclidean metric, skgstat is using pdist under the hood\ndists = pdist(c)\n\nfun = lambda x, *args: like(x[0], x[1], x[2], z=z, dists=dists)\nt3 = time()\nres = minimize(fun, p0, bounds=bounds)\nt4 = time()\nprint(f\"Processing time {np.round(t4 - t3, 2)} seconds\")\nprint('initial guess:     ', p0.round(1))\nprint('optimal parameters:', res.x.round(1))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Finally, you can produce the same plot as before.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\nmod = lambda h: f(h, res.x[0]) * res.x[1] + res.x[2]\n\nx = np.linspace(0, 450, 100)\ny = list(map(mod, x))\ny2 = V.fitted_model(x)\n\nplt.plot(V.bins, V.experimental, '.b', label='experimental')\nplt.plot(x, y, '-g', label='ML fit (Lark, 2000)')\nplt.plot(x, y2, '-b', label='SciKit-GStat default fit')\nplt.legend(loc='lower right')\n#plt.gcf().savefig('compare.pdf', dpi=300)"
      ]
    }
  ],
  "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.9"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}