Main Content

코퓰러를 사용하여 종속 확률 변수 시뮬레이션하기

이 예제에서는 변수 사이에 복잡한 관계가 있거나 개별 변수의 분포가 각기 다른 경우 코퓰러를 사용하여 다변량 분포에서 데이터를 생성하는 방법을 보여줍니다.

MATLAB®은 임의의 입력값이나 잡음을 포함해서 시뮬레이션하기에 이상적인 툴입니다. Statistics and Machine Learning Toolbox™는 여러 일반적인 일변량 분포에 따라 난수 데이터 수열을 생성하는 함수를 제공합니다. 이 툴박스에는 다변량 정규분포 및 다변량 t 분포와 같은 다변량 분포에서 난수 데이터를 생성하는 몇 가지 함수도 포함되어 있습니다. 그러나 모든 주변 분포에 대한 다변량 분포를 생성하는 방법이나 각 변수의 분포가 다를 때 다변량 분포를 생성하는 방법은 제공하지 않습니다.

최근에는 코퓰러를 시뮬레이션 모델에 많이 사용하고 있습니다. 코퓰러는 변수 간의 종속성을 설명하는 함수이며, 상관관계가 있는 다변량 데이터를 모델링하기 위해 분포를 생성하는 방법을 제공합니다. 코퓰러를 사용하면 데이터 분석가가 주변 일변량 분포를 지정하고 변수 간 상관관계 구조를 제공하는 특정 코퓰러를 선택하여 다변량 분포를 생성할 수 있습니다. 이변량 분포뿐만 아니라 더 높은 차원의 분포도 가능합니다. 이 예제에서는 Statistics and Machine Learning Toolbox가 제공하는 코퓰러 함수를 사용하여 MATLAB에서 종속 다변량 난수 데이터를 생성하는 방법에 대해 설명합니다.

시뮬레이션 입력값 간의 종속성

몬테카를로 시뮬레이션을 설계할 때는 임의의 입력값에 대해 어떤 확률 분포를 선택할지 결정해야 합니다. 각각의 개별 변수에 대한 분포를 선택하는 것은 대개 간단하지만, 입력값 간에 어떠한 종속성이 있는지를 결정하는 것은 그렇지 않을 수 있습니다. 이상적으로는 시뮬레이션의 입력 데이터는 모델링하려는 실제 수치 간에 파악된 종속성을 반영해야 합니다. 그러나 시뮬레이션할 때 종속성을 파악할 수 있는 정보가 거의 또는 전혀 없을 수 있으며, 이러한 경우 모델의 민감도를 결정하기 위해 여러 다른 가능성을 시험해 보는 것이 좋습니다.

그렇지만 임의의 입력값이 표준 다변량 분포가 아닌 분포 형태를 띠는 경우 종속성이 있는 입력값을 실제로 생성하기가 어려울 수 있습니다. 더 나아가 표준 다변량 분포 중 일부는 매우 한정된 유형의 종속성만 모델링할 수 있습니다. 입력값을 독립 변수로 만드는 게 항상 가능하며 간단한 선택이긴 하지만 이 방법이 타당하지 않은 경우도 있고 잘못된 결론을 초래할 수도 있습니다.

예를 들어, 재무 위험에 대한 몬테카를로 시뮬레이션에는 보험 손실의 여러 요인을 나타내는 임의의 입력값이 있을 수 있습니다. 이러한 입력값은 로그 정규 확률 변수로 모델링될 수 있습니다. 이 경우 필요한 질문은 두 입력값 간의 종속성이 시뮬레이션 결과에 어떠한 영향을 미치는가입니다. 실제로 ,동일한 임의 조건이 두 요인에 모두 영향을 미친다는 것이 실제 데이터를 통해 확인될 경우 시뮬레이션에서 이를 무시하면 잘못된 결론을 초래할 수 있습니다.

