Модель линейной регрессии степенной функции с двумя членами (y = ax ^ b + cx ^ d) в R

Question

Я пытаюсь найти взаимосвязь между устойчивой скоростью ветра и множителями порыва ветра [множитель порыва = (величина порыва ветра) / (постоянная скорость ветра)]. Сначала я использовал модель мощности (y = ax ^ b) и линейную регрессию для их моделирования, но я также хотел бы попытаться подогнать данные к модели мощности с двумя членами (y = ax ^ b + cx ^ d). .

В модели мощности с одним членом я смог взять логарифм каждой переменной, и затем это позволило мне сохранить коэффициенты (значения a и b) для модели мощности. Я не уверен, как изменить свой код, чтобы иметь возможность создать модель мощности с двумя членами.

model = lm(log(data$GustMult)~log(data$SusWindSpeed))  # not sure how to change for a power model with two terms
coefs = unname(coef(model))
coefa = exp(coefs[1])
coefb = coefs[2]
coefc = exp(coefs[3])  #added for new model
coefd = coefs[4]  #added for new model
rsquared = summary(model)$r.squared

Я попытался запустить код, как указано выше, однако он предоставляет только значения для a и коэффициенты b по-прежнему.

Заранее благодарим вас за любые предложения!

Answer 1

Думаю, Дэни прав, что не существует уникального набора параметров, который подходил бы для этой модели. В любом случае, я думаю, вам нужно будет приспособить это как нелинейную модель.

Я понятия не имею, как выглядят ваши данные, но допустим, что это примерно так:

set.seed(69)

SusWindSpeed <- rgamma(1000, 2, .125)
GustMult  <- (-0.98)*SusWindSpeed^(0.12) + (0.1)*SusWindSpeed^(-0.32) + rnorm(1000, 3, 0.025)
data <- data.frame(SusWindSpeed, GustMult)

plot(data$SusWindSpeed, data$GustMult)

Then you could try fitting a non-linear model like this:

modA <- nls(GustMult ~ var_a * SusWindSpeed^var_b + var_c * SusWindSpeed^var_d + Const,
            start = list(var_a = -1, var_b = 0.1, var_c = 0.1, var_d = -0.3, Const = 3),
            data = data,
            control = list(maxiter = 500))

summary(modA)
#> 
#> Formula: GustMult ~ var_a * SusWindSpeed^var_b + var_c * SusWindSpeed^var_d + 
#>     Const
#> 
#> Parameters:
#>       Estimate Std. Error t value Pr(>|t|)    
#> var_a -1.28314    0.72635  -1.767   0.0776 .  
#> var_b  0.10381    0.04387   2.367   0.0181 *  
#> var_c  0.01506    0.03884   0.388   0.6983    
#> var_d -1.45958    3.32064  -0.440   0.6604    
#> Const  3.38496    0.76287   4.437 1.01e-05 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.02452 on 995 degrees of freedom
#> 
#> Number of iterations to convergence: 22 
#> Achieved convergence tolerance: 3.379e-07

И мы можем показать, что это разумная подгонка:

test_SusWindSpeed <- seq(0, 60, 0.05)

plot.new()
plot(data$SusWindSpeed, data$GustMult)
points(test_SusWindSpeed, 
       predict(modA, newdata = list(SusWindSpeed = test_SusWindSpeed)), 
       col = "red")

However, you might find you struggle to get this model to converge, or that you get different results with different start parameters.

^{Created on 2020-07-09 by the пакет REPEX (v0.3.0)}