Квартет Энскомба

Для подтверждения это мысли рассмотрим такой курьезный пример, как “Квартет Энскомба” – специально подобранные в 1973 году данные английским математиком Ф. Дж. Энскомбом для иллюстрации важности применения графиков для статистического анализа, и влияния выбросов значений на свойства всего набора данных. Эти данные состоят из четырёх пар \(x\) и \(y\) с практически равным средним значением (\(M[x_i] = 9\), \(M[y_i] = 7.5\)) и дисперсией между соответствующими элементами пар (\(D[x_i] = 11\), \(D[y_i]\approx 4.13\)) , а также равным коэффициентом корреляции (\(cor(x_i,y_i) = 0.816\)). Модель линейной регрессии, построенная методом МНК для всех вариантов описывается уравнением \(y = 3.00 + 0.500x\). Графики представлены на рисунке ниже, из которого видно, насколько могут различаться данные, описываемые внешне статистически одинаково.

Код для построения графиков и графики приводятся ниже.

# Загружаем данные и выводим загруженную таблицу 
load(url("https://tushavin.ru/RStudio/Ansc.Rda"))
knitr::kable(Ansc)

Рассчитываем статистику

options(digits=3) # Устанавливаем вывод 3 знаков
apply(Ansc,2,mean) # считаем средние для колонок

##  x1  x2  x3  x4  y1  y2  y3  y4 
## 9.0 9.0 9.0 9.0 7.5 7.5 7.5 7.5

apply(Ansc,2,var)  # считаем дисперсию колонок

##    x1    x2    x3    x4    y1    y2    y3    y4 
## 11.00 11.00 11.00 11.00  4.13  4.13  4.12  4.12

attach(Ansc) # Позволяет обращаться к колонкам по названию столбца
#считаем корелляцию между x и y для каждой пары
cat(cor(x1,y1),cor(x2,y2),cor(x3,y3),cor(x3,y3))

## 0.816 0.816 0.816 0.816

lm(y1~x1)

## 
## Call:
## lm(formula = y1 ~ x1)
## 
## Coefficients:
## (Intercept)           x1  
##         3.0          0.5

Выводим графики

# вывод 4 графиков на лист и смещение границ
# настройки сохраняем
oldpar<-par(mfrow=c(2,2),mar=c(4,4,1,1))
plot(y1~x1,xlab="X",ylab="Y",xlim=c(4,19),ylim=c(4,13),pch=19)
abline(a=3,b=0.5)
plot(y2~x2,xlab="X",ylab="Y",xlim=c(4,19),ylim=c(4,13),pch=19)
abline(a=3,b=0.5)
plot(y3~x3,xlab="X",ylab="Y",xlim=c(4,19),ylim=c(4,13),pch=19)
abline(a=3,b=0.5)
plot(y4~x4,xlab="X",ylab="Y",xlim=c(4,19),ylim=c(4,13),pch=19)
abline(a=3,b=0.5)

Построим регрессионные модели для каждого из четырех случаев

summary(lm(y1~x1)) # Первая пара

## 
## Call:
## lm(formula = y1 ~ x1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.9213 -0.4558 -0.0414  0.7094  1.8388 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)    3.000      1.125    2.67   0.0257 * 
## x1             0.500      0.118    4.24   0.0022 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.24 on 9 degrees of freedom
## Multiple R-squared:  0.667,  Adjusted R-squared:  0.629 
## F-statistic:   18 on 1 and 9 DF,  p-value: 0.00217

summary(lm(y2~x2)) # Вторая пара

## 
## Call:
## lm(formula = y2 ~ x2)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -1.901 -0.761  0.129  0.949  1.269 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)    3.001      1.125    2.67   0.0258 * 
## x2             0.500      0.118    4.24   0.0022 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.24 on 9 degrees of freedom
## Multiple R-squared:  0.666,  Adjusted R-squared:  0.629 
## F-statistic:   18 on 1 and 9 DF,  p-value: 0.00218

summary(lm(y3~x3)) # Третья пара

## 
## Call:
## lm(formula = y3 ~ x3)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -1.159 -0.615 -0.230  0.154  3.241 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)    3.002      1.124    2.67   0.0256 * 
## x3             0.500      0.118    4.24   0.0022 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.24 on 9 degrees of freedom
## Multiple R-squared:  0.666,  Adjusted R-squared:  0.629 
## F-statistic:   18 on 1 and 9 DF,  p-value: 0.00218

summary(lm(y4~x4)) # Четвертая пара

## 
## Call:
## lm(formula = y4 ~ x4)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -1.751 -0.831  0.000  0.809  1.839 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)    3.002      1.124    2.67   0.0256 * 
## x4             0.500      0.118    4.24   0.0022 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.24 on 9 degrees of freedom
## Multiple R-squared:  0.667,  Adjusted R-squared:  0.63 
## F-statistic:   18 on 1 and 9 DF,  p-value: 0.00216

Обратимся к моделям. Все четыре примера проходят по критерию Фишера, коэффициенты значимы, коэффициент детерминации везде 0,67, однако на графиках видно, что нормальная регрессия у нас только в первом случае.

