机器学习基础之线性模型个人文章

回归问题的一般模型如下:$$y=\sumw[i]*x[i]+b$$如下图所示,对于一维数据,线性回归就是根据给定的点$(x_i,y_i)$拟合出一条直线$$y=ax+b$$即求出系数a、b。

importwarningswarnings.filterwarnings('ignore')importpandasaspdimportnumpyasnpimportmatplotlib.pyplotasplt%matplotlibinlineimportseabornassnsimportmglearnmglearn.plots.plot_linear_regression_wave()w[0]:0.393906b:-0.031804

推广到多维数据,线性回归模型的训练过程就是寻找参数向量$\overrightarrow{w}$的过程,只是拟合的目标变为了高维平面,线性回归最常用的两种方法是最小二乘法(OLS)和梯度下降法。

最小二乘法基于这样一个目标,使得数据的实际值$y_i$与预测值$\hat{y_i}$之间的偏差最小,即损失函数最小,OLS使用均方误差(MeanSquareError,MSE)作为损失函数,优化目标为$$min\quadMSE=\frac1n\sum_{i=1}^n(y_i-\hat{y_i})^2$$对于一维数据而言$$MSE=\frac1n\sum_{i=1}^n(y_i-\theta_0-\theta_1x_i)^2$$为求最小值,需要求偏导$$\frac{\partialMSE}{\partial\theta_0}=-\frac2n\sum_{i=1}^n(y_i-\theta_0-\theta_1x_i)=0$$,$$\frac{\partialMSE}{\partial\theta_1}=-\frac2n\sum_{i=1}^n(y_i-\theta_0-\theta_1x_i)x_i=0$$联立可得$$\theta_1=\frac{\sum(y_i-\bary)(x_i-\barx)}{(x_i-\barx)^2}$$,$$\theta_0=\bary-\theta_1\barx$$

对于多元回归同理,下面是矩阵解法,损失函数定义为$$J(\mathbf\theta)=\frac{1}{2}(\mathbf{X\theta}-\mathbf{Y})^T(\mathbf{X\theta}-\mathbf{Y})$$求导$$\frac{\partial}{\partial\mathbf\theta}J(\mathbf\theta)=\mathbf{X}^T(\mathbf{X\theta}-\mathbf{Y})=0$$最终可得$$\mathbf{\theta}=(\mathbf{X^{T}X})^{-1}\mathbf{X^{T}Y}$$

最小二乘原理是通过求导的方式最小化MSE以求得参数$\theta$,下面我们介绍另一种方法梯度下降法。

梯度下降法是一个比较纯粹的计算机编程方法。

如图所示,我们知道,损失函数是系数的函数,一元线性回归有两个参数,组成了损失函数平面,我们首先随机制定一组系数,即在上图平面上随机选取一个初始点,然后同时进行以下变换$$\theta_0=\theta_0-\alpha\frac{\partialJ(\theta)}{\partial\theta_0}$$$$\theta_1=\theta_1-\alpha\frac{\partialJ(\theta)}{\partial\theta_1}$$重复该步骤直到终止。

我们来分析一下发生了什么。首先,偏导数的系数$\alpha$是正数。对于偏导数而言,当偏导大于零时候,$J(\theta)$随$\theta_i$增大而增大,同时,新的$\theta_i$小于旧的$\theta_i$,因此,$J(\theta)$减小;当偏导数小于零的时候,$J(\theta)$随着$\theta_i$增大而减小,同时,新的$\theta_i$大于旧的$\theta_i$,因此,$J(\theta)$还是减小,即每次循环,损失函数都会减小,最终到达一个局部的最小值,即极小值。

但是,我们的损失函数是凸函数,并不是有多个极小值的图形,其真实图形类似下图所示,极小值即为最小值。

算法步骤:

梯度下降法家族

python中可以使用sklearn或者statsmodel进行线性回归,sklearn偏向于机器学习,statsmodel本就是scipy.stats的补充,主要用于线性回归以及方差分析。

