def fig3(x):
  # Figure 3
  # now we create a cumulative histogram of the data
  #
  P.figure()

  n, bins, patches = P.hist(x, 50, normed=1, histtype='step', cumulative=True)

  # add a line showing the expected distribution
  y = P.normpdf(bins, mu, sigma).cumsum()
  y /= y[-1]
  l = P.plot(bins, y, 'k--', linewidth=1.5)

  # create a second data-set with a smaller standard deviation
  sigma2 = 15.
  x = mu + sigma2 * P.randn(10000)

  n, bins, patches = P.hist(x, bins=bins, normed=1, histtype='step', cumulative=True)

  # add a line showing the expected distribution
  y = P.normpdf(bins, mu, sigma2).cumsum()
  y /= y[-1]
  l = P.plot(bins, y, 'r--', linewidth=1.5)

  # finally overplot a reverted cumulative histogram
  n, bins, patches = P.hist(x, bins=bins, normed=1,
                            histtype='step', cumulative=-1)

  P.grid(True)
  P.ylim(0, 1.05)

    def plot(self, normed=True, N=1000, Xmin=None, Xmax=None, bins=50, color='red', lw=2, 
            hist_kw={'color':'#5F9EA0'}, ax=None):

        if ax:
            ax.hist(self.data, normed=normed, bins=bins, **hist_kw)
        else:
            pylab.hist(self.data, normed=normed, bins=bins, **hist_kw)
        if Xmin is None:
            Xmin = self.data.min()
        if Xmax is None:
            Xmax = self.data.max()
        X = pylab.linspace(Xmin, Xmax, N)

        if ax:
            ax.plot(X, [self.model.pdf(x, self.results.x) for x in X], color=color, lw=lw)
        else:
            pylab.plot(X, [self.model.pdf(x, self.results.x) for x in X], color=color, lw=lw)

        K = len(self.results.x)
        # The PIs must be normalised
        for i in range(0, K/3):
            
            mu, sigma, pi_ = self.results.mus[i], self.results.sigmas[i], self.results.pis[i]
            if ax:
                ax.plot(X, [pi_ * pylab.normpdf(x, mu, sigma) for x in X], 'g--', alpha=0.5)
            else:
                pylab.plot(X, [pi_ * pylab.normpdf(x, mu, sigma) for x in X], 'g--', alpha=0.5)

def extraHistogram():
    import pylab as P
    import numpy as N
    
    
    P.figure()
    bins = 25
    min = 23310.
    max = 23455.
    nu, binsu, patchesu = P.hist(telfocusOld, bins=bins, rwidth=0.9, range = (min, max) , 
           label = 'UnCorrected Data', fc ='b', alpha = 0.4, normed = True)
    n, bins, patches = P.hist(telfocusCorrected, bins=bins, rwidth=0.7, range = (min, max),
           label='Corrected Data', fc = 'r', alpha = 0.6, normed = True)
    #P.axvline(medianNew, label = 'Corrected Median', color = 'r', lw = 1.1)
    #P.axvline(medianOld, label = 'UnCorrected Median', color = 'b', lw = 1.1)
    y1 = P.normpdf(binsu, N.mean(telfocusOld), N.std(telfocusOld))
    y2 = P.normpdf(bins, N.mean(telfocusCorrected), N.std(telfocusCorrected))
    P.plot(binsu, y1, 'b-', linewidth = 2.5, label='Gaussian Fit')
    P.plot(bins, y2, 'r-', linewidth = 3., label='Gaussian Fit')
    P.xlim(min,max)
    P.xlabel('Telescope Focus + median Offset')
    P.ylabel('Normed Values')
    P.legend(shadow=True, loc ='best')
    P.savefig('TelFocusHistogram.png')
    P.close()

    def pdf(self, x, params, normalise=True):
        """Expected parameters are


        params is a list of gaussian distribution ordered as mu, sigma, pi, 
        mu2, sigma2, pi2, ...

        """
        assert divmod(len(params), 3)[1] == 0
        assert len(params) >= 3 * self.k
        k = len(params) / 3

        self.k = k

        pis = np.array(params[2::3])

        if any(np.array(pis)<0):
            return 0
        if normalise is True:
            pis /= pis.sum()
        # !!! sum pi must equal 1 otherwise may diverge badly
        data = 0
        for i in range(0, k):
            mu, sigma, pi_ = params[i*3: (i+1)*3]
            pi_ = pis[i]
            if sigma != 0:
                data += pi_ * pylab.normpdf(x, mu, sigma)
        return data

