"This program demonstrates the different ways to optimize a system with 3 unknowns given a set of observations\n",
"and prior information on the PDF of the unknowns. There are currently three options:\n",
"\n",
"(1) a maximum likelihood estimate: this estimate is a simple linear algebra solution to an overconstrained problem. It does not use\n",
" prior information on parameters or their uncertainty as it is assumed that a unique solution exists. It does return an estimate of\n",
" the uncertainty in the optimal parameters. This method is an alternative to the type of minimization routines that come standard with\n",
" many software packages such as IDL, MATLAB, python, etc.\n",
"\n",
"(2) a full Kalman Filter estimate: this method is also suitable when there are more unknowns that observations, or when the solution is\n",
" not uniquely defined. In that case, a-priori information on the probability density (PDF) of the solution is used. The method returns an\n",
" optimal parameter set and its full covariance.\n",
"\n",
"(3) an ensemble Kalman Filter estimate: the Kalman filter equations are again solved, but now using the statistics of an ensemble to\n",
" find the optimum. This solution is in pri nciple inferior to the one above (statistical errors) but sometimes easier to calculate when the structure of the underlying model is complex. The method returns an ensemble of solutions to the problem, from which a mean and covariance\n",
" can be derived that is the same as in the full Kalman filter.\n",
"\n",
"The problem that we are optimizing here is simple. We have an equation of the form:\n",
"\n",
" y = a + b*x + c* x^2\n",
"\n",
"and the question is whether we can find a set [a,b,c] given imperfect observations (y). These imperfect observations are created by taking a\n",
"'perfect' y (outcome of the equation given correct [a,b,c] ) and then perturbing them with some noise. We want to then find the set [a,b,c]\n",
"that was used, despite the noise.\n",
"\n",
"In this problem, we can modify the following settings:\n",
"\n",
" 'method' : Choose which optimization scheme to employ\n",
" '# obs' : modifies the number of observations we have to constrain the parameters\n",
" 'stddev of obs' : modifies the uncertainty (1-sigma) of observations we have to constrain the parameters\n",
" 'setddev of params' : modifies the uncertainty (1-sigma) of a-priori values for the parameters\n",
" '# members' : modifies the number of ensemble members available for the ensemble Kalman filter method\n",
" 'Freeze random' : fixes the random numbers generated so that an experiment becomes repeatable\n"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"scrolled": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/auke/anaconda3/lib/python3.7/site-packages/IPython/html.py:14: ShimWarning: The `IPython.html` package has been deprecated since IPython 4.0. You should import from `notebook` instead. `IPython.html.widgets` has moved to `ipywidgets`.\n",
" \"`IPython.html.widgets` has moved to `ipywidgets`.\", ShimWarning)\n"
"#### Play with the number of observations, and their errors. What is the effect on the estimate?\n",
"\n",
"#### When few observations are present, their location on the curve matters. Verify this.\n",
"\n",
"#### The maximum likelihood method needs more observations than parameters. See what happens if you violate this.\n",
"\n",
"#### Switch to the Kalman Filter solution, and try again to go to only a few observations. What do you see?\n",
"\n",
"#### Try to manipulate the settings such that (a) the solution sticks to the prior, (b) the solution ignores the prior. Think of the Kalman Gain Matrix from the lecture when making this experiment.\n",
"\n",
"#### Now switch to the Ensemble Kalman Filter. Can you make it work really well, and really poorly?"
]
}
],
"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.6"
},
"varInspector": {
"cols": {
"lenName": 16,
"lenType": 16,
"lenVar": 40
},
"kernels_config": {
"python": {
"delete_cmd_postfix": "",
"delete_cmd_prefix": "del ",
"library": "var_list.py",
"varRefreshCmd": "print(var_dic_list())"
},
"r": {
"delete_cmd_postfix": ") ",
"delete_cmd_prefix": "rm(",
"library": "var_list.r",
"varRefreshCmd": "cat(var_dic_list()) "
}
},
"types_to_exclude": [
"module",
"function",
"builtin_function_or_method",
"instance",
"_Feature"
],
"window_display": false
}
},
"nbformat": 4,
"nbformat_minor": 1
}
%% Cell type:markdown id: tags:
This program demonstrates the different ways to optimize a system with 3 unknowns given a set of observations
and prior information on the PDF of the unknowns. There are currently three options:
(1) a maximum likelihood estimate: this estimate is a simple linear algebra solution to an overconstrained problem. It does not use
prior information on parameters or their uncertainty as it is assumed that a unique solution exists. It does return an estimate of
the uncertainty in the optimal parameters. This method is an alternative to the type of minimization routines that come standard with
many software packages such as IDL, MATLAB, python, etc.
(2) a full Kalman Filter estimate: this method is also suitable when there are more unknowns that observations, or when the solution is
not uniquely defined. In that case, a-priori information on the probability density (PDF) of the solution is used. The method returns an
optimal parameter set and its full covariance.
(3) an ensemble Kalman Filter estimate: the Kalman filter equations are again solved, but now using the statistics of an ensemble to
find the optimum. This solution is in pri nciple inferior to the one above (statistical errors) but sometimes easier to calculate when the structure of the underlying model is complex. The method returns an ensemble of solutions to the problem, from which a mean and covariance
can be derived that is the same as in the full Kalman filter.
The problem that we are optimizing here is simple. We have an equation of the form:
y = a + b*x + c* x^2
and the question is whether we can find a set [a,b,c] given imperfect observations (y). These imperfect observations are created by taking a
'perfect' y (outcome of the equation given correct [a,b,c] ) and then perturbing them with some noise. We want to then find the set [a,b,c]
that was used, despite the noise.
In this problem, we can modify the following settings:
'method' : Choose which optimization scheme to employ
'# obs' : modifies the number of observations we have to constrain the parameters
'stddev of obs' : modifies the uncertainty (1-sigma) of observations we have to constrain the parameters
'setddev of params' : modifies the uncertainty (1-sigma) of a-priori values for the parameters
'# members' : modifies the number of ensemble members available for the ensemble Kalman filter method
'Freeze random' : fixes the random numbers generated so that an experiment becomes repeatable
%% Cell type:code id: tags:
``` python
%matplotlibinline
importmpld3
importmatplotlib.pyplotasplt
fromlib_toymodelimport*
mpld3.enable_notebook()
plt.style.use('bmh')
%run-ilib_toymodel
```
%% Output
/Users/auke/anaconda3/lib/python3.7/site-packages/IPython/html.py:14: ShimWarning: The `IPython.html` package has been deprecated since IPython 4.0. You should import from `notebook` instead. `IPython.html.widgets` has moved to `ipywidgets`.
"`IPython.html.widgets` has moved to `ipywidgets`.", ShimWarning)
%% Cell type:code id: tags:
``` python
test.cov
```
%% Output
array([[ 3.08323522, 0.13203644, -1.06827182],
[ 0.13203644, 0.81091986, -0.05830712],
[-1.06827182, -0.05830712, 0.65578418]])
%% Cell type:markdown id: tags:
# Suggested exercises
#### Play with the number of observations, and their errors. What is the effect on the estimate?
#### When few observations are present, their location on the curve matters. Verify this.
#### The maximum likelihood method needs more observations than parameters. See what happens if you violate this.
#### Switch to the Kalman Filter solution, and try again to go to only a few observations. What do you see?
#### Try to manipulate the settings such that (a) the solution sticks to the prior, (b) the solution ignores the prior. Think of the Kalman Gain Matrix from the lecture when making this experiment.
#### Now switch to the Ensemble Kalman Filter. Can you make it work really well, and really poorly?
