Описательная статистика и графический анализ

Описательная статистика и графика

Загрузка данных

Текстовый файл

#прочитаем данные об автомобилях
#разделитель - запятая
cars=read.csv("cars.csv")
str(cars)

## 'data.frame':    205 obs. of  17 variables:
##  $ make           : Factor w/ 22 levels "alfa-romero",..: 1 1 1 2 2 2 2 2 2 2 ...
##  $ fuel.type      : Factor w/ 2 levels "diesel","gas": 2 2 2 2 2 2 2 2 2 2 ...
##  $ num.of.doors   : Factor w/ 2 levels "four","two": 2 2 2 1 1 2 1 1 1 2 ...
##  $ body.style     : Factor w/ 5 levels "convertible",..: 1 1 3 4 4 4 4 5 4 3 ...
##  $ drive.wheels   : Factor w/ 3 levels "4wd","fwd","rwd": 3 3 3 2 1 2 2 2 2 1 ...
##  $ engine.location: Factor w/ 2 levels "front","rear": 1 1 1 1 1 1 1 1 1 1 ...
##  $ wheel.base     : num  88.6 88.6 94.5 99.8 99.4 ...
##  $ length         : num  169 169 171 177 177 ...
##  $ width          : num  64.1 64.1 65.5 66.2 66.4 66.3 71.4 71.4 71.4 67.9 ...
##  $ height         : num  48.8 48.8 52.4 54.3 54.3 53.1 55.7 55.7 55.9 52 ...
##  $ curb.weight    : int  2548 2548 2823 2337 2824 2507 2844 2954 3086 3053 ...
##  $ engine.size    : int  130 130 152 109 136 136 136 136 131 131 ...
##  $ horsepower     : int  111 111 154 102 115 110 110 110 140 160 ...
##  $ peak.rpm       : int  5000 5000 5000 5500 5500 5500 5500 5500 5500 5500 ...
##  $ city.mpg       : int  21 21 19 24 18 19 19 19 17 16 ...
##  $ highway.mpg    : int  27 27 26 30 22 25 25 25 20 22 ...
##  $ price          : int  13495 16500 16500 13950 17450 15250 17710 18920 23875 NA ...

#запись в файл
write.csv(cars,"cars.txt")

#Прочитаем данные о сотрудниках из текстового файла
#разделитель - табуляция
#кодировка указывается для чтения кириллицы
staff=read.table("staff.txt",sep="\t",encoding = "UTF-8")
str(staff)

## 'data.frame':    7 obs. of  6 variables:
##  $ name  : Factor w/ 7 levels "Александр","Алексей",..: 6 2 1 4 7 5 3
##  $ salary: int  500 600 400 650 300 550 2000
##  $ sex   : Factor w/ 2 levels "female","male": 2 2 2 1 1 2 2
##  $ weight: int  82 87 93 NA NA 72 82
##  $ height: num  1.75 1.8 1.9 1.65 1.7 1.72 1.85
##  $ size  : Factor w/ 5 levels "L","M","S","XL",..: 5 3 4 3 3 2 1

#запиcь в файл
write.table(staff,"staff.txt",sep="\t",fileEncoding="UTF-8")

Встроенные данные

#обращаться по имени данных
#для данных из базовой библиотеки datasets
str(iris)

## 'data.frame':    150 obs. of  5 variables:
##  $ Sepal.Length: num  5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
##  $ Sepal.Width : num  3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
##  $ Petal.Length: num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
##  $ Petal.Width : num  0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
##  $ Species     : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...

Управление таблицей данных

#выведем несколько первых записей (6 - по умолчанию)
head(cars)

##          make fuel.type num.of.doors  body.style drive.wheels
## 1 alfa-romero       gas          two convertible          rwd
## 2 alfa-romero       gas          two convertible          rwd
## 3 alfa-romero       gas          two   hatchback          rwd
## 4        audi       gas         four       sedan          fwd
## 5        audi       gas         four       sedan          4wd
## 6        audi       gas          two       sedan          fwd
##   engine.location wheel.base length width height curb.weight engine.size
## 1           front       88.6  168.8  64.1   48.8        2548         130
## 2           front       88.6  168.8  64.1   48.8        2548         130
## 3           front       94.5  171.2  65.5   52.4        2823         152
## 4           front       99.8  176.6  66.2   54.3        2337         109
## 5           front       99.4  176.6  66.4   54.3        2824         136
## 6           front       99.8  177.3  66.3   53.1        2507         136
##   horsepower peak.rpm city.mpg highway.mpg price
## 1        111     5000       21          27 13495
## 2        111     5000       21          27 16500
## 3        154     5000       19          26 16500
## 4        102     5500       24          30 13950
## 5        115     5500       18          22 17450
## 6        110     5500       19          25 15250

