Package ivs :: Package sigproc :: Module evaluate

[frames] | no frames]

Module evaluate

Evaluate model functions.

Most functions accept a domain of interest and a set of parameters, and return a model function (e.g., a sine, gaussian etc...). Parameters are usually a record array containing named fields for the parameters.

Some nonlinear optimizers expect and return parameters as a 1D array of the independent parameters. Therefore, all functions also accept these 1D numpy arrays, and they are converted by check_input to record arrays. The latter function actually accepts both record arrays and 1D arrays, and acts as a translator between the two parameter representations.

Functions

check_input(fctn)
Decorator to check input of evaluating model, and transforming input from record array to 1D numpy array for optimization, and vice versa.

source code

Evaluating model fits

float

RSS(data, model)
Residual sum of squares.

source code

percentage

varred(data, model)
Calculate variance reduction.

source code

float

bic(data, model, k=0, n=None)
BIC for time regression

source code

float

aic(data, model, k=0, n=None)
AIC for time regression

source code

ndarray (shape n)

fishers_method(pvalues)
Use Fisher's test to combine probabilities.

source code

Timeseries

ndarray,ndarray(, ndarray)

phasediagram(time, signal, nu0, D=0, t0=None, forb=None, asini=None, return_indices=False, chronological=False)
Construct a phasediagram, using frequency nu0.

source code

array

sine(times, parameters)
Creates a harmonic function based on given parameters.

source code

array

sine_freqshift(times, parameters, t0=None)
Creates a sine function with a linear frequency shift.

source code

array

sine_orbit(times, parameters, t0=None, nmax=10)
Creates a sine function with a sinusoidal frequency shift.

source code

array

periodic_spline(times, parameters, t0=None)
Evaluate a periodic spline function

source code

array

kepler(times, parameters, itermax=8)
Construct a radial velocity curve due to Kepler orbit(s).

source code

kepler_diffcorr(times, parameters, itermax=8)
Compute coefficients for differential correction method.

source code

binary(times, parameters, n1=None, itermax=8)
Simultaneous evaluation of two binary components.

source code

array

box(times, parameters, t0=None)
Evaluate a box transit model.

source code

spots(times, spots, t0=None, **kwargs)
Spotted model from Lanza (2003).

source code

General purpose

array

gauss(x, parameters)
Evaluate a gaussian fit.

source code

lorentz(x, parameters)
Evaluate Lorentz profile

source code

array

voigt(x, parameters)
Evaluate a Voigt profile.

source code

array

power_law(x, parameters)
Evaluate a power law.

source code

Convert record arrays to flat arrays and back

normal numpy array or record array

sine_preppars(pars)
Prepare sine function parameters in correct form for evaluating/fitting.

source code

normal numpy array or record array

kepler_preppars(pars)
Prepare Kepler orbit parameters in correct form for evaluating/fitting.

source code

normal numpy array or record array

binary_preppars(pars)
Prepare Kepler orbit parameters in correct form for evaluating/fitting.

source code

normal numpy array or record array

box_preppars(pars)
Prepare sine function parameters in correct form for evaluating/fitting.

source code

normal numpy array or record array

gauss_preppars(pars)
Prepare gauss function parameters in correct form for evaluating/fitting.

source code

normal numpy array or record array

lorentz_preppars(pars)
Prepare Lorentz function parameters in correct form for evaluating/fitting.

source code

normal numpy array or record array

voigt_preppars(pars)
Prepare voigt function parameters in correct form for evaluating/fitting.

source code

normal numpy array or record array

power_law_preppars(pars)
Prepare gauss function parameters in correct form for evaluating/fitting.

source code

Helper functions

_complex_error_function(x)
Complex error function

source code

Variables

logger = logging.getLogger('SIGPROC.EVAL')

Function Details

RSS(data, model)

Residual sum of squares.

Parameters:

data (numpy array) - data
model (numpy array) - model

Returns: float

residual sum of squares

varred(data, model)

Calculate variance reduction.

Parameters:

data (numpy array) - data
model (numpy array) - model

Returns: percentage

