项目主页：randomforest github, gitee
C++ implementation of random forests classification, regression, proximity and variable importance.
推荐阅读：Random Forests C++实现：细节，使用与实验

随机森林和数据插补简介

随机森林（Random Forest, RF）作为一种强大的机器学习算法，已被广泛应用于数据（缺失值）插补（imputation）领域。数据插补是指在处理缺失值时，通过某种方法估计这些缺失值的过程。随机森林在数据插补中的应用主要基于其能够处理分类和回归问题的能力，以及其对缺失数据的鲁棒性。
随机森林通过构建多个决策树并进行集成学习来提高预测的准确性和稳定性。每个决策树都是在不同的样本集合上训练的，这些样本集合是通过对原始数据进行有放回的随机抽样得到的。这种方法的一个关键优势是它能够自动识别和利用数据中的非线性关系，从而在处理具有复杂结构的数据时表现出色。
将随机森林应用于数据插补能获得一些令人满意的特性，比如它能适用于不同类型的缺失数据，能够适应非线性数据，也可以适用于较大的数据集。对于随机森林数据插补算法比较完整的分类介绍，以及对应的使用指导，推荐以下综述型论文。

Tang F, Ishwaran H. Random Forest Missing Data Algorithms. Stat Anal Data Min: The ASA Data Sci Journal. 2017; 10: 363–377.

基于随机森林Proximity的缺失数据插补

原理

基于随机森林Proximity的插补算法不仅适用于连续型(continuous /numerical)特征，也适用于类别型(categorical)或离散型特征。算法核心思想是使用随机森林Proximity来计算缺失样本与非缺失样本之间的相似度，然后通过非缺失样本的加权值来估计缺失值。
现有训练集 { x i , y i } i = 1 N \left \{ \mathbf{x} _i,y_i \right \} _{i=1}^N {xi​,yi​}i=1N​，样本${\mathbf{x}}_i$对应类标签或者目标值$y_i$，两个样本之间的相似度为$p(i,j)$。对于连续型变量，样本某维度的缺失值($x_i^d=nan$)通过其他非缺失样本该维度的加权平均计算得到，这里的权值就是样本间的相似度Proximity (the proximity weighted average of non-missing data):

x ˉ i d = ∑ { j : x j d ≠ n a n } p ( i , j ) I ( y i = y j ) x j d ∑ { j : x j d ≠ n a n } p ( i , j ) I ( y i = y j ) \bar{x}_i^d =\frac{ \sum_{\{j:x_j^d \ne nan\} } p\left ( i,j\right ) I\left (y_i=y_j \right )x_j^d}{\sum_{\{j: x_j^d \ne nan\}}p(i,j)I(y_i=y_j)} xˉid​=∑{j:xjd​=nan}​p(i,j)I(yi​=yj​)∑{j:xjd​=nan}​p(i,j)I(yi​=yj​)xjd​​

对于无序类别型(categorical)变量$x_i^d\in S_d$（$S_d$为有限的值集，比如 {0, 1}），缺失值将被插补为该维度上非缺失样本中出现“频度”最高的值，“频度”为累计Proximity (the most frequent non-missing value where frequency is weighted by proximity)：
x ˉ i d = arg ⁡ max ⁡   c ∈ S d ∑ { j : x j d = c } p ( i , j ) I ( y i = y j ) \bar{x}_i^d ={ \underset{c\in S_d}{{\arg\max} \, }}\sum_{ \{ j:x_j^d=c\}} p\left ( i,j \right )I\left ( y_i=y_j \right ) xˉid​=c∈Sd​argmax​{j:xjd​=c}∑​p(i,j)I(yi​=yj​)
说明：对于分类问题，上面的公式还考虑了$x_i$与$x_j$的类别需要相同；对于回归问题，则无需考虑类别。
本文采用RF-GAP （TPAMI2023)）算法计算训练集样本之间相似性度量Proximities。RF-GAP利用out-of-bag样本与in-bag样本在叶子节点中的关系计算得到$P(i,j)$，反映了通过随机森林训练得到了数据集几何结构(Geometry)。根据作者的实验，RF-GAP优于原始的随机森林Proximity方法。

算法流程

步骤如下：
1 . 预处理：如果是分类问题，使用训练样本维度的类内中值(连续型变量)或者该维度类内出现频率最高的数值(类别型变量)插补缺失值；如果是回归问题，使用所有样本在该维度上的中值或者出现频率最高的数值插补缺失值。
2. 训练：使用插补后的数据集训练随机森林$F_0$
3. 计算样本的Proximities: $P(i,j)$
4. 使用上一节的公式，更新缺失值
5. 使用更新后的数据集，重新训练随机森林 $F_t$
6. 重复3-5，直到迭代次数达到预设值

