<jupyter_start><jupyter_text>CMPE 547 Bayesian Statistics and Machine Learning  CMPE 548 Monte Carlo Methods SWE 546 Data Mining  SWE 582 Machine Learning for Data AnalyticsSupplementary NotesBoğaziçi University, Dept. of Computer EngineeringInstructor: <NAME>  Notebook Summary* We review the notation and parametrization of densities of some basic distributions that are often encountered* We show how random numbers are generated using python libraries* We show some basic visualization methods such as displaying histograms   Sampling From Basic Distributions Sampling from basic distribution is easy using the numpy library.  Formally we will write$x \sim p(X|\theta)$where $\theta$ is the _parameter vector_, $p(X| \theta)$ denotes the _density_ of the random variable $X$ and $x$ is a _realization_, a particular draw from the density $p$.The following distributions are building blocks from which more complicated processes may be constructed. It is important to have a basic understanding of these distributions. Continuous Univariate * Uniform $\mathcal{U}$* Univariate Gaussian $\mathcal{N}$* Gamma $\mathcal{G}$* Inverse Gamma $\mathcal{IG}$* Beta $\mathcal{B}$ Discrete * Poisson $\mathcal{P}$* Bernoulli $\mathcal{BE}$* Binomial $\mathcal{BI}$* Categorical $\mathcal{M}$* Multinomial $\mathcal{M}$ Continuous Multivariate (todo)* Multivariate Gaussian $\mathcal{N}$* Dirichlet $\mathcal{D}$ Continuous Matrix-variate (todo)* Wishart $\mathcal{W}$* Inverse Wishart $\mathcal{IW}$* Matrix Gaussian $\mathcal{N}$  Sampling from standard uniform $\mathcal{U}(0,1)$ For generating a single random number in the interval $[0, 1)$ we use the notation$$x_1 \sim \mathcal{U}(x; 0,1)$$In python, this is implemented as<jupyter_code>import numpy as np

x_1 = np.random.rand()

print(x_1)<jupyter_output>0.9944312250274825<jupyter_text>We can also generate an array of realizations $x_i$ for $i=1 \dots N$,$$x_i \sim \mathcal{U}(x; 0,1)$$<jupyter_code>import numpy as np

N = 5
x = np.random.rand(N)

print(x)<jupyter_output>[ 0.21176772  0.20066043  0.26891193  0.5395703   0.62822017]<jupyter_text>For large $N$, it is more informative to display an histogram of generated data:<jupyter_code>%matplotlib inline
import matplotlib.pyplot as plt

# Number of realizations
N = 1000
x = np.random.rand(N)

plt.hist(x, bins=20)
plt.xlabel('x')
plt.ylabel('Count')
plt.show()<jupyter_output><empty_output><jupyter_text>$\newcommand{\indi}[1]{\left[{1}\right]}$$\newcommand{\E}[1]{\left\langle{1}\right\rangle}$We know that the density of the uniform distribution $\mathcal{U}(0,1)$ is$$\mathcal{U}(x; 0,1) = \left\{ \begin{array}{cc}  1 & 0 \leq x < 1 \\ 0 & \text{otherwise} \end{array}  \right.$$or using the indicator notation$$\mathcal{U}(x; 0,1) = \left[ x \in [0,1) \right]$$ Indicator functionTo write and manipulate discrete probability distributions in algebraic expression, the *indicator* function is useful:$$ \left[x\right] = \left\{ \begin{array}{cc}                                  1 & x\;\;\text{is true} \\                                  0 & x\;\;\text{is false}                                \end{array} \right.$$This notation is also known as the Iverson's convention.  How to plot the density and the histogram onto the same plot? In one dimension, the histogram is simply the count of the data points that fall to a given interval. Mathematically, we have$j = 1\dots J$ intervals where $B_j = [b_{j-1}, b_j]$ and $b_j$ are bin boundries such that $b_0 < b_1 < \dots < b_J$. $$h(x) =  \sum_{j=1}^J \sum_{i=1}^N \indi{x \in B_j} \indi{x_i \in B_j}$$This expression, at the first sight looks somewhat more complicated than it really is. The indicator product just encodes the logical condition $x \in B_j$ __and__ $x_i \in B_j$. The sum over $j$ is just a convenient way of writing the result instead of specifying the histogram as a case by case basis for each bin. It is important to get used to such nested sums.When the density $p(x)$ is given, the probability that a single realization is in bin $B_j$ is given by$$\Pr\left\{x \in B_j\right\} = \int_{B_j} dx p(x) = \int_{-\infty}^{\infty} dx \indi{x\in B_j} p(x) = \E{\indi{x\in B_j}}$$In other words, the probability is just the expectation of the indicator.The histogram can be written as follows $$h(x) =  \sum_{j=1}^J \indi{x \in B_j} \sum_{i=1}^N \indi{x_i \in B_j}$$We define the counts at each bin as $$c_j \equiv \sum_{i=1}^N \indi{x_i \in B_j}$$If all bins have the same width, i.e., $b_j - b_{j-1} = \Delta$ for $\forall j$, and if $\Delta$ is sufficiently small we have$$\E{\indi{x\in B_j}} \approx p(b_{j-1}+\Delta/2) \Delta$$i.e., the probability is roughly the interval width times the density evaluated at the middle point of the bin. The expected value of the counts is$$\E{c_j} = \sum_{i=1}^N \E{\indi{x_i \in B_j}} \approx N \Delta p(b_{j-1}+\Delta/2) $$Hence, the density should be roughly $$p(b_{j-1}+\Delta/2) \approx \frac{\E{c_j} }{N \Delta}$$The $N$ term is intuitive but the $\Delta$ term is easily forgotten. When plotting the histograms on top of the corresponding densities, we should scale the normalized histogram ${ c_j }/{N}$ by dividing by $\Delta$.<jupyter_code>N = 1000

# Bin width
Delta = 0.02

# Bin edges
b = np.arange(0 ,1+Delta, Delta)

# Evaluate the density
g = np.ones(b.size)

# Draw the samples
u = np.random.rand(N)

counts,edges = np.histogram(u, bins=b)
plt.bar(b[:-1], counts/N/Delta, width=Delta)
plt.hold(True)
plt.plot(b, g, linewidth=3, color='y')
plt.hold(False)

plt.show()<jupyter_output><empty_output><jupyter_text>The __plt.hist__ function (calling __np.histogram__) can do this calculation automatically if the option normed=True. However, when the grid is not uniform, it is better to write your own code to be sure what is going on.<jupyter_code>N = 1000
Delta = 0.05
b = np.arange(0 ,1+Delta, Delta)
g = np.ones(b.size)

u = np.random.rand(N)

plt.hist(u, bins=b, normed=True)
plt.hold(True)
plt.plot(b, g, linewidth=3, color='y')
plt.hold(False)

plt.show()<jupyter_output><empty_output><jupyter_text>Sampling from Continuous Univariate Distributions* Uniform $\mathcal{U}$* Univariate Gaussian $\mathcal{N}$$${\cal N}(x;\mu, v) =  \frac{1}{\sqrt{2\pi v}} \exp\left(-\frac12 \frac{(x - \mu)^2}{v}\right)  $$ * Gamma $\mathcal{G}$$${\cal G}(\lambda; a, b) =  \frac{b^a \lambda^{a-1}}{\Gamma(a)} \exp( - b \lambda) = \exp((a-1)\log\lambda - b \lambda - \log\Gamma(a) + a\log b)  $$* Inverse Gamma $\mathcal{IG}$$${\cal IG}(v; \alpha, \beta) =  \frac{\beta^\alpha}{\Gamma(\alpha) v^{\alpha+1}}  \exp(- \frac{\beta}{v}) = \exp(-(\alpha+1)\log v - \beta /v - \log\Gamma(\alpha) + \alpha \log \beta)$$* Beta $\mathcal{B}$$${\cal B}(r; \alpha, \beta) =  \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha) \Gamma(\beta) } r^{\alpha-1} (1-r)^{\beta-1}$$ We will illustrate two alternative ways for sampling from continuous distributions. - The first method has minimal dependence on the numpy and scipy libraries. This is initially the preferred method. Only random variable generators and the $\log \Gamma(x)$ (__gammaln__) function is used and nothing more.- The second method uses scipy. This is a lot more practical but requires knowing more about the internals of the library.  Aside: Gamma function $\Gamma(x)$The gamma function $\Gamma(x)$ is the (generalized) factorial.- Defined by $$\Gamma(x) = \int_0^{\infty} t^{x-1} e^{-t}\, dt$$- For integer $x$, $\Gamma(x) = (x-1)!$. Remember that for positive integers $x$, the factorial function can be defined recursively $x! = (x-1)! x $ for $x\geq 1$.- For real $x>1$, the gamma function satisfies$$\Gamma(x+1) = \Gamma(x) x$$- Interestingly, we have $$\Gamma(1/2) = \sqrt{\pi}$$- Hence$$\Gamma(3/2) = \Gamma(1/2 + 1) = \Gamma(1/2) (1/2) = \sqrt{\pi}/2$$- It is available in many numerical computation packages, in python it is available as __scipy.special.gamma__.- To compute $\log \Gamma(x)$, you should always use the implementation as __scipy.special.gammaln__. The gamma function blows up super-exponentially so numerically you should never evaluate $\log \Gamma(x)$ as```pythonimport numpy as npimport scipy.special as spsnp.log(sps.gamma(x))   Don't sps.gammaln(x)         Do```- A related function is the Beta function$$B(x,y) = \int_0^{1} t^{x-1} (1-t)^{y-1}\, dt$$- We have$$B(x,y) = \frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)}$$- Both $\Gamma(x)$ and $B(x)$ pop up as normalizing constant of the gamma and beta distributions. Derivatives of $\Gamma(x)$-   The derivatives of  $\log \Gamma(x)$ pop up quite often when fitting densities. The first derivative has a specific name,  often called the digamma function or the psi function.$$\Psi(x) \equiv \frac{d}{d x} \log \Gamma(x) $$- It is available as __scipy.special.digamma__ or __scipy.special.psi__- Higher order derivatives of the $\log \Gamma(x)$ function (including digamma itself) are available as __scipy.special.polygamma__<jupyter_code>import numpy as np
import scipy.special as sps

