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


19

Cześć można to samo pokazać, aby uzyskać parametr kształtu i skali dla zmodyfikowanej metody największego prawdopodobieństwa


2
Cześć @zaynah i witam na stronie. Nie jestem pewien, czy twoje pytanie brzmi, czy twoje dane są zgodne z rozkładem weibulla, czy jakie parametry rozkładu weibulla opisujące twoje dane. Jeśli zakładasz, że twoje dane są zgodne z rozkładem Weibulla i chcesz znaleźć parametry, możesz użyć fitdistr(mydata, densfun="weibull")w, Raby znaleźć parametry za pośrednictwem MLE. Aby utworzyć wykres, użyj qqPlotfunkcji z carpakietu: qqPlot(mydata, distribution="weibull", shape=, scale=)z parametrami kształtu i skali, które znalazłeś fitdistr.
COOLSerdash

cześć dzięki za szybką odpowiedź, moje dane to średnia miesięczna prędkość wiatru przez 5 lat, jest kompatybilna z weibull. Problem polega na tym, że nie wiem, jak znaleźć k i c, tj. parametry weibulla ... i nie wiem, jak porównać dane eksperymentalne z weibullem ... także co to jest MLE ... :(
Zay

MLE = oszacowanie maksymalnego prawdopodobieństwa. Nie wiem, jakiego oprogramowania używasz, ale w R, który jest swobodnie dostępny, możesz zainstalować i załadować pakiet MASSoraz użyć fitdistrdanych, aby obliczyć szacunkowe wartości k i c. Następnie możesz porównać swoje dane z weibull z oszacowanymi parametrami przy użyciu qqPlotz carpakietu.
COOLSerdash

wielkie dzięki COOlserdash, pobieram oprogramowanie R.
Zay

1
Ok, oto krok po kroku: 1. Pobierz i zainstaluj R. 2. Instalacja pakietów MASSi carprzez wpisywanie: install.packages(c("MASS", "car")). Załaduj paczki, wpisując: library(MASS)i library(car). 3. Zaimportuj daneR , najlepiej za pomocą pliku .txt. 4. Jeżeli dane są nazywane my.datazastosowanie fitdistrw następujący sposób: fitdistr(my.data, distribution="weibull"). 5. Stwórz wykres zgodnie z opisem w pierwszym komentarzu z qqPlot.
COOLSerdash

Odpowiedzi:


28

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.

  1. Pobierz i zainstalujR , to nic nie kosztuje
  2. 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.
  3. Instaluje pakiety MASSi carprzez wpisywanie: install.packages(c("MASS", "car")). Załaduj je, wpisując: library(MASS)i library(car).
  4. 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.
  5. 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.
  6. 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 MLEkc

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

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 n0V1,V1V2,,Vn1Vnf1,f2,,fnp1=f1,p2=f1+f2,,pn=pn1+fny=a+bx

xi=ln(Vi)
yi=ln[ln(1pi)]
ab
c=exp(ab)
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 VmV0.25V0.75 [p(VV0.25)=0.25,p(VV0.75)=0.75]ckc= V m / ln(2 ) 1 / k

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

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


+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?
COOLSerdash
Korzystając z naszej strony potwierdzasz, że przeczytałeś(-aś) i rozumiesz nasze zasady używania plików cookie i zasady ochrony prywatności.
Licensed under cc by-sa 3.0 with attribution required.