Задание 1. Моделирование БСВ
- Примеры
  - Сравнение с распределением цифр числа пи
  - n=1000 БСВ
- Задание 1
  - Варианты
Задание 2. Моделирование дискретных случайных величин
- Примеры
  - Биномиальное распределение
  - Распределение Пуассона
- Задание 2
  - Варианты
Задание 3. Моделирование цепи Маркова
- Пример
- Задание 3
  - Варианты
Задание 4. Моделирование непрерывных случайных величин
- Примеры
  - Нормальное распределение
  - Распределение Хи-квадрат пирсона
- Задание 4
  - Варианты
Задание 5. Статистический анализ данных

Задание 1. Моделирование БСВ

Примеры

Графики распределений

#для отображения двух графиков одновременно
par(mfrow=c(1,2))
#плотность распределения
curve(dunif,0,1,ylab="P(x)",xlim=c(0,1))
title("ПРВ")
#функция распределения
curve(punif,0,1,ylab="F(x)",xlim=c(0,1))
title("ФРВ")

plot of chunk unnamed-chunk-1

#вернуть к отображению одного графика
par(mfrow=c(1,1))

Численные характеристики

#среднее (или МО - математическое ожидание)
(1+0)/2

## [1] 0.5

#дисперсия
(1-0)^2/12

## [1] 0.08333

Сравнение с распределением цифр числа пи

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

## pi2000
##   0   1   2   3   4   5   6   7   8   9 
## 181 213 207 189 195 205 200 197 202 211

barplot(pi2000.digits/2000,col=gray(0.9),main="2000 цифр числа Пи")

plot of chunk unnamed-chunk-3

#равномерное на [1,10]
y=runif(2000,min=0,max=10)
y.cut=cut(y,breaks=seq(from=0,to=10,by=1))
table(y.cut)

## y.cut
##  (0,1]  (1,2]  (2,3]  (3,4]  (4,5]  (5,6]  (6,7]  (7,8]  (8,9] (9,10] 
##    231    187    198    207    191    191    189    206    206    194

barplot(table(y.cut)/2000, col=gray(0.9),main="Равномерное на [0,10]")

plot of chunk unnamed-chunk-3

n=1000 БСВ

par(mfrow=c(1,2))
x=runif(1000)
plot(density(x),xlab="x",ylab="Pn(x)",main="ЭПРВ")
plot(ecdf(x),main="ЭФРВ")

plot of chunk unnamed-chunk-4

#переключиться в режим одного графика
par(mfrow=c(1,1))

Задание 1

1.Написать функцию моделирования БСВ с помощью мультипликативного конгруэнтного метода (МКМ) c параметрами a0,b,M.

Спецификация функции МКМ-датчика:  

n - размерность вектора БСВ  
runif.mkm <- function(n,a0,b,M=2^31){  
    
}

2.Реализовать функцию моделирования БСВ с помощью метода Маклорена-Марсальи (ММ)
c параметром K (размерность таблицы) и с использованием двух МКМ-датчиков из пункта 1 (согласно варианту). использовать функцию runif.mkm.

Спецификация функции ММ-датчика:  

runif.mm <- function(n,K,a01,b1,a02,b2,M=2^31){  
  
}  

Для параметров K,a01,b1,a02,b2  установить значения по умолчанию согласно варианту.

3.Смоделировать выборку из n=1000 реализаций БСВ с помощью функции по методу Маклорена-Марсальи и на основе данной выборки

построить графики

а) ЭПРВ (эмпирическая плотность распределения вероятностей)

б) ЭФРВ (эмпирическая функция распределения вероятностей)

вычислить статистики

а) m - выборочное среднее

б) d - выборочная дисперсия

Сравнить графики и статистики, построеннные по выборке,с теоретическими.

Варианты

Вариант	a01	b1	a02	b2	K
1	1025	1025	715828205	715828205	32
2	71583743	71583743	787410923	787410923	64
3	143166461	143166461	2^4 + 3	858993641	128
4	214749179	214749179	5^3	930576359	256
5	2^28 + 3	286331897	5^11	1002159077	32
6	5^7	357914615	2^8 + 3	1073741795	64
7	5^5	429497333	1025	1025	128
8	2^12 + 3	501080051	71583743	71583743	256
9	572662769	572662769	143166461	143166461	32
10	644245487	644245487	214749179	214749179	64
11	715828205	715828205	2^28 + 3	286331897	128
12	787410923	787410923	5^7	357914615	192
13	2^4 + 3	858993641	5^5	429497333	256
14	5^3	930576359	2^12 + 3	501080051	32

