Estimators

######################################MLR PORTION

#our goal is to write a newton update proceedure to maximize the loglikelihood

#we seek to find the 0 of the first derivative (if our function has only one extrema and it is a maximum #then that 0 for the 1st derivative will be the maximum of our function

#B1 = Bo - inverse(2nd derivative at Bo)*first derivative at Bo #note here Bo is our staring Bvector. B1 is the updated Bvector which becomes Bo in the next step set.seed(99999) trueB0 <- 05 trueB1 <- (-10) trueSd <- 3 sampleSize <- 1000

#Simulate data x = runif(sampleSize) y = trueB0 + trueB1 * x + rnorm(n=sampleSize,mean=0,sd=trueSd)

lm(y~x) bstart=vector(length=3) bstart[1:3]=c(0,0,1) #not informative start

b0=bstart[1] b1=bstart[2] var=bstart[3] bOld=c(b0,b1,var)

change=1 i=1

bstore=matrix(ncol=3,nrow=4000) n=length(y)

while(change>.0001) #while(i<4000) {

#we’ll do the sums ourself to be sure they’re set up right

#sum for 1st element of 1st derivative (note this is literally the sum part, so it will omit any coefficients in front of the summand) sum1=0 sum2=0 sum3=0 for(l in 1:length(y) ) { yi=y[l] xi=x[l] sum1=sum1 + 2yi - 2b0 - 2xib1

  sum2=sum2 + 2*yi*xi - 2*b0*xi - 2*b1*xi^2

 sum3=sum3 + yi^2 - 2*yi*b0 - 2*yi*xi*b1 + b0^2 + 2*b0*xi*b1 + (xi^2)*(b1^2)
}

#apply any coefficients, etc. to the summand
sum1=sum1*(1/(2*var))
sum2=sum2*(1/(2*var))
 sum3=sum3*(1/(2*var^2)) - n/(2*abs(var))

firstDerivVec=c(sum1,sum2,sum3)

#now we need 2nd derivative matrix (hessian) oneOne=0 oneTwo=0 oneThree=0

twoOne=0
twoTwo=0
twoThree=0

threeOne=0
threeTwo=0
threeThree=0


#hessian order is [db0 of db0      db0 of db1     db0 of dvar ]
#                 |db1 of db0      db1 of db1     db1 of dvar |
#                 [dvar of db0    dvar of db1     dvar of dvar]

for(l in 1:length(y))
{
  yi=y[l]
  xi=x[l]
  
   #oneOne= #no sum needed here, it simplified to remove the sum
  oneTwo=oneTwo + xi
  oneThree= oneThree -2*yi + 2*b0 + 2*xi*b1
  
  
  twoOne=twoOne + xi
  twoTwo=twoTwo + xi^2
  twoThree=twoThree + b0*xi + b1*xi^2 - yi*xi

  threeOne=threeOne -2*yi + 2*b0 + 2*xi*b1
  threeTwo=threeTwo + b0*xi + b1*xi^2 - yi*xi
  threeThree=threeThree + yi^2 - 2*yi*b0 - 2*yi*xi*b1 + b0^2 + 2*b0*xi*b1 + (xi^2)*(b1^2)
}

#apply any coefficients, etc. to the sums
oneOne=-n/var
oneTwo=oneTwo*(-1/var)
oneThree=oneThree*(1/(2*var^2))

twoOne=twoOne*(-1/var)
twoTwo=twoTwo*(-1/var)
twoThree=twoThree * (1/(var^2))

threeOne=threeOne*(1/(2*var^2))
threeTwo=threeTwo * (1/(var^2))
threeThree=threeThree*(-1/(abs(var)^3)) +n/(2*var^2)





secondDerivMat=matrix(nrow=3,ncol=3)
secondDerivMat[1,1]=oneOne
secondDerivMat[1,2]=oneTwo
secondDerivMat[1,3]=oneThree

secondDerivMat[2,1]=twoOne
secondDerivMat[2,2]=twoTwo
secondDerivMat[2,3]=twoThree

secondDerivMat[3,1]=threeOne
secondDerivMat[3,2]=threeTwo
secondDerivMat[3,3]=threeThree

bnew=bOld-solve(secondDerivMat)%*%firstDerivVec
change=sum(abs(abs(bOld)-abs(bnew)))
b0=bnew[1]
b1=bnew[2]
var=bnew[3]
bOld=c(b0,b1,var)
i=i+1
bstore[i,]=bnew

}

summary(lm(y~x)) bnew