Hi,
I think it again. You are right. It is all about Covariate matrix of mahalanobis. Which is totally different between Matlab and SAS. Eevn more the function COV() is also a little different between Matlab and SAS. If you want stick with SAS's function , see the following code and it will lead you to the similiar result with Matlab .
Data have1;
Input obs reps id $ x y;
Datalines ;
1
1
00-01
5
3
2
1
00-02
4
6
3
1
00-03
6
4
4
1
00-06
4
6
5
2
00-02
4
6
6
2
00-03
6
4
7
2
00-04
7
5
8
2
00-05
9
2
9
3
00-01
5
3
10
3
00-02
4
6
11
3
00-05
9
2
12
3
00-10
7
4
13
4
00-03
6
4
14
4
00-05
9
2
15
4
00-09
7
3
16
4
00-10
7
4
17
5
00-03
6
4
18
5
00-04
7
5
19
5
00-07
5
8
20
5
00-09
7
3
21
6
00-02
4
6
22
6
00-03
6
4
23
6
00-05
9
2
24
6
00-10
7
4
;
run;
Data have2;
Input obs reps id $ x y;
Datalines ;
25
7
00-01
5
3
26
7
00-08
1
8
27
7
00-09
7
3
28
7
00-10
7
4
;
Run;
proc iml;
use have1 nobs nobs;
read all var {reps};
read all var {x y} into data;
close;
use have2;
read all var {x y} into B;
close;
center=B[:,];
cB=cov(B);
start_end=t(loc(t(reps)^={.}||remove(reps,nobs)))||
t(loc(t(reps)^=remove(reps,1)||{.}));
want=j(nrow(start_end),2,.);
want[,1]=reps[start_end[,1]];
do i=1 to nrow(start_end);
idx=start_end[i,1]:start_end[i,2];
A=data[idx,];
mA=A[:,];
cA=cov(A);
cov=(nrow(A)/(nrow(A)+nrow(B)))#cA +
(nrow(B)/(nrow(A)+nrow(B)))#cB ;
want[i,2]=mahalanobis(mA,center,cov);
end;
create want from want[c={reps distance}];
append from want;
close;
quit;
Sample data ?
OK. Your data is not readable .so I used my original data.
SAS's function code :
Data have1;
Input Sampleid reps id $ x y;
Datalines ;
1 1 00-01 5 3
1 1 00-02 4 6
1 1 00-03 6 4
1 1 00-06 4 6
1 2 00-02 4 6
1 2 00-03 6 4
1 2 00-04 7 5
1 2 00-05 9 2
1 3 00-01 5 3
1 3 00-02 4 6
1 3 00-05 9 2
1 3 00-10 7 4
1 4 00-03 6 4
1 4 00-05 9 2
1 4 00-09 7 3
1 4 00-10 7 4
1 5 00-03 6 4
1 5 00-04 7 5
1 5 00-07 5 8
1 5 00-09 7 3
1 6 00-02 4 6
1 6 00-03 6 4
1 6 00-05 9 2
1 6 00-10 7 4
2 1 00-01 5 3
2 1 00-02 4 6
2 1 00-03 6 4
2 1 00-06 4 6
2 2 00-02 4 6
2 2 00-03 6 4
2 2 00-04 7 5
2 2 00-05 9 2
;
run;
Data have2;
Input obs reps id $ x y;
Datalines ;
25
7
00-01
5
3
26
7
00-08
1
8
27
7
00-09
7
3
28
7
00-10
7
4
;
Run;
proc iml;
use have1 nobs nobs;
read all var {Sampleid reps};
read all var {x y} into data;
close;
use have2;
read all var {x y} into B;
close;
group=catx(' ',Sampleid,reps);
center=B[:,];
cB=cov(B);
start_end=t(loc(t(group)^={' '}||remove(group,nobs)))||
t(loc(t(group)^=remove(group,1)||{' '}));
want=j(nrow(start_end),3,.);
want[,1]=Sampleid[start_end[,1]];
want[,2]=reps[start_end[,1]];
do i=1 to nrow(start_end);
idx=start_end[i,1]:start_end[i,2];
A=data[idx,];
mA=A[:,];
cA=cov(A);
cov=(nrow(A)/(nrow(A)+nrow(B)))#cA +
(nrow(B)/(nrow(A)+nrow(B)))#cB ;
want[i,3]=mahalanobis(mA,center,cov);
end;
create want from want[c={Sampleid reps distance}];
append from want;
close;
quit;
Matlab's function code :
proc iml;
use have1 nobs nobs;
read all var {Sampleid reps};
read all var {x y} into data;
close;
use have2;
read all var {x y} into B;
close;
start new_cov(x);
Xc=x-x[:,];
c=Xc`*Xc/nrow(x);
return (c);
finish;
start new_mahalanobis(A,B);
xDiff=A[:,]-B[:,];
cA=new_cov(A);
cB=new_cov(B);
pC=(nrow(A)/(nrow(A)+nrow(B)))#cA +
(nrow(B)/(nrow(A)+nrow(B)))#cB ;
d=sqrt(xDiff*inv(pC)*xDiff`);
return (d);
finish;
group=catx(' ',Sampleid,reps);
start_end=t(loc(t(group)^={' '}||remove(group,nobs)))||
t(loc(t(group)^=remove(group,1)||{' '}));
want=j(nrow(start_end),3,.);
want[,1]=Sampleid[start_end[,1]];
want[,2]=reps[start_end[,1]];
do i=1 to nrow(start_end);
idx=start_end[i,1]:start_end[i,2];
A=data[idx,];
want[i,3]=new_mahalanobis(A,B);
end;
create want from want[c={Sampleid reps distance}];
append from want;
close;
quit;
Join us for SAS Innovate April 16-19 at the Aria in Las Vegas. Bring the team and save big with our group pricing for a limited time only.
Pre-conference courses and tutorials are filling up fast and are always a sellout. Register today to reserve your seat.
Learn how to run multiple linear regression models with and without interactions, presented by SAS user Alex Chaplin.
Find more tutorials on the SAS Users YouTube channel.