C++PCA-计算协方差矩阵

Question

我正在尝试计算C++中向量的协方差（矩阵）。。。

我进行了以下工作：

std::vector<std::vector<double> > data = { {2.5, 2.4}, {0.5, 0.7} };

然后我计算并减去平均值，得到以下结果：

data = { {0.05, -0.05}, {-0.1, 0.1} }

据我所知，下一步是转置矩阵，将原点相乘，求和，最后除以X-1

我写了以下内容：

void cover(std::vector<std::vector<double> > &d)
{
    double cov = 0.0; 
    for(unsigned i=0; (i < d.size()); i++)
    {
        for(unsigned j=0; (j < d[i].size()); j++)
        {
            cov += d[i][j] * d[j][i] / (d[i].size() - 1);
            std::cout << cov << " ";
        }
        std::cout << std::endl;
    }
}

其中d是从每个点减去平均值后的矢量这给了我一个结果：

0.0025, 0.0075 
0.0125, 0.0225

与matlab相比：

2.0000    1.7000
1.7000    1.4450

有人知道我哪里错了吗？

感谢

Answer 1

此语句：

据我所知，下一步是转置矩阵，将原点相乘，求和，最后除以X-1

这个实现：

cov += d[i][j] * d[j][i] / (d[i].size() - 1);

不要说同样的话。基于此处的定义：

void outer_product(vector<double> row, vector<double> col, vector<vector<double>>& dst) {
    for(unsigned i = 0; i < row.size(); i++) {
        for(unsigned j = 0; j < col.size(); i++) {
            dst[i][j] = row[i] * col[j];
        }
    }
}
//computes row[i] - val for all i;
void subtract(vector<double> row, double val, vector<double>& dst) {
    for(unsigned i = 0; i < row.size(); i++) {
        dst[i] = row[i] - val;
    }
}
//computes m[i][j] + m2[i][j]
void add(vector<vector<double>> m, vector<vector<double>> m2, vector<vector<double>>& dst) {
    for(unsigned i = 0; i < m.size(); i++) {
        for(unsigned j = 0; j < m[i].size(); j++) { 
            dst[i][j] = m[i][j] + m2[i][j];
        }
    }
}
double mean(std::vector<double> &data) {
    double mean = 0.0;
    for(unsigned i=0; (i < data.size());i++) {
        mean += data[i];
    }
    mean /= data.size();
    return mean;
}
void scale(vector<vector<double>> & d, double alpha) {
    for(unsigned i = 0; i < d.size(); i++) {
        for(unsigned j = 0; j < d[i].size(); j++) {
            d[i][j] *= alpha;
        }
    }
}

因此，给定这些定义，我们可以计算协方差矩阵的值。

void compute_covariance_matrix(vector<vector<double>> & d, vector<vector<double>> & dst) {
    for(unsigned i = 0; i < d.size(); i++) {
        double y_bar = mean(d[i]);
        vector<double> d_d_bar(d[i].size());
        subtract(d[i], y_bar, d_d_bar);
        vector<vector<double>> t(d.size());
        outer_product(d_d_bar, d_d_bar, t);
        add(dst, t, dst);
    } 
    scale(dst, 1/(d.size() - 1));
}

Answer 2

我认为对于outer_product中的循环，这是错误的：

void outer_product(vector<double> row, vector<double> col, vector<vector<double>>& dst) {
for(unsigned i = 0; i < row.size(); i++) {
    for(unsigned j = 0; j < col.size(); i++) {
        dst[i][j] = row[i] * col[j];
    }
}

我将更改I++->j++