forked from yordanbabukov/Statistics-and-probabilities-notes
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy path6-new.txt
177 lines (149 loc) · 10.5 KB
/
6-new.txt
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
Bootstrap sampling
Bootstrapping е метод за генериране на събития от множество с данни на базата, на които може да се направи статистически извод.
Идеята е, че като проиграем дадена игра (хвърляне на зар, раздаване на карти) няколко пъти, можем да си изведем извод за очакванията ни от нея.
Процесът представлява многократно генериране на събития от извадка, след което формираме статистика върху настъпилите събития.
Да го приложим за данните faithful.
install.packages('UsingR')
library(UsingR)
data(faithful)
hist(faithful$eruptions)
hist(faithful$eruptions, breaks=25)
Разпределението е странно.
dev.new()
hist(sample(faithful$eruptions,100, replace=TRUE), breaks=25)
Резултатът е подобно разпределение на оригиналните данни, но не същото.
#hist(sample(faithful$eruptions,272), breaks=25)
dev.off()
dev.off()
Пример за приложение на bootstrapping:
set.seed(333); x<-rnorm(30) #така всеки път генерираната поредица числа ще е една и съща
bootMean <- rep(NA,1000) #генерираме 1000 празни стойности в bootMean.
sampledMean <- rep(NA,1000)
for(i in 1:1000){bootMean[i] <- mean(sample(x,replace=TRUE))}
for(i in 1:1000){sampledMean[i] <- mean(rnorm(30))}
plot(density(bootMean))
lines(density(sampledMean), col='red')
Извод:
Имаме данни за 30 човека.
Правим 1000 варианта на тази извадка.
На симулираните резултати, може да извадим извод за съвсем друга извадка.
Създаване на скрипт:
Създаваме текстов файл (с произволен текстов редактор)
Добавяме съдържание, например:
example <- function(max=5)
{
for (i in 1:max)
{
print(i)
}
}
Свързване на скрипта към работната среда:
source("пълно_име_на_файл") или относително име комбинирано с .libPath("...")
Извикване на скрипта:
example(15)
# Section 7 #
---------Simulations---------
Central limit theorem
Централната гранична теорема казва:
Независимо от първоначалното разпределението на данните, с които работим, ако то има дефинирани средно и стандартно отклонение и вземем достатъчно количество извадки (samples) резултатът от някакво тяхно характеризиране ще бъде нормално разпределен.
Пример:
данните за изригванията на гейзер:
data(faithful)
hist(faithful$eruptions, breaks=25)
Очевидно не са нормално разпределени.
Да приложим централната гранична теорема:
samplingMean <- rep(NA,1000)
x <- faithful$eruptions
for(i in 1:1000){samplingMean[i] <- mean(sample(x,replace=TRUE))}
hist(samplingMean, breaks=25)
Да видим дали важи за данни с друго разпределение, например експоненциално:
samplingMean <- rep(NA,1000)
x <- rexp(1000,1/3)
hist(x, breaks=25)
for(i in 1:1000){samplingMean[i] <- mean(sample(x,replace=TRUE))}
hist(samplingMean, breaks=25)
Ако вземем нормално разпределени или равномерно разпределени данни, резултатът ще е същият.
Разпределение получено по този начин се нарича "sampling distribution".
Тълкуване на средното и стандартно отклонение на sampling distribution:
Средното се очаква да е около средното на извадката.
Стандартното отклонение на средните = стандартното на извадката / sqrt(n)
Ако стандартното отклонение е по-голямо, тогава графиката ще е по-плоска. Ако е по-малко, ще е по-висока.
x <-c(1,2,3,5,6,9)
mean(x)
sd(x)
for(i in 1:1000){samplingMean[i] <- mean(sample(x,replace=TRUE))}
mean(samplingMean)
sd(samplingMean)
sd(x)/sqrt(6)
Практическа задача:
Средно човек от мъжки пол приема по 2 литра вода през спортно активен ден, със стандартно отклонение 0.7 литра.
Планираме провеждане на спортно събитие в планината на 50 човека и сме се запасили със 110 литра вода.
Каква е вероятността водата да не ни стигне?
mean = 2
sd = 0.7
Графика със средно 2 и ст. отклонение 0.7
Трябва да изчислим: каква е вероятността водата не ни стигне.
Или
Вероятността да се изпие повече от 110 литра вода.
Или
Вероятността средно човек да изпие повече от колкото е предвидено: 110/50 = 2.2 литра вода.
Да изчислим стандартното отклонение на sampling distribution. sd = 0.7/sqrt(50) = 0.099
За да не ни стигне водата, трябва да пресметнем каква е вероятността да попаднем в графиката над 2.2 литра.
За целта изчисляваме колко стандартни отклонения клоним от средното. (z-score)
2.2 - 2 (стандартно отклонение - средното) = 0.2
0.2/0.099 = 2.02, което означава, че сме 2.02 стандартни отклонения от средното.
Таблица със z-scores: http://media.tumblr.com/tumblr_mb3781q8zL1r80ng4.png
Засичаме по таблицата и резултатът е: 0.9783, но това е вероятността водата да ни стигне! Затова взимаме:
1-0.9783 = 0.0217 или иначе казано: 2.17% водата няма да ни стигне.
#z-score#
z-score (още наричан standard score)
С колко се различаваме от средното.
z-score = 1, значи че бягаме с 1 стандартно отклонение
z-score=2, бягаме с 2 стандартни отклонения
Формула за изчисляване на z-score:
Z = (X - mean(X) ) / sd(X)
z-score = (70 - 60) / 15 = 10/15 #получаваме 70 точки на тест, средно точките на тест са 60, с нормално отклонение 15. Извод: бягам с 2/3 от средното очакване.
#z-score#
Написано на R изглежда така:
Бърз вариант:
1-pnorm(2.2, 2, 0.099)
По-детайлен:
x <- rnorm(50, 2, 0.7)
mean(x)
2.161451
sd(x)
0.5783829 #Това си е наш експеримент, затова показателите са малко по-различни.
samplingMean <- rep(NA,10000)
for(i in 1:10000){samplingMean[i] <- mean(sample(x,replace=TRUE))}
mean(samplingMean)
2.13432
sd(samplingMean)
0.08007483
hist(samplingMean) #камбановидна! ЦГТ работи!
Средното си е около средното.
Стандарното отклонение е: 0.58/sqrt(50) ~ 0.082 (пак е близко, значи до тук всичко е ОК)
y <- 110/50 #средно колко вода трябва да изпие човек, за да ни стигне водата.
1-pnorm(y, mean(samplingMean) ,sd(samplingMean))
0.3111474
Извод: 31% е вероятността да не ни стигне водата.
Начин за демонстриране, дали данните ни са нормално разпределени е чрез така наречените "normal probability plots"
x=rnorm(100,0,1);
qqnorm(x,main='normal(0,1)'); #нарежда данните спрямо квантила, в който се намират.
qqline(x) #това не е резултат от линейна регресия, а права минаваща през 1ви и 3ти квантил на данните.
Разчитане на графиката: Ако несвързания граф прилича на права линия, то данните са приблизително нормално разпределени.
По ординатата са нашите стойности. По абсцисата са теоретичните квантили на нормалното разпределение.
С други думи: под Q1 очаква да има малко данни. Ако са повече става изкривяване.
Под Q2 очаква да има половината данни.
Под Q3 очаква да има 75% от данни
Над Q3 очаква да има малко данни, тъй като е нормално разпределение.
x=runif(100,0,1);
qqnorm(x,main='normal(0,1)'); #нарежда данните спрямо квантила, в който се намират.
qqline(x)
Тук данните не са в права линия, което индикира че не са нормално разпределение. Което е и вярно, тъй като тези данни са равномерно разпределени. Под Q2 и над Q3 има повече данни от очакваното за нормално разпределение.
qqnorm(faithful$eruptions, main='faithful eruptions')
Функция даваща обобщена информация: EDA - exploratory data analysis
library(UsingR) #зареждаме данните! (с install.packages ги зарежда по подразбиране..)
attach(homedata)
simple.eda(y1970)
simple.eda(y2000)
detach(homedata)