clear all
tic
length=200;
o=0.2;
x=1:1:100;
y=1./(x.^3);
y1=y+o;
y2=y-o;
X=[x,x];
Y=[y1,y2];
d(1:100)=1;
d(101:200)=-1;
%求核矩阵
for p=1:length
for i=1:length
%k(p,i)=(X(p)*X(i)+Y(p)*Y(i)+1)^2;
k(p,i)=exp((-((X(p)-X(i))^2+(Y(p)-Y(i))^2))/0.5);
%Q(p,i)=d(p)*d(i)*k(p,i);
end
end
temp=d'*d.*k;
a=pinv(temp)*ones(length,1);
s=find(a>0,1,'first');
b = 1/d(s)-a'*(d'.*k(:,s));
%b=ones(length,1);
%a=zeros(length,1);
%a=inv(Q)*b;
%a=b/Q;
%for p=1:length
% if a(p)>0
% break;
% end
%end
%b=0;
%for i=1:length
% b=b-a(i)*d(i)*(X(i)*X(p)+Y(i)*Y(p)+1)^2;
%end
%b=b+1/d(p);
%syms s
for i=1:100
%for p=1:length
% m(p)=a(p)*d(p)*(X(p)*x(i)+Y(p)*s+1)^2;
%end
%m=@(s)a.*d'.*exp(-((X-x(i)).^2+(Y-s).^2)/0.5)'*ones(1,200)+b;
% l=sum(m)+b;
%solve('l','s')
yp(i)=fzero(@(r) a'.*d.*exp(-((X-x(i)).^2+(Y-r).^2)/0.5)*ones(200,1)+b, y(i) );
end
e=((y-yp).^2)*ones(100,1)
toc
figure
subplot(1,2,1);
plot(x,yp,'xr',x,y,'b');
legend('目标函数','拟合曲线');
xlabel('x轴','Fontsize',15);
ylabel('y轴','Fontsize',15);
h=legend('拟合曲线','实际曲线','location','best');
set(h,'Fontsize',15);
title('y=1/x^3的拟合曲线','Fontsize',15);
subplot(1,2,2);
%figure
xlabel('x轴','Fontsize',15);
ylabel('y轴','Fontsize',15);
plot(yp-y,'r');
title('1/x^3 svm拟合误差','Fontsize',15);
评论0