最佳拟合直线


Time Limit: 1000ms   Memory limit: 65536K  有疑问？点这里^_^

题目描述

在很多情况下，天文观测得到的数据是一组包含很大数量的序列点图象，每一点用x值和y值定义。这就可能需要画一条通过这些点的最佳拟合曲线。
为了避免只对个别数据分析，需要进行最佳曲线拟合。考虑N个数据点，它们的坐标是(X1,Y1)，(X2,Y2)...，(XN,YN)。假设这些值中的X是严格的精确值，Y的值是测量值(含有一些误差)。

         

对于一个给定的X，如X1，对应的值Y1与曲线C上对应的Y值将存在一个差值。我们用D1表示这个差值，有时我们也称这个差值为偏差、误差或残差，它可能是正、负或零。类似的，X2...，XN，对应的差值为D2，....，DN。
 我们用D12 + D22 + ... + DN2 作为衡量曲线C拟合的“最佳”程度，这个值越小越好，越大则越不好。因此，我们做以下定义：任何一种类型的曲线，它们都有一个共同的特性，当ΣDi2最小时，称为最佳拟合曲线。注：∑指“取和”计算。
 一条曲线具有这一特性时，称之为“最小二乘拟合”，这样的曲线称为“最小二乘曲线”。
本次的计算任务是拟合为一条直线，数学上称之为“线性回归”。“回归”一词看起来有点陌生，因为计算最佳曲线没什么好“回归”的，最好的术语就是“曲线似合”，在直线情况下就是“线性曲线拟合”。
你的任务是编写程序用最小二乘法计算出以下线性方程的系数(斜率a以及y轴的截距b)：
 y = a*x + b   (4.1)
 a和b可以使用以下公式计算：

 

式中N是数据点的个数。注意，以上两式具有相同的分母，∑指逐项加法计算(取和)。∑x指对所有的x值求和，∑y指对所以的y值求和，∑(x^2)指对所有x的平方求和。∑xy指对所有的积xy进行取和计算。应注意，∑xy 与 ∑x*∑y是不相同的(“积的和”与“和的积”是不同的)，同样(∑x)^2与∑(x^2)也是不相同的(“和的平方”与“平方的和”是不相同的)。



输入
 n组整数表示xi，yi ，期中|x|<=106，|y|<=106,
 n < 15
输出
 最佳拟合曲线参数a和b，a和b各占一行，a 和b精确到小数点后3位。
示例输入

4
1  6
2  5
3  7
4  10

示例输出

1.400
3.500

提示
 
来源
#include<stdio.h>
#include<math.h>
int n;
double qiu(double *x)
{
	double s=0;
	for(int i=0;i<n;i++)
	s+=x[i];
	return s; 
}
int main()
{
	double x[20],y[20],c[20],d[20],a1,b1,o;
	scanf("%d",&n);
	for(int i=0;i<n;i++)
	{
	 scanf("%lf%lf",&x[i],&y[i]);
	 c[i]=x[i]*y[i];
	 d[i]=x[i]*x[i];
	}
	a1=n*qiu(c)-qiu(x)*qiu(y);
	b1=qiu(y)*qiu(d)-qiu(x)*qiu(c);
	o=n*qiu(d)-pow(qiu(x),2);
	printf("%.3f\n",a1/o);
	printf("%.3f\n",b1/o);
	
}

需要对数据进行函数拟合，首先画一下二维散点图，目测一下大概的分布，
所谓正态分布，就是高斯分布，正态曲线是一种特殊的高斯曲线。
python的scipy.optimize包里的curve_fit函数来拟合曲线，当然还可以拟合很多类型的曲线。scipy.optimize提供了函数最小值(标量或多维)、曲线拟合和寻找等式的根的有用算法。
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import math
#单个高斯模型，如果曲线有多个波峰，可以分段拟合
def func(x, a,u, sig):
return a*np.exp(-(x - u) ** 2 / (2 * sig ** 2)) / (sig * math.sqrt(2 * math.pi))
#混合高斯模型，多个高斯函数相加
def func3(x, a1, a2, a3, m1, m2, m3, s1, s2, s3):
return a1 * np.exp(-((x - m1) / s1) ** 2) + a2 * np.exp(-((x - m2) / s2) ** 2) + a3 * np.exp(-((x - m3) / s3) ** 2)
#正弦函数拟合
#def fmax(x,a,b,c):
#    return a*np.sin(x*np.pi/6+b)+c
#fita,fitb=optimize.curve_fit(fmax,x,ymax,[1,1,1])
#非线性最小二乘法拟合
#def func(x, a, b,c):
#    return a*np.sqrt(x)*(b*np.square(x)+c)
#用3次多项式拟合，可推广到n次多项式，数学上可以证明，任意函数都可以表示为多项式形式
#f1 = np.polyfit(x, y, 3)
#p1 = np.poly1d(f1)
#yvals = p1(x)  #拟合y值
#也可使用yvals=np.polyval(f1, x)
拟合，并对参数进行限制，bounds里面代表参数上下限，p0是初始范围，默认是[1,1,1]
x=np.arange(1,206,1)
num = []<-自己的y值
numhunt = []<-自己的y值
y = np.array(num)
yhunt = np.array(numhunt)
popt, pcov = curve_fit(func3, x, y)
popthunt, pcovhunt = curve_fit(func, x, yhunt,p0=[2,2,2])
ahunt = popthunt[0]
uhunt = popthunt[1]
sighunt = popthunt[2]
a1 = popt[0]
u1 = popt[1]
sig1 = popt[2]
a2 = popt[3]
u2 = popt[4]
sig2 = popt[5]
a3 = popt[6]
u3 = popt[7]
sig3 = popt[8]
yvals = func3(x,a1,u1,sig1,a2,u2,sig2,a3,u3,sig3) #拟合y值
yhuntvals = func(x,ahunt,uhunt,sighunt) #拟合y值
print(u'系数ahunt:', ahunt)
print(u'系数uhunt:', uhunt)
print(u'系数sighunt:', sighunt)
#绘图
plot1 = plt.plot(x, y, 's',label='insect original values')
plot2 = plt.plot(x, yvals, 'r',label='insect polyfit values')
plot3 = plt.plot(x, yhunt, 's',label='predator original values')
plot4 = plt.plot(x, yhuntvals, 'g',label='predator polyfit values')
plt.xlabel('date')
plt.ylabel('Nightly catches log10(N+1)')
plt.legend(loc=4) #指定legend的位置右下角
plt.title('insect/predator')
plt.show()
下图是单个和多个高斯拟合图像
下图是多项式拟合图像
图例的位置可以自定义设置
lower left
upper center
lower right
center
upper left
center left
upper right
lower center
best
center right
right

