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

Для подтверждения это мысли рассмотрим такой курьезный пример, как “Квартет Энскомба” – специально подобранные в 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

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

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

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

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")

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

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

Давайте рассмотрим несколько рецептов на примере данных 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

I. Необходимо установить следующее ПО на домашний компьютер.

1. Язык R
2. RStudio

Богатая библиотека по R и RStudio на русском для самостоятельного изучения доступна вот тут.

II. Наобходимо повторить вот этот курс.  Он соответсвует тому, что было на первом занятии. По результатам курса выполнить контрольный пример. Результаты должны сойтись, если не сошлись — разбирайтесь. В конце-концов, у вас уже есть для этого высшее образование.

Дерево текущей реальности (ДТР)

Не вдаваясь в методику его построения, описанную в приведенной в списке литературы книге, рассмотрим простой пример сделанный на языке dot. Как мы видим нам нужны следующие графические элементы:

• прямоугольники
• овалы
• соединительные линии
• соединительные стрелки
• также дополнительно нужны кружки для многостраничной схемы

Все это делается следующим образом.

• прямоугольники: [shape = box];
• овалы: [shape = oval,label=»»,height=0.1];
• соединительные линии [dir=none];
• соединительные стрелки — и так есть;
• кружки: [shape=circle];

Также ориентируем документ сверху вниз (rankdir=BT) и выставляем шрифт.

Собственно, скрипт:

digraph PS3 {
resize=auto;
charset="utf8";
node [fontname="Verdana"];
rankdir=BT;
node [shape = box];
T102 [label="102. Сотрудники не до конца пони-\nмают значимость внедрения ПС"];
T301 [label="301.У руководства нет опыта\nвнедрения инструметов ПС"];
T302 [label="302. Для успешного внедрения\nнужно, чтобы каждый сотрудник \nбыл вовлечен в процесс "];
T303 [label="303. Для внедрения ПС\nтребуются доп. затраты ре \nсурсов"];
T304 [label="304. Ресурсы, направленные на мероприятия ПС\nограниченны"];
T305 [label="305. Руководство старается вовлечь \nво внедрение ПС\nбольшую часть сотрудников"];
T306 [label="306.Руководство желает \nвнедрить ПС"];
T307 [label="307.Руководство не хочет\nискать истинные \nпроблемы"];
T308 [label="308. Руководство\nобычно не проводит анализ \nпо степени значимости\nинициатив\nпо внедрению ПС "];
T309 [label="309. Отделы боятся\nнесправедливого\nраспределения ресурсов "];
T310 [label="310.Вовлечение\nвсего персонала\nв мероприятия ПС позволят\nраспределить равномерно ре\nсурсы "];
T311 [label="311.Не происходит\nраспределение мероприятий во важности, \nдля получения скорых результатов "];
T311А [label="311А. План мероприятий, \nкоторый не нацелен на\nрешение проблем предприятия, \nне приведет к улучшениям "];
T312 [label="312. Выделяемые ресурсы \nна внедрение ПС\nраспределены равномерно"];
T313 [label="313. Лишь некоторые мероприятия \nпо внедрению ПС \nповышают производительность и эффективность предприятия "];
T314 [label="314. Руководство \nне заостряет внимание\nна существующие проблемы \n(слабые места) "];
T315 [label="315.Нерациональное распределение \nресурсов уменьшает эффект\nпроводимых мероприятий ПС"];
T316 [label="316.Мероприятия, \nкоторые проводятся\nдля внедрения ПС,не\nулучшают производительность \nв целом "];
T317 [label="317. План мероприятий, \nкоторый не нацелен на\nрешение проблем предприятия, \nне приведет к улучшениям "];
T318 [label="318. Улучшения на предприятии\nне происходят, \nпока существуют\nнерешенные проблемы "];
T319 [label="319.Ресурсы тратятся\nна мероприятия, \nне влияющие на поизв-сть "];
T320 [label="320.Устранение проблем\n(слабых мест) оттягивают"];
T212 [label="212. Большинство думает, \nчто реальных результатов либо долго ждать, \nлибо совсем не будет"];
node [shape=circle];
"711"; "211";"112";
node [shape = oval,label="",height=0.1];
Z1;Z2;Z3;Z4;Z5;Z6;Z7;Z8;
T102->T301;
T301->Z1 [dir=none];
T301->T307;
T302->T305;
T303->"711";
T306->Z1 [dir=none];
T305->Z1 [dir=none];
Z1->T308;
Z1->T308;
Z1->T308;
T304->Z2 [dir=none];
T303->Z2 [dir=none];
T305->Z2 [dir=none];
Z2->T310;
Z2->T310;
Z2->T310;
T307->T311;
T308->T311;
T310->Z3 [dir=none];
T309->Z3 [dir=none];
Z3->T312;
Z3->T312;
T311A->Z4 [dir=none];
T311->Z4 [dir=none];
T305->Z4 [dir=none];
Z4->T313;
Z4->T313;
Z4->T313;
T312->T314;
T312->Z5 [dir=none];
T313->T316;
T314->Z6 [dir=none];
T315->Z6 [dir=none];
T313->Z5 [dir=none];
Z5->T319;
Z5->T319;
Z6->T320;
Z6->T320;
T319->"112";
T317->Z7 [dir=none];
T313->Z7 [dir=none];
Z7->T318;
Z7->T318;
T317->T311A;
T320->Z8 [dir=none];
T316->Z8 [dir=none];
T318->Z8 [dir=none];
Z8->T212;
Z8->T212;
Z8->T212;
"211"->T212;
}

