Lloyd Max Quantizer for optimal quantization of a random variable demonstrating a Gaussian PDF

lloyd-max.py

"""
Lloyd Max Quantizer
This program implements a midrise lloyd-max quantizer, for optimal
quantization of a random variable demonstrating a Gaussian PDF.
The clustering mechanism used here equivalent to k-Means clustering.

Author: Dhruv Ilesh Shah
EE308 - Communication Systems, Autumn 2017
Indian Institute of Technology, Bombay
"""

import numpy as np
from scipy.stats import norm
from scipy import integrate


def cluster(levels, mean, variance, max_iter=100):
    """
    This routine clusters the positive portion of the Gaussian PDF of given mean
    and variance into required number of levels using an iterative update routine.
    The output is an array of converged cluster centroids.
    """
    max_int = 10 * variance  # p-Value ~ 1 (< 1e-12)
    intervals = np.linspace(0., max_int, levels + 1)
    centroids = np.zeros(levels)
    intervals[levels] = max_int
    thresh = 1e-5

    # while (update > thresh):
    for i in range(max_iter):
        intervals_prev = np.copy(intervals)
        for j in range(levels):
            centroids[j] = integrate.quad(lambda x: x * norm.pdf(x, mean, variance**0.5), intervals[j], intervals[
                j + 1])[0] / integrate.quad(lambda x: norm.pdf(x, mean, variance**0.5), intervals[j], intervals[j + 1])[0]

        for j in range(levels - 1):
            intervals[j + 1] = (centroids[j + 1] + centroids[j]) / 2.

        if (np.linalg.norm(intervals_prev - intervals) < thresh):
            break

    return intervals, centroids


mean = 0
variance = 4
num_levels = 6

[intervals, clusters] = cluster(num_levels / 2, mean, variance)
clusters = np.append(clusters, -1*clusters)
clusters = np.sort(clusters)

intervals = np.append(intervals, -1*intervals)
intervals = np.unique(intervals)

print clusters, intervals