逆变换法----离散
- 离散均匀分布
直接生成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