Во втором случае у нас явно квадратичная зависимость, а в третьем и четвертом помахи (выбросы) влияют на уровнение регрессии и при их удалении результат будет совершенно другой.

par(oldpar)       # возвращаем сохраненные настройки
options(digits=7) # Значение по умолчанию
detach(Ansc)      # Отсоединяем имена таблиц
rm(Ansc)          # Удаляем таблицу из памяти

Коэффициент корреляции и графики регрессии

Как вы должны помнить, для модели парной линейной регрессии коэффициент детерминации \(R^2\) равен квадрату обычного коэффициента корреляции между y и x.

Напоминаю, что коэффициент корреляции, сокращенно r, — это число от -1 до 1, которое отражает силу линейной связи между двумя числовыми переменными. Например, допустим, вы опросили 30 человек об их весе и росте и изобразили эти 30 пар (вес, рост) на диаграмме рассеяния.

Если все 30 точек данных идеально ложатся на возрастающую линию, то корреляция между этими двумя переменными будет равна r = 1.

Если же общая форма зависимости (вес, рост) возрастающая, но 30 точек данных не идеально ложатся на одну линию, то r будет где-то между 0 и 1; чем ближе точки данных к прямой линии, тем ближе к 1 будет r.

Если зависимость убывающая, то r будет лежать между 0 и -1, а если линейная зависимость между весом и ростом вообще отсутствует, r будет равен 0.

Один коэффициент корреляции может отражать любое количество закономерностей. Коэффициенты корреляции популярны среди исследователей, поскольку они позволяют обобщить взаимосвязь между двумя переменными в одном числе. Однако данный коэффициент корреляции может отражать любое количество закономерностей между двумя переменными, и без дополнительной информации (в идеале в виде диаграммы рассеяния), как видно из предыдущего примера, понять ее невозможно.

Чтобы проиллюстрировать это, Яном Ванховым (Jan Vanhove) была написана функция R, plot_r(), которая принимает на вход коэффициент корреляции и размер выборки и выводит 16 совершенно разных диаграмм рассеяния, которые все характеризуются одним и тем же коэффициентом корреляции.

Для пояснения графиков используется материал указанного автора «What data patterns can lie behind a correlation coefficient?«.

Коэффициент корреляции равен 0.5 (r=0.5)

Построим такие графики с коэффициентом корреляции 0.5 и размером выборки 35.

# Для установки пакета используете две строки ниже без #
# library(devtools)
# install_github("janhove/cannonball")

library(cannonball)
plot_r(r = 0.5, n = 35)

Верхний ряд

Обычно, что когда студенты думают о взаимосвязи с коэффициентом корреляции 0,5, они представляют себе что-то вроде графиков (1) и (2). На обоих графиках базовая зависимость между X и Y линейна, а значения Y нормально распределены относительно наилучшим образом подходящей прямой линии. Небольшое различие между (1) и (2) заключается в том, что для (1) X распределен нормально, а для (2) X распределен равномерно. Эти два графика отражают тот тип отношений, который должен был отразить r.

График (3) отличается от (1) и (2) тем, что переменная X теперь взята из смещенного распределения. В этом случае большинство значений X сравнительно малы, но одно значение X довольно велико. Такое распределение может быть, когда X представляет собой, например, результаты выполнения участниками трудного задания (эффект пола). В этом случае одно или несколько выходящих за рамки, но истинных значений X могут иметь (не будут иметь) “высокий рычаг”, то есть они могут неоправданно повлиять на коэффициент корреляции, вытянув его вверх или вниз.

Проблема в (4) похожа на проблему в (3), но теперь большинство значений X сравнительно большие, а несколько — довольно низкие, возможно, потому что X отражает выполнение участниками задания, которое было для них слишком легким (эффект потолка). Здесь также аутсайдеры могут иметь “высокий рычаг”, т.е. они могут чрезмерно влиять на коэффициент корреляции, так что он не будет точно характеризовать основную часть данных.

Второй ряд

Графики (5) и (6) представляют собой вариации на ту же тему, что и графики (3) и (4): Значения Y не распределены нормально относительно линии регрессии, а смещены. В таких случаях некоторые отклоняющиеся, но истинные значения Y могут (не обязательно) иметь “высокий рычаг”, т.е. они могут тянуть коэффициент корреляции вверх или вниз гораздо сильнее, чем обычные точки данных.

Графики (7) и (8) — это два примера, когда изменчивость значений Y относительно прямой линии увеличивается и уменьшается, соответственно, по мере увеличения X, хотя, конечно, в данном примере это не очень понятно. Это явление известно как гетероскедастичность. Основные проблемы со слепым полаганием на коэффициенты корреляции в присутствии гетероскедастичности, на мой взгляд, заключаются в том, что (а) “r = 0,5” одновременно занижает то, насколько хорошо Y можно оценить по X для низких (высоких) значений X, и завышает то, насколько хорошо Y можно оценить по X для высоких (низких) значений X, и (б) просто сообщая коэффициент корреляции, вы упускаете важный аспект данных. Кроме того, гетероскедастичность может повлиять на вашу инференциальную статистику (инференциальная статистика — это отрасль статистики, которая делает выводы о соответствующей популяции из набора данных, полученных из выборки, подвергнутой случайным, наблюдательным и выборочным вариациям. Как правило, результаты получают из случайной выборки населения, а выводы, полученные из выборки, затем обобщают для представления всей совокупности).

