for (s in 1:S) { bio_sims <- rnorm(K, theta_new_samples_pop[s], sigma_bio_samples[s, 1]) weights_k <- numeric(K) for (k in 1:K) { weights_k[k] <- exp(log(dt((y_new_1 - bio_sims[k]) / sigma_X_samples[s, 1], df = df_samples[s, 1] ) / sigma_X_samples[s, 1]) + log(dt((y_new_2 - bio_sims[k]) / sigma_X_samples[s, 1], df = df_samples[s, 1] ) / sigma_X_samples[s, 1])) } weights_steady[s] <- mean(weights_k, na.rm = T) }
something like: bio_sims <- rnorm((K, N), theta_new_samples, sigma_bio_samples) to sample a matrix of (K,N) normal random values and then plug that into dt(...) and then take the mean over the K axis of the matrix. Or you could at least vectorize the inner loop by first creating a vector of v_diff <- (y_new_1 - bio_sims)/sigma_X_samples and then do exp(log(dt(v_diff))... etc