rm(list=ls()) a<-6.1 b<-1.8 t<-(1:999999)/1000000 ratio<-dbeta(t,a,b) M<-max(ratio) n<-10000 v<-runif(n,0,1) u<-runif(n,0,1) y<-v[u<= (gamma(a+b)*v^(a-1)*(1-v)^(b-1)/(gamma(a)*gamma(b)))/M] y1<-v[u<=dbeta(v,a,b)/M] hist(y) mean(y) a/(a+b) var(y) a*b/((a+b)^2*(a+b+1)) length(y)