Третий ряд

График (9) иллюстрирует, что коэффициенты корреляции выражают силу линейной связи между двумя переменными. Если связь не линейная, то они малоинформативны. В данном случае r = 0,5 сильно занижает силу связи XY, которая оказывается нелинейной (в данном случае квадратичной). То же самое относится и к (10), где r = 0,5 занижает силу связи XY и упускает из виду циклический характер связи.

Графики (11) и (12) иллюстрируют, как единственная помеха, например, из-за технической ошибки, может дать вводящие в заблуждение коэффициенты корреляции. В (11) однин выброс данных дает сильную положительную корреляцию; если бы учитывались только данные слева, наблюдалась бы отрицательная связь (пунктирная красная линия). Слепое использование r = 0,5 неверно характеризует большую часть данных. В (12) зависимость значительно сильнее, чем r = 0,5 для основной массы данных (пунктирная красная линия); выброс снижает коэффициент корреляции. Графики (11) и (12) отличаются от графиков (3) и (4) тем, что в графиках (3) и (4) все значения X были взяты из одного и того же, но смещенного распределения и, как таковые, являются настоящими точками данных; в графиках (11) и (12) выбросы были вызваны механизмом, отличным от других точек данных (например, ошибкой кодирования или техническим сбоем).

Четвертый ряд

В (13) значения Y распределены бимодально относительно линии регрессии. Это говорит о том, что мы упустили из виду какой-то важный аспект данных, например, фактор группировки: возможно, точки данных выше линии регрессии были отобраны из другой популяции, чем точки ниже линии регрессии.

Ситуация на графике (14) похожа на ситуацию в (13), но значительно хуже ее: Набор данных содержит две группы, но в отличие от (13), общая тенденция, отражаемая r = 0,5, скрывает тот факт, что внутри каждой из этих групп связь XY на самом деле отрицательная. График (14) часто, но не всегда, дает такую картину, которая известна как парадокс Симпсона.

На графике (15) показана ситуация, когда исследователи, вместо того чтобы изучать взаимосвязь XY по всему диапазону X, изучали только случаи с самыми экстремальными значениями X. Выборка по крайним значениям завышает коэффициенты корреляции (см. причину № 2, по которой я не очень люблю коэффициенты корреляции). Другими словами, если вы возьмете выборку из 150 случаев XY и посмотрите только на 50 самых экстремальных значений X, вы получите коэффициент корреляции, который с большой вероятностью будет больше, чем тот, который вы наблюдали бы, если бы рассмотрели все 150 случаев.

Наконец, график (16) — это то, что на самом деле представляют собой многие коэффициенты корреляции. Например, данные X и Y неровные, потому что они представляют собой данные подсчета или ответы на вопросы анкеты. Не факт, что коэффициенты корреляции для таких моделей сами по себе обманчивы, но мы явно говорим о другой модели, чем на графиках (1) и (2).

А если коэффициент корреляции равен нулю (r=0)?

plot_r(r = 0, n = 35)

Главное, что видно из этих графиков, это то, что r = 0 не обязательно означает отсутствие связи XY. Это ясно из графиков (9) и (10), которые демонстрируют сильную нелинейную зависимость. Графики (11) и (12) также подчеркивают этот момент: Существует сильная взаимосвязь для большей части данных, но эта тенденция аннулируется одной точкой данных. Иногда, как показано на графике (14), тенденция, присутствующая в двух подгруппах, может быть не видна в агрегированном анализе; однако в данном примере это не так.

А если коэффициент корреляции равен 0,75 (r=0,75)?

Как мы помним, коэффициент детерминации для модели с константой принимает значения от 0 до 1. Чем ближе значение коэффициента к 1, тем сильнее зависимость. При оценке регрессионных моделей это интерпретируется как соответствие модели данным. Для приемлемых моделей предполагается, что коэффициент детерминации должен быть хотя бы не меньше 50 % (в этом случае коэффициент множественной корреляции превышает по модулю 70 %).

Посмотрим на такие графики

plot_r(r = 0.75, n = 35)

Комментарии, надеюсь, излишни. Пояснения к рисунку аналогичны пояснениям к первому примеру с r=0.5

Мораль: мало построить регрессионную модель, надо еще построить ее график.

Разберем задачу.

В интернет-магазин поступает заказ. При получении заказа, если товар есть, то выписывается счет. Если товара нет, то предлагается зарезервировать товар. Если товар зарезервирован, то по поступлению товара на склад выписывается счёт.  После оплаты счета товар отгружается, если оплаты нет 10 дней, то заказ аннулируется.

Для начала, определимся, что является событиями, что действиями, а что условиями (пока только явно описанные в задаче):

События (Events):

• поступление заказа;
• поступление товара;;
• оплата счета;
• отсутствие оплаты 10 дней.