독립 로그 정규 확률 변수에 대한 시뮬레이션은 간단합니다. 가장 간단한 방법은 lognrnd 함수를 사용하는 것입니다. 여기서는 mvnrnd 함수를 사용하여 독립 정규 확률 변수를 n개 쌍 생성한 후 거듭제곱하겠습니다. 여기에 사용되는 공분산 행렬은 대각 행렬로서 Z의 열은 서로 독립적입니다.

n = 1000;
sigma = .5;
SigmaInd = sigma.^2 .* [1 0; 0 1]
SigmaInd =

    0.2500         0
         0    0.2500

ZInd = mvnrnd([0 0], SigmaInd, n);
XInd = exp(ZInd);
plot(XInd(:,1),XInd(:,2),'.');
axis equal;
axis([0 5 0 5]);
xlabel('X1');
ylabel('X2');

종속 이변량 로그 정규 확률 변수도 0이 아닌 비대각선 항을 갖는 공분산 행렬을 사용하여 쉽게 생성할 수 있습니다.

rho = .7;
SigmaDep = sigma.^2 .* [1 rho; rho 1]
SigmaDep =

    0.2500    0.1750
    0.1750    0.2500

ZDep = mvnrnd([0 0], SigmaDep, n);
XDep = exp(ZDep);

두 번째 산점도 플롯은 이 두 이변량 분포 간의 차이를 보여줍니다.

plot(XDep(:,1),XDep(:,2),'.');
axis equal;
axis([0 5 0 5]);
xlabel('X1');
ylabel('X2');

두 번째 데이터 세트에서 X1의 큰 값이 X2의 큰 값과 연관되는 경향이 더 많으며 작은 값들도 마찬가지임을 명확히 알 수 있습니다. 이러한 종속성은 기본 이변량 정규분포의 상관 모수 rho로 결정됩니다. 시뮬레이션에서 도출되는 결론은 X1 및 X2가 종속성을 가지고 생성되었는지 여부에 따라 달라질 수 있습니다.

이 경우에는 이변량 로그 정규분포가 간단한 해결책이며, 더 고차원으로 일반화하거나 주변 분포가 서로 다른 로그 정규분포일 때 일반화하기도 쉽습니다. 다른 다변량 분포도 존재합니다. 예를 들어, 다변량 t 분포 및 디리클레 분포는 각각 종속 t 확률 변수 및 베타 확률 변수를 시뮬레이션하는 데 사용됩니다. 그러나 단순 다변량 분포는 그 목록이 길지 않으며, 주변 성분이 모두 동일한 분포군(또는 정확히 동일한 분포)에 있는 경우에만 적용됩니다. 이는 많은 경우 실질적인 제한이 될 수 있습니다.

종속 이변량 분포를 생성하기 위한 더욱 일반적인 방법

위에서 설명한 이변량 로그 정규분포 생성은 간단하지만 일반적으로 적용 가능한 방법을 보여줍니다. 먼저, 이변량 정규분포에서 값 쌍을 생성합니다. 이 두 변수 간에는 통계적 종속성이 있으며 각 변수는 정규 주변 분포를 가집니다. 다음으로, 변환(지수 함수)을 각 변수에 개별적으로 적용하여 주변 분포를 로그 정규분포로 변경합니다. 변환된 변수는 계속해서 통계적 종속성을 가집니다.

적합한 변환을 찾을 수 있는 경우 이 방법을 일반화하여 다른 주변 분포를 갖는 종속 이변량 확률 벡터를 생성할 수 있습니다. 실제로, 지수(거듭제곱)만큼 간단하지는 않지만 이러한 변환을 생성하는 일반적인 방법이 있습니다.

정의에 따르면, 정규 CDF(여기서는 PHI로 표시됨)를 표준 정규 확률 변수에 적용하면 구간 [0, 1]에서 균일한 확률 변수가 생성됩니다. 이를 확인하기 위해, Z가 표준 정규분포를 가지면 U = PHI(Z)의 CDF는 다음과 같으며

  Pr{U <= u0} = Pr{PHI(Z) <= u0} = Pr{Z <= PHI^(-1)(u0)} = u0,

