你认为这两个代码的相同点和不同点有哪些？第一个代码：install.packages(“fda”)
library(fda)

tobs = seq(0,1,0.01)
nobs = length(tobs)
knots = c(seq(0,1,0.1));
nknots = length(knots);
norder = 4;
nbasis = length(knots) + norder - 2;
basis = create.bspline.basis(c(min(tobs),max(tobs)),nbasis,norder,knots);


# basis values at samplincurv points

basismat = eval.basis(tobs, basis);
dim(basismat)
quartz()
plot(tobs,basismat[,1],type = “l”,col=1,lwd=3)
lines(tobs,basismat[,2],type = “l”,col=2,lwd=3)

quartz()
matplot(tobs,basismat,type=‘l’,lwd=2,lty=1, xlab=‘day’,ylab=‘basis’,cex.lab=1.5,cex.axis=1.5)
for (i in 1:nknots)
{
abline(v=knots[i],type=“l”, lty=2, lwd=3)
}

# Comments: If x(0) = 1, set the coefficient to the first basis function = 1;

# evaluate the first derivative of the basis functions
Dbasismat = eval.basis(tobs, basis,1);
quartz()
matplot(tobs,Dbasismat,type=‘l’,lwd=2,lty=1, xlab=‘day’,ylab=‘basis’,cex.lab=1.5,cex.axis=1.5)
for (i in 1:nknots)
{
abline(v=knots[i],type=“l”, lty=2, lwd=3)
}

# true curve
ytru = (tobs-0.3)^2
plot(tobs,ytru,type = “l”)

# put noise to the true curve and generate noisy data
nobs = length(tobs)
noise = 0.03rnorm(nobs)
yobs = ytru + noise
points(tobs,yobs)

# estimate basis coefficient
Mmat = ginv(t(basismat)%%basismat)%%t(basismat)
chat = Mmat%%yobs

# fitted curve
yhat = basismat%%chat;
lines(tobs,yhat,type = “l”,col=“red”)

# estimate the variance of noise
SSE = t(yhat-yobs)%%(yhat-yobs)
sigma2 = SSE/(nobs-nbasis)
sigma2
sqrt(sigma2)
# estimate the variance of the fitted curve
Smat = basismat%%Mmat
varYhat = diag(Smat%%Smatmatrix(sigma2,nobs,nobs))

# 95% confidence interval

yhat025 = yhat-1.96sqrt(varYhat)
yhat975 = yhat+1.96sqrt(varYhat)

lines(tobs,yhat025,type=“l”, lty=2, lwd=3,col=“blue”)
lines(tobs,yhat975,type=“l”, lty=2, lwd=3,col=“blue”)

###########################
# Smoothing Splines
#########################


# Use quadrature to get integral - Composite Simpson’s Rule

delta <- 0.02
quadpts <- seq(0,1,delta)
nquadpts <- length(quadpts)
quadwts <- as.vector(c(1,rep(c(4,2),(nquadpts-2)/2),4,1),mode=“any”)
quadwts <- c(1,rep(c(4,2),(nquadpts-1)/2))
quadwts[nquadpts] <- 1
quadwts <- quadwtsdelta/3


# Second derivative of basis functions at quadrature points

Q2basismat = eval.basis(quadpts, basis,2);

# estimates for basis coefficients
Rmat = t(Q2basismat)%%(Q2basismat(quadwts%%t(rep(1,nbasis))))

dim(Rmat)
basismat2 = t(basismat)%%basismat;
lambda = 0.05 # smoothing parameter
Bmat = basismat2 + lambdaRmat;
chat = ginv(Bmat)%%t(basismat)%%yobs;

# fitted value
yhat = basismat%%chat;
yhat2 = basismat%%ginv(t(basismat)%%basismat)%%t(basismat)%%yobs;

quartz()
plot(tobs,ytru,type = “l”)
points(tobs,yobs)
lines(tobs,yhat,type = “l”,col=“red”)
lines(tobs,yhat2,type = “l”,col=“blue”)

# degrees of freedom
Mmat = ginv(Bmat)%%t(basismat)
Smat = basismat%%Mmat
df = sum(diag(Smat))
c(df,nbasis)
# estimate the variance of noise
SSE = t(yhat-yobs)%%(yhat-yobs)
sigma2 = SSE/(nobs-df)
sigma2
sqrt(sigma2)
# estimate the variance of the fitted curve
varYhat = diag(Smat%%Smatmatrix(sigma2,nobs,nobs))

# 95% confidence interval

yhat025 = yhat-1.96sqrt(varYhat)
yhat975 = yhat+1.96*sqrt(varYhat)

lines(tobs,yhat025,type=“l”, lty=2, lwd=3,col=“blue”)
lines(tobs,yhat975,type=“l”, lty=2, lwd=3,col=“blue”)第二个代码：#library(readxl)
#library(fda)
#城镇总生育率,15、16......49岁数据
city_total_birthdata_original <- read_excel("C:/Users/86157/Desktop/人口研究/修复数据 ab/总生育率 ab 不分组.xlsx",range = "B1:AE36") 
#数据处理
birthdata_original <- cbind(city_total_birthdata_original)
birthdata_original <- as.matrix(birthdata_original)
total_quantity <- 240 #240条数据:30*4*2
age_group <- seq(1,35) #35个年龄段
# 15-49岁共35个年龄段
Agegroup <- c("15","16","17","18","19","20","21","22","23","24","25","26","27","28","29","30","31","32","33","34","35","36","37","38","39","40","41","42","43","44","45","46","47","48","49")

