[概率论] 和分布的分位数问题 #5

Joe-C-Ding · 2017-11-26T09:45:33Z

题面

概率论的题目，想了一天也没想出来

X, Y 是随机变量，且相互独立。Z 为它们的和分布，即 Z=X+Y。已知 Pr(X<a) = Pr(Y<b) = p (p > 0.5)，则 Pr(Z<a+b) >= p。

这是真命题么？我本来觉得伪的，不过没找到反例呢，现在觉得是真命题的可能性也不小。

背景

背景不太容易说清，但大概可以这么看。p > 0.5 这个条件就是说前面两个事件可以算的上是“大概率”事件。但一般意义上

X < a && Y < b => X+Y < a+b

这样的代数不等式相加，要两个条件同时成立才成立。所以只有

Pr(Z < a+b) > p^2

也就说，即使原来两个条件是“大概率”事件，结果小于的概率也要打不少折扣。而这个命题在表明其实这样的折扣是不存在的。

现有研究结论

X Y 如果都正态分布命题为真，且结论那里取不到等号。如果是均匀分布，正好取等号。

正态分布的

分位数的意思就是 Z < zp 的概率正好是 p。因为 zp < xp+yp，所以 Z < xp+yp 的概率要更大些。

均匀分布

均匀分布就是简单线性叠加，所以分位数不变，概率也还是 p，不细说了。

因为均匀分布的分散性其实最大，所以它取等号。如果是其它分布，我试的一些具体都取不到那个等号。命题里没有 X Y 的具体分布，也不要求同分布。感觉如果是真命题应该从方差那方面着手考虑证明。有人有证明的思路也可以说说看。

Joe-C-Ding · 2017-11-26T10:58:21Z

X Y 是两个相互独立的随机变量，之前忘了交代。

Joe-C-Ding · 2017-11-27T02:34:06Z

以下是一个反例：

现在的结论似乎是：

如果和分布的概率密度为单峰函数，则结论成立。多峰则不成立。

虽然前者也依然不易证明吧？

shouldsee · 2017-11-28T15:55:48Z

增强的推论:
可以考虑证一个更一般的情形 Pr( a*X+b*Y < a*x_p + b*x_p ) >= p ，且当 a=0或b=0时等号成立。

弱化的推论:
另外有一个特殊情况是 p = 0.5 的时候等号是一直成立的

shouldsee · 2017-12-02T12:47:02Z

令 X,Y 为不相关的变量

Median(X)=log2(e) (see wiki_frechet )
Median(Y)=0

根据WolframAlpha计算结果：
P( X + Y < Median(X) + Median(Y)) = 0.104466 + 0.372747 = 0.477213 < 0.5
由此可以猜测原结论需要附加有关对称性的约束。

积分计算：
验证联合概率分布的积分为1

计算P( X + Y < Median(X) + Median(Y))

Joe-C-Ding · 2017-12-02T16:41:32Z

挺有意思的例子，尤其是 Frechet 分布，只有小于 a 阶的矩。

shouldsee · 2017-12-02T22:37:42Z

和分布好像可以用特征函数来更好地刻画。(见稳定分布)

@Joe-C-Ding 我其实是wiki上扒了一个非对称的分布下来，选择Frechet主要因为它的PDF比较好写 :X

shouldsee · 2018-03-19T16:09:00Z

@Joe-C-Ding

的确f(0.5)=0一般不成立的。（见该笔记)

代了几个Gamma分布看起来根的数量也不太好确定(诡异的是f(1)=0看起来并不总是成立，应该跟数值计算有关？)

Joe-C-Ding · 2018-03-19T22:51:54Z

@shouldsee
f(1) = 0 相当于 Pr(X < xp) = 1, 这时 xp 就是 x_max 或者 inf。所以总有 X+Y<xmax+ymax
同理 0 也如此。但一般代进去算都是奇点吧，所以计算有误差很正常。
但你举的最后一个例子确实值得研究一下是什么情况。

由于 0 1 是平凡根，所以要看 (0,1) 之间有几个根。看你举的那几个例子都是一个根的情况。

shouldsee · 2018-03-20T15:09:35Z