以下为代码实现的流程:

具有缺失值的训练集 X = { x 1 , x 2 , . . . . . . , x N } X=\left\{ x_1,x_2,......,x_N \right\} X={x1​,x2​,......,xN​}
以及类标签 Y = { y 1 , y 2 , . . . , y N } Y=\left\{ y_1,y_2,...,y_N \right\} Y={y1​,y2​,...,yN​}
对于每个维度 $d=\mathrm{1......}P$
|  计算该维度上的非缺失样本的中值$m_d$(连续型变量)，
|  或者出现频率最高的变量值$c_d$(类别型变量)
|  对于每个样本 $i=\mathrm{1......}N$
|     如果 $x_i^d$是缺失值，那么$x_i^d\leftarrow m_d（或c_d）$
通过上述插补方法，获得经过初始填充的训练集$X_0$
训练随机森林$F_0\leftarrow fit(X_0,Y)$
迭代： $t=\mathrm{1...}T$
|  对于每个样本 $i=\mathrm{1......}N$
|  |  使用$F_{t-1}$计算样本$i$与其他样本之间的Proximities:
|  |    $P(i,j)$ $,j\in 1,2,...,N$
|  |  对于每个维度 $d=\mathrm{1......}P$
|  |  |  如果 $x_i^d$是缺失值，使用上一节公式计算：${\bar{x}}_i^d$
|  |  |  更新$X_t$: $x_i^d\leftarrow {\bar{x}}_i^d$
|  |  end of $d$
|  end of $i$
|  训练随机森林$F_t\leftarrow fit(X_t,Y)$
end of $t$

实现

上述Random-Forest-Proximity-based数据插补方法计划更新到“项目randomforest”。
算法复杂度为$O(P\ast N^2)$，其中$P$与$N$分别为维数和样本数。另外，当样本数较大时，Proximities矩阵 （$N\times N$） 无法直接计算。如对于mnist数据集（样本数60000），它的完整Proximities矩阵需要占用26G以上 ($60k\times 60k\times 8$)内存。此外，由于RF-GAP的相似度矩阵$P(i,j)\ne P(j,i)$,因此也没有必要将完整的Proximities矩阵计算出，每次计算样本$i$与其它样本的相似度即可。
以下为主干代码：