Результат (кликабельно).

На мой взгляд, досаточно просто и гораздо быстрее, чем пытаться вручную двигать и вырванивать квадратики в графическом редакторе. Что такой graphviz и где его взять можно посмотреть в соответствующей заметке на сайте. По методике построения ДТР рекомендую обратиться к приведенной ниже книге.

Дерево будущей реальности (ДБР)

Построение принципиально не отличается от ДТР. Элементы те же.

Литература

1. Детмер У Теория ограничений Голдратта. Системный подход к непрерывному совершенствованию. М.: Альпина Паблишер, 2017. 855 с.

Отсутствие подтверждения статистических гипотез.

Если у вас раньше брак был 5 изделий из 100, а стал 3 изделия из 102 можно ли делать выводы, что брак уменьшился в 1.7 раза? Нет! Убедимся это простым расчетом в R.

> bad<-c(3,5)
> count<-c(100,102)
> prop.test(bad,count)

2-sample test for equality of proportions with continuity correction

data: bad out of count
X-squared = 0.11036, df = 1, p-value = 0.7397
alternative hypothesis: two.sided
95 percent confidence interval:
-0.08252672 0.04448751
sample estimates:
prop 1 prop 2
0.03000000 0.04901961

Warning message:
In prop.test(bad, count) :
аппроксимация на основе хи-квадрат может быть неправильной

Нулевая гипотеза в двнном случае, что пропорции одинаковы. Она подтверждается (p-value = 0.7397 > 0.05). Warning при расчетах говорит о том, выборка маловата для оценки параметра. Иными словами, говорить в приведенном случае об улучшении при таких то данных не приходится.

Рассмотрим похожий вариант. Все те же цифры, но в 10 раз больше.

> bad<-c(30,50)
> count<-c(1000,1020)
> prop.test(bad,count)

2-sample test for equality of proportions with continuity correction

data: bad out of count
X-squared = 4.3154, df = 1, p-value = 0.03777
alternative hypothesis: two.sided
95 percent confidence interval:
-0.036961241 -0.001077974
sample estimates:
prop 1 prop 2
0.03000000 0.04901961

Как видим, теперь гипотеза об одинаковости пропорций отвергается ( p-value = 0.03777 < 0.05) и налицо улучшение процесса.

Неверно построенные контрольные карты

Как известно, для контроля качества по альтернативному признаку используются следующие карты (см. ГОСТ Р ИСО 7870-2-2015. Статистические методы. Контрольные карты. Часть 2. Контрольные карты Шухарта)

