Вы можете использовать expand.grid()
, чтобы получить все комбинации k
и l
; затем создайте каждый член в виде строки; используйте paste()
с аргументом collapse
, чтобы объединить термины в одну строку; и привести в формулу:
dunham_formula <- function(k, l) {
terms <- with(expand.grid(k = k, l = l), {
glue::glue("y{k}{l} * (v + .5)^{k} * (J * (J + 1))^{l}")
})
as.formula(paste("E ~", paste0(terms, collapse = " + ")))
}
dunham_formula(0:1, 0:1)
#> E ~ y00 * (v + 0.5)^0 * (J * (J + 1))^0 + y10 * (v + 0.5)^1 *
#> (J * (J + 1))^0 + y01 * (v + 0.5)^0 * (J * (J + 1))^1 + y11 *
#> (v + 0.5)^1 * (J * (J + 1))^1
#> <environment: 0x00000000159e58f8>
Давайте проверим это с некоторыми поддельными данными:
set.seed(42)
n <- 50
df <- data.frame(
E = rexp(n),
v = runif(n),
J = runif(n)
)
summary(nls(dunham_formula(0:1, 0:1), data = df))
#> Warning in nls(dunham_formula(0:1, 0:1), data = df): No starting values specified for some parameters.
#> Initializing 'y00', 'y10', 'y01', 'y11' to '1.'.
#> Consider specifying 'start' or using a selfStart model
#>
#> Formula: E ~ y00 * (v + 0.5)^0 * (J * (J + 1))^0 + y10 * (v + 0.5)^1 *
#> (J * (J + 1))^0 + y01 * (v + 0.5)^0 * (J * (J + 1))^1 + y11 *
#> (v + 0.5)^1 * (J * (J + 1))^1
#>
#> Parameters:
#> Estimate Std. Error t value Pr(>|t|)
#> y00 0.11636 1.60105 0.073 0.942
#> y10 0.68614 1.56022 0.440 0.662
#> y01 0.31082 1.61183 0.193 0.848
#> y11 0.05498 1.51326 0.036 0.971
#>
#> Residual standard error: 1.406 on 46 degrees of freedom
#>
#> Number of iterations to convergence: 1
#> Achieved convergence tolerance: 1.81e-07
Поскольку это на самом деле линейная модель, вы можете вместо этого создать функцию, которая возвращает базовую матрицу, а затем использовать lm()
для подгонки модели:
dunham_basis <- function(v, J, k, l) {
dunham_term <- function(k, l) {
(v + .5) ^ k * (J * (J + 1)) ^ l
}
indices <- expand.grid(k = k, l = l)
cols <- with(indices, Map(dunham_term, k, l))
names(cols) <- apply(indices, 1, paste, collapse = ",")
do.call("cbind", cols)
}
df$Y <- with(df, dunham_basis(v, J, k = 0:1, l = 0:1))
summary(lm(E ~ 0 + Y, data = df))
#>
#> Call:
#> lm(formula = E ~ 0 + Y, data = df)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -1.5828 -0.7332 -0.3373 0.2736 5.3756
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> Y0,0 0.11636 1.60105 0.073 0.942
#> Y1,0 0.68614 1.56022 0.440 0.662
#> Y0,1 0.31082 1.61183 0.193 0.848
#> Y1,1 0.05498 1.51326 0.036 0.971
#>
#> Residual standard error: 1.406 on 46 degrees of freedom
#> Multiple R-squared: 0.432, Adjusted R-squared: 0.3826
#> F-statistic: 8.745 on 4 and 46 DF, p-value: 2.453e-05
Создано в 2019-07-02 пакетом Представление (v0.3.0)