float** MissingValuesImputaion(float** data, bool* is_categorical, const float** data_orig, int* label, RandomCForests_info RFinfo, PROXIMITY_TYPE prox_type, int max_iteration, bool verbose, int random_state, int jobs)
{
	int rv = chekcDataSetLabel(label, RFinfo.datainfo.samples_num, RFinfo.datainfo.classes_num);
	if (rv < 0)
	{
		std::cout << ">>>[MissingValuesImputaion] labels are not properly set." << std::endl;
		return NULL;
	}

	const int samples_num = RFinfo.datainfo.samples_num;
	const int var_num = RFinfo.datainfo.variables_num;
	const int class_num = RFinfo.datainfo.classes_num;


	bool todel = false;
	if (is_categorical == NULL)
	{
		todel = true;
		is_categorical = new bool[RFinfo.datainfo.variables_num];
		memset(is_categorical, 0, sizeof(bool) * RFinfo.datainfo.variables_num);
	}

	unsigned char* rowmask = new unsigned char[samples_num];
	memset(rowmask, 0, sizeof(unsigned char) * samples_num);

	// generate mask of the sample row and matrix
	unsigned char** nanmask = new unsigned char* [samples_num];
	for (int n = 0; n < samples_num; n++)
	{
		nanmask[n] = new unsigned char[var_num];
		memset(nanmask[n], 0, sizeof(unsigned char) * var_num);
		for (int j = 0; j < var_num; j++)
		{
			if (std::isnan(data[n][j]))
			{
				nanmask[n][j] = 1;
				rowmask[n] = 1;
			}
		}
	}

	float** data_imp = clone_data(data, samples_num, var_num);
	bool** bzero = MedianFill(data, is_categorical, label, RFinfo, data_imp/*OUT*/);
	float** data_tmp = clone_data(data_imp, samples_num, var_num);

	int method_nrmse = 0;  // 0: RMSE; 1: NRMSE
	float *var_true = NULL; // variance of the ture(original) values along a variable(feature)
	float* mean_true = NULL;
	int* count_var = NULL; // missing values of a variable(feature)
	float* nrmse = NULL; //  the normalized root mean squared error (NRMSE)
	if (data_orig)
	{
		var_true = new float[var_num * 3];
		mean_true = var_true + var_num;
		nrmse = var_true + var_num * 2;
		memset(var_true, 0, sizeof(float) * var_num * 3);

		count_var = new int[var_num];
		memset(count_var, 0, sizeof(int) * var_num);

		PreparePrecisionCalculation(samples_num, var_num, data_orig, nanmask, count_var, mean_true, var_true);

		PrintPrecision(samples_num, var_num, data_imp, data_orig, nanmask, is_categorical, count_var, var_true, method_nrmse);
	}

	if (jobs <= 0)
		jobs = 1;

	if (jobs > 1)
	{
		std::cout << "max_threads: " << omp_get_max_threads() << std::endl;
		jobs = jobs > omp_get_max_threads() ? omp_get_max_threads() : jobs;
		omp_set_num_threads(jobs);
		std::cout << "jobs: " << jobs << std::endl;
	}

	for (int it = 0; it < max_iteration; it++)
	{
		LoquatCForest* forest = NULL;

		if (jobs > 1)
			TrainRandomForestClassifierOMP(data_tmp, label, RFinfo, forest, random_state, jobs, verbose ? 50 : 0);
		else
			TrainRandomForestClassifier(data_tmp, label, RFinfo, forest, random_state, verbose ? 50 : 0);

		LoquatCTreeNode*** leafMatrix = NULL;
		if (prox_type == PROXIMITY_TYPE::PROX_ORIGINAL)
		{
			leafMatrix = createLeafNodeMatrix(forest, data_tmp);
		}

		timeIt(1);
#pragma omp parallel for num_threads(jobs)
		for (int n = 0; n < samples_num; n++)
		{
			if (!rowmask[n])
				continue;

			float* proximity = NULL;
			int rv = 1;
			switch (prox_type)
			{
			case PROXIMITY_TYPE::PROX_ORIGINAL:
				rv = ClassificationForestOrigProximity2(forest, data_tmp[n], leafMatrix, proximity);
				break;
			case PROXIMITY_TYPE::PROX_GEO_ACC:
			default:
				rv = ClassificationForestGAPProximity(forest, data_tmp, n, proximity);
			}
			
			if (rv < 0)
			{
				delete[] proximity;
				continue;
			}

			const int lb = label[n];
			for (int v = 0; v < var_num; v++)
			{
				if (!nanmask[n][v])
					continue;

				bool use_label = !bzero[v][lb];
				if (!is_categorical[v])
				{
					float s = 0.f, estimated = 0.f;

					for (int i = 0; i < samples_num; i++)
					{
						if (!nanmask[i][v] && (!use_label || label[i] == lb))
						{
							s += proximity[i];
							estimated += data_tmp[i][v] * proximity[i];
						}
					}

					if (s > 1e-10f)
					{
						estimated = estimated / s;
						data_imp[n][v] = estimated;
					}

				}
				else
				{
					int k;
					std::unordered_map<int, float> categroy_prox;
					for (int i = 0; i < samples_num; i++)
					{
						if (!nanmask[i][v] && (!use_label || label[i] == lb))
						{
							k = int(data_tmp[i][v]);
							if (categroy_prox.end() == categroy_prox.find(k))
								categroy_prox.emplace(k, proximity[i]);
							else
								categroy_prox[k] += proximity[i];
						}
					}

					int k_max_prox = -1;
					float max_prox = -FLT_MAX;
					for (auto it = categroy_prox.begin(); it != categroy_prox.end(); ++it)
					{
						if (it->second > max_prox)
						{
							max_prox = it->second;
							k_max_prox = it->first;
						}
					}
					data_imp[n][v] = k_max_prox;
				}

			}

			delete[] proximity;

			if (n > 0 && n % 5000 == 0)
				std::cout << "samples " << n << " imputated." << std::endl;
		}

		std::cout << "proximities: " << timeIt(0) << std::endl;

		ReleaseClassificationForest(&forest);

		if (prox_type == PROXIMITY_TYPE::PROX_ORIGINAL)
		{
			for (int n = 0; n < samples_num; n++)
				delete[] leafMatrix[n];
			delete[] leafMatrix;
		}

		if (data_orig)
		{
			PrintPrecision(samples_num, var_num, data_imp, data_orig, nanmask, is_categorical, count_var, var_true, method_nrmse);
		}


		for (int n = 0; n < samples_num; n++)
			memcpy(data_tmp[n], data_imp[n], sizeof(float) * var_num);
	}


	for (int n = 0; n < samples_num; n++)
	{
		delete[] nanmask[n];
		delete[] data_tmp[n];
	}

	delete[] nanmask;
	delete[] data_tmp;
	delete[] rowmask;

	delete[] var_true;
	delete[] count_var;

	for (int v = 0; v < var_num; v++)
		delete[] bzero[v];
	delete[] bzero;

	if (todel)
		delete[] is_categorical;

	return data_imp;
}

