Distribuciones de probabilidad usadas.

Normal

Genera datos en forma de campana

En r > \(norm\), se antepone la letra \(r\) para indicar \(random\), es decir, genera datos aleatorios normales.

Uniforme

Genera datos en forma horizontal (aproximadamente), frecuencia tiende a ser la misma. Ningún valor predomina sobre el otro

En r > \(unif\), se antepone la letra \(r\) para indicar \(random\), es decir, genera datos aleatorios uniformes.

Estandarización (Normalizar)

set.seed(123)
# aceite del cultivo de limonaria en cc de aceite
acc = rnorm(120, 5, 0.25)
hist(acc, col = 'orange',ylim = c(0, 30), main = 'Distribución del contenido de aceite (CC)', 
     xlab = 'Contenido de aceite', breaks = 10)
rug(acc)

# viscosidad (pascal segundo)
vis = runif(120, 2, 3)
hist(vis, col = 'orange',ylim = c(0, 30), main = 'Distribución de la viscosidad del aceite (Pascal/Seg)', 
     xlab = 'Viscosidad de aceite', breaks = 10)
rug(vis)

# grafico dispersión

plot(acc, vis, pch = 19, cex =0.8, main = 'Contenido de aceite vs viscosidad', ylab = 'Viscosidad (Pascal/Seg)', 
     xlab = 'Contenido de aceite (CC)')

Los graficos que presentan unidades diferentes en los ejes son muy engañosos.

Cambian los valores de los ejes, pero la dispersión se mantiene.

Puntuación Z de estandarización

\[Z = \frac{x-\mu_x}{\sigma_x}\]

# funcion para estandarizar
estand = function(x){
  media = mean(x)
  desv = sd(x)
  z = (x - media)/desv
  return(z)
}

# estandarizando acc
acc_z = estand(acc)

# estandarizando vis
vis_z = estand(vis)


par(mfrow = c(2, 2))
hist(acc, main = 'Acc no estandarizada')
hist(acc_z, main = 'Acc estandarizada')
hist(vis, main = 'vis no estandarizada')
hist(vis_z, main = 'vis estandarizada')

Estandarizar > hace que la media de la variable sea 0

(med_acc = mean(acc)) # media sin estandarizar
## [1] 5.00386
(desv_acc = sd(acc)) # desviación sin estandarizar
## [1] 0.223607
(med_acc_z = mean(acc_z)) # media estandarizada = 0
## [1] -4.953991e-16
(desv_acc_z = sd(acc_z)) # desviación estandarizada = 1
## [1] 1
(med_vis = mean(vis)) # media sin estandarizar
## [1] 2.494534
(desv_vis = sd(vis)) # desviación sin estandarizar
## [1] 0.298678
(med_vis_z = mean(vis_z)) # media estandarizada = 0
## [1] 6.907055e-16
(desv_vis_z = sd(vis_z)) # desviación estandarizada = 1
## [1] 1
par(mfrow = c(1, 2))
plot(acc_z, vis_z, pch = 19, cex = 0.8, main = 'Estandarizado')
points(x = mean(acc_z), y = mean(vis_z), col = 'red', pch = 19)
plot(acc, vis, pch = 19, cex = 0.8, main = 'No Estandarizado')
points(x = mean(acc), y = mean(vis), col = 'blue', pch = 19)

Correlación de Pearson entre aceite y viscosidad.

Es una medida de la asociación lineal entre un par de variables, cuando la correlación es 1 es una correlación peefecta, esto indica que una variable esta relacinada estadisticamente con la otra. Oscila entre -1 y 1.

Datos originales

cr = cor(acc, vis, method = 'pearson')
cr
## [1] -0.01905532

Datos estandarizados

cr_z = cor(acc_z, vis_z, method = 'pearson')
cr_z
## [1] -0.01905532
# biomasa de tuberculos
set.seed(1234)
biom = sort(rnorm(48, 3, 0.25))
# tiempo de pesado de tuberculos (segundos)
tiempo = sort(rnorm(48, 30, 2))
plot(y = biom, x =tiempo)

corr = cor(biom, tiempo)
corr
## [1] 0.9806255

Preguntas para entender el fenómeno

  1. ¿Tiene sentido en estudiar la relación la biomasa y lo que se tarde en pesar lo tubérculos?

  2. Según Descartes, en el eje Y debe ir la variable dependiente y en el eje X la variable independiente.

# biomasa de niños al nacer
set.seed(1234)
peso = sort(rnorm(48, 3, 0.25))
# Numero de Cigüeñas contadas
num = sort(floor(rnorm(48, 30, 2)))
plot(x = peso, y = num, ylim = c(0, 36), xlim = c(0,36))