#2 последние записи
tail(cars,2)

##      make fuel.type num.of.doors body.style drive.wheels engine.location
## 204 volvo    diesel         four      sedan          rwd           front
## 205 volvo       gas         four      sedan          rwd           front
##     wheel.base length width height curb.weight engine.size horsepower
## 204      109.1  188.8  68.9   55.5        3217         145        106
## 205      109.1  188.8  68.9   55.5        3062         141        114
##     peak.rpm city.mpg highway.mpg price
## 204     4800       26          27 22470
## 205     5400       19          25 22625

#выведем первые 5 записей по переменной "make" (марка авто)
head(cars$make,5)

## [1] alfa-romero alfa-romero alfa-romero audi        audi       
## 22 Levels: alfa-romero audi bmw chevrolet dodge honda isuzu ... volvo

#подсчитаем средний расход топлива
#with позволяют напрямую обращаться к переменным
#для этого первым аргументом должна быть соотв. таблица
mpg=with(cars,(highway.mpg+city.mpg)/2)

#для непосредственного доступа к переменным
#присоединяем таблицу данных к текщему окружению
attach(cars)

#теперь обращаемся к любой переменной
#выведем первые цены для первых 5 моделей в списке
price[1:5]

## [1] 13495 16500 16500 13950 17450

#по окончании работы с данными не забыть отсоединить таблицу
#присоединение другой таблицы с совпадающими именами переменных:
#новые переменные "перекрывают" предыдущие
detach(cars)

Исключение пропущенных значений

#исключить строки, которые содержат пропущенное значение
#хотя бы по одной переменной
cars.noNA=na.omit(cars)
#сколько наблюдений исключено
nrow(cars)-nrow(cars.noNA)

## [1] 8

#удалим пропущенные значения у переменной
price.noNA=with(cars,price[!is.na(price)])
#сколько пропущенных значений
length(cars$price)-length(price.noNA)

## [1] 4

Статистики

#присоединим данные к текущему окружению
#для непосредственного доступа к переменным
attach(cars)

Простейшие статистики

#минимальный объем двигателя
min(engine.size)

## [1] 61

#максимальный объем двигателя
max(engine.size)

## [1] 326

#диапазон
range(engine.size)

## [1]  61 326

#ранжировать
 rank(c(5,3,2,7,11))

## [1] 3 2 1 4 5

#cумма цен всех моделей
sum(price,na.rm=T)

## [1] 2654633

#кумулятивная сумма: удобно для вероятностей
cumsum(c(0.1,0.2,0.3,0.4))

## [1] 0.1 0.3 0.6 1.0

Характеристики центра

#среднее: средний объем двигателя (куб. дюймов) по всем моделям
mean(engine.size)

## [1] 126.9

#средняя мощность (л.с.)
mean(horsepower) #есть пропущенные значения

## [1] NA

mean(horsepower,na.rm=TRUE) #проп.зн. не участвуют в вычислениях

## [1] 104.3

#медиана: медиана цен
median(engine.size)

## [1] 120

#медиана - частный случай квантиля
quantile(engine.size,probs=0.5)

## 50% 
## 120

Характеристики разброса значений

#стандартное отклонение
sd(engine.size)

## [1] 41.64

#дисперсия
var(engine.size)

## [1] 1734

sd(engine.size)^2 #дисперсия

## [1] 1734

#вычисление дисперсии с помощью векторного выражения
#несмещенная оценка
n=length(engine.size)
x=engine.size
sum((x-mean(x))^2/(n-1))

## [1] 1734

Асиммерия и эксцесс

Коэффициент асимметрии (Sk): мера симметричности распределения СВ.

1. Sk=0: расспределение симметрично относительно математического ожидания (МО).
   
2. Sk>0: распределение скошено влево (левый хвост распределения “тяжелее”).
   
3. Sk<0: распределение скошено вправо.

Коэффициент эксцесса (Ex): мера остроты пика распределения СВ.

1. Ex=0: острота пика нормального распределения.
   
2. Ex>0: пик острее и хвосты “легче”, чем у нормального распределения.
   
3. Ex<0: пик более пологий, а хвосты “тяжелее”, чем у норм. распр.

library(e1071) #для доступа к следующим 2 функциям
skewness(x) #асимметрия: распределение скошено влево

## [1] 1.919

