Repository: StingraySoftware/stingray
Pull request #141 (2016-08-16 13:13:06)

Another go at Uncertainties
:warning:   **NEEDS REVIEW** :warning: 

This PR substitutes #24 and already implements the changes proposed in #140

This implements uncertainties to Ligtcurves, Powserpectrum, AveragedPowerspectrum, Crosspectrum, and AveragedCrosspectrum (For the Crosspectra, I'm not sure if the current implementation is valid)
- In `Lightcurve`:  
  If no errors are given it assumes a poisson error given by: 
  `error = 1 + sqrt(N + 0.75)` (the [Gehrels1986](http://articles.adsabs.harvard.edu/cgi-bin/nph-iarticle_query?1986ApJ...303..336G&amp;data_type=PDF_HIGH&amp;whole_paper=YES&amp;type=PRINTER&amp;filetype=.pdf) approximation, [link](http://stsdas.stsci.edu/cgi-bin/gethelp.cgi?explain_errors) ) , as suggested by @nithinsingh61
- In `Powerspectrum`:  
  The error is given by: 
  `power_err = power / sqrt(M)`  
  Where `M` is the number of frequencies averaged (either by rebining or by averaging segmented ps)

<!-- Reviewable:start -->

---

This change is [<img src="https://reviewable.io/review_button.svg" height="34" align="absmiddle" alt="Reviewable"/>](https://reviewable.io/reviews/stingraysoftware/stingray/141)

<!-- Reviewable:end -->



Patch:
diff --git a/stingray/crossspectrum.py b/stingray/crossspectrum.py
index b5701a76..04e3c238 100644
--- a/stingray/crossspectrum.py
+++ b/stingray/crossspectrum.py
@@ -1,4 +1,4 @@
-from __future__ import division
+from __future__ import division, absolute_import, print_function
 
 import numpy as np
 import scipy
@@ -6,14 +6,18 @@ import scipy.stats
 import scipy.fftpack
 import scipy.optimize
 
-import stingray.lightcurve as lightcurve
-import stingray.utils as utils
+from stingray.lightcurve import Lightcurve
+from stingray.utils import rebin_data, simon
 from stingray.exceptions import StingrayError
 from stingray.gti import cross_two_gtis, bin_intervals_from_gtis, check_gtis
 
 __all__ = ["Crossspectrum", "AveragedCrossspectrum", "coherence"]
 
 
+def _root_squared_mean(array):
+    return np.sqrt(np.sum(array**2)) / len(array)
+
+
 def coherence(lc1, lc2):
     """
     Estimate coherence function of two light curves.
@@ -32,10 +36,10 @@ def coherence(lc1, lc2):
         Coherence function
     """
 
-    if not isinstance(lc1, lightcurve.Lightcurve):
+    if not isinstance(lc1, Lightcurve):
         raise TypeError("lc1 must be a lightcurve.Lightcurve object")
 
-    if not isinstance(lc2, lightcurve.Lightcurve):
+    if not isinstance(lc2, Lightcurve):
         raise TypeError("lc2 must be a lightcurve.Lightcurve object")
 
     cs = Crossspectrum(lc1, lc2, norm='none')
@@ -78,6 +82,13 @@ class Crossspectrum(object):
         power: numpy.ndarray
             The array of cross spectra (complex numbers)
 
+        power_err: numpy.ndarray
+            The uncertainties of `power`.
+            An approximation for each bin given by "power_err= power/Sqrt(m)".
+            Where `m` is the number of power averaged in each bin (by frequency
+            binning, or averaging more than one spectra). Note that for a single
+            realization (m=1) the error is equal to the power.
+
         df: float
             The frequency resolution
 
@@ -112,6 +123,7 @@ class Crossspectrum(object):
             else:
                 self.freq = None
                 self.power = None
+                self.power_err = None
                 self.df = None
                 self.nphots1 = None
                 self.nphots2 = None
@@ -126,10 +138,10 @@ class Crossspectrum(object):
     def _make_crossspectrum(self, lc1, lc2):
 
         # make sure the inputs work!
-        if not isinstance(lc1, lightcurve.Lightcurve):
+        if not isinstance(lc1, Lightcurve):
             raise TypeError("lc1 must be a lightcurve.Lightcurve object")
 
-        if not isinstance(lc2, lightcurve.Lightcurve):
+        if not isinstance(lc2, Lightcurve):
             raise TypeError("lc2 must be a lightcurve.Lightcurve object")
 
         # Then check that GTIs make sense
@@ -179,6 +191,30 @@ class Crossspectrum(object):
         # with the imaginary part still intact.
         self.power = self._normalize_crossspectrum(self.unnorm_power, lc1.tseg)
 
+        if lc1.err_dist.lower() != lc2.err_dist.lower():
+            simon("Your lightcurves have different statistics."
+                  "The errors in the Crossspectrum will be incorrect.")
+        elif lc1.err_dist.lower() != "poisson":
+            simon("Looks like your lightcurve statistic is not poisson."
+                  "The errors in the Powerspectrum will be incorrect.")
+
+        if self.__class__.__name__ in ['Powerspectrum',
+                                       'AveragedPowerspectrum']:
+            self.power_err = self.power / np.sqrt(self.m)
+        elif self.__class__.__name__ in ['Crossspectrum',
+                                         'AveragedCrossspectrum']:
+            # This is clearly a wild approximation.
+            simon("Errorbars on cross spectra are not thoroughly tested. "
+                  "Please report any inconsistencies.")
+            unnorm_power_err = np.sqrt(2) / np.sqrt(self.m) # Leahy-like
+            unnorm_power_err /= (2 / np.sqrt(self.nphots1 * self.nphots2))
+            unnorm_power_err += np.zeros_like(self.power)
+
+            self.power_err = \
+                self._normalize_crossspectrum(unnorm_power_err, lc1.tseg)
+        else:
+            self.power_err = np.zeros(len(self.power))
+
     def _fourier_cross(self, lc1, lc2):
         """
         Fourier transform the two light curves, then compute the cross spectrum.
@@ -225,9 +261,9 @@ class Crossspectrum(object):
         """
 
         # rebin cross spectrum to new resolution
-        binfreq, bincs, step_size = utils.rebin_data(self.freq,
-                                                     self.power, df,
-                                                     method=method)
+        binfreq, bincs, binerr, step_size = \
+            rebin_data(self.freq, self.power, df, self.power_err,
+                       method=method)
 
         # make an empty cross spectrum object
         # note: syntax deliberate to work with subclass Powerspectrum
@@ -240,6 +276,7 @@ class Crossspectrum(object):
         bin_cs.n = self.n
         bin_cs.norm = self.norm
         bin_cs.nphots1 = self.nphots1
+        bin_cs.power_err = binerr
 
         try:
             bin_cs.nphots2 = self.nphots2
@@ -327,6 +364,9 @@ class Crossspectrum(object):
         binpower: numpy.ndarray
             the binned powers
 
+        binpower_err: numpy.ndarray
+            the uncertainties in binpower
+
         nsamples: numpy.ndarray
             the samples of the original periodogram included in each
             frequency bin
@@ -349,6 +389,10 @@ class Crossspectrum(object):
             self.freq.astype(np.double), self.power.astype(np.double),
             statistic="mean", bins=binfreq)
 
+        binpower_err, bin_edges, binno = scipy.stats.binned_statistic(
+            self.freq.astype(np.double), self.power_err.astype(np.double),
+            statistic=_root_squared_mean, bins=binfreq)
+
         # compute the number of powers in each frequency bin
         nsamples = np.array([len(binno[np.where(binno == i)[0]])
                              for i in range(np.max(binno))])
@@ -360,7 +404,7 @@ class Crossspectrum(object):
         # last right bin edge
         binfreq = binfreq[:-1] + df/2
 
-        return binfreq, binpower, nsamples
+        return binfreq, binpower, binpower_err, nsamples
 
     def coherence(self):
         """
@@ -383,7 +427,7 @@ class Crossspectrum(object):
         ps1 = Crossspectrum(self.lc1, self.lc1)
         ps2 = Crossspectrum(self.lc2, self.lc2)
 
-        return self.unnorm_power/(ps1.unnorm_power * ps2.unnorm_power)
+        return self.unnorm_power.real/(ps1.unnorm_power.real * ps2.unnorm_power.real)
 
     def _phase_lag(self):
         """Return the fourier phase lag of the cross spectrum."""
@@ -449,6 +493,13 @@ class AveragedCrossspectrum(Crossspectrum):
         power: numpy.ndarray
             The array of cross spectra
 
+        power_err: numpy.ndarray
+            The uncertainties of `power`.
+            An approximation for each bin given by "power_err= power/Sqrt(m)".
+            Where `m` is the number of power averaged in each bin (by frequency
+            binning, or averaging powerspectrum). Note that for a single
+            realization (m=1) the error is equal to the power.
+
         df: float
             The frequency resolution
 
@@ -485,8 +536,8 @@ class AveragedCrossspectrum(Crossspectrum):
     def _make_segment_spectrum(self, lc1, lc2, segment_size):
 
         # TODO: need to update this for making cross spectra.
-        assert isinstance(lc1, lightcurve.Lightcurve)
-        assert isinstance(lc2, lightcurve.Lightcurve)
+        assert isinstance(lc1, Lightcurve)
+        assert isinstance(lc2, Lightcurve)
 
         if lc1.dt != lc2.dt:
             raise ValueError("Light curves do not have same time binning dt.")
@@ -509,10 +560,14 @@ class AveragedCrossspectrum(Crossspectrum):
         for start_ind, end_ind in zip(start_inds, end_inds):
             time_1 = lc1.time[start_ind:end_ind]
             counts_1 = lc1.counts[start_ind:end_ind]
+            counts_1_err = lc1.counts_err[start_ind:end_ind]
             time_2 = lc2.time[start_ind:end_ind]
             counts_2 = lc2.counts[start_ind:end_ind]
-            lc1_seg = lightcurve.Lightcurve(time_1, counts_1)
-            lc2_seg = lightcurve.Lightcurve(time_2, counts_2)
+            counts_2_err = lc2.counts_err[start_ind:end_ind]
+            lc1_seg = Lightcurve(time_1, counts_1, err=counts_1_err,
+                                 err_dist=lc1.err_dist)
+            lc2_seg = Lightcurve(time_2, counts_2, err=counts_2_err,
+                                 err_dist=lc2.err_dist)
             cs_seg = Crossspectrum(lc1_seg, lc2_seg, norm=self.norm)
             cs_all.append(cs_seg)
             nphots1_all.append(np.sum(lc1_seg.counts))
@@ -523,8 +578,8 @@ class AveragedCrossspectrum(Crossspectrum):
     def _make_crossspectrum(self, lc1, lc2):
 
         # chop light curves into segments
-        if isinstance(lc1, lightcurve.Lightcurve) and \
-                isinstance(lc2, lightcurve.Lightcurve):
+        if isinstance(lc1, Lightcurve) and \
+                isinstance(lc2, Lightcurve):
 
             if self.type == "crossspectrum":
                 self.cs_all, nphots1_all, nphots2_all = \
@@ -568,14 +623,18 @@ class AveragedCrossspectrum(Crossspectrum):
         nphots1 = np.mean(nphots1_all)
 
         power_avg = np.zeros_like(self.cs_all[0].power)
+        power_err_avg = np.zeros_like(self.cs_all[0].power_err)
         for cs in self.cs_all:
             power_avg += cs.power
+            power_err_avg += (cs.power_err)**2
 
         power_avg /= np.float(m)
+        power_err_avg = np.sqrt(power_err_avg) / m
 
         self.freq = self.cs_all[0].freq
         self.power = power_avg
         self.m = m
+        self.power_err = power_err_avg
         self.df = self.cs_all[0].df
         self.n = self.cs_all[0].n
         self.nphots1 = nphots1
@@ -607,9 +666,9 @@ class AveragedCrossspectrum(Crossspectrum):
 
         """
         if self.m < 50:
-            utils.simon("Number of segments used in averaging is "
-                        "significantly low. The result might not follow the "
-                        "expected statistical distributions.")
+            simon("Number of segments used in averaging is "
+                  "significantly low. The result might not follow the "
+                  "expected statistical distributions.")
 
         # Calculate average coherence
         unnorm_power_avg = np.zeros_like(self.cs_all[0].unnorm_power)
@@ -617,7 +676,7 @@ class AveragedCrossspectrum(Crossspectrum):
             unnorm_power_avg += cs.unnorm_power
 
         unnorm_power_avg /= self.m
-        num = np.abs(unnorm_power_avg)**2
+        num = np.absolute(unnorm_power_avg)**2
 
         # this computes the averaged power spectrum, but using the
         # cross spectrum code to avoid circular imports
@@ -626,14 +685,14 @@ class AveragedCrossspectrum(Crossspectrum):
         aps2 = AveragedCrossspectrum(self.lc2, self.lc2,
                                      segment_size=self.segment_size)
 
-        unnorm_powers_avg_1 = np.zeros_like(aps1.cs_all[0].unnorm_power)
+        unnorm_powers_avg_1 = np.zeros_like(aps1.cs_all[0].unnorm_power.real)
         for ps in aps1.cs_all:
-            unnorm_powers_avg_1 += ps.unnorm_power
+            unnorm_powers_avg_1 += ps.unnorm_power.real
         unnorm_powers_avg_1 /= aps1.m
 
-        unnorm_powers_avg_2 = np.zeros_like(aps2.cs_all[0].unnorm_power)
+        unnorm_powers_avg_2 = np.zeros_like(aps2.cs_all[0].unnorm_power.real)
         for ps in aps2.cs_all:
-            unnorm_powers_avg_2 += ps.unnorm_power
+            unnorm_powers_avg_2 += ps.unnorm_power.real
         unnorm_powers_avg_2 /= aps2.m
 
         coh = num / (unnorm_powers_avg_1 * unnorm_powers_avg_2)
@@ -642,3 +701,14 @@ class AveragedCrossspectrum(Crossspectrum):
         uncertainty = (2**0.5 * coh * (1 - coh)) / (np.abs(coh) * self.m**0.5)
 
         return (coh, uncertainty)
+
+    def time_lag(self):
+        """Calculate time lag and uncertainty.
+        
+        Formula from Bendat & Piersol 1986
+        """
+        lag = super(AveragedCrossspectrum, self).time_lag()
+        coh, uncert = self.coherence()
+        dum = (1. - coh) / (2. * coh)
+        lag_err = np.sqrt(dum / self.m) / (2 * np.pi * self.freq)
+        return (lag, lag_err)
diff --git a/stingray/lightcurve.py b/stingray/lightcurve.py
index f92d05fa..e4a4f3f2 100644
--- a/stingray/lightcurve.py
+++ b/stingray/lightcurve.py
@@ -12,12 +12,16 @@ from stingray.exceptions import StingrayError
 from stingray.utils import simon, assign_value_if_none
 from stingray.gti import cross_two_gtis, join_gtis, gti_border_bins
 from stingray.gti import check_gtis
+from astropy.stats import poisson_conf_interval
 
 __all__ = ["Lightcurve"]
 
+valid_statistics = ["poisson", "gauss", None]
+
 
 class Lightcurve(object):
-    def __init__(self, time, counts, input_counts=True, gti=None):
+    def __init__(self, time, counts, err=None, input_counts=True,
+                 gti=None, err_dist='poisson'):
         """
         Make a light curve object from an array of time stamps and an
         array of counts.
@@ -29,8 +33,18 @@ class Lightcurve(object):
 
         counts: iterable, optional, default None
             A list or array of the counts in each bin corresponding to the
-            bins defined in `time` (note: **not** the count rate, i.e.
-            counts/second, but the counts/bin).
+            bins defined in `time` (note: use `input_counts=False` to
+            input the count range, i.e. counts/second, otherwise use
+            counts/bin).
+
+        err: iterable, optional, default None:
+            A list or array of the uncertainties in each bin corresponding to
+            the bins defined in `time` (note: use `input_counts=False` to
+            input the count rage, i.e. counts/second, otherwise use
+            counts/bin). If None, we assume the data is poisson distributed
+            and calculate the error from the average of the lower and upper 
+            1-sigma confidence intervals for the Poissonian distribution with 
+            mean equal to `counts`.
 
         input_counts: bool, optional, default True
             If True, the code assumes that the input data in 'counts'
@@ -43,6 +57,12 @@ class Lightcurve(object):
             They will be used by other methods to have an indication of the
             "safe" time intervals to use during analysis.
 
+        err_dist: str, optional, default=None
+            Statistic of the Lightcurve, it is used to calculate the
+            uncertainties and other statistical values apropriately.
+            Default makes no assumptions and keep errors equal to zero.
+
+
         Attributes
         ----------
         time: numpy.ndarray
@@ -57,12 +77,18 @@ class Lightcurve(object):
         counts: numpy.ndarray
             The counts per bin corresponding to the bins in `time`.
 
+        counts_err: numpy.ndarray
+            The uncertainties corresponding to `counts`
+
         countrate: numpy.ndarray
             The counts per second in each of the bins defined in `time`.
 
+        countrate_err: numpy.ndarray
+            The uncertainties corresponding to `countrate`
+
         meanrate: float
             The mean count rate of the light curve.
-            
+
         meancounts: float
             The mean counts of the light curve.
 
@@ -83,7 +109,13 @@ class Lightcurve(object):
             Good Time Intervals. They indicate the "safe" time intervals
             to be used during the analysis of the light curve.
 
+        err_dist: string
+            Statistic of the Lightcurve, it is used to calculate the
+            uncertainties and other statistical values apropriately.
+            It propagates to Spectrum classes.
+
         """
+
         if not np.all(np.isfinite(time)):
             raise ValueError("There are inf or NaN values in "
                              "your time array!")
@@ -101,18 +133,46 @@ class Lightcurve(object):
             raise StingrayError("A single or no data points can not create "
                                 "a lightcurve!")
 
+        if err is not None:
+            if not np.all(np.isfinite(err)):
+                raise ValueError("There are inf or NaN values in "
+                                 "your err array")
+        else:
+            if err_dist.lower() not in valid_statistics:
+                # err_dist set can be increased with other statistics
+                raise StingrayError("Statistic not recognized."
+                                    "Please select one of these: ",
+                                    "{}".format(valid_statistics))
+            if err_dist.lower() == 'poisson':
+                err_low, err_high = poisson_conf_interval(np.asarray(counts),
+                    interval='frequentist-confidence', sigma=1)
+                # calculate approximately symetric uncertainties
+                err = (np.absolute(err_low)+np.absolute(err_high))/2.0
+                # other estimators can be implemented for other statistics
+            else:
+                simon("Stingray only uses poisson err_dist at the moment, "
+                      "We are setting your errors to zero. "
+                      "Sorry for the inconvenience.")
+                err = np.zeros_like(counts)
+
         self.time = np.asarray(time)
         self.dt = time[1] - time[0]
 
         self.bin_lo = self.time - 0.5 * self.dt
         self.bin_hi = self.time + 0.5 * self.dt
 
+        self.err_dist = err_dist
+
         if input_counts:
             self.counts = np.asarray(counts)
             self.countrate = self.counts / self.dt
+            self.counts_err = np.asarray(err)
+            self.countrate_err = np.asarray(err) / self.dt
         else:
             self.countrate = np.asarray(counts)
             self.counts = self.countrate * self.dt
+            self.counts_err = np.asarray(err) * self.dt
+            self.countrate_err = np.asarray(err)
 
         self.meanrate = np.mean(self.countrate)
         self.meancounts = np.mean(self.counts)
@@ -147,6 +207,7 @@ class Lightcurve(object):
                             gti=self.gti + time_shift)
         new_lc.countrate = self.countrate
         new_lc.counts = self.counts
+        new_lc.counts_err = self.counts_err
         new_lc.meanrate = np.mean(new_lc.countrate)
         new_lc.meancounts = np.mean(new_lc.counts)
         return new_lc
@@ -183,9 +244,24 @@ class Lightcurve(object):
                              "of same dimension and equal.")
 
         new_counts = np.add(self.counts, other.counts)
+
+        if self.err_dist.lower() != other.err_dist.lower():
+            simon("Lightcurves have different statistics!"
+                  "We are setting the errors to zero to avoid complications.")
+            new_counts_err = np.zeros_like(new_counts)
+        elif self.err_dist.lower() in valid_statistics:
+                new_counts_err = np.sqrt(np.add(self.counts_err**2,
+                                                other.counts_err**2))
+            # More conditions can be implemented for other statistics
+        else:
+            raise StingrayError("Statistics not recognized."
+                                " Please use one of these: "
+                                "{}".format(valid_statistics))
+
         common_gti = cross_two_gtis(self.gti, other.gti)
 
-        lc_new = Lightcurve(self.time, new_counts, gti=common_gti)
+        lc_new = Lightcurve(self.time, new_counts,
+                            err=new_counts_err, gti=common_gti)
 
         return lc_new
 
@@ -222,9 +298,24 @@ class Lightcurve(object):
                              "of same dimension and equal.")
 
         new_counts = np.subtract(self.counts, other.counts)
