SAS IML. Использование языка SAS IML. Основы IML (Лекции 2013), страница 2

At least2147483647 more bytes required.• Память выделяется частями (extent). Сначала SAS пытаетсявыделить всё больше и больше extent’ов. Памятьсжимается (избавляемся от промежуточных временныхпеременных). Когда мы достигаем ограничения, выдаетсяошибка.• Освобождение памяти:FREE A B C; *удаляет из памяти матрицы A,B,C;FREE /; *удаляются все матрицы;FREE / A B; *удаляются все матрицы, кроме A и B;25Работа с памятью• Изменение настроек собственно PROC IML для работы спамятью делается опциями worksize= и symsize= (послекоторых указывается лимит памяти в килобайтах):PROC IML <SYMSIZE=n1> <WORKSIZE=n2> ;• Worksize= память для хранения содержимого матриц• Symsize= память для хранения таблицы с именами матриц• Команда для вывода отладочной информации о кол-веextents, памяти, кол-ве сжатий:show space;• Команда для вывода информации о созданных ихранящихся в памяти переменных:show names;26Матричные преобразования• Разложение 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;27Матричные преобразования• 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;28Матричные преобразования• Собственные векторы, собственные значения.• Матрица А имеет только действительные элементы.Полученные с.в., с.з.

имеют или одну (действит.) или две (действит. +мнимую) составляющих.CALL EIGEN( evals, evecs, A ) ;Ax = λxCALL 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;29Матричные преобразования• Если матрица А не симметричная, и собственные значения получаютсякомплексные, evals содержит 2 колонки. Первая – действительная часть,вторая – мнимая.• Если строчки i, i+1 в матрице собственных значений (evals) содержаткомплексно-сопряженные соб.значения, то в матрице соб.векторов(evecs) содержится действительная и мнимая части для координат одногоиз соб. векторов.Иллюстрация результата, полученного на прошлом слайде:1/√2= +i-i/√21/√2Re,Im-i/√21 с.з.2 с.з.1/√2i/√2= -i1/√21 с.в.2 с.в.i/√2Re,Im30Задание• Напишите свой модуль, который вычисляет определитель матрицы(матрица передаётся как параметр).

Нельзя пользоваться функциейdet. Результат задания – программа *.sas• Вычислять можно любым методом, какой вам нравится.• http://stackoverflow.com/questions/2435133/what-is-the-best-matrixdeterminant-algorithm• Оценка от 0 до 3:• 0 – не сделал(а) или программа считает det неправильно• 1 – определитель вычисляется для матрицы фиксированного размера• 2 – п. «1» + сообщения об ошибках (неквадратная матрица, ???)• 3 – п.

«2» + определитель вычисляется для матрицы произвольногоразмера• Сохраните его в хранилище для последующего вызова• Загрузите модуль из хранилища и вызовите его в другой процедуреIML, написав пример(ы) его использованияБольше по этой лекции заданий нет.31Использование языка SAS/IML• Supplementary: Решение разных задач32Статистическое описание данных•Создание отчета с рассчитанными статистиками, в том числе в разбивке поклассифицирующим переменным.

На вход подается набор данных, открытыйоператором USE.SUMMARY <CLASS operand> <VAR operand> <WEIGHT operand> <STAT operand><OPT operand> <WHERE(expression)> ;proc iml;use sashelp.class;SUMMARY CLASS {sex} VAR {age} STAT {mean min max};quit;Analysis variable : AgeSex NobsMEANMINMAX--------------------------------------------F913.2222211.0000015.00000M1013.4000011.0000016.00000All1913.3157911.0000016.00000---------------------------------------------• Расчет матрицы с частотами:CALL TABULATE( levels, freq, x <, method> ) ;• Распечатка открытого набора данных:LIST <range> <VAR operand> <WHERE(expression)> ;33Графика*Чтобы построить график массива , созданного в IML, предлагается сбросить его внабор данных, а затем подать его на вход процедуры gplot:*let's plot sine and cosine;proc iml;do x1=-10 to 10 by 0.1;x=x||x1;y=y||sin(x1);y2=y2||cos(x1);end;A=T(x//y//y2);create work.temp1from A[colname={"x","sine","cosine"}];append from A;quit;proc gplot data=temp1;*goptions vsize=3 hsize=3;SYMBOL1 COLOR=red INTERPOL=JOIN value=x;SYMBOL2 COLOR=green INTERPOL=JOIN value=none;*legend;legend1 label=("Fancy plot") value=("sine plot" "cosine plot");axis1 label=('y=f(x)');plot sine*x cosine*x / overlay legend=legend1 vaxis=axis1;run;quit;* В IML есть собственные графические процедуры, но мы их изучать не будем34ГрафикаПостроение графиков можно сделать меньшим кол-вом кода IML:proc iml;/** define functions **/start Func1(x);return( sin(x));определилиfinish;функцииstart Func2(x);return( cos(x));finish;вектор-строка Хx = do(-5,5,0.1); *make matrix row;y1=Func1(x);без циклов посчиталиy2=Func2(x);вектор-строку УA=T(x//y1//y2);create work.temp1 from A[colname={"x","sine","cosine"}];append from A;quit;35Решение ОДУ••••В IML присутствует только одна процедура для этого: ODEThe ODE subroutine is an adaptive, variable-order, variable-step-size, stiff integrator based onimplicit backward-difference methods.

