在数据集中拟合对数曲线
Fitting logaritmic curve in a dataset
我有一个数据集 archivo
包含自 2003 年以来政府拍卖的每个持续时间的债券利率。前几行是:
Fecha 1 2 3 4 5 6 7 8 9 10 11 12 18 24
2003-01-02 NA NA NA NA NA 44.9999 NA NA 52.0002 NA NA NA NA NA
2003-01-03 NA NA NA NA NA NA NA NA NA NA NA NA NA NA
2003-01-06 NA NA NA NA NA NA NA NA NA NA NA NA NA NA
2003-01-07 NA NA NA NA NA 40.0000 NA NA 45.9900 NA NA NA NA NA
2003-01-08 NA NA NA NA NA NA NA NA NA NA NA NA NA NA
2003-01-09 NA NA NA NA NA 37.0000 NA NA 41.9999 NA NA NA NA NA
名称为 1 到 24 的每一列对应不同的持续时间。 (1 个月、2 个月、...、24 个月)。并非所有期限都在拍卖日售出。这就是为什么我有 NA
s.
我需要用对数拟合曲线计算 NA
s(缺失)率,每一行至少有超过 1 个值。对于包含所有 NA
的行,我可以使用前面构造的曲线。
我知道我可以 运行 像这样的代码:
x<-colnames(archivo[,-1]) # to keep the durations
y<-t(archivo[1,-1])
estimacion<-lm(y ~ log(x))
param<-estimacion$coefficients
并获取第一行的系数。然后 运行 一个循环,对每一行都这样做。
有什么办法可以直接对整个数据集做,不做循环就得到每一行(每个log fitting)的参数吗?
希望问题够清楚。
提前致谢!
尝试:
dat <- as.data.frame(t(archivo[,-1])) ## transpose you data frame
## a function to fit a model `y ~ log(x)` for response vector `y`
fit_model <- function (y) {
non_NA <- which(!is.na(y)) ## non-NA rows index
if (length(non_NA) > 1) {
## there are at least 2 data points, a linear model is possible
lm.fit(cbind(1, log(non_NA)), y[non_NA])$coef
} else {
## not sufficient number of data, return c(NA, NA)
c(NA, NA)
}
}
## fit linear model column-by-column
result <- sapply(dat, FUN = fit_model)
请注意,我正在使用 lm.fit()
,由 lm()
调用的内核拟合例程。如果您不熟悉 ?lm.fit
,请阅读它。它需要 2 个基本参数:
- 首先是模型矩阵。您的模型
y ~ log(x)
的模型矩阵是 matrix(c(rep(1,24), log(1:24)), ncol = 2)
。您也可以通过 model.matrix(~log(x), data = data.frame(x = 1:24))
. 构造它
- 第二个是响应向量。对于您的问题,它是
dat
. 的一列
不像lm()
可以处理NA
,lm.fit()
不能。所以我们需要自己从模型矩阵和响应向量中删除 NA
行。 non_NA
变量就是这样做的。注意,你的模型y ~ log(x)
涉及2个参数/系数,所以拟合至少需要2个数据。如果没有足够的数据,模型拟合是不可能的,我们 return c(NA, NA)
.
最后,我使用 sapply()
逐列拟合线性模型,只保留 $coef
的系数。
测试
我正在使用您在问题中发布的示例行。使用上面的代码,我得到:
# V1 V2 V3 V4 V5 V6
# x1 14.06542 NA NA 13.53005 NA 14.90533
# x2 17.26486 NA NA 14.77316 NA 12.33127
每列给出 dat
的每一列(或 archivo
的每一行)的系数。
更新
最初我在 lm.fit()
中使用 matrix(rep(1,24), log(1:24))[non_NA, ]
作为模型矩阵。但这效率不高。它首先生成完整的模型矩阵,然后删除带有 NA
的行。双重想法表明这样更好:cbind(1, log(non_NA))
.
我有一个数据集 archivo
包含自 2003 年以来政府拍卖的每个持续时间的债券利率。前几行是:
Fecha 1 2 3 4 5 6 7 8 9 10 11 12 18 24
2003-01-02 NA NA NA NA NA 44.9999 NA NA 52.0002 NA NA NA NA NA
2003-01-03 NA NA NA NA NA NA NA NA NA NA NA NA NA NA
2003-01-06 NA NA NA NA NA NA NA NA NA NA NA NA NA NA
2003-01-07 NA NA NA NA NA 40.0000 NA NA 45.9900 NA NA NA NA NA
2003-01-08 NA NA NA NA NA NA NA NA NA NA NA NA NA NA
2003-01-09 NA NA NA NA NA 37.0000 NA NA 41.9999 NA NA NA NA NA
名称为 1 到 24 的每一列对应不同的持续时间。 (1 个月、2 个月、...、24 个月)。并非所有期限都在拍卖日售出。这就是为什么我有 NA
s.
我需要用对数拟合曲线计算 NA
s(缺失)率,每一行至少有超过 1 个值。对于包含所有 NA
的行,我可以使用前面构造的曲线。
我知道我可以 运行 像这样的代码:
x<-colnames(archivo[,-1]) # to keep the durations
y<-t(archivo[1,-1])
estimacion<-lm(y ~ log(x))
param<-estimacion$coefficients
并获取第一行的系数。然后 运行 一个循环,对每一行都这样做。
有什么办法可以直接对整个数据集做,不做循环就得到每一行(每个log fitting)的参数吗?
希望问题够清楚。
提前致谢!
尝试:
dat <- as.data.frame(t(archivo[,-1])) ## transpose you data frame
## a function to fit a model `y ~ log(x)` for response vector `y`
fit_model <- function (y) {
non_NA <- which(!is.na(y)) ## non-NA rows index
if (length(non_NA) > 1) {
## there are at least 2 data points, a linear model is possible
lm.fit(cbind(1, log(non_NA)), y[non_NA])$coef
} else {
## not sufficient number of data, return c(NA, NA)
c(NA, NA)
}
}
## fit linear model column-by-column
result <- sapply(dat, FUN = fit_model)
请注意,我正在使用 lm.fit()
,由 lm()
调用的内核拟合例程。如果您不熟悉 ?lm.fit
,请阅读它。它需要 2 个基本参数:
- 首先是模型矩阵。您的模型
y ~ log(x)
的模型矩阵是matrix(c(rep(1,24), log(1:24)), ncol = 2)
。您也可以通过model.matrix(~log(x), data = data.frame(x = 1:24))
. 构造它
- 第二个是响应向量。对于您的问题,它是
dat
. 的一列
不像lm()
可以处理NA
,lm.fit()
不能。所以我们需要自己从模型矩阵和响应向量中删除 NA
行。 non_NA
变量就是这样做的。注意,你的模型y ~ log(x)
涉及2个参数/系数,所以拟合至少需要2个数据。如果没有足够的数据,模型拟合是不可能的,我们 return c(NA, NA)
.
最后,我使用 sapply()
逐列拟合线性模型,只保留 $coef
的系数。
测试
我正在使用您在问题中发布的示例行。使用上面的代码,我得到:
# V1 V2 V3 V4 V5 V6
# x1 14.06542 NA NA 13.53005 NA 14.90533
# x2 17.26486 NA NA 14.77316 NA 12.33127
每列给出 dat
的每一列(或 archivo
的每一行)的系数。
更新
最初我在 lm.fit()
中使用 matrix(rep(1,24), log(1:24))[non_NA, ]
作为模型矩阵。但这效率不高。它首先生成完整的模型矩阵,然后删除带有 NA
的行。双重想法表明这样更好:cbind(1, log(non_NA))
.