使用 C# 实现对应分析

Miguel Diaz Kusztrich

4.43/5 (5投票s)

2016 年 10 月 31 日

CPOL

7分钟阅读

7936

242

用于添加到您自己的统计类库中的对应分析的完整算法

下载 CADemo 项目的源代码 - 8.4 MB

引言

在本文中，我将向您展示如何执行与统计分析技术（称为对应分析）相关的所有计算。

本文提供的代码是一个 Windows Forms 应用程序，尽管通常您会希望将代码与一些其他关于数学计算的类一起放入类库中。在这种情况下，这样做非常容易，因为几乎所有提供的代码都涉及对应分析计算，所以这只是一个复制粘贴的工作。

为了可读性，代码分为三个步骤，尽管您可以将其全部放在一个函数中。第一步是从 csv 文件加载数据并计算文件和列的质量。接下来，执行计算以获取卡方距离和惯量。最后，将数据投影到二维子空间中，以便在图形图表中绘制。

除了数学代码之外，还提供了一个简单的图形绘制示例。

背景

我假设您知道什么是对应分析。如果不是，您可以访问我的博客（http://software-tecnico-libre.es/en/article-by-topic/software-developing/charp-programming/mathematics-with-csharp/correspondence-analysis-csharp），了解与本文代码相关的该技术的简要说明。

它基本上包括对列联表中的数据进行最优缩放，以便您可以像连续数据一样查看分类数据之间的关系，从而生成图形地图，其中表示不同的行和列，并且您可以发现它们之间一些不明显的关联。

使用代码

首先，让我们回顾一下代码中使用的全局数据。主变量是 _data，它是一个 DataTable，将包含列联表，以及对应分析的行和列计算以及用于在二维中绘制图表的主坐标和标准坐标。

private DataTable _data = null;

列联表的大小在以下变量中

private int _colCount = 0;
private int _rowCount = 0;

为了提高可读性，我定义了两个枚举来标识包含计算值的不同列和行。

public enum ColumnPosition
{
	RowNames, FirstDataColumn, Mass, ChiDist, Inertia, Dim1RowStandard, Dim2RowStandard, Dim1RowPrincipal, Dim2RowPrincipal
}

public enum RowPosition
{
	FirstDataRow, Mass, ChiDist, Inertia, Dim1ColStandard, Dim2ColStandard, Dim1ColPrincipal, Dim2ColPrincipal
}

行和列中的数据类型为 double，除了行名（第一列）为字符串类型。初始列联表由绝对频率组成，值为 0 到 1 之间。

添加到表中的行和列（对应分析对行和列是对称的）用于质量或边缘频率、到质心的卡方距离、惯量以及每行和每列的标准坐标和主坐标。

使用这两个辅助函数，您可以独立于列联表的大小来定位这些行/列中的每一个

private int GetColPosition(ColumnPosition c)
{
	switch (c)
	{
 		case ColumnPosition.RowNames:
                    return 0;
                case ColumnPosition.FirstDataColumn:
                    return 1;
                case ColumnPosition.Mass:
                    return 1 + _colCount;
                case ColumnPosition.ChiDist:
                    return 2 + _colCount;
                case ColumnPosition.Inertia:
                    return 3 + _colCount;
                case ColumnPosition.Dim1RowStandard:
                    return 4 + _colCount;
                case ColumnPosition.Dim2RowStandard:
                    return 5 + _colCount;
                case ColumnPosition.Dim1RowPrincipal:
                    return 6 + _colCount;
                case ColumnPosition.Dim2RowPrincipal:
                    return 7 + _colCount;
        }
	throw new ArgumentException();
}

private int GetRowPosition(RowPosition r)
{
	switch (r)
        {
                case RowPosition.FirstDataRow:
                    return 0;
                case RowPosition.Mass:
                    return _rowCount;
                case RowPosition.ChiDist:
                    return 1 + _rowCount;
                case RowPosition.Inertia:
                    return 2 + _rowCount;
                case RowPosition.Dim1ColStandard:
                    return 3 + _rowCount;
                case RowPosition.Dim2ColStandard:
                    return 4 + _rowCount;
                case RowPosition.Dim1ColPrincipal:
                    return 5 + _rowCount;
                case RowPosition.Dim2ColPrincipal:
                    return 6 + _rowCount;
	}
        throw new ArgumentException();
}

