从0开始学习R语言--Day30--函数型分析
在研究离散变量之间的影响时,我们往往只能获取类似中位数,平均数点来额外数据特点;但如果数据本身具有时间特性的话,我们可以尝试运用函数型分析,将静态的离散点转为动态过程来分析,即若本来是分析离散点对另一个变量的影响,那么转换后就变为研究一条曲线的变化趋势是否对应了另一个变量的变化特点。
比如气温数据分析冰淇淋销量时,只看离散点只能知道气温越高,销量越好,但若看趋势,很可能捕捉到在下午气温会比其他时候更高,卖出更多冰淇淋的特点。以下是一个例子:
library(fda)# 生成模拟数据 -------------------------------------------------
set.seed(123)
n <- 100 # 样本量
time_points <- seq(0, 1, length.out=50) # 观测时间点# 生成随机温度曲线(函数型X)
true_curve <- sin(2 * pi * time_points) * exp(-time_points)
X_curves <- t(replicate(n, true_curve + rnorm(length(time_points), sd=0.5)))# 生成产量Y(与曲线变化率相关)
# Y <- 5 * apply(X_curves, 1, function(x) sum(diff(x))) + rnorm(n, sd=2)
Y <- 5 * apply(X_curves, 1, function(x) sum(abs(diff(x)))) + rnorm(n, sd=2)
# 创建函数型数据对象
basis <- create.bspline.basis(c(0,1), nbasis=15)
X_fd <- smooth.basis(time_points, t(X_curves), basis)$fd# 函数型线性回归
fit <- fRegress(Y ~ X_fd) # 注意这里变量名是X_fd# 预测新数据
new_curve <- true_curve + 0.3 * sin(4 * pi * time_points)
new_fd <- smooth.basis(time_points, new_curve, basis)$fdbeta_fd <- fit$betaestlist[[2]]$fd # 第1个是截距,第2个是X_fd的系数# 计算预测值:Y = 截距 + 积分(X(t) * beta(t) dt)
intercept <- fit$betaestlist[[1]]$fd$coefs # 截距项
pred_integral <- inprod(new_fd, beta_fd) # 函数内积
prediction <- intercept + pred_integralprint(paste("预测产量:", prediction))# 输出结果
plot(new_fd, main="新温度曲线")
print(paste("预测产量:", round(prediction, 2)))
输出:
[1] "预测产量: 153.38"
图像显示温度按照正弦曲线的趋势波动,但在实际应用中,在画图时一般不推荐用标准化的时间轴,会忽略真实的现实信息,像这里便无法判断时间段指的是一个小时还是12个小时。