矩估计说明:
就是把求方程解的Newton法过渡到统计学里面,达到参数估计的目的,其实也是对于方程的求解。
直接上例题
R代码
#矩估计
#下面是之前的Newtons代码
Newtons<-function(fun,x,ep=1e-5,it_max=100){#itmax表示最大迭代次数
index<-0;k<-1
while(k<=it_max){
x1<-x;obj<-fun(x);
x<-x-solve(obj$j,obj$f)
norm<-sqrt((x-x1)%*%(x-x1))
if(norm<ep){
index<-1;break
}
k=k+1
}
obj<-fun(x)
list(root=x,it=k,index=index,Fval=obj$f)
}
#二项分布B(n,p),n*p-mean(x)=0,n*p*(1-p)-sum((x-mean(x))^2/m-1=0
bifun<-function(p){
f=c(p[1]*p[2]-mean(k),p[1]*p[2]*(1-p[2])-sum((k-mean(k))^2)/(length(k)-1))
j=matrix(c(p[2],p[1],p[2]*(1-p[2]),p[1]-2*p[1]*p[2]),nrow = 2,byrow = T)
list(f=f,j=j)
}
k=rbinom(1000,20,0.7)
Newtons(bifun,c(10,0.5))
结果
> #矩估计
> Newtons<-function(fun,x,ep=1e-5,it_max=100){#itmax表示最大迭代次数
+ index<-0;k<-1
+ while(k<=it_max){
+ x1<-x;obj<-fun(x);
+ x<-x-solve(obj$j,obj$f)
+ norm<-sqrt((x-x1)%*%(x-x1))
+ if(norm<ep){
+ index<-1;break
+ }
+ k=k+1
+ }
+ obj<-fun(x)
+ list(root=x,it=k,index=index,Fval=obj$f)
+ }
>
> #二项分布B(n,p),n*p-mean(x)=0,n*p*(1-p)-sum((x-mean(x))^2/m-1=0
> bifun<-function(p){
+ f=c(p[1]*p[2]-mean(k),p[1]*p[2]*(1-p[2])-sum((k-mean(k))^2)/(length(k)-1))
+ j=matrix(c(p[2],p[1],p[2]*(1-p[2]),p[1]-2*p[1]*p[2]),nrow = 2,byrow = T)
+ list(f=f,j=j)
+ }
>
> k=rbinom(1000,20,0.7)
> Newtons(bifun,c(10,0.5))
$`root`
[1] 20.3717921 0.6860467
$it
[1] 6
$index
[1] 1
$Fval
[1] 1.776357e-15 0.000000e+00
结果20.3717921 0.6860467,与20 0.7差不多!!!