Hello,

I have to simulate random deviates within a loop with varying parameters. Unfortunately, I don't see how to do it because parameters are declared with the seed, which means I cannot get a proper random sequence because I would declare the seed several times in the program.

For instance, with an exponential random deviate, the following example is a working one

Int main(){
Int i,n
Doub lambda;

n=100;
lambda=0.4;

Expondev myexp(lambda,seed);
for (i=0;i<n;i++){
cout << myexp.dev()<<"\n";
}
}
The "varying-parameter" version is giving always the same number (because of the seed) :

Int main(){
Int i,n=100;
double lambda;
lambda=0.4;

for (i=0;i<n;i++){
lambda=0.4; // or lambda+=0.001 for varying parameter
Expondev myexp(lambda,seed);
cout << myexp.dev()<<"\n";
}
}

I would like to do something like (the code is not correct):

Int main(){
Int i,n=100;
double lambda;
Expondev myexp(seed);

lambda=0.4;
for (i=0;i<n;i++){
lambda+=0.001
cout << myexp.dev(lambda)<<"\n";
}
}

Do you have any idea of how to program this ? Should I re-program the random deviates code or is there any shortcut ?

Thanks

Fred

Hello,

I have to....I don't see how...is there any shortcut ?

The concept (as well as the code) of using uniformly-distributed deviates to generate exponential deviates (In section 7.3.1) is pretty simple, isn't it? (See Footnote.)

I mean, three lines of code in Expondev::dev() and you want a shortcut?

Why not just make your own version of Expondev::dev() that uses a Ran object (no need to make a new derived class or a new "dev" function for Expondev that takes a parameter unless you just want to):


// Put your #include stuff here

//
// Use deviate from uniformly distributed random variable
// on [0,1] to obtain deviate from exponential distribution
//
Doub expdev(Ran & ran, const Doub & beta)
{
//
// Put the three lines of code from Expondev::dev() here,
// except replace
//
// u = doub();
// with
// u = ran.doub();
//
}


Int main()
{
Int i;
Int n = 10;
double lambda = 0.4;
int seed = 87654321;
//int seed = time(0); // To get different sequences in different runs

Ran myran(seed);

cout << fixed;// << scientific;
for (i = 0; i < n; i++) {
lambda += 0.001; // Or whatever you want to do with lambda
cout << setw(14) << expdev(myran, lambda) << endl;
}

return 0;
}


