Principios fundamentales encontrados en la clase


A continuación se presenta la funcion x^3 y algunas maneras de darle argumentos a la función

cubos <- function(x) x**3
cubos(1:5)
## [1]   1   8  27  64 125
v1 = seq(from = 0, to = 20, by = 2); v1 #Seq: secuencias desde-hasta-paso
##  [1]  0  2  4  6  8 10 12 14 16 18 20
cubos(v1) 
##  [1]    0    8   64  216  512 1000 1728 2744 4096 5832 8000

Luego creamos una función para estandarizar variables a partir de los parametros x, x menos la media, y por ultimo los desvios sobre la desviación estandar. Para provar el optimo funcionamiento de la función se probo con una distribución uniforme y para una distribución normal

#####################
estandarizar <- function(x){
  med = mean(x)
  dv = x - med #Desvios
  z = dv/sd(x)
  par(mfrow=c(1,3))
  hist(x)
  hist(dv)
  hist(z)
  return(list(med.dat = med, med.des = mean(dv),
              var.dat = var(x), var.des = var(dv),
              med.z = mean(z), var.z = var(z),
              z = z))
}
set.seed(123)
##Generacion de datos aleatorios con una distribución normal
estandarizar(rnorm(100, 5, 1))

## $med.dat
## [1] 5.090406
## 
## $med.des
## [1] -2.087642e-16
## 
## $var.dat
## [1] 0.8332328
## 
## $var.des
## [1] 0.8332328
## 
## $med.z
## [1] -2.366325e-16
## 
## $var.z
## [1] 1
## 
## $z
##   [1] -0.71304802 -0.35120270  1.60854170 -0.02179795  0.04259548
##   [6]  1.77983218  0.40589817 -1.48492941 -0.85149566 -0.58726835
##  [11]  1.24195461  0.29513939  0.34000892  0.02221347 -0.70797086
##  [16]  1.85854263  0.44636008 -2.25349176  0.66930255 -0.61698896
##  [21] -1.26885349 -0.33783464 -1.22304003 -0.89754917 -0.78377819
##  [26] -1.94683206  0.81876439  0.06898128 -1.34588242  1.27452758
##  [31]  0.36815564 -0.42229479  0.88157948  0.86296437  0.80101058
##  [36]  0.65537241  0.50778230 -0.16686565 -0.43422620 -0.51585092
##  [41] -0.86009994 -0.32681639 -1.48529653  2.27707482  1.22429519
##  [46] -1.32941869 -0.54040552 -0.61026684  0.75541982 -0.19037243
##  [51]  0.17847258 -0.13031397 -0.14600575  1.40027842 -0.34637532
##  [56]  1.56226982 -1.79571669  0.54141021  0.03664303  0.13752572
##  [61]  0.31685861 -0.64934164 -0.46407310 -1.21490140 -1.27319995
##  [66]  0.23347834  0.39197814 -0.04097396  0.91131364  2.14685001
##  [71] -0.63697081 -2.62876100  1.00275711 -0.87597805 -0.85276181
##  [76]  1.02448422 -0.41101270 -1.43635058  0.09957931 -0.25119772
##  [81] -0.09272595  0.32303830 -0.50510289  0.60688103 -0.34058618
##  [86]  0.26443017  1.10255872  0.37770550 -0.45610238  1.15949090
##  [91]  0.98935390  0.50173432  0.16249260 -0.78691881  1.39156928
##  [96] -0.75663177  2.29720706  1.57995139 -0.35725306 -1.22349625
##Generacion de datos aleatorios con una distribucion uniforme  
estandarizar(runif(100,2,6))