이것이 U(0,1) 확률 변수의 CDF입니다. 시뮬레이션으로 생성된 정규 값과 변환된 값 일부를 히스토그램으로 그려보면 이 사실을 확인할 수 있습니다.

n = 1000;
z = normrnd(0,1,n,1);
hist(z,-3.75:.5:3.75);
xlim([-4 4]);
title('1000 Simulated N(0,1) Random Values');
xlabel('Z');
ylabel('Frequency');

u = normcdf(z);
hist(u,.05:.1:.95);
title('1000 Simulated N(0,1) Values Transformed to U(0,1)');
xlabel('U');
ylabel('Frequency');

이제, 일변량 난수 생성 이론을 빌려서 분포 F의 역 CDF를 U(0,1) 확률 변수에 적용하면 분포가 정확히 F인 확률 변수가 생성됩니다. 이를 역산법이라고 합니다. 이 증명은 근본적으로 앞 사례에서 설명한 증명과 반대입니다. 다음은 감마 분포로의 변환을 보여주는 또 다른 히스토그램입니다.

x = gaminv(u,2,1);
hist(x,.25:.5:9.75);
title('1000 Simulated N(0,1) Values Transformed to Gamma(2,1)');
xlabel('X');
ylabel('Frequency');

두 단계로 구성된 이 변환은 표준 이변량 정규분포의 각 변수에 적용될 수 있으며, 이 경우 임의 주변 분포를 갖는 종속 확률 변수가 생성됩니다. 변환이 각 성분에 대해 개별적으로 동작하기 때문에 결과로 생성되는 두 확률 변수가 동일한 주변 분포를 가지지 않아도 됩니다. 변환은 다음과 같이 정의됩니다.

  Z = [Z1 Z2] ~ N([0 0],[1 rho; rho 1])
  U = [PHI(Z1) PHI(Z2)]
  X = [G1(U1) G2(U2)]

여기서 G1 및 G2는 두 개의 서로 다를 수 있는 분포의 역 CDF입니다. 예를 들어, t(5) 및 Gamma(2,1) 주변 성분을 갖는 이변량 분포에서 확률 벡터를 생성할 수 있습니다.

n = 1000;
rho = .7;
Z = mvnrnd([0 0], [1 rho; rho 1], n);
U = normcdf(Z);
X = [gaminv(U(:,1),2,1) tinv(U(:,2),5)];

이 플롯은 주변 분포와 종속성을 모두 나타내기 위해 산점도 플롯 옆에 히스토그램을 함께 표시했습니다.

[n1,ctr1] = hist(X(:,1),20);
[n2,ctr2] = hist(X(:,2),20);
subplot(2,2,2);
plot(X(:,1),X(:,2),'.');
axis([0 12 -8 8]);
h1 = gca;
title('1000 Simulated Dependent t and Gamma Values');
xlabel('X1 ~ Gamma(2,1)');
ylabel('X2 ~ t(5)');
subplot(2,2,4);
bar(ctr1,-n1,1);
axis([0 12 -max(n1)*1.1 0]);
axis('off');
h2 = gca;
subplot(2,2,1);
barh(ctr2,-n2,1);
axis([-max(n2)*1.1 0 -8 8]);
axis('off');
h3 = gca;
h1.Position = [0.35 0.35 0.55 0.55];
h2.Position = [.35 .1 .55 .15];
h3.Position = [.1 .35 .15 .55];
colormap([.8 .8 1]);

순위 상관 계수

이렇게 생성된 X1과 X2 간의 종속성은 기본 이변량 정규분포의 상관 모수 rho에 의해 결정됩니다. 그러나, X1과 X2의 선형 상관이 rho는 아닙니다. 예를 들어, 원래 로그 정규 사례에서는 해당 상관에 대한 닫힌 형식이 있습니다.

  cor(X1,X2) = (exp(rho.*sigma.^2) - 1) ./ (exp(sigma.^2) - 1)

