文档章节

C++最小二乘法拟合-(线性拟合和多项式拟合)

尘中远
 尘中远
发布于 2016/05/12 23:54
字数 2764
阅读 46
收藏 0
点赞 2
评论 0


在进行曲线拟合时用的最多的是最小二乘法,其中以一元函数(线性)和多元函数(多项式)居多,下面这个类专门用于进行多项式拟合,可以根据用户输入的阶次进行多项式拟合,算法来自于网上,和GSL的拟合算法对比过,没有问题。此类在拟合完后还能计算拟合之后的误差:SSE(剩余平方和),SSR(回归平方和),RMSE(均方根误差),R-square(确定系数)。


1.fit类的实现

先看看fit类的代码:(只有一个头文件方便使用)

#ifndef CZY_MATH_FIT
#define CZY_MATH_FIT
#include <vector>
/*
尘中远,于2014.03.20
主页:http://blog.csdn.net/czyt1988/article/details/21743595
参考:http://blog.csdn.net/maozefa/article/details/1725535
*/
namespace czy{
	///
	/// \brief 曲线拟合类
	///
	class Fit{
		std::vector<double> factor; ///<拟合后的方程系数
		double ssr;                 ///<回归平方和
		double sse;                 ///<(剩余平方和)
		double rmse;                ///<RMSE均方根误差
		std::vector<double> fitedYs;///<存放拟合后的y值,在拟合时可设置为不保存节省内存
	public:
		Fit():ssr(0),sse(0),rmse(0){factor.resize(2,0);}
		~Fit(){}
		///
		/// \brief 直线拟合-一元回归,拟合的结果可以使用getFactor获取,或者使用getSlope获取斜率,getIntercept获取截距
		/// \param x 观察值的x
		/// \param y 观察值的y
		/// \param isSaveFitYs 拟合后的数据是否保存,默认否
		///
		template<typename T>
		bool linearFit(const std::vector<typename T>& x, const std::vector<typename T>& y,bool isSaveFitYs=false)
		{
			return linearFit(&x[0],&y[0],getSeriesLength(x,y),isSaveFitYs);
		}
		template<typename T>
		bool linearFit(const T* x, const T* y,size_t length,bool isSaveFitYs=false)
		{
			factor.resize(2,0);
			typename T t1=0, t2=0, t3=0, t4=0;
			for(int i=0; i<length; ++i)
			{
				t1 += x[i]*x[i];
				t2 += x[i];
				t3 += x[i]*y[i];
				t4 += y[i];
			}
			factor[1] = (t3*length - t2*t4) / (t1*length - t2*t2);
			factor[0] = (t1*t4 - t2*t3) / (t1*length - t2*t2);
			//////////////////////////////////////////////////////////////////////////
			//计算误差
			calcError(x,y,length,this->ssr,this->sse,this->rmse,isSaveFitYs);
			return true;
		}
		///
		/// \brief 多项式拟合,拟合y=a0+a1*x+a2*x^2+……+apoly_n*x^poly_n
		/// \param x 观察值的x
		/// \param y 观察值的y
		/// \param poly_n 期望拟合的阶数,若poly_n=2,则y=a0+a1*x+a2*x^2
		/// \param isSaveFitYs 拟合后的数据是否保存,默认是
		/// 
		template<typename T>
		void polyfit(const std::vector<typename T>& x
			,const std::vector<typename T>& y
			,int poly_n
			,bool isSaveFitYs=true)
		{
			polyfit(&x[0],&y[0],getSeriesLength(x,y),poly_n,isSaveFitYs);
		}
		template<typename T>
		void polyfit(const T* x,const T* y,size_t length,int poly_n,bool isSaveFitYs=true)
		{
			factor.resize(poly_n+1,0);
			int i,j;
			//double *tempx,*tempy,*sumxx,*sumxy,*ata;
			std::vector<double> tempx(length,1.0);

			std::vector<double> tempy(y,y+length);

			std::vector<double> sumxx(poly_n*2+1);
			std::vector<double> ata((poly_n+1)*(poly_n+1));
			std::vector<double> sumxy(poly_n+1);
			for (i=0;i<2*poly_n+1;i++){
				for (sumxx[i]=0,j=0;j<length;j++)
				{
					sumxx[i]+=tempx[j];
					tempx[j]*=x[j];
				}
			}
			for (i=0;i<poly_n+1;i++){
				for (sumxy[i]=0,j=0;j<length;j++)
				{
					sumxy[i]+=tempy[j];
					tempy[j]*=x[j];
				}
			}
			for (i=0;i<poly_n+1;i++)
				for (j=0;j<poly_n+1;j++)
					ata[i*(poly_n+1)+j]=sumxx[i+j];
			gauss_solve(poly_n+1,ata,factor,sumxy);
			//计算拟合后的数据并计算误差
			fitedYs.reserve(length);
			calcError(&x[0],&y[0],length,this->ssr,this->sse,this->rmse,isSaveFitYs);

		}
		/// 
		/// \brief 获取系数
		/// \param 存放系数的数组
		///
		void getFactor(std::vector<double>& factor){factor = this->factor;}
		/// 
		/// \brief 获取拟合方程对应的y值,前提是拟合时设置isSaveFitYs为true
		///
		void getFitedYs(std::vector<double>& fitedYs){fitedYs = this->fitedYs;}