程序界面如下所示

CADemo application

数据文件必须是 csv 格式，第一行包含列名，第一列包含行名，如下图所示

CSV format of the data tables

您可以提供用于分隔列和小数位数的字符，以及用于对结果值进行四舍五入的精度，该精度存储在 _precision 全局变量中。

private int _precision = 6;

工具栏左侧的按钮允许选择并加载一个 csv 数据文件，这是算法的第一步。随着 _data DataTable 的初始化，边际频率被计算并存储在行和列的质量中。

由于在此步骤中计算了质量，因此还初始化了用于缩放图中点的变量，这些变量是

private float _rMassSc = 0f;
private float _cMassSc = 0f;
private float _minRMass = 0f;
private float _minCMass = 0f;

我使用浮点值，因为这些变量将用于图形计算。

对于按质量缩放，还有一个常数，它定义了点将具有的最小相对大小，以避免打印过小的点。

private const float _cMinScale = 0.5f;

代码位于按钮的 Click 事件中。首先，初始化用户界面和局部变量

StreamReader rdr = new StreamReader(OFDlg.FileName);
bsData.DataSource = null;
bChid.Enabled = false;
bProject2D.Enabled = false;
tbWizzard.TabPages.Remove(tpGraph);
int.TryParse(txtPrec.Text, out _precision);
if ((_precision <= 0) || (_precision > 16))
{
	_precision = 6;
	txtPrec.Text = "6";
}

然后，创建数据列和质量列

string line = rdr.ReadLine();
...
string[] columns = line.Split(txtSep.Text[0]);
_colCount = columns.Length - 1;
_rowCount = 0;
_data = new DataTable();
// First Column for Row Names
_data.Columns.Add(new DataColumn("row/col", typeof(string)));
for (int ix = 1; ix < columns.Length; ix++)
{
	_data.Columns.Add(new DataColumn(columns[ix].Replace("\"", ""), typeof(double)));
}
// Add Row Mass Column
_data.Columns.Add(new DataColumn("Mass", typeof(double)));

接下来，从文件中读取数据并计算行质量和比例

DataRow row;
double total = 0;
// Min and Max Mass value
double max = double.MinValue;
double min = double.MaxValue;
while ((line = rdr.ReadLine()) != null)
{
        // Add Rows
        row = _data.NewRow();
        string[] rowdata = line.Split(txtSep.Text[0]);
        // First Element is the Name
        row[GetColPosition(ColumnPosition.RowNames)] = rowdata[0].Replace("\"", "");
        double rtotal = 0;
        for (int ix = GetColPosition(ColumnPosition.FirstDataColumn); ix < rowdata.Length; ix++)
        {
                // Add Values
                double v = Convert.ToDouble(rowdata[ix].Replace(txtDecSep.Text, CultureInfo.CurrentCulture.NumberFormat.NumberDecimalSeparator));
                total += v;
                rtotal += v;
                row[ix] = Math.Round(v, _precision);
        }
        row[GetColPosition(ColumnPosition.Mass)] = Math.Round(rtotal, _precision);
        max = Math.Max(rtotal, max);
        min = Math.Min(rtotal, min);
        _data.Rows.Add(row);
        _rowCount++;
}
// Row Mass scale
_minRMass = (float)min;
_rMassSc = (float)((1 - _cMinScale) / (max - min));

最后计算列质量和列的质量比例。

max = double.MinValue;
min = double.MaxValue;
// Add Column Mass Row
row = _data.NewRow();
row[0] = "Mass";                       
for (int ix = 1; ix < _data.Columns.Count - 1; ix++)
{
        double ctotal = 0;
        for (int ir = 0; ir < _data.Rows.Count; ir++)
        {
                ctotal += (double)_data.Rows[ir][ix];
        }
        row[ix] = Math.Round(ctotal, _precision);
        max = Math.Max(ctotal, max);
        min = Math.Min(ctotal, min);
}
row[GetColPosition(ColumnPosition.Mass)] = Math.Round(total, _precision);
_data.Rows.Add(row);
// Column Mass scale
_minCMass = (float)min;
_cMassSc = (float)((1 - _cMinScale) / (max - min));