x = np.arange(0.1,5,0.01)

f = sps.gammaln(x)
df = sps.psi(x)

# First derivative of the digamma function
ddf = sps.polygamma(1,x)
# sps.psi(x) == sps.polygamma(0,x)

plt.figure(figsize=(8,10))
plt.subplot(3,1,1)
plt.plot(x, f, 'r')
plt.xlabel('x')
plt.ylabel('log Gamma(x)')

plt.subplot(3,1,2)
plt.grid(True)
plt.plot(x, df, 'b')
plt.xlabel('x')
plt.ylabel('Psi(x)')

plt.subplot(3,1,3)
plt.plot(x, ddf, 'k')
plt.xlabel('x')
plt.ylabel('Psi\'(x)')

plt.show()<jupyter_output><empty_output><jupyter_text>Stirling's approximationAn important approximation to the factorial is the famous Stirling's approximation\begin{align}n! \sim \sqrt{2 \pi n}\left(\frac{n}{e}\right)^n\end{align}\begin{align}\log \Gamma(x+1) \approx \frac{1}{2}\log(2 \pi) + x \log(x) - \frac{1}{2} \log(x) \end{align}  Sampling using numpy.random<jupyter_code>%matplotlib inline
import matplotlib.pyplot as plt

import numpy as np
from scipy.special import gammaln


def plot_histogram_and_density(N, c, edges, dx, g, title='Put a title'):
    '''
     N  : Number of Datapoints
     c  : Counts, as obtained from np.histogram function
     edges : bin edges, as obtained from np.histogram
     dx : The bin width
     g  : Density evaluated at the points given in edges
     title : for the plot
    '''
    plt.bar(edges[:-1], c/N/dx, width=dx)
    plt.hold(True)
    plt.plot(edges, g, linewidth=3, color='y')
    plt.hold(False)
    plt.title(title)

def log_gaussian_pdf(x, mu, V):
    return -0.5*np.log(2*np.pi*V) -0.5*(x-mu)**2/V

def log_gamma_pdf(x, a, b):
    return (a-1)*np.log(x) - b*x - gammaln(a) + a*np.log(b)

def log_invgamma_pdf(x, a, b):
    return -(a+1)*np.log(x) - b/x - gammaln(a) + a*np.log(b)

def log_beta_pdf(x, a, b):
    return - gammaln(a) - gammaln(b) + gammaln(a+b) + np.log(x)*(a-1) + np.log(1-x)*(b-1) 


N = 1000

# Univariate Gaussian
mu = 2    # mean
V = 1.2   # Variance
x_normal = np.random.normal(mu, np.sqrt(V), N)

dx = 10*np.sqrt(V)/50
x = np.arange(mu-5*np.sqrt(V) ,mu+5*np.sqrt(V),dx)
g = np.exp(log_gaussian_pdf(x, mu, V))

#g = scs.norm.pdf(x, loc=mu, scale=np.sqrt(V))
c,edges = np.histogram(x_normal, bins=x)

plt.figure(num=None, figsize=(16, 5), dpi=80, facecolor='w', edgecolor='k')
plt.subplot(2,2,1)
plot_histogram_and_density(N, c, x, dx, g, 'Gaussian')


