以下的code是用牛顿算法来找使得目标函数最大化的参数b ，原方程为y[i]=a+b*x[i]+e[i]  假设a=1  ，  e~N（0,1）iid
即通过极大似然估计的方法估计参数b
但是若a也是未知数，能否用该方法同时估计a，b  使得目标函数最大化   


##目标函数
f<-function(x,y,b)
{
sum(-0.5*log(2*pi)-(y-1-b*x)^2/2)
}
##一阶导数
df<-function(x,y,b)
{
sum((y-1-b*x)*x)
}
##二阶导数
d2f<-function(x,y,b)
{
sum(-x^2)
}
##循环找到最好的b
n=10
b=rep(NA,n)
x=matrix(c(1,2,3,4,5),5,1)
y=matrix(c(3,5,7,9,11),5,1)
b[1]=5   #赋给b的初始值
for(i in 1:n)
{
b[i+1]<-b[i]-(d2f(x,y,b[i]))^(-1)*df(x,y,b[i])
}
## 画图 Check the convergence
par(mfrow = c(4, 1))
plot(b,ylim=c(0,5),type = "l", col = "red")

p=f(x,y,b[1])
for(i in 2:n)
p=c(p,f(x,y,b[i]))
plot(p,type="l", col = "red")

q=df(x,y,b[1])
for(i in 2:n)
q=c(q,df(x,y,b[i]))
plot(q,type = "l", col = "red")            

plot(rep(d2f(x,y,b),n), type = "l", col = "red")