iml编程文档
- 格式：doc
- 大小：923.50 KB
- 文档页数：89
ＳＡＳ矩阵运算语言 IML 基础(Interactive Matrix Language) 第一节，矩阵与ＩＭＬ编程 一，矩阵ＩＭＬ定义 １行向量语句
例 x={1 2 3 4 5 6}; 等价于行向量 matrix X:
X 1 2 3 4 5 6 ２列向量语句
y={1,2,3,4,5}; 等价于列向量 matrix Y:
Y 1 2 3 4 5 ３一般矩阵定义
z={1 2, 3 4, 5 6}; 等价于一个 3 ×2 matrix Z:
Z 1 2 3 4 5 6 a={3 -1 2, 2 -2 3, 4 1 -4}; 等价于一个 3 ×３ matrix Ａ: 312223414A
４，生成a,b,…,c为对角元素的对角阵与N阶单位阵 ＢＬＯＣＫ（a,b,…,c） Ｉ（n） PRON IML; RESET PRINT;
B=ＢＬＯＣＫ（1,2,3） C=I(3) PRINT B C;
100020003B 100010001C
5．矩阵各元素仅取同一数值：J（nrow,ncol,数值）
RESET PRINT; A=J(3,3,5） C=J(3,3) PRINT A C;
555555555B 111111111C
The DESIGNF Function 用设计函数DESIGNF建立设计矩阵 DESIGNF( 向量值 );
例为作方差分析作一设计，其中有三水平１，２，３．而它们的容量分别为 n1=3, n2=2, 及 n3=2 d=designf({1,1,1,2,2,3,3});
D 7 rows 2 cols (numeric)
1 0 1 0 1 0 0 1 0 1 -1 -1 -1 -1 二，矩阵运算 加 ＋，减 －，乘 ＊，逆 inv(a)，转置 t(a)，对应元素相乘 ＃． ## ：各元素平方，＜＞ ：对应元素取最大值， ＜＝ ：Ａ的元素满足＜＝Ｂ的对应元素取１，否则取０
例，设有矩阵 A 和 B .
则: 则 (# ：对应元素相乘) : 则 (## ：各元素平方) : 则 (＜＞ ：对应元素取最大值) 则 (＜＝ ：Ａ的元素满足＜＝Ｂ的对应元素取１，否则取０) Ａ的元素与３相除的余数组成的矩阵: 数乘运算k#Ａ与转置运算t(A) 设Ｃ为一个 3 ×2矩阵:
12242610342#68()4812561012TCCtCC则:
x=x#(x>0); Ｘ的元素满足＞０的对应元素取１，否则取０ , 可使用 RESET: PRINTALL 显示
求逆运算inv(A)
3120.0672231.3331.3330.3334140.6670.4672.267A-10.333-0.133则inv(A)==A
proc iml; a={3 -1 2, 2 -2 3, 4 1 -4}; c={8, 2, 9}; x=inv(a); PRINT X; quit;
三，用逆矩阵法求解方程组运算 Ax = c
x = A-1c 例
表为矩阵方程型式 proc iml; a={3 -1 2, 2 -2 3, 4 1 -4}; c={8, 2, 9}; x=inv(a)*c; quit; X 3 rows 1 col (numeric)
3 5 2
四．线性回归 y = Xb + e 其最小二乘法解为： b = (X’X)-1x’ y sas中已有 REG 过程, 现用 PROC IML 编程求之.
proc iml; x={1 1 1, 1 2 4, 1 3 9, 1 4 16, 1 5 25}; y={1,5,9,23,36}; b=inv(t(x)*x)*t(x)*y; yhat=x*b; r=y-yhat; sse=ssq(r); dfe=nrow(x)-ncol(x); mse=sse/dfe;
proc iml; x={1 1 1, 1 2 4, 1 3 9, 1 4 16, 1 5 25}; y={1,5,9,23,36};
start regress; /* 模块开始 */ xpxi=inv(t(x)*x); /* 求 X'X 的逆阵 */ beta=xpxi*(t(x)*y); /* 对回归系数作估计 */ yhat=x*beta; /* 计算预测值 */ resid=y-yhat; /* 求残差 */ sse=ssq(resid); /* 求残差平方和SSE */ n=nrow(x); /* 求容量 */ dfe=nrow(x)-ncol(x); /* ＳＳＥ的自由度 */ mse=sse/dfe; /*求残差均平方和MSE */ ssa=ssq(yhat-sum(y)/n); /* 回归离差 SSa */ cssy=ssq(y-sum(y)/n); /* 总离差 SSＴ */ Ｆ=(cssy-sse); /*求Ｆ值 */ fprob=1-probf(x,ncol(x),dfe); /* p-values */ rsquare=(cssy-sse)/cssy; /*求确定系数值 */
print sse dfe mse fprob rsquare; stdb=sqrt(vecdiag(xpxi)*mse); /* 标准差std 的估计*/
t=beta/stdb; /* 参数的t值 */ prob=1-probf(t#t,1,dfe); /* p-values */ print,"Parameter Estimates",, beta stdb t prob; print,y yhat resid;
finish regress; /* 模块结束 */ reset noprint; run regress; /* 运行模块 */
reset print; /* 进入自动显示 */
covb=xpxi*mse; /* 协方差阵估计 */ s=1/sqrt(vecdiag(covb)); corrb=diag(s)*covb*diag(s); /*相关系数阵 */ run regress;
x1={1,2,3,4,5}; x=orpol(x1,2); reset noprint; run regress; reset print; covb=xpxi*mse; s=1/sqrt(vecdiag(covb)); corrb=diag(s)*covb*diag(s); reset fuzz; corrb=diag(s)*covb*diag(s);
proc iml; start regress; /* Ä£¿é¿ªÊ¼ */ reset noprint; xpxi=inv(t(x)*x); /* Çó X'X µÄÄæÕó */ beta=xpxi*(t(x)*y); /* ¶Ô»Ø¹éϵÊý×÷¹À¼Æ */ yhat=x*beta; /* ¼ÆËãÔ¤²âÖµ */ resid=y-yhat; /* Çó²Ð²î */ sse=ssq(resid); /* Çó²Ð²îƽ·½ºÍSSE */ n=nrow(x); /* ÇóÈÝÁ¿ */ k=ncol(x); dfe=nrow(x)-ncol(x); /* £Ó£Ó£ÅµÄ×ÔÓÉ¶È */ mse=sse/dfe; /*Çó²Ð²î¾ùƽ·½ºÍMSE */ ssa=ssq(yhat-sum(y)/n); /* »Ø¹éÀë²î SSa */ msa=ssa/k; cssy=ssq(y-sum(y)/n); /* ×ÜÀë²î SS£Ô */ F=msa/mse;/*Çó£ÆÖµ */ fprob=1-probf(F#F,k,dfe-1); /* p-values */ rsquare=(cssy-sse)/cssy; /*ÇóÈ·¶¨ÏµÊýÖµ */
print sse dfe mse F fprob rsquare; stdb=sqrt(vecdiag(xpxi)*mse); /* ±ê×¼²îstd µÄ¹À¼Æ*/ t=beta/stdb; /* ²ÎÊýµÄtÖµ */ prob=1-probf(t#t,1,dfe); /* p-values */ print,"Parameter Estimates",, beta stdb t prob; print,y yhat resid; reset noprint; finish regress; /* Ä£¿é½áÊø */ x={1 1 1, 1 2 4, 1 3 9, 1 4 16, 1 5 25}; y={1,5,9,23,36}; reset print; run regress; /* ÔËÐÐÄ£¿é */
五．正交多项式回归 使用 ORPOL function 作二次正交多项式回归.
x1={1,2,3,4,5}; /* second column of X */ x=orpol(x1,2); /* generates orthogonal polynomials */ reset noprint; /* turns off auto printing */ run regress; /* run REGRESS */
Regression Results
SSE DFE MSE RSQUARE 6.4 2 3.2 0.9923518
Parameter Estimates