#使用sklearn进行线性回归fromsklearn.linear_modelimportLinearRegressionfromsklearn.model_selectionimporttrain_test_splitX,y=mglearn.datasets.make_wave(n_samples=60)#导入数据print("数据集的第一个样本点:X:{}y:{}".format(X[0],y[0]))#数据集划分,同一random_state表示对数据集进行相同的划分X_train,X_test,y_train,y_test=train_test_split(X,y,random_state=42)#标准的sklearn风格API,Model().fit(X,y)lr=LinearRegression().fit(X_train,y_train)print('系数:{}'.format(lr.coef_))print('截距:{}'.format(lr.intercept_))print('训练精度:{}'.format(lr.score(X_train,y_train)))print('测试精度:{}'.format(lr.score(X_test,y_test)))数据集的第一个样本点:X:[-0.75275929]y:-1.1807331091906834系数:[0.39390555]截距:-0.031804343026759746训练精度:0.6700890315075756测试精度:0.65933685968637sklearn中使用OLS拟合模型,score是可决系数,可以看出,测试集可决系数只有0.65左右,可以说效果并不好,这是因为原数据为一维数据,当数据维度增加时,线性模型可以变得十分强大。下面,我们再使用statsmodel来训练模型:

对于一元线性回归,当因变量y与x并不成线性关系时,无法直接使用线性回归。根据泰勒定理:令a=0可知,y可以由$x,x^2,x^3...$线性表示,因此,可以将$x^n$看作额外的变量,将一元线性回归转化为多元线性回归,以此来增加模型的准确性。

即通过取对数将原本无线性关系的变量转化为近似线性关系以应用线性回归,如对数线性回归:$$ln\mathbf{Y}=\mathbf{X\theta}$$

岭回归使用L2正则化处理回归模型,其惩罚项为L2范数,惩罚项系数为正数,对应sklearn.Ridge中参数alpha,增大alpha会导致系数趋向于0,从而降低训练集性能,是解决过拟合的一种方法,同样的,sklearn.ridge使用OLS。$$J(\mathbf\theta)=\frac{1}{2}(\mathbf{X\theta}-\mathbf{Y})^T(\mathbf{X\theta}-\mathbf{Y})+\frac{1}{2}\alpha||\theta||_2^2$$

$$\frac{\partial}{\partial\mathbf\theta}J(\mathbf\theta)=\mathbf{X}^T(\mathbf{X\theta}-\mathbf{Y})+\alpha\theta=0$$

$$\theta=(X^TX+\alphaE)^{-1}X^TY$$

#岭回归在sklearn中的实现fromsklearn.linear_modelimportRidgeX,y=mglearn.datasets.load_extended_boston()print('数据规模:{}'.format(X.shape))X_train,X_test,y_train,y_test=train_test_split(X,y,random_state=0)ridge=Ridge(alpha=1).fit(X_train,y_train)LR=LinearRegression().fit(X_train,y_train)print('线性回归精度([训练,测试]):{}'.format([LR.score(X_train,y_train),LR.score(X_test,y_test)]))print('岭回归精度([训练,测试]):{}'.format([ridge.score(X_train,y_train),ridge.score(X_test,y_test)]))数据规模:(506,104)线性回归精度([训练,测试]):[0.952051960903273,0.6074721959665708]岭回归精度([训练,测试]):[0.885796658517094,0.7527683481744755]可以看出,由于boston数据集拥有104个特征,但只有506条数据。线性回归具有十分明显的过拟合,岭回归模型训练精度低于线性回归,但是测试精度高于线性回归。

增加数据量可以缓解过拟合问题,如下图所示,当数据量增大时,线性回归的测试精度与ridge相似。

mglearn.plots.plot_ridge_n_samples()

lasso使用L1正则化,惩罚项是L1范数,但是可以使某特征系数为0,模型更容易解释,也可以呈现模型的重要特征,由于使用绝对值,存在不可导点,因此OLS、梯度下降都不可用。$$J(\mathbf\theta)=\frac{1}{2m}(\mathbf{X\theta}-\mathbf{Y})^T(\mathbf{X\theta}-\mathbf{Y})+\alpha||\theta||_1$$

#lasso实现fromsklearn.linear_modelimportLassolasso=Lasso(alpha=0.1,max_iter=1000000).fit(X_train,y_train)print('训练精度:{}'.format(lasso.score(X_train,y_train)))print('测试精度:{}'.format(lasso.score(X_test,y_test)))print('模型所用特征数:{}'.format(np.sum(lasso.coef_!=0)))训练精度:0.7709955157630054测试精度:0.6302009976110041模型所用特征数:8ElasticNetElasticNet同时使用L1、L2范数进行正则化,$$J(\mathbf\theta)=\frac{1}{2m}(\mathbf{X\theta}-\mathbf{Y})^T(\mathbf{X\theta}-\mathbf{Y})+\alpha\rho||\theta||_1+\frac{\alpha(1-\rho)}{2}||\theta||_2^2$$

