Metadata-Version: 2.1
Name: thetafit
Version: 0.0.1
Summary: Parameter estimation for nonlinear models
Home-page: https://github.com/solbes/thetafit
Author: Antti Solonen
Author-email: antti.solonen@gmail.com
License: MIT
Download-URL: https://github.com/solbes/thetafit/archive/refs/tags/0.0.1.tar.gz
Description: # thetafit
        
        A simple python tool for Bayesian parameter estimation of nonlinear models, inspired by the [mcmcstat MATLAB package](https://github.com/mjlaine/mcmcstat).
        
        Supports MAP estimation (optimization) and MCMC sampling via the Adaptive Metropolis (AM) algorithm.
        
        The user supplies the function for calculating ``-2*log(likelihood)``, which, for Gaussian likelihood, amounts to calculating the sum of squares misfit between data and model.
        
        A simple example is given below, for more examples see the [demos](https://github.com/solbes/thetafit/tree/main/demos) folder.
        
        ## Example: MONOD model
        
        Let us fit the parameters ``theta = [a,b]`` of a simple model ``y = ax/(b+x) + N(0,sig2)``, describing bacterial growth.
        
        ### 1) Define model and data
        
        First, we define the data and the sum of squares function. The function needs to be given in the form ``ssfun(theta,data)``, where ``theta`` is a dict of parameter name-value pairs, and ``data`` is an object of arbitrary form that contains typically the data and all other constants etc. that the function evaluation needs.
        
        ```python
        import numpy as np
        import pandas as pd
        from thetafit import Parameter, optimize, sample
        from thetafit.utils import McmcOptions, jacob
        import matplotlib.pyplot as plt
        
        data = {
            'x': np.array([28., 55., 83., 110., 138., 225., 375.]),
            'y': np.array([0.053, 0.060, 0.112, 0.105, 0.099, 0.122, 0.125]),
            'sig2': 1.0,
        }
        
        modelfun = lambda theta, x: theta['a']*x/(theta['b']+x)
        ssfun    = lambda theta, data: sum((data['y']-modelfun(theta, data['x']))**2/data['sig2'])
        ```
        
        ### 2) MAP estimation
        
        Then, let us first minimize the sum of squares to get the MAP estimate. The initial values (with optional bounds) for the parameters are given as a list of ``Parameter`` objects. Here, we use a helper function for constructing such a list from the ``theta`` dict. 
        
        The optimization is called as ``optimize(ssfun, data, pars, **kwargs)``, where ``**kwargs`` are optional optimization configs fed to ``scipy.optimize``, where one can choose the optimizer and its settings.
        
        ```python
        th0 = {'a': 0.15, 'b': 100}
        
        params = lambda theta: [
            Parameter('a', init=theta['a'], target=True),
            Parameter('b', init=theta['b'], target=True)
        ]
        
        thopt, res = optimize(ssfun, data, params(th0))
        print(thopt)
        ```
        
        which outputs the parameter values (and some information about the optimization process in ``res``):
        
        ```python
        {'a': 0.14541945253897567, 'b': 49.05274598150096}
        ```
        
        ### 3) MCMC
        
        Let us then do MCMC for the parameters, starting from the initial values obtained from the optimization. First, let us fix the error variance ``sig2`` to the mean squared error (MSE) estimate:
        
        ```python
        mse = res.fun/(len(data['y'])-2)
        data['sig2'] = mse
        ```
        
        Then, we calculate the initial proposal covariance matrix via the Gauss-Newton Hessian approximation at the MAP estimate (note that a diagonal initial proposal covariance could be used as well, the adaptation will tune the proposal as the sampling proceeds).
        
        ```python
        J = jacob(modelfun, data['x'], params(thopt))
        qcov = np.linalg.inv(J.T.dot(J))*mse
        ```
        
        Then we define the options for the MCMC run via the ``McmcOptions`` object and perform the sampling:
        
        ```python
        options = McmcOptions(nsimu=4000, qcov=qcov, adaptint=100, printint=500)
        res, chain, sschain = sample(ssfun, data, params(thopt), options)
        ```
        
        which outputs the progress of the sampling:
        
        ```
        Sampling these parameters:
        name	start	[min,max]
        a	0.14541945253897567	[-inf,inf]
        b	49.05274598150096	[-inf,inf]
        i: 500, rejected: 60.8%, out of bounds: 0.0%
        i: 1000, rejected: 62.6%, out of bounds: 0.0%
        i: 1500, rejected: 63.8%, out of bounds: 0.0%
        i: 2000, rejected: 64.3%, out of bounds: 0.0%
        i: 2500, rejected: 64.9%, out of bounds: 0.0%
        i: 3000, rejected: 65.5%, out of bounds: 0.0%
        i: 3500, rejected: 65.8%, out of bounds: 0.0%
        ```
        
        Then, let us visualize the results:
        
        ```python
        pd.plotting.scatter_matrix(chain)
        ```
        ![thetafit_params](https://user-images.githubusercontent.com/6495497/146729259-07466ecf-84ef-4ddf-9d9c-ac7e205cd0b1.png)
        
        ```python
        xplot = np.linspace(0,400)
        preds = np.stack([modelfun(chain.loc[i], xplot) for i in range(2000,4000,10)])
        
        plt.figure(figsize=(7,4))
        plt.plot(xplot, preds.T, 'r-', lw=1, alpha=0.2)
        plt.plot(data['x'], data['y'], 'kx')
        plt.grid(True)
        plt.show()
        ```
        
        ![thetafit_preds](https://user-images.githubusercontent.com/6495497/146729611-5ffb38cc-676e-491f-a3aa-d577b7b894e0.png)
        
Keywords: bayes,nonlinear,parameter estimation
Platform: UNKNOWN
Classifier: Development Status :: 3 - Alpha
Classifier: Intended Audience :: Science/Research
Classifier: Topic :: Scientific/Engineering :: Artificial Intelligence
Classifier: License :: OSI Approved :: MIT License
Classifier: Programming Language :: Python :: 3
Description-Content-Type: text/markdown
Provides-Extra: dev
