SAS/IML

Лекция 1. Основы IML.

Звежинский Дмитрий, SAS Russia/CISdmitry.zvezhinsky@sas.com

1Замечания об ошибках и опечатках просьба направлять лектору.

http://tiny.cc/msu_sas

Что такое IML?• Это процедура языка SAS (PROC IML;), внутри которой можно создавать свои

программы на специальном языке.

• Синтаксис этого языка упрощает работу с матрицами (линейная алгебра, матричные операторы), а также содержит встроенные алгоритмы для решения разных задач:

– Численное решение линейных уравнений, систем уравнений

– Численное решение диф. уравнений и систем диф. уравнений

– Статистический анализ в духе «Большого SAS», т.е. генерирование отчетов с определенной структурой

– Манипуляции наборами данных SAS

– Написание собственных процедур/функций, которые, правда, работают только внутри IML

– Написание программ: циклы, операторы условного перехода, …

– Оптимизация

– Вэйвлет- и Фурье-преобразования

– Вызов процедур и функций из «Большого SAS»

– Выполнение процедур основного языка Base SAS и обмениваться с ним данными

– Выполнение процедур языка R (то есть в том числе позволяет интегрировать R и BaseSAS)

2

Простейшие матричные операции• Основной объект – двумерная матрица.

• Тип матрицы – числовой или символьный (смешивать нельзя).

proc iml;

a={1 2,3 4};

print (a);

quit;

15 proc iml;

NOTE: IML Ready

16 a={1 'a',3 'b'};

ERROR: Mixing character with numeric in matrix literal at line=16 col=8.

17 print (a);

ERROR: Matrix a has not been set to a value.

• Доступ к элементам матрицы – по индексам:proc iml;

a={1 2 3,4 5 6};

print a ,, (a[1,2]) ,, (a[{1 2},2]);

quit;

Индексы можно задавать с помощью матрицы

3

Простейшие матричные операции• Особенность – всё кроме цифр считается символьным

значением:proc iml;

a={1 2,3 4};

b={1 1+1,3 4};

c={1 (1+1),3 4};

d1=1+1;

d2={1 d1,3 4};

print a,,b,,c;

quit;

17 b={1 1+1,3 4};

_

22

200

ERROR: Unequal number of elements per row in literal at line=17 col=16.

4

Матрицы• Можно посмотреть тип матрицы и кол-во строк/столбцов:proc iml;

a={1 2,3 4};

print (a[1,1]);

show names;

a=j(2,2," ");

a[1,1]="PRIVET!";

print (a[1,1]);

show names;

quit;

• Матрицы можно переопределять. Размер элемента символьной матрицы можно поменять, только переопределив её.

• Матрицы можно сохранить в каталоге - хранилище (WORK.IMLSTOR), и пользоваться ими на другом шаге IML:

proc iml;

a={1 2,3 4};

store a;

quit;

proc iml;

load a;

print (a[1,1]);

quit;

proc catalog catalog=work.imlstor;

contents;

run;

5

Простейшие матричные операции• Есть матричные и поэлементные операции (проверяйте размерность!)

• Операции сравнения и минимума/максимума:proc iml;

a={1 2,3 4};

b={0 0,5 6};

print a b ,'a<>b',(a<>b); *set maximum;

b={2,3};

print a b, 'a<b',(a<b);

quit;

= Т( )

6

Subscript reduction operators

• Позволяют преобразовывать матрицы по заданному правилу по указанному измерению:

proc iml;

a={1 2 3, 4 5 6,7 8 9};

print (a[+, ]),,

(a[+,<>]),,

(a[+,<:>]);

quit;

Пример 1: суммируем по строчкам, кол-во

Столбцов остаётся как в исходной матрице.

Пример 2: сначала суммируем по строчкам, потом

ищем максимум в получившейся матрице.

Аналогичный результат дает последовательное

применение reduction operators: (a[+,][,<>])

Пример 3: суммируем по строчкам, и в полученной матрице ищем индекс (номер столбика), где лежит максимум.

Если reduction operators есть сразу в двух измерениях, сначала выполняется первый (для строк), потом второй (для столбиков).

операторы сведения (матриц) по индексам

7

