当前位置:首页 > 网络技术 > 多元正态分布的定义(多元统计分析经典例题)(多元正态分布的例题答案)

多元正态分布的定义(多元统计分析经典例题)(多元正态分布的例题答案)

longge2022-03-24 01:50:09网络技术91

关于本文多元正态分布推断(Inference for a Multivariate Normal Population)理论参见《Applied Multivariate Statistical Analysis and Related Topics with R》第五章内容,书中提供了相关示例的R代码,对于偏爱Python的我,希望通过Python得到同样的结果。

数据集的获取网址:
www.stat.ubc.ca/~lang/text.

示例用到的数据集分别为:class.dat2, consum2000.txt, consum2010.txt.

进行Python编程分析前,先把数据集通过R软件转换下格式,虽然Python也可以读取txt文件,但我更喜欢读取csv格式,所以通过以下代码,将数据集转换为CSV格式并保存本地。

data = read.table('class.dat2', header = T)
write.csv(data, 'class.csv')

示例1

需要用到class数据集,这个示例可以简单概括为期中考试前有两次测试quiz1和quiz2,期中考试后有两次测试quiz3和quiz4,比较quiz1和quiz2之间学生成绩有无进步,以及quiz3和quiz4之间学生成绩有无进步。令μ1 = mean(quiz1 – quiz2), μ2 = mean(quiz3 – quiz4)。

进行多元正态分布推断,编程思路为:

1.导入数据 -> 2.求解样本均值和样本协方差阵 -> 3.计算Hotelling’s T 统计量 -> 4.计算p值,根据p值结合实例分析结果。

代码实现如下:

# 引入第三方库
import pandas as pd
import numpy as np
from scipy import stats


# 读取数据
data = pd.read_csv("class.csv")
df = pd.DataFrame(data)
# 构建矩阵
y = np.c_[df.quiz1-df.quiz2, df.quiz3-df.quiz4]
print(y)
# 计算样本数量n
n = np.shape(y)[0]
# 计算变量数目p
p = np.shape(y)[1]


# 计算样本均值
y = pd.DataFrame(y)
y_bar = y.mean()
print(y_bar)


# 计算样本协方差
S_y = y.cov()
print(S_y)


# 计算 Hotelling's T statistic
T_sq = n * np.dot(np.dot(y_bar.T, np.linalg.inv(S_y)), y_bar)
T_sq2 = ((n - p)/(p * (n - 1))) * T_sq
print('T_sq2:', T_sq2)


# 计算p值
p_value = 1 - stats.f.cdf(T_sq2, p, n-p)
print('p_value:', p_value)

输出结果:

p_value: 0.05442091231270707

由p值可以看出quiz1和quiz2之间、quiz3和quiz4之间存在一些差异,但是这些差异在5%水平不是统计显著的。

示例2

需要用到consum2000, consum2010两个数据集,这个示例可以简单概括为比较2000年和2010年在食品(Food)、衣物(Cloth)、居民数(Resid)、交通(TranC)以及教育(Educ)消费结构有无变化。令

μ1 = mean(Food.2010 – Food.2000);

μ2 = mean(Cloth.2010 – Cloth.2000);

μ3 = mean(Resid.2010 – Resid.2000);

μ4 = mean(TranC.2010 – TranC.2000);

μ5 = mean(Educ.2010 – Educ.2000).

代码实现:

# 引入第三方库
import numpy as np
import pandas as pd
from scipy import stats


# 导入数据
consum00 = pd.read_csv("consum2000.csv")
consum10 = pd.read_csv("consum2010.csv")


# 计算2010年支出份额
data10 = consum10.iloc[:, 1:9]
sum10 = data10.sum(axis=1)
X = data10.div(sum10, axis='rows')
print(X)


# 计算2000年支出份额
data00 = consum00.iloc[:, 1:9]
sum00 = data00.sum(axis=1)
Y = data00.div(sum00, axis='rows')
print(Y)


