close

Вход

Забыли?

вход по аккаунту

Соотечественников ждут здесь, где начинается Россия!;pdf

код для вставкиСкачать
Решение ОДУ
•
•
•
•
В IML присутствует только одна процедура для этого: ODE
The ODE subroutine is an adaptive, variable-order, variable-step-size, stiff integrator based on
implicit backward-difference methods. …
Как и любой численный метод, ODE должен гарантировать нам решение с точностью до
указанного эпсилон (по умолчанию, 10^(-4))… Проверим!
Давайте начнем его проверять с простого ОДУ:
y'(t)=-t*y(t)
•
Начальные условия:
•
Аналитическое решение:
y(t=0)=1
1
Решение ОДУ
/*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;
в
2
Решение ОДУ
• Проверим и построим ошибку (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=green
INTERPOL=JOIN value=none;
plot log10err*t / vref=-4;
run;
quit;
По оси У отложен десятичн.
логарифм ошибки – мы
видим, что она «под
контролем» численного
метода!
3
Решение системы ОДУ
или ОДУ высших порядков, что то же самое
• Пример:
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;
решений
4
Решение системы ОДУ
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;
В этом примере мы не
проверяем ошибку, но
поверим, что она тоже под
контролем (функция и
производная не содержат
каких-то быстрых изменений,
где адаптивный алгоритм
мог бы “наврать”).
5
Оптимизация
• В SAS есть отдельная библиотека процедур для
решения различных задач оптимизации (SAS/OR)
• Посмотрим пример, в котором используются
разные подходы для решения одной и той же
задачи.
• Подберём коэффициенты логистической регрессии
с помощью метода максимального
(логарифмического) правдоподобия. Реализуем его
разными способами. Ответы должны быть
одинаковыми.
• Данные взяты тут:
http://www.jeremymiles.co.uk/regressionbook/data/index.html
6
Оптимизация
SAS/STAT, proc logistic.
proc logistic data=scores /* PLOTS = ALL*/ OUTEST=coef_log;
model pass(EVENT='1')=score exp;
run;
В контексте оптимизации это самая важная часть отчета. Она говорит нам, что численный
метод, который искал максимум функции логарифмического правдоподобия, сошелся, вы
можете доверять результатам анализа, в том числе коэффициентам регрессии. Может ли
быть по другому?
Конечно.
7
Оптимизация
• Quasi-complete separation (перевод этого понятия на
русский язык лектор не знает). При данной обучающей
выборке и модели данные идеально разделяются.
Оказывается, что для логистической регрессии это плохо
т.к. ведёт к бесконечно большим коэффициентам.
• Для нас это означает, что полученным оценкам
параметров мы доверять не можем.
• Оценки параметров для нашей модели (с оптимизацией
всё хорошо)
8
Оптимизация в IML
То же самое в IML:
use scores var _all_;
read all into scores_matrix;
Входные данные берём из
набора данных
Целевая функция для логистической регрессии – ф-ция логарифмического
правдоподобия (l).
start obj_f(p0)
global (scores_matrix);
X1=T(scores_matrix[,1]);
p(x)=1/(1 + exp(−  +  ))
X2=T(scores_matrix[,2]);
Y=T(scores_matrix[,3]);
unit_row=j(1,ncol(Y),1);
p1=1/(exp(-1*(p0[1]#unit_row+p0[2]#X1+p0[3]#X2))
+unit_row);
9
Оптимизация в IML
…
m_l=(p1##Y)#((unit_row-p1)##(unit_row-Y));
l_m_l=log(m_l);
result=-sum(l_m_l);
return(result);
finish obj_f;
10
Оптимизация в IML
Задаём начальную точку:
beta0=J(1,3,1); /* Declares the vector which will be used */
/* as initialization for the optimization routine. */
optn={0 4}; /* Optimization options: 0 to minimize, */
/* 4 to control what will be printed. */
call nlpnra(rc,beta,'obj_f',beta0); /* волшебство!*/
print beta;
quit;
11
Оптимизация SAS/OR
В SAS U её нет, и я с ней никогда не работал, поэтому просто посмотрим:…
proc optmodel;
var x {1..3} init 1;
/* declare variables for coefficients*/
set points;
number pass{points};
number score{points};
number exp{points};
read data scores into points=[n] pass[n]=pass score[n]=score
exp[n]=exp;
/* objective function */
max f = sum{n in points}
if (pass[n]=0) then
log(1-1/(1+exp(-1*(x[1]+x[2]*score[n]+x[3]*exp[n]))))
else log(1/(1+exp(-1*(x[1]+x[2]*score[n]+x[3]*exp[n]))));
/* now run the solver */
solve;
print x;
quit;
12
Оптимизация
• SAS/OR
13
Анализ главных компонент
• В SAS можно делать по разному. А мы будем в IML!
• Цель – нахождение линейных комбинаций переменных,
которые максимально объясняют дисперсию (разброс)
точек.
• Основной технический момент – диагонализация матрицы
ковариаций.
• Тот же анализ (с теми же результатами) можно делать
другим методом - через разложение SVD.
• Одно из популярных приложений – уменьшение
размерности данных (многие статистические модели сами
по себе не умеют бороться с этой проблемой)
14
Анализ главных компонент
Перед выполнением этого анализа рекомендуют приводить переменные к
одному масштабу и центрировать (стандартизовать). Здесь нужно думать, как это
делать в конкретном случае. Мы воспользуемся т.н. Z-scoring, предполагая
равное влияние всех факторов на разброс данных, а также распределение по
каждой переменной близким к нормальному.
x_std=std(x);
x_mean=mean(x);
do i=1 to 3;
X[,i]=(X[,i]-x_mean[i])/x_std[i];
end;
cov_m=cov(x);
Вычислим соб. векторы и соб. значения матрицы ковариаций
call eigen(eVal, eVec, Cov_m);
Обратите внимание на такую особенность этой call-routine:
/* If A is symmetric, the eigenvalues are arranged in descending order. */
Вычислим накопленную пропорцию объясненного разброса (данных) * :
VarExplained = cusum(eval)/sum(eval);
* http://sites.stat.psu.edu/~ajw13/stat505/fa06/16_princomp/04_princomp_coeffic.html
15
Анализ главных компонент
Первые две главные компоненты объясняют 90% разброса данных.
Спроецируем данные на найденные «главные» компоненты. А про третью
можно забыть.
NumPC = 2; /* keep 2 principal components */
scores = x*eVec[,1:NumPC];
/* это и есть проекции, scores */
16
Анализ главных компонент
• В SAS есть процедура Princomp, которая делает тот же
анализ (см. пример). Получаются те же результаты.
PROC PRINCOMP DATA = WORK.TEMP;
VAR COL1 COL2 COL3;
RUN;
17
Факторный анализ
• Цель: выделение «скрытых факторов» - некоторых ненаблюдаемых
характеристик описываемых данных. Соответственно задаются
вопросы – как эти факторы описываются исходными переменными,
насколько хорошо эти новые факторы объясняют данные.
• Одна из техник, которую используют для Ф.А. – метод главных
компонент. Некоторые утверждают, что Ф.А. это, скорее, техника
моделирования, а МГК – описательный метод (см.
http://en.wikipedia.org/wiki/Factor_analysis#Exploratory_factor_analysis
_versus_principal_components_analysis )
• Если упростить, то основная мысль этого анализа – это представить
исходные данные (матрицу) как произведение двух матриц:
10 признаков
100 человек
2 факторных
нагрузки
=
X
Характеристики 100
человек в терминах 2-х
выделенных факторов
*
L
F
18
Факторный анализ
• Оказывается, что матрица L связана с матрицей
ковариаций Х отношением:
S = COV( X ) = L LT + Ψ
• Если «забыть» про Ψ, то это равенство похоже на
спектральное разложение.
S = M Λ M-1 ; где у Λ по диагонали стоят с.з., а у М
по столбцам находятся соб. векторы
• Если ещё вспомнить, что матрица S положительно
определена, то можно воспользоваться свойством
M-1 = MT
• И получить:
Lij = Mij sqrt( Λj )
19
Факторный анализ
Исх. признаки
• При вычислении матрицы L оставим не все λ, а только
самые большие. Если λj< λ*, пусть λj =0
• Таким образом, после стандартизации данных (которая
включает приведение среднего по каждому признаку к 0),
вычисляем всё, что нужно:
Факторные
cov_m=cov(x);
call eigen(eVal, eVec, Cov_m);
NumPC = 2;
B=(j(ncol(x),Numpc,1)#T(sqrt(eval[1:NumPC])));
loadings = eVec[,1:NumPC] # B ;
/* loadings_ij=eigenvector_ij * sqrt(eigenvalue_j) */
psi=(cov_m-loadings*T(loadings)); *specific variances;
communality= loadings[,##];
нагрузки
L
20
Факторный анализ
•
•
•
•
Интерпретация результатов * :
Loadings (факторные нагрузки) - элементы матрицы L - вклад в новые
факторы («главные компоненты») каждого исходного признака. Коэффициент
корреляции между новым фактором и исходным признаком.
Communality (относительная дисперсия) – сумма квадратов факторных
нагрузок для данного признака. Имеют смысл коэффициента детерминации
при предсказании исходных признаков из новых факторов (формула как в
регрессии). Хорошо ли они предсказываются, или не очень (мы что-то
потеряли в при избавлении от «главных компонент»)
Specific variances (??) – информация о том, насколько хорошо оставленные в
модели факторы объясняют данные (чем ближе к нулю – тем лучше).
* http://sites.stat.psu.edu/~ajw13/stat505/fa06/17_factor/07_factor_commun.html
21
Факторный анализ
• SAS approach: PROC FACTOR, результаты те
же:
proc factor method=principal n=2 data=tmp22;
var english--vocab;
run;
22
Факторный анализ
<-Факторные нагрузки
<-Дисперсия данных,
объясненная каждым
фактором (соб. значения
матрицы ковариаций)
<-Относит. дисперсия
для каждого признака
<- матрица Ψ
23
Определение пороговой вероятности для
логистической регрессии
• Для того, чтобы перевести предсказанную вероятность целевого
события p(X) в бинарный отклик, необходимо определить пороговую
вероятность.
Logit(p(X))=a+b*x1+c*x2+….
если p>p*, то предсказываем целевое событие
если p<=p*, то предсказываем альтернативное событие
• Для выбора p* есть много соображений, например, можно применить
подход, основанный на матрице прибыли.
Истинный исход ->
Исход моделиров. |
V
Клиент – хороший
Клиент - жулик
Клиент – хороший,
даём $$$
Прибыль 100$
Потери -200$
Клиент – жулик, не
даём $$$
Потери -25$
Прибыль 0$
24
Определение пороговой вероятности для
логистической регрессии
• Давайте реализуем подсчет функционала «прибыли» на IML
(потенциально – очень гибкое решение, которое, в том числе,
позволяет ввести стратификацию по различным признакам, или более
хитрую логику при вычислении этого функционала)
• Сначала построим модель (вы уже знаете, как это делать!) и
посчитаем для всех наблюдений вер-ть целевого события:
proc logistic data=mycsv;
class laufkont moral;
model kredit(event='1')= laufkont laufzeit moral hoehe
score out=my_super_model;
alter;
25
Определение пороговой вероятности для
логистической регрессии
• Замечание : рекомендуется проводить анализ, показанный ниже, для
выборки, не участвовавшей в построении модели.
/* actually g(1) b(0) */
profit_m={100
-200, /*pred good (1)*/
-25
0};
/*pred bad (0)*/
use my_super_model;
read all var {kredit p_1} into x;
close my_super_model;
a=i(2);
plot=j(1,2,0);
do p_t=0 to 1 by 0.01;
mod_r=(x[,2]>p_t); *result predicted by our model at given p_t;
count_11=sum(((mod_r=1) & (x[,1]=1)));
count_10=sum(((mod_r=1) & (x[,1]=0)));
count_01=sum(((mod_r=0) & (x[,1]=1)));
Считаем правильно/неправильно
count_00=sum(((mod_r=0) & (x[,1]=0)));
классифицированные случаи
a[1,1]=count_11;
a[1,2]=count_10;
a[2,1]=count_01;
a[2,2]=count_00;
Вычислим прибыль и
profit=(a # profit_m)[+,+];
plot=plot//(p_t||profit);
Добавим строку в матрицу (в духе матлаб)
end;
26
Profit!!
Profit
p_t
27
Замечания.
• Не учтено, что обучающая выборка может (в некоторых
случаях - должна) отличаться от генеральной совокупности
в пропорциях «плохих» и «хороших». Соотв., нужно подругому считать предсказанные вероятности для ген.
совокупности.
• В этом случае в формулу для logit(p) надо вводить
поправку *
τ – доля целевых событий в ген. совокупности,
y – их доля в обучающей выборке
• Самый простой способ – опция PRIOREVENT=value в
операторе SCORE proc Logistic (see help). Она задаёт
значение τ.
* напр., см. http://gking.harvard.edu/files/0s.pdf, ф-ла (7)
28
1/--страниц
Пожаловаться на содержимое документа