Действия (Actions)

• выписать счёт;
• предложить зарезервировать товар;
• отгрузить товар;
• аннулировать заказ;

Шлюз/условие (Gateway)

• товар имеется в наличии?

Вроде все просто, правда ведь? А теперь посмотрим на детский рисунок на заданную тему…

Сразу видно, что события перепутаны с действиями. Крайне распространенная ошибка, в большинстве представленных схем она присутствовала, просто не так явно. Например, другая схема:

Как видим, «Поступление заказа» в ней трактуется как событие. Сложнее с другими ошибками, разберем их, пока не придираясь особо к тексту на рисунке 2. Представьте себе мысленно муравьишку, который бежит по нашей бизнес-схеме. Муравьишка запускается по первому событию, спускаем его с цепи и он побежал.

Вариант 1. Товар есть.

На предложенной схеме товар никогда не будет отгружен. Мы будем выписывать счета и практически мгновенно их аннулировать. Возникает вопрос, что не так с условием «Счет оплачен в течение 10 дней»? Дело в том, что для условий по событию нотация BPMN предусматривает несколько иную конструкцию, которая правильно показана на схеме ниже, а именно, т.н. «шлюз, управляемый событиями». Если мы берем в данном случае обычную развилку типа «или», то она обрабатывается без задержки, а в момент визита на развилку муравьишки счёт гарантировано не оплачен, поскольку мы его выписали секунду назад и даже еще не направили покупателю, поэтому наверх к действию «товар отгружается» (кстати, в таком написании текста это однозначно событие, а не действие, еще одна ошибка) муравьишка никогда не попадет. Посочувствуем скотинке.

Спрашивается, а что не так с рисунком 3. Вроде же все хорошо, если не придираться к тексту. А вот и нет. Начнем с самого начала. Два события подряд, первое лишнее. Убираем. Но это так, мелкие придирки. А вот дальше мы проверяем наличие товара и нам оказывается, согласно схеме, всё равно: готов клиент ждать или не готов. Мы в любом случае резервируем товар.  Т.е. действие «Предложить зарезервировать товар» из списка выше просто опущено, а это грубая ошибка и нарушение условий задачи. Ну, а дальше считается, что клиент согласен ждать вечно, когда товар поступит, поскольку событие «Товар поступил» в реальной жизни может наступить очень поздно, либо вообще не наступить. Соглашусь, что это не описано в условиях, но кто обещал, что будет легко? Тем более, что на лекции я несколько раз это рассказывал, что в реальной жизни никто вам подробно все не расскажет. В этом и состоит задача бизнес-аналитика, выявить все такие моменты. Кстати, очень перекликается с управлением рисками, не находите?

Были отдельные граждане, которые умудрились предусмотреть вроде всё, но в процессе наделали другие ошибки.

Во всяком случае, из схемы видно, что разум кипел в процессе рождения этого шедевра.

Перечислю тезисно ошибки схемы, в соответствии с красными цифрами на ней.

1. Дубликат события. Лишняя сущность. Ошибка не грубая, но некрасиво.
2. Развилка. А куда когда идем? Непонятно
3. Условие управляемое событиями. Можно, конечно, и так, но это как в окно выходить при наличии двери. Достаточно было обычного условия «или» в данном случае и не выпендриваться.
4. Такие события очень коварны на схемах. Что тут не так, я описал в предыдущем примере.
5. Действие «выписать счет» должно быть всего одно. Это одно и то же действие, с одним и тем же ответственным и исполнителем. Дубликат в данном случае просто методически неверен с точки зрения нотации.
6. Нотация не предусматривает такое использование шлюза, управляемого событиями. К тому же совершенно лишний элемент на схеме.
7. И снова дубликат действия. Теперь «Аннулировать заказ»
8. Событие «оплата счета» между двумя условиями, с последующей проверкой оплачен ли счет. Полная ерунда.

Итак, как надо было действовать, при решении этой задачи

1. Изучить нотацию и примеры, ссылка на которые была в задании.
2. Определить начало и конец процесса. С этого начинается моделирование любого процесса. Так мы обозначаем рамки, в которых будем работать. Все начинается с события и заканчивается событием.
3. Для начала лучше всего описать линейную последовательность действий: шаг за шагом движение от начала к финальному результату, то, что называется «happy path». Постепенно добавляются ветвления. В таком порядке работать намного проще, чем ставить две или более ветвей одновременно и путаться в стрелках, что откуда и куда идет.
4. Подпроцессов должно быть столько, чтобы избежать ненужной детализации, но не более того. Помните о чувстве меры. Если подпроцессов будет слишком мало, то действия, которые стоило бы спрятать в них, будут находиться в общем процессе, создавая дополнительные объекты, стрелки, ветвления и, как следствие, путаницу. Если вы перестараетесь с желанием убрать все в подпроцессы, то диаграмма потеряет свою информативность, а какие-то изменения в подпроцессе начнут ненаглядно влиять на результаты всего процесса.
5. Все названия процессов должны быть максимально информативны и понятны. Иначе читабельность диаграммы также будет крайне низкой. По методике названия см. выше.
6. Определяем ответственных лиц. До этого мы работали с событиями «в чистом виде». Теперь у них появляются исполнители и ответственные.
7. Добавляем данные, сноски, комментарии, если это необходимо.