Задание 2. Моделирование дискретных случайных величин

Примеры

#установим начальное значение для генератора псевдослучайных чисел
#для воспроизводимости результатов
set.seed(9929)

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

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

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

#количество бракованных изделий в партии из size=100 изделий
#mean=size*prob=5: в среднем попадается 5 бракованных
rbinom(1,100,0.05)

## [1] 0

#n=10000 (количество партий),size=10,prob=0.5
x=rbinom(10000,10,0.5)

#столбиковая диаграмма
barplot(table(x))

plot of chunk unnamed-chunk-6

#mean=size*prob=5 - теорет. матем. ожидание (среднее)
mean(x)

## [1] 5.002

#var=size*prob*(1-prob)=2.5 - теорет. дисперсия
var(x)

## [1] 2.486

t=as.vector(table(x)) #частота каждой слуайной величины (СВ)
m=length(t)-1 #СВ - целые числа в диапазоне от 0 до m
p=dbinom(0:m, size=10,prob=0.5)
#тест Хи-квадрат - проверка на распределение
#p-value > 0.05: на уровне значимости 0.05 
#гипотеза о бином. распр. не отклоняется
chisq.test(t,p=p)

## 
##  Chi-squared test for given probabilities
## 
## data:  t
## X-squared = 7.252, df = 10, p-value = 0.7015

Распределение Пуассона

#пример: количество пользователей на форуме за день
#среднее число пользователей lambda = 5
rpois(n=1,lambda = 5)

## [1] 5

#n=10000 (данные за 10000 дней)
x=rpois(n=10000,lambda = 5)
#максимальное число пользователей не ограничено
#для бином. распр. число столбиков ограничено параметрами
barplot(table(x))

plot of chunk unnamed-chunk-7

#mean=lambda=5 - теоретическое значение
mean(x)

## [1] 4.963

#var=lambda=5
var(x)

## [1] 4.869

t=as.vector(table(x)) #частоты каждой случайной величины (СВ)
m=length(t)-1 #СВ - целые числа в диапазоне от 0 до m
p=dpois(0:m, lambda = 5);
p=p/sum(p) #для того, чтобы sum(p)=1
#тест Хи-квадрат
#simulate.p.value = TRUE: для вычисления p.value методом Монте-Карло
#p-value > 0.05: на уровне значимости 0.05 
#гипотеза о пуассоновском распр. не отклоняется
chisq.test(t,p=p,simulate.p.value = TRUE, B=1000)

## 
##  Chi-squared test for given probabilities with simulated p-value
##  (based on 1000 replicates)
## 
## data:  t
## X-squared = 12.3, df = NA, p-value = 0.7103

Задание 2

Написать функцию моделированию ДСВ с распределением согласно варианту.

Спецификации функций:

rbinom2=function(n,size=10,prob=0.5){

}

rpois2=function(n,lambda=5){

}

Замечания

