code: 
  
  
  import superimport
  
  import numpy as np
  import matplotlib.pyplot as plt
  import pyprobml_utils as pml
  
  import tensorflow as tf
  import tensorflow_datasets as tfds
  
  np.random.seed(0)
  
  ds, info = tfds.load('emnist', split='test', shuffle_files=False, with_info=True) # horribly slow
  print(info)
  
  
  plt.figure(figsize=(10, 10))
  i = 0
  for example in ds:
      image = example["image"]
      label = example["label"]
      plt.subplot(5, 5, i+1)
      plt.xticks([])
      plt.yticks([])
      plt.grid(False)
      plt.imshow(image)
      plt.title(label)
      i += 1
      if i >= 25: break
  
  pml.savefig("emnist-data.pdf")
  plt.show() import superimport
  
  import itertools
  import matplotlib.pyplot as plt
  import numpy as np
  from scipy.linalg import eigh
  from sklearn.cluster import KMeans
  from sklearn.metrics.pairwise import rbf_kernel
  import pyprobml_utils as pml
  
  plt.style.use('classic')
  
  def spectral_clustering_demo():
      np.random.seed(0)
      num_clusters = 2
      for data_type, data in (('circle', sample_circle(num_clusters)),
                              ('spiral', sample_spiral())):
          kmeans = KMeans(n_clusters=num_clusters, random_state=0)
          kmeans.fit(data)
          assignments = kmeans.predict(data)
          plot_data(data, assignments, 'k-means clustering', data_type)
  
          sigma = 0.1
          gamma = 1 / (2 * sigma ** 2)
          W = rbf_kernel(data, gamma=gamma)
          d = np.sum(W, 1, keepdims=True)
          sqrt_d = np.sqrt(d)
  
          normalized_W = (W / sqrt_d) / sqrt_d.T
          paranoid_assert(W, normalized_W, False)
  
          # We select the largest eigen values of normalized_W, rather
          # than the smallest eigenvalues of I - normalized_W.  The two
          # problems are equivalent. The eigen values can be converted
          # between the two problems via `1 - eigen_values`. The eigen
          # vectors are the same between both problems.
          eigen_values, eigen_vectors = eigh(normalized_W,
                                             # Get only the top num_clusters eigenvalues
                                             eigvals=(data.shape[0] - num_clusters, data.shape[0]-1))
          eigen_vectors = eigen_vectors / np.linalg.norm(eigen_vectors, axis=1, keepdims=True)
  
          kmeans.fit(eigen_vectors)
          assignments = kmeans.predict(eigen_vectors)
          plot_data(data, assignments, 'spectral clustering', data_type)
  
          plt.show()
  
  def paranoid_assert(W, normalized_W, enable):
      if not enable:
          return
      D = np.diag(np.sum(W, 1))
      L = D - W
      D_inv_sqrt = np.diag(1 / np.diag(np.sqrt(D)))
      np.testing.assert_almost_equal(np.sum(L, 1), 0, err_msg="Rows of Laplacian must sum to 0.")
      np.testing.assert_allclose(normalized_W, D_inv_sqrt * W * D_inv_sqrt, rtol=0, atol=1)
  
  def sample_circle(num_clusters):
      points_per_cluster = 500
      bandwidth = 0.1
  
      data = np.zeros((num_clusters * points_per_cluster, 2))
      for k, n in itertools.product(range(num_clusters), range(points_per_cluster)):
          theta = 2 * np.pi * np.random.uniform()
          rho = k + 1 + np.random.randn() * bandwidth
          x, y = pol2cart(theta, rho)
          idx = k * points_per_cluster + n
          data[idx, 0] = x
          data[idx, 1] = y
      data = data.reshape((num_clusters * points_per_cluster, 2))
      return data
  
  def pol2cart(theta, rho):
      x = rho * np.cos(theta)
      y = rho * np.sin(theta)
      return(x, y)
  
  def sample_spiral():
      # Only 2 clusters in this case. This is hard-coded.
      points_per_cluster = 500
      bandwidth = 0.1
  
      data = np.empty((points_per_cluster, 2))
  
      w = np.arange(1, points_per_cluster + 1).astype(np.float32) / points_per_cluster
      data[:,0] = (4 * w + 1) * np.cos(2*np.pi * w) + np.random.randn(points_per_cluster) * bandwidth
      data[:,1] = (4 * w + 1) * np.sin(2*np.pi * w) + np.random.randn(points_per_cluster) * bandwidth
      data = np.vstack((data, -data))
  
      return data
  
  def plot_data(data, assignments, title, data_type):
      fig = plt.figure()
      ax = fig.add_subplot(1,1,1)
      ax.plot(data[assignments == 0, 0], data[assignments == 0, 1], 'o', color='r')
      ax.plot(data[assignments == 1, 0], data[assignments == 1, 1], 'o', color='b')
      ax.set_xlabel('x')
      ax.set_ylabel('y')
      ax.axis('square')
      ax.grid(True)
      ax.set_title(title)
      plt.tight_layout()
      pml.savefig(f"{data_type}_{title.replace(' ', '_')}.pdf")
  
  if __name__ == '__main__':
      spectral_clustering_demo()
   # Bayesian Binary logistic regression in 2d for iris flwoers
  
  # Code is based on 
  # https://github.com/aloctavodia/BAP/blob/master/code/Chp4/04_Generalizing_linear_models.ipynb
  
  import superimport
  
  import pymc3 as pm
  import numpy as np
  import pandas as pd
  import theano.tensor as tt
  #import seaborn as sns
  import scipy.stats as stats
  from scipy.special import expit as logistic
  import matplotlib.pyplot as plt
  import arviz as az
  from sklearn.datasets import load_iris
  import pyprobml_utils as pml
  
  iris = load_iris()
  X = iris.data 
  y = iris.target
  
  # Convert to pandas dataframe 
  df_iris = pd.DataFrame(data=iris.data, 
                      columns=['sepal_length', 'sepal_width', 
                               'petal_length', 'petal_width'])
  df_iris['species'] = pd.Series(iris.target_names[y], dtype='category')
  
  
  df = df_iris.query("species == ('setosa', 'versicolor')") 
  
  # We reduce the sample size from 50 to 25 per class,
  # or to 5 + 45 in the unbalanced setting.
  # The latter will increase posterior uncertainty
  unbalanced = False # True
  if unbalanced:
      df = df[45:95]
  else:
      df = df[25:75]
  assert(len(df)==50)
  
  y_1 = pd.Categorical(df['species']).codes 
  x_n = ['sepal_length', 'sepal_width'] 
  x_1 = df[x_n].values
  
  
  with pm.Model() as model_1: 
      α = pm.Normal('α', mu=0, sd=10) 
      β = pm.Normal('β', mu=0, sd=2, shape=len(x_n)) 
       
      μ = α + pm.math.dot(x_1, β) 
      θ = pm.Deterministic('θ', 1 / (1 + pm.math.exp(-μ))) 
      bd = pm.Deterministic('bd', -α/β[1] - β[0]/β[1] * x_1[:,0])
       
      yl = pm.Bernoulli('yl', p=θ, observed=y_1) 
   
      trace_1 = pm.sample(2000, cores=1, chains=2)
      
  varnames = ['α', 'β'] 
  #az.plot_forest(trace_1, var_names=varnames);
  
  idx = np.argsort(x_1[:,0]) 
  bd = trace_1['bd'].mean(0)[idx] 
  
  plt.figure()
  plt.scatter(x_1[:,0], x_1[:,1], c=[f'C{x}' for x in y_1]) 
  plt.plot(x_1[:,0][idx], bd, color='k'); 
  
  az.plot_hdi(x_1[:,0], trace_1['bd'], color='k')
   
  plt.xlabel(x_n[0]) 
  plt.ylabel(x_n[1])
  
  plt.tight_layout()
  if unbalanced:
      pml.savefig('logreg_iris_bayes_2d_unbalanced.pdf', dpi=300)
  else:
      pml.savefig('logreg_iris_bayes_2d.pdf', dpi=300)
          
  plt.show()
apis: matplotlib.pyplot.xticks matplotlib.pyplot.figure matplotlib.pyplot.grid numpy.random.seed matplotlib.pyplot.title matplotlib.pyplot.subplot matplotlib.pyplot.show matplotlib.pyplot.imshow matplotlib.pyplot.yticks
apis: numpy.sum matplotlib.pyplot.tight_layout numpy.random.seed sklearn.cluster.KMeans numpy.vstack matplotlib.pyplot.style.use matplotlib.pyplot.figure numpy.cos scipy.linalg.eigh sklearn.metrics.pairwise.rbf_kernel numpy.random.uniform numpy.zeros numpy.arange numpy.linalg.norm numpy.empty numpy.random.randn matplotlib.pyplot.show numpy.testing.assert_allclose numpy.sqrt numpy.sin
apis: pandas.Series matplotlib.pyplot.figure matplotlib.pyplot.tight_layout pandas.DataFrame numpy.argsort pandas.Categorical matplotlib.pyplot.show matplotlib.pyplot.ylabel matplotlib.pyplot.plot matplotlib.pyplot.xlabel matplotlib.pyplot.scatter sklearn.datasets.load_iris