数据插补没有实际动手写代码前感觉不那么复杂，但实际上需要考虑很多细节，比如要考虑是否某个维度上某个类别对应的数据都缺失的情况。关于归一化（标准化）预处理，如果实际工作中没有真实数据集做参考（结果对比），完全没必要对输入数据集在各自维度上进行标准化。

实验

插补的评价指标

采用NRMSE (normalized root mean squared error)指标评价算法插补数据的准确性，这个指标被其他数据插补的论文广泛使用。

$$NRMSE=\frac{1}{|\mathcal{N}|}\sum _{j\in \mathcal{N}}\sqrt{\frac{mean\left({\left({\mathbf{x}}^j-{\bar{\mathbf{x}}}^j\right)}^2\right)}{var\left({\mathbf{x}}^j\right)}}$$

$$其中,\mathcal{N}=\left\{j:\exists x_i^{j}ismissingvalues\right\}$$

${\mathbf{x}}^j$和 ${\bar{\mathbf{x}}}^j$分别表示训练集维度$j$的真实数据和插补后的数据, 仅针对缺失数据计算对应的$mean$与$var$。

“where $X^{true}$ and $X^{imp}$ are the true and imputed data matrix, respectively, and the mean and variance are computed only over the missing values.” — —
Shangzhi Hong,Henry S. Lynn. Accuracy of random-forest-basedimputation of missing data in the presenceof non-normality, non-linearity, and interaction. BMC Medical Research Methodology, 2020.

由于评价指标NRMSE要计算方差，可能存在不确定性（比如，某个维度缺失值较少方差统计不准确或者接近于0）。如果每个维度的物理属性是一致的情况下(比如字符识别数据集mnist)，或者数据的每个维度进行了归一化，也可以使用RMSE指标来评价插补的效果：
$$RMSE=\frac{1}{|\mathcal{N}|}\sum _{j\in \mathcal{N}}\sqrt{mean\left({\left({\mathbf{x}}^j-{\bar{\mathbf{x}}}^j\right)}^2\right)}$$

原始Proximity与RF-GAP对比

在若干数据集上验证上述算法的数据插补效果，并与原始随机森林数据插补方法进行对比。缺失数据的模拟方法是完全随机的（MCAR, Missing Completely At Random），即以一定概率（对于所有样本和维度是一致的）将原数据值替换为缺失值(nan)。
本文缺失值插补算法的核心是Proximities的计算，将RF-GAP 与原始随机森林Proximity方法的插补方法进行对比。迭代2次，运行3次取平均数。原始Proximities计算方法中不区分袋内(in-bag)袋外(oob)数据，直接统计样本落入到相同叶子节点的概率 $p(i,j)=\frac{1}{T}{\sum }_{t}I(lea{f}_i^t=lea{f}_j^t)$。

数据集 (样本数/特征维数/类别数或输出维度)	随机森林参数	评价方法	缺失百分比	Proximity方法
				Original	GAP
pendigits (7494/16/10)	[200,4, 40, 5]	NRMSE	10%	0.5287	0.4359
			20%	0.5637	0.4553
			30%	0.6002	0.4791
mnist(60000/780/10)	[200,27,40,5]	RMSE	10%	41.32	36.70
			20%	42.68	38.21
			30%	43.4498	39.6647
Sensorless_drive_diagnosis1 (58509/48/11)	[200,6,40,5]	RMSE	10%	0.0436	0.0371
			20%	0.0461	0.0389
			30%	0.0486	0.0417
MAGIC_Gamma_Telescope (19020/10/2)	[200,3,40,5]	NRMSE	10%	0.7707	0.6998
			20%	0.8552	0.7372
			30%	0.9424	0.7992
ionosphere（351,34,2）	[200,6,40,5]	NRMSE	10%	0.9192	0.8479
			20%	0.9278	0.8651
			30%	0.9292	0.8676
