# Run this code! But there is no need to pay much attention to this cell at the first pass through the notebook

#%matplotlib inline 
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as sps


def rand_gauss(n, mu, cov):
    """
    Sample n data points from multivariate Gaussian distribution with mean mu and covariance cov
    """
    return np.atleast_2d(sps.multivariate_normal.rvs(mu, cov, n))

def logpdf_gauss(x, mu, cov):
    """
    Evaluation of the log probability density function for multivariate Gaussian with mean mu and covariance cov
    """
    return sps.multivariate_normal.logpdf(x, mu, cov)
   
def gellipse(mu, cov, n=100, *args, **kwargs):
    """
    Contour plot of 2D Multivariate Gaussian distribution.

    gellipse(mu, cov, n) plots ellipse given by mean vector MU and
    covariance matrix COV. Ellipse is plotted using N (default is 100)
    points. Additional parameters can specify various line types and
    properties. See description of matplotlib.pyplot.plot for more details.
    """
    if mu.shape != (2,) or cov.shape != (2, 2):
        raise RuntimeError('mu must be a two element vector and cov must be 2 x 2 matrix')

    d, v = np.linalg.eigh(4 * cov)
    d = np.diag(d)
    t = np.linspace(0, 2 * np.pi, n)
    x = v @ np.sign(d) @ np.sqrt(np.abs(d)) @ np.array([np.cos(t), np.sin(t)]) + mu[:,np.newaxis]
    return plt.plot(x[0], x[1], *args, **kwargs)

def probit(a):
    from scipy.special import erfinv
    return np.sqrt(2.0) * erfinv(2.0 * a - 1.0)

def plot_det(tar, non, label="",
             axis = [0.2, 40, 0.2, 80],
             xticks = [0.2, 0.5, 1, 2, 5, 10, 20, 35, 50, 65, 80],
             yticks = [0.2, 0.5, 1, 2, 5, 10, 20, 35, 50, 65, 80],
             **kwargs):
        """
        plots DET curve 
        """
        tar = np.array(tar)
        non = np.array(non)
        ntrue=len(tar)
        nfalse=len(non)
        ntotal=ntrue+nfalse

        Pmiss=np.zeros(ntotal+1,np.float32) # 1 more for the boundaries
        Pfa=np.zeros_like(Pmiss)

        scores=np.zeros((ntotal,2),np.float32)
        scores[0:nfalse,0]=non
        scores[0:nfalse,1]=0
        scores[nfalse:ntotal,0]=tar
        scores[nfalse:ntotal,1]=1
        scores=scores[scores[:,0].argsort(),]

        sumtrue=np.cumsum(scores[:,1])
        sumfalse=nfalse - (np.arange(1,ntotal+1)-sumtrue)

        Pmiss[0]=float(ntrue-ntrue) / ntrue
        Pfa[0]=float(nfalse) / nfalse
        Pmiss[1:]=(sumtrue+ntrue-ntrue) / ntrue
        Pfa[1:]=sumfalse / nfalse
        
        idxeer=np.argmin(np.abs(Pfa-Pmiss))
        EER = 0.5*(Pfa[idxeer]+Pmiss[idxeer])*100

        plt.plot(probit(Pfa), probit(Pmiss), label=label + ' EER=%.2f%%' % EER, **kwargs)
        plt.xticks(probit(np.array(xticks)/100), xticks)
        plt.yticks(probit(np.array(yticks)/100), yticks)
        plt.axis(probit(np.array(axis)/100))

        plt.xlabel("FA [%]", fontsize = 12)
        plt.ylabel("Miss [%]", fontsize = 12)
        plt.grid(True)
        plt.legend(loc='upper left', bbox_to_anchor=(1, 1))

#Do not edit this code!
mu_gt =  np.array([2, 1]) #ground truth global mean
Sigma_wc_gt = np.array([[1, -0.8], #ground truth within class covariance
                     [-0.8, 1]])

Sigma_ac_gt = np.array([[10, 2], #ground truth across class covariance
                     [2, 1]])

