title: “Aneirson_Dinâmica” output: html_notebook
library(rsm)
library(lmtest)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(carData)
library(car)
library(nortest)
library(NlcOptim)## Loading required package: MASSlibrary(FrF2)## Loading required package: DoE.base## Loading required package: grid## Loading required package: conf.design## Registered S3 method overwritten by 'DoE.base':
## method from
## factorize.factor conf.design##
## Attaching package: 'DoE.base'## The following objects are masked from 'package:stats':
##
## aov, lm## The following object is masked from 'package:graphics':
##
## plot.design## The following object is masked from 'package:base':
##
## lengths#library(ggthemes)
library(rio)##
## Attaching package: 'rio'## The following object is masked from 'package:conf.design':
##
## factorizeDOE<-bbd(k=~x1+x2+x3,
block = F,
n0=5,
randomize = F,
coding = list(x1~(M-3)/2,
x2~(Mart-3)/1,
x3~(fixa-3)/1))
DOE## run.order std.order M Mart fixa
## 1 1 1 1 2 3
## 2 2 2 5 2 3
## 3 3 3 1 4 3
## 4 4 4 5 4 3
## 5 5 5 1 3 2
## 6 6 6 5 3 2
## 7 7 7 1 3 4
## 8 8 8 5 3 4
## 9 9 9 3 2 2
## 10 10 10 3 4 2
## 11 11 11 3 2 4
## 12 12 12 3 4 4
## 13 13 13 3 3 3
## 14 14 14 3 3 3
## 15 15 15 3 3 3
## 16 16 16 3 3 3
## 17 17 17 3 3 3
##
## Data are stored in coded form using these coding formulas ...
## x1 ~ (M - 3)/2
## x2 ~ (Mart - 3)/1
## x3 ~ (fixa - 3)/1DOE$y<-c(0.48,4.37,0.50,6.85,0,4.95,2.59,7.66,
2.03,1.83,3.51,6.31,3.25,3.44,3.37,3.68,3.51)
DOE## run.order std.order M Mart fixa y
## 1 1 1 1 2 3 0.48
## 2 2 2 5 2 3 4.37
## 3 3 3 1 4 3 0.50
## 4 4 4 5 4 3 6.85
## 5 5 5 1 3 2 0.00
## 6 6 6 5 3 2 4.95
## 7 7 7 1 3 4 2.59
## 8 8 8 5 3 4 7.66
## 9 9 9 3 2 2 2.03
## 10 10 10 3 4 2 1.83
## 11 11 11 3 2 4 3.51
## 12 12 12 3 4 4 6.31
## 13 13 13 3 3 3 3.25
## 14 14 14 3 3 3 3.44
## 15 15 15 3 3 3 3.37
## 16 16 16 3 3 3 3.68
## 17 17 17 3 3 3 3.51
##
## Data are stored in coded form using these coding formulas ...
## x1 ~ (M - 3)/2
## x2 ~ (Mart - 3)/1
## x3 ~ (fixa - 3)/1funcao<-rsm(y~SO(x1,x2,x3),data = DOE) # Função Empírica
summary(funcao)##
## Call:
## rsm(formula = y ~ SO(x1, x2, x3), data = DOE)
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.450000 0.068609 50.2852 3.221e-10 ***
## x1 2.532500 0.054240 46.6907 5.405e-10 ***
## x2 0.637500 0.054240 11.7533 7.308e-06 ***
## x3 1.407500 0.054240 25.9495 3.228e-08 ***
## x1:x2 0.615000 0.076707 8.0175 8.988e-05 ***
## x1:x3 0.030000 0.076707 0.3911 0.707355
## x2:x3 0.750000 0.076707 9.7775 2.481e-05 ***
## x1^2 -0.010000 0.074765 -0.1338 0.897362
## x2^2 -0.390000 0.074765 -5.2164 0.001231 **
## x3^2 0.360000 0.074765 4.8151 0.001932 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Multiple R-squared: 0.9978, Adjusted R-squared: 0.995
## F-statistic: 355.5 on 9 and 7 DF, p-value: 1.8e-08
##
## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## FO(x1, x2, x3) 3 70.408 23.4694 997.1817 1.43e-09
## TWI(x1, x2, x3) 3 3.766 1.2555 53.3445 3.45e-05
## PQ(x1, x2, x3) 3 1.128 0.3759 15.9697 0.001635
## Residuals 7 0.165 0.0235
## Lack of fit 3 0.062 0.0206 0.7994 0.555345
## Pure error 4 0.103 0.0257
##
## Stationary point of response surface:
## x1 x2 x3
## -10.984136 -4.637717 3.333766
##
## Stationary point in original units:
## M Mart fixa
## -18.968271 -1.637717 6.333766
##
## Eigenanalysis:
## eigen() decomposition
## $values
## [1] 0.55064213 0.07547313 -0.66611527
##
## $vectors
## [,1] [,2] [,3]
## x1 0.2587231 0.8802929 0.3976767
## x2 0.4295054 0.2639190 -0.8636387
## x3 0.8652095 -0.3942476 0.3098086funcao<-rsm(y~ FO(x1,x2,x3)+PQ(x2,x3)+TWI(x1,x2)
+TWI(x2,x3),data = DOE)
summary(funcao)##
## Call:
## rsm(formula = y ~ FO(x1, x2, x3) + PQ(x2, x3) + TWI(x1, x2) +
## TWI(x2, x3), data = DOE)
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.445789 0.054414 63.3253 3.081e-13 ***
## x1 2.532500 0.048415 52.3078 1.714e-12 ***
## x2 0.637500 0.048415 13.1673 3.479e-07 ***
## x3 1.407500 0.048415 29.0714 3.287e-10 ***
## x2^2 -0.390526 0.066643 -5.8599 0.0002408 ***
## x3^2 0.359474 0.066643 5.3940 0.0004364 ***
## x1:x2 0.615000 0.068470 8.9821 8.678e-06 ***
## x2:x3 0.750000 0.068470 10.9538 1.668e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Multiple R-squared: 0.9978, Adjusted R-squared: 0.996
## F-statistic: 573.6 on 7 and 9 DF, p-value: 3.411e-11
##
## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## FO(x1, x2, x3) 3 70.408 23.4694 1251.5443 4.133e-12
## PQ(x2, x3) 2 1.127 0.5636 30.0537 0.0001038
## TWI(x1, x2) 1 1.513 1.5129 80.6779 8.678e-06
## TWI(x2, x3) 1 2.250 2.2500 119.9850 1.668e-06
## Residuals 9 0.169 0.0188
## Lack of fit 5 0.066 0.0132 0.5108 0.7605417
## Pure error 4 0.103 0.0257
##
## Stationary point of response surface:
## x1 x2 x3
## -9.117565 -4.117886 2.338022
##
## Stationary point in original units:
## M Mart fixa
## -15.235129 -1.117886 5.338022
##
## Eigenanalysis:
## eigen() decomposition
## $values
## [1] 0.54425291 0.09339524 -0.66870078
##
## $vectors
## [,1] [,2] [,3]
## x1 0.2422866 0.8854479 -0.3965844
## x2 0.4288299 0.2689321 0.8624270
## x3 0.8702885 -0.3790218 -0.3145479persp(funcao,~x1+x3,zlab="Distância",
col = rainbow(50),
contours = ("colors"))
persp(funcao,~x2+x3,zlab="Distância",
col = rainbow(50),
contours = ("colors"))
funcao$residuals # Resíduos## 1 2 3 4 5 6
## -0.020263158 0.034736842 -0.045263158 0.009736842 0.134736842 0.019736842
## 7 8 9 10 11 12
## -0.090263158 -0.085263158 -0.089736842 -0.064736842 0.075263158 0.100263158
## 13 14 15 16 17
## -0.195789474 -0.005789474 -0.075789474 0.234210526 0.064210526curve(dnorm(x,mean(funcao$residuals),sd(funcao$residuals)),
from = min(funcao$residuals),to=max(funcao$residuals),
xlab = "Observações", ylab = "Densidade", lwd=4,
col="turquoise")
lines(density(funcao$residuals),col="blue", lwd=3,pch=2)
legend("topright",legend = c("Real","Estimada"),
col = c("turquoise","blue"),lwd=3,pch=2)
hist(funcao$residuals,lwd=3,col="yellow",main="Resíduos")
curve(dnorm(x,mean(funcao$residuals),sd(funcao$residuals)),
add = T,col="red",lwd=3,pch=5)
shapiro.test(funcao$residuals)##
## Shapiro-Wilk normality test
##
## data: funcao$residuals
## W = 0.97465, p-value = 0.8935residualPlot(funcao)
bptest(funcao)##
## studentized Breusch-Pagan test
##
## data: funcao
## BP = 3.1725, df = 7, p-value = 0.8686dwt(funcao)## lag Autocorrelation D-W Statistic p-value
## 1 0.01500952 1.943119 0.866
## Alternative hypothesis: rho != 0dwtest(funcao)##
## Durbin-Watson test
##
## data: funcao
## DW = 1.9431, p-value = 0.4248
## alternative hypothesis: true autocorrelation is greater than 0bgtest(funcao)##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: funcao
## LM test = 0.0043795, df = 1, p-value = 0.9472predict(funcao, newdata = data.frame(x1=1,x2=0,x3=1),
se.fit = T, interval = "pred")## $fit
## fit lwr upr
## 1 7.745263 7.374278 8.116249
##
## $se.fit
## [1] 0.0902356
##
## $df
## [1] 9
##
## $residual.scale
## [1] 0.1369392obj=function(x){a<-predict(funcao,
newdata = data.frame(x1=x[1],x2=x[2],x3=x[3]))
return(-a)}
x0<-c(0,0,0)
optim(par = x0,fn=obj,lower = -1,
upper = 1,method = c("L-BFGS-B"))## $par
## [1] 1 1 1
##
## $value
## [1] -9.357237
##
## $counts
## function gradient
## 3 3
##
## $convergence
## [1] 0
##
## $message
## [1] "CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL"predict(funcao, newdata = data.frame(x1=1,x2=1,x3=1),
se.fit = T, interval = "pred")## $fit
## fit lwr upr
## 1 9.357237 8.912707 9.801767
##
## $se.fit
## [1] 0.140935
##
## $df
## [1] 9
##
## $residual.scale
## [1] 0.1369392library(Rsolnp)
x0<-rep(0,3) # ponto de partida
otim4<-solnp(pars=x0,
fun = obj,
UB=rep(1,3),
LB=rep(-1,3))##
## Iter: 1 fn: -9.3572 Pars: 1.00000 1.00000 1.00000
## Iter: 2 fn: -9.3572 Pars: 1.00000 1.00000 1.00000
## solnp--> Completed in 2 iterationsotim4## $pars
## [1] 1.0000000 0.9999999 1.0000000
##
## $convergence
## [1] 0
##
## $values
## [1] -3.445789 -9.357237 -9.357237
##
## $lagrange
## [1] 0
##
## $hessian
## [,1] [,2] [,3]
## [1,] 1.914838 -1.0586114 1.382031
## [2,] -1.058611 0.9376314 -1.451680
## [3,] 1.382031 -1.4516802 2.967549
##
## $ineqx0
## NULL
##
## $nfuneval
## [1] 117
##
## $outer.iter
## [1] 2
##
## $elapsed
## Time difference of 0.06175089 secs
##
## $vscale
## [1] 1 1 1 1otim4$pars # Valores dos fatores otimizados## [1] 1.0000000 0.9999999 1.0000000n= 50000
z<-c()
for (j in 1:n) {
x1= 1
x2=1
x3= 1
b0<- runif(1, 3.3391398, 3.5524392 )
b1<-runif(1, 2.4376077, 2.6273923 )
b2<-runif(1, 0.5426077, 0.7323923 )
b3<-runif(1, 1.3126077 , 1.5023923 )
b22<-runif(1, -0.5211450 , -0.2599076 )
b33<-runif(1, 0.2288550 , 0.4900924 )
b12<-runif(1, 0.4808021 , 0.7491979 )
b23<-runif(1, 0.6158021 , 0.8841979 )
y= b0+b1*x1+b2*x2+b3*x3+b22*x2^2+b33*x3^2+b12*x1*x2+b23*x2*x3
#a<- predict(funcao, newdata = data.frame(x1= runif(1,-1,1), x2=runif(1,-1,1),
# x3=runif(1,-1,1)))
z[j]<-y
}
hist(z,freq = F, breaks = 200,xlab = "Distância",
ylab = "Densidade", col="purple", lwd=3,pch=12)
library(fBasics)##
## Attaching package: 'fBasics'## The following object is masked from 'package:car':
##
## densityPlotbasicStats(z)## z
## nobs 50000.000000
## NAs 0.000000
## Minimum 8.633879
## Maximum 10.036809
## 1. Quartile 9.226362
## 3. Quartile 9.487226
## Mean 9.356578
## Median 9.356332
## Sum 467828.897203
## SE Mean 0.000850
## LCL Mean 9.354912
## UCL Mean 9.358244
## Variance 0.036134
## Stdev 0.190089
## Skewness 0.002228
## Kurtosis -0.158210boxplot(z,horizontal = T,col="blue",lwd=4, main="Box-Plot",
xlab="Distância")
library(qcc)## Package 'qcc' version 2.7## Type 'citation("qcc")' for citing this R package in publications.library(SixSigma)
a<-qcc(z,type = c("xbar.one"))
process.capability(a,spec.limits = c(8.8,9.85),
target = 9.36)
##
## Process Capability Analysis
##
## Call:
## process.capability(object = a, spec.limits = c(8.8, 9.85), target = 9.36)
##
## Number of obs = 50000 Target = 9.36
## Center = 9.357 LSL = 8.8
## StdDev = 0.1912 USL = 9.85
##
## Capability indices:
##
## Value 2.5% 97.5%
## Cp 0.9150 0.9094 0.9207
## Cp_l 0.9701 0.9645 0.9757
## Cp_u 0.8600 0.8549 0.8651
## Cp_k 0.8600 0.8539 0.8661
## Cpm 0.9149 0.9092 0.9206
##
## Exp<LSL 0.18% Obs<LSL 0.11%
## Exp>USL 0.49% Obs>USL 0.41%