如何在python中实现EM-GMM？
我已经实现了
算法
对于GMM使用这个
后[GMMs]与最大似然优化
努比](https://towardsdatascience.com/how-to-code-gaussian-mixture-models-
from-scratch-in-python-9e7975df5252）未成功，如下所示：
import numpy as np
def PDF(data, means, variances):
return 1/(np.sqrt(2 * np.pi * variances) + eps) * np.exp(-1/2 * (np.square(data - means) / (variances + eps)))
def EM_GMM(data, k, iterations):
weights = np.ones((k, 1)) / k # shape=(k, 1)
means = np.random.choice(data, k)[:, np.newaxis] # shape=(k, 1)
variances = np.random.random_sample(size=k)[:, np.newaxis] # shape=(k, 1)
data = np.repeat(data[np.newaxis, :], k, 0) # shape=(k, n)
for step in range(iterations):
# Expectation step
likelihood = PDF(data, means, np.sqrt(variances)) # shape=(k, n)
# Maximization step
b = likelihood * weights # shape=(k, n)
b /= np.sum(b, axis=1)[:, np.newaxis] + eps
# updage means, variances, and weights
means = np.sum(b * data, axis=1)[:, np.newaxis] / (np.sum(b, axis=1)[:, np.newaxis] + eps)
variances = np.sum(b * np.square(data - means), axis=1)[:, np.newaxis] / (np.sum(b, axis=1)[:, np.newaxis] + eps)
weights = np.mean(b, axis=1)[:, np.newaxis]
return means, variances
when I run the algorithm on a 1-D time-series dataset, for k equal to 3, it
returns an output like the following:
array([[0.00000000e+000, 0.00000000e+000, 0.00000000e+000,
3.05053810e-003, 2.36989898e-025, 2.36989898e-025,
1.32797395e-136, 6.91134950e-031, 5.47347807e-001,
1.44637007e+000, 1.44637007e+000, 1.44637007e+000,
1.44637007e+000, 1.44637007e+000, 1.44637007e+000,
1.44637007e+000, 1.44637007e+000, 1.44637007e+000,
1.44637007e+000, 1.44637007e+000, 1.44637007e+000,
1.44637007e+000, 2.25849208e-064, 0.00000000e+000,
1.61228562e-303, 0.00000000e+000, 0.00000000e+000,
0.00000000e+000, 0.00000000e+000, 3.94387272e-242,
1.13078186e+000, 2.53108878e-001, 5.33548114e-001,
9.14920432e-001, 2.07015697e-013, 4.45250680e-038,
1.43000602e+000, 1.28781615e+000, 1.44821615e+000,
1.18186109e+000, 3.21610659e-002, 3.21610659e-002,
3.21610659e-002, 3.21610659e-002, 3.21610659e-002,
2.47382844e-039, 0.00000000e+000, 2.09150855e-200,
0.00000000e+000, 0.00000000e+000],
[5.93203066e-002, 1.01647068e+000, 5.99299162e-001,
0.00000000e+000, 0.00000000e+000, 0.00000000e+000,
0.00000000e+000, 0.00000000e+000, 0.00000000e+000,
0.00000000e+000, 0.00000000e+000, 0.00000000e+000,
0.00000000e+000, 0.00000000e+000, 0.00000000e+000,
0.00000000e+000, 0.00000000e+000, 0.00000000e+000,
0.00000000e+000, 0.00000000e+000, 0.00000000e+000,
0.00000000e+000, 0.00000000e+000, 2.14690238e-010,
2.49337135e-191, 5.10499986e-001, 9.32658804e-001,
1.21148135e+000, 1.13315278e+000, 2.50324069e-237,
0.00000000e+000, 0.00000000e+000, 0.00000000e+000,
0.00000000e+000, 0.00000000e+000, 0.00000000e+000,
0.00000000e+000, 0.00000000e+000, 0.00000000e+000,
0.00000000e+000, 0.00000000e+000, 0.00000000e+000,
0.00000000e+000, 0.00000000e+000, 0.00000000e+000,
0.00000000e+000, 1.73966953e-125, 2.53559290e-275,
1.42960975e-065, 7.57552338e-001],
[0.00000000e+000, 0.00000000e+000, 0.00000000e+000,
3.05053810e-003, 2.36989898e-025, 2.36989898e-025,
1.32797395e-136, 6.91134950e-031, 5.47347807e-001,
1.44637007e+000, 1.44637007e+000, 1.44637007e+000,
1.44637007e+000, 1.44637007e+000, 1.44637007e+000,
1.44637007e+000, 1.44637007e+000, 1.44637007e+000,
1.44637007e+000, 1.44637007e+000, 1.44637007e+000,
1.44637007e+000, 2.25849208e-064, 0.00000000e+000,
1.61228562e-303, 0.00000000e+000, 0.00000000e+000,
0.00000000e+000, 0.00000000e+000, 3.94387272e-242,
1.13078186e+000, 2.53108878e-001, 5.33548114e-001,
9.14920432e-001, 2.07015697e-013, 4.45250680e-038,
1.43000602e+000, 1.28781615e+000, 1.44821615e+000,
1.18186109e+000, 3.21610659e-002, 3.21610659e-002,
3.21610659e-002, 3.21610659e-002, 3.21610659e-002,
2.47382844e-039, 0.00000000e+000, 2.09150855e-200,
0.00000000e+000, 0.00000000e+000]])
我认为这是错误的，因为输出是两个向量哪一个
其中一个表示“平均值”，另一个表示“方差”`
价值观。使我对执行产生怀疑的模糊点是
返回“0.00000000e+000”，以显示和显示大多数输出
它不需要真正地将这些输出可视化。顺便说一下，输入数据是
时间序列数据。我检查了所有东西，追踪了好几次，但是
没有虫子出现。
Here are my input data:
[25.31 , 24.31 , 24.12 , 43.46 , 41.48666667,
41.48666667, 37.54 , 41.175 , 44.81 , 44.44571429,
44.44571429, 44.44571429, 44.44571429, 44.44571429, 44.44571429,
44.44571429, 44.44571429, 44.44571429, 44.44571429, 44.44571429,
44.44571429, 44.44571429, 39.71 , 26.69 , 34.15 ,
24.94 , 24.75 , 24.56 , 24.38 , 35.25 ,
44.62 , 44.94 , 44.815 , 44.69 , 42.31 ,
40.81 , 44.38 , 44.56 , 44.44 , 44.25 ,
43.66666667, 43.66666667, 43.66666667, 43.66666667, 43.66666667,
40.75 , 32.31 , 36.08 , 30.135 , 24.19 ]
我想知道是否有一个优雅的方式来实现它通过“numpy”或SciKit学习。任何帮助都将不胜感激。
问题答案：
正如我在评论中提到的，我看到的关键点是“手段”`初始化。遵循[sk]的默认实现
混合物,
我没有随机初始化，而是切换到KMeans。
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
plt.style.use('seaborn')
eps=1e-8
def PDF(data, means, variances):
return 1/(np.sqrt(2 * np.pi * variances) + eps) * np.exp(-1/2 * (np.square(data - means) / (variances + eps)))
def EM_GMM(data, k=3, iterations=100, init_strategy='kmeans'):
weights = np.ones((k, 1)) / k # shape=(k, 1)
if init_strategy=='kmeans':
from sklearn.cluster import KMeans
km = KMeans(k).fit(data[:, None])
means = km.cluster_centers_ # shape=(k, 1)
else: # init_strategy=='random'
means = np.random.choice(data, k)[:, np.newaxis] # shape=(k, 1)
variances = np.random.random_sample(size=k)[:, np.newaxis] # shape=(k, 1)
data = np.repeat(data[np.newaxis, :], k, 0) # shape=(k, n)
for step in range(iterations):
# Expectation step
likelihood = PDF(data, means, np.sqrt(variances)) # shape=(k, n)
# Maximization step
b = likelihood * weights # shape=(k, n)
b /= np.sum(b, axis=1)[:, np.newaxis] + eps
# updage means, variances, and weights
means = np.sum(b * data, axis=1)[:, np.newaxis] / (np.sum(b, axis=1)[:, np.newaxis] + eps)
variances = np.sum(b * np.square(data - means), axis=1)[:, np.newaxis] / (np.sum(b, axis=1)[:, np.newaxis] + eps)
weights = np.mean(b, axis=1)[:, np.newaxis]
return means, variances
This seems to yield the desired output much more consistently:
s = np.array([25.31 , 24.31 , 24.12 , 43.46 , 41.48666667,
41.48666667, 37.54 , 41.175 , 44.81 , 44.44571429,
44.44571429, 44.44571429, 44.44571429, 44.44571429, 44.44571429,
44.44571429, 44.44571429, 44.44571429, 44.44571429, 44.44571429,
44.44571429, 44.44571429, 39.71 , 26.69 , 34.15 ,
24.94 , 24.75 , 24.56 , 24.38 , 35.25 ,
44.62 , 44.94 , 44.815 , 44.69 , 42.31 ,
40.81 , 44.38 , 44.56 , 44.44 , 44.25 ,
43.66666667, 43.66666667, 43.66666667, 43.66666667, 43.66666667,
40.75 , 32.31 , 36.08 , 30.135 , 24.19 ])
k=3
n_iter=100
means, variances = EM_GMM(s, k, n_iter)
print(means,variances)
[[44.42596231]
[24.509301 ]
[35.4137508 ]]
[[0.07568723]
[0.10583743]
[0.52125856]]
# Plotting the results
colors = ['green', 'red', 'blue', 'yellow']
bins = np.linspace(np.min(s)-2, np.max(s)+2, 100)
plt.figure(figsize=(10,7))
plt.xlabel('$x$')
plt.ylabel('pdf')
sns.scatterplot(s, [0.05] * len(s), color='navy', s=40, marker=2, label='Series data')
for i, (m, v) in enumerate(zip(means, variances)):
sns.lineplot(bins, PDF(bins, m, v), color=colors[i], label=f'Cluster {i+1}')
plt.legend()
plt.plot()
最后我们可以看到，纯随机初始化会产生不同的结果；让我们看看结果“means”：
for _ in range(5):
print(EM_GMM(s, k, n_iter, init_strategy='random')[0], '\n')
[[44.42596231]
[44.42596231]
[44.42596231]]
[[44.42596231]
[24.509301 ]
[30.1349997 ]]
[[44.42596231]
[35.4137508 ]
[44.42596231]]
[[44.42596231]
[30.1349997 ]
[44.42596231]]
[[44.42596231]
[44.42596231]
[44.42596231]]
在某些情况下，结果是如何不同的呢是常量，意味着选择了3个相似的值而没有
迭代时会有很大的变化。在EM\GMM中添加一些打印语句我会澄清的。
-
本文向大家介绍python em算法的实现，包括了python em算法的实现的使用技巧和注意事项，需要的朋友参考一下 以上就是python em算法的实现的详细内容，更多关于python em算法的资料请关注呐喊教程其它相关文章！
-
本文向大家介绍Python实现EM算法实例代码，包括了Python实现EM算法实例代码的使用技巧和注意事项，需要的朋友参考一下 EM算法实例 通过实例可以快速了解EM算法的基本思想，具体推导请点文末链接。图a是让我们预热的，图b是EM算法的实例。 这是一个抛硬币的例子，H表示正面向上，T表示反面向上，参数θ表示正面朝上的概率。硬币有两个，A和B，硬币是有偏的。本次实验总共做了5组，每组随机选一个硬
-
问题内容： 我以前的编程中，代码段仅用于调试目的（记录命令等）。通过使用预处理程序指令，可以完全禁用这些语句以进行生产，如下所示： 做类似的事情的最好方法是什么？ 问题答案： 如果只想禁用日志记录方法，请使用该模块。如果日志级别设置为排除调试语句，那么它将非常接近无操作（它仅检查日志级别并返回而不插入日志字符串）。 如果要在特定条件下以字节码编译时实际删除代码块，则唯一的选择是相当神秘的全局变量。
-
问题内容： 如何实现与C＃代码等效的Python？ 这是一个好主意吗？？请在您的答案中举例说明。 问题答案： 正如其他人在这里提到的： 在Python中不需要接口。这是因为Python具有适当的多重继承，还具有鸭式输入法，这意味着 必须 在Java中具有接口的地方，而不必在Python中具有接口。 也就是说，接口还有多种用途。其中一些被Python 2.6中引入的Pythons抽象基类覆盖。如果您
-
问题内容： 我想知道Python 3中的新super是如何实现的。 我做了一个小例子之后，这个问题就浮现在脑海中，我得到了一个奇怪的错误。我正在使用Pyutilib组件体系结构（PCA），并且已经制作了自定义元类来驱动另一个类的创建： 我收到以下错误： 我wonderign究竟超（一样），它引发的错误，而所有的，并存在。另外的“老方法”-正在工作。 问题答案： 如何执行？这是python3.3的代
-
从Udacity的深度学习类来看，y_i的软最大值只是指数除以整个Y向量的指数之和： 其中< code>S(y_i)是< code>y_i的softmax函数，而< code>e是指数，而< code>j是输入向量y中的列数 我尝试了以下方法： 返回： 但建议的解决方案是： 它产生与第一个实现相同的输出，即使第一个实现显式地取每列和max的差值，然后除以总和。 有人能用数学说明为什么吗？一个正确，