便民热线和泊松分布-北太天元学习25
在一段固定的时间内发生了事件的次数往往是一个随机变量。
例如，我们对大城市的便民热线12345的来电数量是一个随机变量。
假设λ为每分钟呼叫的平均次数（假设该数字为固定常数）。
我们想研究某个固定时间内呼叫次数这个随机变量的分布情况，
显然，这个随机变量的取值为非负整数。
为了研究这个问题，假设我们把1分钟平均分成n等分，假设
[0,1/n),[1/2, 2/n), ..., 这个n个小时间区间内，呼叫的概率是相等，
为了是的平均每分钟的呼叫次数是λ, 我们设置每个小时间区间内
发生呼叫的概率是 λ/n, 我们假设这个时间小区间的长度是充分的小，
保证每个小时间区间内要么发生呼叫(而且呼叫最多只发生一次)，要么没有发生呼叫，
发生呼叫的概率是 λ/n , 没有发生呼叫的概率是 1 - λ/n.
再假设不同小区间内是否发生呼叫是相互独立的。 这样我们就把1分钟
内发生呼叫多少次的问题转换成了 n个独立的重复试验的成功与失败的
问题，成功对应一次呼叫发生。 这样，我们就可以利用二项分布的概率
来计算
P(X = k) = C_n^n * p^k * (1-p)^k
其中, p = λ/n, C_n^k 是从n个不同的元素选择k个的组合数。
让 n 趋向于无穷大，我们得到的极限是
P(X = k ) = λ^k / k! exp(-λ).
这个推导过程用到了 (1 - λ/n)^n 在n 趋向于正无穷时的极限是exp(- λ),
如果是高中生看到这儿，暂时忽略如何求极限的问题，直接读下面一段。
我们来比较一下 B(n,p) 和 Poisson(λ) 在满足 λ = n*p 时候的概率, 记X ~ B(n,p), Y ~ Poisson(n*p) 比较 P(X=k) 和 P(Y=k) , 例如 n=1000, p=0.001,
用北太天元的计算如果如下:
北太天元计算二项分布B(1000,0.0010)和Poisson(1.0000)分布的概率
二项分布概率P(X = 1)=0.3677, 泊松分布P(Y=1)=0.3679
二项分布概率P(X = 2)=0.3681, 泊松分布P(Y=2)=0.3679
二项分布概率P(X = 3)=0.1840, 泊松分布P(Y=3)=0.1839
二项分布概率P(X = 4)=0.0613, 泊松分布P(Y=4)=0.0613
二项分布概率P(X = 5)=0.0153, 泊松分布P(Y=5)=0.0153
如果一个随机变量取值为0,1,2,... 的随机变量，且满足
P(X = k ) = λ^k / k! exp(-λ),
那么我们称X 服从参数为λ 的泊松分布，记作 X ~ Poisson(-λ).
这里的λ 有时被称为X的速率参数。X=k 表示在1分钟接到k次电话的时间
发生。 在其它情形下，根据上下文，我们还可以说k是一段时间内事件发生
次数或放射出某种粒子的数量。
假设电话以每分钟λ=5个电话的平均速率呼入便民热线12345。
（a） 找出一分钟内正好有4个电话呼入的概率。
（b） 找出在10分钟的时间内发生少于45个呼叫的概率。
我们假设一分钟呼入电话的次数是服从速率参数是5的泊松分布，
从而得到
P(X=4) = 5^4/4! *exp(−4)
用北太天元计算
Input [6] >> lambda = 5
lambda =
1x1 double
5
Input [7] >> lambda^4/factorial(4)*exp(-lambda)
ans =
1x1 double
0.1755
要回答（b），我们需要找到每10分钟的平均通话次数，即λ=5*10=50。
我们要找到P（X<45）
在北太天元中使用下面的命令计算
Input [16] >> lambda = 5*10;
Input [17] >> k=0:44;
Input [18] >> fk = zeros(size(k));
Input [19] >> for n=0:44
... >> fk(n+1) = factorial(n);
... >> end
Input [20] >> sum( lambda.^k ./ fk * exp(-lambda) )
ans =
1x1 double
0.2210
我们的概率大约是0.2210。
画图的脚本如下：
close all; clear all; clc
n = 1000; % 试验的次数
p = 0.001; % 做一次试验时，事件A发生的概率
f = binomial_p(n,p);
fp = poisson_p(n*p, n);
ind = find(f> 0.01)
bar(ind,[f(ind);fp(ind)]);
xlabel("事件A发生的次数X")
ylabel("概率")
str = sprintf( ...
"北太天元计算二项分布B(%d,%3.4f)和Poisson(%3.4f)分布的概率", n, p,n*p)
title(str)
legend('二项分布','泊松分布');
fprintf("北太天元计算二项分布B(%d,%3.4f)和Poisson(%3.4f)...
分布的概率\n",n,p,n*p);
for i=ind
fprintf("二项分布概率P(X = %d)=%1.4f, 泊松分布P(Y=%d)=%1.4f\n", ...
i, f(i), i, fp(i) );
end
% 伯努利试验，总试验次数时n, 事件A发生的总次数X
% 时一个随机变量，可能的取值为 0,1,..., n
% X=k 表示 总试验次数是n时，事件A总共发生了k次，
% X=k 发生的概率是
% p(X=k) = nchoosek(n,k) * p^k * q^(n-1)
%
function f = binomial_p(n,p)
f = zeros(1,n+1); % 对二项分布的概率f赋初值
for k=0:1:n; % 事件A发生k次,求二项分布f,即成功k次的概率
f(k+1) = nchoosek(n,k) *p^k * (1-p)^(n-k);
end
end
% 从n个不同元素中取出k个元素的组合数
% cnk = n!/(k! * (n-k)!)
% 第二个算法 cnk = 1/(n-1)*beta(k+1,n-k) 这里的beta是beta函数
function cnk = nchoosek(n,k)
if(n==k || k==0)
cnk = 1;
return;
end
%cnk = (factorial(n)/factorial(n-k)/factorial(k))
cnk = 1/(n-k)/beta(k+1,n-k)
end
function f = poisson_p(lambda,n)
f = exp(-lambda)*ones(1,n+1);
for k=1:n
f(k+1) = f(k)*lambda/k;
end
end
北太天元运行这个脚本的截屏图如下：
