Ask Your Question
1

setting seed for Monte Carlo

asked 2024-04-29 17:35:14 +0100

thethinker gravatar image

updated 2024-05-04 19:10:50 +0100

Does sage seriously not have the ability to change the seed for a Monte Carlo execution? I assumed there is at some way to pass the time, or some other seed generator, to the algorithm.

f(x)=ln(x)
for i in range(10):
    I=monte_carlo_integral(f,[0],[1],100000, algorithm='plain')
    print(I[0])

-1.0002232192476281
-1.0002232192476281
-1.0002232192476281
-1.0002232192476281
-1.0002232192476281
-1.0002232192476281
-1.0002232192476281
-1.0002232192476281
-1.0002232192476281
-1.0002232192476281

EDIT: Based on a helpful comment: randstate is apparently not impacting the Monte Carlo routines.

I=monte_carlo_integral(f,[0],[1],100000, algorithm='plain')
print(I[0])
from sage.misc.randstate import randstate
r = randstate(54321)
print(r.seed())
I=monte_carlo_integral(f,[0],[1],100000, algorithm='plain')
print(I[0])
r = randstate()
print(r.seed())
I=monte_carlo_integral(f,[0],[1],100000, algorithm='plain')
print(I[0])

-1.0002232192476281
54321
-1.0002232192476281
305159320933427511743223216889072415483
-1.0002232192476281

Any other thoughts?

edit retag flag offensive close merge delete

Comments

The random component is initialized inside the code of the function, the relevant lines are:

    try:
        # Initialize the random number generator
        gsl_rng_env_setup()
        type_rng = gsl_rng_default
        _rng = gsl_rng_alloc(type_rng)

So in order to have access to the randomizer, we need to go deep into the cython tunnel.

Why is it important to have a random component?

dan_fulea gravatar imagedan_fulea ( 2024-04-29 23:16:04 +0100 )edit

It's one of the basic methods to determine the statistical errors around a MC calculation. Run N times (with N different seeds) to find an average and standard deviation. You can check for Guassianity also, for example. I can't argue for how complicated implementation would be, but assuming Gaussian statistics limits the relevance of the implementation. Mathematica allows seed changes right at the cli, if that matters. Shall I file a feature request?

thethinker gravatar imagethethinker ( 2024-04-30 17:35:24 +0100 )edit

1 Answer

Sort by ยป oldest newest most voted
0

answered 2024-05-02 17:47:18 +0100

Emmanuel Charpentier gravatar image

seed? may enlighten you :

Docstring:     
WARNING: the enclosing module is marked 'needs sage.groups sage.libs.gap sage.libs.ntl sage.libs.pari',
so doctests may not pass.

   The "randstate" class.  This class keeps track of random number
   states and seeds.  Type "sage.misc.randstate?" for much more
   information on random numbers in Sage.
Init docstring:
WARNING: the enclosing module is marked 'needs sage.groups sage.libs.gap sage.libs.ntl sage.libs.pari',
so doctests may not pass.

   Initialize a new "randstate" object with the given seed (which must
   be coercible to a Python long).

   If no seed is given, then a seed is automatically selected using
   "os.urandom()" if it is available, or the current time otherwise.

   EXAMPLES:

      sage: from sage.misc.randstate import randstate
      sage: r = randstate(54321); r

[ Snip... ]

HTH,

edit flag offensive delete link more

Comments

Good suggestion, but apparently sage is not accessing the randstate routines when doing monte carlo - see my edit.

thethinker gravatar imagethethinker ( 2024-05-04 19:11:34 +0100 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2024-04-29 17:35:14 +0100

Seen: 659 times

Last updated: May 04