clc;clear;close all;
load R_shuju.mat
t=37:3:85;
y=151.4513./(1+(151.4513/69.1-1).*exp( -0.0484.*t)); %对总人口进行预测
t=0:25; %令1980年为初始年，25个数据一直到2005年
x=[98.705 100.1 101.654 103.008 104.357 105.851 107.5 109.3 111.026 112.704 114.333 115.823 117.171 118.517 119.85 121.121 122.389 123.626 124.761 125.786 126.743 127.627 128.453 129.227 129.988 130.756]; 
[c,d]=solve( 'c/(1+(c/98.705-1)*exp(-5*d))=105.851' , 'c/(1+(c/98.705-1)*exp(-8*d))=111.026' ,'c' , 'd' );%求初始参数
b0=[ 109.8216, - 0.19157]; %初始参数值
fun=inline( 'b(1)./(1+(b(1)/98.705-1).*exp(-b(2).*t))' , 'b' , 't' )
[b1,r1,j1]=nlinfit(t,x,fun,b0);
R1=r1.^2;
R2=(x-mean(x)).^2; 
R=1-R1/R2'; %可决系数
 W=sum(abs(r1)) ;%残差绝对值之和
%%画图
y= 153.5351./(1+(153.5351/98.705-1).*exp( -0.0477.*t)); %非线性拟合的方程
plot(t,x, '*' ,t,y, '-or' ) %对原始数据与曲线拟合后的值作图
gg=R_shuju(1,28:53)
set(gca,'XTick',0:1:26,'xticklabel',gg);
set(gca,'XTickLabelRotation',45)%使标注旋转角度
% for i=1:26
% text(t(i),y(i),num2str(y(i)))
% end
 %t表示是从1980年开始，t=20表示是2000年
% t=20:3:53;
% y=  153.5351./(1+(153.5351/98.705-1).*exp( -0.0477.*t)); %对总人口进行预测
figure(5);
%令2000年为初始年，30个数据一直到2030年
t2=20:1:50;
y2= 153.5351./(1+(153.5351/98.705-1).*exp( -0.0477.*t2)); %对总人口进行预测
plot(t2,y2, '*') %对原始数据与曲线拟合后的值作图
%%%%%%%%%生成一个年度数据刻度名称
nianfen=30;%需要使用到的所有年数
gg=cell(1,nianfen+1);
u=1;
chushinianfen=2000;%初始年份
for i=chushinianfen:chushinianfen+nianfen
    xx=[num2str(i),'年'];
    gg(1,u)=cellstr(xx)
    u=u+1
end


set(gca,'XTick',20:1:50,'xticklabel',gg);
set(gca,'XTickLabelRotation',45)%使标注旋转角度
for i=1:31
    text(t2(i),y2(i),num2str(y2(i)))
    % line(t2(i),y2(i),'linestyle',':','color','r');
end
grid on
hold on
% figure(6)
mm=(cell2mat(R_shuju(2,48:66)));
mm=double(mm/1000)
plot(20:1:38,mm,'-or' );