Сделать вложенные циклы более эффективными?

Я анализирую большие наборы данных, используя следующий сценарий:

M <- c_alignment 
c_check <- function(x){
    if (x == c_1) {
        1
    }else{
        0
    }
}
both_c_check <- function(x){
    if (x[res_1] == c_1 && x[res_2] == c_1) {
        1
    }else{
        0
    }
}
variance_function <- function(x,y){
    sqrt(x*(1-x))*sqrt(y*(1-y))
}
frames_total <- nrow(M)
cols <- ncol(M)
c_vector <- apply(M, 2, max)
freq_vector <- matrix(nrow = sum(c_vector))
co_freq_matrix <- matrix(nrow = sum(c_vector), ncol = sum(c_vector))
insertion <- 0
res_1_insertion <- 0
for (res_1 in 1:cols){
    for (c_1 in 1:conf_vector[res_1]){
        res_1_insertion <- res_1_insertion + 1
        insertion <- insertion + 1
        res_1_subset <- sapply(M[,res_1], c_check)
        freq_vector[insertion] <- sum(res_1_subset)/frames_total
        res_2_insertion <- 0
        for (res_2 in 1:cols){
            if (is.na(co_freq_matrix[res_1_insertion, res_2_insertion + 1])){
                for (c_2 in 1:max(c_vector[res_2])){
                    res_2_insertion <- res_2_insertion + 1
                    both_res_subset <- apply(M, 1, both_c_check)
                    co_freq_matrix[res_1_insertion, res_2_insertion] <- sum(both_res_subset)/frames_total
                    co_freq_matrix[res_2_insertion, res_1_insertion] <- sum(both_res_subset)/frames_total
                }
            }
        }
    }
}
covariance_matrix <- (co_freq_matrix - crossprod(t(freq_vector)))
variance_matrix <- matrix(outer(freq_vector, freq_vector, variance_function), ncol = length(freq_vector))
correlation_coefficient_matrix <- covariance_matrix/variance_matrix

Входные данные модели будут примерно такими:

1 2 1 4 3
1 3 4 2 1
2 3 3 3 1
1 1 2 1 2
2 3 4 4 2

Я вычисляю биномиальную ковариацию для каждого состояния, найденного в M [, i] с каждым состоянием, найденным в M [, j] . Каждая строка - это состояние, найденное для этого испытания, и я хочу увидеть, как состояния столбцов взаимно меняются.

Уточнение: я нахожу ковариацию двух полиномиальных распределений, но делаю это посредством биномиальных сравнений.

Входными данными является матрица 4200 x 510, а значение c для каждого столбца в среднем составляет около 15. Я знаю, что циклы для ужасно медленные в R, но я не уверен, как здесь можно использовать функцию apply . Если у кого-то есть предложения по правильному использованию здесь apply , я буду очень признателен. Сейчас сценарий занимает несколько часов. Спасибо!