Paramètres de distribution de Weibull et pour les données de vitesse du vent

Parce que @zaynah a posté dans les commentaires que les données sont censées suivre une distribution de Weibull, je vais fournir un bref tutoriel sur la façon d'estimer les paramètres d'une telle distribution en utilisant MLE (estimation de vraisemblance maximale). Il existe un article similaire sur les vitesses du vent et la distribution de Weibull sur le site.

Téléchargez et installezR , c'est gratuit
Facultatif: Téléchargez et installez RStudio , qui est un excellent IDE pour R fournissant une tonne de fonctions utiles telles que la coloration syntaxique et plus encore.
Installez les paquets MASSet caren tapant: install.packages(c("MASS", "car")). Chargez-les en tapant: library(MASS)et library(car).
Importez vos données dansR . Si vous avez vos données dans Excel, par exemple, enregistrez-les en tant que fichier texte délimité (.txt) et importez-les Ravec read.table.
Utilisez la fonction fitdistrpour calculer les estimations du maximum de vraisemblance de votre distribution de Weibull: fitdistr(my.data, densfun="weibull", lower = 0). Pour voir un exemple entièrement élaboré, consultez le lien au bas de la réponse.
Faites un QQ-Plot pour comparer vos données avec une distribution de Weibull avec les paramètres d'échelle et de forme estimés au point 5: qqPlot(my.data, distribution="weibull", shape=, scale=)

Le tutoriel de Vito Ricci sur l'ajustement de la distribution avec Rest un bon point de départ en la matière. Et il y a de nombreux articles sur ce site sur le sujet (voir aussi cet article ).

Pour voir un exemple entièrement élaboré de la façon d'utiliser fitdistr, jetez un œil à cet article .

Regardons un exemple dans R:

# Load packages

library(MASS)
library(car)

# First, we generate 1000 random numbers from a Weibull distribution with
# scale = 1 and shape = 1.5

rw <- rweibull(1000, scale=1, shape=1.5)

# We can calculate a kernel density estimation to inspect the distribution
# Because the Weibull distribution has support [0,+Infinity), we are truncate
# the density at 0

par(bg="white", las=1, cex=1.1)
plot(density(rw, bw=0.5, cut=0), las=1, lwd=2,
xlim=c(0,5),col="steelblue")

Weibull KDE

# Now, we can use fitdistr to calculate the parameters by MLE
# The option "lower = 0" is added because the parameters of the Weibull distribution need to be >= 0

fitdistr(rw, densfun="weibull", lower = 0)

     shape        scale   
  1.56788999   1.01431852 
 (0.03891863) (0.02153039)

Les estimations du maximum de vraisemblance sont proches de celles que nous avons arbitrairement fixées dans la génération des nombres aléatoires. Comparons nos données en utilisant un QQ-Plot avec une distribution de Weibull hypothétique avec les paramètres que nous avons estimés avec fitdistr:

qqPlot(rw, distribution="weibull", scale=1.014, shape=1.568, las=1, pch=19)

QQPlot

Les points sont bien alignés sur la ligne et principalement dans l'enveloppe de confiance à 95%. Nous concluons que nos données sont compatibles avec une distribution de Weibull. Cela était attendu, bien sûr, car nous avons échantillonné nos valeurs à partir d'une distribution de Weibull.

Estimation de (forme) et (échelle) d'une distribution de Weibull sans MLE $k$ $c$

Cet article énumère cinq méthodes pour estimer les paramètres d'une distribution de Weibull pour les vitesses du vent. Je vais en expliquer trois ici.

Des moyennes et de l'écart type

$k$

k = {(\frac{\hat{σ}}{\hat{v}})}^{- 1.086}

$k=\left(\frac{\hat{\sigma}}{\hat{v}}\right)^{-1.086}$

c

$c$

c = \frac{\hat{v}}{Γ (1 + 1 / k)}

$c=\frac{\hat{v}}{\Gamma(1+1/k)}$

\hat{v}

$\hat{v}$

\hat{σ}

$\hat{\sigma}$

Γ

$\Gamma$

Les moindres carrés correspondent à la distribution observée