def plotPLDA(mu, Sigma_ac, Sigma_wc, marker_mu, line_style_Sigmas, name):
    """
    plot of a PLDA model that shows: 
    the global mean mu in green,
    the across-class covariance Sigma_ac in blue, 
    and the within-class covariance Sigma_wc centered around [0,0] in red.
    Possible markers for mu are: '.','o','x' (see https://matplotlib.org/stable/api/markers_api.html)
    Possible line_styles for the countours are: '--', ':', '-'
    name refers to the name of the PLDA model to be displayed in the legend
    
    """
    assert marker_mu in {"." , "o" , "x"}
    assert line_style_Sigmas in {"--" , ":" , "-"}
    
    plt.plot(mu[0], mu[1], 'g'+marker_mu, ms=10, label='Mean '+name) 
    gellipse(mu, Sigma_ac, 100, 'b'+line_style_Sigmas, lw=2, label='Sigma_ac '+name)
    gellipse(np.array([0,0]), Sigma_wc, 100, 'r'+line_style_Sigmas, lw=2, label='Sigma_wc '+name) #for mere visualization purposes, we center the within-class covariance on the origin (0,0)
    plt.axis('equal')
    plt.legend(loc='upper left', bbox_to_anchor=(1, 1))

plotPLDA(mu_gt,Sigma_ac_gt,Sigma_wc_gt,'.','--','GT')

#Do not edit this code!
def sample_from_PLDA(mu, Sigma_ac, Sigma_wc, S, min_samples_per_speaker, max_samples_per_speaker):
    """
    Generate samples from a PLDA distribution, where each speaker can have a different number of samples:
    mu: global mean
    Sigma_ac: across-class covariance matrix
    Sigma_wc: within-class covariance matrix
    S: Number of speakers to be sampled
    min_samples_per_speaker: minimum number of observations to be sampled for a speaker
    max_samples_per_speaker: maximum number of observations to be sampled for a speaker
    
    Returns:
    Z: two dimensional array where each row is a sampled speaker mean
    X: list of two dimensional arrays, where rows of each array are the individual observations per speaker
    """
    N = np.random.randint(min_samples_per_speaker, max_samples_per_speaker+1, S) # Number of observations per speaker
    Z =  rand_gauss(S, mu, Sigma_ac) # speaker means
    X = [] # Collection of all the X_s

    # For each speaker
    for ns, z in zip(N, Z):
        X_s = rand_gauss(ns, z, Sigma_wc)
        X.append(X_s)
    return X, Z   #Recall, that we do not return N, which can be easily derived from X

#Sampled training data for 10 speakers with only one or two samples per speaker:
X, Z = sample_from_PLDA(mu_gt, Sigma_ac_gt, Sigma_wc_gt, S=10, min_samples_per_speaker=1, max_samples_per_speaker=2)

gellipse(mu_gt, Sigma_ac_gt, 100, 'b', lw=2)
for X_s, z in zip(X, Z):
    p = plt.plot(X_s[:,0], X_s[:,1], '.', ms=2)
    c = p[0].get_color()
    plt.plot(z[0], z[1], '.', c=c, ms=10)
    gellipse(z, Sigma_wc_gt, 100, c=c)

def simple_PLDA_estimate(X):
    """
    Estimates the PLDA parameters using the simple maximum-likelihood approach:
    - X: whole set of training embeddings (for all speakers) as a list of two dimensional arrays, 
         where rows of each array are the individual observations per speaker
    
    Returns:
    - mu: global mean
    - Sigma_ac: across-class covariance matrix
    - Sigma_wc: within-class covariance matrix
    """
    
    # your code goes here
    
    return mu, Sigma_wc, Sigma_ac

#Make use of the following variable names
#mu_simple1, Sigma_wc_simple1, Sigma_ac_simple1 = simple_PLDA_estimate(X)



#code for plots goes here

def p_z_given_X(X, mu, Sigma_ac, Sigma_wc):
    """
    For each speaker 's' in the input data X, it calculates the parameters of the posterior distribution 
    p(z_s|X_s)= N(z_s|mu_s,Sigma_s). For each speaker such parameters are mean (mu_s) 
    and covariance matrix (Sigma_s). The function returns the list of such means (one per speaker) 
    and the list of such covariance matrices (one per speaker).
    
    Inputs:
    - X: whole set of embeddings (for all speakers) 
    - mu: global mean
    - Sigma_ac: across-class covariance matrix
    - Sigma_wc: within-class covariance matrix   
    
    Returns:
    - mu_s: list of per-speaker posterior distribution means 
    - Sigma_s: list of per-speaker posterior distribution means covariance matrices
    """
    
    # Your code goes here
    
    return mu_s, Sigma_s
    
    
    