		/// 
		/// \brief 根据x获取拟合方程的y值
		/// \return 返回x对应的y值
		///
		template<typename T>
		double getY(const T x) const
		{
			double ans(0);
			for (size_t i=0;i<factor.size();++i)
			{
				ans += factor[i]*pow((double)x,(int)i);
			}
			return ans;
		}
		/// 
		/// \brief 获取斜率
		/// \return 斜率值
		///
		double getSlope(){return factor[1];}
		/// 
		/// \brief 获取截距
		/// \return 截距值
		///
		double getIntercept(){return factor[0];}
		/// 
		/// \brief 剩余平方和
		/// \return 剩余平方和
		///
		double getSSE(){return sse;}
		/// 
		/// \brief 回归平方和
		/// \return 回归平方和
		///
		double getSSR(){return ssr;}
		/// 
		/// \brief 均方根误差
		/// \return 均方根误差
		///
		double getRMSE(){return rmse;}
		/// 
		/// \brief 确定系数,系数是0~1之间的数,是数理上判定拟合优度的一个量
		/// \return 确定系数
		///
		double getR_square(){return 1-(sse/(ssr+sse));}
		/// 
		/// \brief 获取两个vector的安全size
		/// \return 最小的一个长度
		///
		template<typename T>
		size_t getSeriesLength(const std::vector<typename T>& x
			,const std::vector<typename T>& y)
		{
			return (x.size() > y.size() ? y.size() : x.size());
		}
		/// 
		/// \brief 计算均值
		/// \return 均值
		///
		template <typename T>
		static T Mean(const std::vector<T>& v)
		{
			return Mean(&v[0],v.size());
		}
		template <typename T>
		static T Mean(const T* v,size_t length)
		{
			T total(0);
			for (size_t i=0;i<length;++i)
			{
				total += v[i];
			}
			return (total / length);
		}
		/// 
		/// \brief 获取拟合方程系数的个数
		/// \return 拟合方程系数的个数
		///
		size_t getFactorSize(){return factor.size();}
		/// 
		/// \brief 根据阶次获取拟合方程的系数,
		/// 如getFactor(2),就是获取y=a0+a1*x+a2*x^2+……+apoly_n*x^poly_n中a2的值
		/// \return 拟合方程的系数
		///
		double getFactor(size_t i){return factor.at(i);}
	private:
		template<typename T>
		void calcError(const T* x
			,const T* y
			,size_t length
			,double& r_ssr
			,double& r_sse
			,double& r_rmse
			,bool isSaveFitYs=true
			)
		{
			T mean_y = Mean<T>(y,length);
			T yi(0);
			fitedYs.reserve(length);
			for (int i=0; i<length; ++i)
			{
				yi = getY(x[i]);
				r_ssr += ((yi-mean_y)*(yi-mean_y));//计算回归平方和
				r_sse += ((yi-y[i])*(yi-y[i]));//残差平方和
				if (isSaveFitYs)
				{
					fitedYs.push_back(double(yi));
				}
			}
			r_rmse = sqrt(r_sse/(double(length)));
		}
		template<typename T>
		void gauss_solve(int n
			,std::vector<typename T>& A
			,std::vector<typename T>& x
			,std::vector<typename T>& b)
		{
			gauss_solve(n,&A[0],&x[0],&b[0]);	
		}
		template<typename T>
		void gauss_solve(int n
			,T* A
			,T* x
			,T* b)
		{
			int i,j,k,r;
			double max;
			for (k=0;k<n-1;k++)
			{
				max=fabs(A[k*n+k]); /*find maxmum*/
				r=k;
				for (i=k+1;i<n-1;i++){
					if (max<fabs(A[i*n+i]))
					{
						max=fabs(A[i*n+i]);
						r=i;
					}
				}
				if (r!=k){
					for (i=0;i<n;i++)         /*change array:A[k]&A[r] */
					{
						max=A[k*n+i];
						A[k*n+i]=A[r*n+i];
						A[r*n+i]=max;
					}
				}
				max=b[k];                    /*change array:b[k]&b[r]     */
				b[k]=b[r];
				b[r]=max;
				for (i=k+1;i<n;i++)
				{
					for (j=k+1;j<n;j++)
						A[i*n+j]-=A[i*n+k]*A[k*n+j]/A[k*n+k];
					b[i]-=A[i*n+k]*b[k]/A[k*n+k];
				}
			} 

			for (i=n-1;i>=0;x[i]/=A[i*n+i],i--)
				for (j=i+1,x[i]=b[i];j<n;j++)
					x[i]-=A[i*n+j]*x[j];
		}
	};
}


