Как реализовать модель R в коде C ++

Question

У меня есть модель в R:

> s1 <- toys[1:10000,]
> model <- glm(V11~V2+V3+V5+V7+V8+V9+V10,gaussian,s1)
> model

Call:  glm(formula = V11 ~ V2 + V3 + V5 + V7 + V8 + V9 + V10, family = gaussian, 
    data = s1)

Coefficients:
(Intercept)           V2           V3           V5           V7           V8           V9          V10  
  -0.900106     0.006385    -0.005080     1.006324     0.229282     0.012391    -0.049307    -0.186450  

Degrees of Freedom: 9999 Total (i.e. Null);  9992 Residual
Null Deviance:      11050000 
Residual Deviance: 121200   AIC: 53340 

Теперь, как мне запрограммировать эту модель R как функцию C?(RTFM со ссылкой было бы достаточно)

Возможно, мне просто нужно умножить все коэффициенты из модели R на их соответствующие входные данные и добавить все условия, чтобы получить окончательный результат?

float model(float v2, float v3, ... float v10)
{
  return -0.900106 * v2 + 0.006385 * v3 + .. + (-0.186450) * v10;
}

Мне нужен автономный код, не зависящий от каких-либо внешних источников

Answer 1

Вы запрашиваете модель линейной регрессии (здесь R glm() обозначает обобщенную линейную модель, но поскольку вы используете идентификационную ссылку, вы в конечном итоге получаете линейную регрессию). В C доступно несколько реализаций, например, библиотека apophenia , которая имеет хороший набор статистических функций с привязками для MySQL и Python. Библиотеки GSL и ALGLIB также имеют специальные алгоритмы.

Однако для легкого и почти автономного кода на C я бы посоветовал взглянуть на glm_test.c, доступную в исходном пакете snpMatrix BioC.

---

После обновления вопроса вы, скорее всего, хотите предсказать результат на основе набора параметров регрессии. Тогда, учитывая, что общий вид предполагаемой модели имеет вид y = b0 + b1 * x1 + b2 * x2 + ... + bp * xp, где b0 - точка пересечения, а b1, ..., bp - коэффициенты регрессии ( по оценкам данных), вычисление довольно простое, так как оно составляет взвешенную сумму: возьмите каждое наблюдаемое значение для ваших p-предикторов и умножьте на b (не забывайте термин перехват!).

Вы можете дважды проверить свои результаты с помощью функции R predict(); Вот пример с двумя предикторами, названными V1 и V2, 100 наблюдениями и регулярной сеткой новых значений для прогнозирования результата (вы также можете использовать свои собственные данные):

> df <- transform(X <- as.data.frame(replicate(2, rnorm(100))), 
                                     y = V1+V2+rnorm(100))
> res.lm <- lm(y ~ ., df)
> new.data <- data.frame(V1=seq(-3, 3, by=.5), V2=seq(-3, 3, by=.5))
> coef(res.lm)
(Intercept)          V1          V2 
0.006712008 0.980712578 1.127586352 
> new.data
     V1   V2
1  -3.0 -3.0
2  -2.5 -2.5
...
> 0.0067 + 0.9807*-3 + 1.1276*-3  # with approximation
[1] -6.3182
> predict(res.lm, new.data)[1]
        1 
-6.318185