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;
Are you ready for the spotlight? We're accepting content ideas for SAS Innovate 2025 to be held May 6-9 in Orlando, FL. The call is open until September 25. Read more here about why you should contribute and what is in it for you!
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.