rho가 정확히 1인 경우를 제외하고 rho보다 엄밀히 작습니다. 위의 Gamma/t 생성처럼 더욱 일반적인 사례에서 X1과 X2 간의 선형 상관을 rho로 표현하기는 어렵거나 불가능합니다. 단, 시뮬레이션을 사용하면 동일한 효과가 발생하는 것을 나타낼 수 있습니다.

그 이유는 선형 상관 계수가 확률 변수 간의 선형 종속성을 나타내기 때문이며, 비선형 변환이 이 확률 변수에 적용된 경우 선형 상관이 유지되지 않습니다. 대신, 켄달의 타우(Kendall's tau)와 스피어만의 로(Spearman's rho)와 같은 순위 상관 계수가 더 적합합니다.

간단히 말해, 이러한 순위 상관은 한 확률 변수의 큰 값 또는 작은 값이 다른 확률 변수의 큰 값 또는 작은 값과 연관된 정도를 측정합니다. 그러나 선형 상관 계수와 다르게 오로지 순위의 측면에서 연관성을 측정합니다. 그 결과, 순위 상관이 모든 단조 변환에서 유지됩니다. 특히, 방금 설명한 변환 방법은 순위 상관을 유지합니다. 따라서, 이변량 정규분포 Z의 순위 상관을 파악하면 최종 변환된 확률 변수 X의 순위 상관을 정확히 확인할 수 있습니다. 기본 이변량 정규분포를 모수화하기 위해 rho가 여전히 필요하지만 확률 변수 간의 종속성을 설명할 때는 켄달의 타우나 스피어만의 로가 더 유용합니다. 이 둘은 주변 분포의 선택에 따라 변하지 않기 때문입니다.

이변량 정규분포의 경우 켄달의 타우 또는 스피어만의 로 간에 단순한 일대일 매핑이 있으며 선형 상관 계수 rho는 다음과 같습니다.

  tau   = (2/pi)*arcsin(rho)     or   rho = sin(tau*pi/2)
  rho_s = (6/pi)*arcsin(rho/2)   or   rho = 2*sin(rho_s*pi/6)
subplot(1,1,1);
rho = -1:.01:1;
tau = 2.*asin(rho)./pi;
rho_s = 6.*asin(rho./2)./pi;
plot(rho,tau,'-', rho,rho_s,'-', [-1 1],[-1 1],'k:');
axis([-1 1 -1 1]);
xlabel('rho');
ylabel('Rank correlation coefficient');
legend('Kendall''s tau', 'Spearman''s rho_s', 'location','northwest');

따라서, Z1과 Z2 간의 선형 상관에 올바른 rho 모수 값을 선택하여 주변 분포에 관계없이 X1과 X2 간에 원하는 순위 상관을 쉽게 생성할 수 있습니다.

참고로, 다변량 정규분포의 경우 스피어만의 순위 상관이 선형 상관과 거의 동일합니다. 하지만, 최종 확률 변수로 변환한 후에는 그렇지 않습니다.

코퓰러

위에서 설명한 생성의 첫 번째 단계에서는 코퓰러, 특히 가우스 코퓰러라고 하는 분포를 정의합니다. 이변량 코퓰러는 간단히 말해서 각각 균일한 주변 분포를 가지는 두 개의 확률 변수에 대한 확률 분포입니다. 이 두 변수는 완전히 독립적이거나, 결정적으로 연관되어 있거나(즉, U2 = U1), 아니면 이 둘 사이의 중간일 수 있습니다. 이변량 가우스 코퓰러군은 선형 상관 행렬인 Rho = [1 rho; rho 1]로 모수화됩니다. rho가 +/- 1에 가까워질수록 U1과 U2는 선형 종속성을 띠게 되고 rho가 0에 가까워질수록 완전한 독립성을 띠게 됩니다.

