-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathKood.R
175 lines (113 loc) · 4.18 KB
/
Kood.R
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
set.seed(32)
####n = 100
n = 100
a = 2.0
b = 1.0/3.0
theta = rgamma(n,a,b)
head(theta)
tail(theta)
hist(theta, freq=FALSE)
curve(dgamma(x,a,b), col = "blue", add = TRUE)
sum(theta)/n
mean(theta) #simuleeritud keskmine
a/b #see on antud gamma jaotuse analüütiline keskmine. Meie simuleeritud keskmine pole väga hea
###n = 10000
n = 10000
a = 2.0
b = 1.0/3.0
theta = rgamma(n,a,b)
head(theta)
tail(theta)
hist(theta, freq=FALSE)
curve(dgamma(x,a,b), col = "blue", add = TRUE)
mean(theta)
var(theta) #simuleeritud variatsioon
a/b^2 #analüütiline variatsioon - suht lähedal
ind <- theta < 5.0 #millised väärtused on simulatsioonis alla 5
head(ind)
mean(ind) # simulatsiooni hinnang, et väärtus on alla 5.0
pgamma(5,2,1/3) #analüütiline hinnang, et väärtus on alla 5.0
quantile(theta,.95) #paneb theta vektori numbrid järjekorda ja võtab 95th quantile'i
qgamma(.95,2,1/3) #analüütiline meie gamma quantile
se <- sd(theta)/sqrt(n) #monte carlo hinnangu standardviga
se
2.0*se #usalduspiiride ulatus keskmisele. Ütleb, et oleme 95% kindlad, et tõeline theta väärtus pole
#kaugemal kui see väärtus keskmisest:
mean(theta) - (2.0*se)
mean(theta) + (2.0*se)
se <- sd(ind)/sqrt(n)
se #usalduspiiride ulatus tõenäosusele, et väärtus on alla 5.0
#ehk proportsioonile väärtustele, mis on simulatsioonis alla 5.0.
2*se #Peaks ütlema meile sisuliselt, et 95% tõenäosus, et kui teeme simulatsiooni,
#siis alla 5.0 jäävate väärtuste proportsioon jääb meie mean(ind) väärtusest nii palju
#mõlemale poole (+ ja -)
####HIERHARHILINE MUDEL
##2 sammu:
m = 1e5 #Montecarlo valimi suurus
y <- numeric(m) #tühi numbriline vektor m liikmega
phi <- numeric(m) #---"---
#järgmises funktsioonis me sisuliselt tõmbame ükshaaval phi vektorisse ühe väärtuse
#beta jaotusest ja paneme selle väärtuse binomiaalse jaotuse parameetriks, et tõmmata
#y vektorisse üks väärtus. Teeme seda 1e5 ehk 100k korda.
for (i in 1:m){
phi[i] <- rbeta(1, 2, 2)
y[i] <- rbinom(1,10,phi[i])
}
#Siin teeme sama asja, mis üleval ainult vektoritega, mis on vääääga palju kiirem
#sisuliselt teeme kohe 100k tõmmet phi vektorisse ja seejärel kasutame igat phi
#vektori liiget ükshaaval, et teha igale phi vektori väärtusele vastav tõmmet
#binomiaalsest jaotusest y vektorisse ja nii 100k korda.
phi <- rbeta(m,2,2)
y <- rbinom(m,10,phi)
#Nüüd vaatame y marginaalse jaotuse omadusi. Sisuliselt lihtalt ignoreerime
#Phi-d ja vaatame y-d Aga pea meeles, et y marginaalne jaotus pole binomiaalne
#vaid beta-bionmiaalne, sest parameeter phi ei ole püsiv.
table(y) #Vaatame kui palju me erinevaid väärtuseid saime y vektorisse
plot(table(y)/m) #plotime tõenäosused saada mingit y väärtust
mean(y) # samuti saame leida marginaalne y beta-binomiaalse jaotuse keskmise
m <- 1e5
theta <- numeric(m)
theta <- rbeta(m,5,3)
theta <- theta/(1-theta)
0.625/(1-0.625)
head(theta)
mean(theta>1)
qnorm(0.3,0,1)
monorm <- rnorm(m,0,1)
quantile(monorm, 0.3)
sqrt(5.2/5000)
Q = matrix(c(0.0, 1,
0.3, 0.7),
nrow=2, byrow=TRUE)
n = 100000
x = numeric(n)
x[1] = 1 # fix the state as 1 for time 1
for (i in 2:n) {
x[i] = sample.int(2, size=1, prob=Q[x[i-1],])}
table(x)/n
vec <- c(1, 0)
vec %*% Q %*% Q %*% Q
vec %*% Q
W = matrix(c(.365, .635,
1, 0),
nrow=2, byrow=TRUE)
W30 = W
for (i in 2:30) {
W30 = W30 %*% W
}
round(W30, 3)
install.packages("rjags")
update_mu = function(n, ybar, sig2, mu_0, sig2_0) {
sig2_1 = 1.0 / (n / sig2 + 1.0 / sig2_0)
mu_1 = sig2_1 * (n * ybar / sig2 + mu_0 / sig2_0)
rnorm(n=1, mean=mu_1, sd=sqrt(sig2_1))
}
#siin all on nu_0 ja beta_0 eeljaotuse parameetrid, mille me paika paneme eelnevalt
update_sig2 = function(n, y, mu, nu_0, beta_0) {
nu_1 = nu_0 + n / 2.0
sumsq = sum( (y - mu)^2 ) # vectorized
beta_1 = beta_0 + sumsq / 2.0
out_gamma = rgamma(n=1, shape=nu_1, rate=beta_1) # rate for gamma is shape for inv-gamma kuna r ei võimalda otse gammast tõmmata tõmbame algusest gammast ja siis võtame inverse'i, mis on sama kui otse gammast tõmmata.
1.0 / out_gamma # reciprocal of a gamma random variable is distributed inv-gamma
}
install.packages("rjags")