{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Writing your own fits\n", "``bayesmsd`` can fit any parametric shape of MSD curves; however, utilizing that full flexibility requires a little bit of coding. This tutorial will walk you through the implementation of a new fit model, ``TwoLocusRouseFit``. Note that a production version of this class is available in ``bayesmsd.lib``; after working through the—somewhat streamlined—implementation in this example, you can check its source code to see the last finishing touches.\n", "\n", "For this example we will implement a fit to the following functional expression:\n", "\\begin{equation}\\label{eq:MSD}\n", "\\text{MSD}(\\Delta t) = 2\\sigma^2 + 2\\Gamma\\sqrt{\\Delta t}\\left(1-\\text{e}^{-\\frac{\\tau}{\\Delta t}}\\right) + 2J\\,\\text{erfc}\\,\\sqrt{\\frac{\\tau}{\\Delta t}}\\,,\\tag{1}\n", "\\end{equation}\n", "where $\\tau\\equiv\\frac{1}{\\pi}\\left(\\frac{J}{\\Gamma}\\right)^2$.\n", "\n", "If $x_1(t)$ and $x_2(t)$ are the spatial positions of two loci at fixed positions along a polymer, the above expression is the MSD we should expect to see for the relative position $y(t)\\equiv x_2(t) - x_1(t)$, assuming the polymer dynamics follow the Rouse model.\n", "\n", "So how do we fit something like this? ``bayesmsd`` provides the whole fitting machinery in form of an abstract base class ``Fit``, which we can subclass to implement a new fit model. Upon doing so, our main task is to override the method ``params2msdm()``, which gives the functional form of the MSD dependent on the parameters. Of course we also have to specify which parameters we need and how they relate to each other. Let's get started!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Class definition and constructor\n", "\n", "```py\n", "import numpy as np\n", "from scipy.special import erfc\n", "\n", "import bayesmsd\n", "from noctiluca.analysis import MSD\n", "\n", "class TwoLocusRouseFit(bayesmsd.Fit):\n", " def __init__(self, data, motion_blur_f=0):\n", " super().__init__(data)\n", " self.motion_blur_f = motion_blur_f\n", "...\n", "```\n", "So far, not much happened; we called the constructor of ``Fit``, which will take care of input handling, and also populates the attributes\n", "\n", " + ``self.d``: number of spatial dimensions in the dataset,\n", " + ``self.T``: length of the longest trajectory.\n", " \n", "The ``motion_blur_f`` argument will be used below to take into account finite imaging exposure times; for now we just carry it along.\n", "\n", "Now we have to specify some details about this fit. Specifically: set ``self.ss_order``, populate ``self.parameters``, and potentially ``self.constraints``.\n", "```py\n", "... self.ss_order = 0\n", "```\n", "With the MSD given above, we are always in a steady state of order 0; for other fitting schemes, this attribute might be determined from the input.\n", "```py\n", "...\n", " self.parameters = {\n", " 'log(σ²)' : bayesmsd.Parameter(bounds=(-np.inf, np.inf)),\n", " 'log(Γ)' : bayesmsd.Parameter(bounds=(-np.inf, np.inf)),\n", " 'log(J)' : bayesmsd.Parameter(bounds=(-np.inf, np.inf)),\n", " 'log(τ)' : bayesmsd.Parameter(bounds=(-np.inf, np.inf)),\n", " }\n", "...\n", "```\n", "We initialize all the parameters that appear in the expression \\eqref{eq:MSD} and provide their domains; since we work in log-space, all parameters have the whole real line as domain; so we just set the bounds to infinity.\n", "\n", "Having defined the parameters, we can now fix any pre-determined relationships between them:\n", "```py\n", "...\n", " def fix_tau(params): \n", " return 2*(params['log(J)']-params['log(Γ)']) - np.log(np.pi)\n", " self.parameters['log(τ)'].fix_to = fix_tau\n", "...\n", "```\n", "Note that in this example we could also have omitted $\\tau$ from the parameter list and just hardcoded this identity, since $\\tau$ should never be independent. Sticking with this implementation, we would at least move ``fix_tau()`` to be a static method, instead of defining it at runtime (c.f. summary below).\n", "\n", "Finally, we can list any inequality constraints that the fit parameters should satisfy by default.\n", "```py\n", "... self.constraints = []\n", "```\n", "For the Rouse MSD in this example there are no additional constraints to be satisfied. Note that if you do not reset it, ``self.constraints`` defaults to ``[self.constraint_Cpositive]``. This constraint function is provided by the ``Fit`` baseclass and ensures positive definiteness of the covariance matrix. It is costly to evaluate, so if positive definiteness is guaranteed by the shape of the MSD (as it is here) you can omit it. Other notable examples:\n", "\n", " + a powerlaw MSD with $\\alpha\\in[0, 2)$ is positive definite\n", " + for a spline interpolation we have no way of ensuring positive definiteness, so use ``constraint_Cpositive``" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## MSD definition\n", "Onwards to the key part: how to calculate an MSD from the parameters we set up above.\n", "```py\n", "...\n", " def params2msdm(self, params):\n", " sigma2 = np.exp(params['log(σ²)'])\n", " Gamma = np.exp(params['log(Γ)'])\n", " J = np.exp(params['log(J)'])\n", " tau = np.exp(params['log(τ)'])\n", " \n", " @bayesmsd.deco.MSDfun\n", " @bayesmsd.deco.imaging(noise2=sigma2,\n", " f=self.motion_blur_f,\n", " alpha0=0.5,\n", " )\n", " def msd(dt, G=G, J=J, tau=tau):\n", " return ( 2*G*np.sqrt(dt)*(1-np.exp(-tau/dt))\n", " + 2*J*erfc(np.sqrt(tau/dt)) )\n", " \n", " return self.d*[(msd, 0)]\n", "...\n", "```\n", "That's already the full definition of this function. After assigning local variable names to the parameters (mainly for readability) we define the MSD function and then return it. Let's take a look at some details, starting from the end:\n", "\n", " + ``params2msdm`` should return a separate MSD function for each of the ``self.d`` dimensions in the data set. In this simplified implementation of ``TwoLocusRouseFit`` we assume that all dimensions are identical, so we can just copy the same values to all dimensions.\n", " + Beyond the MSD, ``params2msdm`` should also return the first moment for the Gaussian process, which corresponds to a constant offset (for ``ss_order = 0``, i.e. stationary processes) or a drift term (for ``ss_order = 1``, i.e. increment stationary processes). For the most part this should just be zero; this is why we return ``(msd, 0)`` tuples instead of just the ``msd`` functions for each dimension.\n", " + In defining the ``msd`` function\n", " ```py\n", " def msd(dt, G=G, J=J, tau=tau)\n", " ```\n", " we pass all the parameters as default arguments. This is important to get the scoping right; default arguments are evaluated at definition time, whereas references to external variables (i.e. if we did not use this default argument construction) are evaluated at execution time. By then the values might have changed from what we wanted to use in the definition.\n", " + The ``bayesmsd.deco.imaging()`` decorator adds imaging artifacts to the MSD, specifically motion blur and localization error. Motion blur is determined from the fractional exposure $f\\in[0, 1]$—this is the ratio of shutter time over frame time, i.e. $f=0$ is ideal stroboscopic illumination, $f=1$ is continuous illumination—and the effective short time MSD scaling $\\alpha_0$. The latter is just the log-slope of the \"true\" MSD (no motion blur or localization error) at a time lag of 1 frame, which for our Rouse model example is always 0.5.\n", " + Finally, the ``bayesmsd.deco.MSDfun`` decorator extends the validity of the time lag ``dt`` to the real line. Note how we did not have to take any precautions against division by 0 or square roots of negative numbers in defining the ``msd()`` function. Indeed, the ``MSDfun`` decorator guarantees that the ``dt`` argument is always a numpy array with strictly positive values: anything else will be patched using the analytical identities $\\text{MSD}(-\\Delta t)\\equiv \\text{MSD}(\\Delta t)$ and $\\text{MSD}(\\Delta t) = 0$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Initial values\n", "A ``Fit`` should be able to come up with a rough initial guess for parameters from the data. This can be more or less sophisticated, but of course the better the initial values, the faster the fit will converge, so it's always good to give a best effort estimate here. We will resort to eye-balling the empirical MSD (exactly the methodology that ``bayesmsd`` is set to improve upon).\n", "```py\n", "...\n", " def initial_params(self):\n", " e_msd = MSD(self.data) / self.d\n", " \n", " J = np.nanmean(np.concatenate([traj[:]**2 for traj in data],\n", " axis=0,\n", " ))\n", " G = np.nanmean(e_msd[1:5]/np.sqrt(np.arange(1, 5)))\n", " sigma2 = e_msd[1]/2\n", " \n", " return {\n", " 'log(σ²)' : np.log(sigma2),\n", " 'log(Γ)' : np.log(G),\n", " 'log(J)' : np.log(J),\n", " }\n", "...\n", "```\n", "Note that this implementation is not safe against any of the first 5 MSD points missing. While this is fine in most realistic situations, the library implementation adds a few additional checks here." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Prior\n", "``bayesmsd`` being a Bayesian method, we have to specify a prior over the parameter space. This does not have to be a normalized probability distribution, i.e. can be an improper prior. By default, ``Fit`` assumes a uniform prior over the specified parameter space:\n", "```py\n", "...\n", " def logprior(self, params):\n", " return 0\n", "...\n", "```\n", "Since we are using log parameters, this translates to a $\\frac{1}{x}$ prior over $\\sigma^2$, $\\Gamma$, and $J$, which is appropriate for (*a priori*) scale free, positive parameters. Often, it is more convenient to transform the parameter space and use a uniform prior, than putting some complicated prior on the parameters." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Summary: full implementation\n", "Collecting all the bits from the previous paragraphs, here is the full implementation of ``TwoLocusRouseFit``. Note that the library implementation in [bayesmsd.lib.TwoLocusRouseFit](../bayesmsd.rst#bayesmsd.lib.TwoLocusRouseFit) adds a few bells and whistles to this, most notably dimensionally independent parameters. Still, the below should provide a good baseline for implementing your own fits." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from scipy.special import erfc\n", "\n", "import bayesmsd\n", "from noctiluca.analysis import MSD\n", "\n", "class TwoLocusRouseFit(bayesmsd.Fit):\n", " def __init__(self, data, motion_blur_f=0):\n", " super().__init__(data)\n", " self.motion_blur_f = motion_blur_f\n", " \n", " self.ss_order = 0\n", " \n", " self.parameters = {\n", " 'log(σ²)' : bayesmsd.Parameter(bounds=(-np.inf, np.inf)),\n", " 'log(Γ)' : bayesmsd.Parameter(bounds=(-np.inf, np.inf)),\n", " 'log(J)' : bayesmsd.Parameter(bounds=(-np.inf, np.inf)),\n", " 'log(τ)' : bayesmsd.Parameter(bounds=(-np.inf, np.inf)),\n", " }\n", " \n", " self.parameters['log(τ)'].fix_to = TwoLocusRouseFit.fix_tau\n", " \n", " self.constraints = []\n", " \n", " @staticmethod\n", " def fix_tau(params): \n", " return 2*(params['log(J)']-params['log(Γ)']) - np.log(np.pi)\n", " \n", " def params2msdm(self, params):\n", " sigma2 = np.exp(params['log(σ²)'])\n", " Gamma = np.exp(params['log(Γ)'])\n", " J = np.exp(params['log(J)'])\n", " tau = np.exp(params['log(τ)'])\n", " \n", " @bayesmsd.deco.MSDfun\n", " @bayesmsd.deco.imaging(noise2=sigma2,\n", " f=self.motion_blur_f,\n", " alpha0=0.5,\n", " )\n", " def msd(dt, G=G, J=J, tau=tau):\n", " return ( 2*G*np.sqrt(dt)*(1-np.exp(-tau/dt))\n", " + 2*J*erfc(np.sqrt(tau/dt)) )\n", " \n", " return self.d*[(msd, 0)]\n", " \n", " def initial_params(self):\n", " e_msd = MSD(self.data) / self.d\n", "\n", " J = np.nanmean(np.concatenate([traj[:]**2 for traj in data],\n", " axis=0,\n", " ))\n", " G = np.nanmean(e_msd[1:5]/np.sqrt(np.arange(1, 5)))\n", " sigma2 = e_msd[1]/2\n", "\n", " return {\n", " 'log(σ²)' : np.log(sigma2),\n", " 'log(Γ)' : np.log(G),\n", " 'log(J)' : np.log(J),\n", " }\n", " \n", " # logprior is uniform by default, so can just omit it" ] } ], "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.3" } }, "nbformat": 4, "nbformat_minor": 2 }