Я сознательно не привожу правильную схему, поскольку не все еще выполнили это задание, но, думаю, того что я написал должно хватить для понимания, как не надо делать.

Sapienti Sat.

Если данный пакет не установлен, его необходимо предварительно установить.

Подключаем пакет.

library(sets)

Для примера рассмотрим простую систему терморегулятора на основе нечеткой логики. Возьмем в качестве комфортной температуру 18-24 градуса. Температуру ниже 20 градусов будем считать холодной (COLD), а выше 22 градусов жаркой. Температуру от 18 до 24 будем считать комфортной.

Рассмотрим также состояние клапана регулятора батареи отопления, который может быть закрыт, частично закрыт, частично открыт и открыт. Опишем систему с помощью команд пакета sets.

Для начала задам “вселенную”, т.е. то пространство, на котором будет проводиться вычисление.

sets_options("universe", seq(from = 0, to = 100, by = 0.1))

Теперь опишем обе переменные, используя в первом случае трапецоиды, во втором треугольники.

variables <-set(
    TEMP=fuzzy_variable(COLD=fuzzy_trapezoid(corners = c(-1,0,17,20)),
                    NORMAL=fuzzy_trapezoid(corners = c(18,19,21,24)),
                    HOT=fuzzy_trapezoid(corners = c(22,25,100,101))),
    VALVE = fuzzy_variable(CLOSED=fuzzy_triangular(corners = c(-1,0,25.1)),
                    PART.CLOSED=fuzzy_triangular(corners = c(0,25,75.1)),
                    PART.OPENED=fuzzy_triangular(corners = c(26,75.2,100)),
                    OPENED=fuzzy_triangular(corners = c(75.1,100,101)))
 )

Опишем простые правила регулирования клапана. Если холодно, то открыть, если жарко, то закрыть, если комфортно то не менять. Обратите внимание, что надо предусмотреть все возможные состояния системы. Если мы что-то пропустим, то результат может нас неприятно удивить.

rules <-set(
    fuzzy_rule(TEMP %is% COLD, VALVE %is% OPENED),
    fuzzy_rule(TEMP %is% HOT, VALVE %is% CLOSED),
    fuzzy_rule(TEMP %is% NORMAL &&  VALVE %is% CLOSED ,VALVE %is% CLOSED),
    fuzzy_rule(TEMP %is% NORMAL &&  VALVE %is% PART.CLOSED, VALVE %is% PART.CLOSED),
    fuzzy_rule(TEMP %is% NORMAL &&  VALVE %is% PART.OPENED, VALVE %is% PART.OPENED),
    fuzzy_rule(TEMP %is% NORMAL &&  VALVE %is% OPENED, VALVE %is% OPENED)
    )

Сформируем систему и выведем о ней данные.

system <- fuzzy_system(variables, rules)
print(system)

## A fuzzy system consisting of 2 variables and 6 rules.
## 
## Variables:
## 
## TEMP(COLD, NORMAL, HOT)
## VALVE(CLOSED, PART.CLOSED, PART.OPENED, OPENED)
## 
## Rules:
## 
## TEMP %is% NORMAL && VALVE %is% CLOSED => VALVE %is% CLOSED
## TEMP %is% NORMAL && VALVE %is% OPENED => VALVE %is% OPENED
## TEMP %is% NORMAL && VALVE %is% PART.CLOSED => VALVE %is% PART.CLOSED
## TEMP %is% NORMAL && VALVE %is% PART.OPENED => VALVE %is% PART.OPENED
## TEMP %is% HOT => VALVE %is% CLOSED
## TEMP %is% COLD => VALVE %is% OPENED

plot(system)

Получилось, мягко говоря, не ахти. Выведем отдельный графики с русскими подписями осей.

vars<-as.list(system$variables)
nams <- names(vars)
plot(vars[[2]], main = nams[2],col=rainbow(4),xlab="Температура, град. C",ylab="Степень принадлежности")

Теперь проверим как работает система. Зададим начальные условия.

data<-list(TEMP=5,VALVE=0)
fc1<-fuzzy_inference(system,data)
plot(fc1)

Проведем дефаззификацию, под которой в системах нечеткого вывода понимают процесс перехода от функции принадлежности выходной лингвистической переменной к её четкому (числовому) значению. Обращаю ваше внимание, что это  всего лишь один из возможных методов.

gset_defuzzify(fc1, "centroid")

## [1] 91.73333

Получили результат, что клапан открыт на 91.7%. ПУсть прошло какое-то время и при таком положении клапана воздух в помещении нагрелся до 17 градусов.

data<-list(TEMP=17,VALVE=91.7)
fc1<-fuzzy_inference(system,data)
plot(fc1)

gset_defuzzify(fc1, "centroid")

## [1] 91.73333

Ничего не поменялось. Нагрелось до 18 градусов:

