R 语言中积分的一个问题

shenzhonghua

我现在自己构造了一个密度函数，想生成随机数，因为这个密度函数写不出来，所以我想这么做，密度基本都集中在（-5，5）之间，我想生成均匀（-5，5）上的随机数，比如1000个，带入密度函数的积分F（x)，然后和我在unif（0,1)上生成的随机数（比如说是0.4），进行比较，比如大于0.4的有300个，小于0.4的有700个，那么我近似看这个随机数为[-5+（3/10)*10]

问题1：如果说我对于累积分布函数来说，是（0,X)上f(x)的积分，如果我抽（0,1)上的随机数，然后反解x就能够做出随机数，不知道有什么办法能解，积分的方程

问题2：我在把F(x)和0.4比较的时候好像不能比，也就是说积分出来的值，和一般的值还不能比来比去，说是什么不是倍数，不知道有没有高手能帮我想想，怎么把积分的值，做的能和一般的数据做比较

程序

F<-function(x)

{

return(integrate(fh,-Inf,x,stop.on.error='F'))

}

a<-runif(100,-4,4)

b<-numeric(100)

for(i in 1:100)

{

b<-F(a)

}

这里F(a)根本不能赋值到b上，怎么解决？或者怎么直接比较

zltbaggio

楼主，现在已经有不少办法sample from very complicated distribution。你可以把你的pdf说出来看看。

我的疑问是你写不出f，却可以求得出F？

不过如果你有F，有一个直接的方法sample from the corresponding distribution.

1. 你先generate X from unif(0,1)

2. calculate <bblatex>Y=F^{-1}(X)</bblatex>.

Y 就是你要生成的随机数，from f(x)。反复1000次，你就有1000个随机数了。在第二步的时候，你要求CDF的反函数，可能需要用到uniroot函数。你看一下相关的help文件就可以了

shenzhonghua

回复第2楼的 zltbaggio：

好的，其实是个核密度的非参数密度函数估计，你可以认为是很多非中心的dnorm(x）的求和，你如果把这个函数在（-Inf,x）上积分的话，就是cdf了，不知道你有没有办法解决F(x)=runif(0,1)这个方程的办法，我现在知道积分的值他是一个list，F（x)=0x上f(x)的积分，来用这个等式去求解，不大对，比较麻烦，但是如果有具体的数据，就能得到数字值了，比如as.numeric(F(1)[1])就能把积分值转化为数字F(1)本身是个list,需要转化。

zltbaggio

你如果可以把CDF 告诉我，我可以试着写一下R code，再贴给你。

如果你是mixture of normal, 不建议用上述方法。可以尝试用Gibbs sampling。我也可以试着写一下code。

最好把你的问题完整的叙述一下

shenzhonghua

回复第4楼的 zltbaggio：

fh<-function(x)

{

return(((dnorm((x-Sample[1])/0.4895))+(dnorm((x-Sample[2])/0.4895))+(dnorm((x-Sample[3])/0.4895))+(dnorm((x-Sample[4])/0.4895))+(dnorm((x-Sample[5])/0.4895))+(dnorm((x-Sample[6])/0.4895))+(dnorm((x-Sample[7])/0.4895))+(dnorm((x-Sample[8])/0.4895))+(dnorm((x-Sample[9])/0.4895))+(dnorm((x-Sample[10])/0.4895))+(dunif((x-Sample[11])/0.4895))+

(dnorm((x-Sample[12])/0.4895))+(dnorm((x-Sample[13])/0.4895))+(dnorm((x-Sample[14])/0.4895))+(dnorm((x-Sample[15])/0.4895))+(dnorm((x-Sample[16])/0.4895))+(dnorm((x-Sample[17])/0.4895))+(dnorm((x-Sample[18])/0.4895))+(dnorm((x-Sample[19])/0.4895))+(dnorm((x-Sample[20])/0.4895))+(dnorm((x-Sample[21])/0.4895)))/(21*0.4895))

}

F<-function(x)

{

return(integrate(fh,-Inf,x))

}

zltbaggio

回复第5楼的 shenzhonghua：

你就不能给个cdf的数学表达式？从你给的code里面我无法看出你的distribution到底是什么，抱歉