The following are code examples for showing how to use . They are extracted from open source Python projects. You can vote up the examples you like or vote down the exmaples you don’t like. You can also save this page to your account.
Example 1
def fit_gamma(ti):
mu = ti.mean()
sigma = ti.std()
C = mu**2/sigma**2
def f(x):
return np.exp(x)/(2.*np.expm1(x)/x - 1.)
def f1(x):
y = 2.*np.expm1(x)/x - 1.
z = 2./x**2*(x*np.exp(x) - np.expm1(x))
return np.exp(x)*(1 - z/y)/y
# newton solve
K = 2.*C # initial value
#print "Newton iteration:"
for i in range(10):
dK = -(f(K) - C)/f1(K)
K = K + dK
# print i, "Residual", f(K) - C, "Value K =", K
#print
tau = mu*(1. - np.exp(-K))/K
return K, tau Example 2
def sample_scalar(self, shape, a):
AMAX = 30
if a > AMAX:
return np.random.poisson(a, shape)
k = 1
K = np.full(shape, k)
s = a/np.expm1(a)
S = s
U = np.random.random(shape)
new = S < U
while np.any(new):
k += 1
K[new] = k
s = s*a/float(k)
S = S + s
new = S < U
return K Example 3
def c(vec):
"""Complement function for probabilities in the log-space: robustly computes 1-P(A) in the log-space
Args:
vec: vector of negative numbers representing log-probabilities of an event.
Returns: the log-probabilities of (1-P(A)) were log(P(A)) are given in the vec numpy array
Examples:
>>> c(-1e-200)
-460.51701859880916
# >>> np.log(1 - np.exp(-1e-200)) raises a `RuntimeWarning: divide by zero` error
"""
# return np.log1p(-np.exp(vec)) # Not robust to -1e-200
if np.max(np.array(vec)) > 0:
print('vec', vec)
return np.log(-np.expm1(vec)) Example 4
def c(vec):
"""Complement function for probabilities in the log-space: robustly computes 1-P(A) in the log-space
Args:
vec: vector of negative numbers representing log-probabilities of an event.
Returns: the log-probabilities of (1-P(A)) were log(P(A)) are given in the vec numpy array
Examples:
>>> c(-1e-200)
-460.51701859880916
# >>> np.log(1 - np.exp(-1e-200)) raises a `RuntimeWarning: divide by zero` error
"""
# return np.log1p(-np.exp(vec)) # Not robust to -1e-200
if np.max(np.array(vec)) > 0:
print('vec', vec)
return np.log(-np.expm1(vec)) Example 5
def vmf_logp(x, lon_lat, kappa):
if x[1] < -90. or x[1] > 90.:
raise ZeroProbability
return -np.inf
if kappa < eps:
return np.log(1. / 4. / np.pi)
mu = np.array([np.cos(lon_lat[1] * d2r) * np.cos(lon_lat[0] * d2r),
np.cos(lon_lat[1] * d2r) * np.sin(lon_lat[0] * d2r),
np.sin(lon_lat[1] * d2r)])
test_point = np.transpose(np.array([np.cos(x[1] * d2r) * np.cos(x[0] * d2r),
np.cos(x[1] * d2r) *
np.sin(x[0] * d2r),
np.sin(x[1] * d2r)]))
logp_elem = np.log( -kappa / ( 2. * np.pi * np.expm1(-2. * kappa)) ) + \
kappa * (np.dot(test_point, mu) - 1.)
logp = logp_elem.sum()
return logp Example 6
def _logpdf(self, samples):
if self.theta == 0:
vals = np.zeros(samples.shape[0])
else:
vals = np.log(-self.theta * np.expm1(-self.theta)
* np.exp(-self.theta
* (samples[:, 0] + samples[:, 1]))
/ (np.expm1(-self.theta)
+ np.expm1(-self.theta * samples[:, 0])
* np.expm1(-self.theta * samples[:, 1])) ** 2)
return vals Example 7
def _logcdf(self, samples):
if self.theta == 0:
vals = np.sum(np.log(samples), axis=1)
else:
old_settings = np.seterr(divide='ignore')
vals = np.log(-np.log1p(np.expm1(-self.theta * samples[:, 0])
* np.expm1(-self.theta * samples[:, 1])
/ (np.expm1(-self.theta)))) \
- np.log(self.theta)
np.seterr(**old_settings)
return vals Example 8
def _ccdf(self, samples):
if self.theta == 0:
vals = samples[:, 0]
else:
vals = np.exp(-self.theta * samples[:, 1]) \
* np.expm1(-self.theta * samples[:, 0]) \
/ (np.expm1(-self.theta)
+ np.expm1(-self.theta * samples[:, 0])
* np.expm1(-self.theta * samples[:, 1]))
return vals Example 9
def _ppcf(self, samples):
if self.theta == 0:
vals = samples[:, 0]
else:
vals = -np.log1p(samples[:, 0] * np.expm1(-self.theta)
/ (np.exp(-self.theta * samples[:, 1])
- samples[:, 0] * np.expm1(-self.theta
* samples[:, 1]))) \
/ self.theta
return vals Example 10
def prior_and_posterior_for_category(self, category_to_primitives):
category_to_prior = {}
category_to_posterior = {}
# We want to take the logs of the lengh prior, but some elements
# will be zero; entries in the result being -inf is okay, but we want
# to avoid a potentially alarming warning message
M = np.zeros(self.length_prior.shape)
M[self.length_prior > 0.] = np.log(self.length_prior[self.length_prior > 0.])
M[self.length_prior <= 0.] = -np.inf
lengths = np.arange(len(M))
for category, primitives in category_to_primitives.iteritems():
base_prior = 0.
for primitive in primitives:
base_prior += self.base_prior_by_primitive[primitive]
# Effectively we calculate the probability of never
# choosing this category in a rule with 0 primitives, 1
# primitive, etc, then add those up, scaling by the
# probability of having 0 primitives, 1 primitive, etc.;
# for speed and numerical behavior we do this in a
# vectorized way on a log scale:
inclusion_prior = -1.0*np.expm1(
logsumexp(M + (np.log1p(-1.0*base_prior) * lengths))
)
category_to_prior[category] = inclusion_prior
states_in_category = 0
for state in self.state_ensemble_primitives:
if state.intersection(primitives):
states_in_category += 1
category_to_posterior[category] = (
states_in_category / float(len(self.state_ensemble))
)
return category_to_prior, category_to_posterior Example 11
def __init__(self, daily_returns, benchmark_daily_returns, risk_free_rate, days, period=DAILY):
assert(len(daily_returns) == len(benchmark_daily_returns))
self._portfolio = daily_returns
self._benchmark = benchmark_daily_returns
self._risk_free_rate = risk_free_rate
self._annual_factor = _annual_factor(period)
self._alpha = None
self._beta = None
self._sharpe = None
self._return = np.expm1(np.log1p(self._portfolio).sum())
self._annual_return = (1 + self._return) ** (365 / days) - 1
# self._annual_return = (1 + self._return) ** (self._annual_factor / len(self._portfolio)) - 1
self._benchmark_return = np.expm1(np.log1p(self._benchmark).sum())
self._benchmark_annual_return = (1 + self._benchmark_return) ** (365 / days) - 1
# self._benchmark_annual_return = (1 + self._benchmark_return) ** (self._annual_factor / len(self._portfolio)) - 1
self._max_drawdown = None
self._volatility = None
self._annual_volatility = None
self._benchmark_volatility = None
self._benchmark_annual_volatility = None
self._information_ratio = None
self._sortino = None
self._tracking_error = None
self._annual_tracking_error = None
self._downside_risk = None
self._annual_downside_risk = None
self._calmar = None Example 12
def __init__(self, daily_returns, benchmark_daily_returns, risk_free_rate, days, period=DAILY):
assert(len(daily_returns) == len(benchmark_daily_returns))
self._portfolio = daily_returns
self._benchmark = benchmark_daily_returns
self._risk_free_rate = risk_free_rate
self._annual_factor = _annual_factor(period)
self._daily_risk_free_rate = self._risk_free_rate / self._annual_factor
self._alpha = None
self._beta = None
self._sharpe = None
self._return = np.expm1(np.log1p(self._portfolio).sum())
self._annual_return = (1 + self._return) ** (365 / days) - 1
self._benchmark_return = np.expm1(np.log1p(self._benchmark).sum())
self._benchmark_annual_return = (1 + self._benchmark_return) ** (365 / days) - 1
self._max_drawdown = None
self._volatility = None
self._annual_volatility = None
self._benchmark_volatility = None
self._benchmark_annual_volatility = None
self._information_ratio = None
self._sortino = None
self._tracking_error = None
self._annual_tracking_error = None
self._downside_risk = None
self._annual_downside_risk = None
self._calmar = None
self._avg_excess_return = None Example 13
def pdf_direct(self, tt, N=50):
a = self.K
t = tt/self.tau
S = np.ones_like(t)
s = np.ones_like(t)
for k in range(1, N):
s *= (a*t)/(k*(k+1.))
S += s
return np.exp(-t)*a/np.expm1(a) * S /self.tau Example 14
def pdf_bessel(self, tt):
a = self.K
t = tt/self.tau
return np.exp(-t)*np.sqrt(a/t)*iv(1., 2.*np.sqrt(a*t))/self.tau/np.expm1(a) Example 15
def cdf(self, tt, N=50):
a = self.K
tau = self.tau
gamma = sp.stats.gamma.cdf
S = np.zeros_like(tt)
s = 1.
for k in range(1, N):
s *= a/k
S += s*gamma(tt, k, scale=tau)
return 1./np.expm1(a) * S Example 16
def cfd_vec(self, tt, p, N=50):
# cdf that takes parameter as vector input, for fitting
a = p[:, 0:1]
tau = p[:, 1:2]
gamma = sp.stats.gamma.cdf
S = np.ones((p.shape[0], tt.shape[1]))
s = np.ones((1, tt.shape[0]))
for k in range(1, N):
s = s*a/k
S = S + s*gamma(tt, k, scale=tau)
return 1./np.expm1(a) * S Example 17
def f(x):
return np.exp(x)/(2.*np.expm1(x)/x - 1.) Example 18
def f1(b, N=50):
k = np.arange(1, N)
return np.sum(summand(k, b))*1./np.expm1(b) - np.log(b*np.exp(b)/np.expm1(b)) Example 19
def h(x):
return np.sqrt(g(x)*(2. - x/np.expm1(x))) Example 20
def poisson_from_positiveK(mean):
# solve x/(1 - exp(-x)) == mean
def f(x):
return x/(1. - np.exp(-x))
def f1(x):
return (np.expm1(x) - x)/(2.*np.cosh(x) - 2.)
x = solve(mean, f, f1, mean, n=10)
return x Example 21
def f(x):
return exp(x)/(2.*expm1(x)/x - 1.) Example 22
def f(b, N=50):
k = np.arange(1, N)
return np.sum(summand(k, b))*1./np.expm1(b) - np.log(b*np.exp(b)/np.expm1(b)) Example 23
def poisson_from_positiveK(mean):
# solve x/(1 - exp(-x)) == mean
def f(x):
return x/(1. - np.exp(-x))
def f1(x):
return (np.expm1(x) - x)/(2.*np.cosh(x) - 2.)
x = solve_newton(mean, f, f1, mean, n=10)
return x Example 24
def sample_(self, shape, a):
if np.isscalar(a) or a.size == 1:
return self.sample_scalar(shape, a)
#print shape, a.shape
AMAX = 30
k = 1
K = np.full(shape, k)
# for large a values, just use non-truncated poisson
large = broadcast_mask(a > AMAX)
small = broadcast_mask(a <= AMAX)
K[large] = np.random.poisson(a[large], K[large].shape)
Ksmall = K[small]
a = a[small]
# actual algorithm begins here
s = a/np.expm1(a)
S = s
U = np.random.random(Ksmall.shape)
new = S < U
while np.any(new):
k += 1
Ksmall[new] = k
s = s*a/float(k)
S = S + s
new = S < U
K[small] = Ksmall
return K Example 25
def pdf_(self, x, a):
return stats.poisson.pmf(x, a)*np.exp(a)/np.expm1(a) Example 26
def backward(self, y):
"""
Inverse of the softplus transform:
.. math::
x = \log( \exp(y) - 1)
The bound for the input y is [self._lower. inf[, self._lower is
subtracted prior to any calculations. The implementation avoids overflow
explicitly by applying the log sum exp trick:
.. math::
\log ( \exp(y) - \exp(0)) &= ys + \log( \exp(y-ys) - \exp(-ys)) \\
&= ys + \log( 1 - \exp(-ys)
ys = \max(0, y)
As y can not be negative, ys could be replaced with y itself.
However, in case :math:`y=0` this results in np.log(0). Hence the zero is
replaced by a machine epsilon.
.. math::
ys = \max( \epsilon, y)
"""
ys = np.maximum(y - self._lower, np.finfo(settings.float_type).eps)
return ys + np.log(-np.expm1(-ys)) Example 27
def norm_y_inv(y_bc):
return np.expm1((y_bc * norm_y_lambda + 1)**(1/norm_y_lambda)) Example 28
def expm1(self, out=None):
assert out is None
return self.elemwise(np.expm1) Example 29
def step_math(self, dt, J, spiked, voltage, refractory_time):
# reduce all refractory times by dt
refractory_time -= dt
# compute effective dt for each neuron, based on remaining time.
# note that refractory times that have completed midway into this
# timestep will be given a partial timestep, and moreover these will
# be subtracted to zero at the next timestep (or reset by a spike)
delta_t = (dt - refractory_time).clip(0, dt)
# update voltage using discretized lowpass filter
# since v(t) = v(0) + (J - v(0))*(1 - exp(-t/tau)) assuming
# J is constant over the interval [t, t + dt)
voltage -= (J - voltage) * np.expm1(-delta_t / self.tau_rc)
# determine which neurons spiked (set them to 1/dt, else 0)
spiked_mask = voltage > 1
spiked[:] = spiked_mask / dt
# set v(0) = 1 and solve for t to compute the spike time
t_spike = dt + self.tau_rc * np.log1p(
-(voltage[spiked_mask] - 1) / (J[spiked_mask] - 1))
# set spiked voltages to zero, refractory times to tau_ref, and
# rectify negative voltages to a floor of min_voltage
voltage[voltage < self.min_voltage] = self.min_voltage
voltage[spiked_mask] = 0
refractory_time[spiked_mask] = self.tau_ref + t_spike Example 30
def test_numpy_method():
# This type of code is used frequently by PyMC3 users
x = tt.dmatrix('x')
data = np.random.rand(5, 5)
x.tag.test_value = data
for fct in [np.arccos, np.arccosh, np.arcsin, np.arcsinh,
np.arctan, np.arctanh, np.ceil, np.cos, np.cosh, np.deg2rad,
np.exp, np.exp2, np.expm1, np.floor, np.log,
np.log10, np.log1p, np.log2, np.rad2deg,
np.sin, np.sinh, np.sqrt, np.tan, np.tanh, np.trunc]:
y = fct(x)
f = theano.function([x], y)
utt.assert_allclose(np.nan_to_num(f(data)),
np.nan_to_num(fct(data))) Example 31
def impl(self, x):
# If x is an int8 or uint8, numpy.expm1 will compute the result in
# half-precision (float16), where we want float32.
x_dtype = str(getattr(x, 'dtype', ''))
if x_dtype in ('int8', 'uint8'):
return numpy.expm1(x, sig='f')
return numpy.expm1(x) Example 32
def c_code(self, node, name, inputs, outputs, sub):
(x,) = inputs
(z,) = outputs
if node.inputs[0].type in complex_types:
raise NotImplementedError('type not supported', type)
return "%(z)s = expm1(%(x)s);" % locals() Example 33
def __init__(self, daily_returns, benchmark_daily_returns, risk_free_rate, days, period=DAILY):
assert(len(daily_returns) == len(benchmark_daily_returns))
self._portfolio = daily_returns
self._benchmark = benchmark_daily_returns
self._risk_free_rate = risk_free_rate
self._annual_factor = _annual_factor(period)
self._daily_risk_free_rate = self._risk_free_rate / self._annual_factor
self._alpha = None
self._beta = None
self._sharpe = None
self._return = np.expm1(np.log1p(self._portfolio).sum())
self._annual_return = (1 + self._return) ** (365 / days) - 1
self._benchmark_return = np.expm1(np.log1p(self._benchmark).sum())
self._benchmark_annual_return = (1 + self._benchmark_return) ** (365 / days) - 1
self._max_drawdown = None
self._volatility = None
self._annual_volatility = None
self._benchmark_volatility = None
self._benchmark_annual_volatility = None
self._information_ratio = None
self._sortino = None
self._tracking_error = None
self._annual_tracking_error = None
self._downside_risk = None
self._annual_downside_risk = None
self._calmar = None
self._avg_excess_return = None Example 34
def mape_ln(y,d):
c=d.get_label()
result=np.sum(np.abs(np.expm1(y)-np.abs(np.expm1(c)))/np.abs(np.expm1(c)))/len(c)
return "mape",result Example 35
def transform_back(self, transformed_target):
return np.expm1(transformed_target) Example 36
def process(self, **kwargs):
"""Process module."""
Transform.process(self, **kwargs)
self._kappa = kwargs[self.key('kappa')]
self._kappa_gamma = kwargs[self.key('kappagamma')]
self._m_ejecta = kwargs[self.key('mejecta')]
self._v_ejecta = kwargs[self.key('vejecta')]
self._tau_diff = np.sqrt(self.DIFF_CONST * self._kappa *
self._m_ejecta / self._v_ejecta) / DAY_CGS
self._trap_coeff = (
self.TRAP_CONST * self._kappa_gamma * self._m_ejecta /
(self._v_ejecta ** 2)) / DAY_CGS ** 2
td2, A = self._tau_diff ** 2, self._trap_coeff # noqa: F841
new_lums = np.zeros_like(self._times_to_process)
if len(self._dense_times_since_exp) < 2:
return {self.dense_key('luminosities'): new_lums}
min_te = min(self._dense_times_since_exp)
tb = max(0.0, min_te)
linterp = interp1d(
self._dense_times_since_exp, self._dense_luminosities, copy=False,
assume_sorted=True)
uniq_times = np.unique(self._times_to_process[
(self._times_to_process >= tb) & (
self._times_to_process <= self._dense_times_since_exp[-1])])
lu = len(uniq_times)
num = int(round(self.N_INT_TIMES / 2.0))
lsp = np.logspace(
np.log10(self._tau_diff /
self._dense_times_since_exp[-1]) +
self.MIN_LOG_SPACING, 0, num)
xm = np.unique(np.concatenate((lsp, 1 - lsp)))
int_times = np.clip(
tb + (uniq_times.reshape(lu, 1) - tb) * xm, tb,
self._dense_times_since_exp[-1])
int_te2s = int_times[:, -1] ** 2
int_lums = linterp(int_times) # noqa: F841
int_args = int_lums * int_times * np.exp(
(int_times ** 2 - int_te2s.reshape(lu, 1)) / td2)
int_args[np.isnan(int_args)] = 0.0
uniq_lums = np.trapz(int_args, int_times)
uniq_lums *= -2.0 * np.expm1(-A / int_te2s) / td2
new_lums = uniq_lums[np.searchsorted(uniq_times,
self._times_to_process)]
return {self.key('tau_diffusion'): self._tau_diff,
self.dense_key('luminosities'): new_lums} Example 37
def gpinv(p, k, sigma):
"""Inverse Generalised Pareto distribution function."""
x = np.empty(p.shape)
x.fill(np.nan)
if sigma <= 0:
return x
ok = (p > 0) & (p < 1)
if np.all(ok):
if np.abs(k) < np.finfo(float).eps:
np.negative(p, out=x)
np.log1p(x, out=x)
np.negative(x, out=x)
else:
np.negative(p, out=x)
np.log1p(x, out=x)
x *= -k
np.expm1(x, out=x)
x /= k
x *= sigma
else:
if np.abs(k) < np.finfo(float).eps:
# x[ok] = - np.log1p(-p[ok])
temp = p[ok]
np.negative(temp, out=temp)
np.log1p(temp, out=temp)
np.negative(temp, out=temp)
x[ok] = temp
else:
# x[ok] = np.expm1(-k * np.log1p(-p[ok])) / k
temp = p[ok]
np.negative(temp, out=temp)
np.log1p(temp, out=temp)
temp *= -k
np.expm1(temp, out=temp)
temp /= k
x[ok] = temp
x *= sigma
x[p == 0] = 0
if k >= 0:
x[p == 1] = np.inf
else:
x[p == 1] = -sigma / k
return x Example 38
def test_numpy_ufuncs(self):
# test ufuncs of numpy 1.9.2. see:
# http://docs.scipy.org/doc/numpy/reference/ufuncs.html
# some functions are skipped because it may return different result
# for unicode input depending on numpy version
for name, idx in compat.iteritems(self.indices):
for func in [np.exp, np.exp2, np.expm1, np.log, np.log2, np.log10,
np.log1p, np.sqrt, np.sin, np.cos, np.tan, np.arcsin,
np.arccos, np.arctan, np.sinh, np.cosh, np.tanh,
np.arcsinh, np.arccosh, np.arctanh, np.deg2rad,
np.rad2deg]:
if isinstance(idx, pd.tseries.base.DatetimeIndexOpsMixin):
# raise TypeError or ValueError (PeriodIndex)
# PeriodIndex behavior should be changed in future version
with tm.assertRaises(Exception):
func(idx)
elif isinstance(idx, (Float64Index, Int64Index)):
# coerces to float (e.g. np.sin)
result = func(idx)
exp = Index(func(idx.values), name=idx.name)
self.assert_index_equal(result, exp)
self.assertIsInstance(result, pd.Float64Index)
else:
# raise AttributeError or TypeError
if len(idx) == 0:
continue
else:
with tm.assertRaises(Exception):
func(idx)
for func in [np.isfinite, np.isinf, np.isnan, np.signbit]:
if isinstance(idx, pd.tseries.base.DatetimeIndexOpsMixin):
# raise TypeError or ValueError (PeriodIndex)
with tm.assertRaises(Exception):
func(idx)
elif isinstance(idx, (Float64Index, Int64Index)):
# results in bool array
result = func(idx)
exp = func(idx.values)
self.assertIsInstance(result, np.ndarray)
tm.assertNotIsInstance(result, Index)
else:
if len(idx) == 0:
continue
else:
with tm.assertRaises(Exception):
func(idx) Example 39
def nonnegative_bandpass_filter(data,fa=None,fb=None,
Fs=1000.,order=4,zerophase=True,bandstop=False,
offset=1.0):
'''
For filtering data that must remain non-negative. Due to ringing
conventional fitering can create values less than zero for non-
negative real inputs. This may be unrealistic for some data.
To compensate, this performs the filtering on the natural
logarithm of the input data. For small numbers, this can lead to
numeric underflow, so an offset parameter (default 1) is added
to the data for stability.
Parameters
----------
data (ndarray):
data, filtering performed over last dimension
fa (number):
low-freq cutoff Hz. If none, lowpass at fb
fb (number):
high-freq cutoff Hz. If none, highpass at fa
Fs (int):
Sample rate in Hz
order (1..6):
butterworth filter order. Default 4
zerophase (boolean):
Use forward-backward filtering? (true)
bandstop (boolean):
Do band-stop rather than band-pass
offset (positive number):
Offset data to avoid underflow (1)
Returns
-------
filtered :
Filtered signal
'''
offset -= 1.0
data = np.log1p(data+offset)
filtered = bandpass_filter(data,
fa=fa, fb=fb, Fs=Fs,
order=order,
zerophase=zerophase,
bandstop=bandstop)
return np.expm1(filtered) 