I want to R random draws in Matlab from the following 3D cumulative distribution function (CDF) [This CDF belongs to the GEV family]
To do so, I have tried to adapt this code to my question as follows:
F = @(x0,x1,x2) exp(-exp(-x0) - (exp(-x1/lambda)+exp(-x2/lambda)).^lambda);
lims = [-5,10];
draws = zeros(R,3);
for r=1:R
draws(r,:) = sample_random(F,3,lims);
end
with auxiliary functions
function r = sample_random(F,N,lims)
delta = diff(lims)/10000;
x = linspace(lims(1),lims(2),300);
r = inf(1,N);
for ii = 1:N
marginal = get_marginal(F,r,ii,x,delta);
p = rand * marginal(end);
[~,I] = unique(marginal); % interp1 cannot handle duplicated points, let's remove them
r(ii) = interp1(marginal(I),x(I),p);
end
function marginal = get_marginal(F,r,ii,x,delta)
N = length(r);
marginal = zeros(size(x));
for jj=0:2^(ii-1)-1
rr = flip(dec2bin(jj,N)-'0');
sign = mod(sum(rr,2),2);
if sign == 0
sign = 1;
else
sign = -1;
end
args = num2cell(r - delta * rr);
args{ii} = x;
marginal = marginal + sign * F(args{:});
end
However, when I plot the empirical marginal CDF of the first component of F for different values of lambda, I get almost identical pictures, which I think should not be the case:
clear
rng default
lambda_grid=[0.05 0.4 1];
slambda=size(lambda_grid,2);
store_draws=cell(slambda,1);
R=10^3; %number of draws
for j=1:slambda
lambda=lambda_grid(j);
F = @(x1,x2,x0) exp(-exp(-x0) - (exp(-x1/lambda)+exp(-x2/lambda)).^lambda); %CDF of epsilon0, epsilon1, epsilon2
lims = [-5,10];
draws_temp= zeros(R,3);
for r=1:R
draws_temp(r,:) = sample_random(F,3,lims);
end
store_draws{j}=draws_temp;
end
ecdf(store_draws{1}(:,1))
hold on
ecdf(store_draws{2}(:,1))
hold on
ecdf(store_draws{3}(:,1))
legend('\lambda=0.1', '\lambda=0.5', '\lambda=1')
Therefore, I think I'm doing something wrong. Could you advise?
