Revision 34c2b81061aaf8d77873c5431608970613444ea9 authored by Christiana Kartsonaki on 06 June 2022, 22:10:05 UTC, committed by cran-robot on 06 June 2022, 22:10:05 UTC
1 parent 42f5443
fit.model.R
``````fit.model <-
function (p, q, n, r, starting_values, h_vector, data_true, sim_data,
features, n_iter, print_results = TRUE, variances = TRUE)
{
h <- list()
for (s in 1:p) {
h[[s]] <- h_vector[s]
}

theta_0 <- list()
for (s in 1:p) {
theta_0[[s]] <- starting_values[s]
}

k <- 2^ceiling(log2(p + 1))

# library(survey)
if (p == 1) {
D <- matrix(c(1, -1), nrow = 2, ncol = 1)
}
if (p > 1) {
D <- 2 * survey::hadamard(k - 1) - 1
D <- D[, 2:(p + 1)]
}

sim <- list()
z <- list()

results <- matrix(NA, nrow = n_iter, ncol = p)
k_ones <- matrix(1, nrow = k, ncol = 1)
r_ones <- matrix(1, nrow = r, ncol = 1)
design.points <- list()
z_data <- matrix(NA, nrow = 1, ncol = q)

for (j in 1:q) {
z_data[1, j] <- features(data_true)[j]
}
z_data <- t(z_data)

if (variances == TRUE) {
var_theta <- list()
}

for (iter in 1:n_iter) {

for (d in 1:k) {
sim[[d]] <- list()
z[[d]] <- matrix(NA, nrow = r, ncol = q)
}

for (d in 1:k) {
design.points[[d]] <- list()
for (s in 1:p) {
design.points[[d]][[s]] <- theta_0[[s]] + (D[d, s] * h[[s]])
}
for (i in 1:r) {
sim[[d]][[i]] <- sim_data(n, design.points[[d]])
}
}

for (d in 1:k) {
for (i in 1:r) {
z[[d]][i, ] <- features(sim[[d]][[i]])
}
}

zbar <- list()
for (d in 1:k) {
zbar[[d]] <- (t(r_ones) %*% z[[d]])/r
}

Zbar <- matrix(unlist(zbar), nrow = q, byrow = F)

ztilde <- list()
sigma <- list()
for (d in 1:k) {
ztilde[[d]] <- z[[d]] - (r_ones %*% zbar[[d]])
sigma[[d]] <- (t(ztilde[[d]]) %*% ztilde[[d]])/(r - 1)
}

V <- Zbar %*% D

sigma_vv <- (1/k) * Reduce("+", sigma)
sigma_inv <- solve(sigma_vv)

lambda <- rep(NA, times = p)
for (s in 1:p) {
lambda[s] <- sqrt((t(V[, s]) %*% sigma_inv %*% V[, s]))
}

L <- (1/lambda) * (sigma_inv %*% V)
Ybar <- t(L) %*% Zbar

Y_data <- t(L) %*% z_data

A <- (2/k) * (t(V) %*% L)
for (u in 1:dim(A)[1]) {
for (v in 1:dim(A)[2]) {
A[u, v][u != v] <- 0
}
}

H <- matrix(0, nrow = p, ncol = p)
for (s in 1:p) {
H[s, s] <- h[[s]]
}

theta_hat <- matrix(NA, nrow = p, ncol = 1)
th_0 <- matrix(unlist(theta_0))

theta_hat <- th_0 + (2 * H %*% solve(A) %*% (Y_data - ((Ybar %*% k_ones)/k)))

for (s in 1:p) {
theta_0[[s]] <- theta_hat[s]
}

results[iter, ] <- t(theta_hat)

if (variances == TRUE) {
var_theta[[iter]] <- 4 * H %*% solve(A) %*% t(L) %*% sigma_vv %*% L %*% H %*% solve(A)
}

if (print_results == TRUE) {
print(results[iter, ])
}
}

if (variances == FALSE) {
var_theta <- 4 * H %*% solve(A) %*% t(L) %*% sigma_vv %*% L %*% H %*% solve(A)
}

return(list(estimates = results, var_estimates = var_theta, L = L, sigma = sigma_vv, zbar = Zbar, z_D = z_data, ybar = Ybar, y_D = Y_data))

}

``````

Computing file changes ...