.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_data_analysis/distribution_fitting/plot_estimate_gev_racetime.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_data_analysis_distribution_fitting_plot_estimate_gev_racetime.py: Estimate a GEV on race times data ================================= .. GENERATED FROM PYTHON SOURCE LINES 6-55 In this example, we illustrate various techniques of extreme value modeling applied to the fatest annual race times for the women's 1500 meter event over the period 1975-1992. Readers should refer to [coles2001]_ to get more details. We illustrate techniques to: - estimate a stationary and a non stationary GEV, - estimate a return level, using: - the log-likelihood function, - the profile log-likelihood function. This analyse is particular as we are interested in modeling the annual minimum race times and not the annual maximum race times. In order to transform the minimum modeling into a maximum modeling, we proceeds as follows. We denote by :math:`\tilde{M}_n = \min (X_1, \dots, X_n)` where all the :math:`X_i` are independent and identically distributed variables. We introduce :math:`Y_i = -X_i` for :math:`1 \leq i \leq n`, and :math:`M_n = \min (Y_1, \dots, Y_n)`. Then, we have: .. math:: \tilde{M}_n = - M_n. We can show that if the renormalized maximum :math:`\tilde{M}_n` tends to the GEV distribution parametrized by :math:`(\mu, \sigma, \xi)`, then the renormalized minimum :math:`\tilde{M}_n` tends to the *GEV for minima* distribution parametrized by :math:`(\tilde{\mu}, \tilde{\sigma}, \tilde{\xi})` where: .. math:: :nowrap: \begin{align*} \tilde{\mu} & = -\mu \\ \tilde{\sigma} & = \sigma\\ \tilde{\xi} & = \xi \end{align*} The cumulated distribution function of :math:`\tilde{M}_n`, denoted by :math:`\tilde{G}`, is defined by: .. math:: \tilde{G}(z) = 1-G(-z) = 1-\exp \left( -\left[ 1-\tilde{\xi} \left( \dfrac{z-\tilde{\mu}}{\tilde{\sigma}}\right)\right]^{-1/\tilde{\xi}}\right) for all :math:`z` such that :math:`1-\tilde{\xi}(z-\tilde{\mu})/\tilde{\sigma} > 0`. In that example, we model the :math:`M_n` variable which is the annual maximum of the opposite race times. First, we load the race times dataset. We start by looking at them through time. .. GENERATED FROM PYTHON SOURCE LINES 55-72 .. code-block:: Python import openturns as ot import openturns.viewer as otv import openturns.experimental as otexp from openturns.usecases import coles data = coles.Coles().racetime print(data[:5]) graph = ot.Graph( "Fastest annual women's 1500m race times", "year", "race time (s)", True, "" ) cloud = ot.Cloud(data[:, :2]) cloud.setColor("red") graph.add(cloud) graph.setIntegerXTick(True) view = otv.View(graph) .. image-sg:: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_estimate_gev_racetime_001.png :alt: Fastest annual women's 1500m race times :srcset: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_estimate_gev_racetime_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none [ Year Race time ] 0 : [ 1972 241.371 ] 1 : [ 1973 244.616 ] 2 : [ 1974 242.301 ] 3 : [ 1975 246 ] 4 : [ 1976 235.994 ] .. GENERATED FROM PYTHON SOURCE LINES 73-74 We select the race times column. We transform them into their opposite values. .. GENERATED FROM PYTHON SOURCE LINES 74-76 .. code-block:: Python sample = -1.0 * data[:, 1] .. GENERATED FROM PYTHON SOURCE LINES 77-82 **Stationary GEV modeling via the log-likelihood function** We first assume that the dependence through time is negligible, so we first model the data as independent observations over the observation period. We estimate the parameters of the GEV distribution by maximizing the log-likelihood of the data. .. GENERATED FROM PYTHON SOURCE LINES 82-85 .. code-block:: Python factory = ot.GeneralizedExtremeValueFactory() result_LL = factory.buildMethodOfLikelihoodMaximizationEstimator(sample) .. GENERATED FROM PYTHON SOURCE LINES 86-88 We get the fitted GEV for the variable :math:`M_n` and its parameters :math:`(\hat{\mu}, \hat{\sigma}, \hat{\xi})`. .. GENERATED FROM PYTHON SOURCE LINES 88-93 .. code-block:: Python fitted_GEV = result_LL.getDistribution() desc = fitted_GEV.getParameterDescription() param = fitted_GEV.getParameter() print(", ".join([f"{p}: {value:.3f}" for p, value in zip(desc, param)])) .. rst-class:: sphx-glr-script-out .. code-block:: none mu: -239.260, sigma: 3.643, xi: -0.470 .. GENERATED FROM PYTHON SOURCE LINES 94-96 We get the asymptotic distribution of the estimator :math:`(\hat{\mu}, \hat{\sigma}, \hat{\xi})`. In that case, the asymptotic distribution is normal. .. GENERATED FROM PYTHON SOURCE LINES 96-100 .. code-block:: Python parameterEstimate = result_LL.getParameterDistribution() print("Asymptotic distribution of the estimator : ") print(parameterEstimate) .. rst-class:: sphx-glr-script-out .. code-block:: none Asymptotic distribution of the estimator : Normal(mu = [-239.26,3.64268,-0.469834], sigma = [0.891274,0.778059,0.236181], R = [[ 1 -0.113124 -0.406653 ] [ -0.113124 1 -0.706993 ] [ -0.406653 -0.706993 1 ]]) .. GENERATED FROM PYTHON SOURCE LINES 101-102 We get the covariance matrix and the standard deviation of :math:`(\hat{\mu}, \hat{\sigma}, \hat{\xi})`. .. GENERATED FROM PYTHON SOURCE LINES 102-105 .. code-block:: Python print("Cov matrix = \n", parameterEstimate.getCovariance()) print("Standard dev = ", parameterEstimate.getStandardDeviation()) .. rst-class:: sphx-glr-script-out .. code-block:: none Cov matrix = [[ 0.79437 -0.0784475 -0.0856013 ] [ -0.0784475 0.605376 -0.129919 ] [ -0.0856013 -0.129919 0.0557815 ]] Standard dev = [0.891274,0.778059,0.236181] .. GENERATED FROM PYTHON SOURCE LINES 106-112 At last, we can validate the inference result thanks the 4 usual diagnostic plots: - the probability-probability pot, - the quantile-quantile pot, - the return level plot, - the data histogram and the density of the fitted model. .. GENERATED FROM PYTHON SOURCE LINES 112-117 .. code-block:: Python validation = otexp.GeneralizedExtremeValueValidation(result_LL, sample) graph = validation.drawDiagnosticPlot() view = otv.View(graph) .. image-sg:: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_estimate_gev_racetime_002.png :alt: , Sample versus model PP-plot, Sample versus model QQ-plot, Return level plot, Density :srcset: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_estimate_gev_racetime_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 118-121 **Stationary GEV modeling via the profile log-likelihood function** Now, we use the profile log-likehood function rather than log-likehood function to estimate the parameters of the GEV. .. GENERATED FROM PYTHON SOURCE LINES 121-123 .. code-block:: Python result_PLL = factory.buildMethodOfProfileLikelihoodMaximizationEstimator(sample) .. GENERATED FROM PYTHON SOURCE LINES 124-127 The following graph allows one to get the profile log-likelihood plot. It also indicates the optimal value of :math:`\xi`, the maximum profile log-likelihood and the confidence interval for :math:`\xi` of order 0.95 (which is the default value). .. GENERATED FROM PYTHON SOURCE LINES 127-131 .. code-block:: Python order = 0.95 result_PLL.setConfidenceLevel(order) view = otv.View(result_PLL.drawProfileLikelihoodFunction()) .. image-sg:: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_estimate_gev_racetime_003.png :alt: profile likelihood :srcset: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_estimate_gev_racetime_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 132-137 We can get the numerical values of the confidence interval: it appears to be a bit smaller with the interval obtained from the profile log-likelihood function than with the log-likelihood function. Note that if the order requested is too high, the confidence interval might not be calculated because one of its bound is out of the definition domain of the log-likelihood function. .. GENERATED FROM PYTHON SOURCE LINES 137-143 .. code-block:: Python try: print("Confidence interval for xi = ", result_PLL.getParameterConfidenceInterval()) except Exception as ex: print(type(ex)) pass .. rst-class:: sphx-glr-script-out .. code-block:: none Confidence interval for xi = [-0.791703, -0.156295] .. GENERATED FROM PYTHON SOURCE LINES 144-156 **Return level estimate from the estimated stationary GEV** We estimate the :math:`m`-block return level :math:`z_m`: it is computed as a particular quantile of the GEV model estimated using the log-likelihood function. We just have to use the maximum log-likelihood estimator built in the previous section. The return level of :math:`M_n` and :math:`\tilde{M}_n` have opposite values. As the data are annual sea-levels, each block corresponds to one year: the 10-year return level corresponds to :math:`m=10` and the 100-year return level corresponds to :math:`m=100`. The method provides the asymptotic distribution of the estimator :math:`\hat{z}_m` of :math:`M_n` which mean is the return-level estimate. .. GENERATED FROM PYTHON SOURCE LINES 156-163 .. code-block:: Python zm_10 = factory.buildReturnLevelEstimator(result_LL, 10.0) return_level_10 = zm_10.getMean() print("Maximum log-likelihood function : ") print(f"10-year return level={return_level_10}") return_level_ci10 = zm_10.computeBilateralConfidenceInterval(0.95) print(f"CI={return_level_ci10}") .. rst-class:: sphx-glr-script-out .. code-block:: none Maximum log-likelihood function : 10-year return level=[-234.2] CI=[-235.552, -232.849] .. GENERATED FROM PYTHON SOURCE LINES 164-171 .. code-block:: Python zm_100 = factory.buildReturnLevelEstimator(result_LL, 100.0) return_level_100 = zm_100.getMean() print(f"100-year return level={return_level_100}") return_level_ci100 = zm_100.computeBilateralConfidenceInterval(0.95) print(f"CI={return_level_ci100}") .. rst-class:: sphx-glr-script-out .. code-block:: none 100-year return level=[-232.4] CI=[-235.071, -229.729] .. GENERATED FROM PYTHON SOURCE LINES 172-176 **Return level estimate via the profile log-likelihood function of a stationary GEV** We can estimate the :math:`m`-block return level :math:`z_m` directly from the data using the profile likelihood with respect to :math:`z_m`. .. GENERATED FROM PYTHON SOURCE LINES 176-180 .. code-block:: Python result_zm_10_PLL = factory.buildReturnLevelProfileLikelihoodEstimator(sample, 10.0) zm_10_PLL = result_zm_10_PLL.getParameter() print(f"10 years return level (profile)={zm_10_PLL}") .. rst-class:: sphx-glr-script-out .. code-block:: none 10 years return level (profile)=-234.20091228177077 .. GENERATED FROM PYTHON SOURCE LINES 181-185 We can get the confidence interval of :math:`z_m`: once more, it appears to be a bit smaller than the interval obtained from the log-likelihood function. As for the confidence interval of :math:`\xi`, dependeding on the order requested, the interval might not be calculated. .. GENERATED FROM PYTHON SOURCE LINES 185-194 .. code-block:: Python result_zm_10_PLL.setConfidenceLevel(0.95) try: return_level_ci10 = result_zm_10_PLL.getParameterConfidenceInterval() except Exception as ex: print(type(ex)) pass print("Maximum profile log-likelihood function : ") print(f"CI={return_level_ci10}") .. rst-class:: sphx-glr-script-out .. code-block:: none Maximum profile log-likelihood function : CI=[-235.485, -232.301] .. GENERATED FROM PYTHON SOURCE LINES 195-197 We can also plot the profile log-likelihood function and get the confidence interval, the optimal value of :math:`z_m` and its confidence interval. .. GENERATED FROM PYTHON SOURCE LINES 197-199 .. code-block:: Python view = otv.View(result_zm_10_PLL.drawProfileLikelihoodFunction()) .. image-sg:: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_estimate_gev_racetime_004.png :alt: profile likelihood :srcset: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_estimate_gev_racetime_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 200-209 **Non stationary GEV modeling via the log-likelihood function** If we look at the data carefully, we see that the pattern of variation has not remained constant over the observation period. There is an increase in the data through time. We want to model this trend because a slight increase in extreme sea-levels might have a significant impact on the safety of coastal flood defenses. We still work on the :math:`M_n` variable. First we need to get the time stamps (in years here). .. GENERATED FROM PYTHON SOURCE LINES 209-211 .. code-block:: Python timeStamps = data[:, 0] .. GENERATED FROM PYTHON SOURCE LINES 212-239 Then, we define the functional basis for each parameter of the GEV model. Even if we have the possibility to affect a time-varying model to each of the 3 parameters :math:`(\mu, \sigma, \xi)`, it is strongly recommended not to vary the parameter :math:`\xi` and to let it constant. For numerical reasons, it is strongly recommended to normalize all the data as follows: .. math:: \tau(t) = \dfrac{t-c}{d} where: - the *CenterReduce* method where :math:`c = \dfrac{1}{n} \sum_{i=1}^n t_i` is the mean time stamps and :math:`d = \sqrt{\dfrac{1}{n} \sum_{i=1}^n (t_i-c)^2}` is the standard deviation of the time stamps; - the *MinMax* method where :math:`c = t_1` is the initial time and :math:`d = t_n-t_1` the final time; - the *None* method where :math:`c = 0` and :math:`d = 1`: in that case, data are not normalized. We suppose that :math:`\mu` is linear in time, and that the other parameters remain constant. .. math:: :nowrap: \begin{align*} \mu(t) & = \beta_1 + \beta_2\tau(t) \\ \sigma(t) & = \beta_3 \\ \xi(t) & = \beta_4 \end{align*} .. GENERATED FROM PYTHON SOURCE LINES 239-245 .. code-block:: Python constant = ot.SymbolicFunction(["t"], ["1.0"]) basis_lin = ot.Basis([constant, ot.SymbolicFunction(["t"], ["t"])]) basis_cst = ot.Basis([constant]) # basis for mu, sigma, xi basis_coll = [basis_lin, basis_cst, basis_cst] .. GENERATED FROM PYTHON SOURCE LINES 246-257 We can now estimate the list of coefficients :math:`\vect{\beta} = (\beta_1, \beta_2, \beta_3, \beta_4)` using the log-likelihood of the data. We test the 3 normalizing methods and both initial points in order to evaluate their impact on the results. We can see that: - both normalization methods lead to the same result for :math:`\beta_1`, :math:`\beta_3` and :math:`\beta_4` (note that :math:`\beta_2` depends on the normalization function), - both initial points lead to the same result when the data have been normalized, - it is very important to normalize all the data: if not, the result strongly depends on the initial point and it differs from the result obtained with normalized data. The results are not optimal in that case since the associated log-likelihood are much smaller than those obtained with normalized data. .. GENERATED FROM PYTHON SOURCE LINES 257-279 .. code-block:: Python initiPoint_list = list() initiPoint_list.append("Gumbel") initiPoint_list.append("Static") normMethod_list = list() normMethod_list.append("MinMax") normMethod_list.append("CenterReduce") normMethod_list.append("None") ot.ResourceMap.SetAsUnsignedInteger( "GeneralizedExtremeValueFactory-MaximumEvaluationNumber", 1000000 ) print("Linear mu(t) model: ") for normMeth in normMethod_list: for initPoint in initiPoint_list: print("normMeth, initPoint = ", normMeth, initPoint) # The ot.Function() is the identity function. result = factory.buildTimeVarying( sample, timeStamps, basis_coll, ot.Function(), initPoint, normMeth ) beta = result.getOptimalParameter() print("beta1, beta2, beta3, beta4 = ", beta) print("Max log-likelihood = ", result.getLogLikelihood()) .. rst-class:: sphx-glr-script-out .. code-block:: none Linear mu(t) model: normMeth, initPoint = MinMax Gumbel beta1, beta2, beta3, beta4 = [-242.635,6.25764,2.73173,-0.201355] Max log-likelihood = -51.8940418474017 normMeth, initPoint = MinMax Static beta1, beta2, beta3, beta4 = [-242.635,6.25777,2.73171,-0.201359] Max log-likelihood = -51.894041787764984 normMeth, initPoint = CenterReduce Gumbel beta1, beta2, beta3, beta4 = [-239.506,1.94197,2.73166,-0.201303] Max log-likelihood = -51.89404131401922 normMeth, initPoint = CenterReduce Static beta1, beta2, beta3, beta4 = [-239.506,1.94201,2.73163,-0.201294] Max log-likelihood = -51.894041292150014 normMeth, initPoint = None Gumbel beta1, beta2, beta3, beta4 = [-239.942,-0.000578171,2.80599,0.0958746] Max log-likelihood = -61.18805537945545 normMeth, initPoint = None Static beta1, beta2, beta3, beta4 = [-239.26,1.15331e-06,3.6427,-0.469838] Max log-likelihood = -54.622262531545104 .. GENERATED FROM PYTHON SOURCE LINES 280-282 According to the previous results, we choose the *MinMax* normalization method and the *Gumbel* initial point. This initial point is cheaper than the *Static* one as it requires no optimization computation. .. GENERATED FROM PYTHON SOURCE LINES 282-292 .. code-block:: Python result_NonStatLL = factory.buildTimeVarying( sample, timeStamps, basis_coll, ot.Function(), "Gumbel", "MinMax" ) beta = result_NonStatLL.getOptimalParameter() print("Linear mu(t) model : ") print("beta1, beta2, beta3, beta_4 = ", beta) print(f"mu(t) = {beta[0]:.4f} + {beta[1]:.4f} * tau") print(f"sigma = = {beta[2]:.4f}") print(f"xi = = {beta[3]:.4f}") .. rst-class:: sphx-glr-script-out .. code-block:: none Linear mu(t) model : beta1, beta2, beta3, beta_4 = [-242.635,6.25764,2.73173,-0.201355] mu(t) = -242.6347 + 6.2576 * tau sigma = = 2.7317 xi = = -0.2014 .. GENERATED FROM PYTHON SOURCE LINES 293-294 You can get the expression of the normalizing function :math:`t \mapsto \tau(t)`: .. GENERATED FROM PYTHON SOURCE LINES 294-299 .. code-block:: Python normFunc = result_NonStatLL.getNormalizationFunction() print("Function tau(t): ", normFunc) print("c = ", normFunc.getEvaluation().getImplementation().getCenter()[0]) print("1/d = ", normFunc.getEvaluation().getImplementation().getLinear()[0, 0]) .. rst-class:: sphx-glr-script-out .. code-block:: none Function tau(t): class=LinearFunction name=Unnamed implementation=class=LinearEvaluation name=Unnamed center=[1972] constant=[0] linear=[[ 0.05 ]] c = 1972.0 1/d = 0.05 .. GENERATED FROM PYTHON SOURCE LINES 300-302 You can get the function :math:`t \mapsto \vect{\theta}(t)` where :math:`\vect{\theta}(t) = (\mu(t), \sigma(t), \xi(t))`. .. GENERATED FROM PYTHON SOURCE LINES 302-304 .. code-block:: Python functionTheta = result_NonStatLL.getParameterFunction() .. GENERATED FROM PYTHON SOURCE LINES 305-307 We get the asymptotic distribution of :math:`\vect{\beta}` to compute some confidence intervals of the estimates, for example of order :math:`p = 0.95`. .. GENERATED FROM PYTHON SOURCE LINES 307-327 .. code-block:: Python dist_beta = result_NonStatLL.getParameterDistribution() condifence_level = 0.95 for i in range(beta.getSize()): lower_bound = dist_beta.getMarginal(i).computeQuantile((1 - condifence_level) / 2)[ 0 ] upper_bound = dist_beta.getMarginal(i).computeQuantile((1 + condifence_level) / 2)[ 0 ] print( "Conf interval for beta_" + str(i + 1) + " = [" + str(lower_bound) + "; " + str(upper_bound) + "]" ) .. rst-class:: sphx-glr-script-out .. code-block:: none Conf interval for beta_1 = [-245.03004506706299; -240.23928370487715] Conf interval for beta_2 = [1.4859158023150973; 11.029371006754019] Conf interval for beta_3 = [1.4710369765042874; 3.99242034343222] Conf interval for beta_4 = [-0.7111936015189477; 0.3084839056248185] .. GENERATED FROM PYTHON SOURCE LINES 328-330 In order to compare different modelings, we get the optimal log-likelihood of the data for both stationary and non stationary models. The difference is significant enough to be in favor of the non stationary model. .. GENERATED FROM PYTHON SOURCE LINES 330-334 .. code-block:: Python print("Max log-likelihood: ") print("Stationary model = ", result_LL.getLogLikelihood()) print("Non stationary linear mu(t) model = ", result_NonStatLL.getLogLikelihood()) .. rst-class:: sphx-glr-script-out .. code-block:: none Max log-likelihood: Stationary model = -54.622265554594776 Non stationary linear mu(t) model = -51.8940418474017 .. GENERATED FROM PYTHON SOURCE LINES 335-355 We can draw the mean function :math:`t \mapsto \Expect{\mbox{GEV}(t)}`. Be careful, it is not the function :math:`t \mapsto \mu(t)`. As a matter of fact, the mean is defined for :math:`\xi <1` only and in that case, for :math:`\xi \neq 0`, we have: .. math:: \Expect{\mbox{GEV}(t)} = \mu(t) + \dfrac{\sigma(t)}{\xi(t)} (\Gamma(1-\xi(t))-1) and for :math:`\xi = 0`, we have: .. math:: \Expect{\mbox{GEV}(t)} = \mu(t) + \sigma(t)\gamma where :math:`\gamma` is the Euler constant. We can also draw the function :math:`t \mapsto q_p(t)` where :math:`q_p(t)` is the quantile of order :math:`p` of the GEV distribution at time :math:`t`. Here, :math:`\mu(t)` is a linear function and the other parameters are constant, so the mean and the quantile functions are also linear functions. The graph confirms the increase of the annual maximum sea-levels through time. .. GENERATED FROM PYTHON SOURCE LINES 355-380 .. code-block:: Python graph = ot.Graph( r"Fatest annual race times - Linear $\mu(t)$", "year", "race time (m)", True, "" ) dataModified = data * ot.Point([1.0, -1.0]) graph.setIntegerXTick(True) # data cloud = ot.Cloud(dataModified) cloud.setColor("red") graph.add(cloud) # mean function meandata = [ result_NonStatLL.getDistribution(t).getMean()[0] for t in data[:, 0].asPoint() ] curve_meanPoints = ot.Curve(data[:, 0].asPoint(), meandata) graph.add(curve_meanPoints) # quantile function graphQuantile = result_NonStatLL.drawQuantileFunction(0.95) drawQuant = graphQuantile.getDrawable(0) drawQuant = graphQuantile.getDrawable(0) drawQuant.setLineStyle("dashed") graph.add(drawQuant) graph.setLegends(["data", "mean function", "quantile 0.95 function"]) graph.setLegendPosition("lower right") view = otv.View(graph) .. image-sg:: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_estimate_gev_racetime_005.png :alt: Fatest annual race times - Linear $\mu(t)$ :srcset: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_estimate_gev_racetime_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 381-392 At last, we can test the validity of the stationary model :math:`\mathcal{M}_0` relative to the model with time varying parameters :math:`\mathcal{M}_1`. The model :math:`\mathcal{M}_0` is parametrized by :math:`(\beta_1, \beta_3, \beta_4)` and the model :math:`\mathcal{M}_1` is parametrized by :math:`(\beta_1, \beta_2, \beta_3, \beta_4)`: so we have :math:`\mathcal{M}_0 \subset \mathcal{M}_1`. We use the Likelihood Ratio test. The null hypothesis is the stationary model :math:`\mathcal{M}_0`. The Type I error :math:`\alpha` is taken equal to 0.05. This test confirms that the dependence through time is not negligible: it means that the linear model:math:`\mu(t)` component explains a large variation in the data. .. GENERATED FROM PYTHON SOURCE LINES 392-405 .. code-block:: Python llh_LL = result_LL.getLogLikelihood() llh_NonStatLL = result_NonStatLL.getLogLikelihood() modelM0_Nb_param = 3 modelM1_Nb_param = 4 resultLikRatioTest = ot.HypothesisTest.LikelihoodRatioTest( modelM0_Nb_param, llh_LL, modelM1_Nb_param, llh_NonStatLL, 0.05 ) accepted = resultLikRatioTest.getBinaryQualityMeasure() print( f"Hypothesis H0 (stationary model) vs H1 (linear mu(t) model): accepted ? = {accepted}" ) .. rst-class:: sphx-glr-script-out .. code-block:: none Hypothesis H0 (stationary model) vs H1 (linear mu(t) model): accepted ? = False .. GENERATED FROM PYTHON SOURCE LINES 406-410 We detail the statistics of the Likelihood Ratio test: the deviance statistics :math:`\mathcal{D}_p` follows a :math:`\chi^2_1` distribution. The model :math:`\mathcal{M}_0` is rejected if the deviance statistics estimated on the data is greater than the threshold :math:`c_{\alpha}` or if the p-value is less than the Type I error :math:`\alpha = 0.05`. .. GENERATED FROM PYTHON SOURCE LINES 410-414 .. code-block:: Python print(f"Dp={resultLikRatioTest.getStatistic():.2f}") print(f"alpha={resultLikRatioTest.getThreshold():.2f}") print(f"p-value={resultLikRatioTest.getPValue():.2f}") .. rst-class:: sphx-glr-script-out .. code-block:: none Dp=5.46 alpha=0.05 p-value=0.02 .. GENERATED FROM PYTHON SOURCE LINES 415-425 We can perform the same study with a quadratic model for :math:`\mu(t)`: .. math:: :nowrap: \begin{align*} \mu(t) & = \beta_1 + \beta_2 \tau(t) + \beta_3\tau(t)^2 \\ \sigma(t) & = \beta_4 \\ \xi(t) & = \beta_5 \end{align*} .. GENERATED FROM PYTHON SOURCE LINES 425-439 .. code-block:: Python basis_quad = ot.Basis( [constant, ot.SymbolicFunction(["t"], ["t"]), ot.SymbolicFunction(["t"], ["t^2"])] ) basis_coll_2 = [basis_quad, basis_cst, basis_cst] result_NonStatLL_2 = factory.buildTimeVarying( sample, timeStamps, basis_coll_2, ot.Function(), "Gumbel", "MinMax" ) beta = result_NonStatLL_2.getOptimalParameter() print("Quadratic mu(t) model : ") print("beta1, beta2, beta3, beta4, beta5 = ", beta) print(f"mu(t) = {beta[0]:.4f} + {beta[1]:.4f} * tau + {beta[2]:.4f} * tau^2") print(f"sigma = {beta[3]:.4f}") print(f"xi = {beta[4]:.4f}") .. rst-class:: sphx-glr-script-out .. code-block:: none Quadratic mu(t) model : beta1, beta2, beta3, beta4, beta5 = [-245.692,26.0513,-19.5972,2.28398,-0.182332] mu(t) = -245.6916 + 26.0513 * tau + -19.5972 * tau^2 sigma = 2.2840 xi = -0.1823 .. GENERATED FROM PYTHON SOURCE LINES 440-442 We get the asymptotic distribution of :math:`\vect{\beta}` to compute some confidence intervals of the estimates, for example of order :math:`p = 0.95`. .. GENERATED FROM PYTHON SOURCE LINES 442-461 .. code-block:: Python dist_beta = result_NonStatLL_2.getParameterDistribution() condifence_level = 0.95 for i in range(beta.getSize()): lower_bound = dist_beta.getMarginal(i).computeQuantile((1 - condifence_level) / 2)[ 0 ] upper_bound = dist_beta.getMarginal(i).computeQuantile((1 + condifence_level) / 2)[ 0 ] print( "Conf interval for beta_" + str(i + 1) + " = [" + str(lower_bound) + "; " + str(upper_bound) + "]" ) .. rst-class:: sphx-glr-script-out .. code-block:: none Conf interval for beta_1 = [-248.62771673783448; -242.75550382826742] Conf interval for beta_2 = [9.891070987824822; 42.21154696807771] Conf interval for beta_3 = [-35.56488620024279; -3.629473567407395] Conf interval for beta_4 = [1.4077954878683696; 3.1601684521206312] Conf interval for beta_5 = [-0.9157225901662904; 0.5510589623495372] .. GENERATED FROM PYTHON SOURCE LINES 462-464 In order to compare different modelings, we get the optimal log-likelihood of the data for both stationary and non stationary models. The difference is significant enough to be in favor of the non stationary model. .. GENERATED FROM PYTHON SOURCE LINES 464-467 .. code-block:: Python print("Max log-likelihood = ") print("Non stationary quadratic mu(t) model = ", result_NonStatLL_2.getLogLikelihood()) .. rst-class:: sphx-glr-script-out .. code-block:: none Max log-likelihood = Non stationary quadratic mu(t) model = -48.4684335974957 .. GENERATED FROM PYTHON SOURCE LINES 468-480 At last, we can test the validity of the non stationary model :math:`\mathcal{M}_0` where :math:`\mu(t)` is linear relative to the model :math:`\mathcal{M}_1` where :math:`\mu(t)` is quadratic. The model :math:`\mathcal{M}_0` is parametrized by :math:`(\beta_1, \beta_2, \beta_3, \beta_4)` and the model :math:`\mathcal{M}_1` is parametrized by :math:`(\beta_1, \beta_2, \beta_3, \beta_4, \beta_5)`: so we have :math:`\mathcal{M}_0 \subset \mathcal{M}_1`. We use the Likelihood Ratio test. The null hypothesis is the stationary model :math:`\mathcal{M}_0`. The Type I error :math:`\alpha` is taken equal to 0.05. This test confirms that the dependence through time is not negligible: it means that the :math:`\mu(t)` quadratic model explains even better a large variation in the data. .. GENERATED FROM PYTHON SOURCE LINES 480-489 .. code-block:: Python llh_NonStatLL_2 = result_NonStatLL_2.getLogLikelihood() resultLikRatioTest = ot.HypothesisTest.LikelihoodRatioTest( 4, llh_NonStatLL, 5, llh_NonStatLL_2, 0.05 ) accepted = resultLikRatioTest.getBinaryQualityMeasure() print( f"Hypothesis H0 (linear mu(t) model) vs H1 (quadratic mu(t) model): accepted ? = {accepted}" ) .. rst-class:: sphx-glr-script-out .. code-block:: none Hypothesis H0 (linear mu(t) model) vs H1 (quadratic mu(t) model): accepted ? = False .. GENERATED FROM PYTHON SOURCE LINES 490-491 .. code-block:: Python otv.View.ShowAll() .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 9.892 seconds) .. _sphx_glr_download_auto_data_analysis_distribution_fitting_plot_estimate_gev_racetime.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_estimate_gev_racetime.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_estimate_gev_racetime.py `