## $med.dat
## [1] 3.944207
## 
## $med.des
## [1] -7.103042e-17
## 
## $var.dat
## [1] 1.386103
## 
## $var.des
## [1] 1.386103
## 
## $med.z
## [1] -5.748874e-17
## 
## $var.z
## [1] 1
## 
## $z
##   [1] -0.84029420  1.61826399  0.39178189  0.09845350 -0.28361946
##   [6]  1.33928538 -0.41436110 -0.67207186 -1.07160013 -1.06641378
##  [11] -0.01362095 -0.79191722 -0.91664070  0.63983716 -1.48943276
##  [16]  0.72979231 -0.45582182 -0.26197502  1.13782862  1.47046642
##  [21] -0.69147494  1.61400301  0.82336460  0.68060298 -1.47183253
##  [26] -0.30860214 -0.02788111  0.25210142  0.72098766  1.45968333
##  [31]  0.44949049 -0.19579977  0.19035862 -1.45268904 -0.76510423
##  [36] -0.30203874 -0.97952900  1.17512063 -1.13193336  1.07826062
##  [41]  0.20648254  0.59886732 -1.06802166  0.49944809 -0.59178685
##  [46]  0.81031784 -0.29596441  1.64203806  1.63538557  0.81761637
##  [51] -0.77747158 -0.89784172  0.36351426 -0.74246119  0.15295202
##  [56]  1.01667419 -1.08038077 -0.27741613 -0.04918052  1.29804038
##  [61]  1.49374150  1.34516655  0.63919318  1.57684155  0.10326156
##  [66]  0.30736468 -0.50867853 -0.47132813 -1.58333796  0.05694706
##  [71]  1.30801756 -1.62996390 -1.40655533 -1.09345982  0.96585544
##  [76]  0.84643335  1.65059718 -0.06652113 -1.39864798  0.55300242
##  [81]  0.92596547 -1.18555012 -0.30396636 -0.88697851 -1.45445550
##  [86] -0.30631712 -1.43077566 -0.88391699 -1.46576738  0.62592653
##  [91] -0.63978691 -1.30916723 -1.40707525  1.33994459  0.91120063
##  [96]  1.12306491  1.68547184 -1.29938895 -1.31487422  1.06267643

Para explicar un nuevo concepto fue necesario la creación de dos nueva variables N y P, donde N simboliza el nitrogeno y la P simboliza el Potasio

set.seed(123)
N = sort(rnorm(50, 4, 0.5))
P = sort(rnorm(50, 0.5, 0.01))

Con las nuevas variables creadas P y N creamos una regresión lineal para ver como se comportan las variables

#Regresion lineal simple
mod1 = lm(P~N)
mod1
## 
## Call:
## lm(formula = P ~ N)
## 
## Coefficients:
## (Intercept)            N  
##     0.42369      0.01936
summary(mod1)
## 
## Call:
## lm(formula = P ~ N)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0051859 -0.0008680  0.0003719  0.0008153  0.0021135 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.4236910  0.0016234  260.99   <2e-16 ***
## N           0.0193600  0.0004015   48.22   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.001301 on 48 degrees of freedom
## Multiple R-squared:  0.9798, Adjusted R-squared:  0.9794 
## F-statistic:  2325 on 1 and 48 DF,  p-value: < 2.2e-16

Ahora creamos una grafica para observar graficamente las relaciones existentes entre las variables, ademas de decorar el grafico con lineas, simbolos y texto

plot(N,P, pch = 19, col = 'blue', cex= 1.5,
     ylim= c(0, 5), xlim = c(0,5))
grid(10)
abline(a = 0.423, b = 0.0193, col ='red')
abline(h = 0.423, lty = 2)
abline(v = 0, lty = 2)
text(1,0.423, expression(theta), pos=3)
angulo = atan(0.0193)*180/pi
angulo
## [1] 1.105671
text(1.6,0.423, paste('= ', round(angulo,2)), pos=3)

Para comprobar uno de los principios que se encuentran al inicio de este documento estandarizamos las variables P y N que fueron creadas hace poco (Los analisis estadisticos son invariantes a la estandarización)

#Estandarizando valores N y P
N2 = estandarizar(N)$z

mean(N2); sd(N2)
## [1] 7.006548e-17
## [1] 1
P2 = estandarizar(P)$z

mean(P2); sd(P2)
## [1] -5.144435e-15
## [1] 1

Para demostrar que los analisis estadisticos son invariantes a la estandarización realizamos una nueva regresión lineal para las variables estandarizadas y luego se crea una grafica

#Regresion lineal simple para la variables estandarizadas y una grafico mostrando la interacción entre ellas 
mod2 = lm(P2~N2)
summary(mod2)
## 
## Call:
## lm(formula = P2 ~ N2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.57274 -0.09587  0.04108  0.09004  0.23342 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -5.835e-15  2.032e-02    0.00        1    
## N2           9.898e-01  2.053e-02   48.22   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1437 on 48 degrees of freedom
## Multiple R-squared:  0.9798, Adjusted R-squared:  0.9794 
## F-statistic:  2325 on 1 and 48 DF,  p-value: < 2.2e-16
angulo2 = atan(9.898e-1)*180/pi;angulo2
## [1] 44.7063
plot(P2~N2, col = 'blue', pch = 19)

Se hizo un analisis de correlación para las 4 variables, las dos estandarizadas y las no estandarizadas (o las variables originales) y se muestran los resultados en donde se observa que el resultado es el mismo que el de las variables sisn estadarizar, por lo que probamos que los estadisticos son invariantes a la estandarización 

#Calculo de correlaciones
cor(P,N)
## [1] 0.9898341
cor(P2,N2)
## [1] 0.9898341