corr = cor(peso, num)
corr
## [1] 0.9681666
set.seed(1243)
# N en el tejdio en ppm
N = sort(rnorm(120, 10, 0.8))
# P en el tejido en ppm
P = sort(rnorm(120, 0.1, 0.02))
plot(N, P, pch = 19, cex = 0.6)

Estandarizando

N_z = estand(N)
P_z = estand(P)
par(mfrow = c(1, 2))
plot(N, P, pch = 19, cex = 0.6, main = 'Sin estandarizar')
plot(N_z, P_z, pch = 19, cex = 0.6, xlim = c(-3, 3), 
     ylim = c(-3, 3), main = 'Estandarizado')

plot(N, P, pch = 19, cex = 0.6, main = 'Sin estandarizar', 
     xlim = c(0, 12), ylim = c(0, 12))

Calculando el angulo que forma una linea recta con respecto a su horizontal

Regresión lineal simple

mod = lm(P~N) # regresión lineal
plot(N, P, pch = 19, cex = 0.6)
abline(mod) # recta que mejor se ajusta a los datos

##### Ecuación de la recta \[y = mx+b\]

library(pander)
Summary del modelo
  Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.1638 0.003663 -44.72 7.69e-76
N 0.02655 0.0003637 72.99 4.798e-100
Fitting linear model: P ~ N
Observations Residual Std. Error \(R^2\) Adjusted \(R^2\)
120 0.003158 0.9783 0.9782 
# coeficientes de la recta
b = -0.1638
m = 0.0265
# m es igual a la tangente teta
# tang 0 = 0.0265
# 0 = arctang(0.0265)

\[Pendiente~es~igual~a~tan(\theta)\] \[tan(\theta) = Pendiente\] \[\theta = arctan(pendiente)\]

ang_r = atan(m) # angulo de la recta en radianes
#ang_r
ang_g = ang_r*180/pi # angulo en grados
ang_g
## [1] 1.517983

Para montar en rpubs

  1. Registrar/Simular datos de temperatura y datos de humedad relativa,
  2. Grafico sin estandarizar y grafico estandarizado.
  3. Determinar si el grafico es realista o engañoso.
  4. Graficar los demas modelos de creciiento de Growthmodels
library(growthmodels)
growth <- blumberg(0:100, 10, 2, 0.5)
par(bg = 'white', fg = 'darkred')
plot(growth, pch = 19, cex = 0.5, main = 'Modelo de crecimiento Blumberg', xlab = 'Time', col = 'darkblue')
grid(nx = 10, ny = 10, lwd = 1,col = 'black')

Correlación de Pearson

cor_p = cor(growth, 0:100, method = 'pearson')
cor_p
## [1] 0.7666426

Correlación de Spearman

cor_s = cor(growth, 0:100, method = 'spearman')
cor_s
## [1] 1

TALLER

1. Registrar/Simular datos de temperatura y datos de humedad relativa.

library(readxl)
Acacias_Meta <- read_excel("Acacias_Meta_2011_2019_JOINT.xlsx", 
    sheet = "Acacias_Meta_2011_2019_JOINT", 
    col_types = c("numeric", "text", "text", 
        "text", "numeric", "numeric", "text", 
        "text", "numeric", "text"))

tmedacacias<- Acacias_Meta$Tmed
HRacacias<- Acacias_Meta$RHUM

2.Grafico sin estandarizar y grafico estandarizado.

estand1 = function(x){
media = mean(x)
desv = sd(x)
z = (x - media)/desv
return(z)
}
testand<- estand1(tmedacacias)
HRestand<- estand1(HRacacias)

par(mfrow = c(2, 2))
hist(tmedacacias, main = "Distribución de T \ no Estandarizada", xlab = "Temperatura Media")
hist(HRacacias, main = "Distribucion de la HR\ no Estandarizada", xlab = "Humedad Relativa")
hist(testand, main = "Distribución de T\ Estandarizada", xlab = "Temperatura Media")
hist(HRestand, main = "Distribucion de la HR\ Estandarizada", xlab = "Humedad Relativa")

par(mfrow = c(1,2))
plot(tmedacacias,HRacacias, main = "Sin estandarizar")
points(x = mean(tmedacacias), y = mean(HRacacias), col = "red", pch = 19)
plot(testand,HRestand, main = "Estandarizado")
points(x = mean(testand), y = mean(HRestand), col = "blue", pch = 19)

plot(tmedacacias,HRacacias, xlim=c(0,35), ylim=c(0,100))

### 3.Determinar si el grafico es realista o engañoso.

Analizando las graficas y viendo los resultados de la correlacion de pearson y spearman estandarizando y sin estandarizar se puede concluir que la grafica es realista y la estandarizacion no modifica los resultados.