+
+        if self.err_dist.lower() != other.err_dist.lower():
+            simon("Lightcurves have different statistics!"
+                  "We are setting the errors to zero to avoid complications.")
+            new_counts_err = np.zeros_like(new_counts)
+        elif self.err_dist.lower() in valid_statistics:
+            new_counts_err = np.sqrt(np.add(self.counts_err**2,
+                                            other.counts_err**2))
+            # More conditions can be implemented for other statistics
+        else:
+            raise StingrayError("Statistics not recognized."
+                                " Please use one of these: "
+                                "{}".format(valid_statistics))
+
         common_gti = cross_two_gtis(self.gti, other.gti)
 
-        lc_new = Lightcurve(self.time, new_counts, gti=common_gti)
+        lc_new = Lightcurve(self.time, new_counts,
+                            err=new_counts_err, gti=common_gti)
 
         return lc_new
 
@@ -246,7 +337,8 @@ class Lightcurve(object):
         >>> lc_new.counts
         array([100, 100, 100])
         """
-        lc_new = Lightcurve(self.time, -1*self.counts, gti=self.gti)
+        lc_new = Lightcurve(self.time, -1*self.counts,
+                            err=self.counts_err, gti=self.gti)
 
         return lc_new
 
@@ -348,8 +440,8 @@ class Lightcurve(object):
             tstart = toa[0]
 
         # compute the number of bins in the light curve
-        # for cases where tseg/dt is not integer. 
-        # TODO: check that this is always consistent and that we 
+        # for cases where tseg/dt is not integer.
+        # TODO: check that this is always consistent and that we
         # are not throwing away good events.
 
         if tseg is None:
@@ -401,11 +493,11 @@ class Lightcurve(object):
             raise ValueError("New time resolution must be larger than "
                              "old time resolution!")
 
-        bin_time, bin_counts, _ = utils.rebin_data(self.time,
-                                                   self.counts,
-                                                   dt_new, method)
+        bin_time, bin_counts, bin_err, _ = \
+            utils.rebin_data(self.time, self.counts, dt_new,
+                             yerr=self.counts_err, method=method)
 
-        lc_new = Lightcurve(bin_time, bin_counts)
+        lc_new = Lightcurve(bin_time, bin_counts, err=bin_err)
         return lc_new
 
     def join(self, other):
@@ -461,33 +553,53 @@ class Lightcurve(object):
                         "lc2`.")
 
         new_counts = []
-
+        new_counts_err = []
         # For every time stamp, get the individual time counts and add them.
         for time in new_time:
             try:
                 count1 = self.counts[np.where(self.time == time)[0][0]]
+                count1_err = self.counts_err[np.where(self.time == time)[0][0]]
             except IndexError:
                 count1 = None
+                count1_err = None
 
             try:
                 count2 = other.counts[np.where(other.time == time)[0][0]]
+                count2_err = other.counts_err[np.where(other.time == time)[0][0]]
             except IndexError:
                 count2 = None
+                count2_err = None
 
             if count1 is not None:
                 if count2 is not None:
                     # Average the overlapping counts
                     new_counts.append((count1 + count2) / 2)
+
+                    if self.err_dist.lower() != other.err_dist.lower():
+                        simon("Lightcurves have different statistics!"
+                              "We are setting the errors to zero.")
+                        new_counts_err = np.zeros_like(new_counts)
+                    elif self.err_dist.lower() in valid_statistics:
+                        new_counts_err.append(np.sqrt(((count1_err**2) +
+                                                      (count2_err**2)) / 2))
+                    # More conditions can be implemented for other statistics
+                    else:
+                        raise StingrayError("Statistics not recognized."
+                                            " Please use one of these: "
+                                            "{}".format(valid_statistics))
                 else:
                     new_counts.append(count1)
+                    new_counts_err.append(count1_err)
             else:
                 new_counts.append(count2)
+                new_counts_err.append(count2_err)
 
         new_counts = np.asarray(new_counts)
+        new_counts_err = np.asarray(new_counts_err)
 
         gti = join_gtis(self.gti, other.gti)
 
-        lc_new = Lightcurve(new_time, new_counts, gti=gti)
+        lc_new = Lightcurve(new_time, new_counts, err=new_counts_err, gti=gti)
 
         return lc_new
 
@@ -554,12 +666,13 @@ class Lightcurve(object):
         """Private method for truncation using index values."""
         time_new = self.time[start:stop]
         counts_new = self.counts[start:stop]
+        counts_err_new = self.counts_err[start:stop]
         gti = \
             cross_two_gtis(self.gti,
                            np.asarray([[self.time[start] - 0.5 * self.dt,
                                         time_new[-1] + 0.5 * self.dt]]))
 
-        return Lightcurve(time_new, counts_new, gti=gti)
+        return Lightcurve(time_new, counts_new, err=counts_err_new, gti=gti)
 
     def _truncate_by_time(self, start, stop):
         """Private method for truncation using time values."""
@@ -608,26 +721,34 @@ class Lightcurve(object):
         """
         new_counts = sorted(self.counts, reverse=reverse)
         new_time = []
