Back

# 统计函数

## 随机变量

### 离散随机变量

$P(X=k) = {N\choose k} p^k(1-p)^{n-k}.$ {#eq:stats:binomial_pdf}

let _ =
let toss = Array.make 10 0 in
Array.map (fun _ -> Stats.binomial_rvs 0.3 1) toss
>- : int array = [|0; 0; 0; 0; 0; 0; 0; 0; 1; 0|]


$p_X(k)=P({s \in S | X(s) = k}),$

let x = [|0; 1; 2; 3|]
>val x : int array = [|0; 1; 2; 3|]
let p = Array.map (Stats.binomial_pdf ~p:0.3 ~n:3) x
>val p : float array =
>  [|0.342999999999999916; 0.440999999999999837; 0.188999999999999918;
>    0.0269999999999999823|]
Array.fold_left (+.) 0. p
>- : float = 0.999999999999999778


$p(X\leq~k)=\sum_{i=0}^k{N\choose i} p^k(1-p)^{n-i}.$

let x = [|0; 1; 2; 3|]
>val x : int array = [|0; 1; 2; 3|]
let p = Array.map (Stats.binomial_cdf ~p:0.3 ~n:3) x
>val p : float array = [|0.342999999999999972; 0.784; 0.973; 1.|]


### 连续随机变量

let noise sigma = Stats.gaussian_rvs ~mu:0. ~sigma;;
let x = Array.init 999 (fun _ -> Stats.gaussian_rvs ~mu:1. ~sigma:1.);;
let y = Array.init 999 (fun _ -> Stats.gaussian_rvs ~mu:12. ~sigma:3.);;

(* 将数组转换为矩阵 *)

let x' = Mat.of_array x 1 999;;
let y' = Mat.of_array y 1 999;;

(* 绘制图形 *)

let h = Plot.create ~m:1 ~n:2 "plot_02.png" in

Plot.subplot h 0 0;
Plot.set_ylabel h "频率";
Plot.histogram ~bin:30 ~h x';
Plot.histogram ~bin:30 ~h y';

Plot.subplot h 0 1;
Plot.set_ylabel h "PDF p(x)";
Plot.plot_fun ~h (fun x -> Stats.gaussian_pdf ~mu:1. ~sigma:1. x) (-2.) 6.;
Plot.plot_fun ~h (fun x -> Stats.gaussian_pdf ~mu:12. ~sigma:3. x) 0. 25.;

Plot.output h;;

$p(x\leq~k)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^k~e^{-t^2/2}dt.$

let h = Plot.create "plot_gaussian_cdf.png" in
Plot.set_ylabel h "CDF";
Plot.plot_fun ~h ~spec:[ RGB (66,133,244); LineStyle 1; LineWidth 2.; Marker "*" ] (fun x -> Stats.gaussian_cdf ~mu:1. ~sigma:1. x) (-2.) 6.;
Plot.plot_fun ~h ~spec:[ RGB (219,68,55);  LineStyle 2; LineWidth 2.; Marker "+" ] (fun x -> Stats.gaussian_cdf ~mu:12. ~sigma:3. x) 0. 25.;
Plot.(legend_on h ~position:SouthEast [|"mu=1,sigma=1"; "mu=12, sigma=3"|]);
Plot.output h

### 描述性统计

$E(X) = \frac{1}{n}\sum_{i}x_i,$ {#eq:stats:mean} 其中$$x_i$$是总体中的第$$i$$个元素。方差的定义类似：

$Var(X) = \frac{1}{n}\sum_{i}(x_i - E(X))^2.$ {#eq:stats:variance}

$M_n(X) = \int_x~(x_i - c)^2~f(x_i)dx,$ {#eq:stats:moment}

let data = Array.init 100 (fun _ -> Stats.uniform_rvs 0. 10.);;

Stats.mean data
>- : float = 5.18160409659184573
Stats.std data
>- : float = 2.92844832850280135
Stats.var data
>- : float = 8.57580961271085229
Stats.skew data
>- : float = -0.109699186612116223
Stats.kurtosis data
>- : float = 1.75165078829330856


Stats.central_moment 0 data
>- : float = 1.
Stats.central_moment 1 data
>- : float = -3.13082892944294137e-15
Stats.central_moment 2 data
>- : float = 8.49005151658374224
Stats.central_moment 3 data
>- : float = -2.75496511397836663


Stats.min;;
Stats.max;;
Stats.median;;
Stats.quantile;;
Stats.first_quartile;;
Stats.third_quartile;;
Stats.percentile;;

minmax很容易使用。median是整个样本中排序列表中的中间数字。它有时比mean更能描述数据，因为后者更容易受到异常值的影响。

## 特殊分布

Gaussian distribution $$\frac{1}{\sigma {\sqrt {2\pi }}}e^{-{\frac {1}{2}}\left({\frac {x-\mu }{\sigma }}\right)^{2}}$$
Gamma distribution $$\frac{1}{\Gamma(k)\theta^k}x^{k-1}e^{-x\theta^-{1}}$$
Beta distribution $$\frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)\Gamma(\beta)}x^{\alpha-1}(1-x)^{\beta-1}$$
Cauchy distribution $$(\pi~\gamma~(1 + (\frac{x-x_0}{\gamma})^2))^{-1}$$
Student’s $$t$$-distribution $$\frac{\Gamma((v+1)/2)}{\sqrt{v\pi}\Gamma(v/2)}(1 + \frac{x^2}{v})^{-\frac{v+1}{2}}$$

• gaussian_rvs：随机数生成器。
• gaussian_pdf：概率密度函数。
• gaussian_cdf：累积分布函数。
• gaussian_ppf：百分点函数（CDF的反函数）。
• gaussian_sf：生存函数（1 - CDF）。
• gaussian_isf：逆生存函数（SF的反函数）。
• gaussian_logpdf：对数概率密度函数。
• gaussian_logcdf：对数累积分布函数。
• gaussian_logsf：对数生存函数。

Stats模块支持许多分布。对于每个分布，都有一组以分布名称作为它们共同前缀的相关函数。例如，下面的代码使用gamma_pdf绘制Gamma分布的概率密度函数。结果显示在[@fig:stats:gamma_pdf]中。


module N = Dense.Ndarray.D

let _ =
let x = N.linspace 0. 16. 100 in

let f1 x = Owl_stats.gamma_pdf x ~shape:1. ~scale:2. in
let f2 x = Owl_stats.gamma_pdf x ~shape:2. ~scale:2. in
let f3 x = Owl_stats.gamma_pdf x ~shape:5. ~scale:1. in
let f4 x = Owl_stats.gamma_pdf x ~shape:7.5 ~scale:1. in

let y1 = N.map f1 x in
let y2 = N.map f2 x in
let y3 = N.map f3 x in
let y4 = N.map f4 x in

let h = Plot.create "gamma_pdf.png" in
let open Plot in
set_xlabel h "";
set_ylabel h "";
set_title h "Gamma分布概率密度函数";
plot ~h ~spec:[ RGB (66, 133, 244); LineStyle 1; LineWidth 2. ] x y1;
plot ~h ~spec:[ RGB (219, 68,  55); LineStyle 1; LineWidth 2. ] x y2;
plot ~h ~spec:[ RGB (244, 180,  0); LineStyle 1; LineWidth 2. ] x y3;
plot ~h ~spec:[ RGB (15, 157,  88); LineStyle 1; LineWidth 2. ] x y4;

Plot.(legend_on h ~position:NorthEast [|"k=1, theta=2"; "k=2, theta=2"; "k=5, theta=1"; "k=7.5, theta=1"|]);

output h

## 多变量

$p(xy) = p(X=x \textrm{AND } Y=y) = p(X=x) \times p(Y=y) = p_X(x)p_Y(y).$

$P(X | Y) = \frac{P(X~\cap~Y)}{P(Y)}.$

$P(X|Y) = \frac{P(Y|X)P(X)}{P(Y)}$ {#eq:stats:bayes}

## 抽样

### 无偏估计

$m = \frac{1}{n}\sum_{i=1}^n~x_i,$

$v = \frac{1}{n}\sum_{i=1}^n(x_i - m)^2.$

### 推断总体参数

let f = Stats.gaussian_cdf ~mu:0. ~sigma:1. in
f 1.96 -. f (-1.96)
>- : float = 0.950004209703559


$P( x - 1.96~\frac{\sigma}{\sqrt{n}} \le \mu \le x + 1.96~\frac{\sigma}{\sqrt{n}}).$

let f x = Stats.t_cdf x ~df:16. ~loc:0. ~scale:1. in
f 2.12 -. f (-2.12)
>- : float = 0.950009071286895823


## 假设检验

### 理论

• 零假设 $$H_0$$：两个数据集之间没有关系。
• 备择假设 $$H_1$$：两个数据集之间存在统计显著关系。

### 在假设检验中的正态分布

let data_0 = Array.init 10 (fun _ -> Stats.gaussian_rvs ~mu:0. ~sigma:1.);;
let data_1 = Array.init 10 (fun _ -> Stats.gaussian_rvs ~mu:3. ~sigma:1.);;

Stats.z_test ~mu:0. ~sigma:1. data_0
>- : Owl_stats.hypothesis =
>{Owl.Stats.reject = false; p_value = 0.289340080583773251;
> score = -1.05957041132113083}


type hypothesis = {
reject : bool;
p_value : float;
score : float;
}

Stats.z_test ~mu:0. ~sigma:1. data_1
>- : Owl_stats.hypothesis =
>{Owl.Stats.reject = true; p_value = 5.06534675819424548e-23;
> score = 9.88035435799393547}


### 双样本推断

$\frac{\bar{x} - \bar{y}}{\sqrt{\frac{\sum_{i=1}^a~(x_i - \bar{x})^2 + \sum_{i=1}^b~(y_i - \bar{y})^2}{a + b - 2}(\frac{1}{a} + \frac{1}{b})}},$

### 其他测试类型

• Kolmogorov-Smirnov测试：ks_test ~alpha x f返回一个测试决策，用于零假设，即向量x中的数据来自具有CDF f的分布的独立随机样本。备择假设是x中的数据来自不同的分布。结果(h,p,d)d是Kolmogorov-Smirnov测试统计量。

• 双样本Kolmogorov-Smirnov测试：ks2_test ~alpha x y返回一个测试决策，用于零假设，即向量xy中的数据来自相同分布的独立随机样本。

• 卡方方差测试var_test ~alpha ~side ~variance x返回一个测试决策，用于零假设，即x中的数据来自具有输入variance的正态分布，使用卡方方差测试。备择假设是x来自具有不同方差的正态分布。

• Jarque-Bera测试jb_test ~alpha x返回一个测试决策，用于零假设，即数据x来自具有未知均值和方差的正态分布，使用Jarque-Bera测试。

• Wald–Wolfowitz Runs测试runs_test ~alpha ~v x返回一个测试决策，用于零假设，即数据x是随机排序的，反对它们不是，通过运行Wald–Wolfowitz测试。该测试基于连续值在均值上方或下方的运行的数量。~v是参考值，默认值是x的中位数。

• Mann-Whitney秩测试mannwhitneyu ~alpha ~side x y对样本x和y执行Mann-Whitney秩测试。如果每个样本的长度小于10且没有结合，则使用精确测试，否则使用渐近正态分布。

## 协方差和相关性

x是我们的解释变量，我们从0到10的区间中均匀抽取了50个随机值。yz都是与x线性相关的响应变量。唯一的区别是我们向响应变量添加了不同水平的噪声。噪声值是从高斯分布生成的。

let noise sigma = Stats.gaussian_rvs ~mu:0. ~sigma;;
let x = Array.init 50 (fun _ -> Stats.uniform_rvs 0. 10.);;
let y = Array.map (fun a -> 2.5 *. a +. noise 1.) x;;
let z = Array.map (fun a -> 2.5 *. a +. noise 8.) x;;

(* 将数组转换为矩阵 *)

let x' = Mat.of_array x 1 50;;
let y' = Mat.of_array y 1 50;;
let z' = Mat.of_array z 1 50;;

(* 绘制图形 *)

let h = Plot.create ~m:1 ~n:2 "plot_01.png" in

Plot.subplot h 0 0;
Plot.set_xlabel h "x";
Plot.set_ylabel h "y (sigma = 1)";
Plot.scatter ~h x' y';

Plot.subplot h 0 1;
Plot.set_xlabel h "x";
Plot.set_ylabel h "z (sigma = 8)";
Plot.scatter ~h x' z';

Plot.output h;;

Stats.corrcoef x y
>- : float = 0.991145445979576656
Stats.corrcoef x z
>- : float = 0.692163016204755288