spambase(4601,57,2)	[200,7,40,5]	RMSE	10%	0.0492	0.0476
			20%	0.0504	0.0494
			30%	0.0506	0.0501
letter-recognition(20000,16,26)	[200,4, 40, 5]	NRMSE	10%	0.5959	0.4808
			20%	0.6463	0.5203
			30%	0.6981	0.5658
Concrete_Compressive_Strength (1030,8,1)	[200,3,40,5]	NRMSE	10%	0.6001	0.6388
			20%	0.6695	0.7019
			30%	0.7035	0.7491
airfoil_self_noise(1503,5,1)	[200,2,40,5]	NRMSE	10%	0.6503	0.6946
			20%	0.7136	0.7532
			30%	0.7647	0.7973

1. Sensorless_drive_diagnosis数据集使用NRMSE存在问题，因为其第40维特征的数值分布[-0.90235~3670.8]极其不均，极个别值偏离度很大，使用NRMSE在该维度上结果超过4，并不能反映插补效果。因此将每个特征(维度）均压缩到[0,1]， 使用RMSE。

以上表中，“Original”为原始随机森林Proximity，“GAP”为TPAMI2023方法，数值为原始样本与插补样本间的平均NRMSE或者RMSE，数值越小表示效果越好。此结果与论文TPAMI2023相符。

mnist可视化结果

这里使用mnist数据集，样本数为60000，特征维数780，样本是手写数字的图像，数值在0～255之间（非二值化图像）。模拟缺失数据的方法是，对于图像中非零元素按一定概率将原值替换为缺失值（nan）。这样做的目的是为了验证上述算法对于非完全随机（Missing at Random, MAR）情况下缺失值插补的效果。另一方面，这样也更利于可视化展现插补的效果。
以下随机选取了若干样本插补的效果。第一列为原始样本，第二列是具有缺失值的样本（样本上的黑点就是人为制造的缺失值），第三列是使用“中值”插补预处理的结果，第四列是使用Proximities插补后的结果。
可以观察到，虽然与原样本图像相比还有一定差距，但是Proximities插补后“缺失”情况有明显改善。如果增加一些先验信息，插补效果还能进一步提升。

附录

使用我实现的randomforest算法算法进行缺失值插补的代码如下：

// 读取原始数据集，无缺失值
char filename[500] = "E:\\randomforest\\DataSet\\Classification\\pendigits.tra";
float** data = NULL;
int* label = NULL;
Dataset_info_C datainfo;
InitalClassificationDataMatrixFormFile2(filename, data/*OUT*/, label/*OUT*/, datainfo/*OUT*/);

RandomCForests_info rfinfo;
rfinfo.datainfo = datainfo;
rfinfo.maxdepth = 40;
rfinfo.ntrees = 200;
rfinfo.mvariables = (int)sqrtf(datainfo.variables_num);
rfinfo.minsamplessplit = 5;
rfinfo.randomness = 1;

float** data_out = NULL;
int random_state = 40;

// 模拟制造缺失值
int nan_num = 0, sample_with_mv = 0;
float** data_mv = clone_data(data, datainfo.samples_num, datainfo.variables_num); 
for (int n = 0; n < rfinfo.datainfo.samples_num; n++)
{
	int sample_have_mv = false;
	for (int v = 0; v < rfinfo.datainfo.variables_num; v++)
	{
	    // 缺失比例
		if (rand_freebsd() % 100 < 10 /*&& data[n][v]>0*/ /*mnist！！！*/ )
		{
			sample_have_mv = true;
			data_mv[n][v] = std::nanf("");
			nan_num++;
		}
	}
	if (sample_have_mv)
		sample_with_mv++;
}
	
cout << "sample with missing values: " << sample_with_mv << endl;
cout << "total missing values: " << nan_num << endl;

const float** data_orig = (const float **)data;
// 标志哪些特征是categorical，以下假设全部特征都不是categorical
bool* is_cate = new bool[rfinfo.datainfo.variables_num];
memset(is_cate, 0, sizeof(bool) * rfinfo.datainfo.variables_num);

// 缺失值插补，返回插补后数据集
data_out = MissingValuesImputaion(data_mv, is_cate, data_orig, label, rfinfo, PROX_GEO_ACC, 2/*迭代次数*/, true, random_state, 8/*线程数*/);
	
for (int n = 0; n < datainfo.samples_num; n++)
{
	delete[]data_mv[n];
	delete[]data[n];
	delete[]data_out[n];
}
delete[]data;
delete[] data_mv;
delete[] data_out;
delete[] label;
delete[] is_cate;
return 0;