def svgd_kernel(self, h = -1):
sq_dist = pdist(self.theta)
pairwise_dists = squareform(sq_dist)**2
if h < 0: # if h < 0, using median trick
h = np.median(pairwise_dists)
h = np.sqrt(0.5 * h / np.log(self.theta.shape[0]+1))
# compute the rbf kernel
Kxy = np.exp( -pairwise_dists / h**2 / 2)
dxkxy = -np.matmul(Kxy, self.theta)
sumkxy = np.sum(Kxy, axis=1)
for i in range(self.theta.shape[1]):
dxkxy[:, i] = dxkxy[:,i] + np.multiply(self.theta[:,i],sumkxy)
dxkxy = dxkxy / (h**2)
return (Kxy, dxkxy)
python类log()的实例源码
def saveHintonPlot(self, matrix, num_tests, max_weight=None, ax=None):
"""Draw Hinton diagram for visualizing a weight matrix."""
fig,ax = plt.subplots(1,1)
if not max_weight:
max_weight = 2**np.ceil(np.log(np.abs(matrix).max())/np.log(2))
ax.patch.set_facecolor('gray')
ax.set_aspect('equal', 'box')
ax.xaxis.set_major_locator(plt.NullLocator())
ax.yaxis.set_major_locator(plt.NullLocator())
for (x, y), w in np.ndenumerate(matrix):
color = 'white' if w > 0 else 'black'
size = np.sqrt(np.abs(0.5*w/num_tests)) # Need to scale so that it is between 0 and 0.5
rect = plt.Rectangle([x - size / 2, y - size / 2], size, size,
facecolor=color, edgecolor=color)
ax.add_patch(rect)
ax.autoscale_view()
ax.invert_yaxis()
plt.savefig(self.figures_path + self.save_prefix + '-Hinton.eps')
plt.close()
def _factor_target_indices(self, Y_inds, vocab_size=None, base=2):
if vocab_size is None:
vocab_size = len(self.dp.word_index)
print >>sys.stderr, "Factoring targets of vocabulary size: %d"%(vocab_size)
num_vecs = int(math.ceil(math.log(vocab_size)/math.log(base))) + 1
base_inds = []
div_Y_inds = Y_inds
print >>sys.stderr, "Number of factors: %d"%num_vecs
for i in range(num_vecs):
new_inds = div_Y_inds % base
if i == num_vecs - 1:
if new_inds.sum() == 0:
# Most significant "digit" is a zero. Omit it.
break
base_inds.append(new_inds)
div_Y_inds = numpy.copy(div_Y_inds/base)
base_vecs = [self._make_one_hot(base_inds_i, base) for base_inds_i in base_inds]
return base_vecs
def sampleSize(self):
r = 0.3
alpha = 0.05
# power=0.9
C = 0.5 * np.log((1 + r) / (1 - r))
Za = scipy.stats.norm.ppf(1 - (0.05 / 2))
sizeArray = []
powerArray = []
power = 0.5
for i in range(50, 100, 1):
power = i / 100
powerArray.append(power)
Zb = scipy.stats.norm.ppf(1 - power)
N = abs((Za - Zb) / C)**2 + 3
sizeArray.append(N)
return [powerArray, sizeArray]
def sample(self, sample_size=20, text=None):
"""Sample the documents."""
p = 1
if text != None:
try:
x, word_idxs = self.reader.get(text)
except Exception as e:
print(e)
return
else:
x, word_idxs = self.reader.random()
print(" [*] Text: %s" % " ".join([self.reader.idx2word[word_idx] for word_idx in word_idxs]))
cur_ps = self.sess.run(self.p_x_i, feed_dict={self.x: x})
word_idxs = np.array(cur_ps).argsort()[-sample_size:][::-1]
ps = cur_ps[word_idxs]
for idx, (cur_p, word_idx) in enumerate(zip(ps, word_idxs)):
print(" [%d] %-20s: %.8f" % (idx+1, self.reader.idx2word[word_idx], cur_p))
p *= cur_p
print(" [*] perp : %8.f" % -np.log(p))
def lr_grad_norm_avg(self):
# this is for enforcing lr * grad_norm not
# increasing dramatically in case of instability.
# Not necessary for basic use.
global_state = self._global_state
beta = self._beta
if "lr_grad_norm_avg" not in global_state:
global_state['grad_norm_squared_avg_log'] = 0.0
global_state['grad_norm_squared_avg_log'] = \
global_state['grad_norm_squared_avg_log'] * beta \
+ (1 - beta) * np.log(global_state['grad_norm_squared'] + eps)
if "lr_grad_norm_avg" not in global_state:
global_state["lr_grad_norm_avg"] = \
0.0 * beta + (1 - beta) * np.log(self._lr * np.sqrt(global_state['grad_norm_squared'] ) + eps)
# we monitor the minimal smoothed ||lr * grad||
global_state["lr_grad_norm_avg_min"] = \
np.exp(global_state["lr_grad_norm_avg"] / self.zero_debias_factor() )
else:
global_state["lr_grad_norm_avg"] = global_state["lr_grad_norm_avg"] * beta \
+ (1 - beta) * np.log(self._lr * np.sqrt(global_state['grad_norm_squared'] ) + eps)
global_state["lr_grad_norm_avg_min"] = \
min(global_state["lr_grad_norm_avg_min"],
np.exp(global_state["lr_grad_norm_avg"] / self.zero_debias_factor() ) )
def lr_grad_norm_avg(self):
# this is for enforcing lr * grad_norm not
# increasing dramatically in case of instability.
# Not necessary for basic use.
global_state = self._global_state
beta = self._beta
if "lr_grad_norm_avg" not in global_state:
global_state['grad_norm_squared_avg_log'] = 0.0
global_state['grad_norm_squared_avg_log'] = \
global_state['grad_norm_squared_avg_log'] * beta \
+ (1 - beta) * np.log(global_state['grad_norm_squared'] + eps)
if "lr_grad_norm_avg" not in global_state:
global_state["lr_grad_norm_avg"] = \
0.0 * beta + (1 - beta) * np.log(self._lr * np.sqrt(global_state['grad_norm_squared'] ) + eps)
# we monitor the minimal smoothed ||lr * grad||
global_state["lr_grad_norm_avg_min"] = \
np.exp(global_state["lr_grad_norm_avg"] / self.zero_debias_factor() )
else:
global_state["lr_grad_norm_avg"] = global_state["lr_grad_norm_avg"] * beta \
+ (1 - beta) * np.log(self._lr * np.sqrt(global_state['grad_norm_squared'] ) + eps)
global_state["lr_grad_norm_avg_min"] = \
min(global_state["lr_grad_norm_avg_min"],
np.exp(global_state["lr_grad_norm_avg"] / self.zero_debias_factor() ) )
def sample_dtype(self):
return tf.int32
# WRONG SECOND DERIVATIVES
# class CategoricalPd(Pd):
# def __init__(self, logits):
# self.logits = logits
# self.ps = tf.nn.softmax(logits)
# @classmethod
# def fromflat(cls, flat):
# return cls(flat)
# def flatparam(self):
# return self.logits
# def mode(self):
# return U.argmax(self.logits, axis=1)
# def logp(self, x):
# return -tf.nn.sparse_softmax_cross_entropy_with_logits(self.logits, x)
# def kl(self, other):
# return tf.nn.softmax_cross_entropy_with_logits(other.logits, self.ps) \
# - tf.nn.softmax_cross_entropy_with_logits(self.logits, self.ps)
# def entropy(self):
# return tf.nn.softmax_cross_entropy_with_logits(self.logits, self.ps)
# def sample(self):
# u = tf.random_uniform(tf.shape(self.logits))
# return U.argmax(self.logits - tf.log(-tf.log(u)), axis=1)
def spec_entropy(Rates,time_range=[],bin_w = 5.,freq_range = []):
'''Function to calculate the spectral entropy'''
power,freq,dfreq,dummy,dummy = mypsd(Rates,time_range,bin_w = bin_w)
if freq_range != []:
power = power[(freq>=freq_range[0]) & (freq <= freq_range[1])]
freq = freq[(freq>=freq_range[0]) & (freq <= freq_range[1])]
maxFreq = freq[np.where(power==np.max(power))]*1000*100
perMax = (np.max(power)/np.sum(power))*100
k = len(freq)
power = power/sum(power)
sum_power = 0
for ii in range(k):
sum_power += (power[ii]*np.log(power[ii]))
spec_ent = -(sum_power/np.log(k))
return spec_ent,dfreq,maxFreq,perMax
def __init__(self, nlp, **kwargs):
"""
Initializes a new instance of SpacySimilarityHook class.
:param nlp: `spaCy language object <https://spacy.io/docs/api/language>`_.
:param ignore_stops: Indicates whether to ignore the stop words.
:param only_alpha: Indicates whether only alpha tokens must be used.
:param frequency_processor: The function which is applied to raw \
token frequencies.
:type ignore_stops: bool
:type only_alpha: bool
:type frequency_processor: callable
"""
self.nlp = nlp
self.ignore_stops = kwargs.get("ignore_stops", True)
self.only_alpha = kwargs.get("only_alpha", True)
self.frequency_processor = kwargs.get(
"frequency_processor", lambda t, f: numpy.log(1 + f))
def gaussian_nll(x, mus, sigmas):
"""
NLL for Multivariate Normal with diagonal covariance matrix
See:
wikipedia.org/wiki/Multivariate_normal_distribution#Likelihood_function
where \Sigma = diag(s_1^2,..., s_n^2).
x, mus, sigmas all should have the same shape.
sigmas (s_1,..., s_n) should be strictly positive.
Results in output shape of similar but without the last dimension.
"""
nll = lib.floatX(numpy.log(2. * numpy.pi))
nll += 2. * T.log(sigmas)
nll += ((x - mus) / sigmas) ** 2.
nll = nll.sum(axis=-1)
nll *= lib.floatX(0.5)
return nll
def monotoneTFosc(f):
"""Maps [-inf,inf] to [-inf,inf] with different constants
for positive and negative part.
"""
if np.isscalar(f):
if f > 0.:
f = np.log(f) / 0.1
f = np.exp(f + 0.49 * (np.sin(f) + np.sin(0.79 * f))) ** 0.1
elif f < 0.:
f = np.log(-f) / 0.1
f = -np.exp(f + 0.49 * (np.sin(0.55 * f) + np.sin(0.31 * f))) ** 0.1
return f
else:
f = np.asarray(f)
g = f.copy()
idx = (f > 0)
g[idx] = np.log(f[idx]) / 0.1
g[idx] = np.exp(g[idx] + 0.49 * (np.sin(g[idx]) + np.sin(0.79 * g[idx])))**0.1
idx = (f < 0)
g[idx] = np.log(-f[idx]) / 0.1
g[idx] = -np.exp(g[idx] + 0.49 * (np.sin(0.55 * g[idx]) + np.sin(0.31 * g[idx])))**0.1
return g
def logscale_img(img_array,
cap=255.0,
coeff=1000.0):
'''
This scales the image according to the relation:
logscale_img = np.log(coeff*(img/max(img))+1)/np.log(coeff)
Taken from the DS9 scaling algorithms page at:
http://hea-www.harvard.edu/RD/ds9/ref/how.html
According to that page:
coeff = 1000.0 works well for optical images
coeff = 100.0 works well for IR images
'''
logscaled_img = np.log(coeff*img_array/np.nanmax(img_array)+1)/np.log(coeff)
return cap*logscaled_img
def forward(self, outputs, targets):
"""SoftmaxCategoricalCrossEntropy forward propagation.
.. math:: L_i = - \\sum_j{t_{i,j} \\log(p_{i,j})}
Parameters
----------
outputs : numpy.array
Predictions in (0, 1), such as softmax output of a neural network,
with data points in rows and class probabilities in columns.
targets : numpy.array
Either targets in [0, 1] matching the layout of `outputs`, or
a vector of int giving the correct class index per data point.
Returns
-------
numpy 1D array
An expression for the item-wise categorical cross-entropy.
"""
outputs = np.clip(outputs, self.epsilon, 1 - self.epsilon)
return np.mean(-np.sum(targets * np.log(outputs), axis=1))
def test_GWD(self):
# Compute categorical crossentropy
indices = self.mock_y > 0
selected_log = -np.log(self.mock_x_softmax[indices])
self.loss = 0#np.sum(selected_log) / np.sum(self.mock_y)
# Create keras model with this activation and compile it
model = Sequential()
activation_layer = Lambda(lambda x: x,
input_shape=self.data_shape[1:],
output_shape=self.data_shape[1:]
)
model.add(activation_layer)
model.compile('sgd', loss=gwd)
# Predict data from the model
loss = model.evaluate(self.mock_y, self.mock_y, batch_size=1, verbose=0)
# Assertions
print('Expected loss: {}'.format(self.loss))
print('Actual loss: {}'.format(loss))
self.assertTrue(np.allclose(loss, self.loss),
msg='Categorical cross-entropy loss 3D does not produce the expected results')
def logTickValues(self, minVal, maxVal, size, stdTicks):
## start with the tick spacing given by tickValues().
## Any level whose spacing is < 1 needs to be converted to log scale
ticks = []
for (spacing, t) in stdTicks:
if spacing >= 1.0:
ticks.append((spacing, t))
if len(ticks) < 3:
v1 = int(np.floor(minVal))
v2 = int(np.ceil(maxVal))
#major = list(range(v1+1, v2))
minor = []
for v in range(v1, v2):
minor.extend(v + np.log10(np.arange(1, 10)))
minor = [x for x in minor if x>minVal and x<maxVal]
ticks.append((None, minor))
return ticks
def logTickValues(self, minVal, maxVal, size, stdTicks):
## start with the tick spacing given by tickValues().
## Any level whose spacing is < 1 needs to be converted to log scale
ticks = []
for (spacing, t) in stdTicks:
if spacing >= 1.0:
ticks.append((spacing, t))
if len(ticks) < 3:
v1 = int(np.floor(minVal))
v2 = int(np.ceil(maxVal))
#major = list(range(v1+1, v2))
minor = []
for v in range(v1, v2):
minor.extend(v + np.log10(np.arange(1, 10)))
minor = [x for x in minor if x>minVal and x<maxVal]
ticks.append((None, minor))
return ticks
def summary(self):
s = OrderedDict()
s['source kind'] = self.sourcekind
s['source'] = self.source
if self.params is not None:
for param, value in self.params._asdict().items():
s['parameter: %s' % param] = value
s['log-variance theoretical half-life'] = self.params.logvarhalflife()
s['log-variance theoretical unconditional s.d.'] = np.sqrt(self.params.logvaruncondvar())
s['log-return sample mean'] = np.mean(self.svdf['logreturn'])
s['log-return sample s.d.'] = np.sqrt(np.var(self.svdf['logreturn']))
if 'logvar' in self.svdf.columns:
s['log-variance sample mean'] = np.mean(self.svdf['logvar'])
s['log-variance sample s.d.'] = np.sqrt(np.var(self.svdf['logvar']))
s['correlation timing'] = self.cortiming
s['log-return forward?'] = self.logreturnforward
s['log-return scale'] = self.logreturnscale
return s
def multiclass_log_loss(y_true, y_pred, eps=1e-15):
"""Multi class version of Logarithmic Loss metric.
https://www.kaggle.com/wiki/MultiClassLogLoss
Parameters
----------
y_true : array, shape = [n_samples]
true class, intergers in [0, n_classes - 1)
y_pred : array, shape = [n_samples, n_classes]
Returns
-------
loss : float
"""
predictions = np.clip(y_pred, eps, 1 - eps)
# normalize row sums to 1
predictions /= predictions.sum(axis=1)[:, np.newaxis]
actual = np.zeros(y_pred.shape)
n_samples = actual.shape[0]
actual[np.arange(n_samples), y_true.astype(int)] = 1
vectsum = np.sum(actual * np.log(predictions))
loss = -1.0 / n_samples * vectsum
return loss
def step(self, mode):
if mode == "train" and self.mode == "test":
raise Exception("Cannot train during test mode")
if mode == "train":
theano_fn = self.train_fn
batch_gen = self.train_batch_gen
elif mode == "test":
theano_fn = self.test_fn
batch_gen = self.test_batch_gen
else:
raise Exception("Invalid mode")
data = next(batch_gen)
ys = data[-1]
data = data[:-1]
ret = theano_fn(*data)
return {"prediction": np.exp(ret[0]) - 1,
"answers": ys,
"current_loss": ret[1],
"loss_reg": ret[2],
"loss_mse": ret[1] - ret[2],
"log": ""}
def evaluation(self, X_test, y_test):
# normalization
X_test = self.normalization(X_test)
# average over the output
pred_y_test = np.zeros([self.M, len(y_test)])
prob = np.zeros([self.M, len(y_test)])
'''
Since we have M particles, we use a Bayesian view to calculate rmse and log-likelihood
'''
for i in range(self.M):
w1, b1, w2, b2, loggamma, loglambda = self.unpack_weights(self.theta[i, :])
pred_y_test[i, :] = self.nn_predict(X_test, w1, b1, w2, b2) * self.std_y_train + self.mean_y_train
prob[i, :] = np.sqrt(np.exp(loggamma)) /np.sqrt(2*np.pi) * np.exp( -1 * (np.power(pred_y_test[i, :] - y_test, 2) / 2) * np.exp(loggamma) )
pred = np.mean(pred_y_test, axis=0)
# evaluation
svgd_rmse = np.sqrt(np.mean((pred - y_test)**2))
svgd_ll = np.mean(np.log(np.mean(prob, axis = 0)))
return (svgd_rmse, svgd_ll)
def sample(self, sess, chars, vocab, num, prime, temperature):
state = self.cell.zero_state(1, tf.float32).eval()
for char in prime[:-1]:
x = np.zeros((1, 1))
x[0, 0] = vocab[char]
feed = {self.input_data: x, self.initial_state: state}
[state] = sess.run([self.final_state], feed)
def weighted_pick(a):
a = a.astype(np.float64)
a = a.clip(min=1e-20)
a = np.log(a) / temperature
a = np.exp(a) / (np.sum(np.exp(a)))
return np.argmax(np.random.multinomial(1, a, 1))
char = prime[-1]
for n in range(num):
x = np.zeros((1, 1))
x[0, 0] = vocab[char]
feed = {self.input_data: x, self.initial_state: state}
[probs, state] = sess.run([self.probs, self.final_state], feed)
p = probs[0]
sample = weighted_pick(p)
char = chars[sample]
yield char
def update_qz(self, left = None, right = None):
if left == None:
left = 0
right = self.lc.n
for index in range(left, right):
#if index % 100 == 0:
# print index
# sys.stdout.flush()
qz_1 = np.log(self.theta)
qz_0 = np.log(1 - self.theta)
for (label, worker) in zip(*self.lc.crowd_labels[index]):
if label >=0 and worker in self.quv:
qz_1 += expectation_z(self.quv[worker][0], self.quv[worker][1], label)
qz_0 += expectation_z(self.quv[worker][2], self.quv[worker][3], label)
qz_1 = np.exp(qz_1)
qz_0 = np.exp(qz_0)
temp = qz_1 * 1.0 / (qz_0 + qz_1)
if not math.isnan(temp):
self.total_changes += np.abs(self.qz[index] - temp)
self.qz[index] = temp
def fit_normal(self, a):
"""
fit a normal distribution to [(x, p)]
where p is in log scale
"""
# normalize p:
#print a
list_p = []
for (x, p) in a: list_p.append(p)
ps = scipy.misc.logsumexp(list_p)
s = 0
ss = 0
for (x, p) in a:
s += x * np.exp(p - ps)
ss += x*x * np.exp(p - ps)
var = ss - s*s
ep = 1E-300
if var < ep: var = ep
return (s, var)
def eval_pdf(self, worker, var, x, list_il):
"""
Eval *LOG* of pdf of new qu or qv
"""
if var == 'sen':
res = expectation_binorm('v', self.quv[worker][2], self.quv[worker][3], \
x, self.get_mu(worker), self.get_c(worker))
#print res
for (item, label) in list_il:
item_w = self.get_item_w(item)
res += item_w * self.qz[item] * np.log( self.Ber(S(x), label))
else:
res = expectation_binorm('u', self.quv[worker][0], self.quv[worker][1], \
x, self.get_mu(worker), self.get_c(worker))
for (item, label) in list_il:
item_w = self.get_item_w(item)
res += item_w * (1 - self.qz[item]) * np.log( self.Ber(S(x), label))
return res
def expectation_z(mu, var, L, w = 3):
"""
USE LAPLACE approx
Evaluate the expectation of log[ S(u)^L (1-S(u))^(1-L) ]
when u ~ Normal(mu, var)
"""
if L == 1:
f = lambda u: np.log(S(u))
fpp = lambda x: -np.exp(x) / pow(1 + np.exp(x), 2)
return f(mu) + 0.5* fpp(mu)*var
else:
f = lambda u: np.log(1 - S(u))
fpp = lambda x: -np.exp(x) / pow(1 + np.exp(x), 2)
return f(mu) + 0.5* fpp(mu)*var
# need: E[u~N](f)
return scipy.integrate.quad(f, mu - w*std, mu + w*std)[0]
def expectation_z_quad(mu, var, L, w = 3):
"""
USE QUAD (NUMERICAL INTEGRATION)
Evaluate the expectation of log[ S(u)^L (1-S(u))^(1-L) ]
when u ~ Normal(mu, var)
"""
#U = np.random.normal(mu, np.sqrt(var), num)
if L == 1:
f = lambda u: scipy.stats.norm.pdf(u, loc = mu, scale = np.sqrt(var) ) * np.log(S(u))
else:
f = lambda u: scipy.stats.norm.pdf(u, loc = mu, scale = np.sqrt(var) ) * (np.log(1 - S(u)))
#return f
std = np.sqrt(var)
#return scipy.integrate.quad(f, mu - w*std, mu + w*std, epsabs = 1.0e-4, epsrel = 1.0e-4, limit = 25)[0]
return scipy.integrate.quad(f, mu - w*std, mu + w*std)[0]
def expectation_binorm(rv, mu, var, x, M, C, w = 3):
"""
Evaluate the expectation of log Norm(uv| M, C)
x = u, v ~ Norm(mu, var) rv == 'v'
x = v, u ~ Norm(mu, var) rv == 'u'
"""
#print rv, mu, var, x, M, C
if rv == 'v':
f = lambda v: scipy.stats.norm.pdf(v, loc = mu, scale = np.sqrt(var) ) * \
np.log(scipy.stats.multivariate_normal.pdf([x, v], mean = M, cov = C, allow_singular = True))
else:
f = lambda u: scipy.stats.norm.pdf(u, loc = mu, scale = np.sqrt(var) ) * \
np.log(scipy.stats.multivariate_normal.pdf([u, x], mean = M, cov = C, allow_singular = True))
#return f
#print f(mu)
std = np.sqrt(var)
return scipy.integrate.quad(f, mu - w*std, mu + w*std)[0]
def austerity(log_lik,log_density_prior, X,epsilon,batch_size=30,chain_size=10000, thinning=1, theta_t=np.random.randn()):
A = [theta_t]
N = X.shape[0]
dimension=1
if hasattr(theta_t, "__len__"):
dimension = len(theta_t)
global_evals = 0
for i in range(chain_size*thinning-1):
# if i % 1000 ==0:
# print( 100.0*i/(chain_size*thinning), ' %')
theta_prime = np.random.randn(dimension)+theta_t
u = np.random.rand()
mu_0 = np.log(u)+log_density_prior(theta_t) -log_density_prior(theta_prime)
mu_0 = mu_0/N
accept,evals = approximate_MH_accept(mu_0, log_lik, X, batch_size, epsilon, theta_prime, theta_t, N)
global_evals += evals
if accept:
theta_t = theta_prime
A.append(theta_t)
return np.array(A[::thinning]),global_evals
def evaluate_stance(
self,
testN,
testtimes,
testinfecting_vec,
testinfected_vec,
testeventmemes,
testW,
testT,
):
predictednode_vec = []
for next_event_index in xrange(len(testtimes)):
print 'considering event', next_event_index
words = testW[next_event_index, :]
predictions = []
for label in set(self.node_vec):
loglikelihood_term = 0
loglikelihood_term += np.dot(words,
np.log(self.beta[label, :]))
predictions.append((label, loglikelihood_term))
predictednode_vec.append(max(predictions, key=lambda x: \
x[1])[0])
return predictednode_vec
