Para empezar, si queremos resolver el problema de la catapulta tendremos que usar las siguientes librerías:

library(readxl)
library(SixSigma)
library(sjPlot)
library(ggplot2)
library(tidyverse)
library(DT)
library(qcc)

Cargamos los datos, en nuestra práctica hemos utilizado cuatro variables: Nº de tiro(tiro), Bungee Position(BP), Release Angle(RA) y Distancia(d).

Para iniciar el análisis pasamos BP a factor.

datos <- read_excel("DOE.xlsx")
str(datos)
## tibble [20 x 4] (S3: tbl_df/tbl/data.frame)
##  $ Tiro: num [1:20] 1 2 3 4 5 1 2 3 4 5 ...
##  $ BP  : num [1:20] 150 150 150 150 150 150 150 150 150 150 ...
##  $ RA  : num [1:20] 135 135 135 135 135 160 160 160 160 160 ...
##  $ d   : num [1:20] 124 125 122 124 125 156 156 159 156 154 ...
datos$BP = factor(datos$BP)
summary(datos)
##       Tiro     BP           RA              d        
##  Min.   :1   150:10   Min.   :135.0   Min.   :122.0  
##  1st Qu.:2   200:10   1st Qu.:135.0   1st Qu.:134.8  
##  Median :3            Median :147.5   Median :147.5  
##  Mean   :3            Mean   :147.5   Mean   :150.0  
##  3rd Qu.:4            3rd Qu.:160.0   3rd Qu.:163.2  
##  Max.   :5            Max.   :160.0   Max.   :183.0

Ahora visualizamos los datos obtenidos.

ggplot(datos,aes(x=RA,y=d))+
  geom_point(aes(fill=BP,color=BP))+
  labs(caption="Fig1: Datos registrados de tiros con la catapulta virtual.")

A simple vista podemos ver que a mayor ángulo de tiro y BP mayor distancia recorrerá.

Ajustamos el modelo lineal para predecir la distancia.

fit = lm(d ~ RA*BP,data=datos)
summary(fit)
## 
## Call:
## lm(formula = d ~ RA * BP, data = datos)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##  -3.60  -0.60  -0.10   0.85   3.40 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -49.88000    6.80288  -7.332 1.69e-06 ***
## RA            1.28800    0.04596  28.027 5.00e-15 ***
## BP200       -22.68000    9.62073  -2.357 0.031471 *  
## RA:BP200      0.28800    0.06499   4.431 0.000419 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.817 on 16 degrees of freedom
## Multiple R-squared:  0.9938, Adjusted R-squared:  0.9926 
## F-statistic: 851.8 on 3 and 16 DF,  p-value: < 2.2e-16

Una vez tenemos el modelo ajustado procedemos con la predicción.

newdata<-expand.grid(BP=factor(levels(datos$BP)), RA=as.numeric(datos$RA))
pred.fit<-predict(fit,newdata,se.fit=TRUE, interval="prediction")  
prediccion<-cbind(newdata,pred=pred.fit$fit,sd=pred.fit$se.fit,df=pred.fit$df)

ggplot(prediccion,aes(x=RA,y=pred.fit,group=BP))+
  geom_point(aes(fill=BP,color=BP))+
  geom_errorbar(aes(ymin = pred.lwr, ymax = pred.upr,color=BP))+
  ylab("Distancia")+ggtitle("Predicción de la distancia")+
  labs(caption="Fig2: Ajuste ANOVA.")

Para hayar el ángulo de tiro resolvemos el modelo a continuación:

Y= RA + BP(factor) + RA:BP(factor)

Para predecir el ángulo de tiro para que con BP=150 y BP=200 alcance la distancia de 150+-5 despejamos las ecuaciones tal que así:

Para BP=150 -> Y=-49.88+1.29RA

Para BP=200 -> Y=(-49.88-22.68) + (1.29+0.29)RA -> Y=-72.56+1.58RA

tab_model(fit)
  d
