[SciPy-User] quasi random, Halton sequence

josef.pktd@gmai... josef.pktd@gmai...
Mon Jun 17 22:25:15 CDT 2013


I didn't find any quasi-random sequences in python that is BSD compatible.
The question shows up every few years. Is there anything now?


a quick translation from c to python, (to be translated to cython and
to c (going in a circle))

maybe there is something slightly off (e.g. gap in circle)

Josef

-----------------
# -*- coding: utf-8 -*-
"""

Created on Mon Jun 17 22:12:21 2013

Author: Sebastien Paris
Josef Perktold translation from c

http://www.mathworks.com/matlabcentral/fileexchange/17457-quasi-montecarlo-halton-sequence-generator
"""

#void halton(int dim , int nbpts, double *h  , double *p )
#{
#
#	double lognbpts , d , sum;
#
#	int i , j , n , t , b;
#
#	static int P[11] = {2 ,3 ,5 , 7 , 11 , 13 , 17 , 19 ,  23 , 29 , 31};
#
#
#	lognbpts = log(nbpts + 1);
#
#
#	for(i = 0 ; i < dim ; i++)
#
#	{
#
#		b      = P[i];
#
#		n      = (int) ceil(lognbpts/log(b));
#
#
#		for(t = 0 ; t < n ; t++)
#
#		{
#			p[t] = pow(b , -(t + 1) );
#		}
#
#
#		for (j = 0 ; j < nbpts ; j++)
#
#		{
#
#			d        = j + 1;
#
#			sum      = fmod(d , b)*p[0];
#
#
#			for (t = 1 ; t < n ; t++)
#
#			{
#
#				d        = floor(d/b);
#
#				sum     += fmod(d , b)*p[t];
#
#			}
#
#
#			h[j*dim + i] = sum;
#
#		}
#
#	}
#
#}

from math import log, floor, ceil, fmod
import numpy as np

def halton(dim, nbpts):
    h = np.empty(nbpts * dim)
    h.fill(np.nan)
    p = np.empty(nbpts)
    p.fill(np.nan)
    P = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31]
    lognbpts = log(nbpts + 1)
    for i in range(dim):
        b = P[i]
        n = int(ceil(lognbpts / log(b)))
        for t in range(n):
            p[t] = pow(b, -(t + 1) )

        for j in range(nbpts):
            d = j + 1
            sum_ = fmod(d, b) * p[0]
            for t in range(1, n):
                d = floor(d / b)
                sum_ += fmod(d, b) * p[t]

            h[j*dim + i] = sum_

    return h.reshape(nbpts, dim)

x = halton(2, 5000);
#plot(x(1 , :) , x(2 , :) , '+')
print x[:5]
import matplotlib.pyplot as plt
plt.figure()
plt.plot(x[:500, 0], x[:500, 1], '+')
plt.title('uniform-distribution (500)')

plt.figure()
plt.plot(x[:, 0], x[:, 1], '+')
plt.title('uniform-distribution')

from scipy import stats
plt.figure()
xn = stats.norm._ppf(x)
plt.plot(xn[:, 0], xn[:, 1], '+')
plt.title('normal-distribution')

plt.figure()
plt.plot(stats.t._ppf(x[:, 0], 3), stats.t._ppf(x[:, 1], 3), '+')
plt.title('t-distribution')

plt.figure()
x0 = xn[:100]
x0 /= np.sqrt((x0*x0 + 1e-100).sum(1))[:,None]
plt.plot(x0[:, 0] , x0[:, 1], '+')
plt.xlim(-1.1, 1.1)
plt.ylim(-1.1, 1.1)
plt.title('uniform on circle')

plt.show()
-----------------


More information about the SciPy-User mailing list