@Joe-C-Ding
我发现matlab积分用的是adaptive quadrature条目（见此页reference）。SO上有提到python的quad()函数对python function来说很慢，quadrature()相对会快一点但是精度有时候不够(1E-4)。另外可以考虑写ctypes函数

UPDATE: scipy好像慢在scipy.stats.weibull_min.cdf这里，可以考虑用julia替代（开源且overhead小）notebook

Joe-C-Ding · 2018-03-21T03:49:54Z

@shouldsee
我以前以为 quad() 是自适应积分呢.文档只说了用的是 fotran 的什么包，具体算法没提。quadrature() 是高斯积分，区别就有点大了。

matlab 因为没有名空间类似的技术，所以加载包和确定函数重载这块好像特别慢。导致一个函数第一次运行要比第二次运行慢很多。你跑 python 能比我这快小一半，matlab 反而比我慢一倍多除了有系统启动的时间以外可能还有这个原因。

我是在 python 和 matlab 里面跑的，python 计时也只计了 main 里那部分。不过优化的问题就复杂了，因为很多研究初期要计算的东西都不是那么确定的。原生就快还是有优势。

Joe-C-Ding · 2018-03-21T03:54:27Z

@shouldsee
你后来算的那个还是看不出来是积分慢还是 cdf 本身算的慢吧。这些常用分布的 cdf 都有解析表达式，像 weibull 就是个 exp 的事，底层我觉得都是 libc 里的函数，x86的硬件就有汇编指令能直接算。gamma 的 cdf 可能稍微费力些。

Joe-C-Ding · 2018-03-21T08:57:31Z

@shouldsee 你的程序要不要再检查下？看你后来的一个例子改 loc 会影响 G(p)。我觉得是不会吧，因为 loc 只是平移操作。X < xp 的概率和 X+5 < xp+5 的概率是一样的。所以也不改变 X+Y < xp+yp 的概率。

shouldsee · 2018-03-21T09:12:47Z

本来大早上准备睡个懒觉结果被无数封邮件拉了起来。。。

@Joe-C-Ding
用GammaCDF做测试的原因是这个函数比较好找，三种语言都有原生提供。一开始积的是x^3，然后发现三种都没啥区别，所以推断是scipy.stats里面的overhead比较高(有很多check之类会杀速度)。

benchmark那个笔记里面测了内侧和外侧时间， matlab的速度还是fair的。

改loc为啥会影响G(p)我也是到现在没搞懂。
UPDATE: 改了一个bug，原代码积分的时候写错了定义域导致和分布的pdf并不正规化，进一步导致cdf的值域不是(0,1)。现在G(p)跟loc没啥关系了

UPDATE：python里直接call了gammacdf的解析形式居然快了40倍。然而就目前的需求来说，主要问题还是搞出ppf的形式，不过估计也只能数值逼近了。

Joe-C-Ding · 2018-03-21T09:54:15Z

matlab 同样脚本第一次跑和第二次跑速度能提高至少一两倍。因为改动脚本再执行通常也不会改变执行速度了，所以我觉得不是解释器缓存了语法分析的原因。可能就是第一次加载函数慢的事。
我之前的意思你那种计时方法，不管里外，都算上了解释器解析函数名的时间。不过就结果来说，即使这样看 matlab 也不是很慢。

用 matlab 算了同一个算例，你对对哪不一样。

  Fz = Pr(X + Y < z) = int { Pr(Y = t) * Pr(X < z-t) dt }

用这个公式算的 Z 的 cdf。看你后来的几个算例，Z 的 cdf 都不趋于1，肯定不对吧……

start_tic = tic;
clf;

% p1 =spstat.gamma(loc=-2, a= 7)
% p2 =spstat.gamma(loc=0, a= 20)

func = @(z, t)gampdf(t,20,1) * gamcdf(z+2-t,7,1);
Fz = @(z)integral(@(t)func(z, t), -inf, inf, 'arrayvalued', 1);

