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

19

Hi peut également être montré pour obtenir le paramètre de forme et d'échelle pour la méthode modifiée du maximum de vraisemblance

Zay
la source
2
Salut @zaynah et bienvenue sur le site. Je ne sais pas si votre question est de savoir si vos données sont compatibles avec une distribution weibull ou quels seraient les paramètres d'une distribution weibull décrivant vos données. Si vous supposez que vos données suivent une distribution weibull et que vous souhaitez rechercher les paramètres, vous pouvez utiliser fitdistr(mydata, densfun="weibull")in Rpour rechercher les paramètres via MLE. Pour faire un graphique, utilisez la qqPlotfonction du carpackage: qqPlot(mydata, distribution="weibull", shape=, scale=)avec les paramètres de forme et d'échelle que vous avez trouvés avec fitdistr.
COOLSerdash
salut merci pour une réponse rapide, mes données sont la vitesse moyenne du vent mensuelle pendant 5 ans, elle est compatible avec weibull. le problème est que je ne sais pas comment trouver k et c c'est-à-dire les paramètres de weibull .. et je ne sais pas comment comparer les données expérimentales avec weibull ... aussi ce qui est MLE ... :(
Zay
MLE = Estimation du maximum de vraisemblance. Je ne sais pas quel logiciel vous utilisez, mais dans R, qui est disponible gratuitement, vous pouvez installer et charger le package MASSet utiliser fitdistravec vos données pour calculer les estimations de k et c. Et puis, vous pouvez comparer vos données avec le weibull avec les paramètres estimés à l'aide qqPlotdu carpackage.
COOLSerdash
merci beaucoup COOlserdash, je télécharge le logiciel R.
Zay
1
Ok, voici un tutoriel étape par étape: 1. Téléchargez et installez R. 2. Installez les paquets MASSet caren tapant: install.packages(c("MASS", "car")). Chargez les packages en tapant: library(MASS)et library(car). 3. Importez vos données dansR , de préférence avec un fichier .txt. 4. Si vos données est appelée l' my.datautilisation fitdistrde la manière suivante: fitdistr(my.data, distribution="weibull"). 5. Faites un graphique comme je l'ai décrit dans le premier commentaire avec qqPlot.
COOLSerdash

Réponses:

28

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.

  1. Téléchargez et installezR , c'est gratuit
  2. 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.
  3. Installez les paquets MASSet caren tapant: install.packages(c("MASS", "car")). Chargez-les en tapant: library(MASS)et library(car).
  4. 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.
  5. 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.
  6. 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 c (échelle) d'une distribution de Weibull sans MLEkc

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=(σ^v^)-1.086
c
c=v^Γ(1+1/k)
v^σ^Γ

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

n0-V1,V1-V2,,Vn-1-VnF1,F2,,Fnp1=F1,p2=F1+F2,,pn=pn-1+Fny=une+bX

Xje=ln(Vje)
yje=ln[-ln(1-pje)]
uneb
c=exp(-uneb)
k=b

Vitesse du vent médiane et quartile

VmV0,25V0,75 [p(VV0,25)=0,25,p(VV0,75)=0,75]ck

k=ln[ln(0,25)/ln(0,75)]/ln(V0,75/V0,25)1,573/ln(V0,75/V0,25)
c=Vm/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?
COOLSerdash