剧情提要：
[机器小伟]在[工程师阿伟]的陪同下进入了结丹初期的修炼，
这次要修炼的目标是[函数的应用]。
正剧开始：
星历2016年04月09日 11:38:05, 银河系厄尔斯星球中华帝国江南行省。
[工程师阿伟]正在和[机器小伟]一起研究[函数的应用]。
<span style="font-size:18px;">function funTask(x) {
return Math.log(x)/Math.log(Math.E)+2*x-6;
}
var r = 20;
config.setSector(1,1,1,1);
config.graphPaper2D(0, 0, r);
config.axis2D(0, 0,180, 1);
//坐标轴设定
var scaleX = 2*r, scaleY = 2*r;
var spaceX = 2, spaceY = 2;
var xS = -10, xE = 10;
var yS = -10, yE = 10;
config.axisSpacing(xS, xE, spaceX, scaleX, 'X');
config.axisSpacing(yS, yE, spaceY, scaleY, 'Y');
var transform = new Transform();
//存放函数图像上的点
var a = [], b = [], c = [], d = [], e = [];
//需要显示的函数说明
var f1 = 'y=lnx+2x-6', f2 = 'y = logx', f3 = 'y=x', f4 = 'y=e^x',
f5 = 'y = ln(x)';
//函数描点
for (var x = xS; x <= xE; x+=0.3) {
if (x != 0) {
a.push([x, funTask(x)]);
}
}
//二分法求函数的零点
//区间的最小值和最大值
var minX = 0, maxX = 4;
var y1 = y2 = 0;
x = minX;
y1 = funTask(x);
x = maxX;
y2 = funTask(x);
//如果在给定区间上存在有零点
if (y1 * y2 < 0) {
var epsilon = 0.000001;
while (Math.abs(y1-y2) > epsilon) {
x = minX;
y1 = funTask(x);
x = maxX;
y2 = funTask(x);
x = (minX+maxX)/2;
y = funTask(x);
if (y * y1 < 0) {
maxX = x;
}
else {
minX = x;
}
}
plot.setFillStyle('blue')
.fillText('零点：x = '+ ((minX+maxX)/2).toFixed(3), -200, -60, 120);
}
</span>
这个题，由于[人叫板老师]的精度不够，所以取得的值并不如小伟的准确。
<span style="font-size:18px;"> if (1) {
var r = 20;
config.setSector(1,1,1,1);
config.graphPaper2D(0, 0, r);
config.axis2D(0, 0,180, 1);
//坐标轴设定
var scaleX = 2*r, scaleY = 2*r;
var spaceX = 2, spaceY = 2;
var xS = -10, xE = 10;
var yS = -10, yE = 10;
config.axisSpacing(xS, xE, spaceX, scaleX, 'X');
config.axisSpacing(yS, yE, spaceY, scaleY, 'Y');
var transform = new Transform();
//存放函数图像上的点
var a = [], b = [], c = [], d = [], e = [];
//需要显示的函数说明
var f1 = 'y=2^x+3x-7', f2 = 'y = logx', f3 = 'y=x', f4 = 'y=e^x',
f5 = 'y = ln(x)';
//函数描点
for (var x = xS; x <= xE; x+=0.3) {
if (x != 0) {
a.push([x, funTask(x)]);
}
}
//二分法求函数的零点
//区间的最小值和最大值
var minX = 0, maxX = 2;
var y1 = y2 = 0;
x = minX;
y1 = funTask(x);
x = maxX;
y2 = funTask(x);
//如果在给定区间上存在有零点
if (y1 * y2 < 0) {
var epsilon = 0.000001;
while (Math.abs(y1-y2) > epsilon) {
x = minX;
y1 = funTask(x);
x = maxX;
y2 = funTask(x);
x = (minX+maxX)/2;
y = funTask(x);
if (y * y1 < 0) {
maxX = x;
}
else {
minX = x;
}
}
plot.setFillStyle('blue')
.fillText('零点：x = '+ ((minX+maxX)/2).toFixed(3), -200, -60, 120);
}
//document.write(a+'<br/>');
//存放临时数组
var tmp = [];
//显示变换
if (a.length > 0) {
a = transform.scale(transform.translate(a, 0, 0), scaleX/spaceX, scaleY/spaceY);
//函数1
tmp = [].concat(a);
shape.pointDraw(tmp, 'red');
tmp = [].concat(a);
shape.multiLineDraw(tmp, 'pink');
plot.setFillStyle('red');
plot.fillText(f1, 100, -90, 200);
}
if (b.length > 0) {
b = transform.scale(transform.translate(b, 0, 0), scaleX/spaceX, scaleY/spaceY);
//函数2
tmp = [].concat(b);
shape.pointDraw(tmp, 'blue');
tmp = [].concat(b);
shape.multiLineDraw(tmp, '#22CCFF');
plot.setFillStyle('blue');
plot.fillText(f2, 100, -120, 200);
}
if (c.length > 0) {
c = transform.scale(transform.translate(c, 0, 0), scaleX/spaceX, scaleY/spaceY);
tmp = [].concat(c);
shape.pointDraw(tmp, 'green');
tmp = [].concat(c);
shape.multiLineDraw(tmp, '#CCFF22');
plot.setFillStyle('green');
plot.fillText(f3, 100, -150, 200);
}
if (d.length > 0) {
d = transform.scale(transform.translate(d, 0, 0), scaleX/spaceX, scaleY/spaceY);
tmp = [].concat(d);
shape.pointDraw(tmp, 'orange');
tmp = [].concat(d);
shape.multiLineDraw(tmp, '#CC8800');
plot.setFillStyle('orange');
plot.fillText(f4, 100, -180, 200);
}
if (e.length > 0) {
e = transform.scale(transform.translate(e, 0, 0), scaleX/spaceX, scaleY/spaceY);
tmp = [].concat(e);
shape.pointDraw(tmp, 'purple');
tmp = [].concat(e);
shape.multiLineDraw(tmp, 'purple');
plot.setFillStyle('purple');
plot.fillText(f5, 100, -60, 200);
}
}
}
function funTask(x) {
return Math.pow(2, x)+3*x-7;
}</span>
<span style="font-size:18px;">function funTask(x) {
return Math.pow(x, 3)+1.1*Math.pow(x, 2)+0.9*x-1.4;
}
function funTask(x) {
return 3-Math.log(x)-x;
}
</span>
小伟决定来好好研究一下这个多项式拟合数据点。
这里有一组数据：
<span style="font-size:18px;">//一组点
if (1) {
var r = 20;
config.setSector(1,1,1,1);
config.graphPaper2D(0, 0, r);
config.axis2D(0, 0,180, 1);
//坐标轴设定
var scaleX = 2*r, scaleY = 2*r;
var spaceX = 2, spaceY = 2;
var xS = -10, xE = 10;
var yS = -10, yE = 10;
config.axisSpacing(xS, xE, spaceX, scaleX, 'X');
config.axisSpacing(yS, yE, spaceY, scaleY, 'Y');
var X=[ 1 ,2 ,3 ,4 ,5 ,6],
Y=[ 2.5 ,3.51 ,4.45 ,5.52 ,6.47 ,7.51];
var array = [];
var size = X.length;
for (var i = 0; i < size; i++) {
array.push([X[i], Y[i]]);
}
var transform = new Transform();
var tmp = [];
array = transform.scale(transform.translate(array, 0, 0), scaleX/spaceX, scaleY/spaceY);
tmp = [].concat(array);
shape.pointDraw(tmp, 'green');
//tmp = [].concat(array);
//shape.multiLineDraw(tmp, '#CCFF22');
}</span>
经过这样的拟合计算：
<span style="font-size:18px;">#一次多项式拟合
def linefit(x , y):
N = float(len(x))
sx,sy,sxx,syy,sxy=0,0,0,0,0
for i in range(0,int(N)):
sx += x[i]
sy += y[i]
sxx += x[i]*x[i]
syy += y[i]*y[i]
sxy += x[i]*y[i]
a = (sy*sx/N -sxy)/( sx*sx/N -sxx)
b = (sy - a*sx)/N
r = abs(sy*sx/N-sxy)/math.sqrt((sxx-sx*sx/N)*(syy-sy*sy/N))
return a,b,r
def tmp():
X=[ 1 ,2 ,3 ,4 ,5 ,6]
Y=[ 2.5 ,3.51 ,4.45 ,5.52 ,6.47 ,7.51]
a,b,r=linefit(X,Y)
print("X=",X)
print("Y=",Y)
print("拟合结果: y = %10.5f x + %10.5f , r=%10.5f" % (a,b,r) )
>>>
X= [1, 2, 3, 4, 5, 6]
Y= [2.5, 3.51, 4.45, 5.52, 6.47, 7.51]
拟合结果: y = 1.00000 x + 1.49333 , r= 0.99989</span>
来看一下拟合结果是否正确：
<span style="font-size:18px;"> if (1) {
var r = 20;
config.setSector(1,1,1,1);
config.graphPaper2D(0, 0, r);
config.axis2D(0, 0,180, 1);
//坐标轴设定
var scaleX = 2*r, scaleY = 2*r;
var spaceX = 2, spaceY = 2;
var xS = -10, xE = 10;
var yS = -10, yE = 10;
config.axisSpacing(xS, xE, spaceX, scaleX, 'X');
config.axisSpacing(yS, yE, spaceY, scaleY, 'Y');
var X=[ 1 ,2 ,3 ,4 ,5 ,6],
Y=[ 2.5 ,3.51 ,4.45 ,5.52 ,6.47 ,7.51];
var array = [];
var size = X.length;
for (var i = 0; i < size; i++) {
array.push([X[i], Y[i]]);
}
var transform = new Transform();
var tmp = [];
//绘制拟合曲线
var a = [], b = [], c = [], d = [], e = [];
//需要显示的函数说明
var f1 = 'y=x+1.49333';
//函数描点
for (var x = xS; x <= xE; x+=1) {
if (x != 0) {
a.push([x, funTask(x)]);
}
}
if (a.length > 0) {
a = transform.scale(transform.translate(a, 0, 0), scaleX/spaceX, scaleY/spaceY);
//函数1
tmp = [].concat(a);
shape.pointDraw(tmp, 'red');
tmp = [].concat(a);
shape.multiLineDraw(tmp, 'pink');
plot.setFillStyle('red');
plot.fillText(f1, 100, -90, 200);
}
//原数据点
array = transform.scale(transform.translate(array, 0, 0), scaleX/spaceX, scaleY/spaceY);
tmp = [].concat(array);
shape.pointDraw(tmp, 'green');
//tmp = [].concat(array);
//shape.multiLineDraw(tmp, '#CCFF22');
}
function funTask(x) {
return x+1.49333;
}</span>
好像是没错。
来看看二次多项式拟合。
<span style="font-size:18px;">def polyfit(x, y, degree):
import numpy;
results = {}
coeffs = numpy.polyfit(x, y, degree)
results['polynomial'] = coeffs.tolist()
# r-squared
p = numpy.poly1d(coeffs)
# fit values, and mean
yhat = p(x) # or [p(z) for z in x]
ybar = numpy.sum(y)/len(y) # or sum(y)/len(y)
ssreg = numpy.sum((yhat-ybar)**2) # or sum([ (yihat - ybar)**2 for yihat in yhat])
sstot = numpy.sum((y - ybar)**2) # or sum([ (yi - ybar)**2 for yi in y])
results['determination'] = ssreg / sstot #准确率
return results
def tmp():
x=[ 1 ,2 ,3 ,4 ,5 ,6]
y=[ 2.5 ,3.51 ,4.45 ,5.52 ,6.47 ,7.2]
z1 = polyfit(x, y, 2)
print(z1);
>>>
{'determination': 0.99893014166855909, 'polynomial': [-0.024285714285714542, 1.1257142857142872, 1.3699999999999986]}</span>
<span style="font-size:18px;"> if (1) {
var r = 20;
config.setSector(1,1,1,1);
config.graphPaper2D(0, 0, r);
config.axis2D(0, 0,180, 1);
//坐标轴设定
var scaleX = 2*r, scaleY = 2*r;
var spaceX = 2, spaceY = 2;
var xS = -10, xE = 10;
var yS = -10, yE = 10;
config.axisSpacing(xS, xE, spaceX, scaleX, 'X');
config.axisSpacing(yS, yE, spaceY, scaleY, 'Y');
var X=[ 1 ,2 ,3 ,4 ,5 ,6]
Y=[ 2.5 ,3.51 ,4.45 ,5.52 ,6.47 ,7.2]
var array = [];
var size = X.length;
for (var i = 0; i < size; i++) {
array.push([X[i], Y[i]]);
}
var transform = new Transform();
var tmp = [];
//绘制拟合曲线
var a = [], b = [], c = [], d = [], e = [];
//需要显示的函数说明
var f1 = 'y=-0.024x2+1.126x+1.370';
//函数描点
for (var x = xS; x <= xE; x+=1) {
if (x != 0) {
a.push([x, funTask(x)]);
}
}
if (a.length > 0) {
a = transform.scale(transform.translate(a, 0, 0), scaleX/spaceX, scaleY/spaceY);
//函数1
tmp = [].concat(a);
shape.pointDraw(tmp, 'red');
tmp = [].concat(a);
shape.multiLineDraw(tmp, 'pink');
plot.setFillStyle('red');
plot.fillText(f1, 100, -90, 200);
}
//原数据点
array = transform.scale(transform.translate(array, 0, 0), scaleX/spaceX, scaleY/spaceY);
tmp = [].concat(array);
shape.pointDraw(tmp, 'green');
//tmp = [].concat(array);
//shape.multiLineDraw(tmp, '#CCFF22');
}
</span>
这也不知道它怎么算的，明显有很大偏差啊。
这里顺便贴下对数拟合的代码，这个和小伟没任何关系，
因为不是小伟写的，一时也不好验证对错，先放着看看吧。
<span style="font-size:18px;">对数函数拟合-这个是最难的，baidu上都找不到，google了半天才找到的。指数、幂数拟合啥的，都用这个，把func改写一下就行
from scipy import log as log print pcov
import numpy
from scipy import log
from scipy.optimize import curve_fit
def func(x, a, b):
y = a * log(x) + b
return y
def polyfit(x, y, degree):
results = {}
#coeffs = numpy.polyfit(x, y, degree)
popt, pcov = curve_fit(func, x, y)
results['polynomial'] = popt
# r-squared
yhat = func(x ,popt[0] ,popt[1] ) # or [p(z) for z in x]
ybar = numpy.sum(y)/len(y) # or sum(y)/len(y)
ssreg = numpy.sum((yhat-ybar)**2) # or sum([ (yihat - ybar)**2 for yihat in yhat])
sstot = numpy.sum((y - ybar)**2) # or sum([ (yi - ybar)**2 for yi in y])
results['determination'] = ssreg / sstot
return results
x=[ 1 ,2 ,3 ,4 ,5 ,6]
y=[ 2.5 ,3.51 ,4.45 ,5.52 ,6.47 ,7.51]
z1 = polyfit(x, y, 2)
print z1
</span>
本节到此结束，欲知后事如何，请看下回分解。
时间： 2024-10-09 20:16:10