Цель этой статьи — продемонстрировать различные подходы, которые мы можем использовать для построения одномерной линейной регрессии в R. Мы увидим формулы, используемые для каждого подхода, и сравним их предсказания.
Введение
Линейная регрессия — это метод прогнозирования непрерывной переменной путем нахождения линейной связи между независимыми переменными и зависимыми переменными.
Линейное уравнение:
Линейная регрессия с одной независимой переменной называется одномерной линейной регрессией. Для одномерной линейной регрессии m всегда равно 1 (т. е. m = 1)
Есть два подхода к прогнозированию в LR. Они есть
- Нормальное уравнение
- Градиентный спуск
Нормальное уравнение
Есть два способа добиться предсказания с помощью нормального уравнения.
- Матричная форма
- Ковариация/дисперсия
Матричная форма
Ковариация/дисперсия
Градиентный спуск
Любой вариант градиентного спуска можно использовать для предсказания зависимой переменной.
Функция гипотезы:
Функция стоимости:
Градиентный спуск:
Импорт библиотек и загрузка данных
Установка библиотек
# norm_pkgs <- c('moments','nortest')
# install.packages(norm_pkgs)
# install.packages('matlib')
library(MASS) # for matrix operationsЗагрузка данных
data <- read.csv('input/ex1data1.txt',header=FALSE, col.names = c('population','profit')) head(data)## population profit ## 1 6.1101 17.5920 ## 2 5.5277 9.1302 ## 3 8.5186 13.6620 ## 4 7.0032 11.8540 ## 5 5.8598 6.8233 ## 6 8.3829 11.8860X = as.matrix(data$population) y = as.matrix(data$profit)
Преобразование переменных X и Y в массивы numpy (матрицы/векторы), а также добавление столбца 1 к X (поскольку значение X0 равно 1 в w0. x0, т.е. предполагается, что коэффициент перехвата равен 1)
Нормальное уравнение — ковариация/дисперсия
getCoeffs <- function(X,y) { covariance_xy = cov(X,y) variance_x = var(X)coeffs = covariance_xy/variance_x y_intercept = mean(y) - coeffs * mean(X) return(list('intercept' = y_intercept[1], 'coefficient'=coeffs[1])) }getPred <- function(X,y) { res <- getCoeffs(X,y) print(res) preds = res$intercept +(res$coefficient * X) return(preds) }preds = getPred(X,y)## $intercept ## [1] -3.895781 ## ## $coefficient ## [1] 1.193034## Model EvaluationSSE = sum((preds-y)^2) # Sum of squared error SST = sum((y-mean(y))^2) # sum of squared total n = nrow(X) # number of observations q = 1 # number of coefficientsMSE = SSE/(n-q) # Mean squared error MST = SST/(n-1) # Mean squared total MAPE = sum(abs(preds-y)^2)/n # Mean Absolute Percentage ErrorR_squared = 1-(SSE/SST) # R Square A_rsquare = 1-(MSE/MST) # Adjusted R Square std_err = MSE^1/2 # Standard Error/Root Mean Squared error MSR = (SST-SSE)/(q-1) f_static = MSR/MSE # F static# paste('Mean Squared Percentage Error :',MSE) # paste('Mean Absolute Percentage Error :',MAPE) # print(paste('R Square :',R_squared)) # print(paste('Adjusted R Square :',A_rsquare)) # print(paste('Standard Error :',std_err)) print(paste('F Statistic :',f_static))## [1] "F Statistic : Inf"
Среднеквадратическая ошибка в процентах: 9,047213
Средняя абсолютная ошибка в процентах: 8,9539428
Квадрат R: 0,7020316
Скорректированный квадрат R: 0,7020316
Стандартная ошибка: 4,5236065
F Статистика :
lm() в пакете статистики
lrModel <- lm(profit~population, data=data) print(lrModel)## ## Call: ## lm(formula = profit ~ population, data = data) ## ## Coefficients: ## (Intercept) population ## -3.896 1.193summary(lrModel)## ## Call: ## lm(formula = profit ~ population, data = data) ## ## Residuals: ## Min 1Q Median 3Q Max ## -5.8540 -1.9686 -0.5407 1.5360 14.1982 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) -3.89578 0.71948 -5.415 4.61e-07 *** ## population 1.19303 0.07974 14.961 < 2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 3.024 on 95 degrees of freedom ## Multiple R-squared: 0.702, Adjusted R-squared: 0.6989 ## F-statistic: 223.8 on 1 and 95 DF, p-value: < 2.2e-16
Наблюдение:
P-значение модели является значительным.
Обе оценки, R-квадрат и скорректированный R-квадрат, выше и лучше Независимая переменная «популяция» значима
modsummary <- summary(lrModel) modcoeffs <- modsummary$coefficientsstd.err <- modcoeffs['population','Std. Error'] beta.estim <- modcoeffs['population','Estimate'] t_valve <- modcoeffs['population','t value'] # print(paste('t-statistic :',t_valve)) # print(paste('Std. Error :',std.err)) # print(paste('Beta Estimate/Coefficient :',beta.estim))
t-статистика: 14,9608056
стандарт. Ошибка: 0,0797439
Бета-оценка/коэффициент: 1,1930336
Наблюдение:
t-статистика хорошая (должно быть больше 1,96, чтобы p-значение было меньше 0,05), а стандартная ошибка меньше (близка к 0)
print(paste('AIC :',AIC(lrModel)))## [1] "AIC : 493.907190165596"print(paste('BIC :',BIC(lrModel)))## [1] "BIC : 501.631323101106"
Градиентный спуск
X = cbind(1,matrix(data$population)) y = as.matrix(data$profit)gradDescent <- function(x,y,alpha=0.005,maxIter=100,threshold=0.0001) { thetas <- list(maxIter) costs <- double(maxIter) converged <- FALSE theta <- matrix(c(0,0),nrow = 2) m <- nrow(y) i <- 1 while(converged==FALSE) { preds <- x %*% theta theta <- theta - alpha * (t(x) %*% (preds-y))/m cost = mean((preds-y)^2)/2 if (i > 1) { if ((costs[i-1]-cost) == threshold) { converged <- TRUE }} if (i == maxIter) { converged <- TRUE } costs[i] <- cost thetas[[i]] <- theta i <- i+1 } print(paste('theta :',thetas[which.min(costs)])) print(paste('Cost :',min(costs))) print(paste('Iteration :',which.min(costs))) #return(thetas[which.min(costs)]) }gradDescent(X,y,maxIter=15000)## [1] "theta : c(-3.89577556561726, 1.1930331104718)" ## [1] "Cost : 4.47697137597775" ## [1] "Iteration : 15000"
Вывод
Как мы видели, все подходы линейной регрессии для прогнозирования абсолютно одинаковы (с полным ядром, доступным в Kaggle), и я оставляю профилирование подходов (кодов), чтобы найти оптимизированный подход для будущих целей. Другие сообщения из этой серии с использованием Python и Pyspark
Ссылка:
https://sebastianraschka.com/faq/docs/closed-form-vs-gd.html