data<-list(TEMP=18,VALVE=91.7)
fc1<-fuzzy_inference(system,data)
plot(fc1)

gset_defuzzify(fc1, "centroid")

## [1] 91.03544

Мы прикрыли клапан до 91%. Температура выросла еще на 1 градус.

data<-list(TEMP=19,VALVE=91)
fc1<-fuzzy_inference(system,data)
plot(fc1)

gset_defuzzify(fc1, "centroid")

## [1] 70.63561

Температура стабилизировалась, клапан ещё чуть прикрыт.

data<-list(TEMP=19,VALVE=70.6)
fc1<-fuzzy_inference(system,data)
plot(fc1)

gset_defuzzify(fc1, "centroid")

## [1] 64.90177

Температура падает.

data<-list(TEMP=18.1,VALVE=64.9)
fc1<-fuzzy_inference(system,data)
plot(fc1)

gset_defuzzify(fc1, "centroid")

## [1] 69.54573

Тогда немного приоткроем клапан…

Имитационная модель

Рассмотрим некое абстракное помещение, которое мы обогреваем. Пусть теплопотери будут прямо пропорциональны разнице температур между наружним и внутренним помещением с учетом площади поверхности и коэффициента теплопотерь.

Пусть наружняя температура равна 0, а теплопотери в минуту равны T/30, где Т — температура в помещении. Пусть нагреватель при ста процентах мощности нагревает воздух на 1 градус за 1 минуту.

Промоделируем работу двух систем с использованием различных подходов к дефаззификации.

result<-data.frame(time=1:100,TEMP=rep(0.0,100),K=rep(0.0,100),TEMP1=rep(0.0,100),K1=rep(0.0,100))
for(i in 1:99) {
  fc1<-fuzzy_inference(system,list(TEMP=result$TEMP[i],VALVE=result$K[i]))  
  result$K[i+1]<-round(gset_defuzzify(fc1, "centroid"),1)
  result$TEMP[i+1]<-round(result$TEMP[i]+(1-result$TEMP[i+1]/30)*result$K[i+1]/100,1)
  fc1<-fuzzy_inference(system,list(TEMP=result$TEMP1[i],VALVE=result$K1[i]))  
  result$K1[i+1]<-round(gset_defuzzify(fc1, "meanofmax"),1)
  result$TEMP1[i+1]<-round(result$TEMP1[i]+(1-result$TEMP1[i+1]/30)*result$K1[i+1]/100,1)
  
}
plot(TEMP~time,data=result,col="red",xlab="Время, мин.", ylab="Температура, градусы.",type="l")
lines(result$TEMP1, col="blue")

plot(K~time,data=result,col="red",xlab="Время, мин.", ylab="Открытие вентиля, проценты.",type="l")
lines(result$K1, col="blue")

Обратите внимание, что модель с центроидом не стабилизируется.

Примерно так всё оно и работает и в реальной жизни, поэтому надо достаточно вдумчиво подходить к проектированию таких систем.

1. Основные преобразования таблиц описаны вот здесь.
2. Тест для повторения усвоенного нужно пройти вот здесь.
3. Онлайн редактор формул в формете mathtex есть вот тут. Краткое описание на русcком вот тут.  А вот этот редактор формул для Word я упоминал.

Для закрепления рекомендую набрать формулу:

\[ f_{\mathbf{X}}(\mathbf{x}) = \frac{1}{(2\pi )^{n/2} \vert \Sigma \vert^{1/2}} e^{-\frac{1}{2}(\mathbf{x} — \mathbf{\mu})^{\top} \Sigma^{-1} (\mathbf{x} — \mathbf{\mu})},\; \mathbf{x} \in \mathbb{R}^n \]

Примечание: ответ для сампопроверки содержится в самой формуле, для этого достаточно с помощью правой кнопкой мыши вывести код Tex. Это не единственный правильный вариант, можно написать и несколько иначе.

Давайте рассмотрим несколько рецептов на примере данных mtcars. Напоминаю, что команда head выводит первые несколько строк большой таблицы. Чтобы не перегружать заметку я её буду использовать практически в каждой команде.

data("mtcars")
head(mtcars)

##                    mpg cyl disp  hp drat    wt  qsec vs am gear carb
## Mazda RX4         21.0   6  160 110 3.90 2.620 16.46  0  1    4    4
## Mazda RX4 Wag     21.0   6  160 110 3.90 2.875 17.02  0  1    4    4
## Datsun 710        22.8   4  108  93 3.85 2.320 18.61  1  1    4    1
## Hornet 4 Drive    21.4   6  258 110 3.08 3.215 19.44  1  0    3    1
## Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02  0  0    3    2
## Valiant           18.1   6  225 105 2.76 3.460 20.22  1  0    3    1

head(mtcars[order(mtcars$mpg, mtcars$cyl, mtcars$hp),])

