How do I declare correlation between a random variable with beta distribution and a random variable with normal distribution?

1 回表示 (過去 30 日間)
I have a random variable with beta distribution declared:
aux_alpha = 3;
aux_beta = 813.21;
LB = 41.368542; %lower bound
D = 1337.582858; %difference between upper bound and lower bound
alpha = aux_alpha+1;
beta = aux_beta+1;
aux_fy = betarnd(alpha,beta);
fy = aux_fy*D+LB;
And a random variable with normal distribution:
Ast_mean = Ast0*b*d;
Ast_std = 0.02*Ast_mean;
Ast = normrnd(Ast_mean,Ast_std);
There is a correlation of 0.5 between them. I'm trying to test it with the code:
nsam = 10000;
rhoij = 0.5;
y1 = normrnd(0,1,nsam,1);
y2 = normrnd(0,1,nsam,1);
RY = [1 rhoij;rhoij 1];
Jzy = chol(RY,'lower');
z1 = 0;
z2 = 0;
u1 = 0;
u2 = 0;
for sam = 1:nsam
zaux = Jzy*[y1(sam);y2(sam)];
z1(sam) = zaux(1);
z2(sam) = zaux(2);
u1(sam) = normcdf(z1(sam));
u2(sam) = normcdf(z2(sam));
x1_aux(sam) = betainv(u1(sam),alpha,beta);
x1(sam) = x1_aux(sam)*D+LB;
x2(sam) = norminv(u2(sam),Ast_mean,Ast_std);
end
rho = corrcoef(x1,x2)
But rho is not returning the correct correlation (0.5). Does anyone know how to solve this? I'm in doubt especially in the part where I use betainv().

採用された回答

Jeff Miller
Jeff Miller 2021 年 4 月 24 日
Something like this would work, though you would have to adjust rho by trial and error to get the final corr(x1,x2) value that you want:
% Just example values for some parameters--I was not sure what values you
% want.
alpha = 1.1;
beta = 4.2;
Ast_mean = 23.45;
Ast_std = 6.789;
LB = 41.368542; %lower bound
D = 1337.582858; %difference between upper bound and lower bound
nsam = 100000; % generate a large sample so we can estimate the true correlation accurately
rho = 0.5; % adjust this up or down to get the desired corr(x1,x2)=0.5 (or close enough)
% You shouldn't have to change anything from here
x = mvnrnd([0 0],[1 rho; rho 1],nsam);
xcdf = normcdf(x);
x1_aux = betainv(xcdf(:,1),alpha,beta);
x1 = x1_aux*D+LB;
x2 = x(:,2)*Ast_std + Ast_mean;
corr(x1,x2)

その他の回答 (0 件)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by