-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathbayes-GEV
142 lines (121 loc) · 5.49 KB
/
bayes-GEV
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
######################################################################################################################################
## unnormalized log-posterior function , gradient vector (mu delta xi)' = x = (x[1] x[2] x[3])
## independent normal priors ~ N[0, 100]
######################################################################################################################################
### log-posterior for GEV
######################################################################################################################################
# ll <- function(theta, x) {
# mu = theta[1]
# delta = theta[2]
# xi = theta[3]
# ll <- sum(log(f(x, mu, exp(delta), xi))) + sum(dnorm(theta, mean = 0, sd = 5, log = T))
# ifelse(ll != "-Inf", ll, -100000)
# }
### gradient log-posterior for GEV
######################################################################################################################################
# nll <- function(theta, x) {
# D <- rep(0, 3)
# mu = theta[1]
# delta = theta[2]
# xi = theta[3]
# n = length(x)
# sig = exp(delta)
# d = x - mu
# z = 1 + xi*d/sig
# w = d/sig^2
# q = d/sig
#
# D[1] = 1/sig*sum(z^(-1)*((1+xi) - z^(-1/xi))) - mu/25
# D[2] = (sum((1+xi)*w*z^(-1) - w*z^(-1-1/xi))*sig - n) - delta/25
# D[3] = suppressWarnings(1/xi^(2)*sum(log(z)) - (1/xi+1)*sum(q*z^(-1)) +
# 1/xi*sum(q*z^(-1-1/xi)) - 1/xi^(2)*sum(log(z)*z^(-1/xi)) - xi/25)
# if (suppressWarnings(any(z < 0 | D == "NaN" | D == "Inf" | D == "-Inf"))) rep(-100000, 3)
# else D
# }
### Expected Fisher Information Matrix and its derivatives
######################################################################################################################################
G <- function(theta, x) {
G = matrix(, 3, 3)
mu = theta[1]
delta = theta[2]
sig = exp(theta[2])
xi = theta[3]
n = length(x)
pg <- gamma(2 + xi)
p <- (1 + xi)^(2) * gamma(1 + 2*xi)
q <- pg * (psigamma(1 + xi) + 1 + 1/xi)
G[1, 1] = p/sig^2 + 1/25
G[1, 2] = -1/(sig*xi)*(p - pg)
G[1, 3] = -1/(sig*xi)*(q - p/xi)
G[2, 1] = G[1, 2]
G[2, 2] = 1/(xi^2) * (1 - 2*pg + p) + 1/25
G[2, 3] = -1/(xi^2) * (1 - J + (1 - pg)/xi - q + p/xi)
G[3, 1] = G[1, 3]
G[3, 2] = G[2, 3]
G[3, 3] = 1/(xi^2) * ((pi^2)/6 + (1 - J + 1/xi)^2 - 2*q/xi + p/(xi^2)) + 1/25
return(n*G)
}
### maximum a posteriori
######################################################################################################################################
# map = optim(c(0.1, 0, 0.1), ll, gr = nll, x = y, hessian = T, method = "BFGS", control = list(fnscale = -1, maxit = 10000)); map
#
# G. = (G(map$par, y))
# invG = solve(G.)
### example HMC data port-pirie
#####################################################################
# data(portpirie)
# y = portpirie[, 2]
# n = length(y)
#
# pdf('histex1.pdf', width = 6, height = 5)
# hist(y, prob = T, xlab = "Nível do mar [metros]", ylab = 'Frequência relativa', main = '')
# dev.off()
# pdf('stport.pdf', width = 6, height = 5)
# plot(y, ylab = "Nível do mar [metros]", type = 'l', xlab = 'Tempo', main = '')
# dev.off()
#
# pdf('acfport.pdf', width = 6, height = 5)
# acf(portpirie[, 2], xlab = 'Tempo', main = '')
# dev.off()
# pdf('HMCex1.pdf', width = 7, height = 5)
# par(mfrow=c(3, 1), mar = c(4, 5, 3, 3))
# plot(B$theta[, 1], ylab = expression(mu^(i)), xlab = 'i', type = 'l', main = 'HMC')
# plot(B$theta[, 2], ylab = expression(sigma^(i)), xlab = 'i', type = 'l')
# plot(B$theta[, 3], ylab = expression(xi^(i)), xlab = 'i', type = 'l')
# dev.off()
# pdf('acfHMCex1.pdf', width = 7, height = 5)
# par(mfrow=c(3, 1), mar = c(4, 5, 3, 3))
# acf(B$theta[, 1], lag.max = 100, ylab = expression(r[mu~","~k]), main = 'HMC Autocorrelação', xlab = '')
# acf(B$theta[, 2], lag.max = 100, ylab = expression(r[sigma~","~k]), main = '', xlab = '')
# acf(B$theta[, 3], lag.max = 100, ylab = expression(r[xi~","~k]), xlab = 'defasagem k', main = '')
# dev.off()
#
# pdf('MHex1.pdf', width = 7, height = 5)
# par(mfrow=c(3, 1), mar = c(4, 5, 3, 3))
# plot(as.numeric(m1[, 1]), ylab = expression(mu^(i)), xlab = 'i', type = 'l', main = 'MH')
# plot(as.numeric(m1[, 2]), ylab = expression(sigma^(i)), xlab = 'i', type = 'l')
# plot(as.numeric(m1[, 3]), ylab = expression(xi^(i)), xlab = 'i', type = 'l')
# dev.off()
# pdf('acfMHex1.pdf', width = 7, height = 5)
# par(mfrow=c(3, 1), mar = c(4, 5, 3, 3))
# acf(as.numeric(m1[, 1]), lag.max = 100, ylab = expression(r[mu~","~k]), main = 'MH Autocorrelação', xlab = '')
# acf(as.numeric(m1[, 2]), lag.max = 100, ylab = expression(r[sigma~","~k]), main = '', xlab = '')
# acf(as.numeric(m1[, 3]), lag.max = 100, ylab = expression(r[xi~","~k]), xlab = 'defasagem k', main = '')
# dev.off()
#
# ###########################
# ### return-period inference
# ###########################
# x = seq(min(y)-.5, max(y)+.5, length.out = 100)
# dens.est = matrix(nrow = dim(B$theta)[1], ncol = length(x) )
# quan.est = matrix(nrow = length(x), ncol = 2)
#
# for (i in 1:length(x)) {
# dens.est[ ,i] = f(x[i], B$theta[ ,1], B$theta[ ,2], B$theta[ ,3])
# quan.est[i, ] = quantile(dens.est[ ,i], prob = c(0.05, 0.95))
# }
# hist(y, prob = T)
# lines(x, apply(dens.est, 2, mean), col = 'red')
# lines(x, quan.est[, 1], lty = 2, col = 'red')
# lines(x, quan.est[, 2], lty = 2, col = 'red')
# points(y, rep(0, n), col = 'blue', pch = 4)