def main():
	''' Main Function'''
	# Start and End date of the charts
	dt_start = dt.datetime(2011, 1, 1)
	dt_end = dt.datetime(2012, 12, 31)
	
	#goog = pd.io.data.get_data_yahoo("GOOG",  dt_start, dt_end) # not working
	SPY = DataReader("SPY",  "yahoo", dt_start, dt_end)
	#YHOO = DataReader("YHOO",  "yahoo", dt_start, dt_end)
	
	# normalize prices
	nPrice = sc.normalizedPrice(SPY['Adj Close'].values)
	
	#daily return
	daily_ret = sc.computeDailyReturn(nPrice)
	plt.subplot(1,2,1)
	plt.plot(daily_ret*100, 'b-')
	plt.ylabel('Daily return (%)')
	plt.legend(['SPY-Daily Return based on Adjuested close'])
	
	#daily return histogram
	plt.subplot(1,2,2)
	n, bins, patches = plt.hist(daily_ret, 100, normed=1, facecolor='green', alpha=0.5)
	mean = np.mean(daily_ret)
	sigma = np.std(daily_ret)
	y = P.normpdf( bins, mean, sigma)
	plt.plot(bins, y, 'k--', linewidth=1.5)

	plt.ylabel('Probability')
	plt.legend(['normal distribution approximation','SPY-hostogram of daily return'])
	
	plt.show()

def plot_histogram():
    import numpy as np
    import pylab as P

    # The hist() function now has a lot more options
    #

    #
    # first create a single histogram
    #
    P.figure()
    mu, sigma = 40, 35


    x = abs(np.random.normal(mu, sigma, 1000000))

    # the histogram of the data with histtype='step'
    n, bins, patches = P.hist(x, 100, normed=1, histtype='stepfilled')
    P.setp(patches, 'facecolor', 'g', 'alpha', 0.50)
    P.vlines(np.mean(x), 0, max(n))
    P.vlines(np.median(x), 0, max(n))
    # add a line showing the expected distribution
    y = np.abs(P.normpdf( bins, mu, sigma))
    l = P.plot(bins, y, 'k--', linewidth=1.5)

    P.show()

def load_spectrum_gnirs(file, velScale, resolution):
    """
    Load up a spectrum from the Gemini GNIRS library.
    """
    spec, hdr = pyfits.getdata(file, header=True)

    pixScale = hdr['CD1_1'] * 1.0e-4    # microns

    wavelength = np.arange(len(spec), dtype=float)
    wavelength -= (hdr['CRPIX1']-1.0)   # get into pixels relative to the reference pix
    wavelength *= hdr['CD1_1']          # convert to the proper wavelength scale (Ang)
    wavelength += hdr['CRVAL1']         # shift to the wavelength zeropoint

    # Convert from Angstroms to microns
    wavelength *= 1.0e-4

    deltaWave = 2.21344 / resolution         # microns
    resInPixels = deltaWave / pixScale       # pixels
    sigmaInPixels = resInPixels / 2.355
    psfBins = np.arange(-4*math.ceil(resInPixels), 4*math.ceil(resInPixels))
    psf = py.normpdf(psfBins, 0, sigmaInPixels)
    specLowRes = np.convolve(spec, psf, mode='same')

    # Rebin into equals logarithmic steps in wavelength
    logWave, specNew, vel = log_rebin2(wavelength, specLowRes,
                                       inMicrons=True, velScale=velScale)

    return logWave, specNew

def complicatedHisto(x, dists, var, mean):
    
    '''
    this histogram will take a 2d x with each column being a different
    set of data and then plot each one with its own color.. like in one of the
    histogram matplotlib examples
    
    '''
    fig = plt.figure()
    dist = ['{:.2f}'.format(k) for k in dists]    
    colors = ('g', 'b', 'r','c', 'm','y', 'k','w')    

    n, bins, patches = P.hist(x, 20, normed=1, histtype='bar',
                            color=colors[0:len(dist)],
                            label=dist)
                            
    P.legend()
    P.show()
    
    print "woohoo"
    fig2 = plt.figure(2)
    xbins = np.linspace(-0.5,3,200)
    for sigma, mu,d,c  in zip(var, mean, dists,colors):
        
        y = P.normpdf(xbins,  mu, np.sqrt(sigma))
        l = P.plot(xbins, y, c, label=d, linewidth=1.5)
        P.legend()

def build_hist(ax, dat, title=None) :
    mu, sigma = np.average(dat), np.std(dat)
    n, bins, patches = ax.hist(dat, 50, normed=1, histtype='stepfilled')
    P.setp(patches, 'facecolor', 'g', 'alpha', 0.75)

    y = P.normpdf(bins, mu, sigma)
    l = ax.plot(bins, y, 'k--', linewidth=1.5)
    if title != None :
        ax.set_title(title)

 def histDemo(self):
     mu, sigma = 200, 25
     x = mu + sigma*P.randn(10000)
     n, bins, patches = P.hist(x, 50, normed=1, histtype='stepfilled')
     P.setp(patches, 'facecolor', 'g', 'alpha', 0.75)
     
     y = P.normpdf( bins, mu, sigma)
     l = P.plot(bins, y, 'k--', linewidth=1.5)
     P.savefig(inspect.stack()[0][3]+".png")

