-
Notifications
You must be signed in to change notification settings - Fork 2
/
artigen.m
39 lines (31 loc) · 1022 Bytes
/
artigen.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
clear;
%% predefine
nr = 32;
nc = 32;
p = nr*nc;
n = 1000;
K = 10000;
%% determine the underlying Gaussian process
[Ltt,Lng] = meshgrid((1:nr),(1:nc));
Ltt = Ltt(:);
Lng = Lng(:);
Phi = 8;
S = 0.01*expcor(Ltt,Lng,Phi);
%% sample max-stable process
Z = zeros(n,p);
parfor i = 1:n
ee = exprnd(1,K,1);
pp = cumsum(ee);
W = mvnrnd(ones(1,p),S,K);
Z(i,:) = max(spdiags(1./pp,0,K,K)*W);
end
%% transform the unit Frechet marginals to general GEV marginals
nc = max(Lng)+0.5;
nr = max(Ltt)+0.5;
G=0.15; %0.05+0.001*(Lng-nr/2)-0.001*(Ltt-nc/2)+0.0012*((Lng-nr/2-1).^2-1.2^2)-0.0013*(Lng-nr/2+1).*(Ltt-nc/2-1)-0.0015*((Ltt-nc/2+1).^2-1.2^2);
S=15+0.005*(Lng-nr/2)+0.01*(Ltt-nc/2)-0.008*((Lng-nr/2-0.5).^2-7.5^2)-0.005*(Lng-nr/2-2).*(Ltt-nc/2+3)+0.007*((Ltt-nc/2+0.5).^2-1.5^2);
L=30+0.02*(Lng-nr/2+2)-0.03*(Ltt-nc/2-2)+0.02*((Lng-nr/2).^2-5.5^2)-0.04*(Lng-nr/2+0.5).*(Ltt-nc/2-1.5)-0.03*((Ltt-5).^2-2.5^2);
Ga = G; %repmat(G.',n,1);
Sa = repmat(S.',n,1);
La = repmat(L.',n,1);
XDat = La + Sa./Ga.*(Z.^Ga-1);