##                      mpg cyl disp  hp drat    wt  qsec vs am gear carb
## Cadillac Fleetwood  10.4   8  472 205 2.93 5.250 17.98  0  0    3    4
## Lincoln Continental 10.4   8  460 215 3.00 5.424 17.82  0  0    3    4
## Camaro Z28          13.3   8  350 245 3.73 3.840 15.41  0  0    3    4
## Duster 360          14.3   8  360 245 3.21 3.570 15.84  0  0    3    4
## Chrysler Imperial   14.7   8  440 230 3.23 5.345 17.42  0  0    3    4
## Maserati Bora       15.0   8  301 335 3.54 3.570 14.60  0  1    5    8

head(with(mtcars,mtcars[order(mpg, cyl, hp),]))

##                      mpg cyl disp  hp drat    wt  qsec vs am gear carb
## Cadillac Fleetwood  10.4   8  472 205 2.93 5.250 17.98  0  0    3    4
## Lincoln Continental 10.4   8  460 215 3.00 5.424 17.82  0  0    3    4
## Camaro Z28          13.3   8  350 245 3.73 3.840 15.41  0  0    3    4
## Duster 360          14.3   8  360 245 3.21 3.570 15.84  0  0    3    4
## Chrysler Imperial   14.7   8  440 230 3.23 5.345 17.42  0  0    3    4
## Maserati Bora       15.0   8  301 335 3.54 3.570 14.60  0  1    5    8

library(dplyr)
head(arrange(mtcars, mpg, cyl, hp))

##    mpg cyl disp  hp drat    wt  qsec vs am gear carb
## 1 10.4   8  472 205 2.93 5.250 17.98  0  0    3    4
## 2 10.4   8  460 215 3.00 5.424 17.82  0  0    3    4
## 3 13.3   8  350 245 3.73 3.840 15.41  0  0    3    4
## 4 14.3   8  360 245 3.21 3.570 15.84  0  0    3    4
## 5 14.7   8  440 230 3.23 5.345 17.42  0  0    3    4
## 6 15.0   8  301 335 3.54 3.570 14.60  0  1    5    8

library(tibble)
head(newcars<-rownames_to_column(mtcars))

##             rowname  mpg cyl disp  hp drat    wt  qsec vs am gear carb
## 1         Mazda RX4 21.0   6  160 110 3.90 2.620 16.46  0  1    4    4
## 2     Mazda RX4 Wag 21.0   6  160 110 3.90 2.875 17.02  0  1    4    4
## 3        Datsun 710 22.8   4  108  93 3.85 2.320 18.61  1  1    4    1
## 4    Hornet 4 Drive 21.4   6  258 110 3.08 3.215 19.44  1  0    3    1
## 5 Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02  0  0    3    2
## 6           Valiant 18.1   6  225 105 2.76 3.460 20.22  1  0    3    1

newcars %>% 
  group_by(cyl, am) %>%
  select(mpg, cyl, wt, am) %>%
  summarise(avgmpg = mean(mpg), avgwt = mean(wt)) %>%
  filter(avgmpg > 20)

## # A tibble: 3 x 4
## # Groups:   cyl [2]
##     cyl    am   avgmpg   avgwt
##   <dbl> <dbl>    <dbl>   <dbl>
## 1     4     0 22.90000 2.93500
## 2     4     1 28.07500 2.04225
## 3     6     1 20.56667 2.75500

library(magrittr)
head(X<-newcars %>% filter(cyl == 8))

##              rowname  mpg cyl  disp  hp drat   wt  qsec vs am gear carb
## 1  Hornet Sportabout 18.7   8 360.0 175 3.15 3.44 17.02  0  0    3    2
## 2         Duster 360 14.3   8 360.0 245 3.21 3.57 15.84  0  0    3    4
## 3         Merc 450SE 16.4   8 275.8 180 3.07 4.07 17.40  0  0    3    3
## 4         Merc 450SL 17.3   8 275.8 180 3.07 3.73 17.60  0  0    3    3
## 5        Merc 450SLC 15.2   8 275.8 180 3.07 3.78 18.00  0  0    3    3
## 6 Cadillac Fleetwood 10.4   8 472.0 205 2.93 5.25 17.98  0  0    3    4

library(stringr)
newcars %>% filter(str_detect(rowname,"RX4"))

##         rowname mpg cyl disp  hp drat    wt  qsec vs am gear carb
## 1     Mazda RX4  21   6  160 110  3.9 2.620 16.46  0  1    4    4
## 2 Mazda RX4 Wag  21   6  160 110  3.9 2.875 17.02  0  1    4    4

newcars %>% filter(str_detect(rowname,"Hornet"))

##             rowname  mpg cyl disp  hp drat    wt  qsec vs am gear carb
## 1    Hornet 4 Drive 21.4   6  258 110 3.08 3.215 19.44  1  0    3    1
## 2 Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02  0  0    3    2

newcars %>% filter(str_detect(rowname,"Hornet|RX4"))

##             rowname  mpg cyl disp  hp drat    wt  qsec vs am gear carb
## 1         Mazda RX4 21.0   6  160 110 3.90 2.620 16.46  0  1    4    4
## 2     Mazda RX4 Wag 21.0   6  160 110 3.90 2.875 17.02  0  1    4    4
## 3    Hornet 4 Drive 21.4   6  258 110 3.08 3.215 19.44  1  0    3    1
## 4 Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02  0  0    3    2