def analysis(series_dir):
    files = glob(series_dir + "/*.dicm")
    files.sort()
    series_intensity = [imageIntensity(name) for name in files]
    y = normpdf(series_intensity, min(series_intensity), max(series_intensity))
    plot([i for i in range(len(series_intensity))], series_intensity, 'b')
    plot([i for i in range(len(series_intensity))], deriv2(series_intensity), 'r')
    grid(True)
    show()
    return series_intensity

 def pdf_model(x, p):
     print
     print "pdf_model()"
     print "  x=%s" % x
     print "  p=%s" % (p,)
     mu1, sig1, mu2, sig2, pi_1 = p
     print "  mu1:  %s" % mu1
     print "  sig1: %s" % sig1
     print "  mu2:  %s" % mu2
     print "  sig2: %s" % sig2
     print "  pi_1: %s" % pi_1
     raw1 = py.normpdf(x, mu1, sig1)
     print "  raw1: %s" % raw1
     raw2 = py.normpdf(x, mu2, sig2)
     print "  raw2: %s" % raw2
     ret = pi_1 * raw1 + (1 - pi_1) * raw2
     print "  ret: %s" % ret
     print
     return ret

def stability(paths, show=False, output=None, annotations=None, aspect=2):
  # COLOR = "#548BE3"
  COLOR = "#8CB8FF"
  figure = p.figure()
  pltnum = len(paths)//2 + len(paths)%2
  for i, fname in enumerate(paths):
    with open(fname) as fd:
      values = []
      t = fd.readline().strip("#")
      for l in fd.readlines():
        values += [float(l.strip())]
      p.subplot(pltnum, 2, i+1)

      ## title
      if annotations:
        p.title(annotations[i])
      else:
        p.title(t)
      avg = average(values)
      percents = list(map(lambda x: avg/x-1, values))
      n, bins, patches = p.hist(percents,
        bins=50, normed=True,
        histtype='bar', color=COLOR)

      mu, sigma = norm.fit(percents)
      y = p.normpdf(bins, mu, sigma)
      p.plot(bins, y, 'r-', linewidth=3)
      p.xlim(min(bins), max(bins))

      ## set aspect
      xmin, xmax = p.xlim()
      ymin, ymax = p.ylim()
      xr = xmax - xmin
      yr = ymax - ymin
      aspect = xr/yr/2
      g = p.gca()
      g.set_aspect(aspect)
      # p.figaspect(aspect)


    ## remove y axis
    yaxis = p.gca().yaxis
    yaxis.set_major_locator(MaxNLocator(nbins=4, prune='lower'))
    # yaxis.set_visible(False)

    ## xaxis
    xaxis = p.gca().xaxis
    xaxis.set_major_formatter(to_percent)
    xaxis.set_major_locator(MaxNLocator(nbins=5, prune='lower'))

  p.tight_layout(pad=0.5)
  if output:
    p.savefig(output, bbox_inches='tight')
  if show:
    p.show()

def fig1(x):
  # Figure 1
  # first create a single histogram
  #
  # the histogram of the data with histtype='step'
  P.figure()
  n, bins, patches = P.hist(x, 50, normed=1, histtype='stepfilled')
  P.setp(patches, 'facecolor', 'g', 'alpha', 0.75)

  # add a line showing the expected distribution
  y = P.normpdf(bins, mu, sigma)
  l = P.plot(bins, y, 'k--', linewidth=1.5)

def densityplot(data):
    """
    Plots a histogram of daily returns from data, plus fitted normal density.
    """
    dailyreturns = percent_change(data)
    pylab.hist(dailyreturns, bins=200, normed=True)
    m, M = min(dailyreturns), max(dailyreturns)
    mu = pylab.mean(dailyreturns)
    sigma = pylab.std(dailyreturns)
    grid = pylab.linspace(m, M, 100)
    densityvalues = pylab.normpdf(grid, mu, sigma)
    pylab.plot(grid, densityvalues, 'r-')
    pylab.show()

def probdist1(l):
    '''This function creates an array of random numbers from Gaussian 
    distribution, and then plots a histogram of those random numbers and 
    the exptected probability distribution computed solely from the mean
    and variance of the data points.'''
    lstrand = []
    for i in range(l):
        lstrand.append(numpy.random.randn());
    mu = numpy.mean(lstrand)
    sigma = math.sqrt(numpy.var(lstrand))
    n, bins, patches = pylab.hist(lstrand, 100, normed=1)
    pylab.setp(patches, 'facecolor', 'g', 'alpha', 0.75)
    y = pylab.normpdf(bins, mu, sigma)
    line = pylab.plot(bins, y, 'k--', linewidth=1.5)
    pylab.show()

