######################################## ######################################## ######################################## # Dados da potencia de turbina de aviões ######################################## ######################################## ######################################## # limpar todas as variáveis rm(list = ls(all.names = TRUE)) #library(xtable) #library(TeachingDemos) library(plotrix) library(plyr) library(xtable) library(e1071) library(MASS) # # Carregando os dados m.dados <- scan("C:\\Users\\cnaber\\Trabalho\\windows\\Unicamp\\Disciplinas\\1_semestre_2019\\ME 613\\Dados\\Turbina.prn",what=list(0,0)) # Gerando e/ou nomeando as variáveis veloc<- tempo <- cbind(m.dados[[2]]) grupo<- cbind(m.dados[[1]]) grupofac <- factor(grupo,c("1","2","3","4","5")) # Imprimindo os dados em forma de tabela mdados<-cbind(veloc,grupo) xtable(mdados) #################### #################### # Análise descritiva #################### #################### dataturb<-data.frame(veloc,grupofac) #medidas resumo cturb<-ddply(dataturb,.(grupofac),summarise,media=mean(veloc),dp=sqrt(var(veloc)),vari=var(veloc),cv=100*((sqrt(var(veloc))/mean(veloc))),ca=skewness(veloc),minimo=min(veloc),maximo=max(veloc)) xtable(cturb) # boxplot par(mfrow=c(1,1)) boxplot(veloc~grupofac,names=c("1","2","3","4","5"),xlab="tipo de turbina",ylab="tempo de vida (em milhões de ciclos)",cex=1.2,cex.axis=1.2,cex.lab=1.2) ##################### ##################### # Análise Inferencial ##################### ##################### ################ ################ # Modelo normal fit.model <-result<- lm(veloc~grupofac) summary(fit.model) n <- length(veloc) model.matrix(result) p<-ncol(model.matrix(result)) summary(result) sumres<-summary(result) m.X<-(model.matrix(result)) v.beta <- result$coefficients ep.beta<- sqrt(diag(vcov(result))) library(xtable) quantilt<-qt(0.975,df=n-p) xtable(cbind(sumres$coefficients[,1:2],v.beta-quantilt*ep.beta,v.beta+quantilt*ep.beta,sumres$coefficients[,3],sumres$coefficients[,4]),digits=3) #source("C:\\Users\\cnaber\\Trabalho\\windows\\Unicamp\\Disciplinas\\1_semestre_2019\\ME 613\\Programas\\diag_norm.r") source("C:\\Users\\cnaber\\Trabalho\\windows\\Unicamp\\Disciplinas\\1_semestre_2019\\ME 613\\Programas\\diag2_norm.r") source("C:\\Users\\cnaber\\Trabalho\\windows\\Unicamp\\Disciplinas\\1_semestre_2019\\ME 613\\Programas\\envel_norm.r") diag2norm(fit.model) par(mfrow=c(1,1)) envelnorm(fit.model) #source("C:\\Users\\cnaber\\Trabalho\\windows\\Unicamp\\Disciplinas\\Pos\\2_semestre_2014\\Regressao\\Programas\\alainflu_norm.txt") # R2 summary(result)$r.squared summary(result)$adj.r.squared ################################### ################################### grupofacr <- revalue(grupofac,c("1"="1","2"="2","3"="1","4"="1","5"="5")) fit.model <-result<- lm(veloc~grupofacr) summary(fit.model) n <- length(veloc) model.matrix(result) p<-ncol(model.matrix(result)) summary(result) sumres<-summary(result) m.X<-(model.matrix(result)) v.beta <- result$coefficients ep.beta<- sqrt(diag(vcov(result))) library(xtable) quantilt<-qt(0.975,df=n-p) xtable(cbind(sumres$coefficients[,1:2],v.beta-quantilt*ep.beta,v.beta+quantilt*ep.beta,sumres$coefficients[,3],sumres$coefficients[,4]),digits=3) #source("C:\\Users\\cnaber\\Trabalho\\windows\\Unicamp\\Disciplinas\\1_semestre_2019\\ME 613\\Programas\\diag_norm.r") source("C:\\Users\\cnaber\\Trabalho\\windows\\Unicamp\\Disciplinas\\1_semestre_2019\\ME 613\\Programas\\diag2_norm.r") source("C:\\Users\\cnaber\\Trabalho\\windows\\Unicamp\\Disciplinas\\1_semestre_2019\\ME 613\\Programas\\envel_norm.r") diag2norm(fit.model) par(mfrow=c(1,1)) envelnorm(fit.model) #source("C:\\Users\\cnaber\\Trabalho\\windows\\Unicamp\\Disciplinas\\Pos\\2_semestre_2014\\Regressao\\Programas\\alainflu_norm.txt") # R2 summary(result)$r.squared summary(result)$adj.r.squared # Médias preditas pelo modelo resmupred <- predict(fit.model,type = c("response"),se.fit=TRUE) mupred3 <- unique(resmupred$fit) mupred3 <- c(mupred3[1],mupred3[2],mupred3[1],mupred3[1],mupred3[3]) epmupred3 <- unique(resmupred$se.fit) epmupred3 <- c(epmupred3[1],epmupred3[2],epmupred3[1],epmupred3[1],epmupred3[3]) ez <-qnorm(0.975) plotCI(c(1,2,3,4,5),mupred3,uiw=epmupred3*quantilt,liw=epmupred3*quantilt,pch=19,lwd=2,col=1,xlab="tipo de turbina",ylab="tempo de vida médio",cex=1.2,cex.lab=1.2,cex.axis=1.2,xlim=c(1,5.5)) lines(c(0.9,1.9,2.9,3.9,4.9),cturb$media,pch=15,lwd=2,col=3,type="p",cex=1.2) legend(1.5,14,c("observada","modelo reduzido"),col=c(3,1),pch=c(15,19),bty="n",cex=1.2) xtable(cbind(mupred3,epmupred3,mupred3-quantilt*epmupred3,mupred3+quantilt*epmupred3))