#endif


为了防止重命名,把其放置于czy的命名空间中,此类主要两个函数:

1.求解线性拟合:

///
/// \brief 直线拟合-一元回归,拟合的结果可以使用getFactor获取,或者使用getSlope获取斜率,getIntercept获取截距
/// \param x 观察值的x
/// \param y 观察值的y
/// \param length x,y数组的长度
/// \param isSaveFitYs 拟合后的数据是否保存,默认否
///
template<typename T>
bool linearFit(const std::vector<typename T>& x, const std::vector<typename T>& y,bool isSaveFitYs=false);
template<typename T>
bool linearFit(const T* x, const T* y,size_t length,bool isSaveFitYs=false);


2.多项式拟合:

///
/// \brief 多项式拟合,拟合y=a0+a1*x+a2*x^2+……+apoly_n*x^poly_n
/// \param x 观察值的x
/// \param y 观察值的y
/// \param length x,y数组的长度
/// \param poly_n 期望拟合的阶数,若poly_n=2,则y=a0+a1*x+a2*x^2
/// \param isSaveFitYs 拟合后的数据是否保存,默认是
/// 
template<typename T>
void polyfit(const std::vector<typename T>& x,const std::vector<typename T>& y,int poly_n,bool isSaveFitYs=true);
template<typename T>
void polyfit(const T* x,const T* y,size_t length,int poly_n,bool isSaveFitYs=true);


这两个函数都用模板函数形式写,主要是为了能使用于float和double两种数据类型


2.fit类的MFC示范程序

下面看看如何使用这个类,以MFC示范,使用了开源的绘图控件Hight-Speed Charting,使用方法见 http://blog.csdn.net/czyt1988/article/details/8740500

新建对话框文件,

对话框资源文件如图所示:


加入下面的这些变量:

std::vector<double> m_x,m_y,m_yploy;
	const size_t m_size;
	CChartLineSerie *m_pLineSerie1;
	CChartLineSerie *m_pLineSerie2;

由于m_size是常量,因此需要在构造函数进行初始化,如:

ClineFitDlg::ClineFitDlg(CWnd* pParent /*=NULL*/)
	: CDialogEx(ClineFitDlg::IDD, pParent)
	,m_size(512)
	,m_pLineSerie1(NULL)


初始化两条曲线:

