rand('state',1)
n= 200;
apha=2.5;
e=1e-7;
delta=1e-3;

u1=randn(n+2,1);
u=u1;
mx=abs(max(u));
mi=abs(min(u));
mx=max(mi,mx);
u=2*u/mx;
y1=zeros(n+2,1);
trdata=zeros(100,3);
testy=zeros(100,1);

getx=zeros(100,3);
gety=zeros(100,1);
for i=3:n
    y1(i)=1.5*y1(i-1)*y1(i-2)/(1+y1(i-1)*y1(i-1)+y1(i-2)*y1(i-2))+0.35*sin(y1(i-1)+y1(i-2))+1.2*u(i-1);
    %y1(i)=(0.8-0.5*exp(-y1(i-1)*y1(i-1)))*y1(i-1)-(0.3+0.9*exp(-y1(i-1)*y1(i-1)))*y1(i-2)+u(i);
end    
mx=abs(max(y1));
mi=abs(min(y1));
mx=max(mi,mx);
y=2*y1/mx;
y=y1(3:end);
u=u(3:end);

for i=1:100
    trdata(i,1)=y1(i+1); 
    trdata(i,2)=y1(i+2); 
    trdata(i,3)=u(i);
    testy(i)=y1(i+3);
end
for i=1:100
    getx(i,1)=y1(i+100);
    getx(i,2)=y1(i+101);
    getx(i,3)=u(i+99);
    gety(i)=y1(i+102);
end
[n,m]=size(trdata);
h=zeros(n,n+1);
h(1:n,1)=1;

%y2=zeros(n,1);
for i=1:n
    for j=2:n+1
        dbtmp=norm(trdata(i,1:m)-trdata(j-1,1:m));
        h(i,j)=exp(-0.5*dbtmp*dbtmp/(apha*apha));
    end
end    
sigma=e;
beta=inv(sigma*eye(n+1)+h'*h)*h'*testy;
pbeta=beta;
stime=cputime;
for i=1:10000
    htmp=norm(testy-h*beta);
    sigma=htmp*htmp/n;
    ss=diag(abs(beta));
    beta=ss*inv(sigma*eye(n+1)+ss*h'*h*ss)*ss*h'*testy;
    if i>1
        if norm(beta-pbeta)/norm(pbeta)<delta
        %if sigma<delta
            break
        end    
    end 
    pbeta=beta;
end    
fprintf('Execution time: %4.1f seconds\n',cputime - stime);
hold on;
ind = find((abs(beta)>1e-4));
k=size(ind);
for i=1:k
    if ind(i)>1
        %if y(ind(i))>0.5
            %plot(trdata(ind(i)-1), testy(ind(i)-1), 'rs');
        %else
        %    plot(synthtr(ind(i)-1,1), synthtr(ind(i)-1,2), 'rx');
        %end    
    end    
end    
for i=1:100
    dbtmp=0;
    for j=1:k
        mm=ind(j);
        if mm==1
            dbtmp=dbtmp+beta(mm);
        else
            tmpp=norm(getx(i,1:m)-trdata(mm-1,1:m));
            dbtmp=dbtmp+beta(mm)*exp(-0.5*tmpp*tmpp/(apha*apha));
        end    
    end
    y2(i)=dbtmp;
end    
%plot(x,y,'-k*',x1,y1,'-b',x2,y2,'-r');  
x=1:100;
plot(x,gety,'-*',x,y2,'-.rs');
xlabel('次數');
ylabel('Y(k)');
legend('原始數據','RVM預測')
Error_RVM_single=(gety(1:end)-y2(1:end)');
MAX_relative_rvm=mean(abs(Error_RVM_single))
f2=figure
plot(x,Error_RVM_single);





[nsv, sbeta, bias,H,svi]  = svr(trdata,testy,'rbf',100,'einsensitive',0.02);
%check the effect
%let's see the approximation
%reconstruct the identification result
TST_target= H(:,svi)*sbeta(svi) +bias;
for k=1:100
    XT=getx(k,:);
    svm(k) = svroutput(trdata,XT,'rbf',sbeta,bias);
end

svm_error=sqrt(mse(svm'-gety))
rvm_error=sqrt(mse(y2'-gety))



%Now check the effect
f3=figure;
plot(1:50,gety(1:50),'-o',1:50,y2(1:50),'-.rs',1:50,svm(1:50),':kd')
    xlabel('次數');
    ylabel('值');
    legend('原始數據','稀疏貝葉斯辨識','支持向量機辨識')
Error_SVM_single=(gety(1:end)-svm(1:end)');
MAX_relative_svm=mean(abs(Error_SVM_single))

f4=figure;
plot(1:100,Error_SVM_single,'b-',1:100,Error_RVM_single,'r.-')
legend('one step forecasting result','original data')
xlabel('Time Series');
ylabel('Pressure');