+        new_counts_err = []
         for count in np.unique(new_counts):
             for index in np.where(self.counts == count)[0]:
                 new_time.append(self.time[index])
+                new_counts_err.append(self.counts_err[index])
 
         if reverse:
             new_time.reverse()
+            new_counts_err.reverse()
 
         self.time = np.asarray(new_time)
         self.counts = np.asarray(new_counts)
+        self.counts_err = np.asarray(new_counts_err)
 
-    def plot(self, labels=None, axis=None, title=None, marker='-', save=False,
-             show=False, filename=None):
+    def plot(self, witherrors=False, labels=None, axis=None, title=None,
+             marker='-', save=False, filename=None):
         """
         Plot the Lightcurve using Matplotlib.
 
         Plot the Lightcurve object on a graph ``self.time`` on x-axis and
-        ``self.counts`` on y-axis.
+        ``self.counts`` on y-axis with ``self.counts_err`` optionaly
+        as error bars.
 
         Parameters
         ----------
+        witherrors: boolean, default False
+            Whether to plot the Lightcurve with errorbars or not
+
         labels : iterable, default None
             A list of tuple with xlabel and ylabel as strings.
 
@@ -656,7 +777,11 @@ class Lightcurve(object):
             raise ImportError("Matplotlib required for plot()")
 
         fig = plt.figure()
-        fig = plt.plot(self.time, self.counts, marker)
+        if witherrors:
+            fig = plt.errorbar(self.time, self.counts, yerr=self.counts_err,
+                               fmt=marker)
+        else:
+            fig = plt.plot(self.time, self.counts, marker)
 
         if labels is not None:
             try:
@@ -749,7 +874,8 @@ class Lightcurve(object):
             start = start_bins[i]
             stop = stop_bins[i]
             # Note: GTIs are consistent with default in this case!
-            new_lc = Lightcurve(self.time[start:stop], self.counts[start:stop])
+            new_lc = Lightcurve(self.time[start:stop], self.counts[start:stop],
+                                err=self.counts_err[start:stop])
             list_of_lcs.append(new_lc)
 
         return list_of_lcs
diff --git a/stingray/powerspectrum.py b/stingray/powerspectrum.py
index 9d6b9418..bd8b7a61 100755
--- a/stingray/powerspectrum.py
+++ b/stingray/powerspectrum.py
@@ -154,6 +154,13 @@ class Powerspectrum(Crossspectrum):
             The array of normalized squared absolute values of Fourier
             amplitudes
 
+        power_err: numpy.ndarray
+            The uncertainties of `power`.
+            An approximation for each bin given by "power_err= power/Sqrt(m)".
+            Where `m` is the number of power averaged in each bin (by frequency
+            binning, or averaging powerspectrum). Note that for a single
+            realization (m=1) the error is equal to the power.
+
         df: float
             The frequency resolution
 
@@ -167,6 +174,7 @@ class Powerspectrum(Crossspectrum):
             The total number of photons in the light curve
 
         """
