PROBLEMA.

un hogar desea instalar un generador eólico, para hacer bien la evaluación y no quedar sub o sobre estimado. Le aconsejaron evaluar el proceso por día, en términos de la cantidad de cuartos de hora por día sin intensidad suficiente de viento. Se considera que idealmente no sobrepase de 4 horas al día sin viento suficiente. Se sabe que la intensidad de viento es de un mínimo de 4mts/s, con un máximo de 13mts/s. Para ello se encuentra disponible en DATOS, las mediciones durante un año de la velocidad del viento.

library(fitdistrplus)
library(tidyverse)

Desarrollo

¿qué se requiere averiguar? cuando no hay viento suficiente entra a funcionar las baterías, por lo anterior se sabe que las baterías pueden atender una cantidad de horas con su carga. por tanto, en el dia tengo 96 cuartos de hora y 4 horas corresponden a 16 cuartos.

definiciones

éxito: periodo de 15 min con velocidad del viento inferior a 4 mts/2 o superior a 13 mts/s

X: numero de periodos de 15 min con velocidad del viento inferior a 4 mts/2 o superior a 13 mts/s de 96 cuartos

X: numero de periodos de 15 min con velocidad del viento inferior a 4 mts/2 o superior a 13 mts/s en el dia

x: numero de periodos a revisar hasta encontrar 16 cuartos sin velocidad suficiente.

\[x\sim b(n,p)\]

\[x\sim P(\lambda)\]

\[x\sim bn(r,p)\]

Aunque tengo una idea vaga sobre la función de cuantía que se puede aplicar, para todas se necesita calcular \(p\). Entonces \(p\) se define desde:

Opción 1

sea \(y\) velocidad del viento, entonces \(p = 1- P(4<y<13)= \int_4^{13}f(y)dy\)

Luego, como \(y\) es continua por tanto se debe buscar \(f(y)\).

Analizando comportamiento de y

library(Analitica)
descripYG(datos,vel_viento)

##       n     Mean   Median       SD Kurtosis  Skewness        CV      Min
## 1 35040 8.997743 8.684297 2.979906 3.580125 0.6337527 0.3311837 1.312459
##        Max      P25      P75      IQR Fence_Low Fence_High
## 1 27.18375 6.850357 10.78281 3.932458 0.9516705    16.6815

Si bien hay una gran cantidad de datos atípicos, existe una consistencia en los datos. Como se mide todo el año, es esperable que velocidades del viento superen producto de la temporada de invierno primavera, que suele haber condiciones climáticas favorables a las velocidades encontradas. Por tanto, dado el requerimiento de análisis y modelado de datos, no es necesario replantear el problema y dividir los valores átipicos.

Dado que la forma encontrada es asimétrica, se trabajara desde ahí.

Información disponible

Se tiene conocimiento teórico que la velocidad del viento se puede modelar por una función Rayleigh (que se deriva desde la Gamma).

\[f(x) = \frac{x}{\sigma^2} \exp\left(-\frac{x^2}{2\sigma^2}\right)\] Entonces, toca definir que funciones creemos pueden comportarse de la misma forma.

Funciones Posibles:

  1. Función Gamma \[ f(y) = \frac{y^{a-1}*e^{y/b}}{b^a*\Gamma(a)}\]

  2. Función Weibull \[ f(y) = \frac{a}{b}y^{a-1}e^{-y^a/b}\]

  3. Función Log Normal \[f(x;\mu,\sigma) = \frac{1}{x\sigma\sqrt{2\pi}}\exp\left(-\frac{(\ln x - \mu)^2}{2\sigma^2}\right)\]

Definir en R

dado que la Rayleigh “no esta” en R, se define de la siguiente forma:

fry <-function(x,s){
  x / s^2 * exp(-x^2 / (2 * s^2))
}

drayleigh <- function(x, sigma) ifelse(x >= 0, x / sigma^2 * exp(-x^2 / (2 * sigma^2)), 0) # función densidad
prayleigh <- function(q, sigma) ifelse(q >= 0, 1 - exp(-q^2 / (2 * sigma^2)), 0)           # función de probabilidad
qrayleigh <- function(p, sigma) sigma * sqrt(-2 * log(1 - p))                              # función de cuantiles
rrayleigh <- function(n, sigma) sigma * sqrt(-2 * log(runif(n)))                           # generado de aleatorios