(cr_pse<-cor(tmedacacias,HRacacias, method = 'pearson'))
## [1] -0.6034923
(cr_pe<-cor(testand,HRestand, method = 'pearson'))
## [1] -0.6034923
(cr_sse<-cor(tmedacacias,HRacacias, method = 'spearman'))
## [1] -0.6447906
(cr_se<- cor(testand,HRestand, method = 'spearman'))
## [1] -0.6447906

4. Graficar los demas modelos de crecimiento de Growthmodels

Blumberg growth model

growthblumberg <- blumberg(0:10, 10, 2, 0.5)
# Calculate inverse function
timeblumberg <- blumberg.inverse(growth, 12, 2, 0.5)

par(mfrow = c(1,2))
plot(growthblumberg, main = "Modelo Blumberg", pch = 18)
lines(growthblumberg, lwd = 0.5)
plot(timeblumberg, main = "Modelo Blumberg Inverso", pch = 18)
lines(timeblumberg, lwd = 0.5)

Brody growth model

growthbrody <- brody(0:10, 10, 5, 0.3)
# Calculate inverse function
timebrody <- brody.inverse(growthbrody, 10, 5, 0.3)

par(mfrow = c(1,2))
plot(growthbrody, main = "Modelo Brody", pch = 18)
lines(growthbrody, lwd = 0.5)
plot(timebrody, main = "Modelo Brody Inverso", pch = 18)
lines(timebrody, lwd = 0.5)

Chapman-Richards growth model

growthchapman <- chapmanRichards(0:10, 10, 0.5, 0.3, 0.5)
# Calculate inverse function
timechapman <- chapmanRichards.inverse(growth, 10, 0.5, 0.3, 0.5)

par(mfrow = c(1,2))
plot(growthchapman, main = "Modelo Chapman-Richards", pch = 18)
lines(growthchapman, lwd = 0.5)
plot(timechapman, main = "Modelo Chapman-Richards Inverso", pch = 18)
lines(timechapman, lwd = 0.5)

Generalised Logistic growth model

growthgenlogis <- generalisedLogistic(0:10, 5, 10, 0.3, 0.5, 3)
# Calculate inverse function
timegenlogis <- generalisedLogistic.inverse(growth, 5, 10, 0.3, 0.5, 3)
## Warning in log((alpha - x)/(beta * x)): NaNs produced
par(mfrow = c(1,2))
plot(growthgenlogis, main = "Modelo Generalised Logistic", pch = 18)
lines(growthgenlogis, lwd = 0.5)
plot(timegenlogis, main = "Modelo Generalised Logistic Inverso", pch = 18)
lines(timegenlogis, lwd = 0.5)

Generalised Richard growth model

growthgenrich <- generalisedRichard(0:10, 5, 10, 0.3, 0.5, 1, 3)
timegenrich <- generalisedRichard.inverse(growth, 5, 10, 0.3, 0.5, 1, 3)

par(mfrow = c(1,2))
plot(growthgenrich, main = "Modelo Generalised Richard", pch = 18)
lines(growthgenrich, lwd = 0.5)
plot(timegenrich, main = "Modelo Generalised Richard Inverso", pch = 18)
lines(timegenrich, lwd = 0.5)

Gompertz growth model

growthgompertz <- gompertz(0:10, 10, 0.5, 0.3)
# Calculate inverse function
timegompertz <- gompertz.inverse(growth, 10, 0.5, 0.3)

par(mfrow = c(1,2))
plot(growthgompertz, main = "Modelo Gompertz", pch = 18)
lines(growthgompertz, lwd = 0.5)
plot(timegompertz, main = "Modelo Gompertz Inverso", pch = 18)
lines(timegompertz, lwd = 0.5)

Logistic growth model

growthlogistic <- logistic(0:10, 10, 0.5, 0.3)
# Calculate inverse function
timelogistic <- logistic.inverse(growth, 10, 0.5, 0.3)

par(mfrow = c(1,2))
plot(growthlogistic, main = "Modelo Logistic", pch = 18)
lines(growthlogistic, lwd = 0.5)
plot(timelogistic, main = "Modelo Logistic Inverso", pch = 18)
lines(timelogistic, lwd = 0.5)

Log-logistic growth model

growthloglogistic <- loglogistic(0:10, 10, 0.5, 0.3)
# Calculate inverse function
timeloglogistic <- loglogistic.inverse(growth, 10, 0.5, 0.3)

par(mfrow = c(1,2))
plot(growthloglogistic, main = "Modelo Log-logistic", pch = 18)
lines(growthloglogistic, lwd = 0.5)
plot(timeloglogistic, main = "Modelo Log-logistic Inverso", pch = 18)
lines(timeloglogistic, lwd = 0.5)