前面提到广义线性回归:$$f(y)=X\theta$$使用线性模型进行二分类的想法就自然产生了:只需找一个单调可微函数将分类任务的真实标记y与线性回归模型的预测值联系起来。

如用于二分类的线性模型可以用以下公式预测:$$y=\sumw[i]*x[i]+b>0$$

常用的二分类模型有Logistic回归(LogisticRegression)和线性支持向量机(LinearSupportVectorMachine,线性SVM)。

LogisticRegression是将线性回归所得因变量y进行非线性(Sigmoid函数)变换映射到[0,1]之内,作为分类样本点分属0,1两类的概率。

Sigmoid函数$$g(z)=y=\frac{1}{1+e^{-z}}$$其图像如下:

X=np.linspace(-10,10)y=[]foriinX:y.append(1/(1+np.exp(-i)))plt.plot(X,y)plt.xlabel('z')plt.ylabel('g(z)')在x=0处函数值为0.5,x趋向于无穷时,函数值分别趋向0和1。如果$${z=X\theta}$$那么就把线性回归所得的函数值映射到了0-1之间,$g(z)$可以看作分类为1的概率,越靠近1,被分类为1的概率越大,在临界值0.5处最容易被误分类。

对y进行变形可得$$z=X\theta=ln\frac{y}{1-y}$$

我们将y看作分类为正例的概率,则z就可以看作正反概率之比(称为几率)的对数,因此,模型称为“对数几率回归”。

对于每个样本点$(x_i,y_i)$,$y_i=1,y_i=0$的概率分别为$$P(y_i=1|x_i,\theta)=h_\theta(x_i)$$$$P(y_i=0|x_i,\theta)=1-h_\theta(x_i)$$将其合并为$$P(y_i|x_i,\theta)=h_\theta(x_i)^{y_i}(1-h_\theta(x_i))^{1-y_i}$$假设每个样本点独立同分布,样本数为n,由最大似然法(MLE)构造似然函数得$$L(\theta)=\prod_{i=1}^nP(y_i|x_i,\theta)$$由于似然函数表示的是取得现有样本的概率,应当予以最大化,因此,取似然函数的对数的相反数作为损失函数$$J(\theta)=-lnL(\theta)=-\sum\limits_{i=1}^{m}(y_iln(h_{\theta}(x_i))+(1-y_i)ln(1-h_{\theta}(x_i)))$$

求偏导得$$\frac{\partial}{\partial\theta}J(\theta)=X^T(h_{\theta}(X)-Y)$$

使用梯度下降法$$\theta=\theta-\alphaX^T(h_{\theta}(X)-Y)$$

#乳腺癌数据上使用LogisticRegressionfromsklearn.datasetsimportload_breast_cancerfromsklearn.linear_modelimportLogisticRegressioncancer=load_breast_cancer()X_train,X_test,y_train,y_test=train_test_split(cancer.data,cancer.target,stratify=cancer.target,random_state=42)forC,makerinzip([0.001,1,100],['o','^','v']):logistic=LogisticRegression(C=C,penalty='l2').fit(X_train,y_train)print('训练精度(C={}):{}'.format(C,logistic.score(X_train,y_train)))print('训练精度(C={}):{}'.format(C,logistic.score(X_test,y_test)))plt.plot(logistic.coef_.T,maker,label='C={}'.format(C))plt.xticks(range(cancer.data.shape[1]),cancer.feature_names,rotation=90)plt.xlabel('CoefficientIndex')plt.ylabel('Coefficient')plt.legend()训练精度(C=0.001):0.9507042253521126训练精度(C=0.001):0.9440559440559441训练精度(C=1):0.9413145539906104训练精度(C=1):0.965034965034965训练精度(C=100):0.9507042253521126训练精度(C=100):0.958041958041958

