Python之二维正态分布采样置信椭圆绘制
作者:犹有傲霜枝 发布时间:2021-04-08 06:39:09
标签:二维正态分布采样,置信椭圆,Python
二维正态分布采样后,绘制置信椭圆
假设二维正态分布表示为:
下图为两个二维高斯分布采样后的置信椭圆
和
每个二维高斯分布采样100个数据点,图片为:
代码如下
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
def make_ellipses(mean, cov, ax, confidence=5.991, alpha=0.3, color="blue", eigv=False, arrow_color_list=None):
"""
多元正态分布
mean: 均值
cov: 协方差矩阵
ax: 画布的Axes对象
confidence: 置信椭圆置信率 # 置信区间, 95%: 5.991 99%: 9.21 90%: 4.605
alpha: 椭圆透明度
eigv: 是否画特征向量
arrow_color_list: 箭头颜色列表
"""
lambda_, v = np.linalg.eig(cov) # 计算特征值lambda_和特征向量v
# print "lambda: ", lambda_
# print "v: ", v
# print "v[0, 0]: ", v[0, 0]
sqrt_lambda = np.sqrt(np.abs(lambda_)) # 存在负的特征值, 无法开方,取绝对值
s = confidence
width = 2 * np.sqrt(s) * sqrt_lambda[0] # 计算椭圆的两倍长轴
height = 2 * np.sqrt(s) * sqrt_lambda[1] # 计算椭圆的两倍短轴
angle = np.rad2deg(np.arccos(v[0, 0])) # 计算椭圆的旋转角度
ell = mpl.patches.Ellipse(xy=mean, width=width, height=height, angle=angle, color=color) # 绘制椭圆
ax.add_artist(ell)
ell.set_alpha(alpha)
# 是否画出特征向量
if eigv:
# print "type(v): ", type(v)
if arrow_color_list is None:
arrow_color_list = [color for i in range(v.shape[0])]
for i in range(v.shape[0]):
v_i = v[:, i]
scale_variable = np.sqrt(s) * sqrt_lambda[i]
# 绘制箭头
"""
ax.arrow(x, y, dx, dy, # (x, y)为箭头起始坐标,(dx, dy)为偏移量
width, # 箭头尾部线段宽度
length_includes_head, # 长度是否包含箭头
head_width, # 箭头宽度
head_length, # 箭头长度
color, # 箭头颜色
)
"""
ax.arrow(mean[0], mean[1], scale_variable*v_i[0], scale_variable * v_i[1],
width=0.05,
length_includes_head=True,
head_width=0.2,
head_length=0.3,
color=arrow_color_list[i])
# ax.annotate("",
# xy=(mean[0] + lambda_[i] * v_i[0], mean[1] + lambda_[i] * v_i[1]),
# xytext=(mean[0], mean[1]),
# arrowprops=dict(arrowstyle="->", color=arrow_color_list[i]))
# v, w = np.linalg.eigh(cov)
# print "v: ", v
# # angle = np.rad2deg(np.arccos(w))
# u = w[0] / np.linalg.norm(w[0])
# angle = np.arctan2(u[1], u[0])
# angle = 180 * angle / np.pi
# s = 5.991 # 置信区间, 95%: 5.991 99%: 9.21 90%: 4.605
# v = 2.0 * np.sqrt(s) * np.sqrt(v)
# ell = mpl.patches.Ellipse(xy=mean, width=v[0], height=v[1], angle=180 + angle, color="red")
# ell.set_clip_box(ax.bbox)
# ell.set_alpha(0.5)
# ax.add_artist(ell)
def plot_2D_gaussian_sampling(mean, cov, ax, data_num=100, confidence=5.991, color="blue", alpha=0.3, eigv=False):
"""
mean: 均值
cov: 协方差矩阵
ax: Axes对象
confidence: 置信椭圆的置信率
data_num: 散点采样数量
color: 颜色
alpha: 透明度
eigv: 是否画特征向量的箭头
"""
if isinstance(mean, list) and len(mean) > 2:
print "多元正态分布,多于2维"
mean = mean[:2]
cov_temp = []
for i in range(2):
cov_temp.append(cov[i][:2])
cov = cov_temp
elif isinstance(mean, np.ndarray) and mean.shape[0] > 2:
mean = mean[:2]
cov = cov[:2, :2]
data = np.random.multivariate_normal(mean, cov, 100)
x, y = data.T
plt.scatter(x, y, s=10, c=color)
make_ellipses(mean, cov, ax, confidence=confidence, color=color, alpha=alpha, eigv=eigv)
def main():
# plt.figure("Multivariable Gaussian Distribution")
plt.rcParams["figure.figsize"] = (8.0, 8.0)
fig, ax = plt.subplots()
ax.set_xlabel("x")
ax.set_ylabel("y")
print "ax:", ax
mean = [4, 0]
cov = [[1, 0.9],
[0.9, 0.5]]
plot_2D_gaussian_sampling(mean=mean, cov=cov, ax=ax, eigv=True, color="r")
mean1 = [5, 2]
cov1 = [[1, 0],
[0, 1]]
plot_2D_gaussian_sampling(mean=mean1, cov=cov1, ax=ax, eigv=True)
plt.savefig("./get_pickle_data/pic/gaussian_covariance_matrix.png")
plt.show()
if __name__ == "__main__":
main()
来源:https://blog.csdn.net/qq_41009742/article/details/114871413
0
投稿
猜你喜欢
- 现在我们已经知道了,fixtures是一个非常强大的功能。那么有的时候,我们可能会写一个fixture,而这个fixture所有的测试函数都
- 这次让我们来看看一个真实场景吧:银行转账一、举例银行转账假设现在有一个xuewei的账号里面有 100W。然后有多个任务在转账,转入转出都是
- jQuery的选择器是CSS 1-3,XPath的结合物。jQuery提取这二种查询语言最好的部分,融合后创造出了最终的jQuery表达式查
- 本文实例讲述了Python类装饰器。分享给大家供大家参考,具体如下:编写类装饰器类装饰器类似于函数装饰器的概念,但它应用于类,它们可以用于管
- 如下所示:import arcpy... from arcpy import env... env.workspace="C:\\
- 我打算将WebQQ单独出来运行, 一开始直接拷贝了pyxmpp2的mainloop, 但是跑起来问题多多, 所以我又研究了利用Tornado
- 问题描述最近~ 发现对series里的元素操作挺复杂的,用for loop + Series.iloc[i]会发生卡死的状况,那么,lamb
- Oracle9i中提供强大的迁移功能,可以从多种数据库向Oracle迁移数据。Oracle新发行的迁移工具提供了从Access2
- 本文实例总结了PHP session会话操作技巧。分享给大家供大家参考,具体如下:会话技术session将会话数据存储与服务器端,同时使会话
- 一、装饰器decorator decorator设计模式允许动态地对现有的对象或函数包装以至于修改现有的职责和行为,简单地讲用来动态地扩展现
- 每次抽取后都重新洗牌。计算10000次随机抽取可得到同花的几率。我做的比较复杂,分别累计了四种花色分别出现了几次import randoml
- asp fso type属性取得文件类型代码是用来返回类型指定的文件或文件夹。语法FileObject.Type FolderObject.
- 插入一条记录后,如何得到最新的自动增加ID?我们要用到SQL Server的@@IDENTITY。它能够记录下系统最近使用的一个IDENTI
- 一、请求扩展1.before_request作用: 类比django中间件中的process_request,在请求到来执行路由函数之前先执
- 1、删除Oracal在注册表中的主项:regedit.exe->LocalMachine->Software->Oracl
- 本文实例为大家分享了python批量文件重命名的具体代码,供大家参考,具体内容如下问题描述最近遇到朋友求助,如何将大量文件名前面的某些字符删
- 本文实例为大家分享了使用python的twisted框架编写一个简单的聊天室具体代码,供大家参考,具体内容如下下面是基本架构代码:# -*-
- 1.简述:xlrd是读取excel表格数据;支持 xlsx和xls 格式的excel表格;三方模块安装方式:pip3 install xlr
- 一:命名空间里的namespace关键字和__NAMESPACE__常量的运用PHP支持两种抽象的访问当前命名空间内部元素的方法,__NAM
- 本文实例为大家分享了python实现通讯录系统的具体代码,供大家参考,具体内容如下一、首先,给大家看一下这个系统的主菜单效果如下:二、步骤分