Mitcherlich growth model

growthmitcherlich <- mitcherlich(0:10, 10, 0.5, 0.3)
# Calculate inverse function
timemitcherlich <- mitcherlich.inverse(growth, 10, 0.5, 0.3)

par(mfrow = c(1,2))
plot(growthmitcherlich, main = "Modelo Mitcherlich", pch = 18)
lines(growthmitcherlich, lwd = 0.5)
plot(timemitcherlich, main = "Modelo Mitcherlich Inverso", pch = 18)
lines(timemitcherlich, lwd = 0.5)

Morgan-Mercer-Flodin growth model

growthmorgan <- mmf(0:10, 10, 0.5, 4, 1)
# Calculate inverse function
timemorgan <- mmf.inverse(growth, 10, 0.5, 4, 1)

par(mfrow = c(1,2))
plot(growthmorgan, main = "Modelo Morgan-Mercer-Flodin", pch = 18)
lines(growthmorgan, lwd = 0.5)
plot(timemorgan, main = "Modelo Morgan-Mercer-Flodin Inverso", pch = 18)
lines(timemorgan, lwd = 0.5)

Monomolecular growth model

growthmonom <- monomolecular(0:10, 10, 0.5, 0.3)
# Calculate inverse function
timemonom <- monomolecular.inverse(growth, 10, 0.5, 0.3)

par(mfrow = c(1,2))
plot(growthmonom, main = "Modelo Monomolecular", pch = 18)
lines(growthmonom, lwd = 0.5)
plot(timemonom, main = "Modelo Monomolecular Inverso", pch = 18)
lines(timemonom, lwd = 0.5)

Negative exponential growth model

growthnegative <- negativeExponential(0:10, 1, 0.3)
# Calculate inverse function
timenegative <- negativeExponential.inverse(growth, 10, 0.3)

par(mfrow = c(1,2))
plot(growthnegative, main = "Modelo Negative exponential", pch = 18)
lines(growthnegative, lwd = 0.5)
plot(timenegative, main = "Modelo Negative exponential Inverso", pch = 18)
lines(timenegative, lwd = 0.5)

Richard growth model

growthrichard <- richard(0:10, 10, 0.5, 0.3, 0.5)
timerichard <- richard.inverse(growth, 10, 0.5, 0.3, 0.5)

par(mfrow = c(1,2))
plot(growthrichard, main = "Modelo Richard", pch = 18)
lines(growthrichard, lwd = 0.5)
plot(timerichard, main = "Modelo Richard Inverso", pch = 18)
lines(timerichard, lwd = 0.5)

Schnute growth model

growthschnute <- schnute(0:10, 10, 5, .5, .5)
# Calculate inverse function
timeschnute <- schnute.inverse(growth, 10, 5, .5, .5)
## Warning in log((x^(1/m) - r0)/beta): NaNs produced
par(mfrow = c(1,2))
plot(growthschnute, main = "Modelo Schnute", pch = 18)
lines(growthschnute, lwd = 0.5)
plot(timeschnute, main = "Modelo Schnute Inverso", pch = 18)
lines(timeschnute, lwd = 0.5)

Stannard growth model

growthStannard <- stannard(0:10, 1, .2, .1, .5)
# Calculate inverse function
timeStannard <- stannard.inverse(growth, 1, .2, .1, .5)
## Warning in log((alpha/x)^(1/m) - 1): NaNs produced
par(mfrow = c(1,2))
plot(growthStannard, main = "Modelo Stannard", pch = 18)
lines(growthStannard, lwd = 0.5)

von Bertalanffy growth model

growthvon <- vonBertalanffy(0:10, 10, 0.5, 0.3, 0.5)
# Calculate inverse function
timevon <- vonBertalanffy.inverse(growth, 10, 0.5, 0.3, 0.5)

par(mfrow = c(1,2))
plot(growthvon, main = "Modelo von Bertalanffy", pch = 18)
lines(growthvon, lwd = 0.5)
plot(timevon, main = "Modelo von Bertalanffy Inverso", pch = 18)
lines(timevon, lwd = 0.5)

Weibull growth model

growthWeibull <- weibull(0:10, 10, 0.5, 0.3, 0.5)
# Calculate inverse function
timeWeibull <- weibull.inverse(growth, 10, 0.5, 0.3, 0.5)

par(mfrow = c(1,2))
plot(growthWeibull, main = "Modelo Weibull", pch = 18)
lines(growthWeibull, lwd = 0.5)
plot(timeWeibull, main = "Modelo Weibull Inverso", pch = 18)
lines(timeWeibull, lwd = 0.5)