Простейшие матричные преобразования• inv (A) – обратная матрица

• trace (A) – след матрицы (сумма диаг. эл-тов A*A)

• abs(A) – модуль каждого элемента

• sum(A,B) – сумма всех элементов из всех матриц

• nrow(A), ncol(B) – кол-во строк (столбцов)

• j (n, m, “privet!”) – матрица n строк, m столбцов, состоящая из текстовой константы “privet!”

• diag(A) – возвращает матрицу того же размера, все элементы кроме диагональных заполняются нулями

• i (n) – единичная матрица размером (n x n)

• loc (A) – возвращает индексы ненулевых элементов, в частностиможно найти эл-ты, удовлетворяющие условию.

8

МодулиВстроенные функции:

• Support.sas.com/documentation

• SAS/IML 9.3 User's Guide

Пользовательские (называются модули):

• Можно их определять прямо в proc iml :

START <name> <( arguments )> <GLOBAL( arguments )> ;

…

FINISH <name> ;

• Пользоваться можно тоже только внутри IML.

• Есть в виде подпрограмм (subroutine, запускаются оператором CALL), либо в виде функций (function), которые могут вернуть какое-то значение. Отличаются только наличием оператора Return внутри модуля (если он есть – это функция).

9

Модули• Область видимости для переменных: зависит от того, есть ли в определённом

модуле аргументы

• Из модуля-процедуры без аргументов можно ссылаться на глобальные переменные, которые определятся позднее (из модуля-функции нельзя):

proc iml;

start test ;

print a;

finish;

a=5;

call test;

quit;10

Есть (it depends): Нет (все перем. глобальные):

Модули• Вложенные модули преобразуются так, чтобы у них был один (и тот же)

уровень вложенности.

• Выражения в качестве аргументов – вычисляются, но хранятся как локальная копия аргумента:

proc iml;

start test(in,out);

out=in+1;

print in out;

finish;

inx={1 2 3};

outx={1 2 3};

out=0;

call test(inx[1],outx[1]);

print (outx[1]);

call test(inx[1],out);

print (out);

quit; 11

Значение выражения (которое нужно как-то вычислить) видно внутри модуля, но оно хранится как локальная копия!

Во втором вызове out – это не выражение (его не нужно вычислять), а имя матрицы. Через него можно возвратить значение в глобальную область видимости.

Глобальная или локальная переменная?

12

1. Переменная инициализируется вне модуля – глобальная2. Переменная инициализируется внутри модуля, но у модуля отсутствуют

аргументы – глобальнаяstart test;

output=input+1;

finish;

input=5;

output=0;

call test;

print output;

3. Переменная передаётся как аргумент в модуль – внутри модуля с ней можно работать как с глобальнойstart test(input,output);

output=input+1;

finish;

in=5;

out=0;

call test(in,out);

4. С помощью опции Global указываем переменную – внутри модуля она будет считаться глобальной.5. Если в качестве аргумента передаётся выражение, а не матрица, то переменная, которая в него входит, не считается глобальной переменной.

out

6

output

6

Резюме

Глобальная или локальная переменная?proc iml;

start sum1(x);

return (x[1]+x[2]);

finish;

x={1 2 3 4 5};

do i=1 to 4;

print(sum1(x[i]||x[i+1]));

end;

quit;

proc iml;

x={1 2 3 4 5};

start sum1;

c=x[i]+x[i+1];

finish;

c=0;

do i=1 to 4;

call sum1;

print (c);

end;

quit;13

х – глобальнаяа внутри модуля?(локальная)

В качестве аргумента передаём глоб. переменную

х – глобальнаяа внутри модуля?(как и «i», глобальная)

Глобальная или локальная переменная?proc iml;

start sum1(m) global (c,x);

c=x[m]+x[m+1];

finish;

x={1 2 3 4 5};

c=0;

do i=1 to 4;

call sum1(i);

print (c);

end;

quit;

14

х – глобальная, m-локальная

Создаём вектор-строку х

Вызываем модуль.

Хранение модулей• В IML есть механизм для хранения модулей и матриц

(слайд 5) для их последующего использованияproc iml;

start test;

print "this is a test";

finish;

STORE MODULE= test;

quit;

proc iml;