Creamos una nueva variable a la cual denominamos fertilización, y la agrupamos junto con las otras 4 variables (potasio y sodio estandarizados y las no estandarizadas) 

fert = gl(2, 25, 50, labels = c('control','quimico'))
tbl1 = data.frame(P,N,P2,N2,fert);tbl1
##            P        N          P2          N2    fert
## 1  0.4769083 3.016691 -2.71200476 -2.16123294 control
## 2  0.4845125 3.156653 -1.87218101 -1.85889692 control
## 3  0.4877928 3.367302 -1.50989032 -1.40386871 control
## 4  0.4892821 3.367469 -1.34541195 -1.40350677 control
## 5  0.4897358 3.430932 -1.29530375 -1.26642024 control
## 6  0.4898142 3.438446 -1.28663895 -1.25018864 control
## 7  0.4929080 3.466088 -0.94495739 -1.19047736 control
## 8  0.4931199 3.486998 -0.92155222 -1.14530984 control
## 9  0.4937209 3.635554 -0.85517337 -0.82440817 control
## 10 0.4939974 3.652647 -0.82463985 -0.78748695 control
## 11 0.4949768 3.656574 -0.71647658 -0.77900398 control
## 12 0.4950897 3.687480 -0.70400508 -0.71224128 control
## 13 0.4962934 3.719762 -0.57106399 -0.64250835 control
## 14 0.4966679 3.722079 -0.52970029 -0.63750278 control
## 15 0.4967407 3.763604 -0.52166471 -0.54780365 control
## 16 0.4971523 3.766672 -0.47620808 -0.54117631 control
## 17 0.4976430 3.777169 -0.42201094 -0.51850209 control
## 18 0.4977423 3.798558 -0.41104467 -0.47229999 control
## 19 0.4977951 3.809764 -0.40520842 -0.44809159 control
## 20 0.4986111 3.847019 -0.31509251 -0.36761772 control
## 21 0.4995713 3.852464 -0.20904446 -0.35585454 control
## 22 0.4997145 3.884911 -0.19322498 -0.28576479 control
## 23 0.5000576 3.891013 -0.15533107 -0.27258521 control
## 24 0.5005300 3.896041 -0.10315791 -0.26172231 control
## 25 0.5012385 3.958315 -0.02490926 -0.12720211 control
## 26 0.5018130 3.969044  0.03853920 -0.10402677 quimico
## 27 0.5021594 4.035254  0.07679443  0.03899559 quimico
## 28 0.5023873 4.055341  0.10196449  0.08238648 quimico
## 29 0.5025332 4.064644  0.11807452  0.10248111 quimico
## 30 0.5030353 4.076687  0.17352792  0.12849490 quimico
## 31 0.5033178 4.179907  0.20473164  0.35146434 quimico
## 32 0.5037964 4.200386  0.25758676  0.39570124 quimico
## 33 0.5038528 4.213232  0.26381674  0.42345111 quimico
## 34 0.5043518 4.230458  0.31892883  0.46066150 quimico
## 35 0.5044821 4.248925  0.33331762  0.50055293 quimico
## 36 0.5054840 4.276959  0.44396701  0.56110914 quimico
## 37 0.5058461 4.344320  0.48396580  0.70661834 quimico
## 38 0.5064438 4.350678  0.54996942  0.72035206 quimico
## 39 0.5092227 4.389983  0.85687955  0.80525515 quimico
## 40 0.5099350 4.410791  0.93555492  0.85020310 quimico
## 41 0.5100574 4.418894  0.94906721  0.86770660 quimico
## 42 0.5102557 4.439067  0.97097113  0.91128339 quimico
## 43 0.5109684 4.447563  1.04968101  0.92963605 quimico
## 44 0.5114881 4.603981  1.10707652  1.26751971 quimico
## 45 0.5136065 4.612041  1.34104359  1.28493014 quimico
## 46 0.5136860 4.626907  1.34982360  1.31704386 quimico
## 47 0.5151647 4.779354  1.51313331  1.64634862 quimico
## 48 0.5153261 4.857532  1.53095878  1.81522403 quimico
## 49 0.5205008 4.893457  2.10247091  1.89282472 quimico
## 50 0.5218733 5.084478  2.25405159  2.30545591 quimico

Ahora, realizaremos la prueba T-student, a partir de una función creada por nosotros simulando computacionalmente la ecuación de la curva - T, y utilizaremos los datos contenidos en la tabla que se creo en el paso anterior, es decir que probaremos la funcion calculando el valor para los datos de potasio y sodio estandarizados y las no estandarizadas

