R/MOL.density.R

In eta21/DiffusionRimp: Inference and Analysis for Diffusion Processes via Data Imputation and Method of Lines

MOL.density=function(Xs,s=0,t=2,xlims,N=31,delt=1/100,mu,sig,final.only=FALSE)
{
  xx1=seq(xlims[1],xlims[2],length=N)
  dx1=diff(xx1)[1]
  dxx1=dx1^2
   if((missing(mu)&&missing(sig)))
  {

    namess=c('mu','sig')
    namess2=c('muu','sigg')
    txt =rep('+rep(0,length(X))',2)
    func.list=rep(0,length(namess))
    obs=objects(pos=1)
    muu  =function(X,t){}
    sigg =function(X,t){}
    for(i in 1:length(namess))
    {
      if(sum(obs==namess[i]))
      {
        txt[i] = paste0(body(namess[i])[2],txt[i])
      }
    }
    body(muu)  =  parse(text=txt[1])
    body(sigg) =  parse(text=txt[2])
  }else
  {
    muu  =function(X,t){}
    sigg =function(X,t){}
    body(muu) =  parse(text=paste0(mu,'+rep(0,length(X))'))
    body(sigg)=  parse(text=paste0(sig,'+rep(0,length(X))'))
  }

  
  f=function(U,tme)
  {
    MU1= muu(xx1,tme)*U
    SU1= sigg(xx1,tme)^2*U
    D1 = (MU1[-c(1,2)]-MU1[-c(N-1,N)])/dx1/2
    D2 =  1/2*((SU1[-c(1,2)]-2*SU1[-c(1,N)]+SU1[-c(N-1,N)]))/dxx1
    MMM1=-D1+D2
    return(c(0,MMM1,0))
  }
  
  M=rep(0,N)
  wwwhhh1=min(which(abs(xx1-Xs[1])==min(abs(xx1-Xs[1]))))
  xx1=xx1+(Xs-xx1[wwwhhh1])
  M[wwwhhh1]=1
  N.mesh=round((t-s)/delt)+1
  if(!final.only)
  {
    MM = array(0,dim=c(N,N.mesh))
  }
  ttt=s
  
  t.alpha=
    c(0.000000000000000000000000000000000000000000000000000000000000,
      0.100000000000000000000000000000000000000000000000000000000000,
      0.539357840802981787532485197881302436857273449701009015505500,
      0.809036761204472681298727796821953655285910174551513523258250,
      0.309036761204472681298727796821953655285910174551513523258250,
      0.981074190219795268254879548310562080489056746118724882027805,
      0.833333333333333333333333333333333333333333333333333333333333,
      0.354017365856802376329264185948796742115824053807373968324184,
      0.882527661964732346425501486979669075182867844268052119663791,
      0.642615758240322548157075497020439535959501736363212695909875,
      0.357384241759677451842924502979560464040498263636787304090125,
      0.117472338035267653574498513020330924817132155731947880336209,
      0.833333333333333333333333333333333333333333333333333333333333,
      0.309036761204472681298727796821953655285910174551513523258250,
      0.539357840802981787532485197881302436857273449701009015505500,
      0.100000000000000000000000000000000000000000000000000000000000,
      1.00000000000000000000000000000000000000000000000000000000000)
  
  
  for(i in 2:N.mesh)
  {
    
    x0=M
    fx0=f(x0,ttt)
    x1=x0+delt*(0.1*fx0)
    fx1=f(x1,ttt+t.alpha[2]*delt)
    x2=x0+delt*(-0.915176561375291*fx0  +1.45453440217827*fx1)
    fx2=f(x2,ttt+t.alpha[3]*delt)
    x3=x0+delt*( 0.202259190301118*fx0  +0.606777570903354*fx2)
    fx3=f(x3,ttt+t.alpha[4]*delt)
    x4=x0+delt*( 0.184024714708644*fx0  +0.197966831227192*fx2-0.0729547847313633*fx3)
    fx4=f(x4,ttt+t.alpha[5]*delt)
    x5=x0+delt*( 0.0879007340206681*fx0 +0.410459702520261*fx3+0.482713753678866*fx4)
    fx5=f(x5,ttt+t.alpha[6]*delt)
    x6=x0+delt*(0.085970050490246*fx0   +0.330885963040722*fx3+0.48966295730945*fx4-0.0731856375070851*fx5)
    fx6=f(x6,ttt+t.alpha[7]*delt)
    x7=x0+delt*(0.120930449125334*fx0   +0.260124675758296*fx4+0.0325402621549091*fx5-0.0595780211817361*fx6)
    fx7=f(x7,ttt+t.alpha[8]*delt)
    x8=x0+delt*(0.110854379580391*fx0   -0.0605761488255006*fx5+0.321763705601778*fx6+0.510485725608063*fx7)
    fx8=f(x8,ttt+t.alpha[9]*delt)
    x9=x0+delt*(0.112054414752879*fx0   -0.144942775902866*fx5-0.333269719096257*fx6+0.49926922955688*fx7+0.509504608929686*fx8)
    fx9=f(x9,ttt+t.alpha[10]*delt)
    x10=x0+delt*(0.113976783964186*fx0  -0.0768813364203357*fx5+0.239527360324391*fx6+0.397774662368095*fx7+0.0107558956873607*fx8-0.327769124164019*fx9)
    fx10=f(x10,ttt+t.alpha[11]*delt)
    x11=x0+delt*(0.0798314528280196*fx0 -0.0520329686800603*fx5-0.0576954146168549*fx6+0.194781915712104*fx7+0.145384923188325*fx8-0.0782942710351671*fx9-0.114503299361099*fx10)
    fx11=f(x11,ttt+t.alpha[12]*delt)
    x12=x0+delt*(0.985115610164857*fx0  +0.330885963040722*fx3+0.48966295730945*fx4-1.37896486574844*fx5-0.861164195027636*fx6+5.78428813637537*fx7+3.28807761985104*fx8-2.38633905093136*fx9-3.25479342483644*fx10-2.16343541686423*fx11)
    fx12=f(x12,ttt+t.alpha[13]*delt)
    x13=x0+delt*(0.895080295771633*fx0  +0.197966831227192*fx2-0.0729547847313633*fx3-0.851236239662008*fx5+0.398320112318533*fx6+3.63937263181036*fx7+1.5482287703983*fx8-2.12221714704054*fx9-1.58350398545326*fx10-1.71561608285936*fx11-0.0244036405750127*fx12)
    fx13=f(x13,ttt+t.alpha[14]*delt)
    x14=x0+delt*(-0.915176561375291*fx0+1.45453440217827*fx1+0*fx2+0*fx3-0.777333643644968*fx4+0*fx5-0.0910895662155176*fx6+0.0910895662155176*fx12+0.777333643644968*fx13)
    fx14=f(x14,ttt+t.alpha[15]*delt)
    x15=x0+delt*(0.1*fx0-0.157178665799771*fx2+0.157178665799771*fx14)
    fx15=f(x15,ttt+t.alpha[16]*delt)
    x16=x0+delt*(0.181781300700095*fx0+0.675*fx1+0.34275815984719*fx2+0*fx3+0.259111214548323*fx4-0.358278966717952*fx5-1.04594895940883*fx6+0.930327845415627*fx7+1.77950959431708*fx8+0.1*fx9-0.282547569539044*fx10-0.159327350119973*fx11-0.145515894647002*fx12-0.259111214548323*fx13-0.34275815984719*fx14-0.675*fx15)
    fx16=f(x16,ttt+t.alpha[17]*delt)
    
    M=M+(0.0333333333333333333333333333333333333333333333333333333333333*fx0
         +0.0250000000000000000000000000000000000000000000000000000000000*fx1
         +0.0333333333333333333333333333333333333333333333333333333333333*fx2
         +0.000000000000000000000000000000000000000000000000000000000000*fx3
         +0.0500000000000000000000000000000000000000000000000000000000000*fx4
         +0.000000000000000000000000000000000000000000000000000000000000*fx5
         +0.0400000000000000000000000000000000000000000000000000000000000*fx6
         +0.000000000000000000000000000000000000000000000000000000000000*fx7
         +0.189237478148923490158306404106012326238162346948625830327194*fx8
         +0.277429188517743176508360262560654340428504319718040836339472*fx9
         +0.277429188517743176508360262560654340428504319718040836339472*fx10
         +0.189237478148923490158306404106012326238162346948625830327194*fx11
         -0.0400000000000000000000000000000000000000000000000000000000000*fx12
         -0.0500000000000000000000000000000000000000000000000000000000000*fx13
         -0.0333333333333333333333333333333333333333333333333333333333333*fx14
         -0.0250000000000000000000000000000000000000000000000000000000000*fx15
         +0.0333333333333333333333333333333333333333333333333333333333333*fx16)*delt
    ttt=ttt+delt
    if(!final.only)
    {
      MM[,i]=M
    }
      if(any(is.na(M))){stop(paste("NAs produced at time",ttt,". Perhaps increase NN and decrease delt. Alternatively, widen xlims."))}
  }
  if(final.only)
  {
    MM=M
  }
  
   ret = (list(density=MM/dx1,Xt=xx1,time=seq(s,t,delt)))
   class(ret) = 'MOL.density'
   return(ret)
   
}