C-карта. В таких контрольных картах строится график числа дефектов (в партии, в день, на один станок, в расчете на 100 футов трубы и т.п.). При использовании карты этого типа делается предположение, что дефекты контролируемой характеристики продукции встречаются сравнительно редко, при этом контрольные пределы для данного типа карт рассчитываются на основе свойств распределения Пуассона (распределения редких событий).
U-карта. В карте данного типа строится график относительной частоты дефектов, то есть отношения числа обнаруженных дефектов к n — числу проверенных единиц продукции (здесь n обозначает, например, число футов длины трубы, объем партии изделий). В отличие от C-карты, для построения карты данного типа не требуется постоянство числа единиц проверяемых изделий, поэтому ее можно использовать при анализе партий различного объема.
Np-карта. В контрольных картах этого типа строится график для числа дефектов (в партии, в день, на станок), как и в случае С-карты. Однако, контрольные пределы этой карты рассчитываются на основе биномиального распределения, а не распределения редких событий Пуассона. Поэтому данный тип карт должен использоваться в том случае, когда обнаружение дефекта не является редким событием (например, когда обнаружение дефекта происходит более чем у 5% проверенных единиц продукции). Этой картой можно воспользоваться, например, при контроле числа единиц продукции, имеющих небольшой брак.
P-карта. В картах данного типа строится график процента обнаруженных дефектных изделий (в расчете на партию, в день, на станок и т.д.). График строится так же, как и в случае U-карты. Однако контрольные пределы для данной карты находятся на основе биномиального распределения (для долей), а не распределения редких событий. Поэтому P-карта наиболее часто используется, когда появление дефекта нельзя считать редким событием (если, например, ожидается, что дефекты будут присутствовать в более чем 5% общего числа произведенных единиц продукции).

Поэтому построение U-карты для брака 30% в подавляющем большинстве случаев будет методически неверно. Далее, если у вас различный объем выборки, то, по уму, и контрольные границы будут различными. Рекомендую использовать пакет qcc в R, там все это делается проще.

Для построения контрольных карт используется команда qcc со следующим параметром type:

• «xbar» карта средних значений
• «R» R-карта размахов
• «S» S-карта стандартных отклонений
• «xbar.one» карта индивидуальных значений
• «p» карта пропорций (p-карта)
• «np» карта числа несоответствий (np-карта)
• «c» c-карта
• «u» u-карта
• «g» g-карта (карта редких событий)

Рассмотрим построение карты на примере имеющихся в пакете данных orangejuice

Концентрат замороженного апельсинового сока упакован в картонные банки размером 6 унций. Эти банки формуются на машине из картона с прикрепление металлической нижней панели. Затем проверяют, может ли емкость быть заполнена жидкость без протечек, как на боковом шве, так снизу в местах крепления. Если это происходит, банка считается несоответствующей. Данные собирали в виде 30 образцов по 50 банок каждый с интервалом в полчаса в течение трех смен, когда машина работала непрерывно. Для 15 использовался новый полуфабрикат из картона, который планируется вводить в производство. Образец 23 был получен, когда неопытный оператор был временно назначен на машину. После первых 30 образцов была произведена настройка машины. Затем еще 24 образца были взяты из рабочего процесса.

Данные содержат 54 наблюдения следующих четырех величин:

• sample — номер измерения;
• D — числ дефектов;
• size — размер выборки;
• trial -признак настройка (TRUE) или производство (FALSE)

> library(qcc)
> data(orangejuice)
> head(orangejuice)
sample D size trial
1 1 12 50 TRUE
2 2 15 50 TRUE
3 3 8 50 TRUE
4 4 10 50 TRUE
5 5 4 50 TRUE
6 6 7 50 TRUE
> qcc(orangejuice[orangejuice$trial,2],sizes=orangejuice[orangejuice$trial,3],type="p")

Аналогично можно построить контрольную карту со старыми и новыми значениями.

> qcc(orangejuice[orangejuice$trial,2], sizes=orangejuice[orangejuice$trial,3], newdata=orangejuice[!orangejuice$trial,2], newsizes=orangejuice[!orangejuice$trial,3],type="p")

