牛顿迭代算回归系数

上次作业（求beta^hat）小高的简便做法 然而我还没有研究明白 优化，优化，优化 一个好的优化方法真的很难，但重要的在于知道有这么一种优化算法

#### 2. successful example ####
## data
attach(swiss)
x1 = swiss$Agriculture
x2 = swiss$Education
x3 = swiss$Catholic
y = swiss$Fertility
data = data.frame(x1 , x2, x3, y)

## build matrix (maths work)
## mat is beta vector
f <- function(m = matrix(),
              y = vector(mode = 'numeric'),
              mat = vector(mode = 'numeric')){
  out <- t(y - m%*%mat)%*%(y - m%*%mat)
  return(out)
}

# first-order derivative matrix
## m is independent variable(matrix) dim = 3
## y is dependent variable(vector) dim = 1
d <- function(m = matrix(),
              y = vector(mode = 'numeric'),
              mat = vector(mode = 'numeric')){
  out <- -2*t(m)%*%y+2*t(m)%*%m%*%mat
  return(out)
}

# hession
dd <- function(m = matrix()){
  out <- 2*t(m)%*%m
  return(out)
}

# main function to deal with this problem
mat_esti <- function(m, y){
  MaxIter = 30 # loop times
  m <- as.matrix(m)
  y <- as.matrix(y)
  
  m <- cbind(rep(1,length(m[,1])), m) # add a col of '1'
  mat <- seq(1,length(m[1,]))
  mat1 <- matrix(mat,1,length(mat))
  
  # Newton
  for(i in 1:MaxIter){
    # 'd' & 'dd' are same-wide matrixs
    # a %*% x = b  ->  solve(a,b)=x
    mat <- mat - solve(dd(m), d(m, y, mat)) 
    mat1 <- rbind(mat1, matrix(mat,1,length(mat)))
    if(sum((mat1[i+1,] - mat1[i, ])^2) <1e-6){break}
  }
  # every beta will be stored in mat1
  #
  # write colnames, and make a better output
  names <- c('Intercept','beta 1','beta 2','beta 3')
  dim(mat) <- c(1, length(mat))
  best_mat <- as.data.frame(mat)
  for(i in 1:(length(mat)-1)){
    name <- paste(names[i],i)
  }
  colnames(best_mat) <- names
  out <- data.frame(best_mat)
  return(out)
}

## output

mat_esti(data[,1:3],data[,4])
##   Intercept     beta.1    beta.2    beta.3
## 1  86.22502 -0.2030377 -1.072147 0.1452013
## beta = (XX')^(-1)X'Y

同学的模型

生成数据上值得借鉴

# DGP 数据生成过程
## 给定a，b，构造出y和x来。判断估计出的a^,b^，和真实的相比差距有多大
## 1.
a = 5
b = 3
set.seed(250)
## 2. 生成因变量k
k = rnorm(100,1,0)
## 3. 生成自变量xi
x1 = runif(100)
## 4. 方程
k = a+b*x1
## 5. 估计
# {k,x1} -> b
# 比较估计的b和设定的b之间的误差

总结！！

多元回归
多变量回归（我们通常所说） 多元回归是因变量y有很多列

对于一个模型的三个重要组件

yi = b0 + b1x1 + b2x2 + … + bpxp + e
第一部分：数据（dataset）
知道y，知道xi
第二部分：参数
链接y与x关系的变量（b0~bp）
第三部分：扰动误差项error

三部分都不能少，包括误差项（模型具有不确定性）

建模过程中，误差项往往要假定已知。

高斯-马尔科夫假定（重要七个假定）
约定误差e服从N(0,sigma^2)
遍有了最小二乘的最基本思路sum( e^2)

似然函数

最像的函数，找这样的函数，使得y和x提供的信息最充分
对线性回归而言，有误差项，因此能够得到y的分布
yi~ N(b0+b1x1+b2x2+…+bpxp, sigma^2)
因此可以构建似然函数l
关于yi密度的乘积，观测从1到n

极大似然估计就是在最大化l函数
找到合适的b，使其最大化
到了这一步，我们就可以牛顿迭代了！
其实极大化似然函数，等价于极大化log(l) :b0~bp

如何处理一个统计模型

1.写出模型

2.把我们关心的三个组件写出来

数据，参数，误差项

3.选一个最好的方式，OLS或者极大似然估计等

明白等式左边的y的表达式，分布

4.写出来似然函数，或者对数似然函数

5.极大化时，对log(l)求关于b的偏导数

sum(((yi-b0-b1x1-…-bpxp)/sigma^2)*((-1,x1,…,xp)^(-1)))
i = 1~n

求log(l)对b的二阶导（应该是一个p+1阶的方阵）

homework

写对数一阶导
写对数二阶导
再用dgp法生成数据检验