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

1 vue (au cours des 30 derniers jours)
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().

Réponse acceptée

Jeff Miller
Jeff Miller le 24 Avr 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)

Plus de réponses (0)

Community Treasure Hunt

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

Start Hunting!

Translated by