SWOT анализ

Во всех работах без исключения SWOT-анализ делается просто для того, чтобы занять две-три страницы темы. Смысла в нем практического никакого, делается он неверно, сделанные на основании него выводы (в 80% случаев они просто отсутствуют) настолько тривиальны, что даже говорить о них не хочется. Рецензенты в большинстве случаев отмечают это в замечаниях. Если вы проектируете какую-нибудь инновацию это совершенно бесполезный инструмент, поскольку это инструмент стратегического планирования, а не DFSS. Не позорьтесь, выберите что-нибудь другое в качестве основы своего исследования.

Пусть игрок А обладает m возможными стратегиями А₁,…,А_m, а природа П может находиться в одном из n своих состояний П₁,…,П_n. Предполагается обычно, что игрок А в состоянии оценить результаты выбора им каждой из своих стратегий А_i, i=1,…,m, при каждом состоянии природы П_j, j=1,…,n, количественно выражающиеся действительными числами а_ij.

Эти числа, называемые выигрышами игрока А, можно записать в виде матрицы:

library(knitr)
mtx<-matrix(c(2,1,2,
            3,5,6,
            1,4,2,
            4,3,1.5),ncol=4)
rownames(mtx)<-paste0("A~",1:nrow(mtx),"~")
colnames(mtx)<-paste0("П~",1:ncol(mtx),"~")
kable(mtx)

(w<-c(0.5,0.1,0.2,0.2))

## [1] 0.5 0.1 0.2 0.2

sum(w)

## [1] 1

B<-mtx %*% w
colnames(B)<-"П~B~"
kable(B)

max(B)

## [1] 2.4

(w<-c(0.25,0.25,0.25,0.25))

## [1] 0.25 0.25 0.25 0.25

BL<-mtx %*% w
colnames(BL)<-"П~B~"
kable(BL)

max(BL)

## [1] 3.25

V<-as.matrix(apply(mtx,1,min))
colnames(V)<-"П~W~"
kable(V)

max(V)

## [1] 1.5

kable(HL<-cbind(B,V))

(w<-c(0.75,0.25))

## [1] 0.75 0.25

kable(HLM<-HL %*% w,col.names = "КХЛ")

max(HLM)

## [1] 2.1

M<-as.matrix(apply(mtx,1,max))
colnames(M)<-"П~M~"
kable(M)

max(M)

## [1] 6

kable(HV<-cbind(V,M))

(w<-c(0.5,0.5))

## [1] 0.5 0.5

kable(HVM<-HV %*% w,col.names = "КГ")

max(HVM)

## [1] 3.75

mtx1<- matrix(rep(apply(mtx,2,max),nrow(mtx)),byrow=T,ncol=ncol(mtx))-mtx
kable(mtx1)

S<-as.matrix(apply(mtx1,1,max))
colnames(S)<-"П~S~"
kable(S)

min(S)

## [1] 1

Для начала, откроем эту страницу в режиме просмотра кода. По идее, для динамического построения графиков может использоваться xml или json. Поиск в коде находит нам следующую строку:

<a class="link link_ajax link_theme_normal i-bem" data-bem='{"link":{}}' tabindex="0" href="/quotes/index.html">Ещё</a></div></div><div class="quote i-bem" data-bem='{"quote":{"signValDigit":4,"animation":true,"url":"/quotes/graph_1.json"}}'>

Интересует именно ссылка “/quotes/graph_1.json”. Вспоминаем HTML, понимаем, что это ссылка на данные от корня сайта, иными словами, она имеет вид https://news.yandex.ru/quotes/graph_1.json.

Теперь надо разобраться, как импортировать json в R. Ищем. Первая же ссылка дает ответ.

install.packages("rjson") 
library("rjson") 
json_file <- "https://news.yandex.ru/quotes/graph_1.json" 
json_data <- fromJSON(file=json_file)