Output (for what it's worth):

0.424774
2.079428
0.419669
4.401218
4.167412
0.902654
0.243658
8.969417
4.163782
0.032607


Regards,

Dave

Footnote:
Oh, yeah, I keep forgetting:
"No one was born knowing this stuff, you know."
---davekw7x

I hope this helps.

Thanks Dave ! I have just tried, and it was fine for the exponential deviate.

I mean, three lines of code in Expondev::dev() and you want a shortcut?

Actually, I'm interested in simulating all kinds of deviates : I would mostly like to use the Binomial and Poisson deviates with parameters varying within a loop : I choose the exponential deviate because it was the simplest ! I will try to make my own functions with a Ran object as you suggested, but in the Binomial case for instance, I'm not so sure of how to do it (I should have asked for the Binomial case before) !

Regards,

Fred

(Inadvertent post has been deleted. If you saw it, then: Oops---never mind. If you didn't see it, then: never mind.)

Here are some actual suggestions.

First of all, I'm sure that you have already looked at the code and the methodology described in the text. If you are happy with the NR3 binomial deviates generator concept, and you have tested the code, then here is what I think you are looking for:

The main thing is that you want to get new deviates with different parameters without creating a new object. Therefore, you will be calling the same Ran object successively to get your new values.

So: here's one possibility:

Modify the Binomialdev class by adding a member function that sort of acts like a constructor, but doesn't instantiate a new Ran.

Instead of changing any of the distribution files (I never do that except for posted and approved bug fixes), here's what I might do:

Copy the entire Binomialdev class from deviates.h to a new file. I might call the file something like binomialdevplus.h

Inside the file, I would change the class name to agree with the file name: BinomialdevPlus
(Don't forget to change the name of the constructor to BinomialdevPlus also.)

Now, copy the entire constructor to another block inside the class definition. This will be the new function. I'll call it ChangeParams

Two notes:

It's not a constructor; it is a member function. So you make a few changes (other than the name):

1. Since it is not a constructor, it has a return type. If you can't think of anything it needs to return, make it void.

2. It only has two parameters: nn and ppp

3. Since it is not a constructor, class members are changed by assignment statements, not base class initializers.

So the file binomialdevplus.h looks like:


// Depends on nr3.h, ran.h gamma.h
struct BinomialdevPlus : Ran {
Doub pp,p,pb,expnp,np,glnp,plog,pclog,sq;
Int n,swch;
Ullong uz,uo,unfin,diff,rltp;
Int pbits[5];
Doub cdf[64];
Doub logfact[1024];

//
// Allow change of parameters without creating new object
// Uses previously-created Ran object
void ChangeParams(Int nn, Doub ppp)
{
pp = ppp;
n = nn;
// Everything else in this function is same
// as the constructor.
.
.
.
} // End of ChangeParameters

// The constructor and dev() function are identical to
// those of Binomaldev
BinomialdevPlus(Int nn, Doub ppp, Ullong i) : Ran(i), pp(ppp), n(nn) {
.
.
.
} // End of constructor
Int dev() {
.
.
.
} // End of dev()
}; // End of BinomialdevPlus class


Then your test program could look something like


// Other includes: nr3.h, ran. gamma.h
#include "./binomialdevplus.h"

Int main()
{
Int i;
Int n = 3;
Int num_samples = 10;
Doub p = 0.5;

int seed = 87654321;
//int seed = time(0); // To get different sequences in different runs

BinomialdevPlus mybin(n, p, seed);

for (i = 0; i < num_samples; i++) {
cout << setw(5) << n << ": " << setw(10) << mybin.dev() << endl;
n *= 2;// or whatever you want to do with n or p
mybin.ChangeParams(n, p);
}

return 0;
}


Output (for what it's worth):

3: 2
6: 4
12: 5
24: 16
48: 21
96: 48
192: 96
384: 201
768: 385
1536: 770


If you need something simpler (that is, faster or more efficient in some respect), then evaluate other generators (or look at the various cases covered by Binomialdev and eliminate the ones that you don't need). Create a "normal" class. Test the scaled-down version first, then add the ChangeParameters function.

Regards,

Dave

Footnote: You could have done this with the exponential example also, but I was giving a simple answer to your simple question. (And, by the way, I meant no disrespect by my wisecrack about a shortcut for three lines of code. The fact is that I never use smilies, and so you couldn't see the tone of my voice or the tongue in my cheek when I wrote that.)

Using a function ChangeParams was indeed great. I have written the following code in a deviatesplus.h file :

struct BinomialdevPlus : Ran {
Doub pp,p,pb,expnp,np,glnp,plog,pclog,sq;
Int n,swch;
Ullong uz,uo,unfin,diff,rltp;
Int pbits[5];
Doub cdf[64];
Doub logfact[1024];

//
// Allow change of parameters without creating new object
// Uses previously-created Ran object
void ChangeParams(Int nn, Doub ppp)
{
pp = ppp;
n = nn;
// Everything else in this function is same
// as the constructor.

} // End of ChangeParameters

BinomialdevPlus(Int nn, Doub ppp, Ullong i) : Ran(i), pp(ppp), n(nn) {
Int j;

// following of the program


I admit knowing very little about object oriented programming (I'm new to C++) so I should probably learn more about how to create /manipulate objects and then go back to the code...

There is one thing I'm not sure to get :
// Everything else in this function is same as the constructor.
Is something else needed ? I have also copied the constructor within the ChangeParams block, but it does not seem to change the results !


Thank you very much anyway...

Fred

Using a function ChangeParams was indeed great. I have written...

I tried to make instructions as explicit as possible without violating copyright or licensing provisions.

I'll try again:

After copying the entire Binomialdev class from the original deviates.h to the new file, here's what the new file looks like, except I have put comments to show the line numbers:

struct Binomialdev : Ran { // line 1 of binomialdevplus.h
.
.
.
Binomialdev(Int nn, Doub ppp, Ullong i) : Ran(i), pp(ppp), n(nn) {// line 8 is start of constructor
.
.
}// line 32 is end of constructor
Int dev() { // line 33 is start of dev()
.
.
.
} // Line 80 is end of dev()
}; // Last line is line 81


Now, make the changes I listed:

Change the name of the class and of its constructor:

struct BinomialdevPlus : Ran { // line 1 of binomialdevplus.h
.
.
.
BinomialdevPlus(Int nn, Doub ppp, Ullong i) : Ran(i), pp(ppp), n(nn) {// line 8 is start of BinomialdevPlus constructor
.
.
.
}// line 32 is end of constructor
Int dev() { // line 33 is start of dev()

Now copy the block consisting of lines 8-32 (the entire constructor) into some other place within the class definition. Lines 8-32 remain in place, unchanged from now on.

If you copy them to a place just before the end of the class definition you will now have the following. I have adjusted comments to show the exact line numbers that I have in my file:

.
.
.
} // Line 80 is end of dev()
BinomialdevPlus(Int nn, Doub ppp, Ullong i) : Ran(i), pp(ppp), n(nn) {// line 81 will be changed in the next step
.
.
.
}// line 105
}; // Last line in file is line 106: End of class definition


Now change the new line 81 to the following:

void ChangeParams(Int nn, Doub ppp) {//New line 81: Start of ChangeParams function

Insert the following two lines after line 81:

pp = ppp; // new line 82
n = nn; // new line 83


Now, the last line of this function is:

}// New line 107: End of ChangeParams


And the last line in the file is

}; // Last line is line 108


Try the following test program just for ChangeParams:


// Other includes: nr3.h, ran.h, gamma.h
#include "./binomialdevplus.h"

Int main()
{
Int n = 3;
Doub p = 0.5;
Int seed = 87654321;

BinomialdevPlus mybin(n, p, seed);

cout << "After constructor:" << endl;
cout << " mybin.n = " << mybin.n << endl;
cout << " mybin.pp = " << mybin.pp << endl << endl;

mybin.ChangeParams(15, 0.99);
cout << "After ChangeParams:" << endl;
cout << " mybin.n = " << mybin.n << endl;
cout << " mybin.pp = " << mybin.pp << endl << endl;

return 0;
}


Output should be:

After constructor:
mybin.n = 3
mybin.pp = 0.5

After ChangeParams:
mybin.n = 15
mybin.pp = 0.99



Regards,

Dave

Thank you Dave !

Let me chime in that we (NR authors) recognize that putting parameters only in the constructors, and being derived classes from Ran, was probably not the best design decision for all the deviates classes. So, repackaging the same underlying methods in a more general design is definitely a good thing to do.

Eventually we'll probably re-do all the deviates classes to fix this, but it may be a while!

Cheers,
Bill P.