Changeset 307

Timestamp:

11/08/08 10:58:36 (2 months ago)

Author:

emuller

Message:

Fixed inh_gamma_generator hazard function numerics
Fixed poisson_process ISI buffer overrun
Added stgen module docstring
Added more tests
Added inh_gamma_psth.py example

Files:

• trunk/examples/stgen (added)
• trunk/examples/stgen/inh_gamma_psth.py (added)
• trunk/src/__init__.py (modified) (2 diffs)
• trunk/src/stgen.py (modified) (13 diffs)
• trunk/test/test_stgen.py (modified) (2 diffs)

trunk/src/__init__.py

@@ -6 +6 @@
-
-# The nice thing would be to gathered every non standard
-z
+e
-# messages and the check
-dependencies = {'pylab' : {'website' : 'http://matplotlib.sourceforge.net/', 'is_present' : False, 'check':False},
@@ -16 +16 @@
-                'NeuroTools.facets.hdf5' : {'website' : None, 'is_present' : False, 'check':False},
-                'srblib' : {'website' : 'http://www.sdsc.edu/srb/index.php/Python', 'is_present' : False, 'check':False},
+                'rpy' : {'website' : 'http://rpy.sourceforge.net/', 'is_present' : False, 'check':False},
+
-                ## Add here your extensions ###
-               }

trunk/src/stgen.py