# install.packages("rjson")
library("rjson")
json_file <- "https://news.yandex.ru//quotes/graph_1006.json" 
l <- fromJSON(file=json_file) 
df <- data.frame(matrix(unlist(l$prices), nrow= length(l$prices), byrow=T),stringsAsFactors=FALSE) 
tail(df) 

##                X1    X2
## 6302 1.481748e+12 53.79
## 6303 1.481835e+12 54.19
## 6304 1.481921e+12 55.26
## 6305 1.481927e+12 55.33
## 6306 1.482180e+12 54.76
## 6307 1.482235e+12 55.46

df$Dat=format(as.Date(df$X1/3600/24/1000, origin = "1970-01-01"),"%Y-%m-%d") 
tail(df)

##                X1    X2        Dat
## 6302 1.481748e+12 53.79 2016-12-14
## 6303 1.481835e+12 54.19 2016-12-15
## 6304 1.481921e+12 55.26 2016-12-16
## 6305 1.481927e+12 55.33 2016-12-16
## 6306 1.482180e+12 54.76 2016-12-19
## 6307 1.482235e+12 55.46 2016-12-20

GetYandexData<-function(idx,valuename="Value") {
  json_file <- paste0("https://news.yandex.ru//quotes/graph_",idx,".json")
  l <- fromJSON(file=json_file)
  df <- data.frame(matrix(unlist(l$prices), nrow= length(l$prices), byrow=T),stringsAsFactors=FALSE)
  df$Dat=format(as.Date(df$X1/3600/24/1000, origin = "1970-01-01"),"%Y-%m-%d")
  names(df)[2]<-valuename
  return(df[,c(3,2)]) 
}

z1<-GetYandexData(1006,"Brent")
z2<-GetYandexData(10,"Gold")
z3<-GetYandexData(1,"USD")

mytable<-merge(merge(z1,z2),z3)
tail(mytable)

##             Dat Brent   Gold     USD
## 5107 2016-12-16 55.26 1135.4 61.7515
## 5108 2016-12-16 55.26 1136.8 61.7515
## 5109 2016-12-16 55.33 1135.4 61.7515
## 5110 2016-12-16 55.33 1136.8 61.7515
## 5111 2016-12-19 54.76 1141.0 61.7931
## 5112 2016-12-20 55.46 1134.4 61.7967

graphData <- mytable
graphData$year=as.factor(format(as.Date(graphData$Dat),"%Y"))
graphData$month=as.factor(format(as.Date(graphData$Dat),"%m"))
graphData$wd=as.factor(format(as.Date(graphData$Dat),"%u"))
tail(graphData)

##             Dat Brent   Gold     USD year month wd
## 5107 2016-12-16 55.26 1135.4 61.7515 2016    12  5
## 5108 2016-12-16 55.26 1136.8 61.7515 2016    12  5
## 5109 2016-12-16 55.33 1135.4 61.7515 2016    12  5
## 5110 2016-12-16 55.33 1136.8 61.7515 2016    12  5
## 5111 2016-12-19 54.76 1141.0 61.7931 2016    12  1
## 5112 2016-12-20 55.46 1134.4 61.7967 2016    12  2

