title 'Classical MDS with PROC IML'; data flying(label=intercity flying mileage); input city $1-13 @16 (d1-d10) (10*6.); array d d1-d10; do over d; if d=. then d=0.; end; CARDS; Atlanta 0 Chicago 587 0 Denver 1212 920 0 Houston 701 940 879 0 Los Angeles 1936 1745 831 1374 0 Miami 604 1188 1726 968 2339 0 New York 748 713 1631 1420 2451 1092 0 San Francisco 2139 1858 949 1645 347 2594 2571 0 Seattle 2182 1737 1021 1891 959 2734 2408 678 0 Washington DC 543 597 1494 1220 2300 923 205 2442 2329 0 proc print; proc iml; reset autoname ; use flying ; read all into dist[r=city] ; dist = dist + dist`; * make it symmetric; ndim=2; *--------------cmdsr---------------------+ | | classical mds with ratio data | | input: | dist(n,n) distances | ndim(1,1) # of dimensions | | output: | coor(n,ndim) coordinates | | routines required: | gaussnew | +-----------------------------------------; CMDSR: n= nrow(dist); n2=(n*n-1)/2; ncoor=n*ndim; *-----torgerson mds-----; d=dist##2; d= d[,:]* j(1,n,1)+ j(n,1,1)* d[:,] - d[:,:]* j(n,n,1)-d; call eigen( x, coor, d); free d; coor= coor[,1:ndim] * diag( sqrt( x[1:ndim,])); coor=coor- j( nrow(coor),1,1)* coor[:,]; print "Initial Coordinates"; print coor; *-----initialize other matrices-----; beta= shape(coor,0,1); residual= j(n2,1,0); predict= j(n2,1,0); *-----GAUSS-NEWTON-----; run gaussnew; *-----rotate to principal components-----; coor= shape(beta,0,ndim); coor=coor- j( nrow(coor),1,1)* coor[:,]; call svd( coor, tmp1, tmp2, coor); coor=coor* diag(tmp1); free tmp1 tmp2; print " Final Coordinates"; print coor[r=city]; return; *--------COMPUTE DERIVATIVES-------; start deriv; coor= shape(beta,0,ndim); jacobian= j(n2,ncoor,0); obs=0; do row=2 to n; do col=1 to row-1; obs=obs+1; xx=( coor[col,]- coor[row,])/ predict[obs,1]; xx=xx||(-xx); index=(ndim*(col-1))+(1:ndim) ||(ndim*(row-1))+(1:ndim); jacobian[obs,index]=xx; end; end; finish; *----------END CMDSR--------------; *-------EVALUATE RESIDUALS------; start resid; coor= shape(beta,0,ndim); obs=0; do row=2 to n; do col=1 to row-1; obs=obs+1; predict[obs,1]= sqrt( ssq( coor[row,] - coor[col,])); residual[obs,1]= dist[row,col] - predict[obs,1]; end; end; finish; start GAUSSNEW; *------------------------------------------+ | | let n = # observations | p = # of parameters | | input: | beta(p,1) initial estimates | | output: | beta(p,1) final estimates | | routines required: | resid returns residual(n,1) | deriv returns jacobian(n,p) +-----------------------------------------; print "Gauss-Newton Nonlinear Regression"; run resid; sse= ssq(residual); bprint=(0||sse`); eps=1; *-----gauss-newton iterations-----; do iter=1 to 20 while(eps>{1e-3}); lastsse=sse; run deriv; gradient=jacobian`*residual; hessian=jacobian`*jacobian; delta= sweep(hessian,1: nrow(beta)) *gradient; beta=beta+delta; run resid; sse= ssq(residual); bprint=bprint//(iter||sse`); eps=(lastsse-sse)/lastsse; end; if ( iter>20) then print "*** iteration limit exceeded ***"; *-----print out results-----; names={'Iteration' 'SSE'}; print bprint[ colname=names]; finish; *-----END GAUSSNEW-----; *-------- SAMPLE DATA --------- CLASSICAL MDS BY GAUSS-NEWTON INTERCITY FLYING MILEAGES COOR COL1 COL2 ROW1 1016.48 -202.224 ROW2 540.308 482.02 ROW3 -681.089 35.7584 ROW4 228.348 -810.019 ROW5 -1702.34 -551.685 ROW6 1603.05 -822.941 ROW7 1516.37 734.011 ROW8 -2009.04 -159.225 ROW9 -1897.48 819.875 ROW10 1385.39 474.43 *------------------------------;