# 求X与Y之差
XY_d = np.c_[X.iloc[:, 0:3]-Y.iloc[:, 0:3], X.iloc[:, 5:7]-Y.iloc[:, 5:7]]
XY_d = pd.DataFrame(XY_d, columns=('Food', 'Cloth', 'Resid', 'TranC', 'Educ'))
# 计算样本均值
d_mean = XY_d.mean()
print(d_mean)
# 计算样本协方差阵
d_S = XY_d.cov()


# 计算样本大小
n = np.shape(XY_d)[0]
# 计算变量数
p = np.shape(XY_d)[1]


# 计算 Hotelling's T 统计量
T2 = n * np.dot(np.dot(d_mean.T, np.linalg.inv(d_S)), d_mean)
Tstar2 = ((n-p)/(p*(n-1)))*T2


# 计算p值
p_value = 1 - stats.f.cdf(Tstar2, p, n-p)
print('p_value:', p_value)

输出结果:

p_value: 7.460698725481052e-14

可见p值近似为0,拒绝原假设,说明2000年与2010年的消费结构发生了明显的变化。

免责声明
本站部分资源来源于互联网 如有侵权 请联系站长删除
龙哥网是优质的互联网科技创业资源_行业项目分享_网络知识引流变现方法的平台为广大网友提供学习互联网相关知识_内容变现的方法。#转载请注明出处!
  • 随机文章
  • 热门文章
  • 热评文章

“多元正态分布的定义(多元统计分析经典例题)(多元正态分布的例题答案)” 的相关文章

手机屏幕材质TFT好不好?真的很差吗?闲谈手机屏幕那些事(tft手机屏幕材质怎么样)

手机屏幕材质TFT好不好?真的很差吗?闲谈手机屏幕那些事(tft手机屏幕材质怎么样)

我买手机的一般先关心外观,然后是分辨率、屏幕尺寸、CPU品牌、核心、主频、RAM……屏幕材质往往是不怎么关心的,可能是我对手机要求不高吧。现在的手机厂商在宣传自己产品的时候,总是强调IPS屏幕,可我...

网页右下角弹出广告窗,可以自行关闭的右下角弹窗(如何屏蔽右下角弹窗广告)

网页右下角弹出广告窗,可以自行关闭的右下角弹窗(如何屏蔽右下角弹窗广告)

打开网页右下角弹出窗口,相信经常浏览网页的人已经非常熟悉,这咱弹窗效果最常用的就是广告和网站公告了,一打开网页就会弹出来,广告效果非常好,通过javascript实现的代码: <!D...

wifi是wlan吗?wlan和wifi有什么区别?(WLAN和wifi有什么区别)

wifi是wlan吗?wlan和wifi有什么区别?(WLAN和wifi有什么区别)

以前用诺基亚手机的时间,主界面会显示wlan,当时就以为wlan是wifi,两者是一样的东西,实质上却是不一样的! wifi是什么:wifi是无线保真的缩写,是能将电脑、手机等以无线方式互联的...

健身38种整套私教系统训练(健身私教课)

健身38种整套私教系统训练(健身私教课)

课程介绍想要练出一身的肌肉,需要付出的汗水与努力自然是不言而喻的。可是,网上的健身教练资源却良莠不齐,对于一些健身小白来说,是非常难靠自己来进行整理与分离。本课程带来38套健身方法.让你的健身不再迷茫...

Mybatis Example的高级用法详解_java

Mybatis Example的高级用法详解_java

目录Mybatis Example的高级用法一. mapper接口中的函数及方法二. example实例方法三. 使用案例说说Mybatis Example常见用法一. 说明二. 排序查询三. 查询l...

Java动态脚本Groovy_java(java动态脚本groovy)

Java动态脚本Groovy_java(java动态脚本groovy)

目录1.Groovy特性2.核心涉及3.Java与Groovy转换第一步:引入Groovy依赖第二步:创建interface接口声明方法第三步:在resources目录下创建.groovy文件第四步:...