# "В лесу родилась ёлочка..."
# См. http://www.r-bloggers.com/geom_christmas_tree-a-new-geom-for-ggplot2-v2-0/
# Нижеследующий код нужен для построения елочек
library(ggplot2)
GeomChristmasTree <- ggproto("GeomChristmasTree", Geom,
 required_aes = c("x", "y"),
 default_aes = aes(shape = 19, colour = "black", 
 fill = "green4", size = 3,
 linetype = 1, alpha = 1,
 fontsize = 1),
 draw_key = draw_key_polygon,
 
 draw_panel = function(data, panel_scales, coord) {
 coords <- coord$transform(data, panel_scales)
 
 # each tree has 4*branches + 3 points
 if (length(coords$size) == 1) {
 tsize <- rep(pmax(1, round(coords$size)), length(coords$x))
 theight <- rep(pmax(0, round(coords$size)), length(coords$x))
 } else {
 tsize <- pmax(1, round(coords$size))
 theight <- pmax(0, coords$size)
 }
 
 # scale factors
 r01x <- diff(range(coords$x))/100
 r01y <- diff(range(coords$y))/100
 
 # coords
 longx <- unlist(lapply(seq_along(coords$x), function(i) {
 if (tsize[i] == 1) {
 dx <- -c(0.3, 0.3, 1.2, 0, -1.2, -0.3, -0.3)
 } else {
 dx <- -c(0.3, 0.3, rep(c(1.2,0.3), tsize[i]-1), 1.2, 0, -1.2, rep(c(-0.3,-1.2), tsize[i]-1), -0.3, -0.3)
 }
 r01x*dx + coords$x[i]
 }))
 longy <- unlist(lapply(seq_along(coords$y), function(i) {
 if (tsize[i] == 1) {
 dy <- c(-0.5, 0, 0, theight[i], 0, 0, -0.5)
 } else {
 dy <- c(-0.5, 0, 0, rep(1:(tsize[i]-1), each=2), theight[i], rep((tsize[i]-1):1, each=2), 0, 0, -0.5)
 }
 r01y*dy + coords$y[i]
 }))
 longid <- unlist(sapply(seq_along(coords$y), function(i) {
 rep(i, each=4*tsize[i]+3)
 }))
 
 grid::polygonGrob(
 longx, 
 longy,
 id = longid,
 gp = grid::gpar(col = coords[,"colour"],
 fill = coords[,"fill"],
 fontsize = 10)
 )
 }
)
geom_christmas_tree <- function(mapping = NULL, data = NULL, stat = "identity",
                                position = "identity", na.rm = FALSE, show.legend = NA, 
                                inherit.aes = TRUE, ...) {
  layer(
    geom = GeomChristmasTree, mapping = mapping,  data = data, stat = stat, 
    position = position, show.legend = show.legend, inherit.aes = inherit.aes,
    params = list(na.rm = na.rm, ...)
  )
}
# "В лесу она росла..."

# Строим график зависимости курса от цены на нефть


ggplot(subset(graphData,year=="2016"), aes(x=Brent, y=USD)) + 
  stat_density_2d(aes(color=month))+
    geom_christmas_tree(aes(fill=month,size=wd)) +
  labs(x="Цена на Нефть Brent (ICE.Brent), USD/баррель",y="Курс USD ЦБ РФ, руб.")+
  scale_fill_discrete(name = "Месяц")+
  scale_size_discrete(guide=FALSE)+
  scale_color_discrete(guide=FALSE)+
  theme_bw()+theme(text=element_text(size=14))

## Warning: Using size for a discrete variable is not advised.

# Сохраняем картинку. Размер в дюймах.
ggsave("Pic1.png",width=6,height=6,dpi=300)

1. Скачиваем и устанавливаем на флешку платформу PortableApp.
2. Скачиваем и устанавливаем R Portable.
3. Скачиваем и устанавливаем RStudio Portable.

Все необходимое установлено. Теперь настраиваем. Для этого делаем следующие действия.

1. Запускаем приложение Start  в корне флеш-накопителя. Видим в трее новый значок, если на него нажать открывается окно.

2. Запускаем RStudioPortable. Видим вот такое окно:

Необходимо указать путь к R. Для этого нажимаем кнопку Browse и на флеш-накопителе выбираем папку F:\PortableApps\R-Portable\App\R-Portable\bin (в моём случае это диск F, у вас может быть другой). Приобретает вот такой вид (поскольку у меня 64-битная ОС).

Выбираем ту, что больше нравится, я выбрал 64-битную и нажимаем OK. Видим вот такое окно.

Соглашаемся и получаем запущенную RStudio

Теперь необходимо установить пакеты расширения, которые нам пригодятся в работе. Для этого полностью скопируем команду

source(url("https://tushavin.ru/RStudio/installall.R"))

и вставим после символа приглашения к вводу (значок «>» в левом окне), нажимаем клавишу Enter и идем пить чай, ну, или делаем что-либо другое на компьютере, поскольку все прекрасно работает и в фоновом режиме. Некоторое время будут писаться разные красные и черненькие буковки, не обращаем внимание, значит функционирует. После того как все закончится снова появится символ приглашения ввода.

