|
如何用R语言生成对数凹函数随机数,救救孩子吧 ,看了好久了
R1=rep(0,15)
R1[15]=5
X=c(0.5794672, 0.7091715, 0.8354615 ,0.9360012 ,0.9377477, 0.9993006 ,1.0504122 ,1.0855902,
1.3160316, 1.4271615, 1.5559699, 1.5595734, 1.6737055)
JX=13
RXJ=7
TimX=X[13]
a=1
b=2
#要生成的随机数服从的密度函数,密度函数只在0~0.5794672 有定义 而且在这个范围是对数凹的
Paimiu <- function(miu){
Gmiu <- sum((R1[1:JX]+1)*(X-miu)^2)+RXJ*(TimX-miu)^2
value = prod(X-miu)/((b+Gmiu)^(JX+a))
# cat(value,'\n')
return(value)
}
#对数凹函数算法(帮忙看看有什么问题,拜托了)
M <- 0.01
C <- Paimiu(M)
r <- log(C)
a=TRUE
while (a==TRUE) {
U <- runif(1,0,2)
E <- rexp(1)
if(U<=1){
miuvalue=U
t=-E
}else{
Ex=rexp(1)
miuvalue=1+Ex
t=-E-Ex
}
miuvalue <- M +miuvalue/C
cat(miuvalue,'\n')
if(t<= log(Paimiu(miuvalue))-r){
a=FALSE
}
}
这是算法论文截图
|