다양한 수준의 rho에 대해 시뮬레이션된 난수 값을 산점도 플롯으로 그려보면 여러 범위의 가우스 코퓰러가 가능함을 알 수 있습니다.

n = 500;
Z = mvnrnd([0 0], [1 .8; .8 1], n);
U = normcdf(Z,0,1);
subplot(2,2,1);
plot(U(:,1),U(:,2),'.');
title('rho = 0.8');
xlabel('U1');
ylabel('U2');
Z = mvnrnd([0 0], [1 .1; .1 1], n);
U = normcdf(Z,0,1);
subplot(2,2,2);
plot(U(:,1),U(:,2),'.');
title('rho = 0.1');
xlabel('U1');
ylabel('U2');
Z = mvnrnd([0 0], [1 -.1; -.1 1], n);
U = normcdf(Z,0,1);
subplot(2,2,3);
plot(U(:,1),U(:,2),'.');
title('rho = -0.1');
xlabel('U1');
ylabel('U2');
Z = mvnrnd([0 0], [1 -.8; -.8 1], n);
U = normcdf(Z,0,1);
subplot(2,2,4);
plot(U(:,1),U(:,2),'.');
title('rho = -0.8');
xlabel('U1');
ylabel('U2');

U1과 U2 간의 종속성은 주변 분포 X1 = G(U1) 및 X2 = G(U2)와는 완전히 별개입니다. X1 및 X2에는 임의의 주변 분포가 주어질 수 있으며, 여전히 동일한 순위 상관을 가집니다. 이는 코퓰러의 주요 장점 중 하나입니다. 즉, 코퓰러에서는 이렇게 종속성과 주변 분포를 별도로 지정할 수 있습니다.

t 코퓰러

이변량 t 분포에서 시작하고 그에 대응하는 t CDF를 사용하여 변환을 수행함으로써 다른 코퓰러군을 생성할 수 있습니다. 이변량 t 분포는 선형 상관 행렬 Rho와 자유도 nu로 모수화됩니다. 따라서 예를 들어, t(1) 또는 t(5) 코퓰러를 각각 자유도가 1과 5인 다변량 t를 기반으로 하여 나타낼 수 있습니다.

다양한 수준의 rho에 대해 시뮬레이션된 난수 값을 산점도 플롯으로 그려보면 여러 범위의 t(1) 코퓰러가 가능함을 알 수 있습니다.

n = 500;
nu = 1;
T = mvtrnd([1 .8; .8 1], nu, n);
U = tcdf(T,nu);
subplot(2,2,1);
plot(U(:,1),U(:,2),'.');
title('rho = 0.8');
xlabel('U1');
ylabel('U2');
T = mvtrnd([1 .1; .1 1], nu, n);
U = tcdf(T,nu);
subplot(2,2,2);
plot(U(:,1),U(:,2),'.');
title('rho = 0.1');
xlabel('U1');
ylabel('U2');
T = mvtrnd([1 -.1; -.1 1], nu, n);
U = tcdf(T,nu);
subplot(2,2,3);
plot(U(:,1),U(:,2),'.');
title('rho = -0.1');
xlabel('U1');
ylabel('U2');
T = mvtrnd([1 -.8; -.8 1], nu, n);
U = tcdf(T,nu);
subplot(2,2,4);
plot(U(:,1),U(:,2),'.');
title('rho = -0.8');
xlabel('U1');
ylabel('U2');

t 코퓰러는 가우스 코퓰러와 마찬가지로 U1 및 U2에 대해 균일한 주변 분포를 가집니다. t 코퓰러의 성분 간 순위 상관 tau 또는 rho_s도 가우스 코퓰러와 동일한 rho의 함수입니다. 그러나 이 플롯에서 볼 수 있듯이 t(1) 코퓰러는 성분이 동일한 순위 상관을 갖는 경우에도 가우스 코퓰러와 상당히 다릅니다. 그 차이는 각각의 종속성 구조에 있습니다. 자유도 모수 nu가 커지면 t(nu) 코퓰러는 그에 대응하는 가우스 코퓰러에 가까워집니다.