3. Переходим к проверке работы. Открываем материал по линейному программированию. Внизу видим раздел «Решение с пакетом lpSolve» и прямо под ним код. Копируем его полностью, вставляем после приглашения, нажимаем Enter. Вставляем еще две команды:

result
result$solution

Результат должен совпасть с описанным в заметке.

В предыдущей заметке мы говорили о непрерывных распределениях, однако при работе с качественными признаками мы имеем дело с дискретными распределениями. Как известно, в области качества достаточно часто используется биноминальное распределение (например, Np-карта Шухарта). Рассмотрим пример практической задачи на это распределение.

Постановка задачи

Имеет тест на 45 вопросов по 4 варианта ответа в каждом. Только один вариант верный.Тест разбит на 3 блока по 15 вопросов. За правильные ответы на вопросы из 1-го блока начисляется 1 балл, из 2-го блока – 2 балла, из 3-го блока – 3 балла, и соответственно, максимальный балл, который можно набрать равен 90. Определить вероятность сдачи теста, тыкая ответы наугад, при условии, что для успешной сдачи нужна набрать минимум 61 балл.

> dbinom(0,15,0.25) # Вероятность получить ноль правильных ответов
[1] 0.01336346
> pbinom(0,15,0.25) # Вероятность получить 0 и менее правильных ответов
[1] 0.01336346
> pbinom(0,15,0.25,lower.tail = F) # Вероятность получить более 0
[1] 0.9866365
> dbinom(1,15,0.25) # Вероятность получить ровно один правильный ответ
[1] 0.06681731
> pbinom(1,15,0.25) # Вероятность получить один и менее (0,1)
[1] 0.08018077
> dbinom(0,15,0.25)+dbinom(1,15,0.25) # Сумма вероятностей получить 0 и 1 ответ
[1] 0.08018077
> pbinom(1,15,0.25,lower.tail = F) # вероятность получить более 1 правильного ответа
[1] 0.9198192

# Инициируем переменную для ответа
P<-0
# В Z мы будем накапливать вероятности для исходов теста 0:90
Z<-data.frame(X=0:90,y=rep(0,91)) 

# Перебираем возможные варианты для каждой части теста
 for(i in 0:15) for (j in 0:15) for(k in 0:15) { 
   # считаем вероятность куммулятивно для каждой точки
   Z$y[i+2*j+3*k+1]<-Z$y[i+2*j+3*k+1]+dbinom(i,15,0.25)*dbinom(j,15,0.25)*dbinom(k,15,0.25) 
   # считаем куммулятивно вероятность ответить правильно наугад
   if(i+2*j+3*k>=61)   P<-P+dbinom(i,15,0.25)*dbinom(j,15,0.25)*dbinom(k,15,0.25) 
}
Z$cdf<-cumsum(Z$y)
cat("Ответ на задачу составляет: ",P)

## Ответ на задачу составляет:  1.73944e-08

oldpar<-par(mfrow=c(2,1))
plot(Z$y~Z$X,main="Функция вероятности f(x)",xlab="x",ylab="f(x)",pch=19,cex=0.8,col="blue")
points(Z$y[62:91]~Z$X[62:91],cex=0.9,col="red",pch=19)
text(70,0.02,expression(sum(f[i],i==61,90)==1.73944 %*% 10^{-8}))
plot(Z$cdf~Z$X,main="Функция распределения F(x)",xlab="x",ylab="f(x)",pch=19,cex=0.8,col="blue")
abline(h=Z$cdf[61],v=60,col="darkgrey",lty=2)
points(Z$cdf[61]~Z$X[61],cex=0.9,col="red",pch=19)
text(70,0.8,expression(1-F(60)==1.73944 %*% 10^{-8}))

par(oldpar)

sum(Z$y[62:91]) # sum(f(t),t = 61...90)

## [1] 1.73944e-08

1-Z$cdf[61]  # 1-F(60); 61, а не 60, потому как нумерация в массиве начинается с 1

## [1] 1.73944e-08