断点分析法_局部线性回归_matlab

本文介绍断点分析法在数据作假方面的应用：

题目：

近年来空气质量问题始终是政府、环境保护部门和全国人民关注的热点问题。为了激励城市政府重视空气污染治理，地方官员的政绩考核中往往包括诸如“蓝天数”这样的指标，即全年空气污染指数低于100点的天数。搜集相关空气质量和气候数据，利用数学模型分析其数据的真实性。

思路：

使用不连续回归分析的方法。经验判断，在AQI为100左右最可能发生数据的操作，因为如果数据修改幅度过大，可能会引起公众和其他政府官员的怀疑，因为公众可以获得每日API指数；但是API指数在和左右是难以辨别的，在这里修改数据也不易被发现。因此，在存在数据被修改的情况下可能会存在断点。依据判断空气质量数据的概率密度分布曲线是否在达标指数附近存在不连续，我们可以初步判断数据是否被人为修改。

说明：

断点分析法一般需要在潜在断点左右两侧分别进行局部回归（一般的回归在边界处的估计存在较大偏差，采用局部线性回归可以很好的解决边界问题）；而局部回归多采用局部多项式回归，局部线性回归是局部多项式回归的特例，因为其算法较为简单且性能优越，我们这里采用局部线性回归。

断点分析法步骤：

1.在100左右两侧分别进行频率分直方图的建立（程序中有体现）

2. 取直方图的中点作为自变量，对应的概率密度作为应变量（到这里为止，我们得到拟合所要数据）

3.对左右两侧的数据分别进行局部线性回归

4.以下图为例：

定性分析：在150左侧明显存在不连续，说明很有可能存在数据作假。

定量分析：可以使用150的左右极限的差值的大小定量分析数据作假的可能性

说明1：差值越大，表示越有可能发生数据操作，不代表作假的程度越大

说明2：局部线性回归的带宽与直方图组距的选择很重要（后面有提到）

局部多项式回归：

多项式回归模型如下：

其中 m(x) 表示回归函数，y_i 表示第个 i 采样点 x_i 处的采样值，ε_i表示独立分布零均值噪声

在 m(x) 的形式未确定的情况下，假设是 N 阶局部平滑的，为了估计函数在给定数据下的任意点处的值，我们可以将函数在这一点局部展开。假设 X_i 是 X 附近的采样点，则有N阶泰勒展开式：

首先定义权函数：

其中 为核函数，它以估计点为中心，用来控制各个采样点的权重：距离
X 越近的点，权重越大，h 为带宽（平滑参数），用于控制核的尺度。核函数K( ) 形式不固定，需满足关于 y 轴对称并在零点处取最大值，在这里我们使用高斯核

选择 m(x) 来使得下面的局部加权平方和 Q 最小

估计依赖于目标值
X ，最终有

当N=1，为局部线性回归，因为其算法较为简单且性能优越，所以我们采用局部线性回归分别对临界值左右两侧的数据进行拟合。

局部线性参数的求解：

分别求 的偏导数，并让它们等于零，这里以线性回归进行局部的拟合，所以只需求的偏导数。

引入矩阵

整理关于的线性方程组，使用矩阵表示如下：

若存在，则有：

所以有 y 的估计值：

直方图组距与核函数带宽的选取：

进行局部线性回归进行估计分两步：第一步需确保箱体大到包含足够多的样本使其样本点在临界值两边都比较平滑，但又要小到一定程度使得样本点在临界值处的跳跃能都明显的显现出来，这就需要选择合适的b；第二步以直方图箱体的中点作为观测变量，以对应的概率密度作为结果变量，采用局部线性估计y。但是潜在的不连续点不应包括在箱体中，分别在潜在不连续点左右两侧进行直方图的绘制，获得对应的样本点。

因所获得数据较少，b选择不易过大，但是不宜过小。我们使用，的标准差。带h宽控制着核函数权重变化的速率，即用来控制核的尺度，我们使用McCrary建议的带宽。

说明：组距与带宽的选择有专门的研究，这里不进行深入研究，大家可以根据实际情况进行合适的选择。下图说明：当带宽选择不当会出现过拟合与欠拟合。

以下是matlab代码：

%% 局部线性拟合
PM10 = NO2;
data = PM10(2:482);

cut_off = 120;
x_lim_all = [min(PM10),max(PM10)];

% 获取 150 左侧的样本点
data_low = data(data<cut_off);
bin_width_low = 2*std(data_low)/sqrt(length(data_low));
bin_num_low = (cut_off-min(data_low))/bin_width_low;
[y_low,x_low] = hist(data_low,bin_num_low);
y_low=y_low/length(PM10)/mean(diff(x_low));
% x_low(y_low==0) = [];
% y_low(y_low==0) = [];
% 获取 150 右侧的样本点
data_up = data(data>cut_off);
bin_width_up = 2*std(data_up)/sqrt(length(data_up));
bin_num_up = (max(data_up)-cut_off)/bin_width_up;
[y_up,x_up] = hist(data_up,bin_num_up);
y_up=y_up/length(PM10)/mean(diff(x_up));
% x_up(y_up==0) = [];
% y_up(y_up==0) = [];