+
         Crossspectrum.__init__(self, lc1=lc, lc2=lc, norm=norm, gti=gti)
         self.nphots = self.nphots1
 
@@ -350,6 +358,13 @@ class AveragedPowerspectrum(AveragedCrossspectrum, Powerspectrum):
             The array of normalized squared absolute values of Fourier
             amplitudes
 
+        power_err: numpy.ndarray
+            The uncertainties of `power`.
+            An approximation for each bin given by "power_err= power/Sqrt(m)".
+            Where `m` is the number of power averaged in each bin (by frequency
+            binning, or averaging powerspectrum). Note that for a single
+            realization (m=1) the error is equal to the power.
+
         df: float
             The frequency resolution
 
@@ -394,7 +409,9 @@ class AveragedPowerspectrum(AveragedCrossspectrum, Powerspectrum):
         for start_ind, end_ind in zip(start_inds, end_inds):
             time = lc.time[start_ind:end_ind]
             counts = lc.counts[start_ind:end_ind]
-            lc_seg = lightcurve.Lightcurve(time, counts)
+            counts_err = lc.counts_err[start_ind: end_ind]
+            lc_seg = lightcurve.Lightcurve(time, counts, err=counts_err,
+                                           err_dist=lc.err_dist.lower())
             power_seg = Powerspectrum(lc_seg, norm=self.norm)
             power_all.append(power_seg)
             nphots_all.append(np.sum(lc_seg.counts))
