thouska / spotpy

A Statistical Parameter Optimization Tool
https://spotpy.readthedocs.io/en/latest/
MIT License
254 stars 152 forks source link

Objective function required for MCMC, DE-MCz? #273

Open arthur-e opened 3 years ago

arthur-e commented 3 years ago

An objective function shouldn't be required for MCMC, DE-MCz, as these techniques are sampling directly from the posterior distribution. The choice of parameter sets is guided by the transition probabilities and the acceptance ratio, not some misfit between model predictions and observations. The spotpy documentation acknowledges as much in the section "Rosenbrock:"

"MCMC and DEMCz do not use the objective function and have their own likelihood function (logarithmic probability density)."

As such, in my equivalent of spotpy_setup (really should be capitalized because it is a class), I set the objectivefunction to return None (I wanted to leave the method in place to document the spotpy API in case I should decide to change algorithms later):

def objectivefunction(self, simulation, evaluation, *args, **kwargs):
    return None

And yet, when running DE-MCz, I get errors related to the objectivefunction():

Traceback (most recent call last):
  File "scripts/spot.py", line 495, in <module>
    fire.Fire(CalibrationCLI)
  File "/home/arthur.endsley/python-env/smap-l4c/lib/python3.7/site-packages/fire/core.py", line 141, in Fire
    component_trace = _Fire(component, args, parsed_flag_args, context, name)
  File "/home/arthur.endsley/python-env/smap-l4c/lib/python3.7/site-packages/fire/core.py", line 471, in _Fire
    target=component.__name__)
  File "/home/arthur.endsley/python-env/smap-l4c/lib/python3.7/site-packages/fire/core.py", line 681, in _CallAndUpdateTrace
    component = fn(*varargs, **kwargs)
  File "scripts/spot.py", line 316, in tune_gpp
    algorithm.sample(samples, nChains = 2 * len(PARAMETERS['GPP']))
  File "/home/arthur.endsley/src/spotpy/build/lib/spotpy/algorithms/demcz.py", line 179, in sample
    likelist = self.postprocessing(i, vector, simulations, chains=rep)
  File "/home/arthur.endsley/src/spotpy/build/lib/spotpy/algorithms/_algorithm.py", line 411, in postprocessing
    self.status(like,params,block_print=block_print)
  File "/home/arthur.endsley/src/spotpy/build/lib/spotpy/algorithms/_algorithm.py", line 92, in __call__
    self.compare(objectivefunction, params)
  File "/home/arthur.endsley/src/spotpy/build/lib/spotpy/algorithms/_algorithm.py", line 72, in maximizer
    if objval > self.objectivefunction_max:
TypeError: '>' not supported between instances of 'NoneType' and 'float'

Notably, when running DE-MCz with a proper objective function (one that returns something other than None), the maximum objective function value never changes; it retains its initial value set in the code. Shouldn't there be some sort of check to see which algorithm is used and, for MCMC and DE-MCz, bypass these comparisons with the objective function?

arthur-e commented 3 years ago

In general, the implementation of objectivefunction() is very confusing. The documentation reads:

"In this case we select the Root Mean Squared Error, giving the simulations and the evaluations. As we want to minimize our function, we select a negative objective function at this point."

And the code sample shown is:

def objectivefunction(self,simulation,evaluation):
    objectivefunction= -spotpy.objectivefunctions.rmse(evaluation,simulation)      
    return objectivefunction

However, if this objectivefunction() was truly minimized it would never identify an optimal parameter set because RMSE would diverge toward negative infinity. I believe this code sample was used because many of the algorithms maximize the objective function, not minimize. DE-MCz will maximize the objective function (although, as mentioned, it shouldn't really use an objective function at all); ROPE will maximize the objective function.

The chief issue here is not that users have to be wary about setting up the objectivefunction(); rather, it is that the "optimization_direction" is hard-coded in spotpy, not exposed as an argument to _algorithm sub-classes. As such, users have to look at the source code to figure out what is going on.

thouska commented 3 years ago

Hi @arthur-e thank you for your message. Bayesian algorithms (such as the implemented MCMC and DREAM) use likelihood functions instead of objective functions to determine, to find out if a paricular model run was better or worse. So in the "def objectivefunction" one need to return an appropriate likelihood function, if one wants to use Bayesian optimization. I recommend this tutorial to get used to the general set up in spotpy and this paper for the theoretical background.

Spotpy is designed to allow the user a maximum of flexibility, this comes with the cost that some background knowledge of the user is required. Not every combination of settings is checked or makes sense. However, for those not want to get too much into the theory of model optimization, plenty of examples are given.

For suggestions to improve the documentation, it would help me to get a direct pull request using the source files.

arthur-e commented 3 years ago

Hi @thouska, thanks for your message. I am familiar with MCMC, both likelihood-free approaches and approximate Bayesian computing. I think the central issue is that stochastic gradient-based optimization approaches tend to minimize the objective function while likelihood-based approaches are based on maximizing the likelihood function. Since the optimization direction is hard-coded, not transparent to the user (i.e., not in the function signature or in API documentation), and as the phrase "objective function" is used universally, this can be confusing (DREAM will run smoothly with a function intended to be minimized, it just won't produce good results). I'll try to find time to put a Pull Request together for suggested documentation updates.

thouska commented 3 years ago

Sounds good, I would very much apriciate your valuable time to work on this warmly welcome your contribution!