dat<-data.frame(x=c(2,3,4,5,6),y=c(4,4,6,6,10))

N = length(dat$x)
x<-dat$x-mean(dat$x)
y<-dat$y-mean(dat$y)

ols<-lm(y~x-1)
summary(ols)

min.RSS <- function(pa1, pa2, x, y, N) {
    b0 = pa1
    b1 = pa2
    RSS = vector("numeric",length(b1))
    for (i in 1:length(b1)){
      RSS[i] = sum((pa2[i]*x-y)^2)
    }
    return(RSS)
}

pa1 <- seq(1, 1, by=0)
names(pa1) <- pa1
pa2 <- seq(1, 2, by=0.05)
names(pa2) <- pa2

z.RSS <- outer(X = pa1, Y = pa2, FUN = min.RSS, N = N, x = x, y = y)
min(z.RSS)
zd.RSS<-data.frame(z.RSS)

inds_0 = which(zd.RSS == min(zd.RSS), arr.ind=TRUE)
inds_0
(rnames = rownames(zd.RSS)[inds_0[,1]])
(cnames = colnames(zd.RSS)[inds_0[,2]])


Lik <- function(para1, para2,x, y, N){
  be1 = para1  
  sig2= para2
  Likelihood = vector("numeric",length(be1))
  for (i in 1:length(be1)){
    Likelihood[i] = (2*pi*sig2[i])^(-N/2)*exp((-1/(2*sig2[i]))*sum((y-be1[i]*x)^2))
 }
 return(Likelihood)
}

para1 <- seq(1.1, 2.1, by=0.05)
names(para1) <- para1
para2 <- seq(0.5, 0.9, by=0.01)
names(para2) <- para2

z.Lik <- outer(X = para1, Y = para2, FUN = Lik, N = N, x = x, y = y)
max(z.Lik)
zd.Lik<-data.frame(z.Lik)

inds = which(zd.Lik == max(zd.Lik), arr.ind=TRUE)
inds
(rnames = rownames(zd.Lik)[inds[,1]])
(cnames = colnames(zd.Lik)[inds[,2]])

require(lattice)
#Loading required package: lattice
#wireframe(z, drape=T, xlab="be1", ylab="sig2", zlab="Lik",col.regions=rainbow(100))
wireframe(z.Lik, xlab="be1", ylab="sig2", zlab="Lik")


LogL <- function(param1, param2, x, y, N){
 beta1 = param1
 sigma2 = param2
  LogLikelihood = vector("numeric",length(beta1))
  for (i in 1:length(beta1)){
    LogLikelihood[i] = (N/2)*log(sigma2[i])+1/(2*sigma2[i])*(sum((y-beta1[i]*x)^2))
      }
  return(LogLikelihood)
}

param1 <- seq(1.1, 2.1, by=0.05)
names(param1) <- param1
param2 <- seq(0.5, 0.9, by=0.01)
names(param2) <- param2

z.LogL <- outer(X = param1, Y = param2, FUN = LogL, N = N, x = x, y = y)
min(z.LogL)
zd.LogL<-data.frame(z.LogL)

inds_1 = which(zd.LogL == min(zd.LogL), arr.ind=TRUE)
inds_1
(rnames_1 = rownames(zd.LogL)[inds[,1]])
(cnames_1 = colnames(zd.LogL)[inds[,2]])

#require(lattice)
#Loading required package: lattice
#wireframe(z, drape=T, xlab="beta1", ylab="sigma2", zlab="LogL",col.regions=rainbow(100))
wireframe(z.LogL, xlab="beta1", ylab="sigma2", zlab="LogL")

# test it out
#LogL(param1[3], param2[3],  x = x, y = y, N = N)
#z.LogL[3,3]

library(bbmle)
fn<-function(b, sigma) {
  (N/2)*log(sigma^2)+1/(2*sigma^2)*(sum((y-b*x)^2))
}
mle<-mle2(fn,start=list(b=1.5, sigma=1))
summary(mle)