def main():
    df = pd.read_csv(INDATA)
    df['basin_id'] = np.arange(15006)+1
    df.set_index('basin_id', inplace=True)
    wdf = pd.read_csv(AREA)
    wdf.set_index('basin_id', inplace=True)
    df.join(wdf, inplace=True)
    df.dropna(inplace=True)
    x = np.asarray(df["2040"])
    w = np.asarray(df["F_AREA"])
    sd = np.std(x)
    n, bins, patches = pl.hist(x, 50, normed=1)
    y = pl.normpdf(bins, 0, sd)
    l = pl.plot(bins, y, 'k--', linewidth=1.5)
    pl.show()

def plot_histogram(fig):
    mu, sigma = 100, 15
    x = mu + sigma*randn(10000)

    axes = fig.gca()
    # the histogram of the data
    n, bins, patches = axes.hist(x, 100, normed=1)

    # add a 'best fit' line
    y = normpdf( bins, mu, sigma)
    l = axes.plot(bins, y, 'r--', linewidth=2)

    axes.set_xlim((40, 160))
    axes.set_xlabel('Smarts')
    axes.set_ylabel('P')
    axes.set_title('IQ: mu=100, sigma=15')

def plot_spectrum_componenten():
    """Plot voor componenten pulshoogte spectrum"""

    multiplot = MultiPlot(3, 1, axis='semilogy', width=r'0.5\linewidth')

    pylab.clf()
    subplot0 = multiplot.get_subplot_at(0, 0)
    subplot1 = multiplot.get_subplot_at(1, 0)
    subplot2 = multiplot.get_subplot_at(2, 0)
    x = pylab.linspace(1e-10, 11, 200)
    th_signal = []
    signal = pylab.zeros(x.shape)
    for N in range(1, 15):
        scale = 100. / N ** 3
        pdf = 80 * scale * pylab.normpdf(x, N, pylab.sqrt(N) * 0.35)
        # pylab.plot(x, pdf)
        subplot1.plot(x, pdf, mark=None)
        # subplot1.add_pin('%d MIP' % N, 'above right', x=N, use_arrow=True,
        #                  style='lightgray')
        subplot2.plot(x, pdf, mark=None, linestyle='lightgray')
        signal += pdf
        th_signal.extend(int(100 * scale) * [N])

    gammas = 1e2 * x ** -3
    subplot1.plot(x, gammas, mark=None)
    subplot2.plot(x, gammas, mark=None, linestyle='lightgray')

    signal += gammas
    pylab.plot(x, signal)
    pylab.plot(x, gammas)
    subplot2.plot(x, signal, mark=None)
    pylab.yscale('log')
    pylab.ylim(ymin=1)

    n, bins = pylab.histogram(th_signal, bins=pylab.linspace(0, 11, 100))
    n = pylab.where(n == 0, 1e-10, n)
    subplot0.histogram(n, bins)

    multiplot.show_xticklabels_for_all([(0, 0), (2, 0)])
    multiplot.set_xticks_for_all([(0, 0), (2, 0)], range(20))
    # multiplot.show_yticklabels_for_all([(0, 0), (1, 0), (2, 0)])
    multiplot.set_xlabel('Aantal deeltjes')
    multiplot.set_ylabel('Aantal events')
    multiplot.set_ylimits_for_all(min=1, max=1e5)
    multiplot.set_xlimits_for_all(min=0, max=10.5)

    multiplot.save('spectrum_componenten')

	def test_NormalByCauchyAcceptReject(self, no_of_samples=10000):
		import rpy, math, random, pylab
		sample_list = []
		M = math.pi/(math.sqrt(2*math.pi)) +0.1 #"+0.1" makes it a little bit less efficient, but keeps the density around 0 normal
		for i in range(no_of_samples):
			cauchy_sample = rpy.r.rcauchy(1)
			u = random.random()
			normal_mass = rpy.r.dnorm(cauchy_sample)
			cauchy_mass = 1/(math.pi*(1+cauchy_sample*cauchy_sample))
			if u<= normal_mass/(M*cauchy_mass):
				sample_list.append(cauchy_sample)
		accept_ratio = len(sample_list)/float(no_of_samples)
		print "accept ratio is %s"%(accept_ratio)
		n, bins, patches = pylab.hist(sample_list, 100, normed=1)
		pylab.title("Normal via cauchy, accept_ratio:%s"%accept_ratio)
		y = pylab.normpdf( bins, 0, 1)
		l = pylab.plot(bins, y, 'ro-', linewidth=2)
		pylab.show()