-
Notifications
You must be signed in to change notification settings - Fork 0
/
.Rhistory
194 lines (194 loc) · 6.5 KB
/
.Rhistory
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
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
require(Rcpp)
require(profvis)
require(RcppArmadillo)
require(Matrix)
require(expm)
require(std)
require(rbenchmark)
sourceCpp("RcppPoisson.cpp")
Tsigma2 <- 0.5
Tbeta0 <- log(10)
#Tamanho da amostra
set.seed(123)
Tn = 100
## Simular o efeito aleatório b_i
Tbi <- rnorm(n = Tn, mean = 0, sd = sqrt(Tsigma2))
## lambda
Tlambdai <- exp(Tbeta0 + Tbi)
Tyi <- rpois(n = Tn, lambda = Tlambdai)
#Integrandos
sum_cpp<-cppFunction(
"double sum_cpp(NumericVector x) {
const int n = x.size();
double y ;
for (int i=1; i < n; ++i) {
y = y+x[i];
}
return y;
}")
Tintegrando<-function(Tbi, Tbeta0, Tsigma2, Ty){
a= Exp_dpois(Ty,exp(Tbeta0+Tbi))
b= Ext_dnorm(Tbi, mean = 0, sd = sqrt(Tsigma2))
out <- exp(log(a) + log(b))
return(out)
}
Tintegrando2<-function(Tbi, Tbeta0, Tsigma2, Ty){
dpois(Ty, exp(Tbeta0 + Tbi))*(1/((sqrt(2 * pi * Tsigma2))*exp((Tbi^2)/(2*Tsigma2))))
}
Tintegrando3<-function(Tbi, Tbeta0, Tsigma2, Ty){
Ext_Integrando_NumericVector(Tbi,Tbeta0,Tsigma2,Ty) * dnorm(Tbi, mean = 0, sd = sqrt(Tsigma2))
}
TintegrandoNaMao<-function(Tbi, Tbeta0, Tsigma2, Ty){
lambda <- exp(Tbeta0+Tbi)
x <- log(lambda^Ty)
w <- Tbi^2
v <- (2*Tsigma2)
e <- -lambda - (w/v)
f <- log(sqrt(2*pi*Tsigma2))
z <- log(factorial(Ty))
exp(x + e - z -f)
}
TintegrandoSapply<-function(Tbi, Tbeta0, Tsigma2, Ty){
dpois(Ty, exp(Tbeta0 + Tbi))*(1/((sqrt(2 * pi * Tsigma2))*exp((Tbi^2)/(2*Tsigma2))))
}
integrando <- function(bi, beta0, sigma2, y) {
lambda <- exp(beta0 + bi)
out <- dpois(y, lambda = lambda)*dnorm(bi, mean = 0, sd = sqrt(sigma2))
return(out)
}
Ext_Integrando_Completo(Tbi,Tbeta0,Tsigma2,Tyi[1])
Ext_Integrando(Tbi,Tbeta0,Tsigma2,Tyi[1])
Tintegrando(Tbi,Tbeta0,Tsigma2,Tyi[1])
Tintegrando2(Tbi,Tbeta0,Tsigma2,Tyi[1])
Tintegrando3(Tbi,Tbeta0,Tsigma2,Tyi[1])
TintegrandoNaMao(Tbi,Tbeta0,Tsigma2,Tyi[1])
integrando(Tbi,Tbeta0,Tsigma2,Tyi[1])
Tintegrando(Tbi,Tbeta0,Tsigma2,Tyi[1])
Tintegrando2(Tbi,Tbeta0,Tsigma2,Tyi[1])
Tintegrando3(Tbi,Tbeta0,Tsigma2,Tyi[1])
TintegrandoNaMao(Tbi,Tbeta0,Tsigma2,Tyi[1])
my_ll <- function(par, y) {
integral <- c()
for(i in 1:length(y)) {
integral[i] <- integrate(f = integrando, lower = -Inf, upper = Inf,
beta0 = par[1], sigma2 = exp(par[2]), y = y[i])$value
}
#print(c(round(par, 2)))
ll <- sum(log(integral))
return(-ll)
}
Tmy_ll <- function(par, y){
-sum(sapply(y,function(y) log(integrate(f = Tintegrando, lower = -Inf, upper = Inf, Tbeta0 = par[1], Tsigma2 = exp(par[2]), Ty = y)$value)))
}
Tmy_ll2 <- function(par, y){
integral<- c()
for(i in 1:length(y)){
integral[i] <- integrate(f = Tintegrando2, lower = -Inf, upper = Inf, Tbeta0 = par[1], Tsigma2 = exp(par[2]), Ty = y[i])$value
}
-sum(log(integral))
}
Tmy_ll3 <- function(par, y){
integral<- c()
for(i in 1:length(y)){
integral[i] <- integrate(f = Ext_Integrando_NumericVector, lower = -Inf, upper = Inf, Tbeta0 = par[1], Tsigma2 = exp(par[2]), Ty = y[i])$value
}
-sum_cpp(log(integral))
}
Tmy_ll4 <- function(par, y){
integral<- c()
for(i in 1:length(y)){
integral[i] <- integrate(f = TintegrandoNaMao, lower = 0, upper = Inf, Tbeta0 = par[1], Tsigma2 = exp(par[2]), Ty = y[i])$value
}
-sum_cpp(log(integral))
}
Tmy_ll5 <- function(par, y){
-sum(sapply(-Inf:Inf,function(y) log(sum(Tintegrando2(par[1],exp(par[2])))), Ty = y))
}
Tmy_ll_sapply <- function(par, y){
-sum(sapply(y,function(y) log(integrate(f = Tintegrando, lower = -Inf, upper = Inf, Tbeta0 = par[1], Tsigma2 = exp(par[2]), Ty = y)$value)))
}
temp <- optim(par = c(log(8), log(0.3) ), fn = my_ll, y = yi,method = "Nelder-Mead")
temp <- optim(par = c(log(8), log(0.3) ), fn = my_ll, y = Tyi,method = "Nelder-Mead")
temp
Ttemp <- optim(par = c(log(8), log(0.3)), fn = Tmy_ll, y = Tyi, method = "Nelder-Mead")
Ttemp
Ttemp2 <- optim(par = c(log(8), log(0.3)), fn = Tmy_ll2, y = Tyi, method = "Nelder-Mead")
Ttemp2
Ttemp3 <- optim(par = c(log(8), log(0.3)), fn = Tmy_ll3, y = Tyi, method = "Nelder-Mead")
Ttemp3
Ttemp4 <- optim(par = c(log(8), log(0.3)), fn = Tmy_ll4, y = Tyi, method = "Nelder-Mead")
Ttemp4
Ttemp5 <- optim(par = c(log(8), log(0.3)), fn = Tmy_ll5, y = Tyi, method = "Nelder-Mead")
Ttemp5
Ttemp6 <- optim(par = c(log(8), log(0.3)), fn = Tmy_ll_sapply, y = Tyi, method = "Nelder-Mead")
Ttemp6
benchmark("Tmy_ll" = optim(par = c(log(8), log(0.3)), fn = Tmy_ll, y = Tyi, method = "Nelder-Mead"),
"Tmy_ll2" =optim(par = c(log(8), log(0.3)), fn = Tmy_ll2, y = Tyi, method = "Nelder-Mead"),
"Tmy_ll_sapply" =optim(par = c(log(8), log(0.3)), fn = Tmy_ll_sapply, y = Tyi, method = "Nelder-Mead"),
"my_ll" =optim(par = c(log(8), log(0.3)), fn = Tmy_ll_sapply, y = Tyi, method = "Nelder-Mead"),
replications = 5)
Ext_Integrando_Completo(Tbi,Tbeta0,Tsigma2,Tyi[1])
Ext_Integrando(Tbi,Tbeta0,Tsigma2,Tyi[1])
integrando(Tbi,Tbeta0,Tsigma2,Tyi[1])
Tintegrando(Tbi,Tbeta0,Tsigma2,Tyi[1])
Tintegrando2(Tbi,Tbeta0,Tsigma2,Tyi[1])
Tintegrando3(Tbi,Tbeta0,Tsigma2,Tyi[1])
TintegrandoNaMao(Tbi,Tbeta0,Tsigma2,Tyi[1])
my_ll <- function(par, y) {
integral <- c()
for(i in 1:length(y)) {
integral[i] <- integrate(f = integrando, lower = -Inf, upper = Inf,
beta0 = par[1], sigma2 = exp(par[2]), y = y[i])$value
}
#print(c(round(par, 2)))
ll <- sum(log(integral))
return(-ll)
}
Tmy_ll <- function(par, y){
-sum(sapply(y,function(y) log(integrate(f = Tintegrando, lower = -Inf, upper = Inf, Tbeta0 = par[1], Tsigma2 = exp(par[2]), Ty = y)$value)))
}
temp <- optim(par = c(log(8), log(0.3) ), fn = my_ll, y = Tyi,method = "Nelder-Mead")
temp
Ttemp <- optim(par = c(log(8), log(0.3)), fn = Tmy_ll, y = Tyi, method = "Nelder-Mead")
Ttemp
Ttemp2 <- optim(par = c(log(8), log(0.3)), fn = Tmy_ll2, y = Tyi, method = "Nelder-Mead")
Ttemp2
Ttemp3 <- optim(par = c(log(8), log(0.3)), fn = Tmy_ll3, y = Tyi, method = "Nelder-Mead")
Ext_Integrando_Completo(Tbi,Tbeta0,Tsigma2,Tyi[1])
Ext_Integrando(Tbi,Tbeta0,Tsigma2,Tyi[1])
Ttemp2 <- optim(par = c(log(8), log(0.3)), fn = Tmy_ll2, y = Tyi, method = "Nelder-Mead")
Ttemp2
Ttemp4 <- optim(par = c(log(8), log(0.3)), fn = Tmy_ll4, y = Tyi, method = "Nelder-Mead")
install.packages("irace")
library("irace")
q()
system.file(package = "irace")
install.packages("std")
install.packages("irace")
library("irace")
system.file(package = "irace")
library("irace")
require(Rcpp)
library("irace")
require(profvis)
library("irace")
library("irace")
lybrary("irace")
library("irace")
library("irace")
install.packages("irace")
install.packages("irace")
library("irace")
system.file(package = "irace")
install.packages("irace")
install.packages("irace")
library("irace")
system.file(package = "irace")
install.packages("irace"")
install.packages("irace")
install.packages("irace")
require(irace)
library("irace")
installed.packages("irace")
system.file(packages="irace")
system.file(package = "irace")