下一步是按下卡方距离按钮时执行的，它包括计算行和列与它们的质心之间的卡方距离以及它们的派生度量（称为惯量），该度量类似于方差。

代码在按钮的 Click 事件中。

首先，对于每个数据行，计算每个数据列的实际比例值（通过将绝对频率除以边际频率获得）与质量值（质心的频率）之间的差的平方和，再除以质量值。

此和的平方根是该行的卡方距离，将此值乘以质量，即可得到惯量。

// Total inertia
double inertia = 0;
_data.Columns.Add("Chi-Dist", typeof(double));
_data.Columns.Add("Inertia", typeof(double));
// Calculate for the rows
for (int ir = GetRowPosition(RowPosition.FirstDataRow); ir < GetRowPosition(RowPosition.Mass); ir++)
{
        double chdrow = 0;
        for (int ic = GetColPosition(ColumnPosition.FirstDataColumn); ic < GetColPosition(ColumnPosition.Mass); ic++)
        {
                // (observed (row profile) - expected)^2 / expected
                chdrow += Math.Pow(((double)_data.Rows[ir][ic] / (double)_data.Rows[ir][GetColPosition(ColumnPosition.Mass)]) - (double)_data.Rows[GetRowPosition(RowPosition.Mass)][ic], 2) / (double)_data.Rows[GetRowPosition(RowPosition.Mass)][ic];
        }
        // Chi-distance
        _data.Rows[ir][GetColPosition(ColumnPosition.ChiDist)] = Math.Round(Math.Sqrt(chdrow), _precision);
        // Inertia
        inertia += chdrow * (double)_data.Rows[ir][GetColPosition(ColumnPosition.Mass)];
        _data.Rows[ir][GetColPosition(ColumnPosition.Inertia)] = Math.Round(chdrow * (double)_data.Rows[ir][GetColPosition(ColumnPosition.Mass)], _precision);
}

然后，进行相同的计算以添加两行，其中包含列的卡方距离和惯量。

表的总惯量是行惯量（或列惯量，它们相等）的总和。

DataRow rowc = _data.NewRow();
rowc[0] = "Chi-Dist";
DataRow rowi = _data.NewRow();
rowi[0] = "Inertia";
// Calculate for the columns
for (int ic = GetColPosition(ColumnPosition.FirstDataColumn); ic < GetColPosition(ColumnPosition.Mass); ic++)
{
        double chdcol = 0;
        for (int ir = GetRowPosition(RowPosition.FirstDataRow); ir < GetRowPosition(RowPosition.Mass); ir++)
        {
                // (observed (column profile) - expected)^2 / expected
                chdcol += Math.Pow(((double)_data.Rows[ir][ic] / (double)_data.Rows[GetRowPosition(RowPosition.Mass)][ic]) - (double)_data.Rows[ir][GetColPosition(ColumnPosition.Mass)], 2) / (double)_data.Rows[ir][GetColPosition(ColumnPosition.Mass)];
        }
        // Chi-distance
        rowc[ic] = Math.Round(Math.Sqrt(chdcol), _precision);
        // Inertia
        rowi[ic] = Math.Round(chdcol * (double)_data.Rows[GetRowPosition(RowPosition.Mass)][ic], _precision);
}
// Total inertia
rowi[GetColPosition(ColumnPosition.Inertia)] = Math.Round(inertia, _precision);
_data.Rows.Add(rowc);
_data.Rows.Add(rowi);

对应分析中最复杂的计算是奇异值分解 (SVD)，用于将 n 维剖面空间投影到二维子空间中，以便绘制图形，获得每行和每列的标准坐标和主坐标。

也许您有自己的库来计算此分解，我使用了 NuGet 包仓库中的 Accord.Math 库中的 SingularValueDecomposition 类。您可以轻松更改代码以使用您喜欢的软件。

