GaussianCopula.cpp

Go to the documentation of this file.

#include "GaussianCopula.h"
 
GaussianCopula::GaussianCopula(unsigned int seed, RealMatrix R) : GenericCopula(seed) {
    this->R = nearest_psd(R);
    this->L = cholesky_L(this->R);
    this->logdet = this->R.slogdet().logdet;
    this->R_inv = this->R.inv();
}
 
std::vector<std::vector<double>> GaussianCopula::sample_u(std::size_t n) {
    std::vector<std::vector<double>> U;
 
    for (std::size_t i = 0; i < n; i++) {
        U.emplace_back(sample_u());
    }
 
    return U;
}
 
std::vector<double> GaussianCopula::sample_u() {
    std::size_t d = L.cols();
 
    RealMatrix z (d, 1);
    for (std::size_t j = 0; j < d; j++) {
        z.at(j, 0) = gsl_ran_ugaussian(eng_.get());   // z follows MN(0, 1)
    }
 
    z = L * z; // z follows MN(0, R)
    std::vector<double> u (d, 0.0);
    for (std::size_t j = 0; j < d; j++) {
        u[j] = std::clamp(gsl_cdf_ugaussian_P(z.at(j, 0)), CLIP_U, 1 - CLIP_U);   // u follows C(R)
    }
 
    return u;
}
 
double GaussianCopula::log_density(Vector u) {
    std::size_t d = u.size();
    RealMatrix z (d, 1);
 
    for (size_t i = 0; i < d; i++) {
        z.at(i, 0) = gsl_cdf_ugaussian_Pinv(std::clamp(u[i], 1e-13, 1. - 1e-13));
    }
 
    // return -0.5 * logdet - 0.5 * (z.transpose() * (R_inv - eye(d)) * z).at(0, 0);
    return -0.5 * (z.transpose() * (R_inv - eye(d)) * z).at(0, 0);
}
 
RealMatrix GaussianCopula::dlog_density(std::vector<double> u) {
    std::size_t d = u.size();
    RealMatrix z (d, 1);
 
    for (size_t i = 0; i < d; i++) {
        z.at(i, 0) = gsl_cdf_ugaussian_Pinv(std::clamp(u[i], 1e-13, 1. - 1e-13));
    }
 
    RealMatrix Az = (R_inv - eye(d)) * z;
    RealMatrix dlogc (-Az);
    for (size_t i = 0; i < d; i++) {
        dlogc.at(i, 0) /= gsl_ran_ugaussian_pdf(z.at(i, 0));
    }
    
    return dlogc;
}
 
RealMatrix GaussianCopula::ddlog_density(std::vector<double> u) {
    std::size_t d = u.size();
    RealMatrix z (d, 1);
 
    for (size_t i = 0; i < d; i++) {
        z.at(i, 0) = gsl_cdf_ugaussian_Pinv(std::clamp(u[i], 1e-13, 1. - 1e-13));
    }
 
    RealMatrix A = R_inv - eye(d);
    RealMatrix Az = A * z;
    RealMatrix ddlogc (d, d);
    for (size_t i = 0; i < d; i++) {
        for (size_t j = 0; j < d; j++) {
            double num = -A.at(i, j);
            if (i == j) num -= Az.at(i, 0) * z.at(i, 0); 
            double den = gsl_ran_ugaussian_pdf(z.at(i, 0)) * gsl_ran_ugaussian_pdf(z.at(j, 0));
            ddlogc.at(i, j) = num / den;
        }
    }
    
    return ddlogc;
}
 
LogDensityDiff GaussianCopula::log_c_dc_ddc(std::vector<double> u) {
    std::size_t d = u.size();
    RealMatrix z (d, 1);
    Vector phi_z (d, 0.0);
 
    for (size_t i = 0; i < d; i++) {
        z.at(i, 0) = gsl_cdf_ugaussian_Pinv(std::clamp(u[i], 1e-13, 1. - 1e-13));
        phi_z[i] = gsl_ran_ugaussian_pdf(z.at(i, 0));
    }
 
    RealMatrix A = R_inv - eye(d);
    RealMatrix Az = A * z;
    RealMatrix dlogc (d, 1);
    RealMatrix ddlogc (d, d);
 
    double log_c = -0.5 * (z.transpose() * Az).at(0, 0);
    for (size_t i = 0; i < d; i++) {
        dlogc.at(i, 0) = -Az.at(i, 0) / phi_z[i];
 
        for (size_t j = 0; j < d; j++) {
            double num = -A.at(i, j);
            if (i == j) num -= Az.at(i, 0) * z.at(i, 0); 
            double den = phi_z[i] * phi_z[j];
            ddlogc.at(i, j) = num / den;
        }
    }
    
    return {log_c, dlogc, ddlogc};
}