y_lim_all = [min([min(y_low) min(y_up)]) max([max(y_low) max(y_up)])];

% 拟合
% 绘制样本点
colors = ['m' 'c' 'k' 'r' 'g'];
figure;
hold on;
plot([x_low x_up],[y_low y_up],'.b');
% 局部线性拟合
% 150 左侧拟合
h_bandwidth = [5*bin_width_low 15*bin_width_low 20*bin_width_low]; % bandwidth
[x_low_fit,y_low_fit] = local_linear_fit(x_low',y_low',h_bandwidth);
for k = 1:length(h_bandwidth)
    plot(x_low_fit(:,k),y_low_fit(:,k),colors(k))
end
% 150 右侧拟合
h_bandwidth = [10*bin_width_up 15*bin_width_up 20*bin_width_up]; % bandwidth
[x_up_fit,y_up_fit] = local_linear_fit(x_up',y_up',h_bandwidth);
for k = 1:length(h_bandwidth)
    plot(x_up_fit(:,k),y_up_fit(:,k),colors(k))
end

legend('trainingdata','a=10,','a=15','a=20')
axis([x_lim_all*1.1 y_lim_all*1.1])
%% 鲁棒性带宽过大过小的比较
figure;
hold on;
plot([x_low x_up],[y_low y_up],'.b');
% 局部线性拟合
h_bandwidth = [1*bin_width_low 15*bin_width_low 30*bin_width_low]; % bandwidth
local_linear_fit(x_low',y_low',h_bandwidth)

h_bandwidth = [1*bin_width_up 15*bin_width_up 30*bin_width_up]; % bandwidth
local_linear_fit(x_up',y_up',h_bandwidth)

legend('trainingdata','a=1,','a=15','a=30')
axis([x_lim_all*1.1 y_lim_all*1.1])

以下为使用matlab编写的局部线性回归函数：（matlab中没有该函数，该函数中附有编者邮箱，具体问题可以交流）

function [fit_x,fit_y] = local_linear_fit(x,y,h_bandwidth)
% local_linear_fit:local linear regression
%   [fit_x,fit_y] = local_linear_fit(x,y,h_bandwidth),x: independent
%   variable,y:dependent variable,h_bandwidth:bandwidth
%
% ncf,July,2016
% Email:[email protected]
% log:
% 2016-7-12:Complete

% column matrix
test_sample_x = x;
test_sample_y = y;

test_sample_X = [ones(length(test_sample_x),1) test_sample_x];

num_x = length(test_sample_x);
% Weight_speed
Weight_speed = 1;

%% linear fit
% line_fit_beta = (test_sample_X'*test_sample_X)\(test_sample_X'*test_sample_y);
% line_fit_y = line_fit_beta(1) + line_fit_beta(2)*test_sample_x;
% plot(test_sample_x,line_fit_y,'b')
%% local linear fit
% fit_x and fit_y
fit_x = min(test_sample_x):0.3:max(test_sample_x);
fit_y = zeros(length(test_sample_x),length(h_bandwidth));

% colors = ['m' 'c' 'k' 'r' 'g']; 

for k_bandwidth = 1:length(h_bandwidth)
    h = h_bandwidth(k_bandwidth);
    for k_fit_y = 1:length(fit_x)
        w = zeros(num_x,num_x);
        K_h_all = zeros(num_x,1);
        % compute K_h
        for k_w = 1:num_x
            K_h_all(k_w) = gaussian_kernel((fit_x(k_fit_y)-test_sample_x(k_w))/h,Weight_speed)/h;
        end
        sum_K_h_all = sum(K_h_all);
        % compute w
        for k_w = 1:num_x
           w(k_w,k_w) = K_h_all(k_w)./sum_K_h_all;
        end
        local_beta = (test_sample_X'*w*test_sample_X)\(test_sample_X'*w*test_sample_y);
        fit_y(k_fit_y,k_bandwidth) = local_beta(1)+local_beta(2)*fit_x(k_fit_y);
    end
%     plot(fit_x,fit_y,colors(k_bandwidth))
end
% legend('trainingdata','linear','r=.1','r=.3','r=.8,','r=2','r=10');
fit_x = [fit_x',fit_x',fit_x'];

以下为高斯核函数：

function k = gaussian_kernel(x,Weight_speed)
% Weight_speed 越大，权值变化越慢

k = (1/(sqrt(2*pi)))*exp(-(x^2)/(2*Weight_speed^2));

参考文献：

[1]McCrary, Justin. 2008. “Manipulation of the runningvariable in the regression discontinuity design:A density test.” Journal of Econometrics 142 (2):698–714.

[2]张煜东，颜俊，王水花，吴乐南.非参数估计方法[J].武汉工程大学学报，2010，99-106.

[3]Ghanem,D., & Zhang, J. (2014). ‘Effortlessperfection:’ Do Chinese cities manipulate air pollutiondata?. Journal of Environmental Economics and Management, 68(2), 203-225.

[4]欧祖军，李洪毅.局部多项式估计的带宽选择[J].长春大学学报，2007,17-19.