简单线性回归问题的优化(SGD)R语言

本编博客继续分享简单的机器学习的R语言实现。

今天是关于简单的线性回归方程问题的优化问题

常用方法,我们会考虑随机梯度递降,好处是,我们不需要遍历数据集中的所有元素,这样可以大幅度的减少运算量。

具体的算法参考下面:

首先我们先定义我们需要的参数的Notation

上述算法中,为了避免过拟合,我们采用了L2的正则化,在更新步骤中,我们会发现,这个正则项目,对参数更新的影响

下面是代码部分:

## Load Library
library(ggplot2)
library(reshape2)
library(mvtnorm)

## Function for reading the data
read_data <- function(fname, sc) {
   data <- read.csv(file=fname,head=TRUE,sep=",")
   nr = dim(data)[1]
   nc = dim(data)[2]
   x = data[1:nr,1:(nc-1)]
   y = data[1:nr,nc]
   if (isTRUE(sc)) {
      x = scale(x)  ## Scale x
      y = scale(y)  ## Scale y
   }
   return (list("x" = x, "y" = y))
}

我们定义了一个读取数据的方程,这里,我们会把数据集给scale一下,可以保证进一步提高运算速度

## Matrix Product Function
predict_func <- function(Phi, w){
    return(Phi%*%w)
} 

## Function to compute the cost function
train_obj_func <- function (Phi, w, label, lambda){
    # Cost funtion including the L2 norm regulaztion
    return(.5 * mean((predict_func(Phi, w) - label)^2) + .5 * lambda * t(w) %*% w)
}

## Return the errors for each iteration
get_errors <- function(data, label, W) {
    n_weights = dim(W)[1]
    Phi <- cbind(‘X0‘ = 1, data)
    errors = matrix(,nrow=n_weights, ncol=2)
    for (tau in 1:n_weights) {
       errors[tau,1] = tau
       errors[tau,2] = train_obj_func(Phi, W[tau,],label, 0) ## Get the errors, set the lambda to 0
    }
   return(errors)
}

 同时,我们定义了计算矩阵乘法,计算目标函数以及求误差的方程。

sgd_train <- function(train_x, train_y, lambda, eta, epsilon, max_epoch) {

    ## Prepare the traindata
    ## Attach the 1 for X0
    Phi <- as.matrix(cbind(‘X0‘=1, train.data))

    ## Calculate the max iteration time for the SGD
    train_len = dim(train_x)[1]
    tau_max = max_epoch * train_len

    W <- matrix(,nrow=tau_max, ncol=ncol(Phi))
    set.seed(1234)
    ## Random Generate the start parameter
    W[1,] <- runif(ncol(Phi))

    tau = 1 # counter
    ## Create a dateframe to store the value of cost function for each iteration
    obj_func_val <-matrix(,nrow=tau_max, ncol=1)
    obj_func_val[tau,1] = train_obj_func(Phi, W[tau,],train_y, lambda)

    while (tau <= tau_max){

        # check termination criteria
        if (obj_func_val[tau,1]<=epsilon) {break}

        # shuffle data:
        train_index <- sample(1:train_len, train_len, replace = FALSE)

        # loop over each datapoint
        for (i in train_index) {
            # increment the counter
            tau <- tau + 1
            if (tau > tau_max) {break}

            # make the weight update
            y_pred <- predict_func(Phi[i,], W[tau-1,])
            W[tau,] <- sgd_update_weight(W[tau-1,], Phi[i,], train_y[i], y_pred, lambda, eta)

            # keep track of the objective funtion
            obj_func_val[tau,1] = train_obj_func(Phi, W[tau,],train_y, lambda)
        }
    }
   # resulting values for the training objective function as well as the weights
    return(list(‘vals‘=obj_func_val,‘W‘=W))
}