disp([0:10:50; Fz(0:10:50)].');
% >>>         z           Fz
%             0   1.7961e-21
%            10   0.00013335
%            20      0.16758
%            30      0.83443
%            40      0.99445
%            50      0.99995

p = linspace(0, 1, 101).';
xp = gaminv(p,7,1)-2;
yp = gaminv(p,20,1);
Gp = Fz(xp+yp) - p;
plot(p, Gp);

i = ismembertol(p, 0:0.2:1, 1e-8);
disp([p(i), xp(i), yp(i), Gp(i)]);
% >>>         p           xp           yp           Gp
%             0           -2            0            0
%           0.2       2.7337       16.172    -0.086729
%           0.4       4.0392       18.567    -0.059228
%           0.6       5.3426       20.811     0.011386
%           0.8       7.0754       23.634     0.062772
%             1          Inf          Inf  -1.8874e-15

fprintf('%s elapsed: %f s\n', mfilename, toc(start_tic));

shouldsee · 2018-03-21T11:01:50Z

@Joe-C-Ding

p1 =spstat.gamma(loc=-2, a= 7)
p2 =spstat.gamma(loc=0, a= 20)
main(p1,p2,cutoff=1E-8,num=200)
scdf = make_scdf(p1,p2,cutoff=1E-8)
Gp = lambda p: scdf(p1.ppf(p)+p2.ppf(p))-p
ps = np.arange(0,1.1,0.2)
arr = [ps,p1.ppf(ps),p2.ppf(ps),Gp(ps)]
arr = np.transpose(arr)

from pymisc.util import mat2str

print mat2str(arr,4)

'''
0.0	-2.0	0.0	0.0
0.2	2.7337	16.1725	-0.0869
0.4	4.0392	18.567	-0.0591
0.6	5.3426	20.8111	0.0111
0.8	7.0754	23.6343	0.0625
1.0	inf	inf	-0.0
'''

另外我手动关掉了pdf的error check大概快了5-10倍吧

Joe-C-Ding · 2018-03-21T15:17:40Z

clear all;

need_import = true;
if need_import
    import_time = tic;
    a = linspace(0,1);
    b = gampdf(1,1,1);
    c = gamcdf(2,2,2);
    d = integral(@(x) x.^(1:10), 0, 1, 'arrayvalued', 1);
    e = gaminv(0.5,3,3);
    fprintf('import time: %f s\n', toc(import_time));
end

start_tic = tic;

% p1 =spstat.gamma(loc=-2, a= 7)
% p2 =spstat.gamma(loc=0, a= 20)

func = @(z, t)gampdf(t,20,1) * gamcdf(z+2-t,7,1);
Fz = @(z)integral(@(t)func(z, t), -inf, inf, 'arrayvalued', 1);

p = linspace(0, 1, 101).';
xp = gaminv(p,7,1)-2;
yp = gaminv(p,20,1);
Gp = Fz(xp+yp) - p;
% plot(p, Gp);

fprintf('%s elapsed: %f s\n', mfilename, toc(start_tic));

结果

import = true
import time: 0.110498 s
test elapsed: 0.304823 s

import = false
test elapsed: 0.409313 s

不管 if ... end 之间的语句执不执行花的总时间都差不多。但是如果把最开始的 clear all 去掉，再次执行还能更快。看来 matlab 除了缓存了函数以外，还有其它加速执行的方法。

你这手动关的技术高级了，我都不知道。scipy 我也不是太熟悉，就是哪些包基本能干什么了解个大概。

shouldsee · 2018-03-21T15:26:38Z

这叫为了生存激发的潜能xD。主要我实在实在不想重新用matlab...所以为了继续用python，总不能这么一直龟速下去...幸好python不负我

shouldsee · 2018-03-21T15:28:46Z

侧面也看出来python开源的好处就是碰上这种overhead还有调节的可能性，而matlab碰上一些封装的c就基本上完了个蛋

Joe-C-Ding · 2018-03-21T15:31:04Z

也不是，matlab 只有很小一部分代码是封装在 c 里的。一般会用到的都是 m 文件，有需要能自己改。
我看了它求 pdf 的代码，直接算的就一行，剩下 90% 的代码也是在检查输入的参数是不是正确。

shouldsee · 2018-03-21T15:31:15Z

具体实现在notebook里写了，至于探索流程参见https://zhuanlan.zhihu.com/p/34800439

shouldsee · 2018-03-21T15:36:01Z

Matlab吧...查个源码翻个doc的速度，简直就是灾难大片。然后也没法架jupyter这种网页端，对我这种笔记本轻巧的来说根本无所适从…我们也就处理图像的课程会用matlab，其它大部分时间还是R/Python。综合来讲有三点：

各种overhead，包括启动，查文档。
syntax不习惯，尤其.m定义函数
没有网页端。

不过幸好你对比了一下两者的性能，不然到时候一直在龟速跑都不知道囧。

Joe-C-Ding · 2018-03-21T15:41:34Z

嗯。matlab 作为通用编程语言还是有些烂，只用来做些数学方面一些数值验证还可以。

感觉有些条件是不是还是得查呀，比如 x 如果不在定义域里 pdf 得返回 0 吧。

shouldsee · 2018-03-21T15:44:03Z

你可以考虑一下如何进行快速的error-check。目前因为积分定义域是从ppf里生成的所以不会在定义域外部。

Joe-C-Ding · 2018-03-21T15:52:35Z

我用的那个公式就得看了：

Fz = Pr(X + Y < z) = int { Pr(Y = t) * Pr(X < z-t) dt }

只是简单的从 -inf 到 inf 积的。要用到如果出定义域自动变 0 的功能。

python 的 quad 好像不支持两端都是 inf 的情况，还是像你那样 cutoff 好一些。matlab 不能在脚本里定义函数确实是相当蛋疼的一个缺陷。用函数句柄的话，只能是一个语句。

shouldsee · 2018-03-21T15:59:58Z

python有np.inf，你可以用正态函数试试速度，之前是因为慢才用的cutoff。

理论上你可以somehow把积分域改成随机变量的值域的。

shouldsee · 2018-04-01T19:06:10Z

A somewhat related application of characteristic function and illustration of CLT can be found here

Joe-C-Ding · 2018-04-02T02:11:25Z

CLT 不好用的原因是，当 n 不太大的时候（之前的研究相当于 n=2），CLT 做近似的误差太大。n 稍大时误差通常会很快减小。但像我研究的安全性领域，通常关心拖尾处的概率（比如 p > 0.95 或 p > 0.99），负方向的小误差有时也不能接受。
如果想仔细研究的话，极值分布一般更有参考价值。但通常也要求 n 稍大。咱们之间给的结论是保守的，既 Pr(X+Y<xp+xp) >= p。它成立的话，是正方向的误差。即使大也不要紧。
用 CF 算当然不错，不过通过之前的研究结论，估出的 G(p) = 0 的根，很多时候对分布参数变化不敏感。这样不用估参数，或只用估个大概可以省去很多样本采集的工作。

shouldsee · 2018-04-02T13:12:49Z

@Joe-C-Ding Cool... I was just fascinated by the fact that Cauchy distribution does not converge to Gaussian under CLT. Given the popularity of CLT and its connection to CF
, it's natural to wonder whether any similar conclusion exists for a CDF . This requires finding a linear functional for the indicator function ( For exponential functions, this linear functional is simply the product of two component functions). Unfortunately I have no clue at all at the moment.

Anyway, it seems CF is a more natural choice since we are summing up random variable. If one can transform a CF back into CDF, it might be even possible to derive an analytical form of G(p).

I think you meant G(p)>0 (p>=0.5 in general)? It makes sense to talk about the sample size N, I think the fact that CLT does not work at N=2 is reminiscent of the general finite-size phenomena... The discrete domain is very different from the continuous one, and indeed we need some different tools here.

(Apologies for the English post for I am on my desktop and don't have Chinese input)

Joe-C-Ding · 2018-04-20T03:19:26Z

CLT 背后的极限理论比较复杂（从17xx 年到 19xx 年研究了200多年），我研究极值分布的时候见有文献提过一些但我没再专门研究过。最经典的 CLT 一般要求独立同分布的，但其实不同分布甚至不太独立的随机变量加起来也可能收敛到正态上。CLT 成立的充要条件现在人们好像也有一些结论了，好像和每个变量的方差有比较大的关系，我没具体了解过。

理论上讲 CF 到 CDF 的反向变换积分一般不收敛，通常要用柯西主值(Cauchy principal value)。我可能分析学得也不深吧，一般的教材上对这种积分提的比较少（可能是因为更容易？）。所以能不能用来研究 G(P) 的其它性质我也没仔细想过。

你说样本数 N，不知道指什么。明确一下符号吧，记 k 为随机变量的个数，咱之前讨论的是 X+Y，相当于 k=2。还一个是一般意义上的样本个数 N，就是比如有 {x1,...,xn} iid ~ X, {y1,...,yn} iid ~ Y，那么咱们之前讨论的问题大概是想知道 {z1,...,zn} 这 N 个样本中小于 xp+yp 的样本个数的期望。并且关心这个期望是不是会大于 p * N，或者要求大于 pN 的话对 p 是不是有什么特别的要求。

CLT 是研究当 k->infinity 时和分布的极限，和 N 没关系的。离散分布和连续分布的区别的确是比较大，我也没有深入的研究。

shouldsee · 2018-04-26T12:46:37Z

我对CF到CDF的变换和柯西主值完全没有了解啊，能否举个例子看看？

我感觉平时用的CLT多是考虑N趋近无穷的，当然k和N的区别仅在于对称性，也即若变量在permutation下对称则k退化成N。但是即便要考虑k的话总要对变量的对称性作出稍弱的假设吧。

我还是想强调一下不是所有的CLT都收敛到高斯分布，个人感觉这是一个可能的切入点。

Joe-C-Ding · 2018-04-27T00:32:55Z

柯西主值是这样：

由于是无穷区间上的积分，因此不一定对所有 t 都收敛的。f(x) 绝对可积时这一般不成问题。
如果不是对所有 t 收敛的话，就说这个分布是重尾的 heavy tailed，这时可能不存在有限的方差。重尾分布一般认为多少都是有点病态的，既使在极值理论中也算常见。一般我们用 CLT 的时候是要求有限方差的，我感觉这大概也是收敛到正态的一个必要条件。

反向变换是这个（希望没写错哈）：

问题是：即使正向变换对所有 t 都收敛，也不能保证这个反向变换对 x 都是收敛的。这时就要用柯西主值了。
本来无穷积分是指：

但

可见前者收敛时后者必收敛且值一致，但反过来不一定。反向变换的无穷积分一般是 p.v. 意义下的。Fourier 变换的逆变换也是同理的。

我没太理解你说的 k 和 N 的关系。因为如果 cdf 知道了的话，就没 N 啥事了。只用研究 k 个分布和的问题。当然应用中，要样本 N 大一些才能表现出统计特征。但表现不出来的时候也不能说理论分布是不对的。

CLT 收敛到的分布类好像统称稳定分布 stable distribution。我知道有不止一个，只是没仔细研究过之间有什么关系，收敛到不同的类上的条件之类的。我看 wiki 上有些介绍，蛮复杂的样子。你要感兴趣可以看看，没准得再补充些其它资料才能看懂。

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

[概率论] 和分布的分位数问题 #5

[概率论] 和分布的分位数问题 #5

Joe-C-Ding commented Nov 26, 2017 •

edited

Loading

Joe-C-Ding commented Nov 26, 2017

Joe-C-Ding commented Nov 27, 2017 •

edited

Loading

shouldsee commented Nov 28, 2017 •

edited

Loading

shouldsee commented Dec 2, 2017 •

edited

Loading

Joe-C-Ding commented Dec 2, 2017

shouldsee commented Dec 2, 2017

shouldsee commented Mar 19, 2018 •

edited

Loading

Joe-C-Ding commented Mar 19, 2018 •

edited

Loading

shouldsee commented Mar 20, 2018 •

edited

Loading

Joe-C-Ding commented Mar 21, 2018

Joe-C-Ding commented Mar 21, 2018

Joe-C-Ding commented Mar 21, 2018

shouldsee commented Mar 21, 2018 •

edited

Loading

Joe-C-Ding commented Mar 21, 2018 •

edited

Loading

shouldsee commented Mar 21, 2018

Joe-C-Ding commented Mar 21, 2018

shouldsee commented Mar 21, 2018

shouldsee commented Mar 21, 2018

Joe-C-Ding commented Mar 21, 2018

shouldsee commented Mar 21, 2018

shouldsee commented Mar 21, 2018

Joe-C-Ding commented Mar 21, 2018

shouldsee commented Mar 21, 2018

Joe-C-Ding commented Mar 21, 2018 •

edited

Loading

shouldsee commented Mar 21, 2018 •

edited

Loading

shouldsee commented Apr 1, 2018 •

edited

Loading

Joe-C-Ding commented Apr 2, 2018 •

edited

Loading

shouldsee commented Apr 2, 2018 •

edited

Loading

Joe-C-Ding commented Apr 20, 2018 •

edited

Loading

shouldsee commented Apr 26, 2018

Joe-C-Ding commented Apr 27, 2018 •

edited

Loading

[概率论] 和分布的分位数问题 #5

[概率论] 和分布的分位数问题 #5

Comments

Joe-C-Ding commented Nov 26, 2017 • edited Loading

题面

背景

现有研究结论

正态分布的

均匀分布

Joe-C-Ding commented Nov 26, 2017

Joe-C-Ding commented Nov 27, 2017 • edited Loading

shouldsee commented Nov 28, 2017 • edited Loading

shouldsee commented Dec 2, 2017 • edited Loading

Joe-C-Ding commented Dec 2, 2017

shouldsee commented Dec 2, 2017

shouldsee commented Mar 19, 2018 • edited Loading

Joe-C-Ding commented Mar 19, 2018 • edited Loading

shouldsee commented Mar 20, 2018 • edited Loading

Joe-C-Ding commented Mar 21, 2018

Joe-C-Ding commented Mar 21, 2018

Joe-C-Ding commented Mar 21, 2018

shouldsee commented Mar 21, 2018 • edited Loading

Joe-C-Ding commented Mar 21, 2018 • edited Loading

shouldsee commented Mar 21, 2018

Joe-C-Ding commented Mar 21, 2018

shouldsee commented Mar 21, 2018

shouldsee commented Mar 21, 2018

Joe-C-Ding commented Mar 21, 2018

shouldsee commented Mar 21, 2018

shouldsee commented Mar 21, 2018

Joe-C-Ding commented Mar 21, 2018

shouldsee commented Mar 21, 2018

Joe-C-Ding commented Mar 21, 2018 • edited Loading

shouldsee commented Mar 21, 2018 • edited Loading

shouldsee commented Apr 1, 2018 • edited Loading

Joe-C-Ding commented Apr 2, 2018 • edited Loading

shouldsee commented Apr 2, 2018 • edited Loading

Joe-C-Ding commented Apr 20, 2018 • edited Loading

shouldsee commented Apr 26, 2018

Joe-C-Ding commented Apr 27, 2018 • edited Loading

Joe-C-Ding commented Nov 26, 2017 •

edited

Loading

Joe-C-Ding commented Nov 27, 2017 •

edited

Loading

shouldsee commented Nov 28, 2017 •

edited

Loading

shouldsee commented Dec 2, 2017 •

edited

Loading

shouldsee commented Mar 19, 2018 •

edited

Loading

Joe-C-Ding commented Mar 19, 2018 •

edited

Loading

shouldsee commented Mar 20, 2018 •

edited

Loading

shouldsee commented Mar 21, 2018 •

edited

Loading

Joe-C-Ding commented Mar 21, 2018 •

edited

Loading

Joe-C-Ding commented Mar 21, 2018 •

edited

Loading

shouldsee commented Mar 21, 2018 •

edited

Loading

shouldsee commented Apr 1, 2018 •

edited

Loading

Joe-C-Ding commented Apr 2, 2018 •

edited

Loading

shouldsee commented Apr 2, 2018 •

edited

Loading

Joe-C-Ding commented Apr 20, 2018 •

edited

Loading

Joe-C-Ding commented Apr 27, 2018 •

edited

Loading