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

1 Ansicht (letzte 30 Tage)
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().

Akzeptierte Antwort

Jeff Miller
Jeff Miller am 24 Apr. 2021
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)

Weitere Antworten (0)

Community Treasure Hunt

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

Start Hunting!

Translated by