#Prueba t-student####
#Ho: mu(control(N)) == mu(quimico(N))
#Ha: mu(control(N)) != mu(quimico(N))

#Crear funcion para prueba t####
mi.pt <- function(x1,x2, conf.lev){
  n1 = length(x1)
  n2 = length(x2)
  s1 = sd(x1)
  s2 = sd(x2)
  tc = abs((mean(x1)-mean(x2))/
             sqrt(((((n1-1)*s1^2)+((n2-1)*s2^2))/
                     (n1+n2-2))*((1/n1)+(1/n2))))
  p_v = pt(q = tc, df = (n1+n2-2), lower.tail = FALSE)
  desi = ifelse(p_v < (1-conf.lev)/2, 'Rechazo Ho', 'NO rechazo Ho')
  return(desi)
}

mi.pt(tbl1$P[1:25],tbl1$P[26:50], 0.95)
## [1] "Rechazo Ho"
mi.pt(tbl1$P2[1:25],tbl1$P2[26:50], 0.95)
## [1] "Rechazo Ho"

Ahora compareramos el resultado de la función creada por nosotros con el resultado de la funcion que viene por defecto en RStudio, el problema es que esta función no mostraba un p-valor por lo que tubiemos que hacer un agrreglo en el programa para verlo

t.test(tbl1$P[1:25],tbl1$P[26:50])
## 
##  Welch Two Sample t-test
## 
## data:  tbl1$P[1:25] and tbl1$P[26:50]
## t = -8.6607, df = 47.974, p-value = 2.262e-11
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.01724883 -0.01074889
## sample estimates:
## mean of x mean of y 
## 0.4944647 0.5084635
tbl2 = split(tbl1$P, tbl1$fert)
PC = tbl2$control
PQ = tbl2$quimico
mi.pt(PC,PQ, 0.95)
## [1] "Rechazo Ho"

Esto lo hicimos como forma de aclaración de conceptos para observar como funcionan los datos que se encierran dentro de los parentesis 

(x = 1:10)
##  [1]  1  2  3  4  5  6  7  8  9 10
y = 1:10

Por último realizamos esta nueva función, esta tiene relación con la lue de los gases 

ley.gases <- function(n, R = 0.082, Tp, v, gra = F){ #Argumentos fijos: logicos, numericos, caracteres
  p = n*R*Tp/v
  if(gra) plot3d(Tp, v, p, col = rainbow(100))
  return(p)
}

Evaluaremos la correcta operatividad de la función que acabamos de construir para dos condiciones diferentes ademas se hizo una gráfica para el caso 1 

###Para un único valor 
p1 = ley.gases(n = 1, Tp = 273, v = seq(1,20))
plot(seq(1,20), p1, type = 'l')

###Para un intervalo de valores definido (Se conoce la amplitud y la cantidad de los datos)
p2 = ley.gases(n = 1,
               Tp = rep(seq(273,303,5), each = 20),
               v = rep(seq(1,20), 7), gra = )
#Tamaño de muestra####
set.seed(123)
var = replicate(200, sample(c(0,1),
                           size = 100,
                           replace = T,
                           prob = c(0.2,0.8)))
colSums(var)
##   [1] 82 81 81 81 81 73 78 79 87 79 85 82 78 84 75 85 78 72 77 87 81 79 77
##  [24] 82 84 80 82 78 79 81 77 80 75 77 82 78 80 72 80 74 80 82 84 80 85 85
##  [47] 81 79 89 82 73 79 83 83 79 79 86 83 78 71 82 78 87 85 79 80 79 81 86
##  [70] 78 77 80 77 81 81 85 86 84 75 76 78 81 82 90 88 80 78 82 84 86 82 81
##  [93] 83 80 77 81 87 79 77 74 80 83 86 70 84 84 78 87 76 77 79 81 80 76 78
## [116] 80 77 78 81 81 77 76 84 74 76 81 85 76 77 80 76 77 79 79 79 88 77 73
## [139] 78 82 73 80 83 80 79 83 80 81 80 77 83 81 83 79 85 82 80 78 83 80 84
## [162] 80 76 84 76 79 73 79 80 74 86 76 75 89 79 80 83 77 77 83 80 82 82 82
## [185] 82 84 86 87 74 84 77 85 82 81 84 84 82 79 86 82
hist(colSums(var))

hist(round(runif(200,0,1)))

var2 = replicate(2000,round(runif(200,0,1)))
pg2 = (colSums(var2)/200)*100

hist(pg2)

TM <- function(N,p,e,Z){
  (N*p*(1-p))/(((N-1)*(e/Z))^2+(p*(1-p)))
}