diff --git a/stingray/simulator/simulator.py b/stingray/simulator/simulator.py
index fbe3de82..bbf6ae43 100644
--- a/stingray/simulator/simulator.py
+++ b/stingray/simulator/simulator.py
@@ -9,7 +9,8 @@ import stingray.simulator.models as models
 
 class Simulator(object):
 
-    def __init__(self, dt=1, N=1024, mean=0, rms=1, red_noise=1, random_state=None):
+    def __init__(self, dt=1, N=1024, mean=0, rms=1, red_noise=1,
+                 random_state=None):
         """
         Methods to simulate and visualize light curves.
 
@@ -356,7 +357,8 @@ class Simulator(object):
 
         # Obtain time series
         long_lc = self._find_inverse(real, imaginary)
-        lc = Lightcurve(self.time, self._extract_and_scale(long_lc))
+        lc = Lightcurve(self.time, self._extract_and_scale(long_lc),
+                        err_dist='gauss')
 
         return lc
 
@@ -384,7 +386,8 @@ class Simulator(object):
         a2 = self.random_state.normal(size=len(s))
 
         lc = self._find_inverse(a1*s, a2*s)
-        lc = Lightcurve(self.time, self._extract_and_scale(lc))
+        lc = Lightcurve(self.time, self._extract_and_scale(lc),
+                        err_dist='gauss')
 
         return lc
 
@@ -424,7 +427,8 @@ class Simulator(object):
         # Inverse Fourier transform
         long_lc = np.fft.irfft(arg)
 
-        lc = Lightcurve(self.time, self._extract_and_scale(long_lc))
+        lc = Lightcurve(self.time, self._extract_and_scale(long_lc),
+                        err_dist='gauss')
         return lc
 
 
@@ -471,7 +475,8 @@ class Simulator(object):
             # Inverse Fourier transform
             long_lc = np.fft.irfft(arg)
 
-            lc = Lightcurve(self.time, self._extract_and_scale(long_lc))
+            lc = Lightcurve(self.time, self._extract_and_scale(long_lc),
+                            err_dist='gauss')
             return lc
         else:
             raise ValueError('Model is not defined!')
@@ -511,7 +516,7 @@ class Simulator(object):
             lc = lc[(len(h) - 1):-(len(h) - 1)]
 
         time = self.dt * np.arange(len(lc))
-        return Lightcurve(time, lc)
+        return Lightcurve(time, lc, err_dist='gauss')
 
     def _find_inverse(self, real, imaginary):
         """
