aaa=Sys.time()
   #set.seed(23569)
   n=100
   total=100
   B=50
   h1=0.03
   h2=0.05
   #h1=n^(-1/5)
   #h2=0.5*1.02^40*sqrt(0.2)
   lamda=seq(0,2,0.25)
   N=length(lamda)
   sigmau=0.6
   g1=0
   g2=0
   RASE1=0
   RASE2=0

   for(run in 1:total)
      {
         #x=rnorm(n,0,1)
         x=runif(n,-1,1)
         e=rnorm(n,0,0.4)
         y=x^2+e
         u=rnorm(n,0,sigmau)
         v=matrix(rnorm(n*B,0,sigmau),n,B)
         w=x+u
         #t=sort(x)
         t=seq(-1,1,length=n)

         #SIMEX定义一个lmd的函数arg
         at=array(rep(0,2*n*N*B),dim=c(2,n,N,B))
         gs0=array(rep(0,n*N),dim=c(n,N))
         gs00=array(rep(0,n*N),dim=c(n,N))

            for (j in 1:N)
              {
              for (b in 1:B)
                 {
                 for (i in 1:n)
                   {
                       objfun=function(a)
                           {
                              gs0=a[1]
                              gs00=a[2]
                              #wb[j]=w+sqrt(lamda[j])*v[,b]
                              -sum((1/(n*h))*dnorm((w+sqrt(lamda[j])*v[,b])-t,0,h1)*dnorm(y-gs0-gs00*((w+sqrt(lamda[j])*v[,b])-t),0,h2))
                           
                            }
                       at[,i,j,b]=optim(c(t^2-0.5,2*t),objfun)$par
                       gs0[i,j]=mean(at[1,i,j,])
                       gs00[i,j]=mean(at[2,i,j,])

                   }
                 }
               }



        #外推
        p0=array(rep(0,n),dim=c(n,1))
        p1=array(rep(0,n),dim=c(n,1))
        p2=array(rep(0,n),dim=c(n,1))
        gs0m=array(rep(0,n),dim=c(n,1))
        gs00m=array(rep(0,n),dim=c(n,1))

        for (i in 1:n)
            {   
                 pf=(lamda)^2
                 lm.sol=lm(gs0[i,]~lamda+pf)
                 p0=lm.sol$coefficient[1]
                 p1=lm.sol$coefficient[2]
                 p2=lm.sol$coefficient[3]
                 gs0m=p0-p1+p2
            
                 #pf=(lamda)^2
                 #lm.sol=lm(gs00[i,]~lamda+pf)
                 #p0=lm.sol$coefficient[1]
                 #p1=lm.sol$coefficient[2]
                 #p2=lm.sol$coefficient[3]
                 #gs00m=p0-p1+p2
             }

            

        #naive估计

       bt=array(rep(0,2*n),dim=c(2,n))
       gn=rep(0,n)            
              for (i in 1:n)
                 {
                   nivfun=function(b)
                           {
                              gs1=b[1]
                              gs11=b[2]
                              -sum((1/(n*h))*dnorm(w-t,0,h1)*dnorm(y-gs1-gs11*(w-t),0,h2))
                           
                            }
                       bt[,i]=optim(c(t^2-0.5,2*t),nivfun)$par
                       #gn=matrix(c(1,0),1,2)%*%bt
                       gn=bt[1,i]
                   }
                 
        gt=t^2
        g1=g1+gs0m
        g2=g2+gn
        RASE11=sqrt(sum((gt-gs0m)^2)/n)
        RASE21=sqrt(sum((gt-gn)^2)/n)
        RASE1=RASE1+RASE11
        RASE2=RASE2+RASE21
    }

  gss=g1/total
  gww=g2/total
  RASE=RASE1/total
  RASEn=RASE2/total


  plot(x,y,ylim=c(-0.5,1),xlim=c(-1,1))
  lines(t,gt,col="red")
  lines(t,gss,col="blue")#SIMEX
  lines(t,gww,col="green")#naive

  RASE
  RASEn

  bbb=Sys.time()
  aaa
  bbb