Essa é uma revisão anterior do documento!

---

idados<- 50
adados<-20000
 
#Equação verdadeira de relacionamento (usada para simular Fase I e gerar Modelo de Regressão)
coef0 <- c(10, 2, -3, 2.3)
 
#Equação alterada
coefal<-c(10, 6, -3, 2.3)
 
#Fase I
x1 <-runif(idados, min=-1, max=1)
x2 <-runif(idados, min=-1, max=1)
 
e <-rnorm(idados, mean=0, sd=4)
 
x <-cbind(1, x1, x2, x1*x2)
 
X <- model.matrix(~ x1*x2)
X
 
model.matrix(~ I(x1*x2))
 
 
#invx<-solve(crossprod(x)) #operação (X'X)^-1
da <-data.frame (x1, x2, x1*x2)
 
da <-data.frame (x1=x1, x1=x2)
 
y  <- vector(length=idados)
ii<-1:idados
for (ii in 1:idados) {
  y[ii]<- coef0%*%x[ii,]+e[ii]
}
 
yy <- as.vector(X %*% coef0 + e)
all.equal(y, yy)
 
da$y  <-y
da$y <- as.vector(X %*% coef0 + e)
 
rm(x1, x2, y, yy)
 
mod   <-lm(y~ x1 + x2 + x1*x2, data=da)
 
mod   <-lm(y~ x1 + x2 + I(x1*x2), data=da)
mod   <-lm(y~ x1 + x2 + x1:x2, data=da)
mod   <-lm(y~ x1*x2, data=da)
 
sdmod <-sqrt(sum(residuals(mod)^2)/df.residual(mod))
 
sdmod <- summary(mod)$sigma
 
# Modelo, com a retirada de uma rodada de outliers...
mod  <-lm(y~ x1 + x2 + x1*x2, data=da[abs(residuals(mod))<3*sdmod,])
sdmod<-sqrt(sum(residuals(mod)^2)/df.residual(mod))
 
#Fase II
# Amostra monitorada
xs  <-vector(length=adados)
xs1 <-runif(adados, min=-1, max=1)
xs2 <-runif(adados, min=-1, max=1)
 
es  <-rnorm(adados, mean=0, sd=4)
xs  <-cbind(1, xs1, xs2, xs1*xs2)
 
res <-vector(length=adados)
j   <- 1:adados
res[j] <- 0
ypred <- predict(mod, newdata=data.frame(x1=xs1, x2=xs2))
 
for (j in 1:adados) {
 res[j]<- (coefal%*%xs[j,]+es[j])-ypred[j] 
}
 
 
## 
solve(t(X) %*% X) %*% t(X) %*% y
 
solve(crossprod(X), crossprod(X,y))