CChartAxis *pAxis = NULL; 
pAxis = m_chartCtrl.CreateStandardAxis(CChartCtrl::BottomAxis);
pAxis->SetAutomatic(true);
pAxis = m_chartCtrl.CreateStandardAxis(CChartCtrl::LeftAxis);
pAxis->SetAutomatic(true);
m_x.resize(m_size);
m_y.resize(m_size);
m_yploy.resize(m_size);
for(size_t i =0;i<m_size;++i)
{
	m_x[i] = i;
	m_y[i] = i+randf(-25,28);
	m_yploy[i] = 0.005*pow(double(i),2)+0.0012*i+4+randf(-25,25);
}
m_chartCtrl.RemoveAllSeries();//先清空
m_pLineSerie1 = m_chartCtrl.CreateLineSerie();	
m_pLineSerie1->SetSeriesOrdering(poNoOrdering);//设置为无序
m_pLineSerie1->AddPoints(&m_x[0], &m_y[0], m_size);
m_pLineSerie1->SetName(_T("线性数据"));
m_pLineSerie2 = m_chartCtrl.CreateLineSerie();	
m_pLineSerie2->SetSeriesOrdering(poNoOrdering);//设置为无序
m_pLineSerie2->AddPoints(&m_x[0], &m_yploy[0], m_size);
m_pLineSerie2->SetName(_T("多项式数据"));

rangf是随机数生成函数,实现如下:

double ClineFitDlg::randf(double min,double max)
{
	int minInteger = (int)(min*10000);
	int maxInteger = (int)(max*10000);
	int randInteger = rand()*rand();
	int diffInteger = maxInteger - minInteger;
	int resultInteger = randInteger % diffInteger + minInteger;
	return resultInteger/10000.0;
}

运行程序,如图所示


线性拟合的使用如下:

void ClineFitDlg::OnBnClickedButton1()
{
	CString str,strTemp;
	czy::Fit fit;
	fit.linearFit(m_x,m_y);
	str.Format(_T("方程:y=%gx+%g\r\n误差:ssr:%g,sse=%g,rmse:%g,确定系数:%g"),fit.getSlope(),fit.getIntercept()
		,fit.getSSR(),fit.getSSE(),fit.getRMSE(),fit.getR_square());
	GetDlgItemText(IDC_EDIT,strTemp);
	SetDlgItemText(IDC_EDIT,strTemp+_T("\r\n------------------------\r\n")+str);
	//在图上绘制拟合的曲线
	CChartLineSerie* pfitLineSerie1 = m_chartCtrl.CreateLineSerie();	
	std::vector<double> x(2,0),y(2,0);
	x[0] = 0;x[1] = m_size-1;
	y[0] = fit.getY(x[0]);y[1] = fit.getY(x[1]);
	pfitLineSerie1->SetSeriesOrdering(poNoOrdering);//设置为无序
	pfitLineSerie1->AddPoints(&x[0], &y[0], 2);
	pfitLineSerie1->SetName(_T("拟合方程"));//SetName的作用将在后面讲到
	pfitLineSerie1->SetWidth(2);
}

需要如下步骤:

  • 声明Fit类,用于头文件在czy命名空间中,因此需要显示声明命名空间名称czy::Fit fit;
  • 把观察数据输入进行拟合,由于是线性拟合,可以使用LinearFit函数,此函数把观察量的x值和y值传入即可进行拟合
  • 拟合完后,拟合的相关结果保存在czy::Fit里面,可以通过相关方法调用,方法在头文件中都有详细说明

运行结果如图所示:



多项式拟合的使用如下:

void ClineFitDlg::OnBnClickedButton2()
{
	CString str;
	GetDlgItemText(IDC_EDIT1,str);
	if (str.IsEmpty())
	{
		MessageBox(_T("请输入阶次"),_T("警告"));
		return;
	}
	int n = _ttoi(str);
	if (n<0)
	{
		MessageBox(_T("请输入大于1的阶数"),_T("警告"));
		return;
	}
	czy::Fit fit;
	fit.polyfit(m_x,m_yploy,n,true);
	CString strFun(_T("y=")),strTemp(_T(""));
	for (int i=0;i<fit.getFactorSize();++i)
	{
		if (0 == i)
		{
			strTemp.Format(_T("%g"),fit.getFactor(i));
		}
		else
		{
			double fac = fit.getFactor(i);
			if (fac<0)
			{
				strTemp.Format(_T("%gx^%d"),fac,i);
			}
			else
			{
				strTemp.Format(_T("+%gx^%d"),fac,i);
			}
		}
		strFun += strTemp;
	}
	str.Format(_T("方程:%s\r\n误差:ssr:%g,sse=%g,rmse:%g,确定系数:%g"),strFun
		,fit.getSSR(),fit.getSSE(),fit.getRMSE(),fit.getR_square());
	GetDlgItemText(IDC_EDIT,strTemp);
	SetDlgItemText(IDC_EDIT,strTemp+_T("\r\n------------------------\r\n")+str);
	//绘制拟合后的多项式
	std::vector<double> yploy;
	fit.getFitedYs(yploy);
	CChartLineSerie* pfitLineSerie1 = m_chartCtrl.CreateLineSerie();	
	pfitLineSerie1->SetSeriesOrdering(poNoOrdering);//设置为无序
	pfitLineSerie1->AddPoints(&m_x[0], &yploy[0], yploy.size());
	pfitLineSerie1->SetName(_T("多项式拟合方程"));//SetName的作用将在后面讲到
	pfitLineSerie1->SetWidth(2);
}

