59 lines
2 KiB
Python
59 lines
2 KiB
Python
|
#!/usr/bin/env python3
|
||
|
|
import numpy as np
|
||
|
|
from sipy import integrate
|
||
|
|
from scipy import norm
|
||
|
|
|
||
|
|
M = 8
|
||
|
|
X = np.random.normal(0,1,1000)
|
||
|
|
|
||
|
|
def ddp(x):
|
||
|
|
mean = 0,
|
||
|
|
sigma = 1
|
||
|
|
return norm.pdf(x,mean,sigma)
|
||
|
|
|
||
|
|
def init_thres_vec(M,X):
|
||
|
|
step = (np.max(X)-np.min(X))/M
|
||
|
|
thres_intervals = np.array([])
|
||
|
|
mid = np.mean(X)
|
||
|
|
for i in range(int(M/2)):
|
||
|
|
thres_intervals = np.append(thres_vec,mid+(i+1)*step)
|
||
|
|
thres_intervals = np.insert(thtres_vec,0,mid-(1+1)*step)
|
||
|
|
return thres_intervals
|
||
|
|
|
||
|
|
def quant(x,thres,intervals):
|
||
|
|
thres= np.append(thres, np.inf)
|
||
|
|
thres= np.insert(thres, 0, -np.inf)
|
||
|
|
x_hat_q = np.zeros(np.shape(x))
|
||
|
|
for i in range(len(thres)-1):
|
||
|
|
if i == 0:
|
||
|
|
x_hat_q = np.where(np.logical_and(x > thres[i], x <= thres[i+1]),
|
||
|
|
np.full(np.size(x_hat_q), intervals[i]), x_hat_q)
|
||
|
|
elif i == range(len(thres))[-1]-1:
|
||
|
|
x_hat_q = np.where(np.logical_and(x > thres[i], x <= thres[i+1]),
|
||
|
|
np.full(np.size(x_hat_q), intervals[i]), x_hat_q)
|
||
|
|
else:
|
||
|
|
x_hat_q = np.where(np.logical_and(x > thres[i], x < thres[i+1]),
|
||
|
|
np.full(np.size(x_hat_q), intervals[i]), x_hat_q)
|
||
|
|
return x_hat_q
|
||
|
|
|
||
|
|
|
||
|
|
def LlyodMax(X,intervals, max_iter=1000,eps=1e-5):
|
||
|
|
err_min = np.inf
|
||
|
|
for i in range(max_iter):
|
||
|
|
for j in range(len(x_hat_q)):
|
||
|
|
centroids[i] = integrate.quad(lambda x : x*ddp(x),
|
||
|
|
intervals[j],intervals[j+1])[0]/
|
||
|
|
integrate.quad(lambda x : ddp(x),
|
||
|
|
intervals[j],intervals[j+1])[0]
|
||
|
|
intervals = 0.5*(centroids[1:]+centroids[:-1])
|
||
|
|
x_hat = quant(X,centroids,intervals)
|
||
|
|
err = np.linalg.norm(X-x_hat)
|
||
|
|
if err < err_min:
|
||
|
|
err_min =err
|
||
|
|
intervals_min = intervals
|
||
|
|
centroids_min = centroids
|
||
|
|
if err_min< 1e-5:
|
||
|
|
break
|
||
|
|
best_x_hat = quant(X,centroids_min,intervals_min)
|
||
|
|
return best_x_hat
|