import numpy as np
import [Link] as plt
from [Link] import spearmanr
import pandas as pd
from scipy import stats
import [Link] as opt
from [Link].empirical_distribution import ECDF
def empirical_cdf_inverse(data, u):
"""
Compute the inverse empirical CDF (quantile function).
"""
sorted_data = [Link](data)
quantile_index = [Link](int(u * (len(sorted_data)-1)), 0, len(sorted_data) -
1)
return sorted_data[quantile_index]
def gumbel_copula_inverse(U, theta, p):
"""
Computes the inverse of the Gumbel copula.
Parameters:
U (float): Marginal probability of X (must be between 0 and 1).
theta (float): Copula parameter.
p (float): Probability value for inversion.
Returns:
float: Inverted value V.
"""
# Validate that U is within the valid range (0,1)
if not (0 < U < 1):
raise AssertionError("U must be in the range (0, 1), but got U =
{}".format(U))
# A small constant to prevent log(0) or very small values
epsilon = 1e-6
def root_func(V):
# Ensure U and V are not too small for the logarithm and exponentiation
log_U = -[Link](U + epsilon)
log_V = -[Link](V + epsilon)
# Compute the Gumbel copula inverse equation
return [Link](-((log_U ** theta + log_V ** theta) ** (1 / theta))) - p
# Improved initial guess
initial_guess = U if U < 0.5 else 1 - U # Start from a reasonable region
# Solving for V using fsolve with better settings
V_sol = [Link](root_func, x0=initial_guess, xtol=1e-8, maxfev=10000)[0]
# If fsolve does not converge, use an alternative method (like brentq)
if [Link](V_sol):
V_sol = [Link](root_func, 1e-4, 1 - 1e-4) # Bracket the root between 0
and 1
return [Link](V_sol, 0.001, 0.999)
def compute_xt_empirical(x_t_minus_1, p01, p11, omega_xt, theta, empirical_data):
"""
Computes X_t from X_{t-1} using the Gumbel copula and empirical marginals.
"""
# Step 1: Compute empirical CDF of X_{t-1}
ecdf = ECDF(empirical_data)
U = ecdf(x_t_minus_1)
# Step 2: Calculate argument for inverse copula
inside_copula = ((p01 + p11) * omega_xt - p01) / p11
# Step 3: Apply inverse Gumbel copula
U_t = gumbel_copula_inverse(U, theta, inside_copula)
# Step 4: Transform U_t back to X_t using empirical inverse CDF
X_t = empirical_cdf_inverse(empirical_data, U_t)
return X_t
def trans(M):
n = len(M)
for i in range(n - 1):
for j in range(i + 1, n):
M[i][j], M[j][i] = M[j][i], M[i][j]
# def inverse_ecdf(data, p):
# data_sorted = [Link](data)
# n = len(data_sorted)
# ecdf_values = [Link](1, n + 1) / n
# return [Link](p, ecdf_values, data_sorted)
# Gumbel copula CDF
def gumbel_copula_cdf(u, v, theta):
if theta == 1:
return u * v # independent case
inner_term = ((-[Link](u))**theta + (-[Link](v))**theta)**(1/theta)
return [Link](-inner_term)
# Partial derivative of Gumbel copula with respect to u
def gumbel_copula_partial_u(u, v, theta):
if theta == 1:
return v # independent case
inner_term = (-[Link](u))**theta + (-[Link](v))**theta
c_uv = gumbel_copula_cdf(u, v, theta)
return c_uv * ((-[Link](u))**(theta-1)) / (u * inner_term**(1 - 1/theta))
# Conditional CDF using Gumbel copula for a given x (if x > 0) and rainfall value y
def O_Y_given_X_x_and_X_positive(y, x):
u, v = H_x(x), H_y(y)
c_uv_partial = gumbel_copula_partial_u(u, v, theta)
h_x = np.count_nonzero(rainfall_data[:-1] == x) / len(rainfall_data[:-1]) if
len(rainfall_data[:-1]) > 0 else 0
f_y = np.count_nonzero(rainfall_data[1:] == y) / len(rainfall_data[1:]) if
len(rainfall_data[1:]) > 0 else 0
denom = P_10 * h_x + P_11 * f_y
return (P_10 * h_x + P_11 * f_y * c_uv_partial) / denom if denom > 0 else
H_y(y)
# Conditional CDF for X = 0 (when previous day is dry)
def O_Y_given_X_0(y):
denom = P_00 + P_01
return (P_00 + P_01 * H_y(y)) / denom if denom > 0 else H_y(y)
def convert_strings_to_floats(input_array):
output_array = []
for element in input_array:
converted_float = float(element)
output_array.append(converted_float)
return output_array
# Simulation function
def simulate_markov_chain(T, p00, p01, p11):
xt =[Link](T)# Store results
o1= [Link](T)
pt1= [Link](T)
# # st=[]
# st = [Link](T)
for t in range(1, T-1):
O = [Link]()
# :
if xt[t - 1] == 0 and O > (p00 / (p00 + p01)):
p1 = ((p00+ p01) * O - p00) / p01
xt[t] = empirical_cdf_inverse(rainfall_data, p1)
o1[t]=O
pt1[t]=p1
# st[t]=0
# [Link]("*****")
st ="*****"
# print(p1)
# '
elif xt[t - 1] > 0 and O > (p01 / (p01 + p11)):
xt[t] = compute_xt_empirical(xt[t-1] , p01, p11, O, theta,
rainfall_data[1:])
# print("**1**")
# p2 = ((p01 + p11) * O - p01) / p11
# U=([Link]())
# p3 = gumbel_copula_inverse(U, theta, p2)
# # print(p3)
# xt[t] = empirical_cdf_inverse(rainfall_data, p3) # x is todays rain
o1[t] = O
pt1[t] = 0
# st[t] = 1
st = "**1**"
# [Link]("**1**")
else:
xt[t] = 0
o1[t] = O
pt1[t] = 0
# st[t] = 2
st = "**2**"
# [Link]("**2**")
# print("**2**")
print(t,"{:.3f}".format(xt[t]),"{:.3f}".format(xt[t-
1]),"{:.3f}".format(pt1[t]),"{:.3f}".format(o1[t]),st )
return xt
df1 = pd.read_csv('D:\\pythonProject1\\venv\\Rain_imd\\IMD_DataSET\\[Link]')
rainfall_data = [Link][:31*5, 7].values # Example selection of data
print(rainfall_data)
n = len(rainfall_data) - 1
print(n)
dry_prev = rainfall_data[:-1] == 0
wet_prev = rainfall_data[:-1] > 0
dry_curr = rainfall_data[1:] == 0
wet_curr = rainfall_data[1:] > 0
P_00 = [Link](dry_prev & dry_curr) / n
P_01 = [Link](dry_prev & wet_curr) / n
P_10 = [Link](wet_prev & dry_curr) / n
P_11 = [Link](wet_prev & wet_curr) / n
print("\nTransition Probabilities:")
print(f"P_00 (Dry-Dry): {P_00:.3f}")
print(f"P_01 (Dry-Wet): {P_01:.3f}")
print(f"P_10 (Wet-Dry): {P_10:.3f}")
print(f"P_11 (Wet-Wet): {P_11:.3f}")
H_x= [Link](1, len(rainfall_data[1:]) + 1) / len(rainfall_data[1:])
H_y= [Link](1, len(rainfall_data[:-1]) + 1) / len(rainfall_data[:-1])
wet_indices = (rainfall_data[:-1] > 0) & (rainfall_data[1:] > 0)
wet_y_data = rainfall_data[:-1][wet_indices]
wet_x_data = rainfall_data[1:][wet_indices]
rho, _ = spearmanr(wet_x_data, wet_y_data)
theta = 1 / (1 - rho) # Gumbel copula parameter
T = len(rainfall_data)
xt_values = simulate_markov_chain(T, P_00, P_01, P_11)
simulated_rainfall=xt_values
n = len(simulated_rainfall) - 1
dry_prev = simulated_rainfall[:-1] == 0
wet_prev = simulated_rainfall[:-1] > 0
dry_curr = simulated_rainfall[1:] == 0
wet_curr = simulated_rainfall[1:] > 0
P_00 = [Link](dry_prev & dry_curr) / n
P_01 = [Link](dry_prev & wet_curr) / n
P_10 = [Link](wet_prev & dry_curr) / n
P_11 = [Link](wet_prev & wet_curr) / n
print("\nTransition Probabilities for simulated rainfall:")
print(f"P_00 (Dry-Dry): {P_00:.3f}")
print(f"P_01 (Dry-Wet): {P_01:.3f}")
print(f"P_10 (Wet-Dry): {P_10:.3f}")
print(f"P_11 (Wet-Wet): {P_11:.3f}")
# y=0 # y is previous day rain i.e. xt-1
# if y==0:
# O1=([Link]())
# print(O1)
# f1=P_00/(P_00+P_01)
# if O1>f1:
# p1 = ((P_00+P_01)*O1-P_00)/P_01
# x = inverse_ecdf(rainfall_data[1:], p1) # x is todays rain
# else:
# x =0
# y=x
# elif y>0:
# O2 = ([Link]())
# # O2=0.23
# # print(O2)
# f2 = P_01 / (P_01 + P_11)
# # print(f2)
# if O2 > f2:
# p2 = ((P_01 + P_11) * O2 - P_01) / P_11
# U=([Link]())
# p3 = gumbel_copula_inverse(U, theta, p2)
# # print(p3)
# x = inverse_ecdf(rainfall_data[1:], p3) # x is todays rain
# else:
# x = 0
# y = x
# # print(simulated_rainfall)
# [Link](figsize=(10, 6))
# [Link](rainfall_data, simulated_rainfall)
# print(rainfall_data,simulated_rainfall)
[Link](figsize=(10, 6))
[Link](rainfall_data, label="Observed Rainfall", marker="o", linestyle="-",
alpha=0.7)
[Link](simulated_rainfall, label="Simulated Rainfall", marker="x",
linestyle="--", alpha=0.7)
[Link]("Time (days)")
[Link]("Rainfall (mm)")
[Link]("Observed vs Simulated Rainfall")
[Link]()
[Link](True)
[Link]()