가우스 코퓰러와 마찬가지로, t 코퓰러에도 주변 분포를 적용할 수 있습니다. 예를 들어, 자유도가 1인 t 코퓰러를 사용하여 Gamma(2,1) 및 t(5) 주변 성분을 갖는 이변량 분포에서 확률 벡터를 다시 생성할 수 있습니다.

subplot(1,1,1);
n = 1000;
rho = .7;
nu = 1;
T = mvtrnd([1 rho; rho 1], nu, n);
U = tcdf(T,nu);
X = [gaminv(U(:,1),2,1) tinv(U(:,2),5)];

[n1,ctr1] = hist(X(:,1),20);
[n2,ctr2] = hist(X(:,2),20);
subplot(2,2,2);
plot(X(:,1),X(:,2),'.');
axis([0 15 -10 10]);
h1 = gca;
title('1000 Simulated Dependent t and Gamma Values');
xlabel('X1 ~ Gamma(2,1)');
ylabel('X2 ~ t(5)');
subplot(2,2,4);
bar(ctr1,-n1,1);
axis([0 15 -max(n1)*1.1 0]);
axis('off');
h2 = gca;
subplot(2,2,1);
barh(ctr2,-n2,1);
axis([-max(n2)*1.1 0 -10 10]);
axis('off');
h3 = gca;
h1.Position = [0.35 0.35 0.55 0.55];
h2.Position = [.35 .1 .55 .15];
h3.Position = [.1 .35 .15 .55];
colormap([.8 .8 1]);

앞에서 가우스 코퓰러를 기반으로 생성한 이변량 Gamma/t 분포와 비교했을 때, 여기서 t(1) 코퓰러를 기반으로 생성한 분포는 동일한 주변 분포를 가지며 변수 간 순위 상관도 동일하지만 종속성 구조는 매우 다릅니다. 이는 다변량 분포가 오로지 주변 분포나 상관관계만으로 정의되지 않는다는 사실을 보여줍니다. 실제 관측된 데이터를 기반으로 하여 적용할 특정 코퓰러를 선택하거나, 입력값 분포에 대한 시뮬레이션 결과의 민감도를 결정하기 위해 서로 다른 여러 코퓰러를 사용할 수 있습니다.

고차 코퓰러

가우스 코퓰러와 t 코퓰러를 타원형 코퓰러라고 합니다. 타원형 코퓰러는 더 높은 차원 수로 쉽게 일반화할 수 있습니다. 예를 들어, 다음과 같이 가우스 코퓰러를 사용하여 Gamma(2,1), Beta(2,2), t(5) 주변 성분을 갖는 삼변량 분포에서 데이터를 시뮬레이션할 수 있습니다.

subplot(1,1,1);
n = 1000;
Rho = [1 .4 .2; .4 1 -.8; .2 -.8 1];
Z = mvnrnd([0 0 0], Rho, n);
U = normcdf(Z,0,1);
X = [gaminv(U(:,1),2,1) betainv(U(:,2),2,2) tinv(U(:,3),5)];
plot3(X(:,1),X(:,2),X(:,3),'.');
grid on;
view([-55, 15]);
xlabel('U1');
ylabel('U2');
zlabel('U3');

선형 상관 모수 rho와 예를 들어, 켄달의 타우 간의 관계는 여기에 사용된 상관 행렬 Rho의 각 항목에 적용됩니다. 데이터의 표본 순위 상관이 이론적 값과 대략 같은지 확인할 수 있습니다.

tauTheoretical = 2.*asin(Rho)./pi
tauTheoretical =

    1.0000    0.2620    0.1282
    0.2620    1.0000   -0.5903
    0.1282   -0.5903    1.0000

tauSample = corr(X, 'type','Kendall')
tauSample =

    1.0000    0.2655    0.1060
    0.2655    1.0000   -0.6076
    0.1060   -0.6076    1.0000