def EM_PLDA_estimate(X, mu_init, Sigma_ac_init, Sigma_wc_init, niters):
    """
    Runs expectation-maximization algorithm for PLDA parameter estimation, where the input parameters are:
    - X: whole set of embeddings (for all speakers)
    - mu_init: initialization for the global mean
    - Sigma_ac_init: initialization for the across-class covariance
    - Sigma_wc_init: initialization for the within-class covariance
    - niters: number of EM iterations 
    
    Returns the EM estimates of:
    - mu: global mean
    - Sigma_ac: across-class covariance matrix
    - Sigma_wc: within-class covariance matrix
    """

    # Your code goes here

    return mu, Sigma_wc, Sigma_ac


# Use the following variable names for storing the PLDA model trained with EM in this task:
# mu_EM1, Sigma_wc_EM1, Sigma_ac_EM1 = EM_PLDA_estimate()




#code for plots goes here

#X, Z = #sample_from_PLDA() 1000 speakers, with 1 to 2 samples per speaker

# Use the following variable names for storing the PLDA model trained with simple ML in this task:
# mu_simple2, Sigma_wc_simple2, Sigma_ac_simple2

# Use the following variable names for storing the PLDA model trained with EM in this task:
# mu_EM2, Sigma_wc_EM2, Sigma_ac_EM2

#code for plots goes here

#X, Z = #sample_from_PLDA() 1000 speakers, with 1 to 100 samples per speaker


# Use the following variable names for storing the PLDA model trained with simple ML in this task:
# mu_simple3, Sigma_wc_simple3, Sigma_ac_simple3

# Use the following variable names for storing the PLDA model trained with EM in this task:
# mu_EM3, Sigma_wc_EM3, Sigma_ac_EM3


#code for plots goes here

def LLRs(x_test, X_enrol, mu, Sigma_ac, Sigma_wc):
    """Given a PLDA model, evaluates log-likelihood ratio scores for a set of N verification trials.
    Each trial consists of one test embedding and set of enrollment embeddings.
    - x_test: two dimensional array with N rows, where each row contains a test embedding from one trial 
    - X_enrol: list of N two dimensional numpy arrays. Each element in the list contains a set 
        of enrollment embeddings from one trial. Embeddings from each enrollment set are stored
        in rows of a two dimensional numpy array. 
    - mu: PLDA global mean
    - Sigma_ac: PLDA across-class covariance matrix
    - Sigma_wc: PLDA within-class covariance matrix  
    
    Returns: 
    - one dimensional numpy array containing N log-likelihood ratios.
    """
    #Your code goes here

# Your code for task 9 goes here

#Your code for the plot goes here

#recall, that the contours can be plotted using: 
#gellipse(mu, Sigma, 100, ':', c='b', label='Sigma')
# where ':' stands for for dotted line and 'b' for blue

# Enrollment embedding and trial set generation. Do not edit this field!

# Note, that all generated embeddings are stored in a form compatible with the definition 
# of the LLR given above: the test embeddings are stored in rows of two dimensional arrays 
# and the enrollment embeddings (even the single-enrollment ones) are stored in a list of 
# two dimensional numpy arrays, where each element in the list contains a set of enrollment 
# embeddings.

# we sample two embeddings for 1000 speakers:
X_tar = sample_from_PLDA(mu_gt, Sigma_ac_gt, Sigma_wc_gt, S=1000, min_samples_per_speaker=2, max_samples_per_speaker=2)[0] 

# for enrollment, we use the first embedding per-speaker:
x_enroll     = [xs[:1] for xs in X_tar] 

# the second embedding per-speaker is used as the target test embedding:
x_test_target = np.array([xs[-1] for xs in X_tar]) 

# we generate a non-target embeddings from different random speakers:
x_test_nontar = np.vstack(sample_from_PLDA(mu_gt, Sigma_ac_gt, Sigma_wc_gt, 1000, 1, 1)[0])

# Your code goes here

# Your code goes here

BAYa class Assignment 2025¶

PLDA generative process¶

Joint probability:¶

Handcrafting the PLDA model¶

Sampling training data¶

Simple maximum-likelihood estimate of parameters¶

Task 1¶

PLDA Expectation Maximization training¶

Summary of the EM algorithm¶

E-step¶

Task 2:¶

M-step¶

Task 3:¶

Task 4¶

Task 5¶

Task 6¶

Task 7¶

PLDA Scoring¶

Task 8¶

Task 9¶

Task 10¶

Speaker verification experiments¶

Task 11¶

Task 12¶

Your code goes here¶