Predictors Estimates CI p
(Intercept) -49.88 -64.30 – -35.46 <0.001
RA 1.29 1.19 – 1.39 <0.001
BP [200] -22.68 -43.08 – -2.28 0.031
RA * BP [200] 0.29 0.15 – 0.43 <0.001
Observations 20
R2 / R2 adjusted 0.994 / 0.993 
v1<-(150+49.88)/1.29;v1
## [1] 154.9457
v2<-(150+72.56)/1.58;v2
## [1] 140.8608
ggplot(prediccion,aes(x=RA,y=pred.fit,group=BP))+
  geom_point(aes(fill=BP,color=BP))+
  geom_errorbar(aes(ymin = pred.lwr, ymax = pred.upr,color=BP))+
  ylab("Distancia")+ggtitle("Predicción de la distancia")+
  geom_line(position=position_dodge(width=0.1))+
  geom_hline(yintercept=150,linetype=3)+
  geom_vline(xintercept=c(v1,v2),linetype=3)+
  labs(caption="Fig3: Ajuste ANOVA y líneas de interpolación.")

Creamos una tabla con todos los posibles resultados.

fitr<-lm(d ~ RA*BP,data=datos)

newdatar<-expand.grid(BP=factor(levels(datos$BP)), RA=as.numeric(seq(135,160,by=1)))
pred.fitr<-round(predict(fitr,newdatar,interval = "prediction"),2)
prediccionr<-cbind(newdatar,pred=pred.fitr)
datatable(prediccionr,filter="top",colnames=c("BP","RA","Distancia.pred", "IP.inf","IP.sup"),
          caption="Tabla1: Predicción de la distancia en función de la configuración BP, RA.")

Una vez tenemos los RA localizados para cada BP(BP=150->RA=155 y BP=200->RA=141) verificamos la capacidad del proceso para conseguir una distancia de 150.

datosc <- read_excel("C:/Users/USUARIO/Desktop/UMH/4 ANO/Primer cuatrimestre/Mejora de procesos/DOE2.xlsx")

cap <- qcc.groups(datosc$d,1:length(datosc$d))
q <- qcc(cap, type="xbar.one", nsigmas=3, plot=FALSE)

target= 150
lsl=target-5;usl=target+5
process.capability(q, spec.limits=c(lsl,usl),target=target)

## 
## Process Capability Analysis
## 
## Call:
## process.capability(object = q, spec.limits = c(lsl, usl), target = target)
## 
## Number of obs = 20           Target = 150
##        Center = 150.1           LSL = 145
##        StdDev = 1.446           USL = 155
## 
## Capability indices:
## 
##       Value    2.5%  97.5%
## Cp    1.152  0.7889  1.515
## Cp_l  1.164  0.8299  1.498
## Cp_u  1.141  0.8126  1.469
## Cp_k  1.141  0.7497  1.532
## Cpm   1.152  0.7972  1.505
## 
## Exp<LSL 0.024%    Obs<LSL 0%
## Exp>USL 0.031%    Obs>USL 0%

La capacidad del proceso para cumplir las especificaciones es parcialmente adecuado, Cp=1.15. Así que tendríamos que tener un control estricto.

cap.fun=function(datosc,lsl,usl){
  xbar=mean(datosc)
  s=sd(datosc)
  zu=(usl-xbar)/s
  zl=(lsl-xbar)/s
  pdfe=pnorm(zl)+1-pnorm(zu)
  zbench=qnorm(pdfe)
  zscore=min(c(zu,-zl))
  Pp=(usl-lsl)/(6*s)
  Ppk=min(c(zu/3,-zl/3 ))
  return(data.frame(zl=zl,zu=zu,pdfe=pdfe,zbench=zbench,zscore=zscore,Pp=Pp,Ppk=Ppk))}

cap.fun(datosc$d,lsl,usl)
##          zl       zu       pdfe    zbench   zscore        Pp       Ppk
## 1 -2.866991 2.810219 0.00454737 -2.608472 2.810219 0.9462015 0.9367395

Para concluir la efectividad (capacidad) del ataque al enemigo en términos del nivel sigma hemos calculado zbench que en valor absoluto le sumamos 1.5 para concluir que tenemos un valor de 4.1 sigma (2.6+1.5). Es decir, estimamos que el 99,5% de los tiros acertarán.

Respecto a los costes preveemos un gasto de 2000u.m. de las balas utilizadas para la obtención de los ángulos de tiro más 500u.m. para derribar el castilllo ya que una vez finalizado el análisis anterior, con los nuevos datos, todos los tiros acertarán dentro de los límites establecidos.