Блог

Я хотел бы подогнать (очень) регрессию высокого порядка к набору данных в R, однако poly() функция имеет ограничение порядка 25.

Для этого приложения мне нужен порядок в диапазоне от 100 до 120.

 model <- lm(noisy.y ~ poly(q,50))
# Error in poly(q, 50) : 'degree' must be less than number of unique points
model <- lm(noisy.y ~ poly(q,30))
# Error in poly(q, 30) : 'degree' must be less than number of unique points
model <- lm(noisy.y ~ poly(q,25))
# OK
  

Многочлены и ортогональные многочлены

poly(x) не имеет жестко заданного предела для degree . Однако на практике существуют два числовых ограничения.

• Базисные функции строятся на уникальном расположении x значений. Полином степени k имеет k 1 базис и коэффициенты. poly генерирует базис без термина перехвата, поэтому degree = k подразумевает k базис и k коэффициенты. Если существуют n уникальные x значения, это должно быть удовлетворено k <= n , иначе просто недостаточно информации для построения полинома. Внутри poly() следующая строка проверяет это условие:

   if (degree >= length(unique(x))) 
      stop("'degree' must be less than number of unique points")
    

• Корреляция между x ^ k и x ^ (k 1) становится все ближе и ближе к 1 по мере k увеличения. Такая скорость приближения, конечно, зависит от x значений. poly сначала генерируется обычный полиномиальный базис, затем выполняется QR-факторизация для нахождения ортогонального диапазона. Если между и возникает дефицит численного ранга x ^ k x ^ (k 1) , poly он также остановится и пожалуется:

   if (QR$rank < degree) 
      stop("'degree' must be less than number of unique points")
    

  Но сообщение об ошибке в этом случае не является информативным. Кроме того, это не обязательно должно быть ошибкой; это может быть предупреждение, которое затем poly может быть сброшено degree rank для продолжения. Может быть, R core может улучшить этот бит??

Ваш метод проб и ошибок показывает, что вы не можете построить полином более 25 степени. Вы можете сначала проверить length(unique(q)) . Если у вас степень меньше этой, но все еще вызывающая ошибку, вы точно знаете, что это связано с числовым недостатком ранга.

Но я хочу сказать, что многочлен более 3-5 степени никогда не бывает полезным!Критической причиной является феномен Рунге. С точки зрения статистической терминологии: полином высокого порядка всегда плохо соответствует данным!. Не наивно думайте, что, поскольку ортогональные многочлены численно более стабильны, чем необработанные многочлены, эффект Рунге можно устранить. Нет, полином степени k образует векторное пространство, поэтому какой бы базис вы ни использовали для представления, они имеют одинаковый диапазон!

Сплайны: кусочно-кубические полиномы и их использование в регрессии

Полиномиальная регрессия действительно полезна, но нам часто нужны кусочные полиномы. Самый популярный выбор — кубический сплайн. Подобно тому, что существуют разные представления для полиномов, существует множество представлений для сплайнов:

• усеченный степенной базис
• естественный кубический сплайн-базис
• B-сплайн-базис

Базис B-сплайна является наиболее численно стабильным, поскольку он имеет компактную поддержку. В результате ковариационная матрица X'X является ленточной, поэтому решение нормальных уравнений (X'X) b = (X'y) очень стабильно.

В R мы можем использовать bs функцию из splines пакета (одного из базовых пакетов R), чтобы получить базис B-сплайна. Ибо bs(x) единственным числовым ограничением степени свободы df является то, что мы не можем иметь больше оснований, чем length(unique(x)) .

Я не уверен, как выглядят ваши данные, но, возможно, вы можете попробовать

 library(splines)
model <- lm(noisy.y ~ bs(q, df = 10))
  

Штрафные сглаживающие / регрессионные сплайны

Регрессионный сплайн по-прежнему может переопределять ваши данные, если вы продолжаете увеличивать степень свободы. Лучшая модель, по-видимому, заключается в выборе наилучшей степени свободы.

Отличным подходом является использование штрафного сглаживающего сплайна или штрафного регрессионного сплайна, так что оценка модели и выбор степени свободы (т. Е. «гладкость») интегрированы.

smooth.spline Функция в stats пакете может выполнять оба. В отличие от того, что, по-видимому, предполагает его название, большую часть времени он просто подгоняет сплайн регрессии, а не сглаживает сплайн. Читайте ?smooth.spline дальше. Для ваших данных вы можете попробовать

 fit <- smooth.spline(q, noisy.y)
  

(Обратите внимание, smooth.spline не имеет интерфейса формулы.)

Аддитивные штрафные сплайны и обобщенные аддитивные модели (GAM)

Как только у нас будет более одной ковариации, нам понадобятся аддитивные модели, чтобы преодолеть «проклятие размерности», оставаясь при этом разумными. В зависимости от представления гладких функций, GAM может принимать различные формы. Наиболее популярным, на мой взгляд, является mgcv пакет, рекомендованный R.

Вы все равно можете подогнать одномерный регрессионный сплайн с mgcv :

 library(mgcv)
fit <- gam(noisy.y ~ s(q, bs = "cr", k = 10))