restart;
with(plots):
with(StringTools):
with(LinearAlgebra):
with(DEtools):
##############################################################################################
Region:=`Mosobl`; url:="https://gogov.ru/covid-19/mo#data";
#valp := [12.0545284599666, 40487.7036090492, 0.0878540675828671, 0.119034763474605, 0.208571881538399, 0.163876229453231, #0.0804863589181470, 0.0940935501988715, 0.105057419991328, 0.104244846871936, 0.00441107479765352, 0.00987259141378349, #0.0000931839501779284];
valp:=readdata(cat(Region,`3c.txt`));
##############################################################################################
fdisplay:=proc(f,p)
print(cat(f,`.jpg`)); #print(cat(f,`.eps`));
plotsetup(jpeg,plotoutput=cat(f,`.jpg`),plotoptions=`noborder`); print(display(p));
plotsetup(ps,plotoutput=cat(f,`.eps`),plotoptions=`noborder`); print(display(p));
plotsetup(default,plotoptions=`noborder`): print(display(p));
end:
pr:=proc(x) print(x); x; end:
grad:=(F,V)->map(q->diff(F,q),V):
linsplit:=(F,V)->subs(map(q->q=0,V),[op(grad(F,V)),F]):
corr:=proc(x,y) local i; seq(x[i]=y[i],i=1..nops(x)): end:
ssum:=(F,V)->convert([seq(F,V)],`+`):
pprod:=(F,V)->convert([seq(F,V)],`*`):
Lag:=proc(t,tx,kx) local i,j;
ssum(kx[i]*pprod(piecewise(j=i,1,(t-tx[j])/(tx[i]-tx[j])),j=1..nops(tx)),i=1..nops(tx)):
end:
Lag(t,[ta,tb],[a,b]); Lag(t,[ta,tb,tc],[a,b,c]);
pi:=evalf(Pi);
gM:=evalf(solve((1-x)^2=x,x)[2]):
goldMin:=proc(f,T,epsilon) local a,b,c,d,fa,fb,fc,fd,k;
a:=op(1,T); b:=op(2,T); fa:=f(a); fb:=f(b); k:=0;
c:=a+(b-a)*gM; fc:=f(c); d:=b-(b-a)*gM; fd:=f(d);
while abs(a-b)>epsilon do: k:=k+1;
if fc>fd then a:=c; fa:=fc; c:=d; fc:=fd; d:=b-(b-a)*gM; fd:=f(d);
else b:=d; fb:=fd; d:=c; fd:=fc;+ c:=a+(b-a)*gM; fc:=f(c);
fi;
od: #print(k);
(a+b)/2;
end:
findMin1:=proc(F,V) local f,df,f0,f1,f2,V0,V1,V2,ff,t,dt,i,j;
ff:=V->F(op(evalf(map(exp,V)))); V1:=evalf(map(ln,V)); f1:=F(op(V));
f:=[seq(F(seq(evalf(exp(V1[j]+piecewise(j=i,0.0001,0))),j=1..nops(V))),i=1..nops(V))];
df:=[seq((f[j]-f1)/0.1,j=1..nops(V))];
V0:=V1-0.001*df; f0:=ff(V0); V2:=V1+0.001*df; f2:=ff(V2);
dt:=0.0001; while f0<f1 do: V2:=V1; f2:=f1; V1:=V0; f1:=f0; V0:=V0-dt*df; f0:=ff(V0); dt:=dt*1.1; od;
dt:=0.0001; while f2<f1 do: V0:=V1; f0:=f1; V1:=V2; f1:=f2; V2:=V2+dt*df; f2:=ff(V2); dt:=dt*1.1; od;
t:=goldMin(t->ff(t*V0+(1-t)*V2),0..1,0.0001);
map(exp,t*V0+(1-t)*V2);
end:
findMin:=proc(F,V) local V1,Z1,Z2;
Z2:=pr(F(op(V))); V1:=findMin1(F,V); Z1:=pr(chi2(op(V1)));
while abs(1-Z1/Z2)>0.0001 do; Z2:=Z1; V1:=findMin1(F,V1); Z1:=pr(chi2(op(V1))); end;
V1;
end:
`====================================`;
`VERHULST FItAING`;
`====================================`;
f_:=d->sum(a[j]*d^j,j=0..n); fe_:=d->sum(a[j]*d^j,j=0..ne);
M:='M':
ff:=x->M*(1-1/(exp(x)+1)); ff_:=unapply(solve(y=ff(x),x),y); diff(ff_(x),x); dff_:=unapply(simplify(%,x),x);
ffe:=x->exp(x); ffe_:=unapply(solve(y=ffe(x),x),y); diff(ffe_(x),x); dffe_:=unapply(simplify(%,x),x);
sigma:=x->simplify(sqrt(x));
chi2:=(T,f_)->simplify(sum(evalf(ff_(T[k])-f_(k))^2/dff_(T[k])^2/sigma(T[k])^2,k=1..nops(T)));
chi2e:=(T,f_)->simplify(sum(evalf(ffe_(T[k])-f_(k))^2/dffe_(T[k])^2/sigma(T[k])^2,k=1..nops(T)));
F:=proc(T,chi2,f_) chi2(T,f_);
indets(%); grad(%%,%); subs(solve(%,%%),f_(i)); unapply(%,i);
end:
dig:={"0","1","2","3","4","5","6","7","8","9","0"}: val:=proc() global data,i; local j,f; f:=0;
while not(data[i] in dig) or f=1 and data[i] in {"+"} union dig do:
if f=1 and not(data[i] in dig) then f:=0; else if data[i]="+" then f:=1; fi fi; i:=i+1: od:
j:=i; while (data[i] in dig or data[i] in {"-","+"}) do i:=i+1: od: parse(data[j..i-1]);
end:
` `; Region; status,data,headers:=HTTP:-Get(url): HTTP:-Code(status); i:=Search("<th>",data):
iter:=proc() global i; local r;
r:=val(); if data[i]<>"." then NULL else [r,val(),val(),val(),val(),val()],iter() fi;
end:
[iter()]: tA:=[seq(%[nops(%)+1-i],i=1..nops(%))];
dd:=tA[1][1]+piecewise(tA[1][2]=2,-29,tA[1][2]=4,31,0)-1;
T:=map(q->q[4],tA): #writedata(Region || `-i.txt`,%): # \320\230\320\275\321\204\320\270\321\206\320\270\321\200\320\276\320\262\320\260\320\275\320\275\321\213\321\205
T3:=map(q->q[5],tA): #writedata(Region || `-m.txt`,%): # \320\243\320\274\320\265\321\200\321\210\320\270\321\205
T1:=map(q->q[6],tA): #writedata(Region || `-r.txt`,%): # \320\222\321\213\320\267\320\264\320\276\321\200\320\276\320\262\320\265\320\262\321\210\320\270\321\205
T2:=[seq(T[i]-(T1[i]+T3[i]),i=1..nops(T))]: #writedata(Region || `-h.txt`,%): # \320\221\320\276\320\273\321\214\320\275\321\213\321\205
i:='i':
Region; 'T'=T; 'T1'=T1; 'T2'=T2; 'T3'=T3;
nops(T); [i+dd $ i=1..%];
h:=x->x;
[seq(h(T[i])-h(T[i-1]),i=2..nops(T))]; [seq(%[i]-%[i-1],i=2..nops(%))]; [seq(%[i]-%[i-1],i=2..nops(%))];
[seq([i+dd+1,%%%[i]],i=1..nops(%%%))]: [seq([i+dd+2,%%%[i]],i=1..nops(%%%))]: [seq([i+dd+3,%%%[i]],i=1..nops(%%%))]:
display(
plot([%%%,%%,%],style=point),
plot([%%%,%%,%],legend=[`\320\277\320\265\321\200\320\262\320\260\321\217`,`\320\262\321\202\320\276\321\200\320\260\321\217`,`\321\202\321\200\320\265\321\202\321\214\321\217`]),
title=`\320\240\320\260\320\267\320\275\320\276\321\201\321\202\320\270 \321\200\321\217\320\264\320\260 N[i]`,titlefont=[roman,15],gridlines=true
);
[seq((h(T[i])-h(T[i-5]))/5.,i=6..nops(T))]: [seq((%[i]-%[i-3])/3.,i=4..nops(%))]: [seq((%[i]-%[i-3])/3.,i=4..nops(%))]:
[seq([i+dd+2,%%%[i]],i=1..nops(%%%))]: [seq([i+dd+4,%%%[i]],i=1..nops(%%%))]: [seq([i+dd+6,%%%[i]],i=1..nops(%%%))]:
display(
plot([%%%,%%,%],style=point),
plot([%%%,%%,%],legend=[`\320\277\320\265\321\200\320\262\320\260\321\217`,`\320\262\321\202\320\276\321\200\320\260\321\217`,`\321\202\321\200\320\265\321\202\321\214\321\217`]),
title=`\320\241\320\263\320\273\320\260\320\266\320\265\320\275\320\275\321\213\320\265 \321\200\320\260\320\267\320\275\320\276\321\201\321\202\320\270 \321\200\321\217\320\264\320\260 N[i]`,titlefont=[roman,15],gridlines=true
);
h:=x->ln(x);
[seq(h(T[i])-h(T[i-1]),i=2..nops(T))]; [seq(%[i]-%[i-1],i=2..nops(%))]; [seq(%[i]-%[i-1],i=2..nops(%))];
[seq([i+dd+1,%%%[i]],i=1..nops(%%%))]: [seq([i+dd+2,%%%[i]],i=1..nops(%%%))]: [seq([i+dd+3,%%%[i]],i=1..nops(%%%))]:
display(
plot([%%%,%%,%],style=point),
plot([%%%,%%,%],legend=[`\320\277\320\265\321\200\320\262\320\260\321\217`,`\320\262\321\202\320\276\321\200\320\260\321\217`,`\321\202\321\200\320\265\321\202\321\214\321\217`]),
title=`\320\240\320\260\320\267\320\275\320\276\321\201\321\202\320\270 \321\200\321\217\320\264\320\260 ln(N[i])`,titlefont=[roman,15],gridlines=true
);
[seq((h(T[i])-h(T[i-5]))/5.,i=6..nops(T))]: [seq((%[i]-%[i-3])/3.,i=4..nops(%))]: [seq((%[i]-%[i-3])/3.,i=4..nops(%))]:
[seq([i+dd+2,%%%[i]],i=1..nops(%%%))]: [seq([i+dd+4,%%%[i]],i=1..nops(%%%))]: [seq([i+dd+6,%%%[i]],i=1..nops(%%%))]:
display(
plot([%%%,%%,%],style=point),
plot([%%%,%%,%],legend=[`\320\277\320\265\321\200\320\262\320\260\321\217`,`\320\262\321\202\320\276\321\200\320\260\321\217`,`\321\202\321\200\320\265\321\202\321\214\321\217`]),
title=`\320\241\320\263\320\273\320\260\320\266\320\265\320\275\320\275\321\213\320\265 \321\200\320\260\320\267\320\275\320\276\321\201\321\202\320\270 \321\200\321\217\320\264\320\260 ln(N[i])`,titlefont=[roman,15],gridlines=true
);
h:=x->ln(x);
[seq(h(T[i])-h(T[i-1]),i=2..nops(T))]; [seq(%[i]-%[i-1],i=2..nops(%))]; [seq(%[i]-%[i-1],i=2..nops(%))];
[seq([i+dd+1,%%%[i]],i=1..nops(%%%))]: [seq([i+dd+1,%%%[i]],i=1..nops(%%%))]: [seq([i+dd+1,%%%[i]],i=1..nops(%%%))]:
display(
plot([%%%,%%,%],style=point),
plot([%%%,%%,%],legend=[`\320\277\320\265\321\200\320\262\320\260\321\217`,`\320\262\321\202\320\276\321\200\320\260\321\217`,`\321\202\321\200\320\265\321\202\321\214\321\217`]),
title=`\320\240\320\260\320\267\320\275\320\276\321\201\321\202\320\270 \321\200\321\217\320\264\320\260 ln(N[i])`,titlefont=[roman,15],gridlines=true
);
n:=1: ne:=n: 'f(t)'=Sum(a[j]*t^j,j=0..n);
fM:=proc(x) global M,chi2,F,T,f_; M:=x; chi2(T,F(T,chi2,f_)); end:
` `; `Approximation of the infection schedule by the solution of the Verhulst equation`; ` `;
M:=goldMin(fM,max(T)+2..max(T)*2,1);
nu:=F(T,chi2,f_): f:=unapply(ff(%(t)),t): N(t)=%(t); Chi2:=chi2(T,%%%);
cat(`Next day forecast: `,round(f(nops(T)+1)));
cat(`The level of 0.5 M is reached at `,round(1+fsolve(f(d)=0.5*M,d=30)+dd-31),` apr`);
cat(`The level of 0.85 M is reached at `,round(1+fsolve(f(d)=0.85*M,d=30)+dd-31),` apr`);
` `; `Approximation of the infection schedule by solving the Malthus equation`; ` `;
nue:=F(T,chi2e,f_): fe:=unapply(ffe(%(t)),t): N(t)=%(t);
simplify([diff(nu(d-dd),d),diff(nue(d-dd),d)]): [coeff(%[1],d,i) $ i=0..n-1];
plot(%%,d=1+dd..nops(T)+dd,view=[0..nops(T)+dd,0..0.5],legend=[`\320\244\320\265\321\200\321\205\321\216\320\273\321\214\321\201\321\202`,`\320\234\320\260\320\273\321\214\321\202\321\203\321\201`],
linestyle=[solid,dash],title=`\320\232\320\276\321\215\321\204\321\204\320\270\321\206\320\270\320\265\320\275\321\202 \320\267\320\260\321\200\320\260\320\266\320\265\320\275\320\270\321\217`,titlefont=[roman,20],labels=[t,alpha(t)],gridlines=true);
d1:=fsolve(f(d)=0.5*M,d=30)+dd; K_:=M; alpha_:=coeff(nu(t),t,1);
n:=4: ne:=n: 'f(t)'=Sum(a[j]*t^j,j=0..n);
fM:=proc(x) global M,chi2,F,T,f_; M:=x; chi2(T,F(T,chi2,f_)); end:
` `; `Approximation of the infection schedule by the solution of the Verhulst equation`; ` `;
M:=goldMin(fM,max(T)+2..max(T)*2,1);
nu:=F(T,chi2,f_): f:=unapply(ff(%(t)),t): N(t)=%(t); Chi2:=chi2(T,%%%);
cat(`Next day forecast: `,round(f(nops(T)+1)));
cat(`The level of 0.5 M is reached at `,round(1+fsolve(f(d)=0.5*M,d=30)+dd-31),` apr`);
cat(`The level of 0.85 M is reached at `,round(1+fsolve(f(d)=0.85*M,d=30)+dd-31),` apr`);
` `; `Approximation of the infection schedule by solving the Malthus equation`; ` `;
nue:=F(T,chi2e,f_): fe:=unapply(ffe(%(t)),t): N(t)=%(t);
[seq([i,(
(T[i-dd]-T[i-dd-1]) /(T2[i-dd]+T2[i-dd-1]) /((1-T[i-dd]/M)+(1-T[i-dd-1]/M))
)*4],i=1+dd+1..nops(T)+dd)]: [seq([%[i][1],(%[i-1][2]+%[i][2]+%[i+1][2])/3],i=2..nops(%)-1)]:
Palpha:=display(plot([%],color=blue),plot([%],style=point,symbolsize=8,symbol=solidcircle,color=blue)):
simplify([diff(nu(d-dd),d),diff(nue(d-dd),d)]): [coeff(%[1],d,i) $ i=0..n-1];
plot(%%,d=1+dd..nops(T)+dd,view=[0..nops(T)+dd,0..0.5],legend=[`\320\244\320\265\321\200\321\205\321\216\320\273\321\214\321\201\321\202`,`\320\234\320\260\320\273\321\214\321\202\321\203\321\201`],
linestyle=[solid,dash],title=`\320\232\320\276\321\215\321\204\321\204\320\270\321\206\320\270\320\265\320\275\321\202 \320\267\320\260\321\200\320\260\320\266\320\265\320\275\320\270\321\217`,titlefont=[roman,20],labels=[t,alpha(t)],gridlines=true):
display(Palpha,%);
df:=unapply(diff(f(i),i),i): ddf:=unapply(diff(f(i),i,i),i):
display(
plot([[i+dd,T[i]] $ i=1..nops(T)],style=point,symbolsize=10,symbol=solidcircle),
plot(fe(i-dd),i=1+dd..max(90,dd+nops(T)),color=magenta),
plot(f(i-dd),i=1+dd..max(90,dd+nops(T))),
seq(plot([[i+dd,T[i]+3*sqrt(T[i])],[i+dd,T[i]-3*sqrt(T[i])]],color=blue),i=1..nops(T)),
axis[2]=[mode=log],
view=[1..80,1..M*1.1],labels=[t,N(t)],gridlines=true
);
display(
plot([[i+dd,T[i]] $ i=1..nops(T)],style=point,symbolsize=8,symbol=solidcircle),
plot(fe(i-dd),i=1+dd..max(120,dd+nops(T)),color=magenta),
plot(f(i-dd),i=1+dd..max(120,dd+nops(T))),
# seq(plot([[i+dd,T[i]+3*sqrt(T[i])],[i+dd,T[i]-3*sqrt(T[i])]],color=blue),i=1..nops(T)),
axis[2]=[mode=log],
view=[1..nops(T)+dd+1,1..T[nops(T)]*1.1],labels=[t,N(t)],gridlines=true
);
display(
plot([[i+dd,T[i]] $ i=1..nops(T)],style=point,symbolsize=10,symbol=solidcircle),
plot(fe(i-dd),i=1+dd..max(120,dd+nops(T)),color=magenta),
plot(f(i-dd),i=1+dd..max(dd+nops(T),90)),
plot(10*df(i-dd),i=1+dd..max(dd+nops(T),120),color=black),
plot(100*ddf(i-dd),i=1+dd..max(dd+nops(T),120),color=gray),
seq(plot([[i+dd,T[i]+3*sqrt(T[i])],[i+dd,T[i]-3*sqrt(T[i])]],color=blue),i=1..nops(T)),
view=[1..80,-M*0.3..M*1.1],labels=[t,N(t)],gridlines=true
);
display(
plot([[i+dd,T[i]] $ i=1..nops(T)],style=point,symbolsize=8,symbol=solidcircle),
plot(fe(i-dd),i=1+dd..max(120,dd+nops(T)),color=magenta),
plot(f(i-dd),i=1+dd..max(dd+nops(T),120)),
# seq(plot([[i+dd,T[i]+3*sqrt(T[i])],[i+dd,T[i]-3*sqrt(T[i])]],color=blue),i=1..nops(T)),
view=[1..nops(T)+dd+1,1..T[nops(T)]*1.1],labels=[t,N(t)],gridlines=true
);
dT:=[[i,(T[i-dd]-f(i-dd))/sigma(f(i-dd))] $ i=1+dd..dd+nops(T)]:
display( plot(%), plot(%,style=point,symbolsize=8,symbol=solidcircle),title=`\320\224\320\265\320\262\320\270\320\260\321\206\320\270\321\217`,titlefont=[roman,20] );
`====================================`;
`FORECAST`;
`====================================`;
proc3:=proc(E)
E[1]*convert(map(X->X^coeff(E[2],X,1),M),`*`);
end:
proc2:=proc(X,E)
proc3(E)*(coeff(E[3],X,1)-coeff(E[2],X,1));
end:
proc1:=proc(X)
convert(map(E->proc2(X,E),L),`+`);
end:
A:='A': B:='B': C:='C': M:=[A,B,C];
L:=[
[P[`01`],0,A],
[(B/K)*P[`12`],A,B],
[P[`23`],B,C],
[P[`10`],A,0], [P[`20`],B,0], [P[`30`],C,0]
]: Matrix(%);
eqs:=map(X->Diff(X,t)=proc1(X),M); Vector(%);
v:=M; alpha:='alpha': K:=k0; tA:=[1,15,35,50,58,62,73,nops(T)+dd]; kA:=['k1x||i' $ i=1..nops(tA)];
par:=[d0,k0,op(kA),k2a,k2b,k3];
param:=[
P[`01`]=0, P[`12`]=alpha(t,op(kA)), P[`23`]=beta(t,k2a,k2b),
P[`10`]=0, P[`20`]=k3
];
init:=[ A(-d0)=K, B(-d0)=1, C(-d0)=0 ];
res:=solve(map(rhs,eqs[1..2]),v[1..2]); res:=res[2]: subs(P[`30`]=P[`10`],param,res);
J:=Matrix(subs(res,map(q->grad(rhs(q),v[1..2]),eqs[1..2]))); evalm(%-lambda): collect(Determinant(%),lambda);
subs(P[`30`]=P[`10`],pr(param),%); solve(%,{lambda});
Eqs:=subs(map(q->q=q(t),v),Diff=diff,P[`30`]=P[`10`],param,eqs); #dsolve(%);
N:='N': A:='A': B:='B': C:='C': val:=valp:
#alpha:=unapply(simplify(evalf(piecewise(t<tA[1],kA[1],t<tA[2],Lag(t,tA[1..3],kA[1..3]),
# seq(op([t<tA[i+1],(Lag(t,tA[i-1..i+1],kA[i-1..i+1])+Lag(t,tA[i..i+2],kA[i..i+2]))/2]),i=2..nops(kA)-2),
#t<tA[nops(tA)],Lag(t,tA[nops(tA)-2..nops(tA)],kA[nops(kA)-2..nops(kA)]),
#kA[nops(kA)]))),t,op(kA));
alpha:=unapply(simplify(evalf(piecewise(t<tA[1],kA[1],t<tA[3],Lag(t,tA[1..4],kA[1..4]),
seq(op([t<tA[i+1],Lag(t,tA[i-1..i+2],kA[i-1..i+2])]),i=3..nops(kA)-3),
t<tA[nops(tA)],Lag(t,tA[nops(tA)-3..nops(tA)],kA[nops(kA)-3..nops(kA)]),
kA[nops(kA)]))),t,op(kA));
beta:=(t,k2a,k2b)->piecewise(t<70,k2a,k2b);
EQS:=[op(Eqs),op(init)]:
res:=dsolve(EQS,numeric,map(q->q(t),v),output=listprocedure,parameters=par); assign('v[i]=subs(res,v[i](t))' $ i=1..nops(v)):
chi2a:='chi2a': chi2:=unapply(chi2a(x0,xx,kA,x2a,x2b,x3),x0,xx,op(kA),x2a,x2b,x3):
chi2a:=proc(x0,xx,x1,x2a,x2b,x3) local i; global K_; K_:=xx;
res(parameters=[corr(par,[x0,xx,op(x1),x2a,x2b,x3])]):
sum((T[i]-(K_-A(i+dd)))^2/(K_-A(i+dd)),i=1..nops(T))+
sum((T2[i]-B(i+dd))^2/B(i+dd),i=1..nops(T2))+
sum((T1[i]-C(i+dd))^2/C(i+dd),i=1..nops(T1));
end:
chi2(op(pr(val))); val:=findMin(chi2,val); chi2(op(%));
#plot(map(q->q(t),v),t=0..3.0e4,legend=[`\320\275\320\265\320\270\320\275\321\204\320\270\321\206\320\270\321\200\320\276\320\262\320\260\320\275\320\275\321\213\320\265`,`\320\261\320\276\320\273\321\214\320\275\321\213\320\265`,`\320\277\320\265\321\200\320\265\320\261\320\276\320\273\320\265\320\262\321\210\320\270\320\265`],
#linestyle=[solid,dash,dashdot],gridlines=true);
writedata(cat(Region,`3c.txt`),val);
display(
plot(map(q->q(t),v),t=0..300,legend=[`\320\275\320\265\320\270\320\275\321\204\320\270\321\206\320\270\321\200\320\276\320\262\320\260\320\275\320\275\321\213\320\265`,`\320\261\320\276\320\273\321\214\320\275\321\213\320\265`,`\320\277\320\265\321\200\320\265\320\261\320\276\320\273\320\265\320\262\321\210\320\270\320\265`],
linestyle=[solid,dash,dashdot],gridlines=true),
plot([[seq([i+dd,K_-T[i]],i=1..nops(T))]],style=point,symbolsize=7,symbol=asterisk),
plot([[seq([i+dd,T1[i]],i=1..nops(T1))]],style=point,symbolsize=7,symbol=circle),
plot([[seq([i+dd,T2[i]],i=1..nops(T2))]],style=point,symbolsize=7,symbol=diamond,color=black),
size=[1000,400],legendstyle=[font=[roman,15]]
): fdisplay(cat(Region,`3c`),%);
logplot(map(q->q(t),v),t=0..500,legend=[`\320\275\320\265\320\270\320\275\321\204\320\270\321\206\320\270\321\200\320\276\320\262\320\260\320\275\320\275\321\213\320\265`,`\320\261\320\276\320\273\321\214\320\275\321\213\320\265`,`\320\277\320\265\321\200\320\265\320\261\320\276\320\273\320\265\320\262\321\210\320\270\320\265`],
linestyle=[solid,dash,dashdot],gridlines=true,size=[1000,400],legendstyle=[font=[roman,15]]);
display(
plot([
[[i,(T[i-dd]-(K_-A(i)))/sigma(K_-A(i))] $ i=1+dd..dd+nops(T)],
[[i,(T2[i-dd]-(B(i)))/sigma(B(i))] $ i=1+dd..dd+nops(T)],
[[i,(T1[i-dd]-(C(i)))/sigma(C(i))] $ i=1+dd..dd+nops(T)]
],linestyle=[solid,dash,dashdot],legend=[`\320\275\320\265\320\270\320\275\321\204\320\270\321\206\320\270\321\200\320\276\320\262\320\260\320\275\320\275\321\213\320\265`,`\320\261\320\276\320\273\321\214\320\275\321\213\320\265`,`\320\277\320\265\321\200\320\265\320\261\320\276\320\273\320\265\320\262\321\210\320\270\320\265`]),
plot([
[[i,(T[i-dd]-(K_-A(i)))/sigma(K_-A(i))] $ i=1+dd..dd+nops(T)],
[[i,(T2[i-dd]-(B(i)))/sigma(B(i))] $ i=1+dd..dd+nops(T)],
[[i,(T1[i-dd]-(C(i)))/sigma(C(i))] $ i=1+dd..dd+nops(T)]
],style=point,symbolsize=8,symbol=solidcircle),
size=[1000,300],legendstyle=[font=[roman,15]]
): fdisplay(cat(Region,`3c-dev`),%);
plot([v[1](t),v[2](t),t=0..3.0e4],axis[1]=[mode=log],axis[2]=[mode=log],labels=[v[1],v[2]],labelfont=[roman,15],gridlines=true);
[seq([i,(
(T[i-dd]-T[i-dd-1]) /(T2[i-dd]+T2[i-dd-1]) /((1-T[i-dd]/K_)+(1-T[i-dd-1]/K_))
)*4],i=1+dd+1..nops(T)+dd)]: [seq([%[i][1],(%[i-1][2]+%[i][2]+%[i+1][2])/3],i=2..nops(%)-1)]:
Palpha:=display(plot([%],color=blue),plot([%],style=point,symbolsize=8,symbol=solidcircle,color=blue)):
#display(%,gridlines=true,labels=['t','alpha(t)'],labelfont=[roman,15],view=[0..nops(T)+dd,0..0.9]);
subs(corr(par,val),alpha(t, op(kA)));
plot(%,t=-20..100,gridlines=true,labels=['t',''alpha(t)''],labelfont=[roman,15],view=[-20..100,0..0.24]):
fdisplay(cat(Region,`3c-zar`),%); display([Palpha,%],title=`\320\232\320\276\321\215\321\204\321\204\320\270\321\206\320\270\320\265\320\275\321\202 \320\267\320\260\321\200\320\260\320\266\320\265\320\275\320\270\321\217`,titlefont=[roman,20]);