# updating the weight vector
sgd_update_weight <- function(W_prev, x, y_true, y_pred, lambda, eta) {
   ## Computer the Gradient
   grad = - (y_true-y_pred) * x
   ## Here I just combine the regularisation term with prev w
   return(W_prev*(1-eta * lambda) - eta * grad)
}

  根据上述我们写的计算更新目标函数和参数的方法,讲算法用R实现

下面是实验部分

## Load the train data and train label
train.data <- read_data(‘assignment1_datasets/Task1C_train.csv‘,TRUE)$x
train.label <- read_data(‘assignment1_datasets/Task1C_train.csv‘,TRUE)$y
## Load the testdata and test label
test.data <- read_data(‘assignment1_datasets/Task1C_test.csv‘,TRUE)$x
test.label <- read_data(‘assignment1_datasets/Task1C_test.csv‘,TRUE)$y

# Set MAX EPOCH

max_epoch = 18

## Implement SGD with Ridge regression
options(warn=-1)

## Initilize
## Set the related parameters
epsilon = .001 ## Terimation threshold
eta = .01 ## Learning Rate
lambda= 0.5 ## Regularisation parmater

## Run SGD
## Cost function values
train_res2 = sgd_train(train.data, train.label, lambda, eta, epsilon, max_epoch)
## Calulate the errors
## To be mentioned here, we will only visulisation for the train error to check the converge result
errors2 = get_errors(train.data, train.label, train_res2$W)  

 接着,我们把SGD的error plot给绘制出来

## Visulastion for SGD
plot(train_res2$val, main="SGD", type="l", col="blue",
     xlab="iteration", ylab="training objective function")

  

  

这里我们的方程比较简单,可以看到,目标函数很快就收敛了。

原文地址:https://www.cnblogs.com/chenyusheng0803/p/9657049.html

时间: 2024-10-10 21:30:31

简单线性回归问题的优化(SGD)R语言的相关文章

ArcGIS与R语言的Delaunay 三角网生成法

上一次讲罗马七丘的空间分析的时候,讲了一个Delaunay 三角测量--实际上这个东西也是泰森多边形的基础,泰森多边形就是在这个Delaunay三角上面进化出来的东西,所以很多同学在问我,这个东东怎么做啊?在ArcGIS里面有工具么?还是那句话:没有现成的工具--但是可以通过一系列手段来完成. 还有就是可以通过R语言来做Delaunay三角--这个是一个很神奇的东东--当然,python也行,高手自己写算法也行. 今天简单讲讲如何用ArcGIS和R语言来分别生成Delaunay三角网. 示例数据

R语言:用简单的文本处理方法优化我们的读书体验

前言 延续之前的用R语言读琅琊榜小说,继续讲一下利用R语言做一些简单的文本处理.分词的事情.其实就是继续讲一下用R语言读书的事情啦,讲讲怎么用它里面简单的文本处理方法,来优化我们的读书体验,如果读邮件和读代码也算阅读的话..用的代码超级简单,不涉及其他包 这里讲两个示例,结尾再来吐槽和总结. 1)R-Blogger订阅邮件拆分 2) R代码库快速阅读方法 不在博客园上阅读时才会看到的,这篇博文归 http://www.cnblogs.com/weibaar所有 仅保证在博客园博客上的排版干净利索

R语言解读多元线性回归模型

转载:http://blog.fens.me/r-multi-linear-regression/ 前言 本文接上一篇R语言解读一元线性回归模型.在许多生活和工作的实际问题中,影响因变量的因素可能不止一个,比如对于知识水平越高的人,收入水平也越高,这样的一个结论.这其中可能包括了因为更好的家庭条件,所以有了更好的教育:因为在一线城市发展,所以有了更好的工作机会:所处的行业赶上了大的经济上行周期等.要想解读这些规律,是复杂的.多维度的,多元回归分析方法更适合解读生活的规律. 由于本文为非统计的专业

R语言解读一元线性回归模型