$n$ $0-V_{1},V_{1}-V_{2},\ldots, V_{n-1}-V_{n}$ $f_{1}, f_{2},\ldots,f_{n}$ $p_{1}=f_{1}, p_{2}=f_{1}+f_{2}, \ldots, p_{n}=p_{n-1}+f_{n}$ $y=a+bx$

X_{je} = \ln (V_{je})

$x_{i} = \ln(V_{i})$

y_{je} = \ln [- \ln (1 - p_{je})]

$y_{i} = \ln[-\ln(1-p_{i})]$

a

$a$

b

$b$

c = \exp (- \frac{une}{b})

$c=\exp\left(-\frac{a}{b}\right)$

k = b

$k=b$

Vitesse du vent médiane et quartile

$V_{m}$ $V_{0.25}$ $V_{0.75}$ $\left[p(V\leq V_{0.25})=0.25, p(V\leq V_{0.75})=0.75\right]$ $c$ $k$

k = \ln [\ln (0,25) / \ln (0,75)] / \ln (V_{0,75} / V_{0,25}) \approx 1,573 / \ln (V_{0,75} / V_{0,25})

$k = \ln\left[\ln(0.25)/\ln(0.75)\right]/\ln(V_{0.75}/V_{0.25})\approx 1.573/\ln(V_{0.75}/V_{0.25})$

c = V_{m} / \ln (2)^{1 / k}

$c=V_{m}/\ln(2)^{1/k}$

Comparaison des quatre méthodes

Voici un exemple de Rcomparaison des quatre méthodes:

library(MASS)  # for "fitdistr"

set.seed(123)
#-----------------------------------------------------------------------------
# Generate 10000 random numbers from a Weibull distribution
# with shape = 1.5 and scale = 1
#-----------------------------------------------------------------------------

rw <- rweibull(10000, shape=1.5, scale=1)

#-----------------------------------------------------------------------------
# 1. Estimate k and c by MLE
#-----------------------------------------------------------------------------

fitdistr(rw, densfun="weibull", lower = 0)
shape         scale   
1.515380298   1.005562356 

#-----------------------------------------------------------------------------
# 2. Estimate k and c using the leas square fit
#-----------------------------------------------------------------------------

n <- 100 # number of bins
breaks <- seq(0, max(rw), length.out=n)

freqs <- as.vector(prop.table(table(cut(rw, breaks = breaks))))
cum.freqs <- c(0, cumsum(freqs)) 

xi <- log(breaks)
yi <- log(-log(1-cum.freqs))

# Fit the linear regression
least.squares <- lm(yi[is.finite(yi) & is.finite(xi)]~xi[is.finite(yi) & is.finite(xi)])
lin.mod.coef <- coefficients(least.squares)

k <- lin.mod.coef[2]
k
1.515115
c <- exp(-lin.mod.coef[1]/lin.mod.coef[2])
c
1.006004

#-----------------------------------------------------------------------------
# 3. Estimate k and c using the median and quartiles
#-----------------------------------------------------------------------------

med <- median(rw)
quarts <- quantile(rw, c(0.25, 0.75))

k <- log(log(0.25)/log(0.75))/log(quarts[2]/quarts[1])
k
1.537766
c <- med/log(2)^(1/k)
c
1.004434

#-----------------------------------------------------------------------------
# 4. Estimate k and c using mean and standard deviation.
#-----------------------------------------------------------------------------

k <- (sd(rw)/mean(rw))^(-1.086)
c <- mean(rw)/(gamma(1+1/k))
k
1.535481
c
1.006938

Toutes les méthodes donnent des résultats très similaires. L'approche du maximum de vraisemblance présente l'avantage de donner directement les erreurs standard des paramètres de Weibull.

Utilisation du bootstrap pour ajouter des intervalles de confiance ponctuels au PDF ou au CDF

Nous pouvons utiliser un bootstrap non paramétrique pour construire des intervalles de confiance ponctuels autour du PDF et du CDF de la distribution de Weibull estimée. Voici un Rscript:

#-----------------------------------------------------------------------------
# 5. Bootstrapping the pointwise confidence intervals
#-----------------------------------------------------------------------------

set.seed(123)

rw.small <- rweibull(100,shape=1.5, scale=1)

xs <- seq(0, 5, len=500)


