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;
It's finally time to hack! Remember to visit the SAS Hacker's Hub regularly for news and updates.