Calibration#
R.A. Collenteur, University of Graz
After a model is constructed, the model parameters can be estimated using the ml.solve method. It can (and will) happen that the model fit after solving is not as good as expected. This may be the result of the settings that are used to solve the model or the way the model was constructed. In this notebook common pitfalls and various tips and tricks that may help to improve the calibration of Pastas models are shared.
In general, the following strategy is advised to solve problems with the parameter estimation:
Check the input time series and solve settings
Change the initial parameters,
Change the model structure,
Change the solve method.
import pandas as pd
import pastas as ps
ps.show_versions()
ps.set_log_level("ERROR")
Pastas : 2.0.0
Python : 3.14.6
Numpy : 2.4.6
Pandas : 3.0.3
Scipy : 1.18.0
Matplotlib : 3.11.0
Numba : 0.66.0
Loading the data#
In the following code-block some example data is loaded. It is good practice to visualize all time series before creating the time series model.
head = (
pd.read_csv("data/B32C0639001.csv", parse_dates=["date"], index_col="date")
.squeeze()
.loc["1985":]
)
# Make this millimeters per day
evap = (
pd.read_csv("data/evap_260.csv", index_col=0, parse_dates=[0])
.squeeze()
.loc["1985":"2003"]
)
rain = (
pd.read_csv("data/rain_260.csv", index_col=0, parse_dates=[0])
.squeeze()
.loc["1985":"2003"]
)
ps.plots.series(head, [evap, rain], table=True);
Make a model#
Given the data above we create a Pastas model with a non-linear recharge model (ps.FlexModel) and a constant to simulate the groundwater level. We’ll use this model to show how we may analyse different types of problems and how to solve them.
ml = ps.Model(head)
ps.ArNoiseModel(model=ml)
rch = ps.rch.FlexModel()
rm = ps.RechargeModel(
model=ml, prec=rain, evap=evap, recharge=rch, rfunc=ps.Gamma(), name="rch"
)
Calibrating a model#
In the above code-block a Pastas model was created, but not yet solved. To solve the model we call ml.solve(). This method has quite a few options (see also the docstring of the method) that influence the model calibration, for example:
tmin/tmax: select the time period used for calibrationnoise: use a noise model to model the residuals or notfit_constant: fit the constant as a parameter or notwarmup: length of the warmup periodsolver: the solver that is used to estimate parametersfreq_obs: frequency of the observations to use during calibration
We start without providing any arguments to the solve method.
# ml.solve? ## Run this to see other solve options
ml.solve()
ml.plots.results(figsize=(10, 6));
Fit report head Fit Statistics
==================================================
nfev 45 EVP 71.05
nobs 463 R2 0.71
noise True RMSE 0.11
tmin 1985-01-15 00:00:00 AICc -2550.83
tmax 2005-10-14 00:00:00 BIC -2513.98
freq D Obj nan
freq_obs None ___
warmup 3650 days 00:00:00 Interp. No
Parameters (9 optimized)
==================================================
optimal initial vary
rch_A 0.347578 0.612817 True
rch_n 0.625093 1.000000 True
rch_a 257.719134 10.000000 True
rch_srmax 62.498196 250.000000 True
rch_lp 0.250000 0.250000 False
rch_ks 23.990774 100.000000 True
rch_gamma 2.950897 2.000000 True
rch_kv 0.961464 1.000000 True
rch_simax 2.000000 2.000000 False
constant_d 0.918679 1.356415 True
noise_alpha 95.409021 15.000000 True
The fit report and the Figure above show that the model is not that great. The parameters have large standard errors, the goodness-of-fit metrics are not that high, and the simulated time series shows a very different behavior to the observed groundwater level.
Checking the explanatory time series and solve settings#
A common pitfall is that there is a problem with the explanatory time series (e.g., precipitation, pumping discharge). This should be the first thing to check when the model fit is not as good as expected.
Length of Time Series: The time series should in principle be available for the entire period of calibration,
Warmup Period: For some models it is necessary that the time series are also available before the calibration period, during the warmup period. This is for example the case with the non-linear recharge models (e.g., FlexModel, Berendrecht).
Units of Time Series (1): While Pastas is in principle unitless, the units of the time series can impact the model calibration. For example, a pumping discharge provided in m \(^3\) /day may lead to very small parameter values (‘Gamma_A’) that are harder to estimate. If you end up with very small parameters for the gain parameter, it may help to rescale the input time series.
Units of Time Series (2): The initial parameters and bounds for the non-linear recharge models are set for precipitation and evaporaton time series provided in mm/day. Using these models with time series in m/day will give bad results.
Normalization of Time Series: Sometimes it can help to normalize the expanatory time series. For example, when using a river level that is high above a certain datum (e.g. tens of meters), it may help to subtract the mean water level from the time series first.
In the example model, many of these things are happening. First, the precipitation time series are not available for the entire calibration period. Secondly, because a non-linear model is applied, we need to to have precipitation and evaporation data before the calibration period starts (typically about one year is enough). We should therefore shorten the calibration period by using to 1986-2003. Note that we use 3650 days for the warmup period (warmup=3650 is the default), the last 365 days of which now has real precipitation and evaporation data . For the other 9 years the mean flux is used. Finally, the non-linear model requires the evaporation and precipitation in mm/day (unless we want to manually set all parameter bounds).
ml = ps.Model(head)
ps.ArNoiseModel(model=ml)
rch = ps.rch.FlexModel()
rm = ps.RechargeModel(
model=ml,
prec=rain * 1e3,
evap=evap * 1e3,
recharge=rch,
rfunc=ps.Gamma(),
name="rch",
)
ml.solve(tmin="1986", tmax="2003", report=False)
axes = ml.plots.results(
tmin="1975", figsize=(10, 6)
) # Use tmin=1975 to show warmup period
axes[0].axvline(pd.Timestamp("1986"), c="k", linestyle="--"); # Start of calibration
Changing the explanatory time series and using the correct calibration period definitely improve the model fit in this example. Changing the explanatory time series a bit generally helps to resolve many issues with the calibration. If this does not work, we may try to help the solver a bit.
Improving initial parameters#
Although Pastas tries to set sensible initial parameters when constructing a model, it occurs that the initial parameters set by Pastas are not a great place to start the search for the optimal parameters. In this case, it may be tried to manually adapt the initial parameters using the ml.set_parameter as follows:
ml.set_parameter(
"rch_n", initial=15.0, pmax=100.0
) # Clearly wrong, just for educational purposes
ml.solve(tmin="1986", tmax="2003", report=True)
---------------------------------------------------------------------------
KeyboardInterrupt Traceback (most recent call last)
Cell In[6], line 4
1 ml.set_parameter(
2 "rch_n", initial=15.0, pmax=100.0
3 ) # Clearly wrong, just for educational purposes
----> 4 ml.solve(tmin="1986", tmax="2003", report=True)
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/model.py:1017, in Model.solve(self, tmin, tmax, freq, warmup, solver, report, initial, weights, fit_constant, freq_obs, initialize, reset_settings, noise, **kwargs)
1014 LeastSquares(model=self)
1016 # Solve model
-> 1017 solve_success, result = self.solver.solve(weights=weights, **kwargs)
1018 # Update the parameters with the results from the optimization
1019 for column in result.columns:
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/solver/least_squares.py:830, in LeastSquares.solve(self, weights, **kwargs)
815 bounds = Bounds(
816 lb=pmin,
817 ub=pmax,
818 keep_feasible=True,
819 )
821 objfunction = partial(
822 self.objfunction,
823 noise=noise,
(...) 827 callback=kwargs.pop("callback", None),
828 )
--> 830 self.result = least_squares(
831 fun=objfunction,
832 x0=initial[vary],
833 jac=self.jac,
834 bounds=bounds,
835 method=self.method,
836 ftol=self.ftol,
837 xtol=self.xtol,
838 gtol=self.gtol,
839 x_scale=self.x_scale,
840 loss=self.loss,
841 f_scale=self.f_scale,
842 max_nfev=self.max_nfev,
843 diff_step=self.diff_step,
844 tr_solver=self.tr_solver,
845 **kwargs,
846 )
848 self.pcov = DataFrame(
849 self.get_covariances(
850 self.result.jac,
(...) 856 columns=parameters.index,
857 )
859 # Prepare return values
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/scipy/_lib/_util.py:660, in _workers_wrapper.<locals>.inner(*args, **kwds)
658 with MapWrapper(_workers) as mf:
659 kwargs['workers'] = mf
--> 660 return func(*args, **kwargs)
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/scipy/optimize/_lsq/least_squares.py:1019, in least_squares(fun, x0, jac, bounds, method, ftol, xtol, gtol, x_scale, loss, f_scale, diff_step, tr_solver, tr_options, jac_sparsity, max_nfev, verbose, args, kwargs, callback, workers)
1015 result = call_minpack(vector_fun.fun, x0, vector_fun.jac, ftol, xtol, gtol,
1016 max_nfev, x_scale, jac_method=jac)
1018 elif method == 'trf':
-> 1019 result = trf(vector_fun.fun, vector_fun.jac, x0, f0, J0, lb, ub, ftol, xtol,
1020 gtol, max_nfev, x_scale, loss_function, tr_solver,
1021 tr_options.copy(), verbose, callback=callback_wrapped)
1023 elif method == 'dogbox':
1024 if tr_solver == 'lsmr' and 'regularize' in tr_options:
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/scipy/optimize/_lsq/trf.py:124, in trf(fun, jac, x0, f0, J0, lb, ub, ftol, xtol, gtol, max_nfev, x_scale, loss_function, tr_solver, tr_options, verbose, callback)
120 return trf_no_bounds(
121 fun, jac, x0, f0, J0, ftol, xtol, gtol, max_nfev, x_scale,
122 loss_function, tr_solver, tr_options, verbose, callback=callback)
123 else:
--> 124 return trf_bounds(
125 fun, jac, x0, f0, J0, lb, ub, ftol, xtol, gtol, max_nfev, x_scale,
126 loss_function, tr_solver, tr_options, verbose, callback=callback)
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/scipy/optimize/_lsq/trf.py:376, in trf_bounds(fun, jac, x0, f0, J0, lb, ub, ftol, xtol, gtol, max_nfev, x_scale, loss_function, tr_solver, tr_options, verbose, callback)
372 f_true = f.copy()
374 cost = cost_new
--> 376 J = jac(x)
377 njev += 1
379 if loss_function is not None:
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/scipy/optimize/_differentiable_functions.py:750, in VectorFunction.jac(self, x)
748 def jac(self, x):
749 self._update_x(x)
--> 750 self._update_jac()
751 if hasattr(self.J, "astype"):
752 # returns a copy so that downstream can't overwrite the
753 # internal attribute. But one can't copy a LinearOperator
754 return self.J.astype(self.J.dtype)
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/scipy/optimize/_differentiable_functions.py:719, in VectorFunction._update_jac(self)
716 else:
717 self._njev += 1
--> 719 self.J = self.jac_wrapped(xp_copy(self.x), f0=self.f)
720 self.J_updated = True
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/scipy/optimize/_differentiable_functions.py:454, in _VectorJacWrapper.__call__(self, x, f0, **kwds)
452 self.njev += 1
453 elif self.jac in FD_METHODS:
--> 454 J, dct = approx_derivative(
455 self.fun,
456 x,
457 f0=f0,
458 **self.finite_diff_options,
459 )
460 self.nfev += dct['nfev']
462 if self.sparse_jacobian:
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/scipy/optimize/_numdiff.py:611, in approx_derivative(fun, x0, method, rel_step, abs_step, f0, bounds, sparsity, as_linear_operator, args, kwargs, full_output, workers)
609 with MapWrapper(workers) as mf:
610 if sparsity is None:
--> 611 J, _nfev = _dense_difference(fun_wrapped, x0, f0, h,
612 use_one_sided, method,
613 mf)
614 else:
615 if not issparse(sparsity) and len(sparsity) == 2:
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/scipy/optimize/_numdiff.py:738, in _dense_difference(fun, x0, f0, h, use_one_sided, method, workers)
735 l = next(gen)
736 u = next(gen)
--> 738 f1 = next(f_evals)
739 f2 = next(f_evals)
740 if one_sided:
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/scipy/optimize/_numdiff.py:912, in _Fun_Wrapper.__call__(self, x)
909 if xp.isdtype(x.dtype, "real floating"):
910 x = xp.astype(x, self.x0.dtype)
--> 912 f = np.atleast_1d(self.fun(x, *self.args, **self.kwargs))
913 if f.ndim > 1:
914 raise RuntimeError("`fun` return value has "
915 "more than 1 dimension.")
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/scipy/optimize/_differentiable_functions.py:424, in _VectorFunWrapper.__call__(self, x)
422 def __call__(self, x):
423 self.nfev += 1
--> 424 return np.atleast_1d(self.fun(x))
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/scipy/optimize/_lsq/least_squares.py:263, in _WrapArgsKwargs.__call__(self, x)
262 def __call__(self, x):
--> 263 return self.f(x, *self.args, **self.kwargs)
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/solver/least_squares.py:772, in LeastSquares.objfunction(self, p, noise, weights, initial, vary, callback)
770 par = initial
771 par[vary] = p
--> 772 return misfit(
773 ml=self.model, p=par, noise=noise, weights=weights, callback=callback
774 )
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/solver/objective_function.py:43, in misfit(ml, p, noise, weights, callback, returnseparate)
41 # Get the residuals or the noise
42 if noise:
---> 43 rv = ml.noise(p) * ml._noise_weights(p)
44 else:
45 rv = ml.residuals(p)
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/model.py:706, in Model.noise(self, p, tmin, tmax, freq, warmup)
703 p = self.get_parameters()
705 # Calculate the residuals
--> 706 res = self.residuals(p, tmin, tmax, freq, warmup)
707 p = p[-self.noisemodel.nparam :]
709 # Calculate the noise
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/model.py:611, in Model.residuals(self, p, tmin, tmax, freq, warmup)
606 freq_obs = (
607 freq if self.settings["freq_obs"] is None else self.settings["freq_obs"]
608 )
610 # simulate model
--> 611 sim = self.simulate(
612 p=p, tmin=tmin, tmax=tmax, freq=freq, warmup=warmup, return_warmup=False
613 )
615 # Get the oseries calibration series
616 obs = self.observations(tmin=tmin, tmax=tmax, freq=freq_obs)
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/model.py:509, in Model.simulate(self, p, tmin, tmax, freq, warmup, return_warmup)
501 # Get the simulation index and the time step
502 # Check if the requested index matches the model settings
503 if (
504 tmin == self.settings["tmin"]
505 and tmax == self.settings["tmax"]
506 and freq == self.settings["freq"]
507 and warmup == self.settings["warmup"]
508 ):
--> 509 sim_index = self.sim_index
510 else:
511 # simulate with the requested settings, but do not update
512 # the model settings, since this is just for one time
513 sim_index = _get_sim_index(
514 tmin=tmin - warmup,
515 tmax=tmax,
516 freq=freq,
517 time_offset=self.time_offset,
518 )
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/model.py:1389, in Model.sim_index(self)
1371 @property
1372 def sim_index(self) -> DatetimeIndex:
1373 """Property that returns the simulation index, including the warmup.
1374
1375 Using the tmin, tmax, freq, and warmup from the model
(...) 1383 model is simulated.
1384 """
1385 return _get_sim_index(
1386 tmin=self.settings["tmin"] - self.settings["warmup"],
1387 tmax=self.settings["tmax"],
1388 freq=self.settings["freq"],
-> 1389 time_offset=self.time_offset,
1390 )
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pastas/model.py:1361, in Model.time_offset(self)
1359 mask = t >= base
1360 if np.any(mask):
-> 1361 time_offsets.add(_get_time_offset(t[mask][0], freq))
1362 if len(time_offsets) > 1:
1363 msg = "The time-offset with the frequency is not the same for all stresses."
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pandas/core/indexes/base.py:5401, in Index.__getitem__(self, key)
5398 if is_integer(key) or is_float(key):
5399 # GH#44051 exclude bool, which would return a 2d ndarray
5400 key = com.cast_scalar_indexer(key)
-> 5401 return getitem(key)
5403 if isinstance(key, slice):
5404 # This case is separated from the conditional above to avoid
5405 # pessimization com.is_bool_indexer and ndim checks.
5406 return self._getitem_slice(key)
File ~/checkouts/readthedocs.org/user_builds/pastas/envs/latest/lib/python3.14/site-packages/pandas/core/arrays/datetimelike.py:393, in DatetimeLikeArrayMixin.__getitem__(self, key)
386 """
387 This getitem defers to the underlying array, which by-definition can
388 only handle list-likes, slices, and integer scalars
389 """
390 # Use cast as we know we will get back a DatetimeLikeArray or DTScalar,
391 # but skip evaluating the Union at runtime for performance
392 # (see https://github.com/pandas-dev/pandas/pull/44624)
--> 393 result = cast(Union[Self, DTScalarOrNaT], super().__getitem__(key))
394 if lib.is_scalar(result):
395 return result
File ~/.asdf/installs/python/3.14.6/lib/python3.14/annotationlib.py:317, in ForwardRef.__hash__(self)
316 def __hash__(self):
--> 317 return hash((
318 self.__forward_arg__,
319 self.__forward_module__,
320 id(self.__globals__), # dictionaries are not hashable, so hash by identity
321 self.__forward_is_class__,
322 ( # cells are not hashable as well
323 tuple(sorted([(name, id(cell)) for name, cell in self.__cell__.items()]))
324 if isinstance(self.__cell__, dict) else id(self.__cell__),
325 ),
326 self.__owner__,
327 tuple(sorted(self.__extra_names__.items())) if self.__extra_names__ else None,
328 ))
KeyboardInterrupt:
Often we do not know what good initial parameters are, but we do get a bad fit, like with this initial value for rch_n above. While solving the model with a noise model is recommended, it does make the parameter estimation more difficult and more sensitive to the initial parameter values. One solution that often helps is to first solve the model without a noise model, and then solve the model with a noise model but without re-initializing the parameters.
By default the parameters are initialized upon each solve, such that each time we call solve we obtain the same result. By setting initial=False we prevent the re-initialisation and use the optimal parameters as initial parameters. This can be done as follows:
# First solve without noise model
ml.del_noisemodel()
ml.solve(report=False, tmin="1986", tmax="2003")
# Then solve with noise model, but do not initialize the parameters
ps.ArNoiseModel(model=ml)
ml.solve(initial=False, tmin="1986", tmax="2003", report=True)
axes = ml.plots.results(figsize=(10, 6))
After solving the model without a noise model (providing the solver an easier problem), we solve again with the parameter estimated from the solve without a noise model. This generally works well. We may also choose to fix parameters that are hard to estimate, perhaps because they are correlated to other parameters, to certain values.
Changing the model structure#
At this point, one might start to think that the bad fit has something to do with the model structure. This could off course be an explanatory time series that is missing, but let’s assume that is not the case. One thing that might help is to change the response function. This can either be from a complicated function to a simpler function (e.g., Gamma to Exponential) or the other way around (e.g., Gamma to FourParam). Another option could be to change other parts of the model structure, for example by applying a non-linear recharge model instead of a linear model.
## Example to be added
More advanced solve options#
If all of the above does not work, and we still think we have the right model structure and explanatory time series, we can for example:
Don’t fit the constant. By default the constant (
constant_d) is estimated as a parameter in Pastas. In specific cases it may help to turn this option off (ml.solve(fit_constant=False)).Switch the solver.
ps.solver.LeastSquares()is used by default, butps.solver.Lmfit()provides a lot of different methods for the parameter estimation, from simple least_squares to the use of MCMC.Remove observations from the groundwater level time series . The use of high frequency measurements is known to cause issues when trying to solve a model when using a noise model. See also the example notebook “Reducing Autocorrelation”.
Summary of Tips & Tricks#
In this notebook a variety of methods to improve the calibration result and model fit for Pastas models were shown. Although a specific type of model was used here to demonstrate these methods, the strategy can be applied to other types of time series and model structures as well.
A summary of all tips and tricks that may help to improve the model calibration given below:
Change units of input time series
Normalize the input time series
Change calibration period
Lengthen the warmup period
Solve first without, then with a noise model
Manually change initial parameters
Fix parameters
Change response functions
Fit constant or not
Try a different solve method
Remove observations