diff --git a/stingray/utils.py b/stingray/utils.py
index b156c57e..66ad7a09 100644
--- a/stingray/utils.py
+++ b/stingray/utils.py
@@ -35,7 +35,7 @@ def simon(message, **kwargs):
     warnings.warn("SIMON says: {0}".format(message), **kwargs)
 
 
-def rebin_data(x, y, dx_new, method='sum'):
+def rebin_data(x, y, dx_new, yerr=None, method='sum'):
 
     """Rebin some data to an arbitrary new data resolution. Either sum
     the data points in the new bins or average them.
@@ -48,6 +48,9 @@ def rebin_data(x, y, dx_new, method='sum'):
     y: iterable
         The independent variable to be binned
 
+    yerr: iterable, optional
+        The uncertainties of y, to be propagated during binning.
+
     dx_new: float
         The new resolution of the dependent variable x
 
@@ -63,9 +66,16 @@ def rebin_data(x, y, dx_new, method='sum'):
 
     ybin: numpy.ndarray
         The binned quantity y
+
+    ybin_err: numpy.ndarray
+        The uncertainties of the binned values of y.
+
+    step_size: float
+        The size of the binning step
     """
 
     y = np.asarray(y)
+    yerr = np.asarray(assign_value_if_none(yerr, np.zeros_like(y)))
 
     dx_old = x[1] - x[0]
 
@@ -76,30 +86,39 @@ def rebin_data(x, y, dx_new, method='sum'):
     step_size = dx_new / dx_old
 
     output = []
+    outputerr = []
     for i in np.arange(0, y.shape[0], step_size):
         total = 0
+        totalerr = 0
 
         int_i = int(i)
         prev_frac = int_i + 1 - i
         prev_bin = int_i
         total += prev_frac * y[prev_bin]
+        totalerr += prev_frac * (yerr[prev_bin]**2)
 
         if i + step_size < len(x):
             # Fractional part of next bin:
             next_frac = i + step_size - int(i + step_size)
             next_bin = int(i + step_size)
             total += next_frac * y[next_bin]
+            totalerr += next_frac * (yerr[next_bin]**2)
 
         total += sum(y[int(i+1):int(i+step_size)])
+        totalerr += sum(yerr[int(i+1):int(step_size)]**2)
         output.append(total)
+        outputerr.append(np.sqrt(totalerr))
 
     output = np.asarray(output)
+    outputerr = np.asarray(outputerr)
 
     if method in ['mean', 'avg', 'average', 'arithmetic mean']:
         ybin = output / np.float(step_size)
+        ybinerr = outputerr / np.sqrt(np.float(step_size))
 
     elif method == "sum":
         ybin = output