LogisticRegression也可以使用正则化,方法同样是在损失函数后增加正则化项。sklearn中LogisticRegression默认使用L2正则化,参数penalty可修改正则化方式。上图是不同正则化参数训练所得模型系数,可以看出skleran中正则化项C越小,正则化程度越强,参数的变换范围越小。这是由于取损失函数时取了相反数,-C相当于lasso中的$\alpha$,C越小,-C越大。$$J(\mathbf\theta)=-\sum\limits_{i=1}^{m}(y_iln(h_{\theta}(x_i))+(1-y_i)ln(1-h_{\theta}(x_i)))-\frac{1}{2}C||\theta||_2^2$$

在使用逻辑回归进行预测时,通常要进行特征筛选,以便减少特征数量,增强模型可解释性并加快训练速度。

#导入数据data=pd.read_excel(r'G:\MYFiles\学习\Python数据分析\Python数据分析与挖掘实战\Python数据分析与挖掘实战\chapter5\demo\data\bankloan.xls')data.head(3)年龄教育工龄地址收入负债率信用卡负债其他负债违约041317121769.311.3593925.008608112711063117.31.3622024.000798024011514555.50.8560752.1689250X=data.iloc[:,:-1]y=data.iloc[:,-1]X_train,X_test,y_train,y_test=train_test_split(X,y,stratify=y)fromsklearn.linear_modelimportLogisticRegressionCVlr=LogisticRegressionCV(Cs=[0.0001,0.001,0.01,0.1,1,10],cv=5).fit(X_train,y_train)print("训练精度:{}\n测试精度:{}\n总精度:{}".format(lr.score(X_train,y_train),lr.score(X_test,y_test),lr.score(X,y)))print("正则化系数:%s"%lr.C_)print("变量系数:%s"%lr.coef_)训练精度:0.8285714285714286测试精度:0.7771428571428571总精度:0.8157142857142857正则化系数:[10.]变量系数:[[0.029383150.08359976-0.24002552-0.08184292-0.004432810.100485780.58023866-0.02504046]]线性SVC见SVM

线性判别分析(LinearDiscriminantAnalysis,简称LDA)是一种监督学习的降维方法。

LDA的思想非常朴素。给定训练集,将样例投影到一条直线上,使得同类样例的投影点尽可能接近、异类样例的投影点尽可能远离;在对新样本进行分类时,将其投影到同样的这条直线上,再根据投影点的位置来确定新样本的类别。

一般来说,如果我们的数据是有类别标签的,那么优先选择LDA去尝试降维;当然也可以使用PCA做很小幅度的降维去消去噪声,然后再使用LDA降维。如果没有类别标签,那么肯定PCA是最先考虑的一个选择了。以下是一个使用LDA进行降维的例子:

frommpl_toolkits.mplot3dimportAxes3Dfromsklearn.datasets.samples_generatorimportmake_classificationX,y=make_classification(n_samples=1000,n_features=3,n_redundant=0,n_classes=3,n_informative=2,n_clusters_per_class=1,class_sep=0.5,random_state=10)fig=plt.figure()ax=Axes3D(fig,rect=[0,0,1,1],elev=30,azim=20)ax.scatter(X[:,0],X[:,1],X[:,2],marker='o',c=y)plt.show()

fromsklearn.discriminant_analysisimportLinearDiscriminantAnalysislda=LinearDiscriminantAnalysis(n_components=2).fit(X,y)X_trans=lda.transform(X)plt.scatter(X_trans[:,0],X_trans[:,1],marker='o',c=y)plt.show()

许多线性模型不适用于多分类问题,通常,我们基于一些基本策略,利用二分类学习器来解决多分类问题。多分类学习的基本思路是“拆解法”,即将多分类任务拆为若干个二分类任务求解。具体来说,先对问题进行拆分,然后为拆出的每个二分类任务训练一个分类器;在测试时,对这些分类器的预测结果进行集成以获得最终的多分类结果。有两点关键:

一对一(OvO):对样本进行N分类,将类别两两组合,共$C_N^2=\frac{N(N-1)}2$个二分类器,使用每个分类器对二分类器预测,对预测结果使用多数表决的方法选出最终结果。

一对其余(OVR):对样本进行N分类,对每个类别构造一个二分类器,共N个二分类器,二分类器的结果是“是”或者“不是”该类别。使用N个分类器对每个样本进行分类,被分类为正例的则为该样例的标签,若被预测为多个标签则选择置信度较高的标签作为预测结果。