variance reduction (percentage)

bic(data, model, k=0, n=None)

BIC for time regression

This one is based on the MLE of the variance,

Sigma_k^2 = RSS/n

Where our maximum likelihood estimator is then (RSS/n)^n

Parameters:

data (numpy array) - data
model (numpy array) - model
k (integer) - number of free, independent parameters
n (integer or None) - number of effective observations

Returns: float

BIC

aic(data, model, k=0, n=None)

AIC for time regression

This one is based on the MLE of the variance,

Sigma_k^2 = RSS/n

Where our maximum likelihood estimator is then (RSS/n)^n

Parameters:

data (numpy array) - data
model (numpy array) - model
k (integer) - number of free, independent parameters
n (integer or None) - number of effective observations

Returns: float

AIC

fishers_method(pvalues)

Use Fisher's test to combine probabilities.

http://en.wikipedia.org/wiki/Fisher's_method

>>> rvs1 = distributions.norm().rvs(10000)
>>> rvs2 = distributions.norm().rvs(10000)
>>> pvals1 = scipy.stats.distributions.norm.sf(rvs1)
>>> pvals2 = scipy.stats.distributions.norm.sf(rvs2)
>>> fpval = fishers_method([pvals1,pvals2])

And make a plot (compare with the wiki page).

>>> p = pl.figure()
>>> p = pl.scatter(pvals1,pvals2,c=fpval,edgecolors='none',cmap=pl.cm.spectral)
>>> p = pl.colorbar()

Parameters:

pvalues (ndarray (shape Nprob times n)) - array of pvalues, the first axis will be combined

Returns: ndarray (shape n)

combined pvalues

phasediagram(time, signal, nu0, D=0, t0=None, forb=None, asini=None, return_indices=False, chronological=False)

Construct a phasediagram, using frequency nu0.

Possibility to include a linear frequency shift D.

If needed, the indices can be returned. To 'unphase', just do:

Example usage:

>>> original = np.random.normal(size=10)
>>> indices = np.argsort(original)
>>> inv_indices = np.argsort(indices)
>>> all(original == np.take(np.take(original,indices),inv_indices))
True

Optionally, the zero time point can be given, it will be subtracted from all time points

Joris De Ridder, Pieter Degroote

Parameters:

time (ndarray) - time points [0..Ntime-1]
signal (ndarray) - observed data points [0..Ntime-1]
nu0 (float) - frequency to compute the phasediagram with (inverse unit of time)
t0 (float) - zero time point, if None: t0 = time[0]
D (float) - frequency shift
chronological (boolean) - set to True if you want to have a list of every phase separately
return_indices (boolean) - flag to return indices

Returns: ndarray,ndarray(, ndarray)

phase points (sorted), corresponding observations

sine(times, parameters)

Creates a harmonic function based on given parameters.

Parameter fields: const, ampl, freq, phase.

This harmonic function is of the form

f(p1,...,pm;x1,...,x_n) = ∑ c_i + a_i sin(2 π (f_i+ φ_i))

where p_i are parameters and x_i are observation times.

Parameters can be given as a record array containing columns 'ampl', 'freq' and 'phase', and optionally also 'const' (the latter array is summed up so to reduce it to one number).

pars = np.rec.fromarrays([consts,ampls,freqs,phases],names=('const','ampl','freq','phase'))

Parameters can also be given as normal 1D numpy array. In that case, it will be transformed into a record array by the sine_preppars function. The array you give must be 1D, optionally contain a constant as the first parameter, followed by 3 three numbers per frequency:

pars = [const, ampl1, freq1, phase1, ampl2, freq2, phase2...]

Example usage: We construct a synthetic signal with two frequencies.

>>> times = np.linspace(0,100,1000)
>>> const = 4.
>>> ampls = [0.01,0.004]
>>> freqs = [0.1,0.234]
>>> phases = [0.123,0.234]

The signal can be made via a record array

>>> parameters = np.rec.fromarrays([[const,0],ampls,freqs,phases],names = ('const','ampl','freq','phase'))
>>> mysine = sine(times,parameters)

