rgbeta=function(n=NULL,mu=NULL,sd=NULL,lims=NULL){

# In general, a beta of arbitrary limits can be described in terms of a standard beta, with
# shape factors "a" and "b", like:
  
# x = B(a,b)*(newlims[2] - newlims[1]) - newlims[1]
# where B(a,b) is a vector of random numbers from a Beta with shape factors a & b
# and x is the vector of transformed values.
  
# For a linear transform like:
  
# x = B(a,b)*c + d
# we have:
# mean(x) = mean(B(a,b))*c + d
# and 
# var(x) = (c^2)*var(B(a,b))
  
# Also, from the Beta distribution, we have:
# mean(B(a,b)) = a/(a+b)
# and
# var(B(a,b)) = (a*b)/(((a+b)^2)*(a+b+1))
# Thus we have two constraints:
  
# mean(x) = (a/(a+b))*c + d
# and
# var(x) = (c^2)*((a*b)/(((a+b)^2)*(a+b+1)))
  
# Because we can rearrange and plug in, we have only 1 constraint:
# Letting lambda = (1 - mean(x)/(c) + d)/(mean(x)/c - d)
# We can solve:
#  b = a*lambda
# And our constraint reduces to:
  
# var(x)/(c^2) = lambda*(a^2)/((a^2)*((lambda+1)^2)*(a+lambda*a+1))
# var(x)/(c^2) - lambda*(a^2)/((a^2)*((lambda+1)^2)*(a+lambda*a+1)) = 0
# Thus we choose a such that this difference is close to zero within some tolerance tol. 
# (Here tol=1E-4. SDs are only given with 2 decimal spaces in the OXTR paper for the demo)

s2=sd^2
lambda = (1-(mu/(lims[2]-lims[1]))+lims[1])/((mu/(lims[2]-lims[1]))-lims[1])

s2rat = s2/((lims[2]-lims[1])^2)

betavar=function(a,b){
  numerator=a*b
  denominator=((a+b)^2)*(a+b+1)
  numerator/denominator
}

alpha=1
beta=1*lambda
d=abs(betavar(alpha,beta)-s2rat)
tol=1E-4
while(d>tol){
  alpha=alpha+0.0001
  beta=lambda*alpha
  d=abs(betavar(alpha,beta)-s2rat)
}

xform=function(z,lims){
  z*(lims[2]-lims[1])+lims[1]
}

rgbeta = vector("list",length=3)
names(rgbeta)=c("r","shape1","shape2")
rgbeta$r=xform(rbeta(n,alpha,beta),lims)
rgbeta$shape1=alpha
rgbeta$shape2=beta
return(rgbeta)
}