-
Notifications
You must be signed in to change notification settings - Fork 0
/
data_res.m
60 lines (54 loc) · 1.67 KB
/
data_res.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
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
function data_res(data,CI,nrsamples,resolutions,blocksize,method)
nr_vars=size(data,2);
nr_res=length(resolutions);
upper=zeros(nr_res,nr_vars);
lower=zeros(nr_res,nr_vars);
medians=zeros(nr_res,nr_vars);
for k=1:nr_res
res = resolutions(k);
starts=randi([1 2000001-res*blocksize],1,nrsamples);
MAFs=zeros(nrsamples,nr_vars);
for i=1:nrsamples
data2=data(starts(i):res:starts(i)+(blocksize-1)*res,:);
if method == 'MAF'
Wmaf=MAF(data2);
if Wmaf(1,1)<0
Wmaf(:,1)=Wmaf(:,1)*-1;
end
MAFs(i,:)=Wmaf(:,1);
elseif method == 'PCA'
C=cov(data2);
[V E]=eig(C);
d=diag(E);
idx=find(d==max(d));
PC=V(:,idx);
if PC(1)<0
PC=PC*-1;
end
MAFs(i,:)=PC;
end
end
for j=1:nr_vars
A=sort(MAFs(:,j));
upper(k,j)=A(CI(2));
lower(k,j)=A(CI(1));
medians(k,j)=median(A);
end
end
figure('pos',[100 100 430 300])
hold on
clrs=['r','g','b','k','c'];
clrs=clrs(1:nr_vars);
x=1:length(resolutions);
for k=1:nr_vars
fill([resolutions flip(resolutions)],[lower(:,k)' flip(upper(:,k)')],clrs(k),'LineStyle','none','facealpha',0.15)
plot(resolutions,medians(:,k),'Color',clrs(k))
end
ylim([-1 1])
xlabel('sampling distance')
if method=='MAF'
ylabel('MAF')
elseif method=='PCA'
ylabel('PC')
end
end