original_docstring: 
      Fit model to data using a given celerite kernel. Computes the PSD and generates useful plots.
      
      Note: x, y, and yerr must by astropy Quantities with units!
      x: time [astropy unit quantity]
      y: data [astropy unit quantity]
      yerr: error on data [astropy unit quantity]
      kernel: celerite kernel
      init: 'minimize', 'differential_evolution', or array of user-specified (e.g. previous) initial conditions
      nburn: number of burn-in samples
      nsamp: number of production samples
      target_name: name of target to display in light curve legend
      color: color for plotting
      plot: whether to plot the result
      verbose: whether to print useful messages
      supress_warn: whether to supress warnings
      seed: seed for random number generator
      
      returns: gp, samples, fig (celerite GuassianProcess object, samples array, and figure [None if plot=False])
      
docstring: Fit model to data using a given celerite kernel. Computes the PSD and generates useful plots.
  Note: x, y, and yerr must by astropy Quantities with units.
  
  gp, samples, fig (celerite GuassianProcess object, samples array, and figure [None if plot=False])
docstring_tokens: Fit model to data using a given celerite kernel . Computes the PSD and generates useful plots . Note : x y and yerr must by astropy Quantities with units . gp samples fig ( celerite GuassianProcess object samples array and figure [ None if plot = False ] )
code: def fit_celerite(x, y, yerr, kernel, init="minimize", nburn=500, nsamp=2000, lamb=None, target_name=None, color="#ff7f0e", plot=True, verbose=True, supress_warn=False, seed=None):
      """
      Fit model to data using a given celerite kernel. Computes the PSD and generates useful plots.
      
      Note: x, y, and yerr must by astropy Quantities with units!
      x: time [astropy unit quantity]
      y: data [astropy unit quantity]
      yerr: error on data [astropy unit quantity]
      kernel: celerite kernel
      init: 'minimize', 'differential_evolution', or array of user-specified (e.g. previous) initial conditions
      nburn: number of burn-in samples
      nsamp: number of production samples
      target_name: name of target to display in light curve legend
      color: color for plotting
      plot: whether to plot the result
      verbose: whether to print useful messages
      supress_warn: whether to supress warnings
      seed: seed for random number generator
      
      returns: gp, samples, fig (celerite GuassianProcess object, samples array, and figure [None if plot=False])
      """
      
      # Set seed for reproducability
      np.random.seed(seed)
      
      if supress_warn:
          warnings.filterwarnings("ignore")
      
      gp = celerite.GP(kernel, mean=np.mean(y.value), fit_mean=False)
      gp.compute(x.value, yerr.value)
      if verbose:
          print("Initial log-likelihood: {0}".format(gp.log_likelihood(y.value)))
  
      # Define a cost function
      def neg_log_like(params, y, gp):
          gp.set_parameter_vector(params)
          return -gp.log_likelihood(y)
  
      def grad_neg_log_like(params, y, gp):
          gp.set_parameter_vector(params)
          return -gp.grad_log_likelihood(y)[1]
  
      # Fit for the maximum likelihood parameters
      initial_params = gp.get_parameter_vector()
      bounds = gp.get_parameter_bounds()
      
      def _mle(y, gp, initial_params, bounds):
      
          # MLE solution
          if init == "minimize":
              soln = minimize(neg_log_like, initial_params, jac=grad_neg_log_like,
                              method="L-BFGS-B", bounds=bounds, args=(y.value, gp))
              initial = np.array(soln.x)
              if verbose:
                  print("Final log-likelihood: {0}".format(-soln.fun))
          elif init == "differential_evolution":
              soln = differential_evolution(neg_log_like, bounds=bounds, args=(y.value, gp))
              initial = np.array(soln.x)
              if verbose:
                  print("Final log-likelihood: {0}".format(-soln.fun))
          # Use user-provided initial MLE conditions
          elif np.issubdtype(np.array(init).dtype, np.number):
              initial = init
          else:
              raise ValueError('initial value not recognized!')
              
          return initial
      
      # Find MLE
      initial = _mle(y, gp, initial_params, bounds)   
      gp.set_parameter_vector(initial)
      
      # Filter long-term trends
      if lamb is not None:
          import statsmodels.api as sm
          # Filter on evenly-sampled MLE solution
          mean_cadence = np.mean(np.diff(x.value))
          t = np.arange(np.min(x.value), np.max(x.value), mean_cadence/10)
          mu = gp.predict(y.value, t, return_cov=False)
          cycle, trend = sm.tsa.filters.hpfilter(mu, lamb)
          
          if plot:
              # Light curve & prediction
              fig, ax_lc = plt.subplots(1,1, figsize=(9,5))
              ax_lc.errorbar(x.value, y.value, yerr=yerr.value, c='k', fmt='.', alpha=0.75, elinewidth=1, label=target_name)
              ax_lc.plot(t, trend, c='r', label='trend')
              ax_lc.invert_yaxis()
              ax_lc.set_xlabel("Time (days)", fontsize=20)
              if y.unit == u.mag:
                  ax_lc.set_ylabel("Magnitude", fontsize=20)
              else:
                  ax_lc.set_ylabel("Normalized Flux", fontsize=20)
  
              ax_lc.minorticks_on()
              ax_lc.tick_params('both',labelsize=18)
              ax_lc.tick_params(axis='both', which='both', direction='in')
              ax_lc.tick_params(axis='both', which='major', length=6)
              ax_lc.tick_params(axis='both', which='minor', length=3)
              ax_lc.xaxis.set_ticks_position('both')
              ax_lc.yaxis.set_ticks_position('both')
  
              ax_lc.legend(fontsize=16, loc=1)
              fig.tight_layout()
  
          # Subtract trend at real data
          y = y - np.interp(x.value, t, trend)*(y.unit) + np.median(y)
          
          # Find new MLE
          gp = celerite.GP(kernel, mean=np.mean(y.value), fit_mean=False)
          gp.compute(x.value, yerr.value)
          # Fit for the maximum likelihood parameters
          initial_params = gp.get_parameter_vector()
          bounds = gp.get_parameter_bounds()
      
          initial = _mle(y, gp, initial_params, bounds)   
          gp.set_parameter_vector(initial)
          
      
      # Define the log probablity
      def log_probability(params):
          gp.set_parameter_vector(params)
          lp = gp.log_prior()
          if not np.isfinite(lp):
              return -np.inf
          return gp.log_likelihood(y) + lp
  
      ndim, nwalkers = len(initial), 32
      sampler = emcee.EnsembleSampler(nwalkers, ndim, log_probability)
      
      if verbose:
          print("Running burn-in...")
      p0 = initial + 1e-8 * np.random.randn(nwalkers, ndim)
      p0, lp, _ = sampler.run_mcmc(p0, nburn)
      
      if verbose:
          print("Running production...")
      sampler.reset()
      sampler.run_mcmc(p0, nsamp);
          
      # Get posterior and uncertianty
      samples = sampler.flatchain
      
      s = np.median(samples, axis=0)
      gp.set_parameter_vector(s)
      
      if plot:
          fig = plot_celerite(x, y, yerr, gp, samples, target_name=target_name, color=color)
      else:
          fig = None
      
      return gp, samples, fig