or a normal 1D numpy array

>>> parameters = np.hstack([const,np.ravel(np.column_stack([ampls,freqs,phases]))])
>>> mysine = sine(times,parameters)

Of course, the latter is easier if you have your parameters in a list:

>>> freq1 = [0.01,0.1,0.123]
>>> freq2 = [0.004,0.234,0.234]
>>> parameters = np.hstack([freq1,freq2])
>>> mysine = sine(times,parameters)

Parameters:

times (numpy array) - observation times
parameters (record array or 1D array) - record array containing amplitudes ('ampl'), frequencies ('freq') phases ('phase') and optionally constants ('const'), or 1D numpy array (see above). Warning: record array should have shape (n,)!

Returns: array

sine signal (same shape as times)

Decorators:

@check_input

sine_freqshift(times, parameters, t0=None)

Creates a sine function with a linear frequency shift.

Parameter fields: const, ampl, freq, phase, D.

Similar to sine, but with extra column 'D', which is the linear frequency shift parameter.

Only accepts record arrays as parameters (for the moment).

Parameters:

times (numpy array) - observation times
parameters (record array) - record array containing amplitudes ('ampl'), frequencies ('freq'), phases ('phase'), frequency shifts ('D') and optionally constants ('const')

Returns: array

sine signal with frequency shift (same shape as times)

sine_orbit(times, parameters, t0=None, nmax=10)

Creates a sine function with a sinusoidal frequency shift.

Parameter fields: const, ampl, freq, phase, asini, omega.

Similar to sine, but with extra column 'asini' and 'forb', which are the orbital parameters. forb in cycles/day or something similar, asini in au.

For eccentric orbits, add longitude of periastron 'omega' (radians) and 'ecc' (eccentricity).

Only accepts record arrays as parameters (for the moment).

Parameters:

times (numpy array) - observation times
parameters (record array) - record array containing amplitudes ('ampl'), frequencies ('freq'), phases ('phase'), frequency shifts ('D') and optionally constants ('const')

Returns: array

sine signal with frequency shift (same shape as times)

periodic_spline(times, parameters, t0=None)

Evaluate a periodic spline function

Parameter fields: freq, knots, coeffs, degree, as returned from the corresponding fitting function fit.periodic_spline.

Parameters:

times (numpy array) - observation times
parameters - [frequency,"cj" spline coefficients]

Returns: array

sine signal with frequency shift (same shape as times)

kepler(times, parameters, itermax=8)

Construct a radial velocity curve due to Kepler orbit(s).

Parameter fields: gamma, P, T0, e, omega, K

Parameters:

times (numpy array) - observation times
parameters (record array) - record array containing periods ('P' in days), times of periastron ('T0' in days), eccentricities ('e'), longitudes of periastron ('omega' in radians), amplitude ('K' in km/s) and systemic velocity ('gamma' in km/s) phases ('phase'), frequency shifts ('D') and optionally constants ('const')
itermax (integer) - maximum number of iterations for solving Kepler equation

Returns: array

radial velocity curve (same shape as times)

Decorators:

@check_input

kepler_diffcorr(times, parameters, itermax=8)

Compute coefficients for differential correction method.

See page 97-98 of Hilditch's book "An introduction to close binary stars".

Decorators:

@check_input

binary(times, parameters, n1=None, itermax=8)

Simultaneous evaluation of two binary components.

parameter fields must be 'P','T0','e','omega','K1', 'K2'.

Decorators:

@check_input

box(times, parameters, t0=None)

Evaluate a box transit model.

Parameter fields: cont, freq, ingress, egress, depth

Parameters [[frequency,depth, fractional start, fraction end, continuum]]

Returns: array

Decorators:

@check_input

spots(times, spots, t0=None, **kwargs)


Spotted model from Lanza (2003).
Extract from Carrington Map.