kurtosis(x) #эксцесс: острый пик

## [1] 5.069

Другие характеристики

#квантили
quantile(price.noNA,probs=c(0.25,0.5,0.75))

##   25%   50%   75% 
##  7775 10295 16500

#интерквартильный размах
IQR(engine.size)

## [1] 44

#другой способ
q1=quantile(engine.size,0.25)
q3=quantile(engine.size,0.75)
q3-q1

## 75% 
##  44

Комплексные статистики

#min,q1,median,mean,q3,max,number of NA's
#для одной переменной
summary(price)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    5120    7780   10300   13200   16500   45400       4

#для всех переменных таблицы данных
#для категоральных переменных - только частоты
summary(cars)

##          make      fuel.type   num.of.doors       body.style drive.wheels
##  toyota    : 32   diesel: 20   four:114     convertible: 6   4wd:  9     
##  nissan    : 18   gas   :185   two : 89     hardtop    : 8   fwd:120     
##  mazda     : 17                NA's:  2     hatchback  :70   rwd: 76     
##  honda     : 13                             sedan      :96               
##  mitsubishi: 13                             wagon      :25               
##  subaru    : 12                                                          
##  (Other)   :100                                                          
##  engine.location   wheel.base        length        width     
##  front:202       Min.   : 86.6   Min.   :141   Min.   :60.3  
##  rear :  3       1st Qu.: 94.5   1st Qu.:166   1st Qu.:64.1  
##                  Median : 97.0   Median :173   Median :65.5  
##                  Mean   : 98.8   Mean   :174   Mean   :65.9  
##                  3rd Qu.:102.4   3rd Qu.:183   3rd Qu.:66.9  
##                  Max.   :120.9   Max.   :208   Max.   :72.3  
##                                                              
##      height      curb.weight    engine.size    horsepower     peak.rpm   
##  Min.   :47.8   Min.   :1488   Min.   : 61   Min.   : 48   Min.   :4150  
##  1st Qu.:52.0   1st Qu.:2145   1st Qu.: 97   1st Qu.: 70   1st Qu.:4800  
##  Median :54.1   Median :2414   Median :120   Median : 95   Median :5200  
##  Mean   :53.7   Mean   :2556   Mean   :127   Mean   :104   Mean   :5125  
##  3rd Qu.:55.5   3rd Qu.:2935   3rd Qu.:141   3rd Qu.:116   3rd Qu.:5500  
##  Max.   :59.8   Max.   :4066   Max.   :326   Max.   :288   Max.   :6600  
##                                              NA's   :2     NA's   :2     
##     city.mpg     highway.mpg       price      
##  Min.   :13.0   Min.   :16.0   Min.   : 5118  
##  1st Qu.:19.0   1st Qu.:25.0   1st Qu.: 7775  
##  Median :24.0   Median :30.0   Median :10295  
##  Mean   :25.2   Mean   :30.8   Mean   :13207  
##  3rd Qu.:30.0   3rd Qu.:34.0   3rd Qu.:16500  
##  Max.   :49.0   Max.   :54.0   Max.   :45400  
##                                NA's   :4

#статистики для построения "ящика с усами"
fivenum(engine.size)

## [1]  61  97 120 141 326

Другие статистики

#подсчитать частоты значений для категоральной переменной
#число моделей авто с дизельными двигателем (diesel) и бензиновым (gas)
table(fuel.type)

## fuel.type
## diesel    gas 
##     20    185

Графики

Базовая графическая система graphics

Функция plot

Примеры использования

#одна переменная
plot(engine.size) # для количественной переменной

plot(fuel.type) #для категоральной (фактора)

#две переменные
plot(engine.size,horsepower)

#использование интерфейса "Формула" (y~x: y зависит от x)
plot(horsepower~engine.size)

#заголок (main) и подписи осей (xlab,ylab)
#при явном указании имен параметров порядок не важен
plot(engine.size,horsepower,main="Horsepower vs. Engine size",
     ylab="Horsepower",xlab="Engine size")

#цвет точек (наблюдений) по типу двигателя и легенда
colors=ifelse(fuel.type=="gas","blue","green")
plot(horsepower~engine.size,col=colors)
#легенда накладывается последовательным вызовом данной функции:
legend("top",legend=c("gas","diesel"),fill=c("blue","green"))

добавление элементов на график

plot(engine.size[1:10],ylab=expression(eta)) #expression: мат. символы
title("Figure") # добавить заголовок
lines(engine.size[1:10]) #соединить линиями
points(c(140,150)) #добавить точки
curve(15*x+80,from=2,to=6,add = T) #добавить линию

