Weibull-Verteilungsparameter

Da in den Kommentaren von @zaynah angegeben ist, dass die Daten vermutlich einer Weibull-Verteilung folgen, werde ich ein kurzes Tutorial zur Schätzung der Parameter einer solchen Verteilung mithilfe von MLE (Maximum Likelihood Estimation) geben. Es gibt einen ähnlichen Beitrag über Windgeschwindigkeiten und Weibull-Verteilung auf der Website.

Herunterladen und installierenR , es ist kostenlos
Optional: Laden Sie RStudio herunter und installieren Sie es. RStudio ist eine großartige IDE für R mit zahlreichen nützlichen Funktionen wie Syntaxhervorhebung und vielem mehr.
Installieren Sie die Pakete MASSund cardurch Eingabe: install.packages(c("MASS", "car")). Laden Sie sie durch Eingabe von: library(MASS)und library(car).
Importieren Sie Ihre Daten inR . Wenn Sie Ihre Daten beispielsweise in Excel haben, speichern Sie sie als durch Trennzeichen getrennte Textdatei (TXT) und importieren Sie sie Rmit read.table.
Verwenden Sie die Funktion , fitdistrdie Maximum - Likelihood - Schätzungen Ihrer Weibull - Verteilung zu berechnen: fitdistr(my.data, densfun="weibull", lower = 0). Ein ausführliches Beispiel finden Sie unter dem Link am Ende der Antwort.
Erstellen Sie einen QQ-Plot, um Ihre Daten mit einer Weibull-Verteilung zu vergleichen, wobei die Skalierungs- und Formparameter auf Punkt 5 geschätzt werden: qqPlot(my.data, distribution="weibull", shape=, scale=)

Das Tutorial von Vito Ricci über das Anpassen der Distribution mit Rist ein guter Ausgangspunkt in dieser Angelegenheit. Und es gibt zahlreiche Beiträge auf dieser Seite zu diesem Thema (siehe auch diesen Beitrag ).

fitdistrSchauen Sie sich diesen Beitrag an, um ein vollständig ausgearbeitetes Beispiel für die Verwendung zu sehen .

Schauen wir uns ein Beispiel an 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)

Die maximalen Wahrscheinlichkeitsschätzungen liegen in der Nähe derjenigen, die wir bei der Erzeugung der Zufallszahlen willkürlich festgelegt haben. Vergleichen wir unsere Daten mit einem QQ-Plot mit einer hypothetischen Weibull-Verteilung mit den Parametern, mit denen wir geschätzt haben fitdistr:

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

QQPlot

Die Punkte sind gut auf der Linie ausgerichtet und liegen größtenteils innerhalb des 95% -Vertrauensbereichs. Wir würden daraus schließen, dass unsere Daten mit einer Weibull-Distribution kompatibel sind. Dies wurde natürlich erwartet, da wir unsere Werte aus einer Weibull-Verteilung abgetastet haben.

Abschätzung von (Form) und (Skala) einer Weibull-Verteilung ohne MLE $k$ $c$

In diesem Artikel werden fünf Methoden zur Schätzung der Parameter einer Weibull-Verteilung für Windgeschwindigkeiten aufgeführt. Ich erkläre hier drei davon.

Mittelwert und Standardabweichung

$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$

Die kleinsten Quadrate passen zur beobachteten Verteilung

$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_{ich} = \ln (V_{ich})

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

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

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

a

$a$

b

$b$

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

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

k = b

$k=b$

Median und Quartil Windgeschwindigkeiten

$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}$

Vergleich der vier Methoden

Hier ist ein Beispiel für den RVergleich der vier Methoden:

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

Alle Methoden führen zu sehr ähnlichen Ergebnissen. Der Maximum-Likelihood-Ansatz hat den Vorteil, dass die Standardfehler der Weibull-Parameter direkt angegeben werden.

Verwenden von Bootstrap, um der PDF- oder CDF-Datei punktweise Konfidenzintervalle hinzuzufügen

Wir können einen nicht-parametrischen Bootstrap verwenden, um punktweise Konfidenzintervalle um die PDF- und CDF-Datei der geschätzten Weibull-Verteilung zu erstellen. Hier ist ein RSkript:

#-----------------------------------------------------------------------------
# 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")

Weibull PDF CIs

#-----------------------------------------------------------------------------
# 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")

Weibull CDF CIs

COOLSerdash
quelle

+1, schöne Übersicht. NB, eine leichte Abkürzung könnte sein, qqPlot w / distribution=weibullaus dem Autopaket zu verwenden , das über MLE zu den Parametern passt und das qq-Diagramm in einem Schritt erstellt.

gung - Reinstate Monica

@gung Danke. Mir ist nicht bekannt, dass qqPlot von cardie MLE-Parameter automatisch berechnet. Wenn ich eine Zufallsvariable mit einer Weibull-Distribution ( rweibull) generiere und den Befehl verwende, qqPlot(rw, distribution="weibull")erhalte ich eine Fehlermeldung, die besagt, dass die Parameter shapeund scaleto angegeben werden müssen qqPlot. Vermisse ich etwas?

COOLSerdash

mein Fehler. Offensichtlich werden nur die Parameter einiger Distributionen automatisch geschätzt, und Weibull gehört nicht dazu.

gung - Reinstate Monica

Hallo, ich habe festgestellt, dass nach dem Import von Mydata in R, wenn ich den Befehl "fitdistr" (mydata, densfun = "weibull") ausführe, "mydata" nicht gefunden wurde. Tatsächlich wurden meine Daten in R importiert. Jede Antwort wäre willkommen.

Zay

@zaynah Könnten Sie bitte Ihre Antwort bearbeiten und Ihren Code posten, den Sie zum Importieren der Daten verwenden. Bitte fügen Sie auch die Fehlermeldung hinzu. Könnten Sie die Daten ohne Fehler importieren? Haben Sie überprüft, ob die Daten korrekt importiert wurden?