+        ybinerr = outputerr
 
     else:
         raise ValueError("Method for summing or averaging not recognized. "
@@ -109,11 +128,12 @@ def rebin_data(x, y, dx_new, method='sum'):
 
     if (tseg / dx_new % 1) > 0:
         ybin = ybin[:-1]
+        ybinerr = ybinerr[:-1]
 
     new_x0 = (x[0] - (0.5*dx_old)) + (0.5*dx_new)
     xbin = np.arange(ybin.shape[0]) * dx_new + new_x0
 
-    return xbin, ybin, step_size
+    return xbin, ybin, ybinerr, step_size
 
 
 def assign_value_if_none(value, default):


Test patch:
diff --git a/stingray/tests/test_crossspectrum.py b/stingray/tests/test_crossspectrum.py
index 9872b14a..2edb084d 100644
--- a/stingray/tests/test_crossspectrum.py
+++ b/stingray/tests/test_crossspectrum.py
@@ -85,6 +85,7 @@ class TestCrossspectrum(object):
         assert cs.nphots2 is None
         assert cs.m == 1
         assert cs.n is None
+        assert cs.power_err is None
 
     def test_init_with_one_lc_none(self):
         with pytest.raises(TypeError):
@@ -136,6 +137,10 @@ class TestCrossspectrum(object):
         new_cs = self.cs.rebin(df=1.5)
         assert new_cs.df == 1.5
 
+    def test_rebin_log(self):
+        # For now, just verify that it doesn't crash
+        _ = self.cs.rebin_log(f=0.01)
+
     def test_norm_leahy(self):
         cs = Crossspectrum(self.lc1, self.lc2, norm='leahy')
         assert len(cs.power) == 4999
@@ -165,6 +170,9 @@ class TestCrossspectrum(object):
         assert max(time_lag) <= np.pi
         assert min(time_lag) >= -np.pi
 
+    def test_nonzero_err(self):
+        assert np.all(self.cs.power_err > 0)
+
     def test_timelag_error(self):
         class Child(Crossspectrum):
             def __init__(self):
@@ -298,7 +306,7 @@ class TestAveragedCrossspectrum(object):
             assert len(coh[0]) == 4999
             assert len(coh[1]) == 4999
 
-            assert len(w) == 1
+            assert len(w) == 3
             assert issubclass(w[-1].category, UserWarning)
 
     def test_failure_when_normalization_not_recognized(self):
@@ -306,3 +314,28 @@ class TestAveragedCrossspectrum(object):
             self.cs = AveragedCrossspectrum(self.lc1, self.lc2,
                                             segment_size=1,
                                             norm="wrong")
+
+    def test_timelag(self):
+        from ..simulator.simulator import Simulator
+        simulator = Simulator(0.1, 10000, rms=0.4, mean=200)
+        test_lc1 = simulator.simulate(2)
+        test_lc2 = Lightcurve(test_lc1.time,
+                              np.array(np.roll(test_lc1.counts, 2)),
+                              err_dist=test_lc1.err_dist)
+
+        cs = AveragedCrossspectrum(test_lc1, test_lc2, segment_size=10,
+                                   norm="none")
+
+        time_lag, time_lag_err = cs.time_lag()
+
+        assert np.all(np.abs(time_lag[:10] - 0.1) < 3 * time_lag_err[:10])
+
+    def test_errorbars(self):
+        time = np.arange(10000) * 0.1
+        test_lc1 = Lightcurve(time, np.random.poisson(200, 10000))
+        test_lc2 = Lightcurve(time, np.random.poisson(200, 10000))
+
+        cs = AveragedCrossspectrum(test_lc1, test_lc2, segment_size=10,
+                                   norm="leahy")
+
+        assert np.allclose(cs.power_err, np.sqrt(2/cs.m))
diff --git a/stingray/tests/test_lightcurve.py b/stingray/tests/test_lightcurve.py
index c3f888d7..77093807 100644
--- a/stingray/tests/test_lightcurve.py
+++ b/stingray/tests/test_lightcurve.py
@@ -22,8 +22,8 @@ class TestLightcurve(object):
 
     @classmethod
     def setup_class(cls):
-        cls.times = [1, 2, 3, 4]
-        cls.counts = [2, 2, 2, 2]
+        cls.times = np.array([1, 2, 3, 4])
+        cls.counts = np.array([2, 2, 2, 2])
         cls.dt = 1.0
         cls.gti = [[0.5, 4.5]]
 
@@ -45,9 +45,45 @@ class TestLightcurve(object):
                     "transforms. Please make the input time evenly sampled.")
 
         with warnings.catch_warnings(record=True) as w:
-            lc = Lightcurve(times, counts)
+            lc = Lightcurve(times, counts, err_dist="poisson")
+            assert str(w[0].message) == warn_str
+
+    def test_unrecognize_err_dist_warning(self):
+        """
+        Check if a non-poisson error_dist throws the correct warning.
+        """
+        times = [1, 2, 3, 4, 5]
+        counts = [2, 2, 2, 2, 2]
+        warn_str = ("SIMON says: Stingray only uses poisson err_dist at "
+                    "the moment, We are setting your errors to zero. Sorry for "
+                    "the inconvenience.")
+
+        with warnings.catch_warnings(record=True) as w:
+            lc = Lightcurve(times, counts, err_dist='gauss')
             assert str(w[0].message) == warn_str
 
+    def test_dummy_err_dist_fail(self):
+        """
+        Check if inputting an irregularly spaced time iterable throws out
+        a warning.
+        """
+        times = [1, 2, 3, 4, 5]
+        counts = [2, 2, 2, 2, 2]
+
+        with pytest.raises(StingrayError):
+            lc = Lightcurve(times, counts, err_dist='joke')
+
+    def test_invalid_data(self):
+        times = [1, 2, 3, 4, 5]
+        counts = [2, 2, np.nan, 2, 2]
+        counts_err = [1, 2, 3, np.nan, 2]
+
+        with pytest.raises(ValueError):
+            lc = Lightcurve(times, counts)
+
+        with pytest.raises(ValueError):
+            lc = Lightcurve(times, [2]*5, err=counts_err)
+
     def test_n(self):
         lc = Lightcurve(self.times, self.counts)
         assert lc.n == 4
@@ -187,6 +223,14 @@ class TestLightcurve(object):
 
             lc = lc1 + lc2
 