Графические параметры

#изменим вид точки
par(pch = 2)
plot(engine.size[1:10])

#зададим количество отображаемых графиков: одна строка, два столбца
par(mfrow=c(1,2))
plot(horsepower~engine.size,main="Plot")
hist(horsepower,main="Histogram")

par(mfrow=c(1,1),pch=1) #вернуть по умолчанию

Другие графики

Больше примеров в [Quick-R](http://www.statmethods.net)

Гистограмма

#распределение объемов двигателя для авто из выборки
hist(engine.size) #по оси y - частоты

# по оси y - относительные частоты
hist(engine.size,probability = T) #~ freq=F

Столбиковая диаграмма

(t=table(drive.wheels)) # подсчитаем частоты моделей по типу ведущих колес

## drive.wheels
## 4wd fwd rwd 
##   9 120  76

barplot(t) #выведем диаграмму

(t=table(fuel.type,drive.wheels))

##          drive.wheels
## fuel.type 4wd fwd rwd
##    diesel   0   9  11
##    gas      9 111  65

barplot(t,legend.text = c("diesel","gas"))

“Ящик с усами”

#кружочками показаны аномальные (резко выделяющиеся) значения
boxplot(engine.size,ylab="engine size")

#в разрезе по типу двигателя#аномальные наблюдения только среди бензиновых двигателей
boxplot(engine.size~fuel.type,ylab="engine size",xlab="fuel type")

#по нескольким переменным
boxplot(city.mpg,highway.mpg,ylab="miles per gallon",names=c("city","highway"))

График ядерной оценки плотности

#плотность
d=density(engine.size)
plot(d,main="Распределение для объема двигателя",ylab="Плотность распределения")

График “Квантиль-квантиль”

qqnorm(engine.size)

Диаграмма рассеивания

pairs(iris[1:4])

Графическая система lattice (решетки)

library(lattice) #загрузка графической библиотеки
#построим  гистограммы
histogram(horsepower ~ engine.size |fuel.type)

#загрузим библиотеку MASS
library(MASS)
#используем данные Cars93 из библиотеки MASS
#ящики с усами для максимальной цены в разрезе по количеству цилиндров
#только одно авто имеет роторный двигатель
bwplot( ~ Max.Price | Cylinders , data = Cars93)

#диаграмма рассеяния с добавлением линии регрессии
#зависимость между переменными "Расход топлива на магистрали" и "Объем бака" в разрезе по типу автомобиля
#зададим функцию, строящую график в каждой ячейке
plot.regression = function(x,y) {
 panel.xyplot(x,y)
 panel.abline(lm(y~x))
}

xyplot(MPG.highway ~ Fuel.tank.capacity | Type, 
       data=Cars93,panel = plot.regression)

Использование шкалы для графиков

#загрузим библиотеку
library(MASS)
#построим двумерный график по данным о млекопитающих (mammals)
plot(mammals) #наблюдения

#логарифмируем переменные
attach(mammals)
plot(log(brain)~log(body))

detach(mammals)

Моделирование

Cлучайные величины

Равномерное распределение

Функции

1. runif - генерация случайных величин (СВ)
2. punif - функция распределения вероятностей (ФРВ)
3. dunif - функция плотности РВ
4. qunif - функция вычисления квантилей (обратная ФРВ)

#пример генерации СВ
runif(10) #10 БСВ: на интервале [0,1]

##  [1] 0.08365 0.99914 0.83307 0.13365 0.44708 0.61104 0.22937 0.97055
##  [9] 0.19250 0.44670

runif(10,1,10) #на интервале [1,10]

##  [1] 7.946 7.221 9.431 9.269 8.636 5.207 1.134 1.130 7.409 2.424

#построим гистограмму
x=runif(100)
hist(x,probability=T,col=gray(0.9),main="Равномерное на [0,1]")
curve(dunif(x,0,1),add=T) #добавить теоретическую кривую

(x=runif(5))

## [1] 0.1318 0.7054 0.5508 0.2960 0.7814

#вычислим вероятности распределения x
punif(x)

## [1] 0.1318 0.7054 0.5508 0.2960 0.7814

#вычислим плотность распределения
dunif(x)

## [1] 1 1 1 1 1

Биномиальное распределение

#10 бросаний монеты
rbinom(n=10,size=1,p=0.5)

##  [1] 1 0 1 1 1 0 1 1 0 1

#n=10 партий продукции по size=100 изделий с вероятностью брака p=0.05
#количество бракованных изделий в каждой партии
rbinom(10,100,0.05)

##  [1]  2 10  2  1  6  7  2  5  5  5

#построим гистограммую для n=100,size=50,p=0.25
x=rbinom(100,100,0.25)
hist(x,probability = T)
#наложим теоретическое распределение
points(0:100,dbinom(0:100,100,0.25),type="h",lwd=3)

Нормальное распределение

Теоретическое

#стандартное
x=rnorm(100)
hist(x,freq=F)
curve(dnorm,-4,4,add=T)

#для большего числа наблюдений
hist(rnorm(1000),freq=F)
curve(dnorm,-4,4,add=T)

#характеристики асимметрии и эксцеса: близки к нулю
x=rnorm(1000)
skewness(x)

## [1] 0.07978

kurtosis(x)

## [1] -0.09657

#измxенение среднего и дисперсии
x=seq(from=-10,to=10,by=0.1)
curve(dnorm(x,-3,0.5),from=-10,to=10,col="green",
      ylab="Density",)
curve(dnorm(x),from=-10,to=10,col="red",add=T)
curve(dnorm(x,2,3),from=-10,to=10,col="blue",add=T)
legend("topright",c("N(0,1)","N(2,3)","N(-3,0.5)"),fill=c("red","blue","green"))

#использование функций
pnorm(2) #вероятность что СВ x< 2

## [1] 0.9772

pnorm(2,lower.tail = F) #вероятность, что СВ x>2

## [1] 0.02275

qnorm(0.95) #квантиль уровня 0.95: обратная функция для pnorm

## [1] 1.645

Пример нормального распределения

`Файл [pop1.csv](http://pluto.huji.ac.il/~msby/StatThink/Datasets/pop1.csv) содержит данные о 100000 человек.

pop1=read.csv("pop1.csv")
#построим гистограмму по переменной "height"
hist(pop1$height,freq=F) #нормальное распределение
#параметры
m=mean(pop1$height) #среднее
s=sd(pop1$height) #стандартное отклонение
x=seq(from=min(pop1$height),to=max(pop1$height),by=1)
curve(dnorm(x,m,s),add=T) #нанесем теоретическую кривую

Cлучайная подвыборка: sample

#случайный выбор без повторений: игра в лотерею
sample(1:100,5) #выбрать 5 чисел из 100

## [1] 34 35  8 44 55

sample(1:100,5) #снова выберем 5 чисел: другой результат#случайный выбор #с повторениями: игра в кости

## [1] 24 81 21 33 37

sample(1:6,10,replace=TRUE) #кинем кости 10 раз

##  [1] 3 2 2 4 2 3 5 1 1 4

sample(1:6,10,replace=T) #еще одна серия экспериментов

##  [1] 5 2 1 1 5 2 2 4 1 6

Тесты согласия: проверка распределения

Критерий Хи-квадрат Пирсона

#первый 2000 цифр после запятой для числа пи
pi2000=scan("pi2000.txt")
barplot(table(pi2000))

#проверим с помощью теста гипотезу, 
#что цифры попадаются с одинаковой вероятностью 0.1
chisq.test(table(pi2000)) #гипотеза не отклоняется

## 
##  Chi-squared test for given probabilities
## 
## data:  table(pi2000)
## X-squared = 4.42, df = 9, p-value = 0.8817

Критерий Колмогорова-Смирнова

Модельные данные

#проверим на нормальное распределение
#p-value > 0.05 : гипотеза о нормальном распределении (pnorm - ФРВ норм. распр.)
#подтверждается на уровне 0.05
ks.test(rnorm(1000),pnorm) # гипотеза о нормальном распр. отклоняется

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  rnorm(1000)
## D = 0.027, p-value = 0.4598
## alternative hypothesis: two-sided

ks.test(runif(1000),pnorm)

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  runif(1000)
## D = 0.5004, p-value < 2.2e-16
## alternative hypothesis: two-sided

Реальные данные

#проверим распределение роста на нормальное на выборке из 100 000 человек
#извлечем только неповторяющиеся значения
heights=unique(pop1$height)
#сколько уникальных значений роста
length(heights)

## [1] 94

m=mean(heights)
s=sd(heights)
#указываем выборочные параметры распределения:
ks.test(heights,pnorm,mean=m,sd=s)

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  heights
## D = 0.0605, p-value = 0.8607
## alternative hypothesis: two-sided

Критерий Лиллиефорса (Lilliefors)

Модифицированный критерий Колмогорова-Смирнова

library(nortest)
lillie.test(heights)

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  heights
## D = 0.0605, p-value = 0.5409