逆变换法生成离散、连续随机数

逆变换法----离散

• 离散均匀分布
  直接生成x

rng.uni=function(n,m){
ceiling(m*runif(n))#相当于生成x，因为x在1-m的范围，先0-1再乘m就可了
}
#n是生成的个数，m是均匀分布取值个数
直接用sample函数更方便
sample.int(m,n,replace=TRUE)
例如：在1-6中有放回抽取size个数
sample.int(6,size=20,replace=TRUE)

• 一般离散分布

#x为取值的集合，n是输出的个数，p是取值的概率
rng.dv=function(x,n=1,p=rep(1/length(x),length(x)){
	m=length(x)
	y=numeric(n)
	fvs=c(0,cumsum(prob[-m]))#F的累加和
	for(i in 1:n){
		j=max(which(fvs<=runif(1)))
#判断fvs小于生成的数中最大的那个值是多少，那么这个值就是对应的x中的数
		y[i]=x[j]
	}
	y
}
num=rng.dv(1:3,100,p=c(0.2,0.5,0.3))
prop.table(table(num))*100
#输出
 1  2  3 
18 53 29

• 几何分布随机数
  根据分布函数与k、p的关系求得X的表达式

rng.geom=function(n,p){
	ceiling(log(runif(n)/log(1-p))
}
#例子
x=rng.geom(100,p=0.1)
mean(x)#用均值来检验，几何分布的期望EX=1/p

• 二项分布随机数
  n次实验中成功的次数
  用到了分布函数F(k)与F(k+1)的关系 也就是用递推关系求

rbinom.v=function(num,size=n,p){
	x=numeric(num)
	for(i in 1:num){
		u=runif(1)
		k=0
		cc=p/(1-p)
		a=(1-p)^size
		F=a
		while(u>F){
			a=a*cc*(size-k)/(k+1)
			F=F+a
			k=k+1
		}
		x[i]=k
	}
	x
}

• 泊松分布随机数
  也是递推

rpois.v=function(m,lambda=1){
  x=numeric(m)
  for(i in 1:m){
    U=runif(1)
    k=0
    a=exp(-lambda)#pk
    F=a
    while(U>F){
      a=a*lambda/(k+1)
      F=F+a
      k=k+1
    }
    x[i]=k
  }
  x
}
x=rpois.v(100,lambda=3)
table(x)
mean(x)

连续随机数
比较简单，分布函数与x的关系，求得x的表达式则可

• 三角形分布
  标准的右三角形分布为beta(2,1)
  p(x)=2x F(x)=x^2
  那么x=根号U

rng.rtri=function(n){
  sqrt(runif(n))
}
x=rng.rtri(1000)
hist(x,freq=FALSE)
curve(2*x,0,1,col='red',add=TRUE)

• 指数分布
  x=-log(1-u)/lamba

rng.exp=function(n,lambda=1){
  -log(runif(n))/lambda  
}
x=rng.exp(100,2)
hist(x,freq=FALSE)
curve(2*exp(-2*x),0,10,col='red',add=TRUE)
mean(x)

随机数

• 标准正态随机数

rng.normbm=function(n,mu,sigma){
  n2=ceiling(n/2)
  nn=2*n2
  amp=sqrt(-2*log(runif(n2)))
  ang=2*pi*runif(n2)
  mu+sigma*c(amp*cos(ang),amp*sin(ang))[1:n]
}
x=rng.normbm(1001,10,2)
str(x)
hist(x,freq=FALSE)
curve(dnorm(x,10,2),4,16,col='red',add=TRUE)
#均值为10，方差为2的正态密度曲线

舍选法
第一种
估计出p的上界M，可以想象为用[a,b]*[0,M]的矩形去框住这个f(x)的图

rng.beta=function(n){
  xvec=numeric(n)
  for(i in 1:n){
    repeat{
      U1=runif(1)
      U2=runif(1)
      X=U1
      Y=135/64*U2
      if(Y<=20*X*(1-X)^3) break#跳出repeat记录在xvec中
      
    }
    xvec[i]=X
  }
  xvec
}
x=rng.beta(1000)
hist(x,freq=FALSE,ylim=c(0,2.5))
curve(dbeta(x,2,4),0,1,col='red',add=TRUE)

• 第二种
  找到一个密度函数图像类似的gx，就是说找到一个函数gx去框住fx