Generamos 200 valores aleatorios entre 0 y 1 con distribucion uniforme.
rnd<-runif(200)
hist(rnd)
Llevamos los valores aleatorios a una distribución normal, aplicando la siguiente funcion logaritmica (notamos que tiene un sesgo positivo)
rdm = (log((1+rnd)/(1-rnd)))/1.82
hist(rdm)
Corremos una simulacion con el modelo de valores aleatorios antes descrito, llamamos la modelo y pasamos como parametro un vector con los valores correspondientes a la variabilidad. El boxplot muestra el rendimiento del tuberculo fresco, apreciamo que hay presencia de variabilidad de un escenario a otro lo cual es correcto pero la mediana no debería aumentar.
source("http://inrm.cip.cgiar.org/home/solanum/001_aleatoriedad/Module_PotentialGrowth1.R")
salida<-modelo(c(0.075,0.1,0.15,0.20,0.25,0.30))
boxplot(salida,main="Fresh Tuber Yield (t/ha)",xlab="variability",ylab="t/ha")
muestro la media para cada escenario, se aprecia que le media aumenta de escenario en escenario
apply(salida,2,mean)
## 0.075 0.1 0.15 0.2 0.25 0.3
## 63.26 63.99 66.89 68.84 70.50 72.42
muestro la mediana para cada escenario
apply(salida,2,median)
## 0.075 0.1 0.15 0.2 0.25 0.3
## 62.52 62.99 65.41 67.21 68.73 69.76
Genero el histograma para cada escenario. Las grafica de frecuencias muestra un sesgo positivo similar a los valores aleatorios generados al inicio, el cual no es correcto.
par(mfrow=c(3,2))
hist(salida$"0.075",main="Fresh tuber yield. variability=0.075")
hist(salida$"0.1",main="Fresh tuber yield. variability=0.1")
hist(salida$"0.15",main="Fresh tuber yield. variability=0.15")
hist(salida$"0.2",main="Fresh tuber yield. variability=0.20")
hist(salida$"0.25",main="Fresh tuber yield. variability=0.25")
hist(salida$"0.3",main="Fresh tuber yield. variability=0.30")
El sesgo positivo de los rendimientos no es correcto, deberia mostrar una distribucion normal sin sesgo, para ello aplicaremos un cambio de escala a los valores aleatorios
Generamos 200 valores aleatorios con distribucion uniforme como se muestra en el histograma
rnd<-runif(200)
hist(rnd)
aqui se muestra el rango de los valores aleatorios, su rango esta entre 0 y 1
summary(rnd)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0042 0.2650 0.5070 0.5000 0.7370 0.9880
aplicamos un cambio de escala a los valores aleatorios (VA), primero multiplicamos el VA por 2, asi los VA ya no se moveran en el intervalo de 0 y 1 sino entre 0 y 2, luego hacemos una traslación restando la unidad a los VA, ahora los VA se distribuyen entre -1 y 1. Fijarse en el valor minimo y maximo
rnd<-rnd*2-1
summary(rnd)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.9920 -0.4690 0.0131 -0.0001 0.4730 0.9760
normalizamos los valores aleatorios y ahora obtendremos una distribución normal sin sesgo
rdm = (log((1+rnd)/(1-rnd)))/1.82
hist(rdm)
Este modelo de aleatoriedad lo aplicamos en el modelo potencial obtenemos lo siguiente
source("http://inrm.cip.cgiar.org/home/solanum/001_aleatoriedad/Module_PotentialGrowth2.R")
salida<-modelo(c(0.075,0.1,0.15,0.20,0.25,0.30))
boxplot(salida,main="Fresh Tuber Yield (t/ha)",xlab="variability",ylab="t/ha")
muestro la media para cada escenario, se aprecia que el valor de la media se mantiene
apply(salida,2,mean)
## 0.075 0.1 0.15 0.2 0.25 0.3
## 58.65 60.56 60.16 59.87 62.14 62.96
muestro la mediana, muestra un comportamiento similar ala media
apply(salida,2,median)
## 0.075 0.1 0.15 0.2 0.25 0.3
## 58.17 59.53 59.45 59.92 62.69 61.34
Genero el histograma para cada escenario. Las grafica de frecuencias ya no muestran un sesgo positivo sino mas bien se acercan a la simetria.
par(mfrow=c(3,2))
hist(salida$"0.075",main="Fresh tuber yield. variability=0.075")
hist(salida$"0.1",main="Fresh tuber yield. variability=0.1")
hist(salida$"0.15",main="Fresh tuber yield. variability=0.15")
hist(salida$"0.2",main="Fresh tuber yield. variability=0.20")
hist(salida$"0.25",main="Fresh tuber yield. variability=0.25")
hist(salida$"0.3",main="Fresh tuber yield. variability=0.30")
Para verificar la variabilidad, calculo el coeficiente de variación (%) para cada escenario, el coeficiente de variacion aumenta de un escenario a otro.
co.var<-function(x)((sd(x)/mean(x))*100)
apply(salida,2,co.var)
## 0.075 0.1 0.15 0.2 0.25 0.3
## 7.297 10.845 13.282 19.571 24.005 27.137