+    def test_add_with_different_err_dist(self):
+        lc1 = Lightcurve(self.times, self.counts)
+        lc2 = Lightcurve(self.times, self.counts, err=self.counts / 2,
+                         err_dist="gauss")
+        with warnings.catch_warnings(record=True) as w:
+            lc = lc1 + lc2
+            assert "ightcurves have different statistics" in str(w[0].message)
+
     def test_add_with_same_gtis(self):
         lc1 = Lightcurve(self.times, self.counts, gti=self.gti)
         lc2 = Lightcurve(self.times, self.counts, gti=self.gti)
@@ -230,6 +274,14 @@ class TestLightcurve(object):
 
             lc = lc1 - lc2
 
+    def test_sub_with_different_err_dist(self):
+        lc1 = Lightcurve(self.times, self.counts)
+        lc2 = Lightcurve(self.times, self.counts, err=self.counts / 2,
+                         err_dist="gauss")
+        with warnings.catch_warnings(record=True) as w:
+            lc = lc1 - lc2
+            assert "ightcurves have different statistics" in str(w[0].message)
+
     def test_subtraction(self):
         _counts = [3, 4, 5, 6]
 
@@ -523,4 +575,4 @@ class TestLightcurveRebin(object):
     def test_rebin_equal_numbers(self):
         dt_all = [2, 3, np.pi, 5]
         for dt in dt_all:
-            yield self.rebin_several, dt
+            self.rebin_several(dt)
diff --git a/stingray/tests/test_powerspectrum.py b/stingray/tests/test_powerspectrum.py
index acf0fd2a..fa8b6365 100644
--- a/stingray/tests/test_powerspectrum.py
+++ b/stingray/tests/test_powerspectrum.py
@@ -33,6 +33,7 @@ class TestPowerspectrum(object):
         assert ps.norm == "frac"
         assert ps.freq is None
         assert ps.power is None
+        assert ps.power_err is None
         assert ps.df is None
         assert ps.m == 1
         assert ps.n is None
@@ -41,6 +42,7 @@ class TestPowerspectrum(object):
         ps = Powerspectrum(lc=self.lc)
         assert ps.freq is not None
         assert ps.power is not None
+        assert ps.power_err is not None
         assert ps.df == 1.0 / self.lc.tseg
         assert ps.norm == "frac"
         assert ps.m == 1
@@ -51,6 +53,7 @@ class TestPowerspectrum(object):
         ps = Powerspectrum(lc=self.lc)
         assert isinstance(ps.freq, np.ndarray)
         assert isinstance(ps.power, np.ndarray)
+        assert isinstance(ps.power_err, np.ndarray)
 
     def test_init_with_lightcurve(self):
         assert Powerspectrum(self.lc)
@@ -174,6 +177,7 @@ class TestPowerspectrum(object):
 
         assert bin_ps.freq is not None
         assert bin_ps.power is not None
+        assert bin_ps.power is not None
         assert bin_ps.df == rebin_factor * 1.0 / self.lc.tseg
         assert bin_ps.norm.lower() == "leahy"
         assert bin_ps.m == 2
diff --git a/stingray/tests/test_utils.py b/stingray/tests/test_utils.py
index 38e918bb..25b5ed93 100644
--- a/stingray/tests/test_utils.py
+++ b/stingray/tests/test_utils.py
@@ -15,46 +15,53 @@ class TestRebinData(object):
         cls.counts = 2.0
         cls.x = np.arange(cls.dx/2, cls.dx/2+cls.n*cls.dx, cls.dx)
         cls.y = np.zeros_like(cls.x)+cls.counts
+        cls.yerr = np.sqrt(cls.y)
 
     def test_new_stepsize(self):
         dx_new = 2.0
-        xbin, ybin, step_size = utils.rebin_data(self.x, self.y, dx_new)
+        xbin, ybin, yerr, step_size = utils.rebin_data(self.x, self.y, dx_new,
+                                                       self.yerr)
         assert step_size == dx_new/self.dx
 
     def test_arrays(self):
         dx_new = 2.0
-        xbin, ybin, step_size = utils.rebin_data(self.x, self.y, dx_new)
+        xbin, ybin, yerr, step_size = utils.rebin_data(self.x, self.y, dx_new,
+                                                       self.yerr)
         assert isinstance(xbin, np.ndarray)
         assert isinstance(ybin, np.ndarray)
 
     def test_length_matches(self):
         dx_new = 2.0
-        xbin, ybin, step_size = utils.rebin_data(self.x, self.y, dx_new)
+        xbin, ybin, yerr, step_size = utils.rebin_data(self.x, self.y, dx_new,
+                                                       self.yerr)
         assert xbin.shape[0] == ybin.shape[0]
 
     def test_binned_counts(self):
         dx_new = 2.0
 
-        xbin, ybin, step_size = utils.rebin_data(self.x, self.y, dx_new)
+        xbin, ybin, yerr, step_size = utils.rebin_data(self.x, self.y, dx_new,
+                                                       self.yerr)
 
         ybin_test = np.zeros_like(xbin) + self.counts*dx_new/self.dx
         assert np.allclose(ybin, ybin_test)
 
     def test_uneven_bins(self):
         dx_new = 1.5
-        xbin, ybin, step_size = utils.rebin_data(self.x, self.y, dx_new)
+        xbin, ybin, yerr, step_size = utils.rebin_data(self.x, self.y, dx_new,
+                                                       self.yerr)
         assert np.isclose(xbin[1]-xbin[0], dx_new)
 
     def test_uneven_binned_counts(self):
         dx_new = 1.5
-        xbin, ybin, step_size = utils.rebin_data(self.x, self.y, dx_new)
+        xbin, ybin, yerr, step_size = utils.rebin_data(self.x, self.y, dx_new,
+                                                       self.yerr)
         ybin_test = np.zeros_like(xbin) + self.counts*dx_new/self.dx
         assert np.allclose(ybin_test, ybin)
 
     def test_rebin_data_should_raise_error_when_method_is_different_than_allowed(self):
         dx_new = 2.0
         with pytest.raises(ValueError):
-            utils.rebin_data(self.x, self.y, dx_new,
+            utils.rebin_data(self.x, self.y, dx_new, self.yerr,
                              method='not_allowed_method')