## Gamma
# Shape
a = 1.2
# inverse scale
b = 30

# Generate unit scale first than scale with inverse scale parameter b
x_gamma = np.random.gamma(a, 1, N)/b
dx = np.max(x_gamma)/500
x = np.arange(dx, 250*dx, dx)

g = np.exp(log_gamma_pdf(x, a, b))
c,edges = np.histogram(x_gamma, bins=x)

plt.subplot(2,2,2)
plot_histogram_and_density(N, c, x, dx, g, 'Gamma')


## Inverse Gamma
a = 3.5
b = 0.2

x_invgamma = b/np.random.gamma(a, 1, N)
dx = np.max(x_invgamma)/500
x = np.arange(dx, 150*dx, dx)
g = np.exp(log_invgamma_pdf(x,a,b))
c,edges = np.histogram(x_invgamma, bins=x)

plt.subplot(2,2,3)
plot_histogram_and_density(N, c, x, dx, g, 'Inverse Gamma')

## Beta 
a = 0.5
b = 1
x_beta = np.random.beta(a, b, N)
dx = 0.01
x = np.arange(dx, 1, dx)
g = np.exp(log_beta_pdf(x, a, b))
c,edges = np.histogram(x_beta, bins=x)

plt.subplot(2,2,4)
plot_histogram_and_density(N, c, x, dx, g, 'Beta')

plt.show()<jupyter_output><empty_output><jupyter_text>Sampling using scipy.stats<jupyter_code>%matplotlib inline
import matplotlib.pyplot as plt

import numpy as np
import scipy.stats as scs



N = 2000

# Univariate Gaussian
mu = 2    # mean
V = 1.2   # Variance

rv_normal = scs.norm(loc=mu, scale=np.sqrt(V))
x_normal = rv_normal.rvs(size=N)

dx = 10*np.sqrt(V)/50
x = np.arange(mu-5*np.sqrt(V) ,mu+5*np.sqrt(V),dx)
g = rv_normal.pdf(x)

c,edges = np.histogram(x_normal, bins=x)

plt.figure(num=None, figsize=(16, 5), dpi=80, facecolor='w', edgecolor='k')
plt.subplot(2,2,1)
plot_histogram_and_density(N, c, x, dx, g, 'Gaussian')


## Gamma
a = 3.2
b = 30

# The following is equivalent to our parametrization of gamma, note the 1/b term
rv_gamma = scs.gamma(a, scale=1/b)
x_gamma = rv_gamma.rvs(N)
dx = np.max(x_gamma)/500
x = np.arange(0, 250*dx, dx)
g = rv_gamma.pdf(x)
c,edges = np.histogram(x_gamma, bins=x)

plt.subplot(2,2,2)
plot_histogram_and_density(N, c, x, dx, g, 'Gamma')



## Inverse Gamma
a = 3.5
b = 0.2

# Note the b term
rv_invgamma = scs.invgamma(a, scale=b)
x_invgamma = rv_invgamma.rvs(N)
dx = np.max(x_invgamma)/500
x = np.arange(dx, 150*dx, dx)
g = rv_invgamma.pdf(x)
c,edges = np.histogram(x_invgamma, bins=x)

plt.subplot(2,2,3)
plot_histogram_and_density(N, c, x, dx, g, 'Inverse Gamma')

## Beta 
a = 0.7
b = 0.8

rv_beta = scs.beta(a, b)
x_beta = rv_beta.rvs(N)
dx = 0.02
x = np.arange(0, 1+dx, dx)
g = rv_beta.pdf(x)
c,edges = np.histogram(x_beta, bins=x)

plt.subplot(2,2,4)
plot_histogram_and_density(N, c, x, dx, g, 'Beta')

plt.show()<jupyter_output><empty_output><jupyter_text>Sampling from Discrete Densities* Bernoulli $\mathcal{BE}$$${\cal BE}(r; w)  =  w^r (1-w)^{1-r} \;\; \text{if} \; r \in \{0, 1\} $$* Binomial $\mathcal{BI}$$${\cal BI}(r; L, w)  =   \binom{L}{r, (L-r)} w^r (1-w)^{L-r} \;\; \text{if} \; r \in \{0, 1, \dots, L\} $$Here, the binomial coefficient is defined as $$\binom{L}{r, (L-r)} = \frac{N!}{r!(L-r)!}$$Note that $${\cal BE}(r; w) = {\cal BI}(r; L=1, w)$$* Poisson $\mathcal{PO}$, with intensity $\lambda$$${\cal PO}(x;\lambda) =  \frac{e^{-\lambda} \lambda^x}{x!} = \exp(x \log \lambda - \lambda - \log\Gamma(x+1)) $$  Given samples on nonnegative integers, we can obtain histograms easily using __np.bincount__. ```pythonc = np.bincount(samples)```The functionality is equivalent to the following sniplet, while implementation is possibly different and more efficient.```pythonupper_bound = np.max()c = np.zeros(upper_bound+1)for i in samples:    c[i] += 1```<jupyter_code>%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np

def plot_histogram_and_pmf(N, c, domain, dx, g, title='Put a title'):
    '''
     N  : Number of Datapoints
     c  : Counts, as obtained from np.bincount function
     domain : integers for each c, same size as c
     dx : The bin width
     g  : Density evaluated at the points given in edges
     title : for the plot
    '''
    plt.bar(domain-dx/2, c/N, width=dx)
    plt.hold(True)
    plt.plot(domain, g, 'ro:', linewidth=3, color='y')
    plt.hold(False)
    plt.title(title)

def log_poisson_pdf(x, lam):
    return -lam + x*np.log(lam) - gammaln(x+1)

def log_bernoulli_pdf(r, pr):
    return r*np.log(pr) + (1-r)*np.log(1 - pr)

def log_binomial_pdf(r, pr, L):
    return gammaln(L+1) - gammaln(r+1) - gammaln(L-r+1) + r*np.log(pr) + (L-r)*np.log(1 - pr)

    
N = 100
pr = 0.8

# For plots
bin_width = 0.3

# Bernoulli
L = 1
x_bern = np.random.binomial(n=L, p=pr, size=N)
c = np.bincount(x_bern, minlength=L+1)
g = np.exp(log_bernoulli_pdf(np.arange(L+1), pr))

plt.figure(figsize=(20,4))
plt.subplot(1,3,1)
plot_histogram_and_pmf(N, c, np.arange(L+1), bin_width, g, 'Bernoulli')
plt.xticks([0,1])


# Binomial
L = 10
pr = 0.7
x_binom = np.random.binomial(n=L, p=pr, size=N)
c = np.bincount(x_binom, minlength=L+1)
g = np.exp(log_binomial_pdf(np.arange(L+1), pr, L))

