Ticket #3356: sage-3356-part6.patch

File sage-3356-part6.patch, 3.6 KB (added by cswiercz, 14 years ago)

Added log-normal random distribution.

  • sage/finance/time_series.pyx

    # HG changeset patch
    # User Chris Swierczewski <cswiercz@gmail.com>
    # Date 1213936275 25200
    # Node ID d823dad321bb50a6b23b9a9f717725b244d065aa
    # Parent  db814611f2fa7407a77598ffeaf93006fce59276
    Added _randomize_lognormal to list of available randomizations for TimeSeries class as well as docstring for other functions.
    
    diff -r db814611f2fa -r d823dad321bb sage/finance/time_series.pyx
    a b cdef class TimeSeries: 
    16591659            self._randomize_normal(loc, scale)
    16601660        elif distribution == 'semicircle':
    16611661            self._randomize_semicircle(loc)
     1662        elif distribution == 'lognormal':
     1663            self._randomize_lognormal(loc, scale)
    16621664        else:
    16631665            raise NotImplementedError
    16641666        return self
    16651667
    16661668    def _randomize_uniform(self, double left, double right):
     1669        """
     1670        Generates a uniform random distribution of doubles between left and
     1671        right and stores values in place.
     1672
     1673        INPUT:
     1674            left -- left bound on random distribution
     1675            right -- right bound on random distribution
     1676
     1677        EXAMPLES:
     1678        """
    16671679        if left >= right:
    16681680            raise ValueError, "left must be less than right"
    16691681
    cdef class TimeSeries: 
    16751687
    16761688    def _randomize_normal(self, double m, double s):
    16771689        """
     1690        Generates a normal random distribution of doubles with mean m and
     1691        standard deviation s and stores values in place.
     1692       
    16781693        INPUT:
    16791694            m -- mean
    16801695            s -- standard deviation
    cdef class TimeSeries: 
    17021717                self._values[k] = m + y2*s
    17031718
    17041719    def _randomize_semicircle(self, double center):
     1720        """
     1721        Generates a semicircle random distribution of doubles about center
     1722        and stores values in place. Uses the acceptance-rejection method.
     1723
     1724        INPUT:
     1725            center -- the center of the semicircle distribution
     1726
     1727        EXAMPLES:
     1728        """
    17051729        cdef Py_ssize_t k
    17061730        cdef double x, y, s, d = 2, left = center - 1, z
    17071731        cdef randstate rstate = current_randstate()
    cdef class TimeSeries: 
    17151739                    break
    17161740            self._values[k] = x + center
    17171741
     1742    def _randomize_lognormal(self, double m, double s):
     1743        """
     1744        Generates a log-normal random distribution of doubles with mean m
     1745        and standard deviation s. Usues Box-Muller algorithm and the identity:
     1746        if Y is a random variable with normal distribution then X = exp(Y)
     1747        is a random variable with log-normal distribution.
    17181748
     1749        INPUT:
     1750            m -- mean
     1751            s -- standard deviation
     1752
     1753        EXAMPLES:
     1754        """
     1755        # Ported from http://users.tkk.fi/~nbeijar/soft/terrain/source_o2/boxmuller.c
     1756        # This the box muller algorithm.
     1757        cdef randstate rstate = current_randstate()
     1758        cdef double x1, x2, w, y1, y2
     1759        cdef Py_ssize_t k
     1760        for k from 0 <= k < self._length:
     1761            while 1:
     1762                x1 = 2*rstate.c_rand_double() - 1
     1763                x2 = 2*rstate.c_rand_double() - 1
     1764                w = x1*x1 + x2*x2
     1765                if w < 1: break
     1766            w = sqrt( (-2*log(w))/w )
     1767            y1 = x1 * w
     1768            y2 = x2 * w
     1769            self._values[k] = exp(m + y1*s)
     1770            k += 1
     1771            if k < self._length:
     1772                self._values[k] = exp(m + y2*s)
     1773       
    17191774    def fft(self, bint overwrite=False):
    17201775        """
    17211776        Return the real discrete fast fourier transform of self, as a