Parametry rozkładu Weibulla i dla danych prędkości wiatru

Ponieważ @zaynah opublikował w komentarzach, że uważa się, że dane są zgodne z rozkładem Weibulla, przedstawię krótki samouczek na temat szacowania parametrów takiego rozkładu przy użyciu MLE (szacowanie maksymalnego prawdopodobieństwa). Na stronie jest podobny post dotyczący prędkości wiatru i rozkładu Weibulla.

Pobierz i zainstalujR , to nic nie kosztuje
Opcjonalnie: Pobierz i zainstaluj RStudio , które jest doskonałym IDE dla R, zapewniając mnóstwo przydatnych funkcji, takich jak podświetlanie składni i więcej.
Instaluje pakiety MASSi carprzez wpisywanie: install.packages(c("MASS", "car")). Załaduj je, wpisując: library(MASS)i library(car).
Zaimportuj swoje dane doR . Jeśli masz na przykład swoje dane w programie Excel, zapisz je jako plik tekstowy z ogranicznikami (.txt) i zaimportuj za Rpomocą read.table.
Użyj funkcji fitdistrdo obliczania maksymalnych szacunków prawdopodobieństwa swojej rozkładu Weibulla: fitdistr(my.data, densfun="weibull", lower = 0). Aby zobaczyć w pełni opracowany przykład, zobacz link na dole odpowiedzi.
Wykonaj wykres QQ, aby porównać swoje dane z rozkładem Weibulla z parametrami skali i kształtu oszacowanymi w punkcie 5: qqPlot(my.data, distribution="weibull", shape=, scale=)

Tutorial Vito Ricci na dopasowanie rozkładu ze Rjest to dobry punkt wyjścia w tej sprawie. Na tej stronie znajduje się wiele postów na ten temat (zobacz też ten post ).

Aby zobaczyć w pełni opracowany przykład użycia fitdistr, zobacz ten post .

Spójrzmy na przykład w 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)

Szacunki maksymalnego prawdopodobieństwa są zbliżone do tych, które arbitralnie ustaliliśmy przy generowaniu liczb losowych. Porównajmy nasze dane za pomocą wykresu QQ z hipotetycznym rozkładem Weibulla z parametrami, które oszacowaliśmy fitdistr:

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

QQPlot

Punkty są ładnie wyrównane na linii i głównie w granicach przedziału ufności 95%. Stwierdzilibyśmy, że nasze dane są zgodne z dystrybucją Weibull. Było to oczywiście oczekiwane, ponieważ próbkowaliśmy nasze wartości z rozkładu Weibulla.

Oszacowanie (kształt) (skala) rozkładu Weibulla bez MLE $k$ $c$

W tym artykule wymieniono pięć metod szacowania parametrów rozkładu Weibulla dla prędkości wiatru. Wyjaśnię tutaj trzy z nich.

Od średnich i odchyleń standardowych

Parametr kształtu jest szacowany jako: a parametr skali jest szacowany jako: with to średnia prędkość wiatru, a odchylenie standardowe, a to funkcja Gamma . $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$

Najmniejsze kwadraty pasują do obserwowanego rozkładu

Jeśli obserwowane prędkości wiatru są podzielone na przedział prędkości , o częstotliwości występowania i skumulowane częstotliwości , możesz dopasować regresję liniową postaci do wartości parametrów dla danego Weibulla są powiązane ze współczynnikami liniowymi i o $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_{i} = \ln (V_{i})

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

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

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

a

$a$

b

$b$

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

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

k = b

$k=b$

Średnie i kwartylowe prędkości wiatru

Jeśli nie masz pełnych obserwowanych prędkości wiatru, ale mediana i kwartyle i , a następnie i można obliczyć za pomocą relacji $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}$

Porównanie czterech metod

Oto przykład w Rporównaniu czterech metod:

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

Wszystkie metody dają bardzo podobne wyniki. Podejście oparte na maksymalnym prawdopodobieństwie ma tę zaletę, że standardowe błędy parametrów Weibulla są podawane bezpośrednio.

Używanie bootstrap do dodawania punktowych przedziałów ufności do pliku PDF lub CDF

Możemy użyć nieparametrycznego bootstrapu do skonstruowania punktowych przedziałów ufności wokół PDF i CDF szacowanego rozkładu Weibulla. Oto Rskrypt:

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

Elementy CI Weibull PDF

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

Elementy CI Weibull CDF

COOLSerdash
źródło

+1, niezły przegląd. Uwaga: niewielkim skrótem może być użycie ? QqPlot w / distribution=weibullz pakietu samochodowego, który dopasuje parametry przez MLE i utworzy wykres qq w 1 kroku.

gung - Przywróć Monikę

@gung Dzięki. Nie wiem, czy qqPlot z automatycznie caroblicza parametry MLE. Jeśli wygenerować zmienną losową o rozkładzie Weibulla ( rweibull) i użyć polecenia qqPlot(rw, distribution="weibull")pojawia się komunikat o błędzie, który musi zapewnić parametry shapei scaledo qqPlot. Czy coś brakuje?

COOLSerdash

mój błąd. Oczywiście, tylko automatycznie szacuje parametry z niektórych dystrybucji, a Weibull nie jest jednym z nich.

gung - Przywróć Monikę

cześć, znalazłem, że po zaimportowaniu moich danych do R, kiedy wykonuję polecenie, fitdistr (mydata, densfun = "weibull") mówi komunikat o błędzie, że „mydata” nie znaleziono .. w rzeczywistości moje dane zostały zaimportowane do R. każda odpowiedź byłaby mile widziana.

Zay

@zaynah Czy możesz edytować swoją odpowiedź i opublikować kod, którego używasz do importowania danych? Dodaj również komunikat o błędzie. Czy możesz zaimportować dane bez błędów? Czy sprawdziłeś, czy dane zostały poprawnie zaimportowane?