Parâmetros de distribuição Weibull

19

Oi, o mesmo pode ser mostrado para obter o parâmetro de forma e escala para o método de máxima verossimilhança modificado

Zay
fonte
2
Olá @zaynah e bem-vindo ao site. Não tenho certeza se sua pergunta é se seus dados são compatíveis com uma distribuição weibull ou quais seriam os parâmetros de uma distribuição weibull que descreve seus dados. Se você assumir que os dados seguem uma distribuição Weibull e quer encontrar os parâmetros, você pode usar fitdistr(mydata, densfun="weibull")em Rencontrar os parâmetros via MLE. Para criar um gráfico, use a qqPlotfunção do carpacote: qqPlot(mydata, distribution="weibull", shape=, scale=)com os parâmetros de forma e escala encontrados fitdistr.
COOLSerdash
oi obrigado por uma resposta rápida, meus dados são velocidade do vento mensal média por 5 anos, é compatível com weibull. problema é que eu não sei como encontrar k e c, ou seja, os parâmetros do weibull .. e eu não sei como comparar dados experimentais com o weibull ... também o que é MLE ... :(
Zay
MLE = Estimativa de máxima verossimilhança. Não sei qual software você usa, mas no R, que está disponível gratuitamente, você pode instalar e carregar o pacote MASSe usar fitdistrcom seus dados para calcular as estimativas de k e c. E então, você pode comparar seus dados com o weibull com os parâmetros estimados usando qqPloto carpacote.
COOLSerdash
muito obrigado COOlserdash, estou baixando o software R.
Zay 31/05
1
Ok, aqui está um tutorial passo a passo: 1. Baixe e instale R. 2. Instale os pacotes MASSe cardigitando: install.packages(c("MASS", "car")). Carregue os pacotes digitando: library(MASS)e library(car). 3. Importe seus dados paraR , de preferência com um arquivo .txt. 4. Se os dados é chamado my.datao uso fitdistrda seguinte forma: fitdistr(my.data, distribution="weibull"). 5. Faça um gráfico como eu descrevi no primeiro comentário com qqPlot.
COOLSerdash

Respostas:

28

Como o @zaynah postou nos comentários que os dados devem seguir uma distribuição Weibull, vou fornecer um breve tutorial sobre como estimar os parâmetros dessa distribuição usando o MLE (Estimativa de máxima verossimilhança). Há um post semelhante sobre a velocidade do vento e a distribuição Weibull no site.

  1. Baixe e instaleR , é grátis
  2. Opcional: Faça o download e instale o RStudio , que é um ótimo IDE para R, oferecendo várias funções úteis, como destaque de sintaxe e muito mais.
  3. Instalar os pacotes MASSe cardigitando: install.packages(c("MASS", "car")). Carregue-os digitando: library(MASS)e library(car).
  4. Importe seus dados paraR . Se você tiver seus dados no Excel, por exemplo, salve-os como arquivo de texto delimitado (.txt) e importe-os Rcom read.table.
  5. Use a função fitdistrpara calcular as estimativas de máxima verossimilhança de sua distribuição weibull: fitdistr(my.data, densfun="weibull", lower = 0). Para ver um exemplo completo, consulte o link na parte inferior da resposta.
  6. Faça um QQ-Plot para comparar seus dados com uma distribuição Weibull com os parâmetros de escala e forma estimados no ponto 5: qqPlot(my.data, distribution="weibull", shape=, scale=)

O tutorial de Vito Ricci sobre como ajustar a distribuição Ré um bom ponto de partida sobre o assunto. E há várias postagens neste site sobre o assunto (veja também esta postagem ).

Para ver um exemplo completo de como usar fitdistr, dê uma olhada nesta postagem .

Vejamos um exemplo em 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)

As estimativas de máxima verossimilhança são próximas daquelas que definimos arbitrariamente na geração dos números aleatórios. Vamos comparar nossos dados usando um QQ-Plot com uma distribuição hipotética de Weibull com os parâmetros que estimamos fitdistr:

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

QQPlot

Os pontos estão bem alinhados na linha e, principalmente, dentro do envelope de 95% de confiança. Concluímos que nossos dados são compatíveis com uma distribuição Weibull. Isso era esperado, é claro, já que amostramos nossos valores a partir de uma distribuição Weibull.


Estimativa de (forma) ec (escala) de uma distribuição Weibull sem MLEkc

Este artigo lista cinco métodos para estimar os parâmetros de uma distribuição Weibull para velocidades do vento. Eu vou explicar três deles aqui.

De médias e desvio padrão

k

k=(σ^v^)-1.086
c
c=v^Γ(1+1/k)
v^σ^Γ

Mínimos quadrados ajustados à distribuição observada

n0 0-V1,V1-V2,...,Vn-1-Vnf1,f2,...,fnp1=f1,p2=f1+f2,...,pn=pn-1+fny=uma+bx

xEu=em(VEu)
yEu=em[-em(1-pEu)]
umab
c=exp(ab)
k=b

Velocidades medianas e quartis de vento

VmV0.25V0.75 [p(VV0.25)=0.25,p(VV0.75)=0.75]ck

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

Comparison of the four methods

Here is an example in R comparing the four methods:

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

All methods yield very similar results. The maximum likelihood approach has the advantage that the standard errors of the Weibull parameters are directly given.


Using bootstrap to add pointwise confidence intervals to the PDF or CDF

We can use a the non-parametric bootstrap to construct pointwise confidence intervals around the PDF and CDF of the estimated Weibull distribution. Here's an R script:

#-----------------------------------------------------------------------------
# 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
fonte
+1, nice overview. NB, a slight shortcut might be to use ?qqPlot w/ distribution=weibull from the car package, which will fit the paramaters via MLE & make the qq-plot in 1 step.
gung - Reinstate Monica
@gung Thanks. I'm not aware that qqPlot from car calculates the MLE parameters automatically. If I generate a random variable with a weibull distribution (rweibull) and use the command qqPlot(rw, distribution="weibull") I get an error message saying that must provide the parameters shape and scale to qqPlot. Am I missing something?
COOLSerdash
my mistake. Evidently, it only automatically estimates parameters from some distributions, and Weibull isn't one of those.
gung - Reinstate Monica
hi, i found that after i import mydata into R, when i do the command,fitdistr(mydata, densfun="weibull") it says error messgae that "mydata" not found.. in fact my data hs been imported into R. any answer would be welcome.
Zay
@zaynah Could you please edit your answer and post your code that you use to import the data. Please add the error message too. Could you import the data without errors? Did you check if the data was imported correctly?
COOLSerdash