这些计算在“2D 投影”按钮的 Click 事件中执行。

首先，计算 S 矩阵，在 wdata 局部变量中，公式为 S = D_r^-1/2(P - rc^T)D_c^-1/2，其中 D_r^-1/2 是一个对角矩阵，其对角线由行质量的平方根的倒数组成。D_c^-1/2 对于列相同，P 是列联表，r 是具有行质量的对角矩阵，c 是具有列质量的对角矩阵。

double[,] wdata = new double[_rowCount, _colCount];
for (int ir = GetRowPosition(RowPosition.FirstDataRow); ir < GetRowPosition(RowPosition.Mass); ir++)
{
        for (int ic = GetColPosition(ColumnPosition.FirstDataColumn); ic < GetColPosition(ColumnPosition.Mass); ic++)
        {
                wdata[ir, ic - 1] = Math.Round(((double)_data.Rows[ir][ic] - (double)_data.Rows[ir][GetColPosition(ColumnPosition.Mass)] * (double)_data.Rows[GetRowPosition(RowPosition.Mass)][ic]) /
                        Math.Sqrt((double)_data.Rows[ir][GetColPosition(ColumnPosition.Mass)] * (double)_data.Rows[GetRowPosition(RowPosition.Mass)][ic]), _precision);
        }
}

SVD 将此矩阵分解为另外三个，S = UDV^T，其中 U 是左奇异向量，D 是一个对角矩阵，其对角线包含降序排列的奇异值，V 是一个包含右奇异向量的矩阵。

SingularValueDecomposition svd = new SingularValueDecomposition(wdata);

从这里，您可以从 U 矩阵的前两个左向量除以质量的平方根，得到每个行在二维中的标准坐标。

主坐标的计算方法是将标准坐标乘以最初的两个奇异值。

// Standard coordinates
_data.Columns.Add("Dim1 S.", typeof(double));
_data.Columns.Add("Dim2 S.", typeof(double));
// Principal coordinates
_data.Columns.Add("Dim1 P.", typeof(double));
_data.Columns.Add("Dim2 P.", typeof(double));
// Processing rows
float[,] vmax = new float[2, 4];
for (int ir = GetRowPosition(RowPosition.FirstDataRow); ir < GetRowPosition(RowPosition.Mass); ir++)
{
        for (int iu = 0; iu < 2; iu++)
        {
                // Project using U vectors of SVD divided by the sqrt of mass
                double v = svd.LeftSingularVectors[ir, iu] / Math.Sqrt((double)_data.Rows[ir][GetColPosition(ColumnPosition.Mass)]);
                // Standard coordinates
                _data.Rows[ir][iu + GetColPosition(ColumnPosition.Dim1RowStandard)] = Math.Round(v, _precision);
                // Principal coordinates using singular value
                _data.Rows[ir][iu + GetColPosition(ColumnPosition.Dim1RowPrincipal)] = Math.Round(v * svd.Diagonal[iu], _precision);
                // Store max and min coordinates for each row & dimension
                vmax[iu, 0] = (float)Math.Min(vmax[iu, 0], Math.Round(v, _precision));
                vmax[iu, 1] = (float)Math.Max(vmax[iu, 1], Math.Round(v, _precision));
                vmax[iu, 2] = (float)Math.Min(vmax[iu, 2], Math.Round(v * svd.Diagonal[iu], _precision));
                vmax[iu, 3] = (float)Math.Max(vmax[iu, 3], Math.Round(v * svd.Diagonal[iu], _precision));
        }
}
_rowDimS1 = new SizeF(vmax[0, 0], vmax[0, 1]);
_rowDimS2 = new SizeF(vmax[1, 0], vmax[1, 1]);
_rowDimP1 = new SizeF(vmax[0, 2], vmax[0, 3]);
_rowDimP2 = new SizeF(vmax[1, 2], vmax[1, 3]);

变量 _rowDimS1、_rowDimS2、_rowDimP1 和 _rowDimP2 是绘制图形的辅助变量。它们包含每个维度的最大和最小坐标。

