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 zainstaluj
R
, 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
MASS
i car
przez wpisywanie: install.packages(c("MASS", "car"))
. Załaduj je, wpisując: library(MASS)
i library(car)
.
- Zaimportuj swoje dane do
R
. Jeśli masz na przykład swoje dane w programie Excel, zapisz je jako plik tekstowy z ogranicznikami (.txt) i zaimportuj za R
pomocą read.table
.
- Użyj funkcji
fitdistr
do 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 R
jest 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")
# 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)
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 MLEkdo
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
doc = 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
n0 - V1, V1- V2), … , Vn - 1- Vnfa1, f2), … , Fnp1= f1, p2)= f1+ f2), … , Sn= pn - 1+ fny= a + b x
xja= ln( Vja)
yja= ln[ - ln( 1 - pja) ]
zabc = 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
V.mV.0,25V.0,75 [p(V≤V0.25)=0.25,p(V≤V0.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 R
poró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 R
skrypt:
#-----------------------------------------------------------------------------
# 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")
#-----------------------------------------------------------------------------
# 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")
fitdistr(mydata, densfun="weibull")
w,R
aby znaleźć parametry za pośrednictwem MLE. Aby utworzyć wykres, użyjqqPlot
funkcji zcar
pakietu:qqPlot(mydata, distribution="weibull", shape=, scale=)
z parametrami kształtu i skali, które znalazłeśfitdistr
.