# 拟合生育率 ------------------------------------------------------------------

nod=4
Lfdobj <- int2Lfd(max(0, nod-2))
#权重向量**********************88
wtvec=c(rep(1, 5),rep(1,15),rep(1,15))
# 创建一个长度为 17 的向量 K，用于存储基函数数量的取值。
K = seq(4,29,1)
num = length(K)
#创建一个长度为 num 的向量 SSEsave，用于存储每个基函数数量对应的平方误差（SSE）。
#初始化时，向量中的所有值都是 NA
SSEsave = rep(NA,num)
for (i in 1:num) {
  basis <- create.bspline.basis(range(age_group), nbasis = K[i], norder = nod)
  fdParobj = fdPar(basis, Lfdobj, lambda=0)
  smoothlist = smooth.basis(seq(1,35),birthdata_original,fdParobj,wtvec)
  SSEsave[i] = smoothlist$SSE
}
ak <- min(SSEsave)
bk <- which(SSEsave==ak)
cat('K','=',K[bk],'\n')# 打印最佳基函数数量

plot(0, 0, type = "n", xlim = c(1, length(SSEsave)), ylim = range(SSEsave),xaxt="n", xlab = "基函数数量K", ylab = "SSE")
# 绘制数据点
points(1:length(SSEsave), SSEsave)
# 连接数据点
lines(1:length(SSEsave),SSEsave)
# 设置横坐标刻度和标签
axis(1, at = seq(1,32,1),labels = c(4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35))

#dev.off()  # 如果边界过大，使用关闭当前图形设备
#plot(...)  # 重新运行绘图命令

basis <- create.bspline.basis(c(1,35), nbasis =29,norder = 4)

loglam = seq(-4,2,0.25)
nlam = length(loglam)
dfsave = rep(NA,nlam)
gcvsave = rep(NA,nlam)
for (ilam in 1:nlam) {
  #cat(paste('log10 lambda =',loglam[ilam],'\n'))
  lambda = 10^loglam[ilam]
  fdParobj = fdPar(basis, Lfdobj, lambda)
  smoothlist = smooth.basis(seq(1,35),birthdata_original,fdParobj,wtvec)
  dfsave[ilam] = smoothlist$df
  gcvsave[ilam] = sum(smoothlist$gcv)
}
a <- min(gcvsave)
b <- which(gcvsave==a)
lambda <-  10^loglam[b]
cat('λ','=','10^',loglam[b])


par(mfrow=c(1,1),cex=1)
# 创建一个空白画布
plot(0, 0, type = "n", xlim = c(1, length(gcvsave)), 
     ylim = range(gcvsave),xaxt="n", xlab = "lg(λ)", ylab = "GCV Criterion")
# 绘制数据点
points(1:length(gcvsave), gcvsave)
# 连接数据点
lines(1:length(gcvsave),gcvsave)
# 设置横坐标刻度和标签
axis(1, at = seq(1,26,2),
     labels = c(-4.00 , -3.5 , -3.00,  -2.50 , -2.00, -1.50,  -1.00,
                -0.50 ,0.00, 0.50,1.00,1.5,2))

##===========最终的fdPar=============##
fdParobj = fdPar(basis, Lfdobj, lambda)
##===================================##

#平滑曲线
smoothlist = smooth.basis(seq(1,35),birthdata_original,fdParobj,wtvec)

birthfd <- smoothlist$fd

smoothed_value <- matrix(data=NA, nrow = 20, ncol = 124)

for(i in seq(1,)){
  
  smoothed_value[,i] <- eval.fd(seq(1,20,1),birthfd[i])
}

negative_elements <- smoothed_value[smoothed_value < 0]
cat(negative_elements)

#--------------所有的生育率曲线-------------#
#规范刻度
#birth_scale <- c(0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8)
birth_scale <- seq(0,160,10)
# 创建一个新的图形窗口，分为两行两列
#par(mfrow=c(2,2))
#----------------城镇总和生育率--------------------#
plot(birthfd[1:30,],ylim = c(0,160),bty="o",ann=F,xaxt="n",yaxt="n",col="#44c1f0",lty=1)
#bty用来设置边框形式，默认值为"o"，表示四面边框都画出，其余可选值包括"l"（左下）
axis(side=1,at=seq(1,35,1),labels=Agegroup,las=2,cex.axis=1.2,lwd=2,lwd.ticks = 2)
axis(side=2,at=seq(0,160,10),labels=birth_scale,lwd=2,lwd.ticks = 2,cex.axis=1.2,mgp = c(3, 0.8, 0))

title(main="Urban Totalbirth",cex.main=2)
title(xlab="Age groups",line = 4,cex.lab=1.2)
title(ylab="probability",line=2.5,cex.lab=1.2)
-------------------------------------------------------