二元多元copula优选matlab代码

#t-Copula, Gaussian Copula, Clayton Copula ,Gumbel Copula, Frank Copula函数，五种常见的Copula模型
##BIC进行优选

import matlab.unittest.diagnostics.ConstraintDiagnostic

    % 定义要测试的 copula 类型
    copula_types = {'Gaussian', 't', 'Clayton', 'Gumbel', 'Frank'};
    bic_values = zeros(length(copula_types), 1);

    % 使用 BIC 来选择最优的 copula
    for i = 1:length(copula_types)
        % 拟合 copula 模型
        if strcmp(copula_types{i}, 't')
            [rho_hat, nu_hat] = copulafit('t', [u_SM, u_VPD]);
            copula_model = struct('Rho', rho_hat, 'Nu', nu_hat);
            logL = sum(log(copulapdf('t', [u_SM, u_VPD], rho_hat, nu_hat)));
            k = numel(rho_hat) + 1; % t copula 的参数数量，包括自由度
        else
            copula_model = copulafit(copula_types{i}, [u_SM, u_VPD]);
            logL = sum(log(copulapdf(copula_types{i}, [u_SM, u_VPD], copula_model)));
            k = numel(copula_model); % 其他 copula 的参数数量
        end

        % 计算 BIC
        bic_values(i) = -2 * logL + k * log(length(u_SM));
    end

    % 选择 BIC 最小的模型
    [~, best_idx] = min(bic_values);
    best_copula_type = copula_types{best_idx};

    if strcmp(best_copula_type, 't')
        [rho_hat, nu_hat] = copulafit('t', [u_SM, u_VPD]);
        best_copula_model = struct('Rho', rho_hat, 'Nu', nu_hat);
    else
        best_copula_model = copulafit(best_copula_type, [u_SM, u_VPD]);
    end

    disp(['最佳 Copula 类型: ', best_copula_type]);