步骤如下:

  • 和线性拟合一样,声明Fit变量
  • 输入观察值,同时输入需要拟合的阶次,这里输入2阶,就是2项式拟合,最后的布尔变量是标定是否需要把拟合的结果点保存起来,保存点会根据观察的x值计算拟合的y值,保存结果点会花费更多的内存,如果拟合后需要绘制,设为true会更方便,如果只需要拟合的方程,可以设置为false
  • 拟合完后,拟合的相关结果保存在czy::Fit里面,可以通过相关方法调用,方法在头文件中都有详细说明
代码:
for (int i=0;i<fit.getFactorSize();++i)
{
	if (0 == i)
	{
		strTemp.Format(_T("%g"),fit.getFactor(i));
	}
	else
	{
		double fac = fit.getFactor(i);
		if (fac<0)
		{
			strTemp.Format(_T("%gx^%d"),fac,i);
		}
		else
		{
			strTemp.Format(_T("+%gx^%d"),fac,i);
		}
	}
	strFun += strTemp;
}

是用于生成方程的,由于系数小于时,打印时会把负号“-”显示,而正数时却不会显示正号,因此需要进行判断,如果小于0就不用添加“+”号,如果大于0就添加“+”号
结果如下:



源代码下载:

© 著作权归作者所有

共有 人打赏支持
尘中远
粉丝 1
博文 26
码字总数 47436
作品 0
朝阳
程序员
SP++3.0已发布,欢迎大家使用(同心协力,共创开源)

SP++ (Signal Processing in C++) 是一个关于信号处理与数值计算的开源C++程序库,该库提供了信号处理与数值计算中常用算法的C++实现。SP++中所有算法都以C++类模板方法实现,以头文件形式组...

张明
2011/02/12
0
55
Machine Learning 之Logistic回归算法中最小二乘法的Matlab曲线拟合

Machine Learning 之Logistic回归算法中最小二乘法的Matlab曲线拟合 逻辑回归是机器学习(Machine Learning)中常见的机器学习算法,在处理逻辑回归(Logistic Regression)离散数据点集时,...

zhangphil
2017/12/15
0
0
Python机器学习笔记——监督学习

使用fit() 函数,对训练数据构成的特征X 和标签 y进行训练,调用 predict() 函数,对未知分类样本分类。 K-近邻算法(KNN) 计算待分类数据点与已有数据集中所有数据点的距离。取距离最小的前...

BeerBread134
2017/08/08
0
0
Machine Learning 2 - 非线性回归算法分析

2017-08-02@erixhao 技术极客TechBooster AI 机器学习第二篇 - 非线形回归分析。我们上文深入本质了解了机器学习基础线性回归算法后,本文继续研究非线性回归。 非线性回归在机器学习中并非热...

erixhao
2017/08/07
0
0
7种回归模型

【编者按】回归分析是建模和分析数据的重要工具。本文解释了回归分析的内涵及其优势,重点总结了应该掌握的线性回归、逻辑回归、多项式回归、逐步回归、岭回归、套索回归、ElasticNet回归等七...

NORTHhan
2017/03/09
0
0
【译文】R语言非线性回归入门

作者 Lionel Hertzog   译者 钱亦欣       在一簇散点中拟合一条回归线(即线性回归)是数据分析的基本方法之一。有时,线性模型能很好地拟合数据,但在某些(很多)情形下,变量间的...

上大飞猪钱小莲
2017/01/08
0
0
9.机器学习sklearn-----岭回归及其应用实例

