Nothing
MeanAge=function
(IdotT,
OdotLin,
sol,
times
){
rho=function(a,tk){
startTime=tk-a
startVal=sol(startTime)
if (startTime<0){
print("if")
print(paste("a=",a,"tk=",tk))
return(startVal)
}
times=c(startTime,tk,2*tk)
s=solver(times,OdotLin,startVal)
res=s[2]/sol(tk)
return(res)
}
E=function(tk){
E_integrand=function(Y,a){a*rho(a,tk)}
res=solver(c(0,tk*0.99),E_integrand,0)[2]
return(res)
}
res=sapply(times,E)
return(res)
}
MeanAge2=function
(IdotT,
OdotLin,
sol,
times
){
maxage=max(times)-min(times)
fineTimes=(seq(sqrt(min(times)),sqrt(max(times)),sqrt(maxage)/10000))^2
startval=1
OofT=splinefun(fineTimes,solver(fineTimes,OdotLin,startval))
rho=function(a,tk){
s=OofT(tk)*IdotT(tk-a)/OofT(tk-a)
res=s/sol(tk)
return(res)
}
E=function(tk){
rho_tk=function(a){rho(a,tk)}
E_integrand=function(Y,a){a*rho(a,tk)}
res=solver(c(tk/1000,tk),E_integrand,0)[2]
return(res)
}
res=mcmapply(E,times,mc.cores=16)
return(res)
}
MeanAge3=function
(IdotT,
OdotLin,
sol,
times
){
so=solver(times,OdotLin,1)
deltaT=times[[2]]-times[[1]]
tstart=times[[2]]
startvalue=sol(tstart)
if(startvalue==0){stop("cannot compute the solution since we are in a fixed point")}
Id=function(t){IdotT(t)}
l=100
pdf(file="runit.MeanAge.Densities.pdf",paper="a4r")
OofT=splinefun(times,so)
s=function(a,tk){
OofT(tk)*Id(tk-a)/OofT(tk-a)
}
rho=function(a,tk){
res=s(a,tk)/sol(tk)
return(res)
}
E=function(tk){
l=100
ages=seq(0,tk-times[[2]],length=l)
s_tk=function(a){s(a,tk)}
ss_tk=splinefun(ages,mapply(s_tk,ages))
one=integrate(ss_tk,lower=0,upper=tk-times[[2]])$value/sol(tk)
rho_0=(1-one)/(deltaT)
Integrand=function(a){a*rho(a,tk)}
spIntegrand=splinefun(ages,sapply(ages,Integrand))
res=integrate(spIntegrand,lower=0,upper=max(ages))$value+rho_0*(tk-(deltaT)/2)
return(res)
}
times[[1]]<-times[[2]]
dev.off()
res=mclapply(times,E,mc.cores=16)
return(res)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.