@@ +1 @@
+"""
+NeuroTools.stgen
+================
+
+A collection of tools for stochastic process generation.
+
+
+Classes
+-------
+
+StGen - Object to generate stochastic processes of various kinds
+        and return them as SpikeTrain or AnalogSignal objects.
+
+
+Functions
+---------
+
+shotnoise_fromspikes - Convolves the provided spike train with shot decaying exponential.
+
+gamma_hazard - Compute the hazard function for a gamma process with parameters a,b.
+"""
+
+
-# TODO : needs refactoring
-# TODO: needs python gsl wrapper / convert to numpy
@@ -9 +32 @@
-
-
-
+_scipy
-    """
-    Compute the hazard function for a gamma process with parameters a,b
@@ -21 +44 @@
-    """
-
-
+
+
-    # Try: 
-    # plot(stgen.gamma_hazard(arange(0,1000.0,0.1),10.0,1.0/50.0))
@@ -39 +63 @@
-    else:
-        return val
+
+
+def gamma_hazard(x, a, b, dt=1e-4):
+    """
+    Compute the hazard function for a gamma process with parameters a,b
+    where a and b are the parameters of the gamma PDF:
+    y(t) = x^(a-1) \exp(-x/b) / (\Gamma(a)*b^a)
+
+    Inputs:
+        x   - in units of seconds
+        a   - dimensionless
+        b   - in units of seconds
+    """
+
+    # Used by inh_gamma_generator
+
+    # Ideally, I would like to see an implementation which does not depend on RPy
+    # but the gamma_hazard_scipy above using scipy exhibits numerical problems, as it does not
+    # support directly returning the log.
+
+    if not check_dependency('rpy'):
+        raise ImportError("gamma_hazard requires RPy (http://rpy.sourceforge.net/)")
+
+    from rpy import r
+
+    # scipy.special.gammaincc has numerical problems
+    #Hpre = -log(scipy.special.gammaincc(a,(x-dt)/b))
+    #Hpost = -log(scipy.special.gammaincc(a,(x+dt)/b))
+
+    # reverting to the good old r.pgamma
+    Hpre = -r.pgamma(x-dt,shape=a,scale=b,lower=False,log=True)
+    Hpost = -r.pgamma(x+dt,shape=a,scale=b,lower=False,log=True)
+    val =  0.5*(Hpost-Hpre)/dt
+
+    return val
+
+
-    
-
@@ -45 +106 @@
-    def __init__(self, rng=None, seed=None):
-        """ 
-
-
-
+
+
+
+
+
+
-
-        Inputs:
-Seed for the random number generator
+The random number generator state object (optional)
-                  a numpy.random.RandomState object, or an object with the same 
-                  interface.
+
+            seed - A seed for the rng (optional).
-
-        If rng is not None, the provided rng will be used to generate random numbers, 
@@ -60 +126 @@
-        Examples:
-            >> x = StGen()
+
+
+
+        StGen Methods:
+
+        Spiking point processes:
+        ------------------------
+ 
+        poisson_generator - homogeneous Poisson process
+        inh_poisson_generator - inhomogeneous Poisson process (time varying rate)
+        inh_gamma_generator - inhomogeneous Gamma process (time varying a,b)
+
+        Continuous time processes:
+        --------------------------
+
+        OU_generator - Ohrnstein-Uhlenbeck process
+        
+
+        See also:
+          shotnoise_fromspikes
+
-        """
-
@@ -75 +162 @@
-
-
-
+,debug=False
-        """
-        Returns a SpikeList whose spikes are a realization of a Poisson process
@@ -98 +185 @@
-        """
-
-
+
+
+
+
+
+
+
+
-        if number > 0:
-            isi = self.rng.exponential(1.0/rate, number)*1000.0
@@ -111 +205 @@
-        i = numpy.searchsorted(spikes, t_stop)
-
+        extra_spikes = []
-        if i==len(spikes):
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
-
@@ -270 +384 @@
-    # TODO: provide optimized C/weave implementation if possible
-
-
+
+
+
+
+
+
+
+
-
-
@@ -432 +553 @@
-    OU_generator = _OU_generator_python2
-
-
+optimized 
-
-
@@ -439 +560 @@
-# TODO fix shotnoise stuff below  ... and write tests
-
-    def shotnoise_generator
-    
-    
-
-    
-    
-    
-    
-
-    
-    
-
-    
-    
-
-    
-    
-    
-
-    
+def shotnoise_fromspikes
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
-
@@ -466 +587 @@
-
-
-exp_conv(poisson_train,t,q,tau,
+shotnoise_fromspikes(spike_train,q,tau,dt,array=False, 
-    """ 
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
-    # time of vanishing significance
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
-

trunk/test/test_stgen.py

@@ -41 +41 @@
-        assert isinstance(st, numpy.ndarray)
-
+        st = stg.poisson_generator(rate,0.0,t_stop,array=True,debug=True)
+        
+        assert isinstance(st[0], numpy.ndarray)
+        assert isinstance(st[1], list)
+
+        st = stg.poisson_generator(rate,0.0,t_stop,debug=True)
+        
+        assert isinstance(st[0], signals.SpikeTrain)
+        assert isinstance(st[1], list)
+
-
-    def testStatsPoisson(self):
+
+        # this is a statistical test with non-zero chance of failure
+
+        def test_poisson(rate,t_start,t_stop):
+            stg = stgen.StGen()
+            dt = t_stop-t_start
+            N = rate*dt/1000.0
+
+            st = stg.poisson_generator(rate,t_start=t_start,t_stop=t_stop,array=True)
+
+            if len(st) in (0,1,2,3):
+                assert N<15
+                return
+
+
+
+            assert st[-1] < t_stop
+            assert st[0] > t_start
+
+
+            # last spike should not be more than 4 ISI away from t_stop
+            err = """
+    Last spike should not be more than 4 ISI behind t_stop.
+    There is a non-zero chance for this to occur during normal operation.
+    Re-run the test to see if the error persists."""
+
+            if st[-1] < t_stop-4.0*1.0/rate*1000.0:
+                raise StatisticalError(err)
+
+
+            # first spike should not be more than 4 ISI away from t_start
+            err = """
+    First spike should not be more than 4 ISI in front of t_start.
+    There is a non-zero chance for this to occur during normal operation.
+    Re-run the test to see if the error persists."""
+
+            if st[0] > t_start+4.0*1.0/rate*1000.0:
+                raise StatisticalError(err)
+
+            err = """
+    Number of spikes should be within 3 standard deviations of mean.
+    There is a non-zero chance for this to occur during normal operation.
+    Re-run the test to see if the error persists."""
+
+
+            if len(st) > N+3.0*numpy.sqrt(N) or len(st) < N-3.0*numpy.sqrt(N):
+                raise StatisticalError(err)
+
+
+        # high rates
+
+        test_poisson(100.0,500.0,1500.0)
+
+        # high rates, short time
+
+        test_poisson(100.0,500.0,550.0)
+
+        # low rates, short time
+
+        test_poisson(2.0,500.0,550.0)
+
+        # low rates, long time
+
+        test_poisson(5.0,500.0,50500.0)
+
+
+
+
+    def testStatsInhPoisson(self):
-
-        # this is a statistical test with non-zero chance of failure
@@ -52 +131 @@
-        t_stop = 1500.0 # milliseconds
-
-        st = stg.poisson_generator(rate,t_start=t_start,t_stop=t_stop,array=True)
-
-        assert st[-1] < t_stop
-        assert st[0] > t_start
-
-
-        # last spike should not be more than 4 ISI away from t_stop
-        err = """
-Last spike should not be more than 4 ISI behind t_stop.
-There is a non-zero chance for this to occur during normal operation.
-Re-run the test to see if the error persists."""
-
-        if st[-1] < t_stop-4.0*1.0/rate*1000.0:
-            raise StatisticalError(err)
-
-
-        # first spike should not be more than 4 ISI away from t_start
-        err = """
-First spike should not be more than 4 ISI in front of t_start.
-There is a non-zero chance for this to occur during normal operation.
-Re-run the test to see if the error persists."""
-
-        if st[0] > t_start+4.0*1.0/rate*1000.0:
-            raise StatisticalError(err)
-
-        err = """
-Number of spikes should be within 3 standard deviations of mean.
-There is a non-zero chance for this to occur during normal operation.
-Re-run the test to see if the error persists."""
-
-        # time interval is one second
-
-        if len(st) > rate+3.0*numpy.sqrt(rate) or len(st) < rate-3.0*numpy.sqrt(rate):
-            raise StatisticalError(err)
-
-
-    def testStatsInhPoisson(self):
-
-        # this is a statistical test with non-zero chance of failure
-
-        stg = stgen.StGen()
-
-        rate = 100.0 #Hz
-        t_start = 500.0
-        t_stop = 1500.0 # milliseconds
-
-        st = stg.inh_poisson_generator(numpy.array([rate]),numpy.array([t_start]),t_stop,array=True)
-