1.基本概念 对于一般地线性回归问题,参数的求解采用的是最小二乘法,其目标函数如下: 参数w的求解,也可以使用如下矩阵方法进行: 对于矩阵X,若某些列线性相关性较大(即训练样本中某些属...

bxg1065283526
04/22
0
0
【18-03-24】Matlab 数据分析

[18-03-24] Matlab 数据分析 1° 数据插值 一般地,从各种试验得来的数据总是有一定的数量,而利用插值技术能够从有限的数据中获取系统整体的状态,因此,数据插值在各行各业,特别是信号处理...

千阳Weston
03/25
0
0
机器学习之线性回归的最小二乘法求解

机器学习之线性回归的最小二乘法求解 假设现在一个普通的一阶线性方程,y=2x+2t。t是随机噪音,生成的散列点(x,y)会沿直线y=2*x上下摆动。利用最小二乘法做一次简单的一阶“曲线”拟合。用m...

zhangphil
2017/12/18
0
0
Matlab 线性拟合 & 非线性拟合

使用Matlab进行拟合是图像处理中线条变换的一个重点内容,本文将详解Matlab中的直线拟合和曲线拟合用法。 关键函数: fittype Fit type for curve and surface fitting Syntax ffun = fitty...

SibylY
2014/01/05
0
0

没有更多内容

加载失败,请刷新页面

加载更多

下一页

Hbase增删查改工具类

package cn.hljmobile.tagcloud.service.data.repository;import java.util.ArrayList;import java.util.HashMap;import java.util.List;import java.util.Map;import java.util......

gulf
8分钟前
0
0
详解机器学习中的梯度消失、爆炸原因及其解决方法

前言 本文主要深入介绍深度学习中的梯度消失和梯度爆炸的问题以及解决方案。本文分为三部分,第一部分主要直观的介绍深度学习中为什么使用梯度更新,第二部分主要介绍深度学习中梯度消失及爆...

tantexian
9分钟前
0
0
JavaMail 发送邮件

参考 https://www.cnblogs.com/xdp-gacl/p/4216311.html 发送html格式邮件 package com.example.stumgr;import java.util.Properties;import javax.mail.Message;import javax.mail......

阿豪boy
11分钟前
0
0
Mongodb安装教程

MongoDB是一个基于分布式文件存储的数据库,是一个介于关系数据库和非关系数据库之间的产品,是非关系数据库当中功能最丰富,最像关系数据库的。它支持的数据结构非常松散,是类似json的bso...

木筏笔歆
12分钟前
0
0
Hadoop之YARN命令

概述 YARN命令是调用bin/yarn脚本文件,如果运行yarn脚本没有带任何参数,则会打印yarn所有命令的描述。 使用: yarn [--config confdir] COMMAND [--loglevel loglevel] [GENERIC_OPTIONS] [...

舒运
13分钟前
0
0
个推数据统计产品(个数)iOS集成实践

最近业务方给我们部门提了新的需求,希望能一站式统计APP的几项重要数据。这次我们尝试使用的是个推(之前专门做消息推送的)旗下新推出的产品“个数·应用统计”,根据官方的说法,个推的数...

个推
14分钟前
0
0
Git 修改提交的用户名和邮箱名字

在通过git提交代码时,发现提交的用户名是自己mac的账户名,想要修改为其他名字和邮箱。 首先可以通过以下命令查看当前配置下的信息,包括用户名和邮箱: > git config --list 针对单项目的相...

edwardGe
17分钟前
0
0
Object.defineProperty()

Object.defineProperty(obj, props)方法直接在一个对象上定义新的属性或修改现有属性,并返回该对象。 obj 在其上定义或修改属性的对象 props 要定义其可枚举属性或修改的属性描述符的对象 ...

litCabbage
18分钟前
0
0
JEESZ分布式框架--单点登录集成方案(三)

多项目集成单点登录配置 当sso验证完成之后,客户端系统需要接收sso系统返回的结果时,需要定义一个过滤器获取返回结果,然后针对返回结果做相关处理.如果不需要做处理时,此处Filter也可以不...

明理萝
19分钟前
0
1
超简单的利用plist 查看ipa包名及其它信息

1.下载ipa安装包 2.用rar等工具打开 3.将iTunesMetadata.plist文件解压出来 4.用http://www.atool.org/plist_reader.php在线反编译工具 5.在其中中找到softwareVersionBundleId 就是包名...

xiaogg
19分钟前
0
0

没有更多内容

加载失败,请刷新页面

加载更多

下一页

返回顶部
顶部