library(spuRs) library(locfit) epanker <- function(x,h){ y <- ifelse(abs(x/h)<1,15*(1-(x/h)^2)^2/16/h,0) return(y) } n <- 200 m <- 150 u <- runif(n,0,1) h <- kdeb(u, h0 = 0.01 * sd(u), h1 = sd(u), meth ="LSCV", kern ="gauss", gf = 2.5) z <- rnorm(n,0,1) x <- rnorm(n,0,1) t <- runif(n,0,1) beta.y <- 2 beta <- 2.5 epsi <- rnorm(n,0,0.64) xt <- function(a){ f <- sin(2*pi*a) return(f) } eta <- xt(t)+epsi alpha <- function(a){ f <- 2.5*a^2+2 return(f) } y <- x*alpha(u)+(z*beta.y)^2+epsi ybo <- y-(z*beta)^2 ksi <- c(1,0) l <- rep(1,n) vt <- function(t0){ return(cbind(l,t-t0)) } w <- function(t0){ wl <- rep(0,n) for(j in 1:n){ wl[j] <- epanker((t[j]-t0),h) } return(diag(wl)) } gj1 <- rep(0,n) for(i in 1:n){ t0 <- t[i] L <- solve(t(vt(t0))%*%w(t0)%*%vt(t0)) gj1[i] <- t(ksi)%*%L%*%t(vt(t0))%*%w(t0)%*%eta } W <- gj1 w <- function(u0){ wl <- rep(0,n) for(j in 1:n){ wl[j] <- epanker((u[j]-u0),h) } return(diag(wl)) } eye <- diag(rep(1,4))/n^3 Wu <- function(u0){ ul <- c() for(j in 1:n){ ul[j] <- (u[j]-u0)/h } wu <- cbind(W,ul*W) return(wu) } gj <- matrix(0,2,n) gjl <- c() S <- matrix(0,n,n) Q <- matrix(0,n,n) for(i in 1:n){ u0 <- u[i] L <- t(Wu(u0))%*%w(u0)%*%Wu(u0) LN <- solve(L) gj[,i]<- LN%*%t(Wu(u0))%*%w(u0)%*%ybo S[i,] <- c(eta[i],0)%*%LN%*%t(Wu(u0))%*%w(u0) Q[i,] <- c(1,0)%*%LN%*%t(Wu(u0))%*%w(u0) } gjl <- c(1,0)%*%gj all <- gjl RASE <- sum((all-alpha(u))^2) f <- function(beta){ ybo <- y-S%*%y gbo <- (z*beta)^2-S%*%(z*beta)^2 g1 <- 2*(z*beta)*z g1bo <- g1-S%*%g1 f1 <- c() df1 <- c() for(i in 1:n){ f1[i] <-(ybo[i]-gbo[i])*g1bo[i]-t(g1)%*%Q[i,]* Q[i,]%*%(y-(z*beta)^2) df1[i] <- 0-(g1bo[i])^2+(ybo[i]-gbo[i])*2*z[i]^2-t(2*z^2)%*% Q[i,]*Q[i,]%*%(y-(z*beta)^2)+t(g1)%*%Q[i,]*Q[i,]%*%g1 } fx <- sum(f1)/n dfx <- sum(df1)/n return(c(fx,dfx)) } newtonraphson(f,beta,1e-10)