.. 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_gpd_rain.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_gpd_rain.py: Estimate a GPD on the daily rainfall data ========================================= .. GENERATED FROM PYTHON SOURCE LINES 6-20 In this example, we illustrate various techniques of extreme value modeling applied to the daily rainfall accumulations in south-west England, over the period 1914-1962. Readers should refer to [coles2001]_ to get more details. We illustrate techniques to: - estimate a stationary and a non stationary GPD, - estimate a return level, using: - the log-likelihood function, - the profile log-likelihood function. .. GENERATED FROM PYTHON SOURCE LINES 20-25 .. code-block:: Python import openturns as ot import openturns.experimental as otexp import openturns.viewer as otv from openturns.usecases import coles .. GENERATED FROM PYTHON SOURCE LINES 26-27 First, we load the Rain dataset. We start by looking at it through time. .. GENERATED FROM PYTHON SOURCE LINES 27-40 .. code-block:: Python dataRain = coles.Coles().rain print(dataRain[:10]) graph = ot.Graph( "Daily rainfall accumulations SW England", "day", "level (mm)", True, "" ) days = ot.Sample([[i] for i in range(len(dataRain))]) cloud = ot.Cloud(days, dataRain) cloud.setColor("red") cloud.setPointStyle(",") graph.add(cloud) graph.setIntegerXTick(True) view = otv.View(graph) .. image-sg:: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_estimate_gpd_rain_001.png :alt: Daily rainfall accumulations SW England :srcset: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_estimate_gpd_rain_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none [ x ] 0 : [ 0 ] 1 : [ 2.3 ] 2 : [ 1.3 ] 3 : [ 6.9 ] 4 : [ 4.6 ] 5 : [ 0 ] 6 : [ 1 ] 7 : [ 1.5 ] 8 : [ 1.8 ] 9 : [ 1.8 ] .. GENERATED FROM PYTHON SOURCE LINES 41-51 In order to select a threshold upon which the GPD model will be fitted, we draw the mean residual life plot with approximate :math:`95\%` confidence interval. It appears that the graph is linear from the threshold around :math:`u=30`. Then, it decays sharply although with a linear trend. We should be tempted to choose :math:`u=60` but there are only 6 exceedances of the threshold :math:`u=60`. So it is not enough to make meaningful inferences. Moreover, the graph is not reliable for large values of :math:`u` due to the limited amount of data on which the estimate and the confidence interval are based. For all these reasons, it appears preferable to chose :math:`u=30`. .. GENERATED FROM PYTHON SOURCE LINES 51-55 .. code-block:: Python factory = ot.GeneralizedParetoFactory() graph = factory.drawMeanResidualLife(dataRain) view = otv.View(graph) .. image-sg:: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_estimate_gpd_rain_002.png :alt: Mean residual life plot :srcset: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_estimate_gpd_rain_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 56-61 To support that choice, we draw the parameter stability plots. We see that the perturbations appear small relatively to sampling errors. We can see that the change in pattern observed in the mean residual life plot is still apparent here for high thresholds. Hence, we choose the threshold :math:`u=30`. .. GENERATED FROM PYTHON SOURCE LINES 61-65 .. code-block:: Python u_range = ot.Interval(0.5, 50.0) graph = factory.drawParameterThresholdStability(dataRain, u_range) view = otv.View(graph, figure_kw={"figsize": (6.0, 6.0)}) .. image-sg:: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_estimate_gpd_rain_003.png :alt: , Modified scale threshold stability, Shape threshold stability :srcset: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_estimate_gpd_rain_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 66-83 **Stationary GPD 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 consider the model :math:`\mathcal{M}_0` defined by: .. math:: :nowrap: \begin{align*} \sigma(t) & = \sigma\\ \xi(t) & = \xi \end{align*} We estimate the parameters of the GPD distribution by maximizing the log-likelihood of the data for the selecte threshold :math:`u=30`. .. GENERATED FROM PYTHON SOURCE LINES 83-86 .. code-block:: Python u = 30 result_LL = factory.buildMethodOfLikelihoodMaximizationEstimator(dataRain, u) .. GENERATED FROM PYTHON SOURCE LINES 87-88 We get the fitted GPD and its parameters of :math:`(\hat{\sigma}, \hat{\xi}, u)`. .. GENERATED FROM PYTHON SOURCE LINES 88-94 .. code-block:: Python fitted_GPD = result_LL.getDistribution() desc = fitted_GPD.getParameterDescription() param = fitted_GPD.getParameter() print(", ".join([f"{p}: {value:.3f}" for p, value in zip(desc, param)])) print("log-likelihood = ", result_LL.getLogLikelihood()) .. rst-class:: sphx-glr-script-out .. code-block:: none sigma: 7.446, xi: 0.184, u: 30.000 log-likelihood = -485.0937375341342 .. GENERATED FROM PYTHON SOURCE LINES 95-97 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 97-101 .. 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 : BlockIndependentDistribution(Normal(mu = [7.44573,0.184112], sigma = [0.959381,0.101172], R = [[ 1 -0.675368 ] [ -0.675368 1 ]]), Dirac(point = [30])) .. GENERATED FROM PYTHON SOURCE LINES 102-103 We get the covariance matrix and the standard deviation of :math:`(\hat{\sigma}, \hat{\xi}, \hat{\xi})`. .. GENERATED FROM PYTHON SOURCE LINES 103-106 .. 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.920412 -0.0655531 0 ] [ -0.0655531 0.0102358 0 ] [ 0 0 0 ]] Standard dev = [0.959381,0.101172,0] .. GENERATED FROM PYTHON SOURCE LINES 107-108 We get the marginal confidence intervals of order 0.95. .. GENERATED FROM PYTHON SOURCE LINES 108-113 .. code-block:: Python order = 0.95 for i in range(2): # exclude u parameter (fixed) ci = parameterEstimate.getMarginal(i).computeBilateralConfidenceInterval(order) print(desc[i] + ":", ci) .. rst-class:: sphx-glr-script-out .. code-block:: none sigma: [5.56538, 9.32608] xi: [-0.0141821, 0.382406] .. GENERATED FROM PYTHON SOURCE LINES 114-121 At last, we can validate the inference result thanks to 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 121-125 .. code-block:: Python validation = otexp.GeneralizedParetoValidation(result_LL, dataRain) graph = validation.drawDiagnosticPlot() view = otv.View(graph) .. image-sg:: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_estimate_gpd_rain_004.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_gpd_rain_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 126-131 **Stationary GPD modeling via the profile log-likelihood function** Now, we use the profile log-likehood function with respect to the :math:`\xi` parameter rather than log-likehood function to estimate the parameters of the GPD. .. GENERATED FROM PYTHON SOURCE LINES 131-133 .. code-block:: Python result_PLL = factory.buildMethodOfXiProfileLikelihoodEstimator(dataRain, u) .. GENERATED FROM PYTHON SOURCE LINES 134-137 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 137-141 .. 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_gpd_rain_005.png :alt: profile likelihood :srcset: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_estimate_gpd_rain_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 142-147 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 147-153 .. 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.0135616, 0.276116] .. GENERATED FROM PYTHON SOURCE LINES 154-162 **Return level estimate from the estimated stationary GPD** We evaluate the :math:`T`-year return level which corresponds to the :math:`m`-observation return level, where :math:`m = T*n_y` with :math:`n_y` the number of observations per year. Here, we have daily observations, hence :math:`n_y = 365`. As we assumed that the observations were independent, the extremal index is :math:`\theta=1` which is the default value. The method also provides the asymptotic distribution of the estimator :math:`\hat{z}_m`. .. GENERATED FROM PYTHON SOURCE LINES 162-178 .. code-block:: Python ny = 365 T10 = 10 zm_10 = factory.buildReturnLevelEstimator(result_LL, dataRain, T10 * ny) 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}") T100 = 100 zm_100 = factory.buildReturnLevelEstimator(result_LL, dataRain, T100 * ny) 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 Maximum log-likelihood function : 10-year return level = [65.9517] CI = [55.6696, 76.2339] 100-year return level = [106.284] CI = [65.4932, 147.075] .. GENERATED FROM PYTHON SOURCE LINES 179-183 **Return level estimate via the profile log-likelihood function of a stationary GPD** We can estimate the :math:`m`-observation return level :math:`z_m` directly from the data using the profile likelihood with respect to :math:`z_m`. .. GENERATED FROM PYTHON SOURCE LINES 183-189 .. code-block:: Python result_zm_100_PLL = factory.buildReturnLevelProfileLikelihoodEstimator( dataRain, u, T100 * ny ) zm_100_PLL = result_zm_100_PLL.getParameter() print(f"100-year return level (profile) = {zm_100_PLL}") .. rst-class:: sphx-glr-script-out .. code-block:: none 100-year return level (profile) = 106.32802943042479 .. GENERATED FROM PYTHON SOURCE LINES 190-192 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. .. GENERATED FROM PYTHON SOURCE LINES 192-197 .. code-block:: Python result_zm_100_PLL.setConfidenceLevel(0.95) return_level_ci100 = result_zm_100_PLL.getParameterConfidenceInterval() print("Maximum profile log-likelihood function : ") print(f"CI={return_level_ci100}") .. rst-class:: sphx-glr-script-out .. code-block:: none Maximum profile log-likelihood function : CI=[80.8575, 184.988] .. GENERATED FROM PYTHON SOURCE LINES 198-200 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 200-202 .. code-block:: Python view = otv.View(result_zm_100_PLL.drawProfileLikelihoodFunction()) .. image-sg:: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_estimate_gpd_rain_006.png :alt: profile likelihood :srcset: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_estimate_gpd_rain_006.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 203-238 **Non stationary GPD modeling via the log-likelihood function** Now, we want to check whether it is necessary to model the time dependency over the observation period. We have to define the functional basis for each parameter of the GPD model. Even if we have the possibility to affect a time-varying model to each of the 2 parameters :math:`(\sigma, \xi)`, it is strongly recommended not to vary the shape parameter :math:`\xi`. We suppose that :math:`\sigma` is linear with time, and that the other parameters remain 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. This method is the default one; - the *None* method where :math:`c = 0` and :math:`d = 1`: in that case, data are not normalized. We consider the model :math:`\mathcal{M}_1` defined by: .. math:: :nowrap: \begin{align*} \sigma(t) & = \beta_1 + \beta_2\tau(t) \\ \xi(t) & = \beta_3 \end{align*} .. GENERATED FROM PYTHON SOURCE LINES 238-244 .. code-block:: Python constant = ot.SymbolicFunction(["t"], ["1.0"]) basis = ot.Basis([ot.SymbolicFunction(["t"], ["t"]), constant]) # basis for mu, sigma, xi sigmaIndices = [0, 1] # linear xiIndices = [1] # stationary .. GENERATED FROM PYTHON SOURCE LINES 245-246 We need to get the time stamps (in days here). .. GENERATED FROM PYTHON SOURCE LINES 246-248 .. code-block:: Python timeStamps = ot.Sample([[i + 1] for i in range(len(dataRain))]) .. GENERATED FROM PYTHON SOURCE LINES 249-251 We can now estimate the list of coefficients :math:`\vect{\beta} = (\beta_1, \beta_2, \beta_3)` using the log-likelihood of the data. .. GENERATED FROM PYTHON SOURCE LINES 251-260 .. code-block:: Python result_NonStatLL = factory.buildTimeVarying( dataRain, u, timeStamps, basis, sigmaIndices, xiIndices ) beta = result_NonStatLL.getOptimalParameter() print(f"beta = {beta}") print(f"sigma(t) = {beta[1]:.4f} * tau(t) + {beta[0]:.4f}") print(f"xi = {beta[2]:.4f}") print(f"Max log-likelihood = {result_NonStatLL.getLogLikelihood()}") .. rst-class:: sphx-glr-script-out .. code-block:: none beta = [0.882226,7.01759,0.18162] sigma(t) = 7.0176 * tau(t) + 0.8822 xi = 0.1816 Max log-likelihood = -484.8061273999882 .. GENERATED FROM PYTHON SOURCE LINES 261-263 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 263-282 .. 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 = [0.8822255895863612; 0.8822255895863612] Conf interval for beta_2 = [7.017592488770449; 7.017592488770449] Conf interval for beta_3 = [0.18161982259931925; 0.18161982259931925] .. GENERATED FROM PYTHON SOURCE LINES 283-284 You can get the expression of the normalizing function :math:`t \mapsto \tau(t)`: .. GENERATED FROM PYTHON SOURCE LINES 284-289 .. 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=[1] constant=[0] linear=[[ 5.70451e-05 ]] c = 1.0 1/d = 5.7045065601825445e-05 .. GENERATED FROM PYTHON SOURCE LINES 290-292 You can get the function :math:`t \mapsto \vect{\theta}(t)` where :math:`\vect{\theta}(t) = (\sigma(t), \xi(t))`. .. GENERATED FROM PYTHON SOURCE LINES 292-294 .. code-block:: Python functionTheta = result_NonStatLL.getParameterFunction() .. GENERATED FROM PYTHON SOURCE LINES 295-298 In order to compare different modelings, we get the optimal log-likelihood of the data for both stationary and non stationary models. The difference seems to be non significant enough, which means that the non stationary model does not really improve the quality of the modeling. .. GENERATED FROM PYTHON SOURCE LINES 298-302 .. code-block:: Python print("Max log-likelihood: ") print("Stationary model = ", result_LL.getLogLikelihood()) print("Non stationary linear sigma(t) model = ", result_NonStatLL.getLogLikelihood()) .. rst-class:: sphx-glr-script-out .. code-block:: none Max log-likelihood: Stationary model = -485.0937375341342 Non stationary linear sigma(t) model = -484.8061273999882 .. GENERATED FROM PYTHON SOURCE LINES 303-322 In order to draw some diagnostic plots similar to those drawn in the stationary case, we refer to the following result: if :math:`Y_t` is a non stationary GPD model parametrized by :math:`(\sigma(t), \xi(t), u)`, then the standardized variables :math:`\hat{Y}_t` defined by: .. math:: \hat{Y}_t = \dfrac{1}{\xi(t)} \log \left[1+ \xi(t)\left( \dfrac{Y_t-u}{\sigma(t)} \right)\right] have the Exponential distribution which is the GPD model with :math:`(\sigma, \xi, u) = (1, 0, 0)`. As a result, 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. using the transformed data compared to the Exponential model. We can see that the adequation seems similar to the graph of the stationary model. .. GENERATED FROM PYTHON SOURCE LINES 322-325 .. code-block:: Python graph = result_NonStatLL.drawDiagnosticPlot() view = otv.View(graph) .. image-sg:: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_estimate_gpd_rain_007.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_gpd_rain_007.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 326-335 We can draw the mean function :math:`t \mapsto \Expect{\mbox{GPD}(t)}`, defined for :math:`\xi <1` only: .. math:: \Expect{\mbox{GPD}(t)} = u + \dfrac{\sigma(t)}{1 - \xi(t)} 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 GPD distribution at time :math:`t`. Here, :math:`\sigma(t)` is a linear function and the other parameters are constant, so the mean and the quantile functions are also linear functions. .. GENERATED FROM PYTHON SOURCE LINES 335-363 .. code-block:: Python graph = ot.Graph( r"Maximum rain - Linear $\sigma(t)$", "day", "level (mm)", True, "", ) graph.setIntegerXTick(True) # data cloud = ot.Cloud(timeStamps, dataRain) cloud.setColor("red") graph.add(cloud) # mean function meandata = [ result_NonStatLL.getDistribution(t).getMean()[0] for t in timeStamps.asPoint() ] curve_meanPoints = ot.Curve(timeStamps.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_gpd_rain_008.png :alt: Maximum rain - Linear $\sigma(t)$ :srcset: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_estimate_gpd_rain_008.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 364-374 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)` and the model :math:`\mathcal{M}_1` is parametrized by :math:`(\beta_1, \beta_2, \beta_3)`: 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 there is no evidence of a linear trend for :math:`\mu`. .. GENERATED FROM PYTHON SOURCE LINES 374-386 .. code-block:: Python llh_LL = result_LL.getLogLikelihood() llh_NonStatLL = result_NonStatLL.getLogLikelihood() modelM0_Nb_param = 2 modelM1_Nb_param = 3 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 ? = True .. GENERATED FROM PYTHON SOURCE LINES 387-391 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 391-396 .. 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=0.58 alpha=0.05 p-value=0.45 .. GENERATED FROM PYTHON SOURCE LINES 397-406 We can test a linear trend in the log-scale parameter for :math:`\sigma(t)`: .. math:: :nowrap: \begin{align*} \sigma(t) & = exp(\beta_1 + \beta_2\tau(t)) \\ \xi(t) & = \beta_3 \end{align*} .. GENERATED FROM PYTHON SOURCE LINES 406-423 .. code-block:: Python sigmaLink = ot.SymbolicFunction("x", "exp(x)") result_NonStatLL_Link = factory.buildTimeVarying( dataRain, u, timeStamps, basis, sigmaIndices, xiIndices, sigmaLink ) beta = result_NonStatLL_Link.getOptimalParameter() print(f"beta = {beta}") print(f"sigma(t) = exp({beta[1]:.4f} * tau(t) + {beta[0]:.4f})") print(f"xi = {beta[2]:.4f}") print(f"Max log-likelihood = {result_NonStatLL.getLogLikelihood()}") # % # The maximized log-likelihood we obtain with the log-linear model is very similar # to the one we obtained with the linear model. Hence, there is no evidence of a time trend. # We draw the diagnostic plots which are similar to the previous ones. graph = result_NonStatLL_Link.drawDiagnosticPlot() view = otv.View(graph) .. image-sg:: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_estimate_gpd_rain_009.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_gpd_rain_009.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none beta = [0.343272,1.80443,0.197766] sigma(t) = exp(1.8044 * tau(t) + 0.3433) xi = 0.1978 Max log-likelihood = -484.8061273999882 .. GENERATED FROM PYTHON SOURCE LINES 424-425 .. code-block:: Python otv.View.ShowAll() .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 28.593 seconds) .. _sphx_glr_download_auto_data_analysis_distribution_fitting_plot_estimate_gpd_rain.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_gpd_rain.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_estimate_gpd_rain.py `