Q-Q plot Multivariate Statistical Analysis 学习笔记 1
发布时间
阅读量：
阅读量
2020年11月4日
Q-Q plot,用于检测样本数据是否符合某一分布，定义链接：https://ww2.mathworks.cn/help/stats/qqplot.html
就以正态分布为例，进行试验。步骤如下：
正态分布测试I（已知数据为正态分布）
- 生成正态分布的随机数列 X，长度为N（下文取N=1e5），作为样本
- sort（X）
- 取n=20个观测点，导出对应Xi作为纵坐标
- 计算X的平均值和方差，导出对应分位值作为横坐标
- 作图
安慰剂测试II（数据不是正态分布，但比较像）
- 生成随机数列 X，长度为N（下文取N=1e5），作为样本
- sort（X）
- 取n=20个观测点，导出对应Xi作为纵坐标
- 计算X的平均值和方差，导出对应分位值作为横坐标
- 作图
正态分布测试I
生成正态分布的随机数列 X
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
//生成一维的正态分布随机数
//pow(e,n)=exp(n)
double p=70,d=10;
const int N=100000;
const double e=2.718281828;
const double pi=3.1415926535;
double exp(double k){
return pow(e,k);
}
vector<double>X;
int main(){
X.clear();
freopen("正态分布随机数 .txt","w",stdout);
default_random_engine E;
int cnt=0;
double y_lim1,y_lim2;
double x_lim=3*d;
y_lim1=(1/(d * (sqrt(2*pi)) ));
y_lim2=(1/(d * (sqrt(2*pi)) )) * exp(-(x_lim*x_lim) / (2*d*d));
//在3d外分类随机，提高效率
while(cnt<N){
uniform_int_distribution<unsigned>x_k(0,100);
if( x_k(E) <99){
uniform_real_distribution<double>xx(p-x_lim,p+x_lim);
double x=xx(E);
double y_up=(1/(d * (sqrt(2*pi)) )) * exp(-( (x-p) * (x-p) ) / (2*d*d));
uniform_real_distribution<double>y(0,y_lim1);
if(y(E)<=y_up){
X.push_back(x);
cnt++;
}
}else{
uniform_real_distribution<double>xx(p-100*d,p+100*d);
double x=xx(E);
double y_up=(1/(d * (sqrt(2*pi)) )) * exp(-( (x-p) * (x-p) ) / (2*d*d));
uniform_real_distribution<double>y(0,y_lim2);
if(y(E)<=y_up){
X.push_back(x);
cnt++;
}
}
}
for(int i=0;i<X.size();i++){
printf("%lf ",X[i]);
}
return 0;
}此处取数学期望为70，标准差为10的正态分布（以学生体育成绩为原型题主就是那个40-的，悲 ）
此处的生成函数其实和真正的正态随机变量还是有区别的。为了测试其相似性，特意写了一个mat的scatter。
导出离散化的统计值：
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=100000;
double in;
int t;
int cnt[14500];
int main(){
memset(cnt,0,sizeof(cnt));
freopen("正态分布随机数 .txt","r",stdin);
freopen("loadin.txt","w",stdout);
for(int i=0;i<N;i++){
scanf("%lf",&in);
t=in*100;
cnt[t]++;
}
cout<<"x=[";
for(int i=3000;i<11000;i++){
t=i/100;
printf("%d ",i);
}
cout<<"]; \n y=[";
for(int i=3000;i<11000;i++){
printf("%d ",cnt[i]);
}
cout<<"]";
return 0;
}加入matlab，生成散点。
（横坐标被拉大了100倍 懒得改了 ）
导出纵坐标
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=100000;
double ave,var;
double x[N+50];
int main(){
memset(x,0,sizeof(x));
ave=var=0;
freopen("正态分布随机数 .txt","r",stdin);
freopen("sample_operation.txt","w",stdout);
for(int i=0;i<N;i++){
scanf("%lf",&x[i]);
ave+=x[i]/N;
}
for(int i=0;i<N;i++){
var+=( (x[i]-ave)*(x[i]-ave) )/N;
}
//printf("ave=%lf \nvar=%lf \n",ave,var);
sort(x,x+N);
for(int i=0;i<20;i++){
printf("%lf ",x[2499+5000*i]);
}
return 0;
}导出横坐标
二分法找出目标拟合正态分布的分位点的值，作为横坐标
x=0:0:20;
ave=70.007051;
var=98.074006;
dev=sqrt(var);
for i=1:20
targ=1/40+1/20*(i-1);
l=ave-5*dev;
r=ave+5*dev;
m=0;
while(r-l>1e-4)
m=(l+r)/2;
%syms x;
%fun=(1/dev/sqrt(2*pi))*exp(-(x-ave)*(x-ave)/2/var);
%F=int(fun,-inf,m);
fun = @(x)(1/dev/sqrt(2*pi))*exp(-(x-ave).*(x-ave)/2/var);
F = integral(fun,-inf,m);
if F>targ
r=m;
else l=m;
end
end
x(i)=m;
end
y=[50.563635 55.626153 58.511703 60.680574 62.443618 64.019945 65.443353 66.805756 68.107271 69.360707 70.634768 71.886007 73.182585 74.567182 76.011793 77.618535 79.363309 81.510039 84.332060 89.503564 ];
scatter(x,y);结果如下：
可以看出，这是一个正态分布（这tm不是废话吗 ）
安慰剂测试II
正在施工
全部评论 (0)
还没有任何评论哟~