Анализ данных измерения массы грецких орехов
В результате случайно отобранных из мешка грецких орехов, были произведены измерения их массы на точных электронных весах. Необходимо провести анализ полученных данных на выбросы, и проверить подчиняются ли измеренные массы орешков закону нормального распределения, так же определить среднюю массу орешков и доверительный интервал с точностью 95 процентов.
Данные по измерению массы орешков (в граммах): 11.95, 10.95, 12.56, 9.16, 10.08, 10.72, 9.29, 6.15, 5.12, 11.16, 9.51, 10.73, 10.44, 10.10, 10.10, 11.35, 13.05, 12.09, 11.09, 7.86, 5.12, 6.41, 8.76, 10.40, 11.36, 6.07, 11.11, 9.29, 10.25, 11.72, 9.82, 8.07, 9.20, 9.65, 9.87, 12.39, 10.72, 11.84, 9.81, 10.10, 11.04, 10.08, 12.91, 10.30, 12.80, 9.40, 8.65, 10.21, 9.50, 10.97, 11.65, 8.27, 9.90, 9.18, 9.83, 10.41, 10.24, 11.26, 12.37, 9.73, 6.10, 5.84, 12.00, 8.70, 10.15, 10.62, 5.21, 11.34, 9.39, 9.96, 4.80, 10.97, 11.51, 12.44, 11.17, 12.43, 9.41, 10.78, 11.05, 10.85, 10.68, 9.02, 12.37, 10.43, 11.96, 10.93, 10.59, 10.76, 8.27, 9.09, 9.32, 9.36, 10.51, 11.70, 11.76, 10.72, 7.23, 9.34, 9.39, 11.53, 9.85, 10.68, 11.19, 11.29, 10.14, 12.88, 8.85, 10.88, 11.17, 13.30, 10.45, 10.13, 12.07, 10.57, 10.20, 9.74, 8.09, 10.62, 11.12, 13.98, 11.87, 13.48, 11.00, 5.86, 11.24, 10.14, 10.58, 11.04, 9.03, 9.67, 10.73, 10.21, 9.99, 10.98, 5.49, 10.28, 10.70, 13.00, 11.98, 10.04, 10.79, 9.65, 10.05, 13.12, 9.41, 4.32, 10.77, 11.40
library(sm)
library(car)
library(DescTools)
massa <- c(11.95, 10.95, 12.56, 9.16, 10.08, 10.72, 9.29, 6.15, 5.12,
11.16, 9.51, 10.73, 10.44, 10.10, 10.10, 11.35, 13.05,
12.09, 11.09, 7.86, 5.12, 6.41, 8.76, 10.40, 11.36, 6.07, 11.11,
9.29, 10.25, 11.72, 9.82, 8.07, 9.20, 9.65, 9.87, 12.39, 10.72,
11.84, 9.81, 10.10, 11.04, 10.08, 12.91, 10.30, 12.80, 9.40, 8.65,
10.21, 9.50, 10.97, 11.65, 8.27, 9.90, 9.18, 9.83, 10.41, 10.24,
11.26, 12.37, 9.73, 6.10, 5.84, 12.00, 8.70, 10.15, 10.62, 5.21,
11.34, 9.39, 9.96, 4.80, 10.97, 11.51, 12.44, 11.17, 12.43, 9.41,
10.78, 11.05, 10.85, 10.68, 9.02, 12.37, 10.43, 11.96, 10.93, 10.59,
10.76, 8.27, 9.09, 9.32, 9.36, 10.51, 11.70, 11.76, 10.72, 7.23,
9.34, 9.39, 11.53, 9.85, 10.68, 11.19, 11.29, 10.14, 12.88, 8.85,
10.88, 11.17, 13.30, 10.45, 10.13, 12.07, 10.57, 10.20, 9.74, 8.09,
10.62, 11.12, 13.98, 11.87, 13.48, 11.00, 5.86, 11.24, 10.14, 10.58,
11.04, 9.03, 9.67, 10.73, 10.21, 9.99, 10.98, 5.49, 10.28, 10.70,
13.00, 11.98, 10.04, 10.79, 9.65, 10.05, 13.12, 9.41,
4.32, 10.77, 11.40)
hist(massa)
grid(col='black')
boxplot(massa)
# создадим новый вектор данных без выбросов
# Выбросы хранятся тут: boxplot.stats(x)$out
# Заберём индексы точек выбросов
ind <- which(massa %in% boxplot.stats(massa)$out)
massa_outliers <- massa[-ind]
boxplot(massa_outliers)
# Полная описательная статистика с выбросами
summary(massa)
# Полная описательная статистика без выбросов
summary(massa_outliers)
# проверка на нормальное распределение
#Библиотека car предоставляет функцию
#qqPlot(...), которая по умолчанию добавляет
#точечную доверительную огибающую
#к обычному qq-графику:
qqPlot(massa_outliers)
shapiro.test(massa_outliers)
hist(massa_outliers,21, freq = F)
lines(density(massa_outliers), col='red', lwd=2)
grid(col='black')
print('доверительный интервал с выбросами')
MeanCI(massa, conf.level = 0.95)
print('доверительный интервал после удаления выбросов')
MeanCI(massa_outliers, conf.level = 0.95)
sm.density(massa_outliers, model = "Normal", xlab = "Имитированная выборка",
ylab = "Функция плотности распределения")
Первым делом постараемся ответить на первый вопрос задачи: подчинены ли массы орехов закону нормальному распределению?. Для этого просто построим гистограмму по данным: hist(massa). При рассмотрении гистограммы видно, что распределение похоже на нормальное, но все же в левой части наблюдаем множество столбиков: от 4 до 7, которые «размывают» купол нормального распределения. Можем предположить, что это могут быть выбросы. А следовательно, необходимо провести специальные расчеты на выявление и удаление выбросов (если они есть). В скрипте с помощью «ящика с усами» в строке boxplot(massa) мы наглядно видим, что у нас в наборе данных data действительно содержится большое количество данных с «выбросами». Эти данные можно вывести отдельно в консоль командой:
boxplot.stats(massa)$out и в результате мы получим данные «выбросов»: «[1] 6.15 5.12 5.12 6.41 6.07 6.10 5.84 5.21 4.80 13.98 5.86 5.49 4.32». Выбросы принимают минимальное значение от 4.3 и до 6.15 заканчивая максимальным значением 13.98. Поэтому далее скрипт автоматически удалит выбросы, и мы получим в итоге новый набор данных massa_outliers без выбросов. И вот теперь снова построим гистограмму данных без выбросов с помощью команды:
hist(massa_outliers,21, freq = F)
lines(density(massa_outliers), col='red', lwd=2)
grid(col='black')
В итоге мы получим новую гистограмму без выбросов. Если мы сравним первую (с выбросами) и вторую (без выбросов) гистограммы, то уже сможем увидеть картину уже более напоминающий купол нормального распределения. Это дает нам основание предположить что данные по массе грецких орехов подчиняются закону нормального распределения. Для полной уверенности в этом проведем QQ тест и shapiro.test(massa_outliers)
Shapiro-Wilk normality test data: massa_outliers W = 0.99291, p-value = 0.7374
Значение p-value = 0.7374 подтверждает нулевую гипотезу (предположение что данные взяты из генеральной совокупности подчиненной нормальному распределению) в высшей степени.
Осталось наконец ответить и на последний вопрос задачи: определить среднюю массу орешков и доверительный интервал с точностью 95 процентов. Это можно сделать сразу с помощью одной команды: MeanCI(massa_outliers, conf.level = 0.95). В результате ее выполнения мы получим в консоли
> MeanCI(massa_outliers, conf.level = 0.95) mean lwr.ci upr.ci 10.57296 10.36618 10.77975
Таким образом среднее значение массы орешков составило 10.57 грамм, и с точностью 95% процентов мы можем утверждать, что среднее значение массы орешков находится от 10.37 до 10.78 грамм).
Отличной альтернативой QQ графику, построенные с помощью функции qqPlot(), является практическое применение функции sm.densety из пакета sm. Установите параметр model = "Normal", и на графике получим доверительную полосу относительно теоретической функции плотности распределения. так же на графике будет отражена кривая ядерной плотности для эмпирических данных, благодаря чему можно оценить правдоподобность предположения о нормальности распределения и в каких диапазонах анализируемых данных наблюдаются отклонения.
sm.density(massa_outliers, model = "Normal", xlab = "Имитированная выборка",
ylab = "Функция плотности распределения")
