Skip to content

Pseudo-random numbers between two numbers

6 messages · guox at ucalgary.ca, (Ted Harding), Duncan Murdoch +1 more

guox at ucalgary.ca
# 
Please forget the last email I sent with the same subject.
=================
I would like to generate pseudo-random numbers between two numbers using
R, up to a given distribution,
for instance, norm. That is something like rnorm(HowMany,Min,Max,mean,sd)
over rnorm(HowMany,mean,sd).
I am wondering if

qnorm(runif(HowMany, pnorm(Min,mean,sd), pnorm(Max,mean, sd)), mean, sd)

is good. Any idea? Thanks.
-james
(Ted Harding)
#
On 10-Mar-09 23:01:45, guox at ucalgary.ca wrote:
It would certainly work as intended! For instance:

truncnorm<-function(HowMany,Min,Max,mean=0,sd=1){
  qnorm(runif(HowMany, pnorm(Min,mean,sd), pnorm(Max,mean, sd)),
        mean, sd)}
Sample <- truncnorm(1000,-1,2.5)
hist(Sample,breaks=100)

Hoping this helps,
Ted.

--------------------------------------------------------------------
E-Mail: (Ted Harding) <Ted.Harding at manchester.ac.uk>
Fax-to-email: +44 (0)870 094 0861
Date: 10-Mar-09                                       Time: 23:39:13
------------------------------ XFMail ------------------------------
(Ted Harding)
# 
I have modified my example to make it more convincing! See at end.
On 10-Mar-09 23:39:17, Ted Harding wrote:
truncnorm<-function(HowMany,Min,Max,mean=0,sd=1){
  qnorm(runif(HowMany, pnorm(Min,mean,sd), pnorm(Max,mean, sd)),
        mean, sd)}
Sample <- truncnorm(100000,-1,2.5)
hist(Sample,breaks=100)
u0<-(-0.975)+0.05*(0:69)
yy<-dnorm(u0)
du<-min(diff(H$breaks))
Y <- 100000*yy*du/(pnorm(2.5)-pnorm(-1.0))
lines(u0,Y)

Ted.

--------------------------------------------------------------------
E-Mail: (Ted Harding) <Ted.Harding at manchester.ac.uk>
Fax-to-email: +44 (0)870 094 0861
Date: 11-Mar-09                                       Time: 00:05:23
------------------------------ XFMail ------------------------------
Duncan Murdoch
#
guox at ucalgary.ca wrote:
It depends on what Min and Max are.  If you get far out in the tails, 
rounding error will kill you.  For example, pnorm(x) is exactly 1 for x 
bigger than 10 or so, so this approach would fail if Min and Max were 
both bigger than 10.  The solution is to switch to lower=FALSE in the 
upper tail, and possibly switch to a log scale if you want to be really 
extreme.

Duncan Murdoch
(Ted Harding)
# 
Because of a mistake I made in copying code into email, there has
been confusion (expressed in some private emails). Hence the
corrected version is appended at the end.
On 11-Mar-09 00:05:26, Ted Harding wrote:
Corrected version:
------------------

truncnorm<-function(HowMany,Min,Max,mean=0,sd=1){
  qnorm(runif(HowMany, pnorm(Min,mean,sd), pnorm(Max,mean, sd)),
        mean, sd)}
Sample <- truncnorm(100000,-1,2.5)
H <- hist(Sample,breaks=100)
u0<-(-0.975)+0.05*(0:69)
yy<-dnorm(u0)
du<-min(diff(H$breaks))
Y <- 100000*yy*du/(pnorm(2.5)-pnorm(-1.0))
lines(u0,Y)


Apologies for the error.
Ted.

--------------------------------------------------------------------
E-Mail: (Ted Harding) <Ted.Harding at manchester.ac.uk>
Fax-to-email: +44 (0)870 094 0861
Date: 11-Mar-09                                       Time: 19:48:12
------------------------------ XFMail ------------------------------
Berwin A Turlach
# 
On Wed, 11 Mar 2009 06:57:02 -0400
Duncan Murdoch <murdoch at stats.uwo.ca> wrote:

            

      
      
      
    

        
This is the approach taken in 
http://www.jstatsoft.org/v16/c02
http://www.jstatsoft.org/v16/c02/paper

which you may want to consult.  They give general routines to truncate
arbitrary distributions.

            

      
      
      
    

        
I agree with this assessment and guess that similar problem will arise
with other truncated distributions created with the code from the paper
above.

And, note, that if sd is small, you can easily be in the situation
where you evaluate pnorm() at extreme values.

            

      
      
      
    

        
Or use package truncnorm.  Though their c code seem to sample with
rejection from a normal proposal which would be quite inefficient for
"extreme" truncations, e.g., standard normal truncated to [8,10].  

Or use package msm, whose rtnorm() implements a more efficient
algorithm for simulating from a truncated normal distribution.

Cheers,

	Berwin