newcars %>% filter(str_detect(rowname,"Hornet|RX4") & cyl==6)

##          rowname  mpg cyl disp  hp drat    wt  qsec vs am gear carb
## 1      Mazda RX4 21.0   6  160 110 3.90 2.620 16.46  0  1    4    4
## 2  Mazda RX4 Wag 21.0   6  160 110 3.90 2.875 17.02  0  1    4    4
## 3 Hornet 4 Drive 21.4   6  258 110 3.08 3.215 19.44  1  0    3    1

head(Y<-newcars %>% 
  group_by(cyl, am) %>%
  select(mpg, cyl,am) %>%
  summarise(avgmpg = mean(mpg)))

## # A tibble: 6 x 3
## # Groups:   cyl [3]
##     cyl    am   avgmpg
##   <dbl> <dbl>    <dbl>
## 1     4     0 22.90000
## 2     4     1 28.07500
## 3     6     0 19.12500
## 4     6     1 20.56667
## 5     8     0 15.05000
## 6     8     1 15.40000

library(tidyr)

## 
## Attaching package: 'tidyr'

## The following object is masked from 'package:magrittr':
## 
##     extract

Y %>% spread(am,avgmpg)

## # A tibble: 3 x 3
## # Groups:   cyl [3]
##     cyl    `0`      `1`
## * <dbl>  <dbl>    <dbl>
## 1     4 22.900 28.07500
## 2     6 19.125 20.56667
## 3     8 15.050 15.40000

4. Ёще несколько полезных примеров

newcars %>% group_by(cyl, am) %>% summarise_all(c("mean", "sd"))

## Warning in mean.default(rowname): argument is not numeric or logical:
## returning NA

## Warning in var(if (is.vector(x) || is.factor(x)) x else as.double(x), na.rm
## = na.rm): в результате преобразования созданы NA

## # A tibble: 6 x 22
## # Groups:   cyl [?]
##     cyl    am rowname_mean mpg_mean disp_mean   hp_mean drat_mean  wt_mean
##   <dbl> <dbl>        <dbl>    <dbl>     <dbl>     <dbl>     <dbl>    <dbl>
## 1     4     0           NA 22.90000  135.8667  84.66667  3.770000 2.935000
## 2     4     1           NA 28.07500   93.6125  81.87500  4.183750 2.042250
## 3     6     0           NA 19.12500  204.5500 115.25000  3.420000 3.388750
## 4     6     1           NA 20.56667  155.0000 131.66667  3.806667 2.755000
## 5     8     0           NA 15.05000  357.6167 194.16667  3.120833 4.104083
## 6     8     1           NA 15.40000  326.0000 299.50000  3.880000 3.370000
## # ... with 14 more variables: qsec_mean <dbl>, vs_mean <dbl>,
## #   gear_mean <dbl>, carb_mean <dbl>, rowname_sd <dbl>, mpg_sd <dbl>,
## #   disp_sd <dbl>, hp_sd <dbl>, drat_sd <dbl>, wt_sd <dbl>, qsec_sd <dbl>,
## #   vs_sd <dbl>, gear_sd <dbl>, carb_sd <dbl>

newcars %>% summarise_if(is.numeric,c("mean", "sd"))

##   mpg_mean cyl_mean disp_mean  hp_mean drat_mean wt_mean qsec_mean vs_mean
## 1 20.09062   6.1875  230.7219 146.6875  3.596563 3.21725  17.84875  0.4375
##   am_mean gear_mean carb_mean   mpg_sd   cyl_sd  disp_sd    hp_sd
## 1 0.40625    3.6875    2.8125 6.026948 1.785922 123.9387 68.56287
##     drat_sd     wt_sd  qsec_sd     vs_sd     am_sd   gear_sd carb_sd
## 1 0.5346787 0.9784574 1.786943 0.5040161 0.4989909 0.7378041  1.6152

mtcars %>% group_by(cyl) %>% mutate(rank = min_rank(desc(mpg)))

## # A tibble: 32 x 12
## # Groups:   cyl [3]
##      mpg   cyl  disp    hp  drat    wt  qsec    vs    am  gear  carb  rank
##    <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
##  1  21.0     6 160.0   110  3.90 2.620 16.46     0     1     4     4     2
##  2  21.0     6 160.0   110  3.90 2.875 17.02     0     1     4     4     2
##  3  22.8     4 108.0    93  3.85 2.320 18.61     1     1     4     1     8
##  4  21.4     6 258.0   110  3.08 3.215 19.44     1     0     3     1     1
##  5  18.7     8 360.0   175  3.15 3.440 17.02     0     0     3     2     2
##  6  18.1     6 225.0   105  2.76 3.460 20.22     1     0     3     1     6
##  7  14.3     8 360.0   245  3.21 3.570 15.84     0     0     3     4    11
##  8  24.4     4 146.7    62  3.69 3.190 20.00     1     0     4     2     7
##  9  22.8     4 140.8    95  3.92 3.150 22.90     1     0     4     2     8
## 10  19.2     6 167.6   123  3.92 3.440 18.30     1     0     4     4     5
## # ... with 22 more rows