plt.subplot(1,3,2)
plot_histogram_and_pmf(N, c, np.arange(L+1), bin_width, g, 'Binomial')
plt.xticks(np.arange(L+1))

# Poisson
intensity = 10.5
x_poiss = np.random.poisson(intensity, size =N )
c = np.bincount(x_poiss)

x = np.arange(len(c))
g = np.exp(log_poisson_pdf(x, intensity))

plt.subplot(1,3,3)
plot_histogram_and_pmf(N, c, x, bin_width, g, 'Poisson')<jupyter_output><empty_output><jupyter_text>MathJax.Hub.Config({ TeX: { extensions: ["AMSmath.js"] }}); Bernoulli, Binomial, Categorical and Multinomial DistributionsThe Bernoulli and Binomial distributions are quite simple and well known distributions on small integers, so it may come as a surprise that they have another, less obvious but arguably more useful representation as discrete multivariate densities. This representation makes the link to categorical distributions where there are more than two possible outcomes. Finally, all Bernoulli, Binomial or Categorical distributions are special cases of Multinomial distribution. BernoulliRecall the Bernoulli distribution $r \in \{0, 1\}$$${\cal BE}(r; w)  =  w^r (1-w)^{1-r}  $$We will define $\pi_0 = 1-w$ and $\pi_1 = w$, such that $\pi_0 + \pi_1 = 1$. The parameter vector is $\pi = (\pi_0, \pi_1)$We will also introduce a positional encoding such that \begin{eqnarray}r = 0 & \Rightarrow & s = (1, 0) \\r = 1 & \Rightarrow & s = (0, 1) \end{eqnarray}In other words $s = (s_0, s_1)$ is a 2-dimensional vector where $$s_0, s_1 \in \{0,1\}\;\text{and}\; s_0 + s_1 = 1$$ We can now write the Bernoulli density $$p(s | \pi)  =  \pi_0^{s_0} \pi_1^{s_1}$$ BinomialSimilarly, recall the Binomial density where $r \in \{0, 1, \dots, L\}$$${\cal BI}(r; L, w)  =   \binom{L}{r, (L-r)} w^r (1-w)^{L-r}  $$We will again define $\pi_0 = 1-w$ and $\pi_1 = w$, such that $\pi_0 + \pi_1 = 1$. The parameter vector is $\pi = (\pi_0, \pi_1)$\begin{eqnarray}r = 0 & \Rightarrow & s = (L, 0) \\r = 1 & \Rightarrow & s = (L-1, 1)\\r = 2 & \Rightarrow & s = (L-2, 2)\\\dots \\r = L & \Rightarrow & s = (0, L)\end{eqnarray}where $s = (s_0, s_1)$ is a 2-dimensional vector where $$s_0, s_1 \in \{0,\dots,L\} \;\text{and}\; s_0 + s_1  = L$$ We can now write the Binomial density as$$p(s | \pi)  =  \binom{L}{s_0, s_1} \pi_0^{s_0} \pi_1^{s_1}$$ CategoricalOne of the advantages of this new notation is that we can write the density even if the outcomes are not numerical. For example, the result of a single coin flip experiment when $r \in \{$ 'Tail', 'Head' $\}$ where the probability of 'Tail' is $w$ can be written as$$p(r | w)  =  w^{\indi{r=\text{'Tail'}}} (1-w)^{\indi{r=\text{'Head'}}}$$We define $s_0 = \indi{r=\text{'Head'}}$ and $s_1 = \indi{r=\text{'Tail'}}$, then the density can be written in the same form as$$p(s | \pi)  = \pi_0^{s_0} \pi_1^{s_1}$$where $\pi_0 = 1-w$ and $\pi_1 = w$.More generally, when $r$ is from a set with $K$ elements, i.e., $r \in R = \{ v_0, v_1, \dots, v_{K-1} \}$ with probability of the event $r = v_k$ given as $\pi_k$, we define $s = (s_0, s_1, \dots, s_{K-1})$ for $k=0,\dots, K-1$$$s_k = \indi{r=v_k}$$Note that by construction, we have $\sum_k s_k = 1$. The resulting density, known as the Categorical density, can be writen as $$p(s|\pi) = \pi_0^{s_0} \pi_1^{s_1} \dots \pi_{K-1}^{s_{K-1}}$$ Multinomial When drawing from a categorical distribution, one chooses a single category from $K$ options with given probabilities. A standard model for this is placing a single ball into $K$ different bins. The vector $s = (s_0, s_1, \dots,s_k, \dots, s_{K-1})$ represents how many balls eack bin $k$ contains.Now, place $L$ balls instead of one into $K$ bins with placing each ball idependently into bin $k$ where $k \in\{0,\dots,K-1\}$ with the probability $\pi_k$. The multinomial is the joint distribution of $s$ where $s_k$ is the number of balls placed into bin $k$. The density will be denoted as $${\cal M}(s; L, \pi)  =  \binom{L}{s_0, s_1, \dots, s_{K-1}}\prod_{k=0}^{K-1} \pi_k^{s_k} $$Here $\pi \equiv [\pi_0, \pi_2, \dots, \pi_{K-1} ]$ is the probability vector and $L$ is referred as the _index parameter_.Clearly, we have the normalization constraint $ \sum_k \pi_k = 1$ and realization of the counts $s$ satisfy $ \sum_k s_k = L $. Here, the _multinomial_ coefficient is defined as $$\binom{L}{s_0, s_1, \dots, s_{K-1}} = \frac{L!}{s_0! s_1! \dots s_{K-1}!}$$Binomial, Bernoulli and Categorical distributions are all special cases of the Multinomial distribution, with a suitable representation.The picture is as follows:Balls/Bins | $2$ Bins | $K$ Bins  -------- | -------- | ---------$1$ Ball | Bernoulli ${\cal BE}$ | Categorical  ${\cal C}$$L$ Balls | Binomial ${\cal BI}$ | Multinomial ${\cal M}$Murphy calls the categorical distribution ($1$ Ball, $K$ Bins) as the Multinoulli. This is non-standard but logical (and somewhat cute).It is common to think of Bernoulli and Binomial as scalar random variable. However, when we think of them as special case of a Multinomial it is better to think of them as bivariate, albeit degenerate, random variables, as illustrated in the following cell along with an alternative visualization.<jupyter_code># The probability parameter
pr = 0.3

fig = plt.figure(figsize=(16,50), edgecolor=None)

maxL = 12
plt.subplot(maxL-1,2,1)
plt.grid(False)
# Set up the scalar binomial density as a bivariate density  
for L in range(1,maxL):
    r = np.arange(L+1)
    p = np.exp(log_binomial_pdf(r, pr=pr, L=L))
    A = np.zeros(shape=(13,13))
    for s in range(L):
        s0 = s
        s1 = L-s
        A[s0, s1] = p[s]

    #plt.subplot(maxL-1,2,2*L-1)