多对多(MvM):每次将若干个类作为正类,若干个其他类作为反类。下面是一个多对多的常用方法:

纠错输出码(ErrorCorrectingOutputCodes,简称ECOC)工作分为两步:

如图所示,是两种编码方式。对于左侧的二元码,构造五个二分类器,也就是对数据集进行五次划分,如$f_1$将$C_2$类数据划分为正例,其他划分为反例。对测试样例使用五个分类器分别分类,得到五个分类标签,作为分类向量,计算分类向量与各个类别分类向量的距离,即可衡量测试样例与五个类别的差异程度,选择差异最小的类别作为最终分类类别。右侧分类器只是划分时分为三个类别。左边分类结果是$C_3$,右侧是$C_2$。

ECOC编码对分类器的错误有一定的容忍和修正能力。例如图(a)中对测试示例的正确预测编码是(-1,+1,+1,-1,+1),在预测时某个分类器出错了,例如f2出错从而导致了错误编码(-1,-1,+1,-1,+1),但基于这个编码仍能产生正确的最终分类结果$C_3$.一般来说对同一个学习任务,ECOC编码越长,纠错能力越强。但是,编码越长,意味着所需训练的分类器越多,计算、存储开销都会增大;另一方面,对有限类别数,可能的组合数目是有限的,码长超过一定范围后就失去了意义。

sklearn中,线性SVC使用的就是一对其余方法进行多分类,数据如下所示:

fromsklearn.datasetsimportmake_blobsX,y=make_blobs(random_state=42)mglearn.discrete_scatter(X[:,0],X[:,1],y)plt.xlabel('Feature0')plt.ylabel('Featrue1')plt.legend(['Class0','Class1','Class2'])plt.show()

对于以上数据使用线性SVC进行分类:

fromsklearn.svmimportLinearSVCLSVC=LinearSVC().fit(X,y)print("系数矩阵:{}\n截距:{}\n评分:{}".format(LSVC.coef_,LSVC.intercept_,LSVC.score(X,y)))系数矩阵:[[-0.174924820.23141171][0.47621845-0.06937277][-0.18913971-0.20400358]]截距:[-1.077455270.13140594-0.08604909]评分:1.0可以看出,LinearSVC输出了三条直线,下面将其可视化,三条直线将其分为7个区域,交叉区域平均分配。

mglearn.plots.plot_2d_classification(LSVC,X,fill=True,alpha=.7)mglearn.discrete_scatter(X[:,0],X[:,1],y)line=np.linspace(-10,10)forcoef,intercept,colorinzip(LSVC.coef_,LSVC.intercept_,['b','r','g']):plt.plot(line,-(line*coef[0]+intercept)/coef[1])plt.xlabel('Feature0')plt.ylabel('Featrue1')plt.legend(['Class0','Class1','Class2','Lineclass0','Lineclass1','Lineclass2'],loc=(1.01,0))

类别不平衡(classimbalance)就是指分类任务中不同类别的训练样例数目差别很大的情况。

比如,一个样本,正例2个,反例998个,只要分类器预测全部为反例,就能得到99.8%的精度,但这个精度是没有意义的。

再缩放:判定正例的临界条件改为$$\frac{y}{1-y}\times\frac{m^-}{m^+}>1$$

在使用分类器进行分类时,分类器给出一个预测值y,将y与阈值进行比较,如通常在$y>0.5$的时候预测正例,否则反例,y反映的是正例的可能性,设置阈值为0.5,则是认为正例与反例的可能性相等。

然而,当训练集中正、反例的数目不同时,令$m^+$表示正例数目,$m^-$表示反例数目,则观测几率是$\frac{m^+}{m^-}$;,由于我们通常假设训练集是真实样本总体的无偏采样,因此观测几率就代表了真实几率。

在实际中,针对类别不平衡问题通常有以下几种处理方法:

欠采样法若随机丢弃反例,可能丢失一些重要信息。欠采样法的代表性算法EasyEnsemble则是利用集成学习机制,将反例划分为若干个集合供不同学习器使用,这样对每个学习器来看都进行了欠采样,但在全局来看却不会丢失重要信息。

过采样法不能简单地对初始正例样本进行重复采样,否则会招致严重的过拟合,过采样法的代表性算法SMOTE是通过对训练集里的正例进行插值来产生额外的正例。

