*===== 7章：行列計算と数値積分 ;
*--- 7.1 行列の作成 ;
proc iml ;
X={1 2,
3 4} ;
y={1 2 3 4} ;
print X y ; * 行列の表示(横) ;
print X,y ; * 行列の表示(縦) ;
quit ;
proc iml ;
E = i(3) ; * 3×3の単位行列 ;
F = do(1, 5, 2) ; * (1, 3, 5) なるベクトル;
J = j(2, 3, 1) ; * 2×3の「1行列」;
K = {1 2, 3 4} ;
L = repeat(K, 1, 2) ; * 行列Kを横に2つ並べる ;
print E F J K L ;
quit ;
proc iml ;
call randseed(7777) ;
X = J(1, 10, .) ; * 乱数を格納する空箱を用意 ;
Y = J(1, 10, .) ; * 乱数を格納する空箱を用意 ;
* 超幾何分布 ;
call randgen(X,'HYPERGEOMETRIC', 40, 3, 20) ;
* テーブル分布 ;
p = {0.1, 0.2, 0.3, 0.4} ;
call randgen(Y, 'TABLE', p) ;
print X, Y ;
quit ;
proc iml ;
call randseed(777) ;
X = J(1, 5, .) ; * 乱数を格納する空箱を用意 ;
call randgen(X, "NORMAL", 0, 1) ; * 平均0, 分散1 ;
print X ;
quit ;
proc iml ;
MEAN = {0, 0} ; * 平均(0,0) ;
COV = {1 0.5, 0.5 1} ; * 分散(1,1), 共分散0.5 ;
call vnormal(X, MEAN, COV, 5, 777) ; * 5組生成 ;
print X ;
quit ;
proc iml ;
call randseed(777) ;
df = 20 ;
sigma = {1 0.5, 0.5 1} ;
X = randwishart(1, df, sigma) ;
Y = shape(X, 2, 2) ; * X を 2×2 行列に変換 ;
print X Y ;
quit ;
data RAW ;
do I=1 to 2 ;
X = I ; Y=I+2 ;
output ;
end ;
run ;
proc iml ;
use work.RAW ;
read all var {X Y} into Z ;
close work.RAW ;
print Z ;
quit ;
proc iml ;
X = {1 2 3, 4 5 6} ;
Y = {7 8 9 0 1 2} ;
C = {"A" "B" "C"} ;
create MYDATA1 from X[colname=C] ;
append from X ;
close X ;
create MYDATA2 var {X Y} ;
append ;
close X Y ;
quit ;
*--- 7.2 行列へのアクセス ;
proc iml ;
X = {1 2,
3 4} ;
Z1 = X[1,2] ; * 1行2列目の要素 ;
Z2 = X[1, ] ; * 1行目の行ベクトル ;
Z3 = X[ ,2] ; * 2列目の列ベクトル ;
Z4 = X[+, ] ; * 行方向（縦方向）の和を計算した行ベクトル ;
Z5 = X[ ,#] ; * 列方向（横方向）の積を計算した列ベクトル ;
Z6 = X[<>,] ; * 行方向（縦方向）の最大値を計算した行ベクトル ;
Z7 = X[><,] ; * 行方向（縦方向）の最小値を計算した行ベクトル ;
Z8 = X[:, ] ; * 行方向（縦方向）の平均を計算した行ベクトル ;
Z9 = X[##] ; * 全ての要素に対して平方和を計算した結果の値 ;
print X, Z1 Z2 Z3, Z4 Z5 Z6, Z7 Z8 Z9 ;
quit ;
proc iml ;
X = {1 2 3,
4 5 6,
7 8 9} ;
Z = X[1:2,2:3] ; * 行列Xの「1～2行目」「2～3列目」にアクセス ;
print X, Z ;
quit ;
proc iml ;
X = {1 2,
3 4} ;
Y = {5 6,
7 8} ;
Z1 = X//Y ; * 縦結合 ;
Z2 = X||Y ; * 横結合 ;
print X Y, Z1 Z2 ;
quit ;
proc iml ;
X = {1 2 2 3} ;
Y = {1 3 5 7} ;
Z1 = unique(X) ; * 重複削除 ;
Z2 = union(X, Y) ; * 和集合 ;
Z3 = xsect(X, Y) ; * 積集合 ;
Z4 = setdif(X, Y); * 差集合 ;
print X, Y, Z1, Z2, Z3, Z4 ;
quit ;
*--- 7.3 行列計算 ;
proc iml ;
X = {1 2,
3 4} ;
Y = {5 6,
7 8} ;
Z1 = X + Y ; * 和 ;
Z2 = X - Y ; * 差 ;
Z3 = X * Y ; * 積 ;
Z4 = X # Y ; * いわゆる行列の積ではない！ ;
Z5 = X**2 ; * 行列のべき乗 ;
print X Y, Z1 Z2 Z3, Z4 Z5 ;
quit ;
proc iml ;
X = {1 3,
2 4} ;
Z1 = det(X) ;
Z2 = diag(X) ;
print X, Z1 Z2 ;
quit ;
proc iml ;
X = {1 3,
2 4} ;
Z1 = eigval(X) ;
Z2 = eigvec(X) ;
print X, Z1 Z2 ;
quit ;
proc iml ;
X = {1 2,
2 4} ;
Z = inv(X) ;
print X, Z ;
quit ;
proc iml ;
X = {1 2,
2 4} ;
Z = ginv(X) ;
print X, Z ;
quit ;
proc iml ;
X = {1 2,
2 4} ;
Z = half(X) ;
*Z = root(X) ;
Z2 = t(Z)*Z ;
print X, Z Z2 ;
quit ;
proc iml ;
X = {1 2,
2 4} ;
call qr(Q, R, PIV, LINDEP, X) ;
RANK= ncol(X) - LINDEP ; * ランク ;
print X, Q R RANK ;
quit ;
proc iml ;
X = {1 3,
2 4} ;
call svd(U, Q, V, X);
print X, U Q V ;
quit ;
proc iml ;
X = {1 3,
2 4} ;
Y = {5,
6} ;
Z = solve(X,y) ;
print X Y, Z ;
quit ;
proc iml ;
X = {1 3,
2 4} ;
Z1 = t(X) ;
Z2 = trace(X) ;
Z3 = vecdiag(X) ;
print X, Z1 Z2 Z3 ;
quit ;
data MYDATA ;
input X Y @@ ;
cards ;
1 8
1 4
0 7
0 5
0 3
;
run ;
proc glm data=MYDATA ;
model Y = X / solution ;
run ;
proc iml ;
y = {8, 4, 7, 5, 3} ;
X = {1 1,
1 1,
1 0,
1 0,
1 0} ;
b = inv( t(X)*X )*t(X)*y ;
print y X, b ;
quit ;
proc iml ;
y = {8, 4, 7, 5, 3} ;
X = {1 1 0,
1 1 0,
1 0 1,
1 0 1,
1 0 1} ;
b = ginv( t(X)*X )*t(X)*y ;
print y X, b ;
quit ;
*--- 7.4 プログラミング ;
proc iml ;
X = 1 ;
if (X > 1) then Y = 1 ;
else Y = 2 ;
print X Y ;
quit ;
proc iml ;
X = 1 ;
Y = choose(X > 1, 1, 2) ;
print X Y ;
quit ;
proc iml ;
X = J(1, 10, 0) ;
print X ;
do i=1 to ncol(X) ;
X[i] = i ;
end ;
print X ;
quit ;
proc iml ;
X = J(1000, 1, .) ;
call streaminit(777) ;
do I=1 to 1000 ;
X[I] = rand("TABLE",1/6,1/6,1/6,1/6,1/6,1/6)
+rand("TABLE",1/6,1/6,1/6,1/6,1/6,1/6)
+rand("TABLE",1/6,1/6,1/6,1/6,1/6,1/6)
+rand("TABLE",1/6,1/6,1/6,1/6,1/6,1/6) ;
end ;
MEAN = X[:] ;
print MEAN ;
quit ;
*--- 7.5 関数の定義と数値積分 ;
proc iml ;
start f(x) ;
return( x+1 ) ;
finish ;
y = f(2) ;
print y ;
run ;
proc iml ;
start f(x) ;
fx = x + 1 ;
return(fx) ;
finish ;
call quad(RESULT, "f", 0 || 1 ) ;
print RESULT ;
quit ;
proc iml ;
start f(x) global(a) ;
fx = x + a ;
return(fx) ;
finish;
a = 1 ;
call quad(RESULT, "f", 0 || a );
print RESULT ;
quit ;
proc iml ;
start f1(x) global(y) ;
fx = 2*x + y ;
return(fx) ;
finish ;
start f2(z) global(y) ;
y = z ;
call quad(gy, "f1", 0 || 1 ) ;
return(gy) ;
finish ;
call quad(RESULT, "f2", 0 || 2 ) ;
print RESULT ;
quit ;
proc iml ;
start f1(x) global(y) ;
fx = 2*x + y ;
return(fx) ;
finish ;
start f2(z) global(y) ;
y = z ;
call quad(gy, "f1", 0 || y/2 ) ;
return(gy) ;
finish ;
call quad(RESULT, "f2", 0 || 2 ) ;
print RESULT ;
quit ;
ods graphics on ;
proc seqdesign plots=all ;
Pocock1: design method=poc nstages=2 alt=upper alpha=0.025 ;
run ;
ods graphics off ;
proc iml ;
start f1(u2) global(u1) ;
f = exp(-u1**2/2)*exp(-(u2-u1)**2/2)/(2*constant('pi')) ;
return(f) ;
finish ;
start f2(x) global(u1, c) ;
u1 = x ;
call quad(f, "f1", -10 || (c#sqrt(2)) ) ;
return(f) ;
finish ;
c = 2.17827 ;
call quad(f, "f2", -10 || c ) ;
alpha = 1 - f ;
print alpha ;
quit ;
*************************************************;
*--- 7章 演習問題 ;
* (1)-(3) ;
proc iml ;
X = {1 1 1,
2 2 2,
3 3 3} ;
Y = {0 1 2,
2 3 4,
4 5 6} ;
Z = X + Y ;
T = ginv( Z ) ; * 関数 inv() だと行列式が0であるため計算できない ;
print X, Y, Z, T ;
quit ;
* (4) ;
proc iml ;
start f(x) ;
fx = exp(-x**2/2)/sqrt(2*3.14) ;
return(fx) ;
finish ;
call quad(RESULT, "f", .M || .P ) ;
print RESULT ;
quit ;