#    plt.bar(r-0.25, p, width=0.5)
#    ax.set_xlim(-1,maxL)
#    ax.set_xticks(range(0,maxL))
    if True:
        plt.subplot(maxL-1,2,2*L-1)
        plt.barh(bottom=r-0.25, width=p, height=0.5)
        ax2 = fig.gca()
        pos = ax2.get_position()
        pos2 = [pos.x0, pos.y0, 0.04, pos.height]
        ax2.set_position(pos2)
        ax2.set_ylim(-1,maxL)
        ax2.set_yticks(range(0,maxL))
        ax2.set_xlim([0,1])
        ax2.set_xticks([0,1])
        plt.ylabel('s1')
        ax2.invert_xaxis()
    
    plt.subplot(maxL-1,2,2*L)
    plt.imshow(A,  interpolation='nearest', origin='lower',cmap='gray_r',vmin=0,vmax=0.7)
    plt.xlabel('s0')
        
    ax1 = fig.gca()
    
    pos = ax1.get_position()
    pos2 = [pos.x0-0.45, pos.y0, pos.width, pos.height]
    ax1.set_position(pos2)

    ax1.set_ylim(-1,maxL)
    ax1.set_yticks(range(0,maxL))
    ax1.set_xlim(-1,maxL)
    ax1.set_xticks(range(0,maxL))


    
plt.show()<jupyter_output><empty_output><jupyter_text>The following cell illustrates sampling from the Multinomial density.<jupyter_code># Number of samples
sz = 3

# Multinomial
p = np.array([0.3, 0.1, 0.1, 0.5])
K = len(p) # number of Bins
L = 20     # number of Balls

print('Multinomial with number of bins K = {K} and Number of balls L = {L}'.format(K=K,L=L))
print(np.random.multinomial(L, p, size=sz))

# Categorical 
L = 1     # number of Balls

print('Categorical with number of bins K = {K} and a single ball L=1'.format(K=K))
print(np.random.multinomial(L, p, size=sz))

# Binomial
p = np.array([0.3, 0.7])
K = len(p) # number of Bins = 2
L = 20     # number of Balls

print('Binomial with two bins K=2 and L={L} balls'.format(L=L))
print(np.random.multinomial(L, p, size=sz))

# Bernoulli
L = 1     # number of Balls
p = np.array([0.3, 0.7])
K = len(p) # number of Bins = 2

print('Bernoulli, two bins and a single ball')
print(np.random.multinomial(L, p, size=sz))<jupyter_output>Multinomial with number of bins K = 4 and Number of balls L = 20
[[ 5  4  3  8]
 [ 7  4  1  8]
 [ 5  3  0 12]]
Categorical with number of bins K = 4 and a single ball L=1
[[1 0 0 0]
 [0 0 0 1]
 [0 0 0 1]]
Binomial with two bins K=2 and L=20 balls
[[10 10]
 [ 5 15]
 [ 8 12]]
Bernoulli, two bins and a single ball
[[0 1]
 [0 1]
 [0 1]]<jupyter_text>Probability tables and the categorical distributionThe following cell illustrates drawing from a categorical distribution with on an alphabet, not necessarly $0\dots K-1$.<jupyter_code># Sampling from a Categorical Distribution
a = np.array(sorted(['blue', 'red', 'black', 'yellow']))
pr = np.array([0.2, 0.55, 0.15, 0.1])
N = 100
x = np.random.choice(a, size=N, replace=True, p=pr)


print('Symbols:')
print(a)
print('Probabilities:')
print(pr)
print('{N} realizations:'.format(N=N))
print(x)<jupyter_output>Symbols:
['black' 'blue' 'red' 'yellow']
Probabilities:
[ 0.2   0.55  0.15  0.1 ]
100 realizations:
['blue' 'blue' 'blue' 'yellow' 'black' 'blue' 'yellow' 'blue' 'red' 'blue'
 'blue' 'black' 'blue' 'yellow' 'blue' 'blue' 'red' 'black' 'blue' 'blue'
 'yellow' 'black' 'blue' 'red' 'blue' 'black' 'red' 'red' 'black' 'blue'
 'yellow' 'blue' 'blue' 'blue' 'blue' 'blue' 'red' 'blue' 'blue' 'blue'
 'blue' 'red' 'blue' 'blue' 'blue' 'blue' 'blue' 'blue' 'blue' 'black'
 'yellow' 'blue' 'blue' 'black' 'blue' 'red' 'blue' 'blue' 'blue' 'red'
 'red' 'blue' 'blue' 'blue' 'blue' 'blue' 'red' 'blue' 'blue' 'blue' 'blue'
 'blue' 'blue' 'yellow' 'black' 'red' 'blue' 'black' 'red' 'blue' 'blue'
 'blue' 'blue' 'red' 'blue' 'blue' 'blue' 'black' 'blue' 'blue' 'blue'
 'blue' 'blue' 'black' 'red' 'yellow' 'blue' 'black' 'blue' 'blue']<jupyter_text>Often we need to opposite of the above process, that is given a list of elements, we need to count the number of occurences of each symbol. The following method creates such a statistics.<jupyter_code>import collections

c = collections.Counter(x)
print(c.most_common())

counts = [e[1] for e in c.most_common()]
symbols = [e[0] for e in c.most_common()]

print('Sorted according to counts')
print(counts)
print(symbols)

# If we require the symbols in sorted order with respect to symbol names, use:
counts = [e[1] for e in sorted(c.most_common())]
symbols = [e[0] for e in sorted(c.most_common())]

print('Sorted according to symbols')
print(counts)
print(symbols)<jupyter_output>[('blue', 64), ('red', 15), ('black', 13), ('yellow', 8)]
Sorted according to counts
[64, 15, 13, 8]
['blue', 'red', 'black', 'yellow']
Sorted according to symbols
[13, 64, 15, 8]
['black', 'blue', 'red', 'yellow']<jupyter_text>Counting letter bigrams in several languages<jupyter_code>%matplotlib inline
from collections import defaultdict
from urllib.request import urlopen
import string
import numpy as np
import matplotlib.pyplot as plt

# Turkish
#"ç","ı","ğ","ö","ş","ü",'â'
# German
#"ä","ß","ö","ü"
# French
#"ù","û","ô","â","à","ç","é","è","ê","ë","î","ï","æ"
tr_alphabet = ['•','a','b','c','ç','d','e','f',
                'g','ğ','h','ı','i','j','k','l',
                'm','n','o','ö','p','q','r','s','ş',
                't','u','ü','w','v','x','y','z']