boot.pdf <- sapply(1:1000, function(i) {
  xi <- sample(rw.small, size=length(rw.small), replace=TRUE)
  MLE.est <- suppressWarnings(fitdistr(xi, densfun="weibull", lower = 0))  
  dweibull(xs, shape=as.numeric(MLE.est[[1]][13]), scale=as.numeric(MLE.est[[1]][14]))
}
)

boot.cdf <- sapply(1:1000, function(i) {
  xi <- sample(rw.small, size=length(rw.small), replace=TRUE)
  MLE.est <- suppressWarnings(fitdistr(xi, densfun="weibull", lower = 0))  
  pweibull(xs, shape=as.numeric(MLE.est[[1]][15]), scale=as.numeric(MLE.est[[1]][16]))
}
)   

#-----------------------------------------------------------------------------
# Plot PDF
#-----------------------------------------------------------------------------

par(bg="white", las=1, cex=1.2)
plot(xs, boot.pdf[, 1], type="l", col=rgb(.6, .6, .6, .1), ylim=range(boot.pdf),
     xlab="x", ylab="Probability density")
for(i in 2:ncol(boot.pdf)) lines(xs, boot.pdf[, i], col=rgb(.6, .6, .6, .1))

# Add pointwise confidence bands

quants <- apply(boot.pdf, 1, quantile, c(0.025, 0.5, 0.975))
min.point <- apply(boot.pdf, 1, min, na.rm=TRUE)
max.point <- apply(boot.pdf, 1, max, na.rm=TRUE)
lines(xs, quants[1, ], col="red", lwd=1.5, lty=2)
lines(xs, quants[3, ], col="red", lwd=1.5, lty=2)
lines(xs, quants[2, ], col="darkred", lwd=2)
#lines(xs, min.point, col="purple")
#lines(xs, max.point, col="purple")

CI PDF Weibull

#-----------------------------------------------------------------------------
# Plot CDF
#-----------------------------------------------------------------------------

par(bg="white", las=1, cex=1.2)
plot(xs, boot.cdf[, 1], type="l", col=rgb(.6, .6, .6, .1), ylim=range(boot.cdf),
     xlab="x", ylab="F(x)")
for(i in 2:ncol(boot.cdf)) lines(xs, boot.cdf[, i], col=rgb(.6, .6, .6, .1))

# Add pointwise confidence bands

quants <- apply(boot.cdf, 1, quantile, c(0.025, 0.5, 0.975))
min.point <- apply(boot.cdf, 1, min, na.rm=TRUE)
max.point <- apply(boot.cdf, 1, max, na.rm=TRUE)
lines(xs, quants[1, ], col="red", lwd=1.5, lty=2)
lines(xs, quants[3, ], col="red", lwd=1.5, lty=2)
lines(xs, quants[2, ], col="darkred", lwd=2)
lines(xs, min.point, col="purple")
lines(xs, max.point, col="purple")

CI CDF de Weibull

COOLSerdash
la source

+1, belle vue d'ensemble. NB, un léger raccourci pourrait être d'utiliser ? QqPlot w / distribution=weibullfrom the car package, qui s'adaptera aux paramètres via MLE et fera le qq-plot en 1 étape.

gung - Rétablir Monica

@gung Merci. Je ne suis pas au courant que qqPlot de carcalcule automatiquement les paramètres MLE. Si je génère une variable aléatoire avec une distribution weibull ( rweibull) et utilise la commande, qqPlot(rw, distribution="weibull")j'obtiens un message d'erreur disant que doit fournir les paramètres shapeet scaleà qqPlot. Suis-je en train de manquer quelque chose?

COOLSerdash

mon erreur. De toute évidence, il estime automatiquement uniquement les paramètres de certaines distributions, et Weibull n'en fait pas partie.

gung - Rétablir Monica

salut, j'ai trouvé qu'après avoir importé mes données dans R, quand je fais la commande fitdistr (mydata, densfun = "weibull"), il dit erreur messgae que "mydata" est introuvable .. en fait mes données ont été importées dans R. toute réponse serait la bienvenue.

Zay

@zaynah Pourriez-vous s'il vous plaît modifier votre réponse et publier votre code que vous utilisez pour importer les données. Veuillez également ajouter le message d'erreur. Pourriez-vous importer les données sans erreur? Avez-vous vérifié si les données ont été importées correctement?