code_tokens: def fit_celerite ( x , y , yerr , kernel , init = "minimize" , nburn = 500 , nsamp = 2000 , lamb = None , target_name = None , color = "#ff7f0e" , plot = True , verbose = True , supress_warn = False , seed = None ) : np . random . seed ( seed ) if supress_warn : warnings . filterwarnings ( "ignore" ) gp = celerite . GP ( kernel , mean = np . mean ( y . value ) , fit_mean = False ) gp . compute ( x . value , yerr . value ) if verbose : print ( "Initial log-likelihood: {0}" . format ( gp . log_likelihood ( y . value ) ) ) def neg_log_like ( params , y , gp ) : gp . set_parameter_vector ( params ) return - gp . log_likelihood ( y ) def grad_neg_log_like ( params , y , gp ) : gp . set_parameter_vector ( params ) return - gp . grad_log_likelihood ( y ) [ 1 ] initial_params = gp . get_parameter_vector ( ) bounds = gp . get_parameter_bounds ( ) def _mle ( y , gp , initial_params , bounds ) : if init == "minimize" : soln = minimize ( neg_log_like , initial_params , jac = grad_neg_log_like , method = "L-BFGS-B" , bounds = bounds , args = ( y . value , gp ) ) initial = np . array ( soln . x ) if verbose : print ( "Final log-likelihood: {0}" . format ( - soln . fun ) ) elif init == "differential_evolution" : soln = differential_evolution ( neg_log_like , bounds = bounds , args = ( y . value , gp ) ) initial = np . array ( soln . x ) if verbose : print ( "Final log-likelihood: {0}" . format ( - soln . fun ) ) elif np . issubdtype ( np . array ( init ) . dtype , np . number ) : initial = init else : raise ValueError ( 'initial value not recognized!' ) return initial initial = _mle ( y , gp , initial_params , bounds ) gp . set_parameter_vector ( initial ) if lamb is not None : import statsmodels . api as sm mean_cadence = np . mean ( np . diff ( x . value ) ) t = np . arange ( np . min ( x . value ) , np . max ( x . value ) , mean_cadence / 10 ) mu = gp . predict ( y . value , t , return_cov = False ) cycle , trend = sm . tsa . filters . hpfilter ( mu , lamb ) if plot : fig , ax_lc = plt . subplots ( 1 , 1 , figsize = ( 9 , 5 ) ) ax_lc . errorbar ( x . value , y . value , yerr = yerr . value , c = 'k' , fmt = '.' , alpha = 0.75 , elinewidth = 1 , label = target_name ) ax_lc . plot ( t , trend , c = 'r' , label = 'trend' ) ax_lc . invert_yaxis ( ) ax_lc . set_xlabel ( "Time (days)" , fontsize = 20 ) if y . unit == u . mag : ax_lc . set_ylabel ( "Magnitude" , fontsize = 20 ) else : ax_lc . set_ylabel ( "Normalized Flux" , fontsize = 20 ) ax_lc . minorticks_on ( ) ax_lc . tick_params ( 'both' , labelsize = 18 ) ax_lc . tick_params ( axis = 'both' , which = 'both' , direction = 'in' ) ax_lc . tick_params ( axis = 'both' , which = 'major' , length = 6 ) ax_lc . tick_params ( axis = 'both' , which = 'minor' , length = 3 ) ax_lc . xaxis . set_ticks_position ( 'both' ) ax_lc . yaxis . set_ticks_position ( 'both' ) ax_lc . legend ( fontsize = 16 , loc = 1 ) fig . tight_layout ( ) y = y - np . interp ( x . value , t , trend ) * ( y . unit ) + np . median ( y ) gp = celerite . GP ( kernel , mean = np . mean ( y . value ) , fit_mean = False ) gp . compute ( x . value , yerr . value ) initial_params = gp . get_parameter_vector ( ) bounds = gp . get_parameter_bounds ( ) initial = _mle ( y , gp , initial_params , bounds ) gp . set_parameter_vector ( initial ) def log_probability ( params ) : gp . set_parameter_vector ( params ) lp = gp . log_prior ( ) if not np . isfinite ( lp ) : return - np . inf return gp . log_likelihood ( y ) + lp ndim , nwalkers = len ( initial ) , 32 sampler = emcee . EnsembleSampler ( nwalkers , ndim , log_probability ) if verbose : print ( "Running burn-in..." ) p0 = initial + 1e-8 * np . random . randn ( nwalkers , ndim ) p0 , lp , _ = sampler . run_mcmc ( p0 , nburn ) if verbose : print ( "Running production..." ) sampler . reset ( ) sampler . run_mcmc ( p0 , nsamp ) ; samples = sampler . flatchain s = np . median ( samples , axis = 0 ) gp . set_parameter_vector ( s ) if plot : fig = plot_celerite ( x , y , yerr , gp , samples , target_name = target_name , color = color ) else : fig = None return gp , samples , fig
short_docstring: Fit model to data using a given celerite kernel.
short_docstring_tokens: Fit model to data using a given celerite kernel .
comment: """
      Fit model to data using a given celerite kernel. Computes the PSD and generates useful plots.
      
      Note: x, y, and yerr must by astropy Quantities with units!
      x: time [astropy unit quantity]
      y: data [astropy unit quantity]
      yerr: error on data [astropy unit quantity]
      kernel: celerite kernel
      init: 'minimize', 'differential_evolution', or array of user-specified (e.g. previous) initial conditions
      nburn: number of burn-in samples
      nsamp: number of production samples
      target_name: name of target to display in light curve legend
      color: color for plotting
      plot: whether to plot the result
      verbose: whether to print useful messages
      supress_warn: whether to supress warnings
      seed: seed for random number generator
      
      returns: gp, samples, fig (celerite GuassianProcess object, samples array, and figure [None if plot=False])
      """$SEP$# Set seed for reproducability$SEP$# Define a cost function$SEP$# Fit for the maximum likelihood parameters$SEP$# MLE solution$SEP$# Use user-provided initial MLE conditions$SEP$# Find MLE$SEP$# Filter long-term trends$SEP$# Filter on evenly-sampled MLE solution$SEP$# Light curve & prediction$SEP$# Subtract trend at real data$SEP$# Find new MLE$SEP$# Fit for the maximum likelihood parameters$SEP$# Define the log probablity$SEP$# Get posterior and uncertianty