# Union of Frequent letters in French, Turkish, German and English
my_alphabet = ['•','a','â','ä',"à","æ",'b','c','ç','d','e',"é","è","ê","ë",'f',
                'g','ğ','h','ı','i',"î",'ï','j','k','l',
                'm','n','o','œ',"ô",'ö','p','q','r','s','ş',
                't','u','ù',"û",'ü','w','v','x','y','z','ß']
# Only ascii characters
ascii_alphabet = list('•'+string.ascii_lowercase)
# Reduction table from my alphabet to ascii
my2ascii_table = {
    ord('â'):"a",
    ord('ä'):"ae",
    ord("à"):"a",
    ord("æ"):"ae",
    ord('ç'):"c",
    ord("é"):"e",
    ord("è"):"e",
    ord("ê"):"e",
    ord("ë"):"e",
    ord('ğ'):"g",
    ord('ı'):"i",
    ord("î"):"i",
    ord('ï'):"i",
    ord('œ'):"oe",
    ord("ô"):"o",
    ord('ö'):"o",
    ord('ş'):"s",
    ord('ù'):"u",
    ord("û"):"u",
    ord('ü'):"u",
    ord('ß'):"ss"
    }
# Reduction table from my alphabet to frequent letters in turkish text
my2tr_table = {
    ord('â'):"a",
    ord('ä'):"ae",
    ord("à"):"a",
    ord("æ"):"ae",
    ord("é"):"e",
    ord("è"):"e",
    ord("ê"):"e",
    ord("ë"):"e",
    ord("î"):"i",
    ord('ï'):"i",
    ord('œ'):"oe",
    ord("ô"):"o",
    ord('ù'):"u",
    ord("û"):"u",
    ord('ß'):"ss"
    }


def count_transitions(fpp, alphabet, tab):

    #ignore punctuation
    tb = str.maketrans(".\t\n\r ","•••••", '0123456789!"\'#$%&()*,-/:;<=>?@[\\]^_`{|}~+')
    #replace other unicode characters with a bullet (alt-8)
    tbu = {
        ord("İ"):'i',
        ord(u"»"):'•', 
        ord(u"«"):'•', 
        ord(u"°"):'•', 
        ord(u"…"):'•',
        ord(u"”"):'•',
        ord(u"’"):'•',
        ord(u"“"):'•',
        ord(u"\ufeff"):'•',
        775: None}

    # Character pairs 
    D = defaultdict(int)

    for line in fpp:
        s = line.decode('utf-8').translate(tb).lower()
        s = s.translate(tbu)
        s = s.translate(tab)
        #print(s)

        if len(s)>1:
            for i in range(len(s)-1):
                D[s[i:i+2]]+=1
    
    M = len(alphabet)
    a2i = {v: k for k,v in enumerate(alphabet)}
    DD = np.zeros((M,M))

    ky = sorted(D.keys())
    for k in D.keys():
        i = a2i[k[0]]
        j = a2i[k[1]]
        DD[i,j] = D[k]

    return D, DD, alphabet<jupyter_output><empty_output><jupyter_text>Count and display occurences of letters in text<jupyter_code>local = 'file:///Users/cemgil/Dropbox/Public/swe546/data/'
#local = 'https://dl.dropboxusercontent.com/u/9787379/swe546/data/'
#files = ['starwars_4.txt', 'starwars_5.txt', 'starwars_6.txt', 'hamlet.txt', 'hamlet_deutsch.txt', 'hamlet_french.txt', 'juliuscaesar.txt','othello.txt', 'sonnets.txt', 'antoniusandcleopatra.txt']
files = ['hamlet_turkce.txt','hamlet_deutsch.txt', 'hamlet_french.txt', 'hamlet.txt','starwars_4.txt', 'starwars_5.txt','juliuscaesar.txt','othello.txt']

plt.figure(figsize=(16,18))

i = 0
for f in files:
    url = local+f
    data = urlopen(url) 
    #D, DD, alphabet = count_transitions(data, my_alphabet, {})
    D, DD, alphabet = count_transitions(data, ascii_alphabet, my2ascii_table)
    #D, DD, alphabet = count_transitions(data, tr_alphabet, my2tr_table)
    M = len(alphabet)
    # Ignore space, space transitions
    DD[0,0] = 1
    
    i+=1
    plt.subplot(len(files)/2,2,i)
        
    S = np.sum(DD,axis=0)
    #Subpress spaces
    S[0] = 0
    S = S/np.sum(S)

    plt.bar(np.arange(M)-0.5, S, width=0.7)
    plt.xticks(range(M), alphabet)
    plt.gca().set_ylim((0,0.2))
    plt.title(f)

plt.show()<jupyter_output><empty_output><jupyter_text>Counting Bigrams<jupyter_code>local = 'file:///Users/cemgil/Dropbox/Public/swe546/data/'
#local = 'https://dl.dropboxusercontent.com/u/9787379/swe546/data/'

#files = ['starwars_4.txt', 'starwars_5.txt', 'starwars_6.txt', 'hamlet.txt', 'hamlet_deutsch.txt', 'hamlet_french.txt', 'juliuscaesar.txt','othello.txt', 'sonnets.txt', 'antoniusandcleopatra.txt']
files = ['hamlet_turkce.txt','hamlet_deutsch.txt', 'hamlet_french.txt', 'hamlet.txt','starwars_4.txt', 'starwars_5.txt','juliuscaesar.txt','othello.txt']

plt.figure(figsize=(17,2*17))

i = 0
for f in files:
    url = local+f
    data = urlopen(url) 
    #D, DD, alphabet = count_transitions(data, my_alphabet, {})
    #D, DD, alphabet = count_transitions(data, ascii_alphabet, my2ascii_table)
    D, DD, alphabet = count_transitions(data, tr_alphabet, my2tr_table)
    M = len(alphabet)
    DD[0,0] = 1
    
    i+=1
    plt.subplot(len(files)/2,2,i)
    plt.imshow(DD, interpolation='nearest', vmin=0)
    plt.xticks(range(M), alphabet)
    plt.xlabel('x(t)')
    plt.yticks(range(M), alphabet)
    plt.ylabel('x(t-1)')
    ax = plt.gca()
    ax.xaxis.tick_top()
    #ax.set_title(f, va='bottom')
    plt.xlabel('x(t) '+f)<jupyter_output><empty_output><jupyter_text>Normalized probability table of $p(x_t|x_{t-1})$<jupyter_code>def normalize(A, axis=0):
    Z = np.sum(A, axis=axis,keepdims=True)
    idx = np.where(Z == 0)
    Z[idx] = 1
    return A/Z

local = 'file:///Users/cemgil/Dropbox/Public/swe546/data/'
#local = 'https://dl.dropboxusercontent.com/u/9787379/swe546/data/'

