-
Notifications
You must be signed in to change notification settings - Fork 4
/
Copy pathTime-dependent_GEV_fit_with_NAO.py
136 lines (121 loc) · 5.7 KB
/
Time-dependent_GEV_fit_with_NAO.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
"""Fit a time- and NAO-dependent GEV distribution to monthly maxima (MM).
NAO stands for the index of the North Atlantic Oscillation.
See tools_climate.py for information how to obtain this data.
Written by Markus Reinert, June 2020–July 2021.
"""
import numpy as np
from scipy import optimize
from advanced_GEV_analysis import negative_log_likelihood, Modifiers
from advanced_GEV_analysis import get_month_selection
from advanced_GEV_analysis import compute_amplitude_and_phase, check_significance
from tools_surge import load_data, Timeseries
from tools_climate import load_NOAA_data, ClimateIndex
data = load_data("Brest", Timeseries.SKEW_SURGE)
data_cli = load_NOAA_data(ClimateIndex.NAO)
# Get the maximum value in every month and its date-time,
# as well as the NAO index of this month
h_MM = []
t_MM = []
cli_MM = []
cli_years = np.array(np.round(data_cli["t_years"], 2), dtype=int)
for year in range(data["year_start"], data["year_end"] + 1):
for month in range(1, 13):
sel = get_month_selection(year, month, data["t"])
if np.any(sel):
i_max = np.argmax(data["h"][sel])
h_MM.append(data["h"][sel][i_max])
t_MM.append(data["t"][sel][i_max])
# Check if there is climate data for this year
if year in cli_years:
i_year = np.where(cli_years == year)[0]
# There must be exactly 12 monthly values (with possibly some being NaN)
assert len(i_year) == 12, "not exactly 12 monthly values of climate index"
cli_MM.append(data_cli["index"][i_year[month - 1]])
else:
cli_MM.append(np.nan)
h_MM = np.array(h_MM)
t_MM = np.array(t_MM)
cli_MM = np.array(cli_MM)
# Replace all NaNs in the climate data by 0s, because adding zero does not change anything
cli_MM[np.isnan(cli_MM)] = 0
# Fit time-independent GEV to the extreme values for a first estimate of the fit parameters
print("")
result = optimize.minimize(negative_log_likelihood, [10, 15, -0.1], args=(h_MM,))
if not result.success:
print("Warning:", result.message)
params_const = result.x
errors_const = np.sqrt(np.diag(result.hess_inv))
# Compute the (positive) log-likelihood of the model
LL_const = -negative_log_likelihood(params_const, h_MM)
print("\nEstimated GEV parameters for time-independent model:")
for i, name in enumerate(["mu", "sigma", "xi"]):
print("{:5} = {:7.4f} ± {:.4f}".format(name, params_const[i], errors_const[i]))
print("Parameter values of mu and sigma in cm, xi dimensionless.")
# Fit GEV model with linear trends and annual cycles in the parameters
print("")
args = (
h_MM,
t_MM,
[Modifiers.LINEAR_TREND, Modifiers.SEASONAL_OSCILLATION],
[Modifiers.LINEAR_TREND, Modifiers.SEASONAL_OSCILLATION],
)
mu, sigma, xi = params_const
params_initial = [round(mu), 0.1, 1, 1, round(sigma), 0.1, 1, 1, xi]
result = optimize.minimize(negative_log_likelihood, params_initial, args, options={"gtol": 5e-4})
if not result.success:
print("Warning:", result.message)
params_trends = result.x
errors_trends = np.sqrt(np.diag(result.hess_inv))
LL_trends = -negative_log_likelihood(params_trends, *args)
print("\nEstimated GEV parameters for model with linear trends and annual cycles:")
for i, name in enumerate([
"mu", "mu_trend", "mu_cos", "mu_sin",
"sigma", "sigma_trend", "sigma_cos", "sigma_sin", "xi",
]):
print("{:11} = {:7.4f} ± {:.4f}".format(name, params_trends[i], errors_trends[i]))
print("Trends in cm per century, xi dimensionless, other parameter values in cm.")
print("\nCompare time-dependent GEV model with time-independent model.")
D = 2 * (LL_trends - LL_const)
k = len(params_trends) - len(params_const)
check_significance(D, k)
# Fit GEV model with linear trends, annual cycles, and linear dependencies on climate index
print("")
args = (
h_MM,
t_MM,
[Modifiers.LINEAR_TREND, Modifiers.SEASONAL_OSCILLATION, Modifiers.CLIMATE_INDEX],
[Modifiers.LINEAR_TREND, Modifiers.SEASONAL_OSCILLATION, Modifiers.CLIMATE_INDEX],
None, # value of shape parameter xi, None means “best fit”, 0.0 means Gumbel distribution
cli_MM,
)
params_initial = [round(mu), 0.1, 1, 1, 0.1, round(sigma), 0.1, 1, 1, 0.1, xi]
result = optimize.minimize(negative_log_likelihood, params_initial, args, options={"gtol": 5e-4})
if not result.success:
print("Warning:", result.message)
params_cli = result.x
errors_cli = np.sqrt(np.diag(result.hess_inv))
LL_cli = -negative_log_likelihood(params_cli, *args)
print("\nEstimated GEV parameters for model with linear trends, annual cycles, "
"and linear {}-dependencies:".format(data_cli["name"]))
for i, name in enumerate([
"mu", "mu_trend", "mu_cos", "mu_sin", "mu_" + data_cli["name"],
"sigma", "sigma_trend", "sigma_cos", "sigma_sin", "sigma_" + data_cli["name"], "xi",
]):
print("{:11} = {:7.4f} ± {:.4f}".format(name, params_cli[i], errors_cli[i]))
print("Trends in cm per century, {} coefficients in cm per unit, "
"xi dimensionless, other parameter values in cm.".format(data_cli["name"]))
print("\nCompare time-dependent GEV model with linear {0}-dependencies "
"to model without {0}-dependencies.".format(data_cli["name"]))
D = 2 * (LL_cli - LL_trends)
k = len(params_cli) - len(params_trends)
check_significance(D, k)
# Convert coefficients of cos and sin to amplitude and phase
print("")
for name, i_cos in zip(["mu", "sigma"], [2, 7]):
print("Annual cycle in {}:".format(name))
amp, phase, std_amp, std_phase = compute_amplitude_and_phase(
params_cli[i_cos], params_cli[i_cos + 1],
errors_cli[i_cos], errors_cli[i_cos + 1],
)
print(" Amplitude = ({:.1f} ± {:.1f}) cm".format(amp, std_amp))
print(" Phase = ({:.0f} ± {:.0f})°".format(np.rad2deg(phase), np.rad2deg(std_phase)))