前言 在我们的日常生活中,存在大量的具有相关性的事件,比如大气压和海拔高度,海拔越高大气压强越小:人的身高和体重,普遍来看越高的人体重也越重.还有一些可能存在相关性的事件,比如知识水平越高的人,收入水平越高:市场化的国家经济越好,则货币越强势,反而全球经济危机,黄金等避险资产越走强. 如果我们要研究这些事件,找到不同变量之间的关系,我们就会用到回归分析.一元线性回归分析是处理两个变量之间关系的最简单模型,是两个变量之间的线性相关关系.让我们一起发现生活中的规律吧. 由于本文为非统计的专业文章,所

R语言与数据分析之七:时间序列简单指数平滑

上篇我们对时间序列数列有了整体的认识并将时间序列进行了分解,今天和小伙伴们分享常用预测算法中相对最简单的:简单指数平滑法.简单指数平滑适用于可用相加模型描述,并且处于恒定水平和没有季节变动的时间序列地短期预测. 简单指数平滑法提供了一种方法估计当前时间点上的水平.为了更加准确的估计当前时间的水平,我们使用alpha参数来控制平滑,alpha的取值在0-1之间.当alpha越接近0,临近预测的观测值在预测中的权重就越小. 我们采用伦敦1813年到1912年全部的每年每英尺降雨量来做分析对象,首先读

多元线性回归公式推导及R语言实现

多元线性回归 多元线性回归模型 实际中有很多问题是一个因变量与多个自变量成线性相关,我们可以用一个多元线性回归方程来表示. 为了方便计算,我们将上式写成矩阵形式: Y = XW 假设自变量维度为N W为自变量的系数,下标0 - N X为自变量向量或矩阵,X维度为N,为了能和W0对应,X需要在第一行插入一个全是1的列. Y为因变量 那么问题就转变成,已知样本X矩阵以及对应的因变量Y的值,求出满足方程的W,一般不存在一个W是整个样本都能满足方程,毕竟现实中的样本有很多噪声.最一般的求解W的方式是最小

R语言确定聚类的最佳簇数:3种聚类优化方法

原文链接:http://tecdat.cn/?p=7275 确定数据集中最佳的簇数是分区聚类(例如k均值聚类)中的一个基本问题,它要求用户指定要生成的簇数k. 一个简单且流行的解决方案包括检查使用分层聚类生成的树状图,以查看其是否暗示特定数量的聚类.不幸的是,这种方法也是主观的. 我们将介绍用于确定k均值,k medoids(PAM)和层次聚类的最佳聚类数的不同方法. 这些方法包括直接方法和统计测试方法: 直接方法:包括优化准则,例如簇内平方和或平均轮廓之和.相应的方法分别称为弯头方法和轮廓方法

R语言中最简单的向量赋值方法

R语言中最简单的向量赋值方法简介: 1. 生成等差数列的向量x x <- 1:10 #将x向量赋值为1 2 3 4 5 6 7 8 9 10 结果为 > x [1] 1 2 3 4 5 6 7 8 9 10 2. 将x的值全部修改成0 x[] <- 0 #非常简洁的赋值方法,建议使用 x[1:length(x)] <- 0 #不建议使用的赋值方法 结果为: > x[] <- 0 > x [1] 0 0 0 0 0 0 0 0 0 0 3.使用seq函数 x <

用R语言实现一个求导的简单例子

在学习导数的时候,老师都会这样解释,假设y=f(x),对点(x0,y0)的求导过程如下:设dx是一个很小的数=> y0+dy=f(x0+dx)=> dy=f(x0+dx)-y0则在这一点的导数a=dy/dx这样假设的前提是dx很小,所以x0至(x0+dx)很短,可以近似为一条直线,则dy/dx可以看成是点(x0,y0)和点(x0+dx,y0+dy)连成直线的斜率因此关于某个点的导数,其意义也可以解释成在该点的变化率 下面用R语言实现一个简单例子:假设:y=1/x,求点(1,1)的导数 x<