file = 'hamlet_turkce.txt'
data = urlopen(local+file) 
D, DD, alphabet = count_transitions(data, tr_alphabet, my2tr_table)


plt.figure(figsize=(9,9))

T = normalize(DD, axis=1)


plt.imshow(T, interpolation='nearest', vmin=0)
plt.xticks(range(M), alphabet)
plt.yticks(range(M), alphabet)
plt.gca().xaxis.tick_top()

plt.show()<jupyter_output><empty_output><jupyter_text>Is Markov(0), Markov(1) or Markov(2) a better model for English letters in plain text ?  Continuous Multivariate (todo)  The Multivariate Gaussian distribution\begin{align}\mathcal{N}(x; \mu, \Sigma) &= |2\pi \Sigma|^{-1/2} \exp\left( -\frac{1}{2} (x-\mu)^\top \Sigma^{-1} (x-\mu) \right) \\& = \exp\left(-\frac{1}{2} x^\top \Sigma^{-1} x + \mu^\top \Sigma^{-1} x   - \frac{1}{2} \mu^\top \Sigma^{-1} \mu   -\frac{1}{2}\log \det(2\pi \Sigma) \right) \\\end{align}Draw a vector $x \in \mathbf{R}^N$ where each element $x_i \sim \mathcal{N}(x; 0, 1)$ for $i = 1\dots N$.$\newcommand{\E}[1]{\left\langle1\right\rangle}$Construct \begin{align}y = Ax\end{align}The expectation and the variance are obtained by\begin{align}\E{y} = \E{Ax} = 0\end{align}\begin{align}\E{y y^\top} = A \E{x x^\top} A^\top = A A^\top\end{align}So\begin{align}y \sim \mathcal{N}(y; 0, A A^\top)\end{align} In two dimensions, a bi-variate Gaussian is conveniently represented by an ellipse. The ellipse shows a contour of equal probability. In particular, if we plot the $3\sigma$ ellipse, $99 \%$ of all the data points should be inside the ellipse.<jupyter_code>%matplotlib inline

def ellipse_line(A, mu, col='b'):
    '''
    Creates an ellipse from short line segments y = A x + \mu 
    where x is on the unit circle.
    '''
    N = 18
    th = np.arange(0, 2*np.pi+np.pi/N, np.pi/N)
    X = np.array([np.cos(th),np.sin(th)])
    
    Y = np.dot(A, X)
    ln = plt.Line2D(mu[0]+Y[0,:],mu[1]+Y[1,:],markeredgecolor='k', linewidth=1, color=col)
    return ln

N = 100

A = np.random.randn(2,2)
mu = np.zeros(2)
X = np.random.randn(2,N)
Y = np.dot(A,X)

plt.cla()
plt.axis('equal')
ax = plt.gca()
ax.set_xlim(-8,8)
ax.set_ylim(-8,8)

col = 'b'
ln = ellipse_line(3*A, mu, col)
ax.add_line(ln)
plt.hold(True)
plt.plot(mu[0]+Y[0,:],mu[1]+Y[1,:],'.'+col)

plt.show()
np.dot(A,A.T)<jupyter_output><empty_output><jupyter_text>When the covariance matrix $\Sigma$ is given, as is typically the case, we need a factorization of $\Sigma = W W^\top$. The Cholesky factorization is such a factorization. (Another possibility, whilst computationally more costly, is the matrix square root.)<jupyter_code>Sigma = np.dot(A, A.T)

W = np.linalg.cholesky(Sigma)
X = np.random.randn(2,N)
Y = np.dot(W,X)

plt.cla()
plt.axis('equal')
ax = plt.gca()
ax.set_xlim(-8,8)
ax.set_ylim(-8,8)

col = 'b'
ln = ellipse_line(3*W, mu, col)
ax.add_line(ln)
plt.hold(True)
plt.plot(mu[0]+Y[0,:],mu[1]+Y[1,:],'.'+col)

plt.show()<jupyter_output><empty_output><jupyter_text>The numpy function __numpy.random.multivariate_normal__ generates samples from a multivariate Gaussian with the given mean and covariance.<jupyter_code>N = 100
Sig = np.dot(A, A.T)
x = np.random.multivariate_normal(mu, Sig, size=N)

plt.cla()
plt.axis('equal')
ax = plt.gca()
ax.set_xlim(-8,8)
ax.set_ylim(-8,8)

plt.plot(x[:,0], x[:,1], 'b.')
ln = ellipse_line(3*A,mu,'b')
plt.gca().add_line(ln)

plt.show()<jupyter_output><empty_output><jupyter_text>Evaluation of the multivariate Gaussian densityThe log-density of the multivariate Gaussian has the following exponential form\begin{align}\log \mathcal{N}(x; \mu, \Sigma) &= -\frac{1}{2}\log \det(2\pi \Sigma) -\frac{1}{2} (x-\mu)^\top \Sigma^{-1} (x - \mu)  \end{align}It is tempting to implement these expression as written -- indeed it is useful to do so for debugging purposes. However, this direct method is both inefficient and numerically not very stable. This will be a problem when the dimension of $x$ is high. A direct implementation might be as follows:<jupyter_code>def log_mvar_gaussian_inefficient(x, mu, Sig):
    return -0.5*np.log(np.linalg.det(2*np.pi*Sig)) - 0.5*np.sum((x-mu)*np.dot(np.linalg.inv(Sig), x-mu),axis=0)<jupyter_output><empty_output><jupyter_text>The evaluation seemingly requires the following steps:- Evaluation of the log of the determinant of the covariance matrix $\Sigma$- Inversion of the covariance matrix $\Sigma$ - Evaluation of the quadratic form $(x-\mu)^\top \Sigma^{-1} (x - \mu)$A more efficient implementation uses the following observations:- The covariance matrix $\Sigma$ is positive semidefinite and has a __Cholesky__ factorization \begin{align}\Sigma = W W^\top\end{align}where $W$ is a lower triangular matrix- The determinant satisfies the following identity\begin{align}\det(\Sigma) & = \det(W W^\top) = \det(W) \det(W^\top) = \det(W)^2\end{align}- The determinant of the triangular matrix $W$ is simply the product of its diagonal elements $W_{i,i}$ so \begin{align}\log \det(\Sigma) & = 2 \log \det(W) & = 2 \sum_i \log W_{i,i} \end{align}- The quadratic form can be evaluated by the inner product $(x - \mu)^\top u$ where $u = \Sigma^{-1} (x - \mu)$.Finding $u$ is equivalent to the solution of the linear system $$\Sigma u = (x - \mu) $$and the solution is equivalent to $$ u = (W^\top)^{-1}W^{-1} (x - \mu) $$and can be solved efficiently by backsubstitution as $W$ is triangular. This can be implemented as follows<jupyter_code>import scipy as sc