Proceso de simulación

Simulación por función

f1<-fitdist(datos$vel_viento,"rayleigh",
            start = list(sigma = sqrt(mean(datos$vel_viento^2) / 2)),
            lower = 1e-6)
f2<-fitdist(datos$vel_viento,"gamma")
f3<-fitdist(datos$vel_viento,"weibull")
f4<-fitdist(datos$vel_viento,"lnorm")

Análisis de simulaciones

plot.legend<-c("Rayleigh","Gamma","Weibull","LogNormal")
ppcomp(list(f1,f2,f3,f4),legendtext = plot.legend)

qqcomp(list(f1,f2,f3,f4),legendtext = plot.legend)

denscomp(list(f1,f2,f3,f4),legendtext = plot.legend)

cdfcomp(list(f1,f2,f3,f4),legendtext = plot.legend)

Claramente se descarta la f Rayleighy podría estar entre Gamma y LogNormal

Usamos Criterio Akaike o Bayes

gofstat(list(f1,f2,f3,f4))
## Goodness-of-fit statistics
##                              1-mle-rayleigh 2-mle-gamma 3-mle-weibull
## Kolmogorov-Smirnov statistic      0.1800289 0.005910395    0.04038657
## Cramer-von Mises statistic      393.6999171 0.150036143   23.56282021
## Anderson-Darling statistic     2237.8368959 0.790460577  166.31341735
##                              4-mle-lnorm
## Kolmogorov-Smirnov statistic  0.02683717
## Cramer-von Mises statistic    8.60147060
## Anderson-Darling statistic   52.62625274
## 
## Goodness-of-fit criteria
##                                1-mle-rayleigh 2-mle-gamma 3-mle-weibull
## Akaike's Information Criterion       186700.9    173523.7      175834.1
## Bayesian Information Criterion       186709.4    173540.7      175851.0
##                                4-mle-lnorm
## Akaike's Information Criterion    174292.9
## Bayesian Information Criterion    174309.8

Ambos criterios coinciden (AIC y BIC), siendo los menos indices Gamma y Log NOrmal.

Validando función definida.

ks.test(datos$vel_viento,"pgamma",shape=f2$estimate["shape"],rate=f2$estimate["rate"])
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  datos$vel_viento
## D = 0.0059104, p-value = 0.1728
## alternative hypothesis: two-sided

dado H0: los datos presentan un comportamiento aproximadamente Gamma, NO existe evidencia para rechazar esta afirmación. Entonces, podemos usar esta función para el calculo de \(p\).

entonces \(p=1-P(4<y<13)=\int_4^{13}f(y)dy\)

en r, calcularemos \(1-(P(y<13)-P(y<4))\)

1-pgamma(13,shape=f2$estimate["shape"],rate=f2$estimate["rate"])+
  pgamma(4,shape=f2$estimate["shape"],rate=f2$estimate["rate"])
## [1] 0.1199858

Luego de cada 100 horas hay 12 horas, que no tienen suficiente viento para el generador

Entonces, como tengo el valor de p, ahora es que distribución ocupo.

Modelar la variable discreta.

Descartando funciones

en la Binomial y en la Poisson la distribución de los cuartos de hora es aleatoria, y requiero que el proceso sea continuado, ya que con 4 horas de uso continuo se descarga la batería y me quedo sin energía. En la Binomial, pueden existir espacios de carga de batería. Luego el que resuelve esto es la Binomial Negativa.

entonces,

si definimos \(X\) como: el numero de cuartos de hora hasta encontrar 17 sin velocidad suficiente (ya que no se produciría energía por el generador y no habría energía en las baterías.)

luego en 96 cuartos, pueden darse varios periodos consecutivos.

Modelo de probabilidad \(P(X\le 96)=\sum_{x=17}^{96-17}P(x)\)

p <- 0.12; r <- 17

pnbinom(96 - r, r, p)       # cuanto sin energía
## [1] 0.06444134
1 - pnbinom(96 - r,r,p)     # cuanto con energía
## [1] 0.9355587

Entonces

Respuesta

se lo lleva a días, aproximadamente 24 días al año, no habría viento suficiente para generar energía. En términos de horas, se esperaría que de cada 24 horas, aprox 1.5 horas no haya velocidad de viento suficiente.