#############Median ##Preparations ##Moments of K~Bin(m,r/sqrt(m)) Binmoment:=proc(i) local KKi,p,x,KKis; KKi:=x^i*binomial(n,x)*p^x*(1-p)^(n-x); KKis:=collect(subs(p=r/n^(1/2),simplify(subs((-1/(-1+p))^n=1/(1-p)^n,simplify(sum(KKi,x=0..n))))),n); RETURN(eval(KKis)) end; Binmoment(4); ### (A.19)--(A.21) for i to 4 do K[i] := Binmoment(i); od; #transformation and backtransformation of m,j,k: # +to allow MAPLE to work with rational functions: m -> m2^2 (or m2 = m^(1/2)) # +to tell MAPLE that j,k are O((sqrt(m)): j -> j m2, k -> k m2 # +to ease readability: j -> k/2 - j1 transf:=proc(x) subs(j=k/2-j1,subs(m=m2^2,u=u/m2,j1=j1*m2,j=j*m2,k=k*m2,x)); end; btransf:=proc(x) subs(m2=m^(1/2),subs(u=u*m2,j1=j1/m2,k=k/m2,x)); end; ##Asymptotics for Binomialcoefficients --- gamma n j k ### see (A.36)--(A.38) asbinom:=proc(upbinom,downbinom,i,pot,withprint) local j1,upb,dpb1,dpb2,n,n2,sqnS,sqnnom,sqndenom1,sqndenom2,sqn0,sqn,sqnN,SQN,bSQN,sqnN0,SQN0,bSQN0; # ## sqrt-Term in stirling approximation sqnS:=subs(n2=sqrt(n),convert(asympt((n2^2)!*exp(n2^2)/((n2^2)^(n2^2)),n2,i+2),polynom)); # #transformation of m,j,k: upb:=transf(upbinom); dpb1:=transf(downbinom); dpb2:=upb-dpb1; # #plugging in sqnS: sqnnom:=asympt(collect(subs(n=upb,sqnS),m2),m2,i+2); sqndenom1:=asympt(collect(subs(n=dpb1,sqnS),m2),m2,i+2); sqndenom2:=asympt(collect(subs(n=dpb2,sqnS),m2),m2,i+2); # #calculate the sqrt-Term x factors before integrand * sqrt(2pi) sqn0:=collect(simplify(asympt(collect(sqnnom/sqndenom1/sqndenom2,m2),m2,i+2)),m2); sqn:=asympt((2*Pi)^(1/2)*(upb+k*m2+1)^(pot/2)*(upb+1)*sqn0,m2,i+2); #control: if withprint=1 then print(btransf(asympt(sqn,m2,i))); fi; #gives common factor of 2^(5/2)*m^(3/2) -- to be separated from rest sqnN0:=asympt(sqn/4*(1/2)^(1/2)*(1/m2)^3,m2,i); sqnN:=asympt(sqn,m2,i); # SQN:=collect(simplify(collect(simplify(sqnN),m2)),m2); SQN0:=collect(simplify(collect(simplify(sqnN0),m2)),m2); #control: backtransformation if withprint=1 then bSQN:=btransf(SQN); print(bSQN); bSQN0:=btransf(SQN0); print(bSQN0); fi; # RETURN(eval(convert(SQN,polynom))) end; #Tests #asbinom(2*m-k,m-j,4,2,1); #asbinom(2*m-1-k,m-j,4,2,1); #asbinom(2*m-1-k,m-j-1,4,2,1); #asbinom(2*m-k,m-j,4,1,1); ##calculation of asmpytotic fractions: asfrac:=proc(nom,denom,i,withtransf,withprint) local an,ad,frn,FRN,bFRN; # #with transformation into m2, j1, k... ? if withtransf=1 then an:=transf(nom); ad:=transf(denom); else an:=nom; ad:=denom; fi; # frn:=collect(collect(asympt(collect(convert(an/ad,polynom),m2),m2,i),k),m2); FRN:=collect(simplify(collect(simplify(frn),m2)),m2); #control: backtransformation if withprint=1 then bFRN:=btransf(FRN); print(bFRN); fi; # RETURN(eval(FRN)) end; ###Tests (compare (A.34)) #asfrac(m-j,2*m-k,6,1,1); #collect(asfrac((2*m-k)^2,2*(m-j),6,1,0)+asfrac((2*m-k)^2,2*(m+j-k),6,1,0),m2); #btransf(collect(asfrac((2*m-k)^2,2*(m-j),6,1,0)+asfrac((2*m-k)^2,2*(m+j-k),6,1,0),m2)); #collect(asfrac((2*m-k)^3,3*(m-j)^2,9,1,0)-asfrac((2*m-k)^3,3*(m+j-k)^2,9,1,0),m2); #btransf(collect(asfrac((2*m-k)^3,3*(m-j)^2,9,1,0)-asfrac((2*m-k)^3,3*(m+j-k)^2,9,1,0),m2)); ###quantile ##Determining xnjk: # #asquantile:=proc(Fx,nom,denom,i,withprint) asquantile:=proc(Fx,i,withprint,withctrl) local xnjk1,xnjk0,t,G,G1,Gt,G1t,ddfn,Xnjk,ctrl,ctrl0,bxnjk,cx; # # one step Newton: # F(xnjk)!=Fx; # F(xnjk)=1/2+f0*xnjk+f1*xnjk^2/2+f2*xnjk^3/6 # G(xnjk)=1/2+f0*xnjk+f1*xnjk^2/2+f2*xnjk^3/6-Fxnjk # xnjk0=(Fx-1/2)/f0 # xnjk=xnjk0-G(xnjk0)/G'(xnjk0) # xnjk0:=(Fx-1/2)/f0; G:=1/2+f0*t+f1*t^2/2+f2*t^3/6-Fx; G1:=f0+f1*t+f2*t^2/2; Gt:=convert(asympt(collect(simplify(subs(t=xnjk0,G)),m2),m2,4),polynom); G1t:=convert(asympt(collect(simplify(subs(t=xnjk0,G1)),m2),m2,2),polynom); ddfn:=convert(asfrac(eval(Gt),eval(G1t),i,0,0),polynom); xnjk1:=xnjk0-ddfn; Xnjk:=convert(collect(collect(simplify(asympt(xnjk1,m2,i)),j1),m2),polynom); if withctrl=1 then ##control G(xnjk0), G(Xnjk) #G(xnjk0) cx:=collect(simplify(collect(subs(t=xnjk0,G),m2)),m2); ctrl0:=collect(simplify(asympt(cx,m2,i)),m2); print(ctrl0); ### O(m^-1) #G(Xnjk) cx:=collect(simplify(collect(subs(t=Xnjk,G),m2)),m2); ctrl:=collect(simplify(asympt(cx,m2,i+1)),m2); print(ctrl); ### O(m^-2) fi; # if withprint=1 then bxnjk:=btransf(Xnjk); print(bxnjk); fi; # RETURN(eval(Xnjk)) end; ##Tests # xnjk:=asquantile(Fxnjk,m-j,2*m-k,4,0); # xnjknX:=asquantile(FxnjknX,3,0,0); # for control purposes: f(xnjk) asfxnjk:=proc(x,i,withprint) local fxnjk0,fxnjk,bfnjk; fxnjk0:=collect(collect(f0+f1*x+f2*x^2/2,f0),m2); fxnjk:=collect(collect(simplify(collect(asympt(fxnjk0,m2,i),f0)),f0),m2); if withprint=1 then bfnjk:=btransf(fxnjk); print(bfnjk); fi; RETURN(eval(fxnjk)) end; ##Tests ## see (A.35) and display below #Fxnjk:=convert(asfrac(m-j,2*m-k-1,4,1,0),polynom); #xnjk:=asquantile(Fxnjk,m-j,2*m-k,4,1); Fxnjk:=convert(asfrac(m-j,2*m-k-1,4,1,0),polynom); xnjk:=asquantile(Fxnjk,4,1,1); fxnjk:=asfxnjk(xnjk,3,1); ##Taylor for log and exp logx:=x-x^2/2+x^3/3-x^4/4; expx:=1+x+x^2/2+x^3/6+x^4/24; ##for the terms in the exponential: dfLogf:=proc(nom,denom,i,withprint) local dflog,dflog1,blog,dnl,nl; nl:=transf(nom); dnl:=transf(denom); dflog1:=collect(subs(x=nl/dnl*DF,logx),m2); dflog:=collect(collect(collect(convert(asympt(collect(simplify(collect(dnl*dflog1,k)),m2),m2,i),polynom),k),m2),DF); if withprint=1 then blog:=btransf(dflog); print(blog); fi; RETURN(eval(dflog)) end; # ##Tests #dfLp:=dfLogf(2*m-k-1,m-j,9,1); #dfLm:=dfLogf(-(2*m-k-1),m-k-1+j,9,1); # #dfexpn:=collect(convert(asympt(collect(collect(collect(dfLp+dfLm,k),j1),m2),m2,2),polynom),DF); ##dF (later: =dfn) #expansion of F(t) around 0: ## calculates hyjkn (5.50) asdf:=proc(xx,yy,z,i,withprint) local dfn0,dfn1,dfn2,dfn3,dfn,dF1,dF2,dF3,dF4, dfexpl,exptu,expu,expy,expty1,expty2,a,sigman; # dfn0:=1/2+f0*a+f1*a^2/2+f2*a^3/6; # t=u+xnjk # tell MAPLE u ~ O(1/sqrt(n)) (up to log....) dfn:=asympt(collect(collect(subs(a=u/m2+xx,dfn0)-yy,u),m2),m2,i); #print(dfn); ## poowers of dfn in minimal order of approx: dF1:=collect(collect(collect(dfn,m2),k),u); dF2:=collect(collect(convert(asympt(dF1^2,m2,i+2),polynom),k),u); dF3:=collect(collect(convert(asympt(dF1^3,m2,i+2),polynom),k),u); dF4:=collect(collect(convert(asympt(dF1^4,m2,i+2),polynom),k),u); # dfexpl:=collect(collect(z,k),DF); #print(dfexpl); # exptu:=convert(asympt(collect(subs(DF^2=dF2,DF^3=dF3,DF^4=dF4,dfexpl),m2),m2,i+2),polynom); # sigman:=sqrt(8)*m2*f0; #substitute u=y*sqrt(n)/sigman expty1:=collect(collect(subs(u=y*m2/sigman,exptu),k),m2); # substract -y^2/2 (for Gaussian density!) and reduce order expty2:=convert(asympt(collect(expty1+y^2/2,m2),m2,i),polynom); #print(expty1);print(subs(x=expty2,expx)); # the rest out of the exponent by taylor of exp #to avoid too large terms do it by hand: expty22:=convert(asympt(collect(expty2^2/2,m2),m2,i),polynom); expty23:=convert(asympt(collect(expty2^3/6,m2),m2,i-1),polynom); expty24:=convert(asympt(collect(expty2^4/24,m2),m2,i-2),polynom); expy:=asympt(collect(collect(1+expty2+expty22+expty23+expty24,j1),m2),m2,i); #print(expty3); # if withprint=1 then bexpt:=btransf(expy); print(bexpt); fi; RETURN(eval(expy)) end; ##Tests (see (A.40); P is just the polynomial for factor m2^{-1}) #hyjkn:=asdf(xnjk,Fxnjk,dfexpn,3,1); #btransf(collect(collect(collect(convert(asympt(hyjkn,m2,3),polynom),j1),y),m2)); # expand f(t)/sigman: [(A.39)] asft:=proc(xn, i,withprint) local ft1,ft2,ft,bft,sigman; sigman:=sqrt(8)*m2*f0; ft1:=(f0+f1*t+f2*t^2/2)/sigman; # substitute according to ... ft2:=collect(collect(subs(u=y/sigman,subs(t=xn+u,ft1)),y),m2); # reduce the order of approx to necessary degree ft:=convert(asympt(ft2,m2,i),polynom); if withprint=1 then bft:=btransf(ft); print(bft); fi; RETURN(eval(ft)) end; # expand G(t)=t^2: ## for displays after (A.43) astint:=proc(xn,G, i,withprint) local t21,t2,bt2,sigman; sigman:=sqrt(8)*m2*f0; # substitute according to ... t21:=collect(collect(subs(u=y/sigman,subs(t=xn+u,G)),y),m2); # reduce the order of approx to necessary degree t2:=convert(asympt(t21,m2,i),polynom); if withprint=1 then bt2:=btransf(t2); print(bt2); fi; RETURN(eval(t2)) end; ##Tests #aft:=asft(xnjk,3,1); #atint:=astint(xnjk,t^2,5,1); # determines integrand getasintegrand:=proc(binup,binlow,i,Gi,pot,withprint,withprintall) local inte,inte0, xnjkn,Fxnjkn,cnjkn,aftn,atintn,hyn,binte,dfLpn,dfLmn,dfexpN,hyn0; ### (A.43) cnjkn:=asbinom(binup,binlow,i-1,pot,withprintall); Fxnjkn:=convert(asfrac(binlow,binup,i,1,withprintall),polynom); xnjkn:=asquantile(Fxnjkn,i,0,withprintall); # print(xnjkn); #dfs: # dfLpn:=dfLogf(binup,binlow,i+6,withprintall); dfLmn:=dfLogf(-binup,binup-binlow,i+6,withprintall); # dfexpN:=collect(convert(asympt(collect(collect(collect(dfLpn+dfLmn,k),j1),m2),m2,2),polynom),DF); # # hyn0:=asdf(xnjkn,Fxnjkn,dfexpN,i+1,withprintall);print(hyn0); hyn0:=asdf(xnjkn,Fxnjkn,dfexpN,i+1,withprintall); hyn:=convert(asympt(hyn0,m2,i),polynom); aftn:=asft(xnjkn,i+1,withprintall); atintn:=astint(xnjkn,Gi,i+2,withprintall); if withprint=2 then # print(cnjkn); print(hyn); print(aftn); print(atintn); fi; inte0:=cnjkn*hyn*aftn*atintn*(2*f0)^pot; inte:=collect(collect(collect(collect(convert(asympt(collect(inte0,m2),m2,i),polynom),f0),k),j1),m2); if withprint=1 then binte:=btransf(inte); print(binte); fi; RETURN(eval(inte)) end; ##Tests #aa:=getasintegrand(2*m-k-1,m-j,3,t^2,2,2,0); #intjk:=convert(asympt(getasintegrand(2*m-k-1,m-j,3,t^2,2,1,0),m2,3),polynom); ## => gives P1njk P2njk #intkp:=collect(collect(subs(j1=k/2,intjk),y),m2); #intkm:=collect(collect(subs(j1=-k/2,intjk),y),m2); # Q2 coefficient to m2^{-1} of [cf (A.45)] ## subs(y1 = y, subs(y = 0, subs(y^8 = y1^8, y^6 = y1^6, y^4 = y1^4, y^2 = y1^2, collect(collect(intkp, y), m2)))); # integrates out y inteyout:=proc(ing,withprint) local tt1,tt2,tt3,tt4,tt,btt; tt1:=subs(y^8=105,collect(collect(ing,m2),y)); tt2:=subs(y^6=15,tt1); tt3:=subs(y^4=3,tt2); tt4:=subs(y^2=1,tt3); tt:=collect(collect(subs(y=0,tt4),k),m2); if withprint=1 then btt:=btransf(tt); print(btt); fi; RETURN(eval(tt)) end; ##Tests #int0kp:=inteyout(intkp,1); #int0km:=inteyout(intkm,1); # integrates out k intekout:=proc(nn,mm,ing,i) local j0,ing0,ing1,ing2,ing3,ing4,ing5,ing6,KK; ing0:=collect(subs(k=k/m2,ing),k); for j0 to 4 do KK := subs(n=nn,Binmoment(5-j0)); ing0:=subs(k^(5-j0)=KK,ing0); od; ing1:=convert(asympt(collect(ing0,m2),m2,i+6),polynom); ing2:=subs(m2=mm^(1/2),ing1); ing3:=convert(asympt(collect(ing2,n),n,i+2),polynom); ing4:=convert(asympt(collect(subs(n=n2^2,ing3),n2),n2,i),polynom); ing5:=collect(collect(ing4,r),f0); ing6:=collect(collect(collect(subs(n2=n^(1/2),collect(ing5,n2)),r),f0),n); RETURN(eval(ing6)) end; ##Tests #intp:=intekout(2*m2^2,n/2,int0kp,3,1); #intp:=intekout(2*m2^2+1,(n-1)/2,int0kp,3,1); #intm:=intekout(2*m2^2+1,(n-1)/2,int0km,3,1); getasrisk:=proc(binup,bindown,i,Gi,pot,pm,withprintl,asyj,withprintall) local intjkn,intk0,intk1,erg; intjkn:=convert(asympt(getasintegrand(binup,bindown,i,Gi,pot,withprintall,withprintall),m2,i),polynom); # decide how to choose j: if withprintl=1 then print(asympt(intjkn,m2,asyj)); fi; if pm=1 then intk0:=collect(collect(subs(j1=k/2,intjkn),y),m2); else intk0:=collect(collect(subs(j1=-k/2,intjkn),y),m2); fi; intk1:=inteyout(intk0,withprintall); if transf(binup+k+1)=2*m2^2+1 then erg:=intekout(2*m2^2+1,(n-1)/2,intk1,3,0); else erg:=intekout(2*m2^2,n/2,intk1,3,0); fi; RETURN(eval(erg)); end; getasrisk(2*m-k,m-j,3,t^2,2,1,1,1,0); getasrisk(2*m-k,m-j,3,t^2,2,-1,1,1,0); getasrisk(2*m-k-1,m-j,3,t^2,2,1,1,1,0); getasrisk(2*m-k-1,m-j,3,t^2,2,-1,1,1,0); getallrisks:=proc(binup,bindown,i,pm,withprintall,withr0,withf00) local asMSE,asBias,asBias2,asVar; asMSE:=getasrisk(binup,bindown,i,t^2,2,pm,0,0,withprintall); asBias:=getasrisk(binup,bindown,i,t,1,pm,0,0,withprintall); asBias2:=collect(collect(collect(convert(subs(n2=n^(1/2),asympt(subs(n=n2^2,asBias^2),n2,3)),polynom),r),f0),n); asVar:=collect(collect(collect(asMSE-asBias2,r),f0),n); if withf00=1 then if withr0=1 then print(subs(f1=0,r=0,asMSE));print(subs(f1=0,r=0,asBias));print(subs(f1=0,r=0,asBias2));print(subs(f1=0,r=0,asVar)); else print(subs(f1=0,asMSE));print(subs(f1=0,asBias));print(subs(f1=0,asBias2));print(subs(f1=0,asVar)); fi; else if withr0=1 then print(subs(r=0,asMSE));print(subs(r=0,asBias));print(subs(r=0,asBias2));print(subs(r=0,asVar)); else print(asMSE);print(asBias);print(asBias2);print(asVar); fi; fi; end; getallrisks(2*m-k,m-j,3,1,0,1,1); getallrisks(2*m-k-1,m-j,3,-1,0,0,0); #################################### #### Midpoint estimator #################################### #transformation and backtransformation of m,j,k,u, and [new!] s: # +to allow MAPLE to work with rational functions: m -> m2^2 # +to tell MAPLE that j,k are O((sqrt(m)): j -> j m2, k -> k m2 # +to ease readability: j -> k/2 - j1 transfs:=proc(x) collect(subs(j=k/2-j1,subs(m=m2^2,u=u/m2,j1=j1*m2,j=j*m2,k=k*m2,s=s/m2,x)),m2); end; btransfs:=proc(x) subs(m2=m^(1/2),collect(subs(u=u*m2,s=s*m2,j1=j1/m2,k=k/m2,x),m2)); end; ## Calculation Delta(s,u) getDelta:=proc(withprint) local Pu,PP1,Delta0,Delta,bDelta; #expansion for the cdf Pu:=1/2+f0*u+f1*u^2/2+f2*u^3/6; #expansion for the 1/cdf PP1:=convert(taylor(1/Pu,u,4),polynom); Delta0:=(f0*(s-u)+f1*(s^2-u^2)/2+f2*(s^3-u^3)/6)*PP1; Delta:=convert(asympt(collect(transfs(Delta0),m2),m2,6),polynom); bDelta:=btransfs(Delta); if withprint=1 then print(bDelta); fi; RETURN(eval(Delta)) end; getlogDelta:=proc(withprint,asy) local DE, logx,logdel,blogdel,mj; mj:=transf(m-j-1); DE:=getDelta(0); logx:=x-x^2/2+x^3/3-x^4/4; logdel:=convert(asympt(collect(mj*subs(x=DE,logx),m2),m2,asy),polynom); if withprint=1 then blogdel:=btransfs(logdel); print(blogdel); fi; RETURN(eval(logdel)) end; getlogDelta(1); evalExp:=proc(x,asy,del) ## plug in s=u, cut off negligible terms and ## do sorting of terms for evaluating the derivatives collect(collect(collect(collect(convert(asympt(collect(subs(exp(0)=1,subs(s=u,x)),m2),m2,asy),polynom)/del,f0),k),u),m2); end; ## calculating H(u) (A.53) getHu:=proc(withprint,pot,asy) local logH,p,logH1,logH21,H22,H21,H2,H2d1,H2d2,H2d3,H2d4,delf,HH0,HH1,HH2,HH3,HH4,HH,Hu,Pu,PP1, erg, cnjkn; # get log(1-Delta(s,u)) logH:=getlogDelta(withprint,asy-2); # factor of H1... delf:=2*f0*m2; ### logH1:=delf*(s-u); logH21:=logH-logH1; #c.f. (A.54) #expansion for the density # print(logH21); p:=f0+f1*t+f2*t^2/2; H22:=transfs((((u+s)/2)^pot)*subs(t=s,p)); H21:=exp(logH21); H2:=H21*H22; # ### derviatves of H2 H2d1:=diff(H2,s); H2d2:=diff(H2d1,s); H2d3:=diff(H2d2,s); H2d4:=diff(H2d3,s); HH0:=evalExp(H2,asy,delf); HH1:=evalExp(H2d1,asy,delf^2); HH2:=evalExp(H2d2,asy,delf^3); HH3:=evalExp(H2d3,asy,delf^4); HH4:=evalExp(H2d4,asy,delf^5); HH:=collect(HH0-HH1+HH2-HH3+HH4,m2); ## (c.f. (A.55)) Pu:=1/2+f0*u+f1*u^2/2+f2*u^3/6; #expansion for the 1/cdf PPa1:=transf((m-j)*convert(taylor(1/Pu,u,4),polynom))/m2;#print(PP1); Hu:=collect(collect(collect(collect(convert(asympt(collect(HH*PPa1,m2),m2,asy),polynom),j1),k),u),m2); if withprint=1 then print(Hu); print(btransfs(Hu)); fi; if withprint=3 then cnjkn:=collect(collect(subs(j1=-k/2,asbinom(2*m-k-1,m-j,3,2,0)/4/2^(1/2)/m2^3),k),m2); print(collect(collect(collect(collect(convert(asympt(collect(subs(j1=-k/2,Hu)*cnjkn,m2),m2,6),polynom),k),f0),u),m2)); fi; erg:=collect(subs(u=t*m2,Hu),m2); # RETURN(eval(erg)) end; Hux:=getHu(3,2,5); ### corresponding function for getasintegrand # getasintegrandmid:=proc(i,pot,withprint,withprintall, withHu) local inte,inte0, xnjkn,Fxnjkn,cnjkn,aftn,atintn,hyn,binte,dfLpn,dfLmn,dfexpN,hyn0,Hu; # cnjkn:=asbinom(2*m-k-1,m-j,i-1,pot,withprintall); Fxnjkn:=convert(asfrac(m-j,2*m-k-1,i,1,withprintall),polynom); xnjkn:=asquantile(Fxnjkn,i,0,withprintall); # dfLpn:=dfLogf(2*m-k-1,m-j,i+6,withprintall); dfLmn:=dfLogf(-2*m+k+1,m-k-1+j,i+6,withprintall); # dfexpN:=collect(convert(asympt(collect(collect(collect(dfLpn+dfLmn,k),j1),m2),m2,2),polynom),DF); # hyn0:=asdf(xnjkn,Fxnjkn,dfexpN,i+1,withprintall); hyn:=convert(asympt(hyn0,m2,i),polynom); aftn:=asft(xnjkn,i+1,withprintall); # # new: replace t^2 by H(t).... [need higher accuracy than in getasintegrand!] Hu:=getHu(withprint,pot,i+2); if withHu=1 then print(Hu); fi; # as in getasintegrand but with Hu instead of t^2 atintn:=astint(xnjkn,Hu,i+2,withprintall); #[need higher accuracy than in getasintegrand!] # if withprint=2 then print(hyn); print(aftn); print(atintn); fi; # rest remains the same... inte0:=cnjkn*hyn*aftn*atintn*(2*f0)^pot; inte:=collect(collect(collect(collect(convert(asympt(collect(inte0,m2),m2,i),polynom),f0),k),j1),m2); if withprint=1 then binte:=btransf(inte); print(binte); fi; RETURN(eval(inte)) end; getasriskmid:=proc(i,pot,pm,withprint,withprintall) local intjkn,intk0,intk1,erg; intjkn:=convert(asympt(getasintegrandmid(i,pot,withprintall,withprintall,withprintall),m2,i),polynom); # decide how to choose j: if withprint=1 then print(asympt(intjkn,m2,1)); fi; if pm=1 then intk0:=collect(collect(subs(j1=k/2,intjkn),y),m2); else intk0:=collect(collect(subs(j1=-k/2,intjkn),y),m2); fi; intk1:=inteyout(intk0,withprintall); erg:=intekout(2*m2^2,n/2,intk1,3,0); RETURN(eval(erg)); end; getasriskmid(3,2,1,1,0); getallrisksmid:=proc(i,pm,withprintall,withr0,withf00) local asMSE,asBias,asBias2,asVar; asMSE:=getasriskmid(i,2,pm,0,withprintall); asBias:=getasriskmid(i+1,1,pm,0,withprintall); asBias2:=collect(collect(collect(convert(subs(n2=n^(1/2),asympt(subs(n=n2^2,asBias^2),n2,3)),polynom),r),f0),n); asVar:=collect(collect(collect(asMSE-asBias2,r),f0),n); if withf00=1 then if withr0=1 then print(subs(f1=0,r=0,asMSE));print(subs(f1=0,r=0,asBias));print(subs(f1=0,r=0,asBias2));print(subs(f1=0,r=0,asVar)); else print(subs(f1=0,asMSE));print(subs(f1=0,asBias));print(subs(f1=0,asBias2));print(subs(f1=0,asVar)); fi; else if withr0=1 then print(subs(r=0,asMSE));print(subs(r=0,asBias));print(subs(r=0,asBias2));print(subs(r=0,asVar)); else print(asMSE);print(asBias);print(asBias2);print(asVar); fi; fi; end; getallrisksmid(3,1,0,1,1);