import scipy.linalg as la

def log_mvar_gaussian_pdf(x, mu, Sig):
    W = np.linalg.cholesky(Sig)
    z = -np.sum(np.log(2*np.pi)/2 + np.log(np.diag(W))) - 0.5* np.sum((x-mu)*la.cho_solve((W,True), x-mu),axis=0)
    return z

# Dimension of the problem 
N = 2

# Generate K points to evaluate the density at
K = 10
x = np.random.randn(N,K)

# Generate random parameters
mu = np.random.randn(N,1) 
R = np.random.randn(N,N)
Sig = np.dot(R, R.T)

z1 = log_mvar_gaussian_pdf(x, mu, Sig)
z2 = log_mvar_gaussian_inefficient(x, mu, Sig)

print(z1)
print(z2)<jupyter_output>[-8.72265771 -6.31543527 -5.4738194  -1.80930712 -2.15963976 -5.54244251
 -3.82657297 -4.16592483 -2.09530621 -8.63078085]
[-8.72265771 -6.31543527 -5.4738194  -1.80930712 -2.15963976 -5.54244251
 -3.82657297 -4.16592483 -2.09530621 -8.63078085]<jupyter_text>For the solution of $\Sigma u = b$ where $\Sigma = WW^\top$, we use the implementation in _scipy.linalg.cho_solve_.<jupyter_code>import scipy.linalg as la

N = 2
# Construct a random positive definite matrix
R = np.random.randn(N,N)
Sig = np.dot(R, R.T)
b = np.random.randn(N,1)

# Direct implementation -- inefficient
u_direct = np.matrix(np.linalg.inv(Sig))*b

# Efficient implementation
W = np.linalg.cholesky(Sig)
u_efficient = la.cho_solve((W,True), b)

# Verify that both give the same result
print(u_direct)
print(u_efficient)<jupyter_output>[[-0.40649925]
 [ 0.39775871]]
[[-0.40649925]
 [ 0.39775871]]<jupyter_text>Dirichlet DistributionThe Dirichlet distribution is a distribution over probability vectors.$$\mathcal{D}(w_{1:N}; \alpha_{1:N}) = \frac{\Gamma(\sum_i \alpha_i)}{\prod_i \Gamma(\alpha_i) } \prod w_i^{\alpha_i-1}$$<jupyter_code>np.random.dirichlet(0.01*np.array([1,2,3]))<jupyter_output><empty_output><jupyter_text>Some useful functions from np.randomWe can get a specific random number state and generate data from it.<jupyter_code>import numpy as np

u = np.random.RandomState()
print(u.permutation(10))

lam = 3;
print(u.exponential(lam))

%matplotlib inline
import matplotlib as mpl
import matplotlib.pyplot as plt

A = np.random.randn(4,5)

plt.figure(figsize=(5,4))
plt.imshow(A, interpolation='nearest')

plt.show()

ax = plt.gca()
ax.set_xlim(-4,4)
ax.set_ylim(-4,4)

A = np.array([[1, 0.3, 0],[0, 1, 0],[0,0,1]])
for i in range(5):
#    circ = plt.Circle(np.random.randn(2),radius=np.random.rand(1)*4,alpha=0.1)
    circ = mpl.patches.Ellipse(np.random.randn(2),width=np.random.rand(1)*4,height=np.random.rand(1)*4, angle=60, alpha=0.1)
    ax.add_patch(circ)

plt.show()<jupyter_output><empty_output><jupyter_text>Importance Sampling Example<jupyter_code>import numpy as np

L = 200
lam = 3.0
rho = 100*L

lp = -L/lam

def log_w(x, lam, rho):
    return np.log(rho)-np.log(lam) - x/lam + x/rho

N = 10000

Inf = np.float('inf')

logS = -Inf
S2 = 0
for i in range(N):
    u = np.random.rand(1)
    # Sample from the proposal
    # Generate x ~ Exp(rho)
    x = -rho*np.log(u)
    x2 = -lam*np.log(u)
    
    if x>L:
        lw = log_w(x, lam, rho)
        m = np.max((logS, lw))
        logS = m + np.log(np.exp(logS-m) + np.exp(lw-m) )

    if x2>L:
        S2 += 1
        
lp_est = logS - np.log(N)
mc_est = np.log(S2) - np.log(N)

print('Ground Truth :', lp)
print('IS Estimator :', lp_est)
print('NI Estimator :', mc_est)<jupyter_output>Ground Truth : -66.66666666666667
IS Estimator : [-66.44972943]
NI Estimator : -inf<jupyter_text>Importance sampling without knowing the normalizing constant<jupyter_code>import numpy as np

L = 200
lam = 3.0
rho = L
k = 1.1

lp = -L/lam


def log_w(x, lam, rho, k):
    return np.log(rho) - np.power((x/lam),k) + x/rho

N = 10000

Inf = np.float('inf')

logS = -Inf
logSW = -Inf
for i in range(N):
    u = np.random.rand(1)
    # Sample from the proposal
    # Generate x ~ Exp(rho)
    x = -rho*np.log(u)
    x2 = -lam*np.log(u)
    
    lw = log_w(x, lam, rho, k)
    m = np.max((logSW, lw))
    logSW = m + np.log(np.exp(logSW-m) + np.exp(lw-m) )
        
    if x>L:
        m = np.max((logS, lw))
        logS = m + np.log(np.exp(logS-m) + np.exp(lw-m) )

lp_est = logS - logSW

print('IS Estimator :', lp_est)

np.max((1,2))
%connect_info<jupyter_output>{
  "ip": "127.0.0.1",
  "stdin_port": 56571,
  "hb_port": 56573,
  "signature_scheme": "hmac-sha256",
  "transport": "tcp",
  "control_port": 56572,
  "iopub_port": 56570,
  "shell_port": 56569,
  "key": "<KEY>"
}

Paste the above JSON into a file, and connect with:
    $> ipython <app> --existing <file>
or, if you are local, you can connect with just:
    $> ipython <app> --existing kernel-0612b3c2-77a2-4d55-b991-57a41d566995.json 
or even just:
    $> ipython <app> --existing 
if this is the most recent IPython session you have started.
{'fraction_non_alphanumeric': 0.17926405369771065, 'fraction_numerical': 0.03224259858564066, 'mean_word_length': 4.000119860961285, 'pattern_counts': {'":': 9, '<': 111, '<?xml version=': 0, '>': 115, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '27', 'dir': 'jupyter-structured-clean-dedup', 'id': '546987', 'license': 'MIT', 'n_tokens_mistral': 17324, 'n_tokens_neox': 16253, 'n_words': 4366, 'path': 'Sampling.ipynb', 'repo_name': 'tugbatugbatugba/data-mining', 'size': '41756'}