原博文
2019-06-02 14:35 −
需要对数据进行函数拟合，首先画一下二维散点图，目测一下大概的分布， 所谓正态分布，就是高斯分布，正态曲线是一种特殊的高斯曲线。 python的scipy.optimize包里的curve_fit函数来拟合曲线，当然还可以拟合很多类型的曲线。scipy.optimize提供了函数最小值(标量或多维)、...
0
23537
相关推荐
2019-12-09 15:04 −
举个例子就知道了
class test():   def __init__(self, *a, **b):     print(a)     print(b)     print(b.get('test'))
tester = test(1, 2, 3, test='abc')
...
2019-12-20 15:12 −
装饰器无参和有参案例：
...
2019-12-20 15:20 −
* 接受参数类型为元组
In [16]: def func(a,*pargs):   ...:   print(a,pargs)   ...: In [17]: func(1,2)
1 (2,) In [18]: func(1,2,3,4)
1 (2, 3, 4)
*...
2019-12-25 09:48 −
原文稿：https://www.cnblogs.com/xiaoxiao-niao/p/7765627.html
每日更新的list 英文版本  https://peter.sh/experiments/chromium-command-line-switches/
1 --allow-o...
2019-12-05 15:20 −
#!/usr/bin/python  #导入sys模块 import sys  print("the first argv:",sys.argv[0],"\n") print("the second argv:",sys.argv[1],"\n") print("the third argv:",s...
2019-12-08 10:15 −
需求:
我要用python通过api,启动这个job,并且启动这个job需要1个参数
安装依赖:
pipenv install python-jenkins
熟悉API的使用方法:
了解一个API的最好办法,就是先看它的 官方文档,
从中找到自己需要的方法 ,弄明白要给方法...
2019-12-20 16:36 −
一、python参数问题：1）位置参数（必需参数）：调用函数时根据函数定义的参数位置来传递参数，调用时的数量必须和声明时的一样： eg：  　def fun(name,sex):    pass  　fun(name1,sex1) #两个参数的顺序必须是一一对应的，自然参数个数也是要一样的了，
2）...
2019-12-21 21:29 −
**1.给定验证码长度n，生成随机验证码，验证码由数字、字母组成（参考chr()内置方法）**
程序代码如下： ```python
import random
def create_check_code(n):   check_code = '' # 创建一个变量用以接收随机出来的验证码字符   f...
2019-11-13 13:16 −
在 python函数的声明和调用 中我们简单的了解了函数的相关使用，然而在函数传递参数的时候，我们埋下了一个坑，关于不定长参数的传递我们还没有讲，今天这篇文章主要就是讲解这个问题。
一.函数不定长参数语法
函数的不定长参数包含*args...

1.将线性回归得到的两个参数打印出来，并且在图像中画出该拟合曲线。

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd


data=pd.read_csv( "finalData.csv", thousands=',')
data.plot(kind='scatter', x="GDP per capita", y='Life satisfaction', figsize=(5,3))


plt.axis([0, 60000, 0, 10])
X=np.linspace(0, 60000, 1000)
plt.plot(X, t0 + t1*X, "b")
plt.text(5000, 3.1, r"$\theta_0 = 4.85$", fontsize=14, color="b")
plt.text(5000, 2.2, r"$\theta_1 = 4.91 \times 10^{-5}$", fontsize=14, color="b")

plt.show()


from sklearn import linear_model
lin1 = linear_model.LinearRegression()
Xsample = np.c_[data["GDP per capita"]]
ysample = np.c_[data["Life satisfaction"]]
lin1.fit(Xsample, ysample)
# 输出 intercept_是截距, coef_是斜率系数
t0, t1 = lin1.intercept_[0], lin1.coef_[0][0]
print('打印t0',t0)
print('打印t1',t1)


2.结果展示



输出 intercept_是截距, coef_是斜率系数 ，这些也被成为回归系数。