MCMC Estimation of Copula CAN'T WORK

未确定任何仓库文件或失败的测试。首先缩减围绕 joint_log_prob、construct_hmc 和 run_chain 的粘贴示例，并使用显示的初始值确定确切错误。完成要求包括：可复现的失败、记录在案的预期结果，以及经过验证且可正常工作的 MCMC copula 估计路径。

The samples is the initial values.

import numpy as np
import os
import re
os.environ["CUDA_VISIBLE_DEVICES"] ="-1"
from scipy.special import gamma
import scipy.stats

import tensorflow as tf
from tensorflow.keras.models import Model
from tensorflow.keras.layers import Dense
from tensorflow.keras.optimizers import Adam
from math import pi

import tensorflow_probability as tfp
tfd = tfp.distributions
tfb = tfp.bijectors
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import seaborn as sns
from scipy.integrate import quad

print("TFP Version", tfp.version)
print("TF Version", tf.version)

np.random.seed(1704)

import scipy
import seaborn
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import nquad

data = np.array(value_lists)[:, :, 0]
NUM_POSTERIOR_SAMPLES = 10000
NUM_BURNIN_ITERATIONS = int(0.25 * NUM_BURNIN_ITERATIONS)
NUM_ADAPTATION = int(0.5 * NUM_BURNIN_ITERATIONS)

class GaussianCopulaTriL(tfd.TransformedDistribution):
    """Takes a location, and lower triangular matrix for the Cholesky factor."""

    def __init__(self, loc, scale_tril):
        super(GaussianCopulaTriL, self).__init__(
            distribution=tfd.MultivariateNormalTriL(
                loc=loc,
                scale_tril=scale_tril),
            bijector=tfb.NormalCDF(),
            validate_args=False,
            name="GaussianCopulaTriLUniform")

    def _parameter_properties(self, dtype, num_classes=None):
        return dict(
            # Annotations may optionally specify properties, such as `event_ndims`,
            # `default_constraining_bijector_fn`, `specifies_shape`, etc.; see
            # the `ParameterProperties` documentation for details.
            loc=tfp.util.ParameterProperties(),
            scale_tril=tfp.util.ParameterProperties())


class WarpedGaussianCopula(tfd.TransformedDistribution):
    """Application of a Gaussian Copula on a list of target marginals.

    This implements an application of a Gaussian Copula. Given [x_0, ... x_n]
    which are distributed marginally (with CDF) [F_0, ... F_n],
    `GaussianCopula` represents an application of the Copula, such that the
    resulting multivariate distribution has the above specified marginals.

    The marginals are specified by `marginal_bijectors`: These are
    bijectors whose `inverse` encodes the CDF and `forward` the inverse CDF.

    block_sizes is a 1-D Tensor to determine splits for `marginal_bijectors`
    length should be same as length of `marginal_bijectors`.
    See tfb.Blockwise for details
    """

    def __init__(self, loc, scale_tril, marginal_bijectors, block_sizes=None):
        super(WarpedGaussianCopula, self).__init__(
            distribution=GaussianCopulaTriL(loc=loc, scale_tril=scale_tril),
            bijector=tfb.Blockwise(bijectors=marginal_bijectors,
                                   block_sizes=block_sizes),
            validate_args=False,
            name="GaussianCopula")

    def _parameter_properties(self, dtype, num_classes=None):
        return dict(
            # Annotations may optionally specify properties, such as `event_ndims`,
            # `default_constraining_bijector_fn`, `specifies_shape`, etc.; see
            # the `ParameterProperties` documentation for details.
            loc=tfp.util.ParameterProperties(),
            scale_tril=tfp.util.ParameterProperties(),
            marginal_bijectors = tfp.util.ParameterProperties(),
            block_sizes = tfp.util.ParameterProperties())

def joint_log_prob(data,
                       correlation,
                       ijv_loc,
                       ijv_scale,
                       cca_loc,
                       cca_scale):
    """
    Estimate the mean and variance of a Normal distribution.

    Parameters
    ----------
      @data: the observed data, binary (0, 1)
      @mu:   the mean parameter. This is "learnable"
      @sigma: the variance parameter. This is "learnable"
    """

    correlation_prior = tfd.Normal(loc=0., scale=1.)
    ijv_loc_prior     = tfd.Normal(loc=400., scale=200.)
    ijv_scale_prior   = tfd.Normal(loc=100., scale=50.)
    cca_loc_prior     = tfd.Normal(loc=50., scale=50.)
    cca_scale_prior   = tfd.Normal(loc=20., scale=10.)
    zero_prior        = tfd.Normal(loc=0., scale=1.)


    if True:

        rv_data = WarpedGaussianCopula(
                loc=[0., 0.],
                scale_tril=[[1., 0.], [correlation, tf.sqrt(1. - correlation ** 2)]],
                #These encode the marginals we want. In this case we want X_0 has
                #Kumaraswamy marginal, and X_1 has Gumbel marginal.

            marginal_bijectors=[
                tfb.Invert(tfb.KumaraswamyCDF(ijv_loc, ijv_scale)),
                tfb.Invert(tfb.GumbelCDF(loc=cca_loc, scale=cca_scale))])


    return (
            correlation_prior.log_prob(correlation) +
            ijv_loc_prior.log_prob(ijv_loc) +
            ijv_scale_prior.log_prob(ijv_scale) +
            cca_loc_prior.log_prob(cca_loc) +
            cca_scale_prior.log_prob(cca_scale) +
            tf.reduce_sum(rv_data.log_prob(data))
    )


def prediction(correlation_,
               ijv_loc_,
               ijv_scal_,
               cca_loc_,
               cca_scale_):

    sampler = lambda z: np.random.choice(z, 1, replace=True)

    correlation = np.array(sampler(correlation_))[0]
    ijv_loc     = np.array(sampler(ijv_loc_))[0]
    ijv_scale   = np.array(sampler(ijv_scal_))[0]
    cca_loc     = np.array(sampler(cca_loc_))[0]
    cca_scale   = np.array(sampler(cca_scale_))[0]

    print("correlation, ijv_loc, ijv_scale, cca_loc, cca_scale",
          correlation, ijv_loc, ijv_scale, cca_loc, cca_scale)
    y_hat = WarpedGaussianCopula(
        loc=[0., 0.],
        scale_tril=[[1., 0.], [correlation, tf.sqrt(1. - correlation ** 2)]],

        marginal_bijectors=[
            tfb.Invert(tfb.KumaraswamyCDF(ijv_loc, ijv_scale)),
            tfb.Invert(tfb.GumbelCDF(loc=cca_loc, scale=cca_scale))])

    return y_hat


data = np.array(value_lists)[:, 1:, 0]

unnormalized_posterior = lambda correlation, ijv_loc, ijv_scale, cca_loc, cca_scale: \
            joint_log_prob(data, correlation, ijv_loc, ijv_scale, cca_loc, cca_scale)


[correlation_, ijv_loc_, ijv_scale_, cca_loc_, cca_scale_] = run_chain(
    construct_hmc(unnormalized_posterior, adaptation_steps=NUM_ADAPTATION),
    inits=[0., 300., 200., 50., 20.],
    iters=[NUM_POSTERIOR_SAMPLES, NUM_BURNIN_ITERATIONS]
)