…Как и любой численный метод, ODE должен гарантировать нам решение с точностью доуказанного эпсилон (по умолчанию, 10^(-4))… Проверим!Давайте начнем его проверять с простого ОДУ:y'(t)=-t*y(t)•Начальные условия:•Аналитическое решение:y(t=0)=136Решение ОДУ/*plot ode solve*/Модуль, где (система) ОДУ задаетсяproc iml;виде вектор-столбца dy/dt:start fun(t,y);v = -t*y; *the equation is y'(t)=-t*y(t);return(v);finish;/* Call ODE */c= 1; *vector of initial values;Начальное условие для y(t(1))t= do(0,10,0.1); *start time, end time; вектор t, для которых мы получим решениеh= {1E-12 1 1E-5}; * the minimum allowable step size, hmin;* the maximum allowable step size, hmax;* and the initial step size to start the integration process,hinit;call ode(r1, "FUN", c, t, h);…магия…TAB=T(t//(c||r1));решение не содержит начальное условие, мы его добавимcreate work.temp1 from TAB[colname={"t","solution"}];append from TAB;сбросим решение в наборquit;/*plot*/данных и построим его сproc gplot data=temp1;помощью gplot.goptions reset;SYMBOL1 COLOR=green INTERPOL=JOIN value=none;plot solution*t;run;quit;в37Решение ОДУ• Проверим и построим ошибку (EPS=1E-4)/*-=-check error-=-*/data errors;set temp1;analyt=exp(-t*t/2);err=analyt-solution;log10err=log(abs(err))/log(10);run;/*plot*/proc gplot data=errors;goptions reset;SYMBOL1 COLOR=greenINTERPOL=JOIN value=none;plot log10err*t / vref=-4;run;quit;По оси У отложен десятичн.логарифм ошибки – мывидим, что она «подконтролем» численногометода!38Решение системы ОДУили ОДУ высших порядков, что то же самое• Пример:y(t)''= exp(-y(t)' /10) * sin(2*t)Нач.

условия: y(0)=1, y'(0)=-1• Перепишем ур-е со второй производной как систему двух уравнений с первойпроизводной:y(t)’ =z(t)z(t)’ =exp(-z/10)*sin(2*t) = y(t)’’• Как будет выглядеть модуль для функции ODE?start ode2ord(t,y);y1=y[2];y2=exp(-y1/10)*sin(2*t);y=y1//y2;return(y);finish;•Вызовем функцию ODE:t-независ. переменная,y – вектор-столбец:[ y (t) ][ y’(t) ]на выходе – тоже вектор-столбец:[ y’(t) = f1() ][ y’’(t) = f2() ]c= {1, -1}; *vector of initial values;t= do(0,2,0.1); *start time, end time;h= {1E-12 1 1E-5};порядокcall ode(r1, "ODE2ORD", c, t, h);TAB=T(t//(c||r1));create work.temp1 from TAB[colname={"t","sol1","sol2"}];append from TAB;решений39Решение системы ОДУgoptions reset;goptions hsize=500pt vsize=250pt;SYMBOL1 COLOR=red INTERPOL=JOIN value=none;SYMBOL2 COLOR=blue INTERPOL=JOIN value=none;Legend1 value=('function y(t)');Legend2 value=('derivative y''(t)');proc gplot data=temp1;plot sol1*t / legend=legend1;plot sol2*t=2 / legend=legend2 ;run;quit;В этом примере мы непроверяем ошибку, ноповерим, что она тоже подконтролем (функция ипроизводная не содержаткаких-то быстрых изменений,где адаптивный алгоритммог бы “наврать”).40Решение жесткого (stiff) ОДУ• Уравнение, которое дает резкое возрастание погрешности принедостаточно малом шаге – хорошая проверка численных методов!• Пример с аналитическим решением *:dy/dt = y2- y3y(0) = δ0 < t < 2/δ• Модуль, задающий ОДУ + программа:proc iml;Проверяем, как работаетstart fun(t,y);v = y##2-y##3;адаптивный контрольreturn(v);шага.

Выводимfinish;последовательные шаги в/* Call ODE */набор данных steps:c= 1E-4; *vector of initial values;t= 0 || 2/c; *start time, end time;h= {1E-15 1 1E-5};call ode(r1, "FUN", c, t ,h) data="work.steps";*The data set is used to save the successful independent and dependentvariables of the integration at each step. The keyword for this option isDATA;quit;* http://www.mathworks.com/company/newsletters/articles/stiff-differential-equations.html41Жесткое ОДУ• Построим график решения у(t):goptions reset vsize=250pt hsize=500pt;SYMBOL1 COLOR=green INTERPOL=none value=x;proc gplot data=steps;plot y1*time;run;OOOPS!http://www.mathworks.com/cmsimages/64812_wl_may03_ode23_large.gif42Решение жесткого ОДУ• Ошибка явно вышла из-под контроля числ.