基于移动最小二乘法的曲线曲面拟合(python语言实现)
1.移动最小二乘法
上篇论文采用最小二乘法来拟合曲线,如果离散数据量比较大,形状复杂,还需要分段拟合和平滑化,因此采用移动最小二乘法进行曲线拟合,可以克服上面的缺点,还具有一些优点;
移动最小二乘法与传统的最小二乘法相比,有两个比较大的改进:
( 1)拟合函数的建立不同。这种方法建立拟合函数不是采用传统的多项式或其它函数,而是由一个系数向量 a(x)和基函数 p(x)构成, 这里 a(x)不是常数,而是坐标 x 的函数。
( 2)引入紧支( Compact Support)概念,认为点 x 处的值 y 只受 x
附近子域内节点影响,这个子域称作点 x 的影响区域, 影响区域外的节点对 x的取值没有影响。在影响区域上定义一个权函数w(x), 如果权函数在整个区域取为常数, 就得到传统的最小二乘法。参考自《基于移动最小二乘法的曲线曲面拟合-曾清红》
2.拟合函数的建立
在拟合区域的一个局部子域上, 拟合函数 f (x)表示为:
式中 为待求系数,它是坐标x的函数。称为基函数。它是一个k阶完备的多项式,m是基函数的项数,
对于一维问题 :
基函数可以为 p(x)=[1,x,x2,,,,xm][1,x,x2,,,,xm][1,x,x^2,,,,x^m]
二维问题可以为: 线性基 p(x)=[1.x.y]T[1.x.y]T[1.x.y]^T, m=3 二次基 p(x) =[1,x,y,x2,xy,y2]T[1,x,y,x2,xy,y2]T[1,x,y,x^2,xy,y^2]^T m=6
这是为在阅读文献时的疑惑,因为我解决的是一维问题,所以不需要二维的基函数。
在移动最小二乘近似中, 系数 ai(x)ai(x)a_i(x) 是通过令近似函数 u(x) 在点 x 的邻域 内各节点误差的加权平方和为最小来确定的
式中 n 为点 x 的邻域 内所包含的节点数.,w(x)=(x−xI)w(x)=(x−xI)w(x) = (x -x_I) 称为节点 xIxIx_I 处的权函数, 它在节点 xI 周围的一个有限区域中大于零, 而在该区域外为零 . 权函数的定义表明, 只有在节点 xI 的影响域范围内的节点才对该点的近似函数产生影响.
这里对支撑域进行说明:
如图:
将整个x范围划分为若干个区域,每个区域包含若干个x点,那么并且规定其中一点为标准点,其他点为参考点。
参考点与标准点的距离作为权函数的参数。得出权重。
3.权函数
权函数在移动最小二乘法中起着非常重要的作用。移动最小二乘法中的权函数 w(x−xI)w(x−xI)w(x-x_I)应该具有紧支性,也就是权函数在 x
的一个子域内不等于零, 在这个子域之外全为零, 这个子域称为权函数的支持域(即 x 的影响区域)。一般选择圆形作为权函数的支持域(见图其半径记为 smaxsmaxs_{max}。 由于权函数的紧支性,只有这些包含在影响区域内的数据点对点 x 的取值有影响权函数 w(x−xI)w(x−xI)w(x-x_I)应该是非负的,并且随着||x−xi||2||x−xi||2||x - x_i||^2 的增加单调递减。权函数还应具有一定的光滑性,因为拟合函数会继承权函数的连续性:如果权函数w(x−xI)w(x−xI)w(x-x_I)是 C1 阶连续的,则拟合函数也是 C1 阶连续的。常用的权函数是样条函数
4
3 法方程的推导
对于任意函数 h(x) 和 g(x), 引入记号:
那么:
公式4可以写为:
写成矩阵形式:
由上面的法方程, 解得 a(x).
然后求解得出A(x),B(x) 求解得出αα\alpha(x)
4.拟合流程
这里说明一下为什么要网格化,网格化主要是选取标准点,并以标准点来划分支撑域,确定支撑域半径和支撑域内的节点
x。
我仍然以上篇最小二乘法的数据点为例,通过代码编写移动最小二乘法的方法:
#主题部分
X=np.arange(-0.9,0.9,0.05)
# 数据点x个数
M=len(x)
# 基函数个数
N=2
p=np.zeros((M,2))
Y=[]
for XX in X:w = np.zeros((M,1))d=0.1 # 影响区域的半径for i in range(0,M):w[i]=W_fun(d, x[i], XX)p[i][0]=1p[i][1]=x[i]A=fun_A(x,w,p)B=fun_B(y,w,p)a=np.linalg.solve(A,B)Y.append(a[0]+a[1]*XX)----------#其他函数部分
# 权函数
def W_fun(d,x,X):s=abs(x-X)/dif (s<=0.5):return (2/3)-4*s**2+4*s**3elif(s<=1):return (4/3)-4*s+4*s**2-(4/3)*s**3else:return 0
# 权函数记号(pm,pn)的计算
def pm_pn(w,x,p,m,n):# x为数据点,w为支撑域的权重,M为数据点个数 p1,p2为传入的数值pmn=0M=len(x)# i代表数据点,m n代表(pm,pn)的下标for i in range(M):pmn=pmn+w[i]*p[i][m]*p[i][n]return float(pmn)
# B矩阵的建立
def fun_B(u,w,p):pumI=0M=len(u) #数据点个数m=p.shape[1] # 基函数个数B=[]for j in range(m):for i in range(M):pumI=pumI+w[i]*p[i][j]*u[i]B.append(float(pumI))return B# A矩阵的建立
def fun_A(x,w,p):M=len(x)函数m=p.shape[1]A=[]for mm in range(m):matA=[]for nn in range(m):pmn=pm_pn(w,x,p,mm,nn)matA.append(pmn)A.append(matA)return A
结果
5.结果
绿色为移动最小二乘法,红色为最小二乘法。
参考文献:
1基于移动最小二乘法的曲线曲面拟合-曾清红
2移动最小二乘法在多功能传感器数据重构中的应用-刘丹
3 移动最小二乘法(MLS)曲线曲面拟合C++代码实现
https://blog.csdn.net/liumangmao1314/article/details/54179526
基于移动最小二乘法的曲线曲面拟合(python语言实现)相关推荐
- 基于移动最小二乘法的曲线曲面拟合论文阅读笔记
基于移动最小二乘法的曲线曲面拟合论文阅读笔记 论文地址:http://www.cnki.com.cn/Article/CJFDTotal-GCTX200401016.htm 一.Problem Sta ...
- 基于强化学习的坦克大战python语言实现
这个项目是基于一个人工智能算法(基于全连接线性BP网络的增强学习模型)代理玩经典游戏坦克大战. 在个游戏中,机器控制的坦克主要目标是射击敌方坦克并保卫自家的根据地. 游戏中坦克的动作空 ...
- 移动最小二乘法(MLS)曲线曲面拟合C++代码实现
移动最小二乘法(MLS)曲线曲面拟合 曲线曲面拟合有很多种方法,Beizer,B样条等,曲面拟合移动最小二乘法是一个很好的选择,本文会详细讲解一下移动最小二乘法方法拟合曲面,并给出C++代码实现. 本 ...
- Python语言更适合低代码开发平台
今年低代码开发真的很火,钉钉.企业微信都相继推出了自己的低代码平台.为什么低代码平台这么火呢? 这也主要取决于企业管理系统的供求关系.有资料统计出,未来5年国内至少需要开发5亿个新应用,才能满足中国企 ...
- livechart 只显示 y 值_基于Python语言的SEGY格式地震数据读取与显示编程
敬请关注<地学新视野> 摘要:本文简单介绍了SEG-Y地震数据文件格式,以及如何用Python语言编写读写SEG-Y格式的地震数据并绘制地震剖面,其中用到了Segyio和matplotli ...
- 「实战案例」基于Python语言开发的信用评分卡
信用风险计量模型可以包括跟个人信用评级,企业信用评级和国家信用评级.人信用评级有一系列评级模型组成,常见是A卡(申请评分卡).B卡(行为模型).C卡(催收模型)和F卡(反欺诈模型). 今天我们展示的是 ...
- python语言deLong‘s test:通过统计学的角度来比较两个ROC曲线、检验两个ROC曲线的差异是否具有统计显著性
python语言deLong's test:通过统计学的角度来比较两个ROC曲线.检验两个ROC曲线的差异是否具有统计显著性 目录
- 如何使用Python语言将测量数据曲线动起来? matplotlib转换成 GIF文件
简 介: 本文就如何使用的Python中的Matplotlib绘制的图运动起来进行讨论.实际上是生成多个图片,将它们合成一个动图.编写了PlotGIF函数,将matplotlib绘制的图片存储,然后生 ...
- 基于python的界面自动化测试-基于Python语言的自动化测试实战第二章(上)
原标题:基于Python语言的自动化测试实战第二章(上) 测试环境搭建 2.1 Windows 下的环境搭建 如果想要学习一门编程语言,对于新手来说只需到其官方网站上去下载最新版本安装即可,但对于想要 ...
最新文章
- pandas使用nunique函数计算dataframe每个数据列的独特值的个数(count number of unique values in each column of dataframe)
- 10进制与16进制之间的转换 delphi
- 中国互联网大佬江湖拼什么?拼财力拼出身拼前景拼造势
- Goreplay来做流量回放
- java.lang.ClassNotFoundException: org.jaxen.JaxenException
- 安卓获取Java端连接
- sssp-springmvc+spring+spring-data-jpa问题总结
- Java基础,使用双循环嵌套,实现输出1-1000之间所有的素数,并统计有多少个
- 计算机专业可以转英语吗,计算机专业英语词汇转.doc
- 360浏览器Linux版内核,360安全浏览器国产操作系统版下载,附特色功能介绍
- python文本自动伪原创_现在有哪些好用的伪原创工具?
- 百度竞价后台操作技巧
- HTTP压测工具之wrk
- java api微盘_新浪微盘的API Python版封装
- 大数据治理包括哪几个方面
- Python + 批量裁剪图片
- kafka auto.offset.reset值失效 earliest/latest详解
- 快速搭建一个本地的FTP服务器
- abb式c语言,ABB机器人是用什么语言编程的? ——ABB机器人
- 【TA-霜狼_may-《百人计划》】美术2.7 Metallic 与 Speculer流程
热门文章
- maven install时报错 dependencies.dependency.version' for cn.e3mall:e3-manager-pojo:jar is missing.
- Matlab怎么计算信号的能量,Matlab小波包分解后如何求各频带信号的能量值? [转]...
- 简单servlet第三篇---使用servlet读取配置文件的内容
- 电磁场与磁路理论基础 (总结收藏)
- 安徽省16市什么企业适合做知识产权贯标申报好处和什么材料汇编整理
- 【数据库】广州Y行数据库工程师笔试题
- vulnhub View2akill
- 【PAT乙级】狼人杀-简单版(暴力逻辑)
- 有关The Last Promise的汉化①国外高手的FE7 hack版
- 经纬度标定及大地坐标系相关概念(一)