#  Some random number generation in Splus

# generate 100 standard normals
x_rnorm(100)

# generate 100 standard normals constrained to lie
# above xlo and below xhi

# a function that uses rejection

rnormtrunc_function(n,xlo,xhi){
  x_rep(NA,n)
  for(i in 1:n){
    can_rnorm(1)
    while((can < xlo) | (can > xhi)) can_rnorm(1)
    x[i]_can
  }
  return(x)
}

# an example of its use

x_rnormtrunc(10,0,1)

 
# a function that uses the cdf
rnormtrunc2_function(n,xlo,xhi){
  ulo_pnorm(xlo)
  uhi_pnorm(xhi)
  u_runif(n,ulo,uhi)
  x_qnorm(u)
  return(x)
}

# an example of its use

x_rnormtrunc2(10,0,1)


# Now a function that generates a multivariate normal

# generate multivariate Normals
rmultnorm_function(n,mu,sigma){
   # returns n rows of multinormal mu,sigma random vectors
   # ie: returns a n x length(mu) matrix
  p_length(mu)
  z_matrix(rnorm(n * p),nrow=n)
  ch_chol(sigma,pivot=T)
  piv_attr(ch,"pivot")
  zz_(z%*%ch)
  zzz_0*zz
  zzz[,piv]_zz
  zzz + matrix(mu,nrow=n,ncol=p,byrow=T)
}

# generate 100 copies of a bivariate normal...
mean_c(0,0)
covmat_matrix(c(1,.8,.8,1),nrow=2,ncol=2)
z_rmultnorm(100,mean,covmat)
plot(z[,1],z[,2])


# now, a pretty square plot

par(mfrow=c(2,2),pty="s",oma=c(1,1,3,1),mar=c(5,5,2,1))

covmat_matrix(c(1,.3,.3,1),nrow=2,ncol=2)
z_rmultnorm(100,mean,covmat)
plot(z[,1],z[,2],xlab="1st component",ylab="second component",
     xlim=c(-3,3),ylim=c(-3,3),pch=1)
mtext("corr = .3",side=3,line=0,cex=1.0,outer=F)

covmat_matrix(c(1,.6,.6,1),nrow=2,ncol=2)
z_rmultnorm(100,mean,covmat)
plot(z[,1],z[,2],xlab="1st component",ylab="second component",
     xlim=c(-3,3),ylim=c(-3,3),pch=1)
mtext("corr = .6",side=3,line=0,cex=1.0,outer=F)

covmat_matrix(c(1,.9,.9,1),nrow=2,ncol=2)
z_rmultnorm(100,mean,covmat)
plot(z[,1],z[,2],xlab="1st component",ylab="second component",
     xlim=c(-3,3),ylim=c(-3,3),pch=1)
mtext("corr = .9",side=3,line=0,cex=1.0,outer=F)

covmat_matrix(c(1,.98,.98,1),nrow=2,ncol=2)
z_rmultnorm(100,mean,covmat)
plot(z[,1],z[,2],xlab="1st component",ylab="second component",
     xlim=c(-3,3),ylim=c(-3,3),pch=1)
mtext("corr = .98",side=3,line=0,cex=1.0,outer=F)

mtext("Bivariate normal realization with various correlations",
      side=3,line=0,cex=2.0,outer=T)



################################################
################################################

# now a function that generates from an arbitrary
# discrete, finite distribution - use the Splus function sample()

# for example, 10 draws from the distribution over the values
# {1,2,3,4} given by p={.1,.5,.3,.1} can be obtained with the
# commands below

values_c(1,2,3,4)
p_c(.1,.5,.3,.1)

x_sample(values,size=10,replace=T,prob=p)