绘制 weibull 的可能性
Plot the likelihood of weibull
我想绘制一个大小为 1000 的 weibull 样本的似然函数,其形状参数为 theta。我使用了标准化的 weibull,所以比例 lambda 是 1。但是输出是一条水平直线。
n<-1000
lik <- function(theta, x){
K<- length(theta)
n<- length(x)
out<- rep(0,K)
for(k in 1:K){
out[k] <- prod(dweibull(x, shape= theta[k], scale=1))
}
return(out)
}
theta<-seq(0.01, 10, by = 0.01)
x <- rweibull(n, shape= 0.5, scale= 1)
plot(theta, lik(theta, x), type="l", lwd=2)
您所做的确实没有错,但计算机很难计算许多小数的乘积,因此最终可能为零(甚至 0.99^1000
= 4^-5)。因此 log
转换然后 sum
更容易。 (由于log变换是单调递增函数最大化log-likelihood和最大化likelihood是一样的)因此改变
prod(dweibull(x, shape= theta[k], scale=1))
到
sum(dweibull(x, shape= theta[k], scale=1, log=TRUE))
另一个小的变化是在 theta
的合理范围内绘制可能性,以便
你可以看到曲线。
工作代码:
set.seed(1)
n<-1000
lik <- function(theta, x){
K <- length(theta)
n <- length(x)
out <- rep(0,K)
for(k in 1:K){
out[k] <- sum(dweibull(x, shape= theta[k], scale=1, log=TRUE))
}
return(out)
}
popTheta = 0.5
theta = seq(0.01, 1.5, by = 0.01)
x = rweibull(n, shape=popTheta, scale= 1)
plot(theta, lik(theta, x), type="l", lwd=2)
abline(v=popTheta)
theta[which.max( lik(theta, x))]
我想绘制一个大小为 1000 的 weibull 样本的似然函数,其形状参数为 theta。我使用了标准化的 weibull,所以比例 lambda 是 1。但是输出是一条水平直线。
n<-1000
lik <- function(theta, x){
K<- length(theta)
n<- length(x)
out<- rep(0,K)
for(k in 1:K){
out[k] <- prod(dweibull(x, shape= theta[k], scale=1))
}
return(out)
}
theta<-seq(0.01, 10, by = 0.01)
x <- rweibull(n, shape= 0.5, scale= 1)
plot(theta, lik(theta, x), type="l", lwd=2)
您所做的确实没有错,但计算机很难计算许多小数的乘积,因此最终可能为零(甚至 0.99^1000
= 4^-5)。因此 log
转换然后 sum
更容易。 (由于log变换是单调递增函数最大化log-likelihood和最大化likelihood是一样的)因此改变
prod(dweibull(x, shape= theta[k], scale=1))
到
sum(dweibull(x, shape= theta[k], scale=1, log=TRUE))
另一个小的变化是在 theta
的合理范围内绘制可能性,以便
你可以看到曲线。
工作代码:
set.seed(1)
n<-1000
lik <- function(theta, x){
K <- length(theta)
n <- length(x)
out <- rep(0,K)
for(k in 1:K){
out[k] <- sum(dweibull(x, shape= theta[k], scale=1, log=TRUE))
}
return(out)
}
popTheta = 0.5
theta = seq(0.01, 1.5, by = 0.01)
x = rweibull(n, shape=popTheta, scale= 1)
plot(theta, lik(theta, x), type="l", lwd=2)
abline(v=popTheta)
theta[which.max( lik(theta, x))]