Parameter distribusi Weibull

Karena @zaynah memposting di komentar bahwa data dianggap mengikuti distribusi Weibull, saya akan memberikan tutorial singkat tentang cara memperkirakan parameter distribusi tersebut menggunakan MLE (Estimasi kemungkinan maksimum). Ada posting serupa tentang kecepatan angin dan distribusi Weibull di situs.

Unduh dan pasangR , gratis
Opsional: Unduh dan instal RStudio , yang merupakan IDE bagus untuk R yang menyediakan banyak fungsi berguna seperti penyorotan sintaks dan banyak lagi.
Instal paket MASSdan cardengan mengetik: install.packages(c("MASS", "car")). Muat dengan mengetik: library(MASS)dan library(car).
Impor data Anda keR . Jika Anda memiliki data di Excel, misalnya, simpan sebagai file teks terbatas (.txt) dan impor Rdengan read.table.
Gunakan fungsi fitdistruntuk menghitung perkiraan kemungkinan maksimum distribusi weibull Anda: fitdistr(my.data, densfun="weibull", lower = 0). Untuk melihat contoh yang sepenuhnya berhasil, lihat tautan di bagian bawah jawaban.
Buat QQ-Plot untuk membandingkan data Anda dengan distribusi Weibull dengan parameter skala dan bentuk yang diperkirakan pada poin 5: qqPlot(my.data, distribution="weibull", shape=, scale=)

The tutorial dari Vito Ricci pada pas distribusi dengan Radalah titik awal yang baik tentang masalah tersebut. Dan ada banyak posting di situs ini tentang subjek (lihat posting ini juga).

Untuk melihat contoh cara penggunaan yang sepenuhnya berhasil fitdistr, lihat pos ini .

Mari kita lihat contoh di 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)

Perkiraan kemungkinan maksimum dekat dengan yang kami tetapkan secara sewenang-wenang dalam pembuatan angka acak. Mari kita bandingkan data kami menggunakan QQ-Plot dengan distribusi Weibull hipotetis dengan parameter yang telah kami perkirakan fitdistr:

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

QQPlot

Poin-poinnya selaras dengan baik pada garis dan sebagian besar dalam amplop kepercayaan-95%. Kami akan menyimpulkan bahwa data kami kompatibel dengan distribusi Weibull. Ini diharapkan, tentu saja, karena kami telah mengambil sampel nilai-nilai kami dari distribusi Weibull.

$k$ $c$

Makalah ini berisi daftar lima metode untuk memperkirakan parameter distribusi Weibull untuk kecepatan angin. Saya akan menjelaskan tiga dari mereka di sini.

Dari cara dan standar deviasi

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

Kuadrat-terkecil cocok dengan distribusi yang diamati

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

Kecepatan angin median dan kuartil

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

Perbandingan keempat metode tersebut

Berikut adalah contoh dalam Rmembandingkan empat metode:

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

Semua metode menghasilkan hasil yang sangat mirip. Pendekatan kemungkinan maksimum memiliki keuntungan bahwa kesalahan standar dari parameter Weibull diberikan secara langsung.

Menggunakan bootstrap untuk menambah interval kepercayaan titik ke PDF atau CDF

Kita dapat menggunakan bootstrap non-parametrik untuk membangun interval kepercayaan searah di sekitar PDF dan CDF dari perkiraan distribusi Weibull. Ini sebuah Rskrip:

#-----------------------------------------------------------------------------
# 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
sumber

+1, ikhtisar yang bagus. NB, sedikit jalan pintas mungkin untuk menggunakan ? QqPlot w / distribution=weibulldari paket mobil, yang akan sesuai dengan paramaters melalui MLE & membuat plot-qq dalam 1 langkah.

gung - Reinstate Monica

@ung, terima kasih. Saya tidak mengetahui bahwa qqPlot dari carmenghitung parameter MLE secara otomatis. Jika saya menghasilkan variabel acak dengan distribusi weibull ( rweibull) dan menggunakan perintah qqPlot(rw, distribution="weibull")saya mendapatkan pesan kesalahan yang mengatakan bahwa harus memberikan parameter shapedan scaleuntuk qqPlot. Apakah saya melewatkan sesuatu?

COOLSerdash

kesalahanku. Jelas, itu hanya secara otomatis memperkirakan parameter dari beberapa distribusi, dan Weibull bukan salah satunya.

gung - Reinstate Monica

hai, saya menemukan bahwa setelah saya mengimpor mydata ke R, ketika saya melakukan perintah, fitdistr (mydata, densfun = "weibull") ia mengatakan kesalahan messgae bahwa "mydata" tidak ditemukan .. sebenarnya data saya hs telah diimpor ke R. jawaban apa pun akan diterima.

Zay

@zaynah Bisakah Anda mengedit jawaban Anda dan memposting kode Anda yang Anda gunakan untuk mengimpor data. Silakan tambahkan pesan kesalahan juga. Bisakah Anda mengimpor data tanpa kesalahan? Apakah Anda memeriksa apakah data itu diimpor dengan benar?