THE END
1.数学优化算法最小二乘法? 原理:共轭梯度法结合了最速下降法和牛顿法的优点,它利用一阶导数信息并沿着共轭方向进行搜索,以加快收敛速度。 ? 优点:所需存储量小,具有步收敛性,稳定性高。 四、应用 最小二乘法优化算法广泛应用于各种领域,如回归分析、曲线拟合、参数估计、图像处理等。特别是在机器学习和深度学习中,最小二乘法作为https://blog.csdn.net/xioayanran123/article/details/144229105
2.普通最小二乘法(OLS)(空间统计)—ArcMap文档执行全局“普通最小二乘法 (OLS)”线性回归可生成预测,也可为一个因变量针对它与一组解释变量关系建模。 可从结果窗口获取此工具的结果(包括可选报表文件)。如果禁用了后台处理,结果也将被写入进度对话框。 注: 此工具的功能包含在ArcGIS Pro 2.3新增的广义线性回归工具中。广义线性回归工具支持其他模型。 https://desktop.arcgis.com/zh-cn/arcmap/latest/tools/spatial-statistics-toolbox/ordinary-least-squares.htm
3.常用算法分析——最小二乘法普通最小二乘法(OLS) OLS实现 广义最小二乘法(GLS)简介 1、引言 最小二乘法应该是我们最早接触的一种数值估计算法。它的特殊形式——一元线性回归,被广泛地应用于多种数值统计分析场合。例如,在验证欧姆定律(U=IR)时,通常的实验方法是分别测量出多个不同电压Ui下,通过电阻的电流值Ii,然后将这些(Ui,Ii)观测https://www.jianshu.com/p/3c058de103bf
4.SPSS普通最小二乘法大家好,我现在正在写毕业论文,用到普通最小二乘法,但是我不知道要怎么用,求救大家。 公式: ⅠDAⅠ = β0+β1*M+β2*M²+β3*M³+β4*REWARD+β5*SIZE+β6*DEBT+ξ 因变量:ⅠDAⅠ为操控性应计利润绝对值;自变量:M为高管持股比例; https://bbs.pinggu.org/jg/huiji_huijiku_965982_1.html
5.计量经济学中的普通最小二乘法(OLS)的4个基本假设条件是什么什么是广义最小二乘法GLS?与普通最小二乘法OLS有什么区别? 计量经济学中的OLS是什么意思 如何用SPSS17.0进行普通最小二乘法分析数据?是有分析--回归里么?里面只有二阶最小二乘法.没有OLS啊. 特别推荐 热点考点 2022年高考真题试卷汇总 2022年高中期中试卷汇总 2022年高中期末试卷汇总 2022年高中月考试卷汇总https://www.zybang.com/question/c923cf08f1a87b927b91ebffd14f2aa3.html
6.§广义最小二乘法.doc矩阵的估计为 二、广义最小二乘法的示例 湖北省病虫灾成灾面积与受灾面积对应关系的分析 病虫灾成灾面积与受灾面积的对应关系的研究对于指导抗灾、救灾有着重大的意义。从统计分析的角度出发,利用逐年的统计资料将病虫灾成灾面积数据看成时间序列,病虫灾受灾面积数据看成时间序列,应用普通最小二乘法可以建立线性模型给https://max.book118.com/html/2017/0805/126114905.shtm
7.2012年1月计量经济学自考试题A.普通最小二乘法 B.广义差分法 C.间接最小二乘法 D.阿尔蒙多项式法 18.当替代弹性σ→1,替代参数ρ→0时,CES生产函数趋于( ) A.线性生产函数 B.C—D生产函数 C.投入产出函数 D.其它 19.进行宏观经济模型的总体设计时,首先需确定( ) A.模型的结构 B.函数形式 https://www.hbzkw.com/exam/20120217170032.html
8.科学网—线性回归最小二乘法和梯度下降法mwy线性回归-最小二乘法和梯度下降法mwy 什么是一元线性模型呢?监督学习中,如果预测的变量是离散的,我们称其为分类(如决策树,支持向量机等),如果预测的变量是连续的,我们称其为回归。回归分析中,如果只包括一个自变量和一个因变量,且二者的关系可用一条直线近似表示,这种回归分析称为一元线性回归分析。如果回归分析https://blog.sciencenet.cn/blog-3413658-1177223.html