summaryrefslogtreecommitdiff
path: root/numpy/oldnumeric
diff options
context:
space:
mode:
Diffstat (limited to 'numpy/oldnumeric')
-rw-r--r--numpy/oldnumeric/__init__.py22
-rw-r--r--numpy/oldnumeric/compat.py2
-rw-r--r--numpy/oldnumeric/linear_algebra.py84
-rw-r--r--numpy/oldnumeric/mlab.py72
-rw-r--r--numpy/oldnumeric/olddefaults.py39
-rw-r--r--numpy/oldnumeric/random_array.py268
-rw-r--r--numpy/oldnumeric/rng.py137
-rw-r--r--numpy/oldnumeric/rng_stats.py35
8 files changed, 617 insertions, 42 deletions
diff --git a/numpy/oldnumeric/__init__.py b/numpy/oldnumeric/__init__.py
index 67874f814..14895910f 100644
--- a/numpy/oldnumeric/__init__.py
+++ b/numpy/oldnumeric/__init__.py
@@ -1,33 +1,25 @@
+# Don't add these to the __all__ variable
from numpy import *
+
+# Add these
from compat import *
from olddefaults import *
-from typeconv import *
from functions import *
-import numpy
import compat
import olddefaults
-import typeconv
import functions
+import numpy
__version__ = numpy.__version__
+del numpy
__all__ = ['__version__']
-__all__ += numpy.__all__
__all__ += compat.__all__
-__all__ += typeconv.__all__
-for name in olddefaults.__all__:
- if name not in __all__:
- __all__.append(name)
-
-for name in functions.__all__:
- if name not in __all__:
- __all__.apend(name)
+__all__ += olddefaults.__all__
+__all__ += functions.__all__
-del name
-del numpy
del compat
del olddefaults
del functions
-del typeconv
diff --git a/numpy/oldnumeric/compat.py b/numpy/oldnumeric/compat.py
index 3c8cf9b70..fa6ea94c4 100644
--- a/numpy/oldnumeric/compat.py
+++ b/numpy/oldnumeric/compat.py
@@ -54,7 +54,7 @@ LittleEndian = (sys.byteorder == 'little')
# backward compatible names from old Precision.py
-Character = 'S1'
+Character = 'c'
UnsignedInt8 = _dt_(nt.uint8)
UInt8 = UnsignedInt8
UnsignedInt16 = _dt_(nt.uint16)
diff --git a/numpy/oldnumeric/linear_algebra.py b/numpy/oldnumeric/linear_algebra.py
new file mode 100644
index 000000000..90d6f7138
--- /dev/null
+++ b/numpy/oldnumeric/linear_algebra.py
@@ -0,0 +1,84 @@
+
+"""Backward compatible with LinearAlgebra from Numeric
+"""
+# This module is a lite version of the linalg.py module in SciPy which contains
+# high-level Python interface to the LAPACK library. The lite version
+# only accesses the following LAPACK functions: dgesv, zgesv, dgeev,
+# zgeev, dgesdd, zgesdd, dgelsd, zgelsd, dsyevd, zheevd, dgetrf, dpotrf.
+
+
+__all__ = ['LinAlgError', 'solve_linear_equations',
+ 'inverse', 'cholesky_decomposition', 'eigenvalues',
+ 'Heigenvalues', 'generalized_inverse',
+ 'determinant', 'singular_value_decomposition',
+ 'eigenvectors', 'Heigenvectors',
+ 'linear_least_squares'
+ ]
+
+from numpy.core import transpose
+import numpy.linalg as linalg
+
+# Linear equations
+
+LinAlgError = linalg.LinAlgError
+
+def solve_linear_equations(a, b):
+ return linalg.solve(a,b)
+
+# Matrix inversion
+
+def inverse(a):
+ return linalg.inv(a)
+
+# Cholesky decomposition
+
+def cholesky_decomposition(a):
+ return linalg.cholesky(a)
+
+# Eigenvalues
+
+def eigenvalues(a):
+ return linalg.eigvals(a)
+
+def Heigenvalues(a, UPLO='L'):
+ return linalg.eigvalsh(a,UPLO)
+
+# Eigenvectors
+
+def eigenvectors(A):
+ w, v = linalg.eig(A)
+ return w, transpose(v)
+
+def Heigenvectors(A):
+ w, v = linalg.eigh(A)
+ return w, transpose(v)
+
+# Generalized inverse
+
+def generalized_inverse(a, rcond = 1.e-10):
+ return linalg.pinv(a, rcond)
+
+# Determinant
+
+def determinant(a):
+ return linalg.det(a)
+
+# Linear Least Squares
+
+def linear_least_squares(a, b, rcond=1.e-10):
+ """returns x,resids,rank,s
+where x minimizes 2-norm(|b - Ax|)
+ resids is the sum square residuals
+ rank is the rank of A
+ s is the rank of the singular values of A in descending order
+
+If b is a matrix then x is also a matrix with corresponding columns.
+If the rank of A is less than the number of columns of A or greater than
+the number of rows, then residuals will be returned as an empty array
+otherwise resids = sum((b-dot(A,x)**2).
+Singular values less than s[0]*rcond are treated as zero.
+"""
+ return linalg.lstsq(a,b,rcond)
+
+def singular_value_decomposition(A, full_matrices=0):
+ return linalg.svd(A, full_matrices)
diff --git a/numpy/oldnumeric/mlab.py b/numpy/oldnumeric/mlab.py
new file mode 100644
index 000000000..409238e57
--- /dev/null
+++ b/numpy/oldnumeric/mlab.py
@@ -0,0 +1,72 @@
+# This module is for compatibility only. All functions are defined elsewhere.
+
+from numpy.oldnumeric import *
+
+__all__ = numpy.oldnumeric.__all__
+
+__all__ += ['rand', 'tril', 'trapz', 'hanning', 'rot90', 'triu', 'diff', 'angle', 'roots', 'ptp', 'kaiser', 'randn', 'cumprod', 'diag', 'msort', 'LinearAlgebra', 'RandomArray', 'prod', 'std', 'hamming', 'flipud', 'max', 'blackman', 'corrcoef', 'bartlett', 'eye', 'squeeze', 'sinc', 'tri', 'cov', 'svd', 'min', 'median', 'fliplr', 'eig', 'mean']
+
+import linear_algebra as LinearAlgebra
+import random_array as RandomArray
+from numpy import tril, trapz as _Ntrapz, hanning, rot90, triu, diff, \
+ angle, roots, ptp as _Nptp, kaiser, cumprod as _Ncumprod, \
+ diag, msort, prod as _Nprod, std as _Nstd, hamming, flipud, \
+ amax as _Nmax, amin as _Nmin, blackman, bartlett, corrcoef as _Ncorrcoef,\
+ cov as _Ncov, squeeze, sinc, median, fliplr, mean as _Nmean
+
+from numpy.linalg import eig, svd
+from numpy.random import rand, randn
+
+from typeconv import oldtype2dtype as o2d
+
+def eye(N, M=None, k=0, typecode=None):
+ """ eye returns a N-by-M 2-d array where the k-th diagonal is all ones,
+ and everything else is zeros.
+ """
+ dtype = o2d[typecode]
+ if M is None: M = N
+ m = nn.equal(nn.subtract.outer(nn.arange(N), nn.arange(M)),-k)
+ if m.dtype != dtype:
+ return m.astype(dtype)
+
+def tri(N, M=None, k=0, typecode=None):
+ """ returns a N-by-M array where all the diagonals starting from
+ lower left corner up to the k-th are all ones.
+ """
+ dtype = o2d[typecode]
+ if M is None: M = N
+ m = nn.greater_equal(nn.subtract.outer(nn.arange(N), nn.arange(M)),-k)
+ if m.dtype != dtype:
+ return m.astype(dtype)
+
+def trapz(y, x=None, axis=-1):
+ return _Ntrapz(y, x, axis=axis)
+
+def ptp(x, axis=0):
+ return _Nptp(x, axis)
+
+def cumprod(x, axis=0):
+ return _Ncumprod(x, axis)
+
+def max(x, axis=0):
+ return _Nmax(x, axis)
+
+def min(x, axis=0):
+ return _Nmin(x, axis)
+
+def prod(x, axis=0):
+ return _Nprod(x, axis)
+
+def std(x, axis=0):
+ return _Nstd(x, axis)
+
+def mean(x, axis=0):
+ return _Nmean(x, axis)
+
+def cov(m, y=None, rowvar=0, bias=0):
+ return _Ncov(m, y, rowvar, bias)
+
+def corrcoef(x, y=None):
+ return _Ncorrcoef(x,y,0,0)
+
+
diff --git a/numpy/oldnumeric/olddefaults.py b/numpy/oldnumeric/olddefaults.py
index 356f5f00c..384261437 100644
--- a/numpy/oldnumeric/olddefaults.py
+++ b/numpy/oldnumeric/olddefaults.py
@@ -1,45 +1,32 @@
-__all__ = ['ones', 'empty', 'identity', 'zeros', 'eye', 'tri']
+__all__ = ['ones', 'empty', 'identity', 'zeros']
import numpy.core.multiarray as mu
import numpy.core.numeric as nn
-
-def ones(shape, dtype=int, order='C'):
+from typeconv import convtypecode
+
+def ones(shape, typecode='l', savespace=0, dtype=None):
"""ones(shape, dtype=int) returns an array of the given
dimensions which is initialized to all ones.
"""
- a = mu.empty(shape, dtype, order)
+ dtype = convtypecode(typecode,dtype)
+ a = mu.empty(shape, dtype)
a.fill(1)
return a
-def zeros(shape, dtype=int, order='C'):
+def zeros(shape, typecode='l', savespace=0, dtype=None):
"""zeros(shape, dtype=int) returns an array of the given
dimensions which is initialized to all zeros
"""
- return mu.zeros(shape, dtype, order)
+ dtype = convtypecode(typecode,dtype)
+ return mu.zeros(shape, dtype)
-def identity(n,dtype=int):
+def identity(n,typecode='l', dtype=None):
"""identity(n) returns the identity 2-d array of shape n x n.
"""
+ dtype = convtypecode(typecode, dtype)
return nn.identity(n, dtype)
-def eye(N, M=None, k=0, dtype=int):
- """ eye returns a N-by-M 2-d array where the k-th diagonal is all ones,
- and everything else is zeros.
- """
- if M is None: M = N
- m = nn.equal(nn.subtract.outer(nn.arange(N), nn.arange(M)),-k)
- if m.dtype != dtype:
- return m.astype(dtype)
-
-def tri(N, M=None, k=0, dtype=int):
- """ returns a N-by-M array where all the diagonals starting from
- lower left corner up to the k-th are all ones.
- """
- if M is None: M = N
- m = nn.greater_equal(nn.subtract.outer(nn.arange(N), nn.arange(M)),-k)
- if m.dtype != dtype:
- return m.astype(dtype)
-
-def empty(shape, dtype=int, order='C'):
+def empty(shape, typecode='l', dtype=None):
+ dtype = convtypecode(typecode, dtype)
return mu.empty(shape, dtype, order)
diff --git a/numpy/oldnumeric/random_array.py b/numpy/oldnumeric/random_array.py
new file mode 100644
index 000000000..061465d02
--- /dev/null
+++ b/numpy/oldnumeric/random_array.py
@@ -0,0 +1,268 @@
+# Backward compatible module for RandomArray
+
+__all__ = ['ArgumentError','F','beta','binomial','chi_square', 'exponential', 'gamma', 'get_seed',
+ 'mean_var_test', 'multinomial', 'multivariate_normal', 'negative_binomial',
+ 'noncentral_F', 'noncentral_chi_square', 'normal', 'permutation', 'poisson', 'randint',
+ 'random', 'random_integers', 'seed', 'standard_normal', 'uniform']
+
+ArgumentError = ValueError
+
+import numpy.random.mtrand as mt
+import numpy as Numeric
+
+from types import IntType
+
+def seed(x=0, y=0):
+ if (x == 0 or y == 0):
+ mt.seed()
+ else:
+ mt.seed((x,y))
+
+def get_seed():
+ raise NotImplementedError, \
+ "If you want to save the state of the random number generator.\n"\
+ "Then you should use obj = numpy.random.get_state() followed by.\n"\
+ "numpy.random.set_state(obj)."
+
+def random(shape=[]):
+ "random(n) or random([n, m, ...]) returns array of random numbers"
+ if shape == []:
+ shape = None
+ return mt.random_sample(shape)
+
+def uniform(minimum, maximum, shape=[]):
+ """uniform(minimum, maximum, shape=[]) returns array of given shape of random reals
+ in given range"""
+ if shape == []:
+ shape = None
+ return mt.uniform(minimum, maximum, shape)
+
+def randint(minimum, maximum=None, shape=[]):
+ """randint(min, max, shape=[]) = random integers >=min, < max
+ If max not given, random integers >= 0, <min"""
+ if not isinstance(minimum, IntType):
+ raise ArgumentError, "randint requires first argument integer"
+ if maximum is None:
+ maximum = minimum
+ minimum = 0
+ if not isinstance(maximum, IntType):
+ raise ArgumentError, "randint requires second argument integer"
+ a = ((maximum-minimum)* random(shape))
+ if isinstance(a, Numeric.ArrayType):
+ return minimum + a.astype(Numeric.Int)
+ else:
+ return minimum + int(a)
+
+def random_integers(maximum, minimum=1, shape=[]):
+ """random_integers(max, min=1, shape=[]) = random integers in range min-max inclusive"""
+ return randint(minimum, maximum+1, shape)
+
+def permutation(n):
+ "permutation(n) = a permutation of indices range(n)"
+ return mt.permutation(n)
+
+def standard_normal(shape=[]):
+ """standard_normal(n) or standard_normal([n, m, ...]) returns array of
+ random numbers normally distributed with mean 0 and standard
+ deviation 1"""
+ if shape == []:
+ shape = None
+ return mt.standard_normal(shape)
+
+def normal(mean, std, shape=[]):
+ """normal(mean, std, n) or normal(mean, std, [n, m, ...]) returns
+ array of random numbers randomly distributed with specified mean and
+ standard deviation"""
+ if shape == []:
+ shape = None
+ return mt.normal(mean, std, shape)
+
+def multivariate_normal(mean, cov, shape=[]):
+ """multivariate_normal(mean, cov) or multivariate_normal(mean, cov, [m, n, ...])
+ returns an array containing multivariate normally distributed random numbers
+ with specified mean and covariance.
+
+ mean must be a 1 dimensional array. cov must be a square two dimensional
+ array with the same number of rows and columns as mean has elements.
+
+ The first form returns a single 1-D array containing a multivariate
+ normal.
+
+ The second form returns an array of shape (m, n, ..., cov.shape[0]).
+ In this case, output[i,j,...,:] is a 1-D array containing a multivariate
+ normal."""
+ if shape == []:
+ shape = None
+ return mt.multivariate_normal(mean, cov, shape)
+
+def exponential(mean, shape=[]):
+ """exponential(mean, n) or exponential(mean, [n, m, ...]) returns array
+ of random numbers exponentially distributed with specified mean"""
+ if shape == []:
+ shape = None
+ return mt.exponential(mean, shape)
+
+def beta(a, b, shape=[]):
+ """beta(a, b) or beta(a, b, [n, m, ...]) returns array of beta distributed random numbers."""
+ if shape == []:
+ shape = None
+ return mt.beta(a, b, shape)
+
+def gamma(a, r, shape=[]):
+ """gamma(a, r) or gamma(a, r, [n, m, ...]) returns array of gamma distributed random numbers."""
+ if shape == []:
+ shape = None
+ return mt.gamma(a, r, shape)
+
+def F(dfn, dfd, shape=[]):
+ """F(dfn, dfd) or F(dfn, dfd, [n, m, ...]) returns array of F distributed random numbers with dfn degrees of freedom in the numerator and dfd degrees of freedom in the denominator."""
+ if shape == []:
+ shape == None
+ return mt.f(dfn, dfd, shape)
+
+def noncentral_F(dfn, dfd, nconc, shape=[]):
+ """noncentral_F(dfn, dfd, nonc) or noncentral_F(dfn, dfd, nonc, [n, m, ...]) returns array of noncentral F distributed random numbers with dfn degrees of freedom in the numerator and dfd degrees of freedom in the denominator, and noncentrality parameter nconc."""
+ if shape == []:
+ shape = None
+ return mt.noncentral_f(dfn, dfd, nconc, shape)
+
+def chi_square(df, shape=[]):
+ """chi_square(df) or chi_square(df, [n, m, ...]) returns array of chi squared distributed random numbers with df degrees of freedom."""
+ if shape == []:
+ shape = None
+ return mt.chisquare(df, shape)
+
+def noncentral_chi_square(df, nconc, shape=[]):
+ """noncentral_chi_square(df, nconc) or chi_square(df, nconc, [n, m, ...]) returns array of noncentral chi squared distributed random numbers with df degrees of freedom and noncentrality parameter."""
+ if shape == []:
+ shape = None
+ return mt.noncentral_chisquare(df, nconc, shape)
+
+def binomial(trials, p, shape=[]):
+ """binomial(trials, p) or binomial(trials, p, [n, m, ...]) returns array of binomially distributed random integers.
+
+ trials is the number of trials in the binomial distribution.
+ p is the probability of an event in each trial of the binomial distribution."""
+ if shape == []:
+ shape = None
+ return mt.binomial(trials, p, shape)
+
+def negative_binomial(trials, p, shape=[]):
+ """negative_binomial(trials, p) or negative_binomial(trials, p, [n, m, ...]) returns
+ array of negative binomially distributed random integers.
+
+ trials is the number of trials in the negative binomial distribution.
+ p is the probability of an event in each trial of the negative binomial distribution."""
+ if shape == []:
+ shape = None
+ return mt.negative_binomial(trials, p, shape)
+
+def multinomial(trials, probs, shape=[]):
+ """multinomial(trials, probs) or multinomial(trials, probs, [n, m, ...]) returns
+ array of multinomial distributed integer vectors.
+
+ trials is the number of trials in each multinomial distribution.
+ probs is a one dimensional array. There are len(prob)+1 events.
+ prob[i] is the probability of the i-th event, 0<=i<len(prob).
+ The probability of event len(prob) is 1.-Numeric.sum(prob).
+
+ The first form returns a single 1-D array containing one multinomially
+ distributed vector.
+
+ The second form returns an array of shape (m, n, ..., len(probs)).
+ In this case, output[i,j,...,:] is a 1-D array containing a multinomially
+ distributed integer 1-D array."""
+ if shape == []:
+ shape = None
+ return mt.multinomial(trials, probs, shape)
+
+def poisson(mean, shape=[]):
+ """poisson(mean) or poisson(mean, [n, m, ...]) returns array of poisson
+ distributed random integers with specified mean."""
+ if shape == []:
+ shape = None
+ return mt.poisson(mean, shape)
+
+
+def mean_var_test(x, type, mean, var, skew=[]):
+ n = len(x) * 1.0
+ x_mean = Numeric.sum(x)/n
+ x_minus_mean = x - x_mean
+ x_var = Numeric.sum(x_minus_mean*x_minus_mean)/(n-1.0)
+ print "\nAverage of ", len(x), type
+ print "(should be about ", mean, "):", x_mean
+ print "Variance of those random numbers (should be about ", var, "):", x_var
+ if skew != []:
+ x_skew = (Numeric.sum(x_minus_mean*x_minus_mean*x_minus_mean)/9998.)/x_var**(3./2.)
+ print "Skewness of those random numbers (should be about ", skew, "):", x_skew
+
+def test():
+ obj = mt.get_state()
+ mt.set_state(obj)
+ obj2 = mt.get_state()
+ if (obj2[1] - obj[1]).any():
+ raise SystemExit, "Failed seed test."
+ print "First random number is", random()
+ print "Average of 10000 random numbers is", Numeric.sum(random(10000))/10000.
+ x = random([10,1000])
+ if len(x.shape) != 2 or x.shape[0] != 10 or x.shape[1] != 1000:
+ raise SystemExit, "random returned wrong shape"
+ x.shape = (10000,)
+ print "Average of 100 by 100 random numbers is", Numeric.sum(x)/10000.
+ y = uniform(0.5,0.6, (1000,10))
+ if len(y.shape) !=2 or y.shape[0] != 1000 or y.shape[1] != 10:
+ raise SystemExit, "uniform returned wrong shape"
+ y.shape = (10000,)
+ if Numeric.minimum.reduce(y) <= 0.5 or Numeric.maximum.reduce(y) >= 0.6:
+ raise SystemExit, "uniform returned out of desired range"
+ print "randint(1, 10, shape=[50])"
+ print randint(1, 10, shape=[50])
+ print "permutation(10)", permutation(10)
+ print "randint(3,9)", randint(3,9)
+ print "random_integers(10, shape=[20])"
+ print random_integers(10, shape=[20])
+ s = 3.0
+ x = normal(2.0, s, [10, 1000])
+ if len(x.shape) != 2 or x.shape[0] != 10 or x.shape[1] != 1000:
+ raise SystemExit, "standard_normal returned wrong shape"
+ x.shape = (10000,)
+ mean_var_test(x, "normally distributed numbers with mean 2 and variance %f"%(s**2,), 2, s**2, 0)
+ x = exponential(3, 10000)
+ mean_var_test(x, "random numbers exponentially distributed with mean %f"%(s,), s, s**2, 2)
+ x = multivariate_normal(Numeric.array([10,20]), Numeric.array(([1,2],[2,4])))
+ print "\nA multivariate normal", x
+ if x.shape != (2,): raise SystemExit, "multivariate_normal returned wrong shape"
+ x = multivariate_normal(Numeric.array([10,20]), Numeric.array([[1,2],[2,4]]), [4,3])
+ print "A 4x3x2 array containing multivariate normals"
+ print x
+ if x.shape != (4,3,2): raise SystemExit, "multivariate_normal returned wrong shape"
+ x = multivariate_normal(Numeric.array([-100,0,100]), Numeric.array([[3,2,1],[2,2,1],[1,1,1]]), 10000)
+ x_mean = Numeric.sum(x)/10000.
+ print "Average of 10000 multivariate normals with mean [-100,0,100]"
+ print x_mean
+ x_minus_mean = x - x_mean
+ print "Estimated covariance of 10000 multivariate normals with covariance [[3,2,1],[2,2,1],[1,1,1]]"
+ print Numeric.dot(Numeric.transpose(x_minus_mean),x_minus_mean)/9999.
+ x = beta(5.0, 10.0, 10000)
+ mean_var_test(x, "beta(5.,10.) random numbers", 0.333, 0.014)
+ x = gamma(.01, 2., 10000)
+ mean_var_test(x, "gamma(.01,2.) random numbers", 2*100, 2*100*100)
+ x = chi_square(11., 10000)
+ mean_var_test(x, "chi squared random numbers with 11 degrees of freedom", 11, 22, 2*Numeric.sqrt(2./11.))
+ x = F(5., 10., 10000)
+ mean_var_test(x, "F random numbers with 5 and 10 degrees of freedom", 1.25, 1.35)
+ x = poisson(50., 10000)
+ mean_var_test(x, "poisson random numbers with mean 50", 50, 50, 0.14)
+ print "\nEach element is the result of 16 binomial trials with probability 0.5:"
+ print binomial(16, 0.5, 16)
+ print "\nEach element is the result of 16 negative binomial trials with probability 0.5:"
+ print negative_binomial(16, 0.5, [16,])
+ print "\nEach row is the result of 16 multinomial trials with probabilities [0.1, 0.5, 0.1 0.3]:"
+ x = multinomial(16, [0.1, 0.5, 0.1], 8)
+ print x
+ print "Mean = ", Numeric.sum(x)/8.
+
+if __name__ == '__main__':
+ test()
+
+
diff --git a/numpy/oldnumeric/rng.py b/numpy/oldnumeric/rng.py
new file mode 100644
index 000000000..77a90c694
--- /dev/null
+++ b/numpy/oldnumeric/rng.py
@@ -0,0 +1,137 @@
+# This module re-creates the RNG interface from Numeric
+# Replace import RNG with import numpy.oldnumeric.rng as RNG
+#
+# It is for backwards compatibility only.
+
+
+__all__ = ['CreateGenerator','ExponentialDistribution','LogNormalDistribution','NormalDistribution',
+ 'UniformDistribution', 'error', 'default_distribution', 'random_sample', 'ranf',
+ 'standard_generator']
+
+import numpy.random.mtrand as mt
+import math
+
+class error(Exception):
+ pass
+
+class Distribution(object):
+ def __init__(self, meth, *args):
+ self._meth = meth
+ self._args = args
+
+ def density(self,x):
+ raise NotImplementedError
+
+ def __call__(self, x):
+ return self.density(x)
+
+ def _onesample(self, rng):
+ return getattr(rng, self._meth)(*self._args)
+
+ def _sample(self, rng, n):
+ kwds = {'size' : n}
+ return getattr(rng, self._meth)(*self._args, **kwds)
+
+
+class ExponentialDistribution(Distribution):
+ def __init__(self, lambda_):
+ if (lambda_ <= 0):
+ raise error, "parameter must be positive"
+ Distribution.__init__(self, 'exponential', lambda_)
+
+ def density(x):
+ if x < 0:
+ return 0.0
+ else:
+ lambda_ = self._args[0]
+ return lambda_*exp(-lambda_*x)
+
+class LogNormalDistribution(Distribution):
+ def __init__(self, m, s):
+ m = float(m)
+ s = float(s)
+ if (s <= 0):
+ raise error, "standard deviation must be positive"
+ Distribution.__init__(self, 'lognormal', m, s)
+ sn = math.log(1.0+s*s/(m*m));
+ self._mn = math.log(m)-0.5*sn
+ self._sn = math.sqrt(sn)
+ self._fac = 1.0/math.sqrt(2*math.pi)/self._sn
+
+ def density(x):
+ m,s = self._args
+ y = (math.log(x)-self._mn)/self._sn
+ return self._fac*exp(-0.5*y*y)/x
+
+
+class NormalDistribution(Distribution):
+ def __init__(self, m, s):
+ m = float(m)
+ s = float(s)
+ if (s <= 0):
+ raise error, "standard deviation must be positive"
+ Distribution.__init__(self, 'normal', m, s)
+ self._fac = 1.0/math.sqrt(2*math.pi)/s
+
+ def density(x):
+ m,s = self._args
+ y = (x-m)/s
+ return self._fac*exp(-0.5*y*y)
+
+class UniformDistribution(Distribution):
+ def __init__(self, a, b):
+ a = float(a)
+ b = float(b)
+ width = b-a
+ if (width <=0):
+ raise error, "width of uniform distribution must be > 0"
+ Distribution.__init__(self, 'uniform', a, b)
+ self._fac = 1.0/width
+
+ def density(x):
+ a, b = self._args
+ if (x < a) or (x >= b):
+ return 0.0
+ else:
+ return self._fac
+
+default_distribution = UniformDistribution(0.0,1.0)
+
+class CreateGenerator(object):
+ def __init__(self, seed, dist=None):
+ if seed <= 0:
+ self._rng = mt.RandomState()
+ elif seed > 0:
+ self._rng = mt.RandomState(seed)
+ if dist is None:
+ dist = default_distribution
+ if not isinstance(dist, Distribution):
+ raise error, "Not a distribution object"
+ self._dist = dist
+
+ def ranf(self):
+ return self._dist._onesample(self._rng)
+
+ def sample(self, n):
+ return self._dist._sample(self._rng, n)
+
+
+standard_generator = CreateGenerator(-1)
+
+def ranf():
+ "ranf() = a random number from the standard generator."
+ return standard_generator.ranf()
+
+def random_sample(*n):
+ """random_sample(n) = array of n random numbers;
+
+ random_sample(n1, n2, ...)= random array of shape (n1, n2, ..)"""
+
+ if not n:
+ return standard_generator.ranf()
+ m = 1
+ for i in n:
+ m = m * i
+ return standard_generator.sample(m).reshape(*n)
+
+
diff --git a/numpy/oldnumeric/rng_stats.py b/numpy/oldnumeric/rng_stats.py
new file mode 100644
index 000000000..60248247e
--- /dev/null
+++ b/numpy/oldnumeric/rng_stats.py
@@ -0,0 +1,35 @@
+
+__all__ = ['average', 'histogram', 'standardDeviation', 'variance']
+
+import numpy.oldnumeric as Numeric
+
+def average(data):
+ data = Numeric.array(data)
+ return Numeric.add.reduce(data)/len(data)
+
+def variance(data):
+ data = Numeric.array(data)
+ return Numeric.add.reduce((data-average(data))**2)/(len(data)-1)
+
+def standardDeviation(data):
+ data = Numeric.array(data)
+ return Numeric.sqrt(variance(data))
+
+def histogram(data, nbins, range = None):
+ data = Numeric.array(data, Numeric.Float)
+ if range is None:
+ min = Numeric.minimum.reduce(data)
+ max = Numeric.maximum.reduce(data)
+ else:
+ min, max = range
+ data = Numeric.repeat(data,
+ Numeric.logical_and(Numeric.less_equal(data, max),
+ Numeric.greater_equal(data,
+ min)))
+ bin_width = (max-min)/nbins
+ data = Numeric.floor((data - min)/bin_width).astype(Numeric.Int)
+ histo = Numeric.add.reduce(Numeric.equal(
+ Numeric.arange(nbins)[:,Numeric.NewAxis], data), -1)
+ histo[-1] = histo[-1] + Numeric.add.reduce(Numeric.equal(nbins, data))
+ bins = min + bin_width*(Numeric.arange(nbins)+0.5)
+ return Numeric.transpose(Numeric.array([bins, histo]))
View Lorry Controller Status