接下来，执行相同的计算以获得列的坐标，但在这种情况下使用来自 V 矩阵的右奇异向量。

// Same for columns
DataRow row1s = _data.NewRow();
row1s[0] = "Dim1 S.";
DataRow row2s = _data.NewRow();
row2s[0] = "Dim2 S.";
DataRow row1p = _data.NewRow();
row1p[0] = "Dim1 P.";
DataRow row2p = _data.NewRow();
row2p[0] = "Dim2 P.";
vmax = new float[2, 4];
for (int ic = GetColPosition(ColumnPosition.FirstDataColumn); ic < GetColPosition(ColumnPosition.Mass); ic++)
{
        double v1 = svd.RightSingularVectors[ic - 1, 0] / Math.Sqrt((double)_data.Rows[GetRowPosition(RowPosition.Mass)][ic]);
        row1s[ic] = Math.Round(v1, _precision);
        row1p[ic] = Math.Round(v1 * svd.Diagonal[0], _precision);
        vmax[0, 0] = (float)Math.Min(vmax[0, 0], Math.Round(v1, _precision));
        vmax[0, 1] = (float)Math.Max(vmax[0, 1], Math.Round(v1, _precision));
        vmax[0, 2] = (float)Math.Min(vmax[0, 2], Math.Round(v1 * svd.Diagonal[0], _precision));
        vmax[0, 3] = (float)Math.Max(vmax[0, 3], Math.Round(v1 * svd.Diagonal[0], _precision));
        double v2 = svd.RightSingularVectors[ic - 1, 1] / Math.Sqrt((double)_data.Rows[GetRowPosition(RowPosition.Mass)][ic]);
        row2s[ic] = Math.Round(v2, _precision);
        row2p[ic] = Math.Round(v2 * svd.Diagonal[1], _precision);
        vmax[1, 0] = (float)Math.Min(vmax[1, 0], Math.Round(v2, _precision));
        vmax[1, 1] = (float)Math.Max(vmax[1, 1], Math.Round(v2, _precision));
        vmax[1, 2] = (float)Math.Min(vmax[1, 2], Math.Round(v2 * svd.Diagonal[1], _precision));
        vmax[1, 3] = (float)Math.Max(vmax[1, 3], Math.Round(v2 * svd.Diagonal[1], _precision));
}
_colDimS1 = new SizeF(vmax[0, 0], vmax[0, 1]);
_colDimS2 = new SizeF(vmax[1, 0], vmax[1, 1]);
_colDimP1 = new SizeF(vmax[0, 2], vmax[0, 3]);
_colDimP2 = new SizeF(vmax[1, 2], vmax[1, 3]);
_data.Rows.Add(row1s);
_data.Rows.Add(row2s);
_data.Rows.Add(row1p);
_data.Rows.Add(row2p);

最后的计算是将表的总惯量投影到二维空间。总惯量与投影惯量之和之间的差值可以让你了解维度缩减所造成的惯量损失。

// Calculate total inertias for each dimension (row standard coordinates * contingency table * column standard coordinates)
for (int iv = 0; iv < 2; iv++)
{
        double itotal = 0;
        for (int ic = GetColPosition(ColumnPosition.FirstDataColumn); ic < GetColPosition(ColumnPosition.Mass); ic++)
        {
                double rtmp = 0;
                for (int ir = 0; ir < _data.Rows.Count - 7; ir++)
                {
                        rtmp += (double)_data.Rows[ir][ic] * (double)_data.Rows[ir][iv + GetColPosition(ColumnPosition.Dim1RowStandard)];
                }
                itotal += rtmp * (double)_data.Rows[iv + GetRowPosition(RowPosition.Dim1ColStandard)][ic];
        }
        // Put in the inertia row and column
        _data.Rows[GetRowPosition(RowPosition.Inertia)][iv + GetColPosition(ColumnPosition.Dim1RowStandard)] = Math.Round(itotal * itotal, _precision);
        _data.Rows[iv + GetRowPosition(RowPosition.Dim1ColStandard)][GetColPosition(ColumnPosition.Inertia)] = Math.Round(itotal * itotal, _precision);
}