LOAD MODULE= test;

call test;

quit;

• Можно просматривать название модулей: show storage;

15

3568 STORE MODULE= test;NOTE: Opening storage library WORK.IMLSTOR3569 quit;NOTE: Exiting IML.NOTE: Storage library WORK.IMLSTOR closed.NOTE: PROCEDURE IML used (Total process time):

real time 0.00 secondscpu time 0.00 seconds

3570 proc iml;NOTE: IML Ready3571 LOAD MODULE= test;NOTE: Opening storage library WORK.IMLSTOR

Элементы языка IML• Циклы:

16

скобочки в условии обязательны!

скобочки в условии обязательны!

Элементы языка IML• Оператор условного перехода:

proc iml;

if . then do;

print ". is true";

end;

if 0 then do;

print "0 is true";

end;

quit;

• Условие должно возвращать число. Истина – это любое значение кроме нуля и пропущенного значения.

• Если в условии один оператор, то группу do; .. end; можно не добавлять.

• есть дополняющий оператор “else … ;”17

Программа ничего не печатает.

Попробуйте добавить оператор отрицания «^» перед условием.

Элементы языка IML• Оператор PRINT – вывод в отчет результатов и отладочной

информации. Общий вид вот такой:PRINT <matrices> (expression) "message" pointers <[options]> ;

• Если нужно вычислить выражение (вывести элементы матрицы, арифметика, функции) – указывайте его в скобочках!

• pointers – управляет разбивкой по строкам ( , ) и по страницам ( / )

• Можно подписывать строки и столбцы при выводе на экран:proc iml;

x = {45.125 50.500,

75.375 90.825};

r = {"Div A" "Div B"};

c = {"Amount" "Net Pay"};

print x[rowname=r colname=c format=12.2];

quit;

18

Работа с набором данных SAS

• USE – открыть набор данных для чтения• EDIT – открыть набор данных для чтения и записи• CREATE – создать набор данных, при этом

открытый/созданный набор данных становится активным

• SETIN – выбрать текущий набор для ввода• SETOUT – выбрать текущий набор для вывода.

Используются, например, если нужно только открыть набор данных и сделать активным, ничего не считывая и не записывая:

• SHOW – отобразить информацию о наборах данных

19

Интеграция IML и SASЧтение наборов данных:

1) Открыть набор данных для чтения:proc iml;

use test;

2) Прочитать нужные столбцы в матрицы:read all var {a} into char_test;

read all var {b} into num_test;

print char_test num_test;

quit;

20

Интеграция IML и SAS• Мы помним, что в матрице можно хранить данные только

одного типа. Размер для хранения данных берётся из дескриптора набора данных:

data test;

length a $ 3 text $ 200;

input a $ b text $ ;

datalines;

aaa 10 blahblahblah

bbb 20 blahblahblah

ccc 30 blahblahblah

ddd 40 blahblahblah

;

run;

proc iml;

use test;

show datasets;

read all var {a} into char_test;

read all var {b} into num_test;

read all var {text} into text_test;

print (char_test || text_test);

show names;

quit;

21

SYMBOL ROWS COLS TYPE SIZE

------ ------ ------ ---- ------

char_test 4 1 char 3

num_test 4 1 num 8

text_test 4 1 char 200

Интеграция IML и SASЗапись наборов данных

• 1) Создать набор данных с помощью оператора Create,

• 2) заполнить набор данных из матрицы – оператор Append:proc iml;

x=repeat(do(1,3,1),3,1);

y=repeat(T(do(1,3,1)),1,3);

print x y;

create test1 from x;

append from x;

create test2 from y;

append from y;

quit;

3849 quit;

NOTE: Exiting IML.

NOTE: The data set WORK.TEST1 has 3 observations and 3

variables.

NOTE: The data set WORK.TEST2 has 3 observations and 3

variables.

22

Интеграция IML и SASОператор Use – открыть набор данных и сделать его текущим для ввода/вывода (то есть тем набором, который будет использован операторами Read, Edit, Append, …)

Setin/Setout – выбрать текущий набор для ввода/выводаcreate test1 from x;

create test2 from y;

setout test1;

append from x;

setout test2;

append from y;

