Simulation/simulation5_1/MCMC1/spline_basis.R

In shijiaw/smcDE: Adaptive semiparametric Bayesian differential equations via sequential Monte Carlo methods

tobs = times
#knots = tobs
#knots = c(1,45,90)
knots = seq(times[1], max(times),by = 4)
#knots = c(times[1], seq(2, max(times),by=2))
#knots = c(1,seq(5,90,by=5))
nknots   = length(knots)
norder   = 4
nbasis   = length(knots) + norder - 2
bsbasis = create.bspline.basis(range(knots),nbasis,norder,knots)


# basis values at sampling points
basismat   = eval.basis(tobs, bsbasis)

# square of basis matrix (Phi). It is used frequent for optimization,
# so we calculate in advance
basismat2  = t(basismat)%*%basismat

# values of the first derivative of basis functions at sampling points
Dbasismat  = eval.basis(tobs, bsbasis, 1)

# values of the second derivative of basis functions at sampling points
D2basismat = eval.basis(tobs, bsbasis, 2)

# number of quadrature points (including two knots) between two knots
nquad                        = 5;

# set up quadrature points and weights
#[sinbasis, quadpts, quadwts] 
quadre = quadset(nquad, bsbasis)
quadpts=quadre$quadvals[,1]
quadwts=quadre$quadvals[,2]

# values of the second derivative of basis functions at quadrature points
D0quadbasismat = eval.basis(quadpts, bsbasis, 0)
D1quadbasismat = eval.basis(quadpts, bsbasis, 1)
#D2quadbasismat = eval.basis(quadpts, bsbasis, 2)

#Rmat  = t(D2quadbasismat)%*%(D2quadbasismat*((quadwts)%*%matrix(1,1,nbasis))) 


psi <- basismat  # bsplineS(norder=4,nderiv=0)
psi1 <- Dbasismat # bsplineS(norder=4,nderiv=1)
#psi2 <- D2basismat # bsplineS(norder=4,nderiv=2)