当所有计算完成后，程序会显示一个新选项卡，您可以在其中看到对应分析的图形，如下图所示

The graph of the correspondence analysis

您可以使用工具栏按钮显示使用列、行或两者的主坐标的数据，采用对称版本。

使用“质量”按钮，您可以按与点质量成比例的大小显示这些点。

通过 SetGraphScales，根据所选选项从辅助变量中选择图形边界。

private void SetGraphScales()
{
        // Min and Max coordinates
        if (bSymm.Checked)
        {
                _sizeX = new SizeF(Math.Min(_rowDimP1.Width, _colDimP1.Width), Math.Max(_rowDimP1.Height, _colDimP1.Height));
                _sizeY = new SizeF(Math.Min(_rowDimP2.Width, _colDimP2.Width), Math.Max(_rowDimP2.Height, _colDimP2.Height));
        }
        else if (bRow.Checked)
        {
                _sizeX = new SizeF(Math.Min(_rowDimP1.Width, _colDimS1.Width), Math.Max(_rowDimP1.Height, _colDimS1.Height));
                _sizeY = new SizeF(Math.Min(_rowDimP2.Width, _colDimS2.Width), Math.Max(_rowDimP2.Height, _colDimS2.Height));
        }
        else
        {
                _sizeX = new SizeF(Math.Min(_rowDimS1.Width, _colDimP1.Width), Math.Max(_rowDimS1.Height, _colDimP1.Height));
                _sizeY = new SizeF(Math.Min(_rowDimS2.Width, _colDimP2.Width), Math.Max(_rowDimS2.Height, _colDimP2.Height));
        }
}

还有一些辅助函数，可以根据当前设置获取行和列坐标以及质量。

private float RowCoordinate(int dim, int row)
{
        if (bRow.Checked || bSymm.Checked)
        {
                // Rows in principal coordinates
                return dim == 0 ? Convert.ToSingle(_data.Rows[row][GetColPosition(ColumnPosition.Dim1RowPrincipal)]) : Convert.ToSingle(_data.Rows[row][GetColPosition(ColumnPosition.Dim2RowPrincipal)]);
        }
        else
        {
                // Rows in standard coordinates
                return dim == 0 ? Convert.ToSingle(_data.Rows[row][GetColPosition(ColumnPosition.Dim1RowStandard)]) : Convert.ToSingle(_data.Rows[row][GetColPosition(ColumnPosition.Dim2RowStandard)]);
        }
}
private float ColCoordinate(int dim, int col)
{
        if (bCol.Checked || bSymm.Checked)
        {
                // Columns in principal coordiantes
                return dim == 0 ? Convert.ToSingle(_data.Rows[GetRowPosition(RowPosition.Dim1ColPrincipal)][col]) : Convert.ToSingle(_data.Rows[GetRowPosition(RowPosition.Dim2ColPrincipal)][col]);
        }
        else
        {
                // Columns in standard coordinates
                return dim == 0 ? Convert.ToSingle(_data.Rows[GetRowPosition(RowPosition.Dim1ColStandard)][col]) : Convert.ToSingle(_data.Rows[GetRowPosition(RowPosition.Dim2ColStandard)][col]);
        }
}
private float RowMass(int row, float max)
{
        return max * (bMass.Checked ? _cMinScale + (Convert.ToSingle(_data.Rows[row][GetColPosition(ColumnPosition.Mass)]) - _minRMass) * _rMassSc : 1f);
}
private float ColMass(int col, float max)
{
        return max * (bMass.Checked ? _cMinScale + (Convert.ToSingle(_data.Rows[GetRowPosition(RowPosition.Mass)][col]) - _minCMass) * _cMassSc : 1f);
}

绘制图的代码位于“Graph”选项卡中面板的 Paint 事件中。它并不复杂，我认为它是不言自明的。它只是一个示例，可以改进以将其包含在您自己的图形库中。

这就是示例项目所包含的全部内容。我使用 Visual Studio 2013 版和 .NET Framework 4.5 版编写了它，但我认为可以使用任何其他版本生成它。

感谢阅读，希望您觉得这些代码有用。