Close .. ; – закрыть набор данных, открытый операторами Use, Create, Edit

Show datasets; - статус наборов данных, использованных в текущем сеансе IML.

23

Встроенные подпрограммы• В программах SAS IML вы можете использовать

большинство функций, доступных в Base SAS.

• Перечисленные ниже функции Base SAS либо не доступны из SAS/IML, либо их поведение отличается от функций Base SAS с тем же названием:

• http://support.sas.com/documentation/cdl/en/imlug/66845/HTML/default/viewer.htm#imlug_langref_sect503.htm

• В IML есть и свои подпрограммы, самые интересные мы обсудим далее.

24

Преобразования типов

• Символы в число: NUMc = {"1" "2" "3"};

reset print; /*display values and type of matrices*/

m = num(c);

• Числа в символы: CHAR( matrix <, w> <, d> ) ;

a = {-1.1 0 3.1415 4};

reset print; /*display values and type of matrices*/

m = char(a, 4, 1);

25

Вызов процедур SAS из IMLSUBMIT <parameters> </ options> ;

language statements

ENDSUBMIT;

• Пример :proc iml;

x=do(-10,10,0.5);

eps=j(1,ncol(x),1);

call randgen(eps, "normal", 0,2);

y=3#x+5+eps;

data_iml=(T(x//y));

create test1 from data_iml;

append from data_iml;

close test1;

SUBMIT;

proc reg data=test1;

model col2=col1;

quit;

ENDSUBMIT;

quit;

26

!!! Наличие пробела между “T” и “;” приводит к зависанию IML

Любопытно, что внутрь SUBMIT можно передавать переменные из IML:

proc iml;

VarName = "Sex";

submit VarName;

proc freq data=Sashelp.Class;

table &VarName / out=OutFreq;

run;

endsubmit;

quit;

Работа с памятью• Для размещения матриц используется выделяемая в

памяти область (Workspace)

• Что будет, если мы всю её займём?

ERROR: (execution) Unable to allocate sufficient memory. At least 2147483647 more bytes required.

• Память выделяется частями (extent). Сначала SAS пытается выделить всё больше и больше extent’ов. Память сжимается (избавляемся от промежуточных временных переменных). Когда мы достигаем ограничения, выдается ошибка.

• Освобождение памяти:FREE A B C; *удаляет из памяти матрицы A,B,C;

FREE /; *удаляются все матрицы;

FREE / A B; *удаляются все матрицы, кроме A и B;

27

Работа с памятью• Изменение настроек собственно PROC IML для работы с

памятью делается опциями worksize= и symsize= (после которых указывается лимит памяти в килобайтах):

PROC IML <SYMSIZE=n1> <WORKSIZE=n2> ;

• Worksize= память для хранения содержимого матриц

• Symsize= память для хранения таблицы с именами матриц

• Команда для вывода отладочной информации о кол-ве extents, памяти, кол-ве сжатий:

show space;

• Команда для вывода информации о созданных и хранящихся в памяти переменных:

show names;

28

Матричные преобразования• Разложение SVD:proc iml;

A={1 2 3,

4 5 6,

8 9 10,

11 12 13};

CALL SVD( u, q, v, A ) ;

/*check decomposition*/

print A ,,(u*diag(q)*T(v)) ;

quit;

• Разложение Холецкого:Матрица А – симметричная и положительно определённая.proc iml;

A={10 0 0,

0 1 2,

0 2 9};

u=ROOT(A) ;

/*check decomposition*/

print A ,,(T(u)*u) ;

quit;

29

U, V = унитарные матрицы q – столбец, содержащий диагональные элементы матрицы с сингулярными собств. числами

where is upper triangular.

Матричные преобразования• QR- разложение:

CALL QR( q, r, piv, lindep, a <, ord> <, b> ) ;

proc iml;

A={1 2 3,

4 5 6,

7 8 9};

call qr(q,r,piv,lindep,A) ;

/*check decomposition*/

print q ,,(q*T(q)),r,,(q*r) ;

quit;

30

Матричные преобразования• Собственные векторы, собственные значения.

• Матрица А имеет только действительные элементы. Полученные с.в., с.з. имеют или одну (действит.) или две (действит. + мнимую) составляющих.

CALL EIGEN( evals, evecs, A ) ;

Ax = λx

CALL GENEIG( eval-M, evecs-E, A, B ) ;

EIGVEC( A ) ;

EIGVAL( A ) ;

proc iml;

M={0 -1,1 0};

call eigen(evals, evecs, M);

print evals ,, evecs;

quit;

31

1/√2

-i/√2

Матричные преобразования• Если матрица А не симметричная, собственные значения получаются

комплексные, evals содержит 2 колонки. Первая – действительная часть, вторая – мнимая.

• Собственные векторы кодируются в матрице (evecs) по следующему правилу. Если строчки i, i+1 в матрице собственных значений (evals)содержат комплексно-сопряженные соб.значения, то в матрице соб.векторов (evecs) содержится действительная и мнимая части для координат комплексно сопряженных соб. векторов

32

1/√2

-i/√2= +i

1/√2

i/√2

1/√2

i/√2= -i

1 с.з.2 с.з.

Две компоненты к.с. с.в.

Re,Im

Re,Im

1/√2

-i/√2

К.С.

Матричные преобразования• Если матрица А симметричная, и собственные значения действительные:

proc iml;

M1={1 2,2 1};

call eigen(evals, evecs, M1);

print evals ,, evecs;

quit;

1 и 2 с.з.

1 и 2 соб. вектор

(Дискретное) Фурье-преобразование

• Оно есть!• Для примера посмотрим на фильтрацию суммы синусов:%let fd=10; *discretization frequency;

proc iml;

t=do(-25.5,25.6,1/&fd);

pi=3.1415926;

y=sin(2*pi*0.2#t)+0.5*sin(2*pi*0.5#t);

y1=y;

do i=1 to ncol(t);

y1[i]=y1[i]+RAND('NORMAL',0,0.5);

end;

PLT1=T(t)||T(y)||T(y1);

create plot1 from plt1;

append from plt1;

quit;

goptions reset;

goptions hsize=500pt vsize=400pt;

SYMBOL1 COLOR=red INTERPOL=JOIN value=none;

SYMBOL2 COLOR=blue INTERPOL=JOIN value=none;

proc gplot data=plot1;

plot col2*col1 col3*col1/ overlay;

run;

quit;

34

Фурье-преобразованиеproc iml;

use plot1;

read all into plt1;

NFFT=nrow(plt1); *points in signal;

fft_y=fft(plt1[,2])/(NFFT /2 );

*set norm for fourier transform;

harm=sqrt(fft_y[,1]##2+fft_y[,2]##2);

*amplitude of harmonic;

fd=&fd ; *discretization frequency;

f=do(0,1,2./NFFT);

*frequencies for fourier-transformed

signal;

Теперь можно построить:

• Гармоники исходного сигнала

• и зашумленного

35

Фурье-фильтр• Пара строчек:

fil_fft_y=fft_y;

fil_fft_y[loc(harm<0.1),]=0;

• И обратное Фурье-преобразование:

/*and inverse fourier transform*/

fil_y=ifft(fil_fft_y)/2;

• Теперь сбрасываем всё что нужно в наборы данных, и строим графики:

fourier=T(fd/2#f)||harm||harm_1;

create plot2 from fourier;

append from fourier;

fil_y_tab= plt1[,1]||fil_y;

create plot3 from fil_y_tab;

append from fil_y_tab;

quit;

36

Фурье-фильтр• Объединим наборы данных для построения исходного, зашумленного и

отфильтрованного сигнала:data merged;

set plot3(rename=(col1=t1 col2=yfilt));

set plot1(rename=(col1=t col2=y col3=noisy));

run;

• И построим график:SYMBOL1 COLOR=green INTERPOL=none value=dot;

SYMBOL2 COLOR=red INTERPOL=JOIN value=none;

SYMBOL3 COLOR=blue INTERPOL=JOIN value=none;

proc gplot data=merged;

plot noisy*t y*t yfilt*t / overlay;

run;

quit;

37

Полезные ссылки

• http://support.sas.com/documentation/cdl/en/imlug/67502/HTML/default/titlepage.htm

• http://blogs.sas.com/content/iml/files/2011/10/IMLTipSheet.pdf

38