Included:
    - presence of faculae through Q-factor (set constant to 5 or see 
    Chapman et al. 1997 or Solank & Unruh for a scaling law, since this
    factor decreases strongly with decreasing rotation period.
    - differential rotation through B-value (B/period = 0.2 for the sun)
    - limb-darkening through coefficients ap, bp and cp
    - non-equidistant time series
    - individual starspot evolution and lifetime following a Gaussian law,
    (see Hall & Henry (1993) for a relation between stellar radius and
    sunspot lifetime), through an optional 3rd (time of maximum) and
    4th parameter (decay time)

Not included:
    - distinction between umbra's and penumbra
    - sunspot structure

Parameters of global star:
    - i_sun: inclination angle, i=pi/2 is edge-on, i=0 is pole on.
    - l0: 'epoch angle', if L0=0, spot with coordinate (0,0)
    starts in center of disk, if L0=0, spot is in center of
    disk after half a rotation period
    - period: equatorial period
    - b: differential rotation b = Omega_pole - Omega_eq where Omega is in radians per time unit (angular velocity)
    - ap,bp,cp: limb darkening

Parameters of spots:
    - lambda: 'greenwich' coordinate (longitude) (pi means push to back of star if l0=0)
    - theta: latitude (latitude=0 means at equator, latitude=pi/2 means on pole (B{not
    colatitude})
    - A: size
    - maxt0 (optional): time of maximum size
    - decay (optional): exponential decay speed
    

Note: alpha = (Omega_eq - Omega_p) / Omega_eq
      alpha = -b /2*pi * period

Use for fitting!

gauss(x, parameters)

Evaluate a gaussian fit.

Parameter fields: 'ampl','mu','sigma' and optionally 'const'.

Evaluating one Gaussian:

>>> x = np.linspace(-20,20,10000)
>>> pars = [0.5,0.,3.]
>>> y = gauss(x,pars)
>>> p = pl.plot(x,y,'k-')

Evaluating multiple Gaussian, with and without constants:

>>> pars = [0.5,0.,3.,0.1,-10,1.,.1]
>>> y = gauss(x,pars)
>>> p = pl.plot(x,y,'r-')

Parameters:

x (ndarray) - domain
parameters - record array.

Returns: array

fitted signal

Decorators:

@check_input

lorentz(x, parameters)

Evaluate Lorentz profile

P(f) = A / ( (x-mu)**2 + gamma**2)

Parameters fields: ampl,mu,gamma and optionally const.

Evaluating one Lorentzian:

>>> x = np.linspace(-20,20,10000)
>>> pars = [0.5,0.,3.]
>>> y = lorentz(x,pars)
>>> p = pl.figure()
>>> p = pl.plot(x,y,'k-')

Evaluating multiple Lorentzians, with and without constants:

>>> pars = [0.5,0.,3.,0.1,-10,1.,.1]
>>> y = lorentz(x,pars)
>>> p = pl.plot(x,y,'r-')

Decorators:

@check_input

voigt(x, parameters)

Evaluate a Voigt profile.

Field parameters: ampl, mu, gamma, sigma and optionally const

Function:

   z = (x + gamma*i) / (sigma*sqrt(2))
   V = A * Real[cerf(z)] / (sigma*sqrt(2*pi))

>>> x = np.linspace(-20,20,10000)
>>> pars1 = [0.5,0.,2.,2.]
>>> pars2 = [0.5,0.,1.,2.]
>>> pars3 = [0.5,0.,2.,1.]
>>> y1 = voigt(x,pars1)
>>> y2 = voigt(x,pars2)
>>> y3 = voigt(x,pars3)
>>> p = pl.figure()
>>> p = pl.plot(x,y1,'b-')
>>> p = pl.plot(x,y2,'g-')
>>> p = pl.plot(x,y3,'r-')

Multiple voigt profiles:

>>> pars4 = [0.5,0.,2.,1.,0.1,-3,0.5,0.5,0.01]
>>> y4 = voigt(x,pars4)
>>> p = pl.plot(x,y4,'c-')

Returns: array

Voigt profile

Decorators:

@check_input

power_law(x, parameters)

Evaluate a power law.

Parameter fields: A, B, C, optionally const

Function:

P(f) = A / (1+ Bf)**C + const

>>> x = np.linspace(0,10,1000)
>>> pars1 = [0.5,1.,1.]
>>> pars2 = [0.5,1.,2.]
>>> pars3 = [0.5,2.,1.]
>>> y1 = power_law(x,pars1)
>>> y2 = power_law(x,pars2)
>>> y3 = power_law(x,pars3)
>>> p = pl.figure()
>>> p = pl.plot(x,y1,'b-')
>>> p = pl.plot(x,y2,'g-')
>>> p = pl.plot(x,y3,'r-')

Multiple power laws: >>> pars4 = pars1+pars2+pars3 >>> y4 = power_law(x,pars4) >>> p = pl.plot(x,y4,'c-')

Parameters:

x (ndarray) - domain
parameters - record array.

Returns: array

fitted signal

Decorators:

@check_input

sine_preppars(pars)

Prepare sine function parameters in correct form for evaluating/fitting.

If you input a record array, this function will output a 1D numpy array containing only the independent parameters for use in nonlinear fitting algorithms.

If you input a 1D numpy array, it will output a record array.

Parameters:

pars (record array or normal numpy array) - input parameters in record array or normal array form

Returns: normal numpy array or record array

input parameters in normal array form or record array

kepler_preppars(pars)

Prepare Kepler orbit parameters in correct form for evaluating/fitting.

If you input a record array, this function will output a 1D numpy array containing only the independent parameters for use in nonlinear fitting algorithms.

If you input a 1D numpy array, it will output a record array.

Parameters:

pars (record array or normal numpy array) - input parameters in record array or normal array form

Returns: normal numpy array or record array

input parameters in normal array form or record array

binary_preppars(pars)

Prepare Kepler orbit parameters in correct form for evaluating/fitting.

If you input a record array, this function will output a 1D numpy array containing only the independent parameters for use in nonlinear fitting algorithms.

If you input a 1D numpy array, it will output a record array.

Parameters:

pars (record array or normal numpy array) - input parameters in record array or normal array form

Returns: normal numpy array or record array

input parameters in normal array form or record array

box_preppars(pars)

Prepare sine function parameters in correct form for evaluating/fitting.

If you input a record array, this function will output a 1D numpy array containing only the independent parameters for use in nonlinear fitting algorithms.

If you input a 1D numpy array, it will output a record array.

Parameters:

pars (record array or normal numpy array) - input parameters in record array or normal array form

Returns: normal numpy array or record array

input parameters in normal array form or record array

gauss_preppars(pars)

Prepare gauss function parameters in correct form for evaluating/fitting.

If you input a record array, this function will output a 1D numpy array containing only the independent parameters for use in nonlinear fitting algorithms.

If you input a 1D numpy array, it will output a record array.

Parameters:

pars (record array or normal numpy array) - input parameters in record array or normal array form

Returns: normal numpy array or record array

input parameters in normal array form or record array

lorentz_preppars(pars)

Prepare Lorentz function parameters in correct form for evaluating/fitting.

If you input a record array, this function will output a 1D numpy array containing only the independent parameters for use in nonlinear fitting algorithms.

If you input a 1D numpy array, it will output a record array.

Parameters:

pars (record array or normal numpy array) - input parameters in record array or normal array form

Returns: normal numpy array or record array

input parameters in normal array form or record array

voigt_preppars(pars)

Prepare voigt function parameters in correct form for evaluating/fitting.

If you input a record array, this function will output a 1D numpy array containing only the independent parameters for use in nonlinear fitting algorithms.

If you input a 1D numpy array, it will output a record array.

Parameters:

pars (record array or normal numpy array) - input parameters in record array or normal array form

Returns: normal numpy array or record array

input parameters in normal array form or record array

power_law_preppars(pars)

Prepare gauss function parameters in correct form for evaluating/fitting.

If you input a record array, this function will output a 1D numpy array containing only the independent parameters for use in nonlinear fitting algorithms.

If you input a 1D numpy array, it will output a record array.

Parameters:

pars (record array or normal numpy array) - input parameters in record array or normal array form

Returns: normal numpy array or record array

input parameters in normal array form or record array