4930.sagews

from sage.stats.hmm.distributions import GaussianMixtureDistribution
mean = [2,4,6,8]
prob = [.1,.1,.6,.2]
vari = [.25,.2,.3,.25]
r = len(M)
n = 200
V = GaussianMixtureDistribution((prob[i], mean[i], vari[i]) for i in range(r)).sample(n)
show(V.plot_histogram(bins=100,normalize=False))
V = sorted(list(V))

W = []
for i in range(n):
    s = []
    for j in range(n):
        s.append(float(exp(-(V[i]-V[j])**2/2)))
    W.append(s)
W = matrix(W)

D = []
for i in range(n):
    s = 0
    for j in range(n):
        s += W[i][j]
    D.append(s)
D = diagonal_matrix(D)

L = D - W
e = L.right_eigenvectors()
e.sort(key=lambda x: x[0])
e = [(x,v) for x,vs,_ in e for v in vs]
show(points((i,e[i][0]) for i in range(r*3)))
U = matrix(list(list(x[1]) for x in e[:r])).transpose()
show(matrix_plot(U.transpose()))

DD = diagonal_matrix(list(1/x for x in D.diagonal()))

Lrw = DD * L
e = Lrw.right_eigenvectors()
e.sort(key=lambda x: x[0])
e = [(x,v) for x,vs,_ in e for v in vs]
show(points((i,e[i][0]) for i in range(r*3)))
U = matrix(list(list(x[1]) for x in e[:r])).transpose()
show(matrix_plot(U.transpose()))

DD = diagonal_matrix(list(1/sqrt(x) for x in D.diagonal()))

Lsym = DD * L * DD
e = Lsym.right_eigenvectors()
e.sort(key=lambda x: x[0])
e = [(x,v) for x,vs,_ in e for v in vs]
show(points((i,e[i][0]) for i in range(r*3)))
U = matrix(list(list(x[1]) for x in e[:r])).transpose()
show(matrix_plot(U.transpose()))

T = []
for i in range(n):
    t = 0
    for j in range(r):
        t += U[i][j] * U[i][j]
    t = sqrt(t)
    s = []
    for j in range(r):
        s.append(U[i][j] / t)
    T.append(s)
T = matrix(T)
show(matrix_plot(T.transpose()))