n <- c(100, 200, 500, 1000, 2000, 5000, 10000, 20000) # sample size
mu <- 1 # mean
sigma <- 1 # standard deviation
C <- 1 # constant
theta_true <- C - (sigma^2 + mu^2)
a <- 0.4 # stepsize constant
g <- function(X, theta) {
2 * mean(X^2 + theta - C)
}
f <- function(X, theta){
mean(X^2 + theta - C)^2
}Via the Robbins–Monro algorithm
rm_algo <- function(size) {
X <- rnorm(size, mu, sigma)
theta <- 1 # initial guess
temp <- Inf
k <- 0
while (abs(temp - theta) > 1e-7) {
k <- k + 1
a_k <- a/k
temp <- theta
theta <- theta - a_k * g(X, theta)
theta <- max(-3, min(theta, 1))
}
theta
}
part1_sa <- function(size) {
thetas <- replicate(50, rm_algo(size))
rmse <- sqrt(mean((thetas-theta_true)^2))
error <- sd(thetas)/sqrt(size)
rbind(rmse, error)
}
result1_sa <- sapply(n, part1_sa)
colnames(result1_sa) <- n
rownames(result1_sa) <- c("RMSE", "Standard Error")
result1_sa## 100 200 500 1000 2000
## RMSE 0.27796433 0.1915456 0.118876987 0.07459586 0.049090415
## Standard Error 0.02741368 0.0136813 0.005344907 0.00238286 0.001070532
## 5000 10000 20000
## RMSE 0.0329527543 0.0248970424 0.0189169760
## Standard Error 0.0004700945 0.0002514948 0.0001350654saa <- function(size) {
thetas <- seq(-3, 1, 0.05)
X <- rnorm(size, mu, sigma)
f_theta <- sapply(thetas, function(t){mean((X^2 + t - C)^2)})
thetas[which.min(f_theta)]
}
part1_saa <- function(size) {
thetas <- replicate(50, saa(size))
rmse <- sqrt(mean((thetas-theta_true)^2))
error <- sd(thetas)/sqrt(size)
rbind(rmse, error)
}
result1_saa <- sapply(n, part1_saa)
colnames(result1_saa) <- n
rownames(result1_saa) <- c("RMSE", "Standard Error")
result1_saa## 100 200 500 1000 2000
## RMSE 0.22616366 0.1565248 0.117898261 0.082158384 0.063245553
## Standard Error 0.02280821 0.0110887 0.005319199 0.002624259 0.001427857
## 5000 10000 20000
## RMSE 0.0393700394 0.0353553391 0.0212132
## Standard Error 0.0005607939 0.0003500729 0.0001500