---
# System prepended metadata

title: 'Parameters finding using Maximum Likelihood method '
tags: ['#modmesh', '#minimization', '#likelihood', '#optimization']

---

# Parameter finding prospective
(Modmesh) parameter finding: an application estimates the parameters of the defined functions on some numerical data. We use maximum likelihood method as an estimation, and finding maximum/minimum is an optimization problem.

Log-maximum likelihood method can be represented as minimized of the form 

$F(\pmb{\theta})=-\sum_{k=1}^{K}\ln f(x=x_k;\pmb{\theta})$
where $x_k$ is a data sample, and $\pmb\theta$ is parameter vector. 


## Code information

Parameters finding is part of modmesh. Parts of gradient calculator and iterator/stopping will be impleted by c++. Part of plot will be implemented by pyside6.
The source code of modmesh is here. 
https://github.com/solvcon/modmesh

1. Drawing Plots of results
2. Iterator and stopping for minimization
3. Gradient calculator

## Problem to solve
Optimization problem have many solutions and there is no universal solution can handle all situations. In our prospestive, gradient methods are adopted as they can cover many defined functions and expand the dimension of data. 

As we used the defined PDFs, their derivatives (first order) is also analytically defined. For example, radioactive decay follows an exponential decay, that is $N(t)=ae^{-\lambda t}$. Its derivative value is $\frac{dN(t)}{d\tau}=-at e^{\lambda t}$.


Some gradient methods, such as Newton method and conjugate directions, will use second derivative matrix $G$ (hessian matrix) and its inverse matrix $V$ (covariance matrix). 

Here are some gradient methods which are used in scipy.optimize.minimize 
* Unconstrained minimization
    * BFGS: quasi-Newton method of Broyden, Fletcher, Goldfarb, and Shanno (BFGS). It uses the first derivatives only. It is the default method. Ref. 3 p.136
    * CG: a nonlinear conjugate gradient algorithm. Ref 3 p.120-122
    * Newton-CG: the truncated Newton method. Suitable for large-scale problems. Ref. 3 p.168
    * trust-ncg: the Newton conjugate gradient trust-region algorithm for unconstrained minimization. This algorithm requires the gradient and either the Hessian or a function that computes the product of the Hessian with a given vector. 

* Bound-constrained minimization
    * L-BFGS-B: the L-BFGS-B algorithm. Default method. ref. [Fortran 77 code](https://dl.acm.org/action/downloadSupplement?doi=10.1145%2F279232.279236&file=778.gz), [paper](https://dl.acm.org/doi/10.1145/279232.279236). I saw someone wrap the fortran code with c++ as [here](https://github.com/colinfang/lbfgsb_cpp/tree/master).
    * Powell: a modification of Powell’s method, which is a conjugate direction method. Ref. 4 Section 10.7. Direction Set (Powell's) Methods in Multidimensions

* MIGRAD: Its strategy to find the local minimum is simply to vary the set of parameters $p$, by small (variable-sized) steps, in a direction which causes $F$ to decrease until it finds the point $\hat{p}$ from which $F$ increases in all allowed directions. Although not needed, if the numerical values of the derivative $\frac{\partial F(p)}{\partial p}$ at any point $p$ are known, they can be provided to MIGRAD to help in the minimization. [from Belle 2 document](https://software.belle2.org/development/sphinx/zzz-fitting/doc/minuit/doc/background.html)


## Prospective Users

Stastical 
 
## System Architecture

c++ in the bottom 
Python API

For the coding side, I have little experiences about programming architecture. 

Some objects need to implete as below
* Defined functions:
    * Gaussian function
    * Exponential function
    * Polynomial
    * Chebyshev polynomial
    * Poisson function
    * Crystal ball function:  The Crystal ball function is a Gaussian with a tail on the low side (or both side).
    * User defined function: it should be the same with scipy's user defined function
    * Stack up of above functions (or Call "ADD")
    * 
* First derivative vector of defined functions
* Hessian matrix of defined functions
* Constrains (parameters bound): for example $0<\mu<20.0$ and $0.5<\sigma<3.0$ for Gaussian function
* Log-likelihood value 
* Minimization methods
    * Iterator
    * Stop 
    * Linear programming for constrained minimization
* Plots


However, I found a good project based on MINUIT and provide a python API as *iminuit* [homepage](https://iminuit.readthedocs.io/en/stable/index.html).
## Extension
* Chi-square methods


## Reference
1. [scipy.optimize.minimize](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html#scipy.optimize.minimize): Minimization of scalar function of one or more variables.
2. [MINUIT Tutorial](https://root.cern.ch/download/minuit.pdf): MINUIT is conceived as a tool to find the minimum value of a multi-parameter function and analyze the shape of the function around the minimum. 
3. Numerical Optimization: book, Jorge Nocedal Stephen J. Wright, [pdf online](https://www.ime.unicamp.br/~pulino/MT404/TextosOnline/NocedalJ.pdf)
4. Numerical Recipes: book, William H. Press, Saul A. Teukolsky, William T. Vetterling, Brian P. Flannery [online book](http://numerical.recipes/book.html)