코퓰러 및 경험적 주변 분포

코퓰러를 사용하여 종속 다변량 데이터를 시뮬레이션하려면 다음을 지정해야 한다는 것을 확인했습니다.

  1) the copula family (and any shape parameters),
  2) the rank correlations among variables, and
  3) the marginal distributions for each variable

두 개의 주식수익률 데이터 세트가 있으며 이 데이터와 동일한 분포를 따르는 입력값을 사용해서 몬테카를로 시뮬레이션을 실행한다고 가정하겠습니다.

load stockreturns
nobs = size(stocks,1);
subplot(2,1,1);
hist(stocks(:,1),10);
xlabel('X1');
ylabel('Frequency');
subplot(2,1,2);
hist(stocks(:,2),10);
xlabel('X2');
ylabel('Frequency');

(이 두 데이터 벡터는 길이가 같지만 이는 중요하지 않습니다.)

각 데이터 세트별로 모수적 모델을 피팅하고 이렇게 얻은 추정값을 주변 분포로 사용할 수 있습니다. 하지만, 모수적 모델이 충분히 유연하지 않을 수 있습니다. 대신, 주변 분포에 경험적 모델을 사용할 수 있습니다. 역 CDF를 계산할 수만 있으면 됩니다.

이 데이터 세트에 대한 경험적 역 CDF는 1/nobs, 2/nobs, ...1 값에 계단이 있는 계단 함수입니다. 1. 계단 높이는 정렬된 데이터입니다.

invCDF1 = sort(stocks(:,1));
n1 = length(stocks(:,1));
invCDF2 = sort(stocks(:,2));
n2 = length(stocks(:,2));
subplot(1,1,1);
stairs((1:nobs)/nobs, invCDF1,'b');
hold on;
stairs((1:nobs)/nobs, invCDF2,'r');
hold off;
legend('X1','X2');
xlabel('Cumulative Probability');
ylabel('X');

시뮬레이션을 위해 여러 다른 코퓰러와 상관관계를 실험해볼 수 있습니다. 여기서는 상당히 큰 음의 상관 모수와 함께 이변량 t(5) 코퓰러를 사용하겠습니다.

n = 1000;
rho = -.8;
nu = 5;
T = mvtrnd([1 rho; rho 1], nu, n);
U = tcdf(T,nu);
X = [invCDF1(ceil(n1*U(:,1))) invCDF2(ceil(n2*U(:,2)))];

[n1,ctr1] = hist(X(:,1),10);
[n2,ctr2] = hist(X(:,2),10);
subplot(2,2,2);
plot(X(:,1),X(:,2),'.');
axis([-3.5 3.5 -3.5 3.5]);
h1 = gca;
title('1000 Simulated Dependent Values');
xlabel('X1');
ylabel('X2');
subplot(2,2,4);
bar(ctr1,-n1,1);
axis([-3.5 3.5 -max(n1)*1.1 0]);
axis('off');
h2 = gca;
subplot(2,2,1);
barh(ctr2,-n2,1);
axis([-max(n2)*1.1 0 -3.5 3.5]);
axis('off');
h3 = gca;
h1.Position = [0.35 0.35 0.55 0.55];
h2.Position = [.35 .1 .55 .15];
h3.Position = [.1 .35 .15 .55];
colormap([.8 .8 1]);

시뮬레이션된 데이터의 주변 히스토그램은 원래 데이터와 거의 일치하며 더 많은 값 쌍을 시뮬레이션할수록 더 일치하게 됩니다. 값을 원래 데이터에서 추출했고 각 데이터 세트에 관측값이 단지 100개 있었기 때문에 시뮬레이션된 데이터가 다소 "이산적"임을 알 수 있습니다. 이를 해결하는 한 방법은 시뮬레이션된 최종 값에 가능하면 정규분포된 불규칙 변동을 약간 추가하는 것입니다. 이는 경험적 역 CDF의 평활화된 버전을 사용하는 것과 동일합니다.