1. В качестве БСВ использовать функцию runif.
2. Для моделирование вектора СВ вместо цикла использовать функцию sapply. Код для моделирования СВ оформить в виде анонимной функции:  
vectorOfRandNumb=sapply(1:n,function(x){#алгоритм моделирования})  
3. Алгоритмы моделирования СВ взять из следующей книги:

Вадзинский, Р.Н. Справочник по вероятностным распределениям.

Смоделировать n=10000 реализаций ДСВ (c пом. функции rbinom2 или rpois2 - согласно варианту).

Для воспроизводимости результатов установить начальное значение генератора СВ (для БСВ). Использовать функции set.seed(). В качестве аргумента использовать простое число.

Провести следующий анализ:

вычислить среднее, дисперсию и сравнить с теоретическими;
построить столбиковую диаграмму;
с помощью функции chisq.test проверить точность моделирования

Варианты

Вариант	Распределение	Параметры	Номер алгоритма
1	Bi	size=4,prob=0.25	1
2	Bi	size=4,prob=0.25	2
3	Bi	size=4,prob=0.25	3
4	Bi	size=10,prob=0.5	1
5	Bi	size=10,prob=0.5	2
6	Bi	size=10,prob=0.5	3
7	Bi	size=5,prob=0.4	1
8	Bi	size=5,prob=0.4	2
9	Bi	size=5,prob=0.4	3
10	P	lambda=5	1
11	P	lambda=5	2
12	P	lambda=2	1
13	P	lambda=2	2
14	P	lambda=1	1

Замечания

Bi - биномиальное распределение, P - Пуассона.
Номер алгоритма: для соответствующего распределения( P - c.45, Bi - c.51)
В книге для параметров распределений следующее соответствие:
size - параметр n, prob - p, lambda = \(\mu\)
При моделировании пуассоновских СВ с помощью алгоритма 2 произведение a формировать из разных реализаций БСВ (вызывать функцию runif(1) на каждой итерации цикла: в книге - ошибка!)

Задание 3. Моделирование цепи Маркова

Теорию по цепи Маркова можно посмотреть здесь.

Пример

Моделирование погоды

Пусть в некотором городе в году 100 солнечных дней, 155 - облачных, 110 - дождливых.

#Зададим параметры
#распределение вероятностей по состояниям в начальный момент времени t=0
#относительные частоты разных состояний погоды в течение года
p0=c(100,155,110)/365
state.names=c("shiny","cloudy","rainy")
names(p0)=state.names

p0

##  shiny cloudy  rainy 
## 0.2740 0.4247 0.3014

#матрица переходных вероятностей
#вероятности смены состояний погоды (оценивается за год)
#byrow = TRUE: матрица заполняется числами построчно
P=matrix(c(0.6,0.4,0,0.1,0.6,0.3,0.1,0.3,0.6),ncol=3,byrow = TRUE)
colnames(P)=state.names
rownames(P)=state.names

P

##        shiny cloudy rainy
## shiny    0.6    0.4   0.0
## cloudy   0.1    0.6   0.3
## rainy    0.1    0.3   0.6

#Объединим параметры в одну структуру данных - список
Params=list(p0=p0,P=P)

#теперь можно обращаться к элементам списка
Params$p0

##  shiny cloudy  rainy 
## 0.2740 0.4247 0.3014

Моделирование однородной цепи Маркова (ОЦМ)

#загрузим функцию моделирования из файла
source("markovchain.R")

#смоделируем одну цепь Маркова длины 30, Params - список параметров
mc=markovchain(365,Params)
#результат - фактор, спец. структура данных для катег. перем.
#ОЦМ - ДСВ: посл-ть 
str(mc)

##  Factor w/ 3 levels "shiny","cloudy",..: 2 3 3 2 2 1 1 2 3 3 ...

#10 первых наблюдений (выводятся метки: 1="shiny",2="cloudy",3="rainy")
mc[1:10]

##  [1] cloudy rainy  rainy  cloudy cloudy shiny  shiny  cloudy rainy  rainy 
## Levels: shiny cloudy rainy

#найдем частоту каждого состояния погоды
table(mc)

## mc
##  shiny cloudy  rainy 
##     65    159    141

Стационарное распределение ОЦМ

#получим 1000 реализаций
#lapply повторяет выполнение функции markovchain
#возвращает список (l=list)
mc.list=lapply(rep(365,1000),markovchain,Params=Params)

#извлечем последнее 365-ое наблюдение по каждой реализации
#для этого ипользуем анонимную функцию
#sapply перебирает индексы и из соотв. цепи Маркова извлекает посл. знач.
#sapply возвр. вектор (от s=simplify)
obs365=sapply(1:1000,function(i){mc.list[[i]][365]})

#вычислим стационарное распределение
#нужно решить СЛАУ x * P=x
#x - стац. распр. вероятн. (вектор), P - матр. пер. вероятн.
(freqInf=c(7,16,12)/35)

## [1] 0.2000 0.4571 0.3429

#эмпирическое распр. вероятностей состояний погоды на 365-ый день 
freq365=table(obs365)
#эмпирическое распределение приближается к стационарному (см. как данное распределение изменилось по отношению к p0)
(freq365=freq365/1000)

## obs365
##  shiny cloudy  rainy 
##  0.190  0.452  0.358

Задание 3

Написать функцию моделирование однородной цепи Маркова (ОЦМ) со следующей спецификацией:

markovchain=function(n,Params){  
 
 p0=Params$p0;P=Params$P
 s=rep(0,n)#результирующая последовательность = ОЦМ
 
 
 mc=as.factor(s)  
 levels(s)=c("shiny","cloudy","rainy")  
 return(s)  

}

В качестве БСВ использовать функцию runif. Функцию markovchain включить в отчет.

Смоделировать 1 ОЦМ длины 365 (параметры p0,P - по варианту) и сделать следующее:

вывести первые 10 наблюдений
подсчитать частоты каждого состояния погоды

Смоделировать 1000 реализаций ОЦМ, каждая длины 365 с одинаковыми параметрами (параметры p0,P - по варианту).

Вычислить freq365 - относительные частоты (доли) каждого состояния на момент времени t=365.
Вычислить стационарное распределение freqInf и сравнить его с freq365: вероятности должны иметь близкие значения. Для вычисления freqInf решается СЛАУ с матрицей P.

Моделирование однородной цепи Маркова, а также СЛАУ, которое нужно решить можно посмотреть здесь.

Варианты

Параметры цепи Маркова

начальное распределение вероятностей: \[{{\pi }_{i}}=\text{P}\left\{ {{x}_{1}}=i \right\},\,i\in S\]
матрица вероятностей одношаговых переходов: \[P={{\left( {{p}_{i,j}} \right)}_{i,j\in S}};\quad {{p}_{i,j}}=\text{P}\left\{ \left. {{x}_{t+1}}=j \right|x{}_{t}=i \right\};\ \sum\nolimits_{j\in S}{{{p}_{i,j}}}=1,\ i\in S\].

Замечание

Параметру p0 соответствует \(\pi\).

Варианты

\(\pi =\left( \begin{matrix} 0.33333 \\ 0.33333 \\ 0.33334 \\ \end{matrix} \right)\), \(P=\left( \begin{matrix} 0.1 & 0.2 & 0.7 \\ 0.3 & 0.4 & 0.3 \\ 0.6 & 0.2 & 0.2 \\ \end{matrix} \right)\).
\(\pi =\left( \begin{matrix} 0.33333 \\ 0.33333 \\ 0.33334 \\ \end{matrix} \right)\), \(P=\left( \begin{matrix} 0.4 & 0.2 & 0.4 \\ 0.3 & 0.4 & 0.3 \\ 0.5 & 0 & 0.5 \\ \end{matrix} \right)\).
\(\pi =\left( \begin{matrix} 0.33333 \\ 0.33333 \\ 0.33334 \\ \end{matrix} \right)\), \(P=\left( \begin{matrix} 0.2 & 0.2 & 0.6 \\ 0.2 & 0.5 & 0.3 \\ 0.6 & 0.3 & 0.1 \\ \end{matrix} \right)\).
\(\pi =\left( \begin{matrix} 0.33333 \\ 0.33333 \\ 0.33334 \\ \end{matrix} \right)\), \(P=\left( \begin{matrix} 0.75 & 0 & 0.25 \\ 0.3 & 0.5 & 0.2 \\ 0.2 & 0.2 & 0.6 \\ \end{matrix} \right)\).
\(\pi =\left( \begin{matrix} 0.5 \\ 0 \\ 0.5 \\ \end{matrix} \right)\), \(P=\left( \begin{matrix} 0.1 & 0.2 & 0.7 \\ 0.2 & 0.3 & 0.5 \\ 0.3 & 0.4 & 0.3 \\ \end{matrix} \right)\).
\(\pi =\left( \begin{matrix} 1 \\ 0 \\ 0 \\ \end{matrix} \right)\), \(P=\left( \begin{matrix} 0.1 & 0.2 & 0.7 \\ 0.3 & 0.4 & 0.3 \\ 0.6 & 0.2 & 0.2 \\ \end{matrix} \right)\).
\(\pi =\left( \begin{matrix} 0.2 \\ 0.3 \\ 0.5 \\ \end{matrix} \right)\), \(P=\left( \begin{matrix} 0.4 & 0.2 & 0.4 \\ 0.3 & 0.4 & 0.3 \\ 0.5 & 0 & 0.5 \\ \end{matrix} \right)\).
\(\pi =\left( \begin{matrix} 0.25 \\ 0.5 \\ 0.25 \\ \end{matrix} \right)\), \(P=\left( \begin{matrix} 0.2 & 0.2 & 0.6 \\ 0.2 & 0.5 & 0.3 \\ 0.6 & 0.3 & 0.1 \\ \end{matrix} \right)\).
\(\pi =\left( \begin{matrix} 0.5 \\ 0 \\ 0.5 \\ \end{matrix} \right)\), \(P=\left( \begin{matrix} 0.75 & 0 & 0.25 \\ 0.3 & 0.5 & 0.2 \\ 0.2 & 0.2 & 0.6 \\ \end{matrix} \right)\).
\(\pi =\left( \begin{matrix} 1 \\ 0 \\ 0 \\ \end{matrix} \right)\), \(P=\left( \begin{matrix} 0.1 & 0.2 & 0.7 \\ 0.2 & 0.3 & 0.5 \\ 0.3 & 0.4 & 0.3 \\ \end{matrix} \right)\).
\(\pi =\left( \begin{matrix} 0.2 \\ 0.3 \\ 0.5 \\ \end{matrix} \right)\), \(P=\left( \begin{matrix} 0.1 & 0.2 & 0.7 \\ 0.3 & 0.4 & 0.3 \\ 0.6 & 0.2 & 0.2 \\ \end{matrix} \right)\).
\(\pi =\left( \begin{matrix} 0.25 \\ 0.5 \\ 0.25 \\ \end{matrix} \right)\), \(P=\left( \begin{matrix} 0.4 & 0.2 & 0.4 \\ 0.3 & 0.4 & 0.3 \\ 0.5 & 0 & 0.5 \\ \end{matrix} \right)\).
\(\pi =\left( \begin{matrix} 0.5 \\ 0 \\ 0.5 \\ \end{matrix} \right)\), \(P=\left( \begin{matrix} 0.2 & 0.2 & 0.6 \\ 0.2 & 0.5 & 0.3 \\ 0.6 & 0.3 & 0.1 \\ \end{matrix} \right)\).
\(\pi =\left( \begin{matrix} 1 \\ 0 \\ 0 \\ \end{matrix} \right)\), \(P=\left( \begin{matrix} 0.75 & 0 & 0.25 \\ 0.3 & 0.5 & 0.2 \\ 0.2 & 0.2 & 0.6 \\ \end{matrix} \right)\).

Нахождение стационарного распределения (\(x=({{x}_{1}},{{x}_{2}},{{x}_{3}})\)): решить систему \[\left\{ \begin{matrix} {P}'x=x \\ {{x}_{1}}+{{x}_{2}}+{{x}_{3}}=1 \\ \end{matrix} \right.\]

Замечания

При решении системы исключить одно уравнение из \({P}'x=x\).
Стационарному распределению вероятностей freqInf (в примере) соответстует решение СЛАУ х.

Задание 4. Моделирование непрерывных случайных величин

Примеры

#установим начальное значения датчика псевдослучайных чисел
#для воспроизводимости результатов моделирования
#устанавливать перед первым вызовом функции генерации СВ
set.seed(19)

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

#сгенерируем выборку норм. распр. СВ с МО=5 и Ст.откл.=2
x=rnorm(n=10000,mean = 5,sd = 2)
#построим гистограмму
hist(x,probability = T)

plot of chunk unnamed-chunk-12

#выведем числовые характеристики
mean(x) #математическое ожидание (МО)

## [1] 4.988

sd(x) #стандртное отклонение (СО)

## [1] 1.984

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

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

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

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

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

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

## [1] 0.01293

kurtosis(x) #эксцесс: острота пика прибл. соотв. таковой для норм. распр.

## [1] 0.01875

Тест согласия Колмогорова-Смирнова

Тест Колмогорова-Смирнова:  на выборке СВ х проверяется гипотеза о соответствии указанному закону распределения вероятностей. Требуется указать функцию раcпределения вероятностей (ФРВ,distribution function) для проверяемого закона и значения параметров распределения.

#тест Колмогорова-Смирнова
#указана ФРВ норм. закона (pnorm) и параметры
#p-value > 0.05: на уровне значимости гипотеза не отклоняется
ks.test(x,pnorm,mean=5,sd=2)

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  x
## D = 0.0075, p-value = 0.6188
## alternative hypothesis: two-sided

Распределение Хи-квадрат пирсона

#параметры распределения: df - число степеней свободы (degrees of freedom)
k=5
x=rchisq(n=10000,df = k)
hist(x,probability = T)

plot of chunk unnamed-chunk-15

#числовые характеристики: через точку с запятой указать
#эмпирические и теоретичиские характеристики соответственно
mean(x);k #математическое ожидание

## [1] 4.994

## [1] 5

sd(x);sqrt(2*k) #стандартное отклонение

## [1] 3.149

## [1] 3.162

skewness(x);sqrt(8/k) #коэфф. асимметрии

## [1] 1.242

## [1] 1.265

kurtosis(x);12/k #коэфф. эксцесса

## [1] 2.2

## [1] 2.4

#проверка теста согласия
ks.test(x,pchisq,df=k)

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  x
## D = 0.0049, p-value = 0.9698
## alternative hypothesis: two-sided

Выводы о форме распределения:

Распределение асимметрично: имеет тяжелый левый хвост.
Распределение имеет более острый пик относительно нормального распределения.

Задание 4

Реализовать функцию генерации вектора СВ из непрерывного закона распределения вероятностей (НРВ). Закон РВ и параметры - согласно варианту. Общая спецификация функции:

rdistr2=function(n,par1,par2){

}

Замечания

1. Вместо подстроки "distr" в названии функции использовать принятое в R обозначение для закона распределения СВ (norm, chisq). Принятые обозначения можно посмотреть в разделе справки "Distributions" (Вкладка Help). В качестве аргументов par1, par2 использовать названия параметров согласно варианту (если у распределения другое количество парамеметров, то количество аргументов может изменяться).
2. В качестве датчика БСВ (стандартное равномерной СВ) использовать функцию **runif** (r=runif(1) дает одну реализацию БСВ), а датчика стандартной нормальной величины - функцию rnorm (u=rnorm(1)).  
3. Алгоритм моделирования взять из книги  
Вадзинский, Р.Н. Справочник по вероятностным распределениям.

Смоделировать выборку из n=10000 СВ. Провести следующий анализ:

построить гистограмму
вычислить МО, дисперсию, коэфф-ты асимметрии и эксцесса (если они существуют для данного РВ) и сравнить их с теоретическими значениями.
проверить моделируемый закон распределения с помощью теста Колмогорова-Смирнова (в качестве теорет. ФРВ применить встроенную функцию pdistr, где distr - обозначение для соотв. закона РВ, например norm, тогда ФРВ будет pnorm).

В конце отчета сделать выводы относительно формы распределения на основе коэффициентов асимметрии и эксцесса.

Варианты

Вариант	Распределение	Обозначение	Алгоритм	Параметры в R	Соотв. в книге
1	Нормальное	norm	2	mu=-3,sd=2	\(\mu,\sigma\)
2	Коши	cauchy	1	location=-2,scale=0.5	\(\mu,\lambda\)
3	Логистическое	logis	1	location=2,scale=0.64	\(\mu,\lambda\)
4	Экспоненциальное	exp	1	rate=0.4	\(\lambda\)
5	Гамма	gamma	1	shape=4,rate=3	\(\alpha,\lambda\)
6	Вейбулла-Гнеденко	weibull	1	shape=2.5,scale=2	c,a
7	Логнормальное	lnorm	1	meanlog=0,sdlog=0.5	\(\mu,a\)
8	Бета	beta	1	shape1=4,shape2=2	u,v
9	\({{\chi }^{2}}\)-Пирсона	chisq	1	df=10	n
10	Релея	rayleigh	1	scale=2	a
11	Максвелла	maxwell	1	a=1	a
12	Парето	pareto	1	location=1,shape=2	\(x_{0},\alpha\)
13	Лапласа	laplace	1	location=2,scale=0.8	\(\mu,\lambda\)
14	\({{\chi }^{2}}\)-Пирсона	chisq	2	df=5	n

Замечание

При проверке гипотезы согласия с помощью теста Колмогорова-Смирнова для вариантов 10-13 нужны теоретические ФРВ prayleigh, pmaxwell,…, которые могут быть найдены в библиотеке VGAM. Для использования данных функций необходимы следующие действия:

Установить библиотеку VGAM: install.packages(“VGAM”).
Загрузить данную библиотеку: library(VGAM).

Задание 5. Статистический анализ данных

Указания

1. Для решения данного задания воспользоваться шаблоном "ex5_template.Rmd".   
2. Для справки использовать примеры кода из Разделов 1-3.   
3. Использовать свой набор данных из папки "DATA" (на "\\ Fpmi-stud\Subfaculty\Кафедра ММАД\Новопольцев А.Ю\МСАД &MCЭ"). Заполнить шаблон необходимым кодом для решения задач и скомпилировать в HTML или WORD (см. инструкцию в "Лабораторный практикум по МСАДиМСЭ.docx").

Управление данными

Ввести данные из текстового файла согласно варианту (функция “read.table”). Сохранить в переменную qc. Для этого поместить файл “QC#.txt” из папки “DATA” в рабочую директорию (узнать путь можно вызовом getwd() в консоли).

Распечатать таблицу данных в сокращенном виде (str)

## 'data.frame':    100 obs. of  9 variables:
##  $ V1: num  50.8 54.7 57.7 56.4 52.7 ...
##  $ V2: int  2 2 5 1 6 2 2 2 3 0 ...
##  $ V3: int  2 3 1 1 2 5 2 4 1 3 ...
##  $ V4: int  1 1 1 1 2 1 1 2 1 2 ...
##  $ V5: int  4 4 4 3 3 2 4 4 3 3 ...
##  $ V6: num  31.2 31.9 31 30.4 32.8 ...
##  $ V7: num  16.6 18.2 19.5 18.8 17.8 ...
##  $ V8: num  4.57 4.55 5.76 4.72 4.2 ...
##  $ V9: int  4 3 3 5 3 6 3 3 2 3 ...

Вывести размерность таблицы (dim)

## [1] 100   9

Вывести первые 3 записи (head)

##      V1 V2 V3 V4 V5    V6    V7    V8 V9
## 1 50.75  2  2  1  4 31.18 16.64 4.565  4
## 2 54.70  2  3  1  4 31.88 18.21 4.555  3
## 3 57.66  5  1  1  4 30.96 19.54 5.765  3

Изменить имена следующих переменных (применить функцию names к переменной qc): “б”-“kit”, “п”-“nbug”,“№изг”-“maker”,“№пост”-“vendor”,“№вд”-“bug” “н” - последов. в строке “с1”, “c2”,“c3”,“c4” Вывести часть таблицы функцией head.

##      c1 kit nbug maker vendor    c2    c3    c4 bug
## 1 50.75   2    2     1      4 31.18 16.64 4.565   4
## 2 54.70   2    3     1      4 31.88 18.21 4.555   3
## 3 57.66   5    1     1      4 30.96 19.54 5.765   3
## 4 56.43   1    1     1      3 30.42 18.78 4.725   5
## 5 52.69   6    2     2      3 32.76 17.76 4.200   3
## 6 55.26   2    5     1      2 32.56 18.68 5.365   6

Вычисление статистик

Вывести общую статистику по переменным (summary)

##        c1            kit            nbug          maker     
##  Min.   :49.9   Min.   :0.00   Min.   :0.00   Min.   :1.00  
##  1st Qu.:53.5   1st Qu.:2.00   1st Qu.:1.00   1st Qu.:1.00  
##  Median :54.9   Median :2.00   Median :2.00   Median :1.00  
##  Mean   :54.7   Mean   :2.64   Mean   :2.09   Mean   :1.48  
##  3rd Qu.:56.3   3rd Qu.:4.00   3rd Qu.:3.00   3rd Qu.:2.00  
##  Max.   :58.0   Max.   :6.00   Max.   :7.00   Max.   :2.00  
##      vendor           c2             c3             c4      
##  Min.   :1.00   Min.   :29.5   Min.   :16.5   Min.   :3.64  
##  1st Qu.:2.00   1st Qu.:31.9   1st Qu.:17.8   1st Qu.:4.72  
##  Median :3.00   Median :33.1   Median :18.3   Median :5.04  
##  Mean   :2.87   Mean   :33.1   Mean   :18.2   Mean   :4.99  
##  3rd Qu.:4.00   3rd Qu.:34.3   3rd Qu.:18.7   3rd Qu.:5.31  
##  Max.   :4.00   Max.   :37.4   Max.   :19.7   Max.   :5.95  
##       bug      
##  Min.   :2.00  
##  1st Qu.:3.00  
##  Median :3.00  
##  Mean   :3.66  
##  3rd Qu.:4.25  
##  Max.   :7.00

Преобразовать категоральные переменные maker и vendor в переменные-факторы (c помощью as.factor).
Снова воспользоваться summary и сравнить результаты для преобразованных переменных.

##        c1            kit            nbug      maker  vendor       c2      
##  Min.   :49.9   Min.   :0.00   Min.   :0.00   1:52   1:14   Min.   :29.5  
##  1st Qu.:53.5   1st Qu.:2.00   1st Qu.:1.00   2:48   2:23   1st Qu.:31.9  
##  Median :54.9   Median :2.00   Median :2.00          3:25   Median :33.1  
##  Mean   :54.7   Mean   :2.64   Mean   :2.09          4:38   Mean   :33.1  
##  3rd Qu.:56.3   3rd Qu.:4.00   3rd Qu.:3.00                 3rd Qu.:34.3  
##  Max.   :58.0   Max.   :6.00   Max.   :7.00                 Max.   :37.4  
##        c3             c4            bug      
##  Min.   :16.5   Min.   :3.64   Min.   :2.00  
##  1st Qu.:17.8   1st Qu.:4.72   1st Qu.:3.00  
##  Median :18.3   Median :5.04   Median :3.00  
##  Mean   :18.2   Mean   :4.99   Mean   :3.66  
##  3rd Qu.:18.7   3rd Qu.:5.31   3rd Qu.:4.25  
##  Max.   :19.7   Max.   :5.95   Max.   :7.00

Присоединить таблицу данных qc к списку текущих переменных. Далее можно обращаться непосредственно к именам переменных.

attach(qc)

Для переменной c1 вывести выборочный квантиль уровня 0.9 (quantile)

##   90% 
## 56.99

Для переменной с1 вычислить дисперсию (var)

## [1] 3.414

Вычислить дисперсию (var) для переменной c1*2. Сравнить результат с предыдущим

## [1] 13.66

Для переменных с1,с3 вычислить ковариацию (cov)

## [1] 1.132

Вывести корреляционную матрицу по всем числовым характеристикам с1-c4 (cor). Для выбора нужных столбцов таблицы воспользоваться оператором “[]”.
Между какими переменными есть значимая корреляция?

##          c1        c2       c3        c4
## c1 1.000000  0.003348  0.87507  0.539350
## c2 0.003348  1.000000 -0.01003 -0.009796
## c3 0.875068 -0.010034  1.00000  0.583419
## c4 0.539350 -0.009796  0.58342  1.000000

Представить матрицу корреляции в более удобном виде(symnum, в качестве аргумента данной функции использовать результат cor)

##    c1 c2 c3 c4
## c1 1          
## c2    1       
## c3 +     1    
## c4 .     .  1 
## attr(,"legend")
## [1] 0 ' ' 0.3 '.' 0.6 ',' 0.8 '+' 0.9 '*' 0.95 'B' 1

Построить диаграмму рассеяния с линией регрессии для пары переменных с наибольшей корреляцией (последовательно вызвать функции plot, abline)
plot of chunk unnamed-chunk-29

Анализ распределения вероятностей

Проверить распределение у переменной c1 Построить гистограмму для переменной c1 (hist) plot of chunk unnamed-chunk-30

Вывести общую статистику по переменной c1 (summary)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    49.9    53.5    54.9    54.7    56.3    58.0

#загрузка библиотеки: предварительно установить! install.packages("e1071")
library(e1071) #загрузка необходимой библиотеки

Вывести асимметрию для переменной с1 (skewness)

## [1] -0.4038

Вывести эксцесс для переменной с1 (kurtosis)

## [1] -0.4987

Применить тест Колмогорова-Смирнова (ks.test) для проверки переменной c1 на нормальное распределение. (первый параметр - проверяемая переменная, второй параметр - ФРВ проверяемого закона, mean,sd - выборочные оценки параметров)

## Warning: ties should not be present for the Kolmogorov-Smirnov test

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  c1
## D = 0.061, p-value = 0.8509
## alternative hypothesis: two-sided

Выводы:

Распределение слегка скошено вправо (Асимметрия Sk<0) и имеет более тяжелые хвосты ( Эксцесс Ex<0) по сравнению с нормальным распределением, имеющим в качестве значений параметров их выборочные оценки.
Согласно статистическим тестам, гипотеза о нормальном распределении не отклоняется на общепринятых уровнях значимости (0.05, 0.1 и т.д.). Вероятно, наблюдаемые отрицательные значения асимметрии и эксцесса связаны с недостаточно большим объемом выборки.

Отсоединить таблицу данных qc от списка текущих переменных (работа с данными закончена).

detach(qc)

Задания

Alexander Novopoltsev

Thursday, September 18, 2014

Задание 1. Моделирование БСВ

Примеры

Сравнение с распределением цифр числа пи

n=1000 БСВ

Задание 1

Варианты

Задание 2. Моделирование дискретных случайных величин

Примеры

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

Распределение Пуассона

Задание 2

Варианты

Задание 3. Моделирование цепи Маркова

Пример

Задание 3

Варианты

Задание 4. Моделирование непрерывных случайных величин

Примеры

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

Распределение Хи-квадрат пирсона

Задание 4

Варианты

Задание 5. Статистический анализ данных

Управление данными

Вычисление статистик

Анализ распределения вероятностей