如何基于Archimedean Copula由U与theta推导并发变量V？

我有两组变量，其中一组是另一组的函数。先使用copulas包拟合Copula，再对因变量做极值分析（Extreme Value Analysis, EVA）。最终希望通过因变量对应值U和Archimedean Copula的theta参数，获取并发的自变量对应值V，但这一步遇到困难。试过ChatGPT的方案，效率较低，寻求高效解决思路或帮助。

def check_and_coupla_fit(df, ice_dist, sped_dist):
    """
    check the copula and return the 
    """

    ice = ice_dist.cdf(df.jones)
    sped = sped_dist.cdf(df.sped)


    tau_values = []
    rho_values = []
    theta_values = []
    for copula_type in ['clayton', 'frank', 'gumbel']:
        copula = getattr(copulas, copula_type.capitalize() + 'Copula')
        copula_instance = copula()
        copula_instance.fit([ice, sped])
        
        # Compute tau and rho
        u, v = copula_instance.sample(1000)
        tau, rho = compute_tau_and_rho(u, v)
        
        tau_values.append(tau)
        rho_values.append(rho)
        theta_values.append(copula_instance.compute_theta())
        
    # Determine the best copula based on tau and rho
    best_idx = np.argmax(np.array(tau_values) + np.array(rho_values))
    
    # Create the result array
    result = np.full((3, 1), np.nan)
    result[best_idx] = theta_values[best_idx]
    
    return result

def compute_tau_and_rho(u, v):
    """
    Compute Kendall's tau and Spearman's correlation coefficient.
    """
    tau, _ = copulas.stats.kendall_tau(u, v)
    rho, _ = spearmanr(u, v)
    return tau, rho

def get_stat(df):
    """
    gets distribution parameters for both ice accretion and wind speed.
    """
    
    # Extract the year from the 'valid' column
    df['year'] = df['valid'].dt.year

    # Group by year and get the maximum annual ice accretion values from the dataset
    max_ice_values = df.groupby(['year'])['jones'].max().reset_index()

    df_max_with_sped = pd.merge(max_ice_values, df,
                                on = ['jones', 'year'],
                                how = 'inner')[['jones', 'sped']]
    
    c_sped,loc_sped, scale_sped = weibull_min.fit(df_max_with_sped.sped)
    c_ice, loc_ice, scale_ice = genextreme.fit(df_max_with_sped.jones)
    
    return c_ice, loc_ice, scale_ice, c_sped,loc_sped, scale_sped
def inv_gumbel(u, theta):
    """
    Inverse of the Gumbel copula for given u and theta.
    """
    return np.exp(-(-np.log(u))**(1/theta))


def inv_frank(u, theta):
    """
    Inverse of the Frank copula for given u and theta.
    """
    term = np.log(1 + (np.exp(-theta*u) - 1) * (np.exp(-theta) - 1) / (np.exp(-theta*u) - np.exp(-theta) - 1))
    return -1/theta * np.log(1 - term)


def inv_clayton(u, theta):
    """
    Inverse of the Clayton copula for given u and theta.
    """
    return (u**(-theta) + 1)**(-1/theta)
def get_concurrent_wind(arr, U):
    """
    transforms U to V
    arr: vector of thetas
    U: is standardized ice value
    """
    V = 0
    if np.isnan(arr[0]):
        # Clayton copula
        V = inv_clayton(U, arr[0])
        
    elif np.isnan(arr[1]):
        # Frank copula
        V = inv_frank(U, arr[1])
        
    elif np.isnan(arr[2]):
        # Gumbel copula
        V = inv_gumbel(U, arr[2])
        
    else:
        raise ValueError('Size of arr should be 3, not more!');

    return V
stations = df.station.unique()
# unique stations

Results = np.array([])
counter_rows = 0
proceeed_statistical_analysis = True
for counter_rows, station in bar(enumerate(stations), total=203):
    
    station_data = df.loc[df.station == station].copy()
    # EVA for each station
    station_data = station_data[['valid','jones']]
    # only accounts for conditions in which EVA can be implemented! 
    try:
        station_data.set_index('valid',
                               inplace = True)
        
        station_data = station_data.squeeze()
        
        model = EVA(station_data)
        
        model.get_extremes(
            method="BM",
            extremes_type="high",
            block_size="365.2425D",
            errors="ignore",
        )

        # Get max values for different return periods (50 and 500 yrs).
        model.fit_model()
        summary = model.get_summary(
        return_period=[50, 500],
        alpha=0.95,
        n_samples=1_000)
        summary = summary.iloc[:, 0]
        
        jones_50 = summary[50]
        jones_500 = summary[500]    

        
                                       
    except:
        jones_50 = jones_500 = sped_50 = sped_500 = np.NaN
        arr = np.zeros([1, 3])
        proceeed_statistical_analysis = False

    if proceeed_statistical_analysis:
        # Obtaining concurrent wind speed for 50 and 500 years return periods.
        station_data = df.loc[df.station == station].copy()
        # EVA for each station
        station_data = station_data[['valid','jones']]
        
        c_ice, loc_ice, scale_ice, c_sped,loc_sped, scale_sped = get_stat(data)
    
        ice_dist = genextreme(c = c_ice,
                              loc = loc_ice,
                              scale = scale_ice)
    
        sped_dist = weibull_min(c = c_sped,
                                loc = loc_sped,
                                scale = scale_sped)
        
        # acquire theta and the corresponding copula function for the observations in the ith station
        arr = check_and_coupla_fit(df = station_data,
                                   ice_dist = ice_dist,
                                   sped_dist = sped_dist)
    
        # computes concurrent 50 years wind speed
        sped_50 = get_concurrent_wind(arr, sped_dist.cdf(jones_50))
    
        # computes concurrent 500 years wind speed
        sped_500 = get_concurrent_wind(arr, sped_dist.cdf(jones_500))

    Temp1 = df.loc[df.station == station].copy()
    lon, lat,elev = Temp1.iloc[0,2], Temp1.iloc[0,3], Temp1.iloc[0,4]

    Temp = np.zeros(shape = [1, 10])
    Temp[0, 0] = lon
    Temp[0, 1] = lat
    Temp[0, 2] = elev
    Temp[0, 3] = jones_50
    Temp[0, 4] = jones_500
    Temp[0, 5] = sped_50
    Temp[0, 6] = sped_500
    Temp[0, 7:] = arr
    
    if counter_rows == 0:
        Results = Temp.copy()
    elif counter_rows >0:
        
        Results = np.concatenate([Results, Temp],
                                axis=0)
    else:
        raise ValueError('Something happened!');

Results = pd.DataFrame(Results,
                      columns=['lon',
                               'lat',
                               'elev',
                               'jones_50',
                               'jones_500',
                               'sped_50',
                               'sped_500',
                               'Clayton',
                               'Frank',
                               'Gumbel'])

Results = Results.loc[~np.isnan(Results.jones_50)]

更新8-01
尝试多种方法后明白，仅靠theta和U值无法直接得到V。不过从拟合的Copula函数（如Clayton Copula）生成足够多的样本，大概率能得到接近目标值的结果；若定义为优化问题，经若干迭代也能得到合理结果，耗时取决于计算机内存。

内容的提问来源于stack exchange，提问作者matt