设数列 \(x_0=0, x_1=1\) ， 后续值按如下公式递推计算： x_n = x_{n-2} + x_{n-1}, \ n=2,3,\dots 这样的数列叫做斐波那契数列。 希望编写R函数， 输入 \(n\) ，返回计算的 \(x_n\) 的值。

设向量 x 的各个元素为某个集合的元素。 想要列出 x 的元素的所有不同排列。 比如，如果 x = 1:3 , 所有排列为

1 2 3
1 3 2
2 1 3
2 3 1
3 1 2
3 2 1

共 \(3!=6\) 种不同的排列。

设有 \(n\) 个编号卡片， 分别有号码 \(1, 2, \dots, n\) 。 从中有放回地抽取 \(m\) 个并记录每次的号码， 穷举 \(m\) 个号码中多少个1，多少个2， \(\dots\) , 多少个 \(n\) 这样的结果。

例如，有3个编号卡片， 随机有放回第抽取2次。 用 \((x_1, x_2, x_3)\) 表示每一种个数组合， \(x_1\) 表示2次抽取中号码1的个数， \(x_2\) 表示2次抽取中号码2的个数， \(x_3\) 表示2次抽取中号码3的个数。 问题就是列出 \((x_1, x_2, x_3)\) 的所有不同值。

考虑如下的数据：

如果输入的序列中有重复数值， 将一个升连定义为：序列中连续的一段， 长度至少有两个数， 后一个数大于等于前一个数， 两个升连之间是下降的。

国际象棋的棋盘是 8×8 的格子。 国际象棋手马克斯·贝瑟尔于 1848 年提出如下问题： 在 8×8 格的国际象棋上摆放八个皇后， 使其不能互相攻击， 即任意两个皇后都不能处于同一行、同一列或同一斜线上， 问有多少种摆法。

高斯认为有 76 种方案。 1854 年在柏林的象棋杂志上不同的作者发表了 40 种不同的解， 后来有人用图论的方法解出 92 种结果。

动态规划是最优化问题的一种算法。 最优化涉及最优决策问题， 此决策问题经常可以分解为若干步骤， 要解决许多子问题， 而这些子问题往往是重复的， 记住已经解决的子问题， 会使得决策问题大大简化， 这是动态规划问题的主要思想。

成语说：“智者千虑，必有一失；愚者千虑，必有一得”。 设智者作判断的准确率为 \(p_1 = 0.99\) , 愚者作判断的准确率为 \(p_2=0.01\) ， 计算智者做1000次独立的判断至少犯一次错误的概率， 与愚者做1000次独立判断至少对一次的概率。

在一张圆桌上用餐时， \(n\) 对夫妇随机入座。 计算没有任何一位妻子和她的丈夫相邻的概率。 通过推导可得此概率为 p_n = 1 + \sum_{k=1}^{n} (-1)^k C_{n}^k \frac{(2n-k-1)! 2^k}{(2n-1)!}. \qquad(1) 例如 \(p_2 = 1/3\) , \(p_3 = 4/15 = 0.2667\) , \(p_4 = 0.2952\) , \(p_5 = 0.3101\) 。

分别用上面的理论公式以及直接穷举验证的方法， 对 \(n=2, 3, 4, 5\) 的情形进行验证。

假设有 \(n\) 个城市，编号为 \(1,2,\dots,n\) 。 已知其中的部分城市之间有高速公路连通， 每对连通城市记为 \((F_i, T_i)\) , 其中 \(F_i, T_i \in \{1,2,\dots,n \}\) 且 \(F_i < T_i\) , \(i=1,2,\dots,m\) 。 除了这些直接连通的城市以外， 其它的任意两个城市只能途经别的城市连通， 或者根本不能靠高速公路连通。 用一个 \(m\times 2\) 的R矩阵 \(M\) 可以输入这些连通情况， 矩阵的每行是一对 \((F_i, T_i)\) 值。

要求编写一个R函数， 输入直接连通情况 \(M\) 后， 输出一个 \(n \times n\) 矩阵 \(R\) , R[i,i]=0, R[i,j]=1表示直接相连， R[i,j]= \(k\) ( \(k \geq 2\) )表示城市i与城市j至少需要经过 \(k\) 段高速公路连通， R[i,j]=Inf表示城市i与城市j不能靠高速公路连通。 \(R\) 的元素值仅考虑途经的高速公路段数而不考虑具体里程。 如果从一个城市通过高速公路移动到直接相连城市叫做移动一步， \(R\) 的 \((i,j)\) 元素是从第 \(i\) 城市通过高速公路到第 \(j\) 城市需要移动的步数。

这个问题也可以作为“相识”问题的模型。 设 \(n\) 个人中有些人是直接相识的， 如果两个不相识的人想认识， 假设必须经过相识的人引荐， 问最少需要多少个引荐人。

本问题必须使用循环， 很难向量化，属于R比较不擅长的问题。 可以考虑使用Rcpp把程序用C++语言实现， 可以大大加速。 如果这个程序仅需要执行不多的次数， 用R就足够了。

这个例子主要展示了不用每次计算函数值而是尽可能从已经计算并储存的函数值中查找的技巧。 程序中用了比较多的循环， 如果有需要， 可考虑用Rcpp转换成C++代码以提高效率。

小波是重要数学工具, 在图像处理、信号处理等方面有广泛应用。 小波中一个重要的函数叫做尺度函数(scale function), 它满足所谓双尺度方程： \phi (x) = \sqrt{2} \sum_k {h_k \phi (2x - k)}

一种特殊的尺度函数是只在有限区间上非零的, 叫做紧支集的。 紧支集尺度函数可以在给定 \(\{h_k\}\) 后用以下迭代公式生成： \begin{aligned} \eta _0 (x) =& I_{[ - 0.5,0.5]} (x) \\ \eta _{n + 1} (x) =& \sqrt 2 \sum\limits_{k = 0}^{2N - 1} {h_k \eta _n (2x - k)} \end{aligned} 其中 \(N\) 是正整数, \(N\) =2时 \(h_0\) =0.482962913145, \(h_1\) =0.836516303738, \(h_2\) =0.224143868042, \(h_3=-0.129409522551\) 。 已知 \(\phi(x)\) 的支集(不为零的区间)为 \([0,2N-1]\) , \(\eta_n(x)\) 的支集包含于 \([-0.5, 2N-1]\) 中。

某个房屋带有天花板，天花板与屋顶之间有一定的空间。 用壁炉保持房间温度。 为了研究房间内与天花板上方的温度变化， 建立了如下的微分方程组： \begin{aligned} \frac{d x_1}{dt} =& 0.35 \left( -9.7 \sin \frac{(t+3)\pi}{12} + 8.3 - x_1(t) \right) + 0.46( x_2(t) - x_1(t) ) + 11.1 \\ \frac{d x_2}{dt} =& 0.28 \left( -9.7 \sin \frac{(t+3)\pi}{12} + 8.3 - x_2(t) \right) + 0.46( x_1(t) - x_2(t) ) \end{aligned} 其中 \(x_1(t)\) 是房间在 \(t\) 时刻的温度（单位：摄氏度）， \(x_2(t)\) 是天花板上方在 \(t\) 时刻的温度， \(t\) 是单位为小时的时间。 设 \(t=0\) 时 \(x_1(t)=x_2(t)=4\) 。

（1） 用每一秒钟重新计算的方法计算24小时内的房间温度与天花板上温度， 绘制求解出的这两个温度的时间序列曲线图。

（2）计算7天的温度，作图，查看周期性。 当温度循环变化差距小于0.1度时认为开始周期变化了。

考虑核回归问题。核回归是非参数回归的一种, 假设变量 \(Y\) 与变量 \(X\) 之间的关系为： Y = f(X) + \varepsilon 其中函数 \(f\) 未知。 观测到X和Y的一组样本 \(X_i, Y_i\) , \(i\) =1,…, \(n\) 后, 对 \(f\) 的一种估计为： \hat f(x) = \frac{\sum\limits_{i = 1}^n {K\left( {\frac{{x - X_i }}{h}} \right)Y_i } } {\sum\limits_{i = 1}^n {K\left( {\frac{{x - X_i }}{h}} \right)} } 其中 \(h>0\) 称为窗宽， 窗宽越大， 得到的密度估计越平滑。 \(K\) 叫做核函数, 一般是一个非负的偶函数, 原点处的函数值最大, 在两侧迅速趋于零。 例如正态密度函数, 或所谓双三次函数核： K(x) = \left\{ {\begin{array}{*{20}c} {\left( {1 - \left| x \right|^3 } \right)^3 } & {\left| x \right| \leq 1} \\ 0 & \mbox{其它} \\ \end{array}} \right.

与核回归类似，可以用核平滑方法估计总体分布密度。 设样本为 \(Y_1, Y_2, \dots, Y_n\) , 密度估计公式为 \hat f(x) = \frac{1}{n} \sum\limits_{i = 1}^n \frac{1}{h} {K\left( {\frac{{x - Y_i }}{h}} \right) } 其中 \(K(x)\) 满足 \(\int_{-\infty}^\infty K(x) \,dx = 1\) 。 \(h\) 是窗宽，窗宽越大， 估计的曲线越光滑。 \(h\) 的一种建议公式为 \(h = 1.06 S n^{-1/5}\) , \(S\) 为样本标准差。

对以上两个问题进行编程，其中窗宽 \(h\) 由用户输入。

设二元函数 \(f(x,y)\) 定义如下 \begin{aligned} & f(x,y) \\ =& \exp\left\{ -45(x + 0.4)^2 - 60(y-0.5)^2 \right\} \\ & + 0.5 \exp\left\{ -90(x-0.5)^2 - 45(y+0.1)^4 \right\} \end{aligned} (注意其中有一个4次方。) 求如下二重定积分 I = \int_{-1}^1 \int_{-1}^1 f(x,y) \,dx\,dy \(f(x,y)\) 有两个分别以 \((-0.4, 0.5)\) 和 \((0.5, -0.1)\) 为中心的峰， 对积分有贡献的区域主要集中在 \((-0.4,0.5)\) 和 \((0.5, -0.1)\) 附近， 在其他地方函数值很小，对积分贡献很小。 可以用随机模拟方法估计 \(I\) 的值。

第一种估计方法是平均值法。 设 \(X_i, i=1,2,\dots, N\) 是均匀分布 \(U(-1,1)\) 的随机数， \(Y_i, i=1,2,\dots, N\) 也是均匀分布 \(U(-1,1)\) 的随机数， 两者独立，可估计 \(I\) 为 \hat I_1 = \frac{4}{N} \sum_{i=1}^N f(X_i, Y_i).

第二种估计方法是重要抽样法。 设正态分布 \(\text{N}(\mu, \sigma^2)\) 的密度记为 \(p(x; \mu, \sigma^2)\) ， \begin{aligned} & g(x,y) \\ =& 0.5358984 p(x;-0.4,90^{-1}) p(y;0.5,120^{-1}) \\ & + 0.4641016 p(x;0.5,180^{-1}) p(y;-0.1,20^{-1}), \\ & -\infty < x < \infty, -\infty < y < \infty, \end{aligned} 这是一个二元随机向量的密度， 是两个二元正态密度的混合分布。 设 \(K_i, i=1,2,\dots,N\) 是取值于 \(\{1,2\}\) 的随机数， \(P(K_i = 1) = 0.5358984\) , \(P(K_i = 2) = 1 - P(K_i = 1)\) 。 当 \(K_i=1\) 时，取 \(X_i\) 为N( \(-0.4,90^{-1}\) )随机数， \(Y_i\) 为N( \(0.5,120^{-1}\) )随机数； 当 \(K_i=2\) 时，取 \(X_i\) 为N( \(0.5,180^{-1}\) )随机数， \(Y_i\) 为N( \(-0.1,20^{-1}\) )随机数， 这样得到的 \((X_i, Y_i), i=1,2,\dots,N\) 是 \(g(x,y)\) 的随机数。 \hat I_2 = \frac{1}{N} \sum_{i=1}^n \frac{f(X_i, Y_i)}{g(X_i, Y_i)}, 称为 \(I\) 的重要抽样法估计。

注意尽量用向量化编程。

设时间序列 \(\{ y_t \}\) 有如下模型: y_t = \sum_{k=1}^m A_k \cos(\lambda_k t + \phi_k) + x_t, \ t=1,2,\dots 其中 \(x_t\) 为线性平稳时间序列， \(\lambda_k \in (0,\pi)\) , \(k=1,2,\dots,m\) 。 这样的模型称为潜周期模型。 如果有 \(\{ y_t \}\) 的一组样本 \(y_1, y_2, \dots, y_n\) , 可以定义周期图函数 P(\omega) = \frac{1}{2\pi n} \left| \sum_{t=1}^n y_t e^{-it\omega} \right|^2, \ \omega \in [0, \pi]. 这里 \(\omega\) 是角频率。 对于潜周期数据， 在 \(\lambda_j\) 的对应位置 \(P(\omega)\) 会有尖峰， 而且当 \(n\to\infty\) 时尖峰高度趋于无穷。 下面是一个样例图形。

如下算法可以在 \(n\) 较大时估计 \(m\) 和 \(\{\lambda_k \}\) ： 首先，对 \(\omega_j=\pi j / n\) , \(j=1,2,\dots,n\) 计算 \(h_j = P(\omega_j)\) ， 求 \(\{ h_j, j=1,2,\dots, n \}\) 的3/4分位数记为 \(q\) 。 令 \(C = q n^{0.25}\) ，以 \(C\) 作为分界线， 设 \(\{ h_j \}\) 中大于 \(C\) 的下标 \(j\) 的集合为 \(J\) , 当 \(J\) 非空时，把 \(J\) 中相邻点分入一组， 但是当两个下标的差大于等于 \(n^{0.6}\) 时就把后一个点归入新的一组。 在每组中，以该组的 \(h_j\) 的最大值点对应的角频率 \(j \pi / n\) 作为潜频率 \(\{ \lambda_k \}\) 中的一个的估计。

用如下R程序可以模拟生成一组 \(\{ y_t \}\) 的观测数据:

考虑如下的零均值高斯ARMA(1,1)模型：

若 \(x_0, e_0\) 已知，从 \(x_1, \dots, x_T\) 可递推地计算 e_t = x_t - a x_{t-1} - b e_{t-1}, \ t=1,2,\dots, T 记 \(\boldsymbol{\theta} = (a, b, \sigma_e, x_0, e_0)^T\) ， 有 \(X_t | x_{t-1}, \dots, x_1, \boldsymbol{\theta}\) 服从 N( \(a x_{t-1} + b e_{t-1}, \sigma_e^2\) )条件分布。 于是用联合密度的乘积公式可得 \[\begin{aligned} f(x_1, \dots, x_T | \boldsymbol{\theta}) =& f(x_1 | \boldsymbol{\theta}) f(x_2 | x_1, \boldsymbol{\theta}) \dots f(x_T | X_{T-1}, \dots, x_1, \boldsymbol{\theta}) \\ =& \prod_{t=1}^T \text{dnorm}(x_t, a x_{t-1} + b e_{t-1}, \sigma_e) \\ =& \prod_{t=1}^T (2\pi)^{-1/2} \sigma_e^{-1} \exp\left\{ -\frac12 \frac{1}{\sigma_e^2} (x_t - a x_{t-1} - b e_{t-1})^2 \right\} \\ =& (2\pi)^{-T/2} \sigma_e^{-T} \exp\left\{ -\frac12 \frac{1}{\sigma_e^2} \sum_{t=1}^T e_t^2 \right\} \end{aligned}\] 其中 \(\text{dnorm}(x, \mu, \sigma)\) 表示N( \(\mu, \sigma^2\) )的分布密度。

取 \(x_0=0, e_0=0\) ，给定 \(a, b\) 后递推计算 \(\{ e_t \}\) 序列， 对数似然函数为 l(a, b, \sigma_e) = -T \ln\sigma_e - \frac12 \frac{1}{\sigma_e^2} \sum_{t=1}^T e_t^2 注意其中的 \(\{e_t \}\) 也与待估参数 \(a, b\) 有关。

可以用数值优化算法估计求如上的问题的最大值点。 显然 \(a, b\) 的最大值点不依赖于 \(\sigma_e\) 的取值， 所以可以先求 \(\sum_{t=1}^T e_t^2\) 的最小值点。 这基本是一个最小二乘估计问题， 但是 \(e_{t-1}\) 也依赖于 \(a, b\) 的值所以不能直接用线性最小二乘求解。

称 \(k\) 元时间序列 \(\boldsymbol r_t\) 服从一个VAR( \(p\) )模型， \[\begin{align} \boldsymbol r_t = \boldsymbol\phi_0 + \boldsymbol\Phi_1 \boldsymbol r_{t-1} + \dots + \boldsymbol\Phi_p \boldsymbol r_{t-p} + \boldsymbol a_t \tag{58.1} \end{align}\] 其中 \(\boldsymbol\phi_0\) 和 \(\{ \boldsymbol a_t \}\) 同VAR(1)的规定， \(\boldsymbol\Phi_j\) 是 \(k\) 阶方阵（ \(k=1,2,\dots,p\) ）。 利用向后推移算子（滞后算子） \(B\) 可以将模型写成 \[\begin{aligned} (\boldsymbol I - \boldsymbol\Phi_1 B - \dots - \boldsymbol\Phi_p B^p) \boldsymbol r_t = \boldsymbol\phi_0 + \boldsymbol a_t \end{aligned}\] \[\begin{align} = \boldsymbol I - \boldsymbol\Phi_1 z - \dots - \boldsymbol\Phi_p z^p \tag{58.2} \end{align}\] 这是一个从复数 \(z\) 到 \(k\) 阶方阵 \(P(z)\) 的变换， \(P(z)\) 的每个元素是关于 \(z\) 的阶数不超过 \(p\) 的多项式。 称一元多项式函数 \(\text{det}(P(z))\) （或记为 \(|P(z)|\) ）为模型的（逆序）特征多项式。

设某种设备的贮存寿命为随机变量 \(X\) , 服从指数分布Exp( \(\theta\) ), \(E X = \theta > 0\) 。 假设有 \(n\) 台此种设备分别在时间 \(t_1, t_2, \dots, t_n\) 进行了试验， 第 \(i\) 台试验成功用 \(\delta_i=1\) 表示， 试验失效用 \(\delta_i=0\) 表示。把 \(n\) 次试验的结果写成 Z_n = \left(\begin{array}{cc} t_1 & \delta_1 \\ t_2 & \delta_2 \\ \vdots & \vdots \\ t_n & \delta_n \end{array}\right) \tag{1} 这里 \(\delta_1, \delta_2, \dots, \delta_n\) 认为是随机变量， \(t_1, t_2, \dots, t_n\) 是固定的时间。 \(\delta_i\) 服从两点分布 \(\text{B}(1, \exp(-t_i / \theta))\) 。 \(Z_n\) 取某个特定组合的概率为 P(Z_n) = \prod_{i=1}^n \left\{ \left[ \exp(-t_i / \theta) \right]^{\delta_i} \left[ 1 - \exp(-t_i / \theta) \right]^{1 - \delta_i} \right\}. \tag{2}

要利用观测数据 \(Z_n\) 估计 \(\theta\) 的置信度为 \(1-\alpha\) 的置信下限， 可以用陈家鼎《生存分析与可靠性》中的样本空间排序法。 注意到 \(Z_n\) 中每个 \(\delta_i\) 可以取1或0， 所以 \(Z_n\) 的所有不同取值共有 \(2^n\) 个， 记这些所有不同取值为 \(\mathcal D\) ， 在 \(\mathcal D\) 中定义如下的序： 设 \(Z_n'\) 与 \(Z_n\) 都是 \(\mathcal D\) 中的试验结果， 称 \(Z_n'\) 不次于 \(Z_n\) , 并记作 \(Z_n' \gtrsim Z_n\) ， 如果如下两个条件之一成立：

当所有 \(n\) 次试验都失效时， 恒有 \(G(\theta)=1\) , 方程(4)无解， 取 \(\underline{\theta}=0\) 。

当 \(n\) 次试验都没有失效， 即失效数 \(n - \sum_{i=1}^n \delta_i = 0\) 时, \exp \left( -\frac{1}{\theta} \sum_{i=1}^n t_i \right) - \alpha = 0, \underline{\theta} = \frac{\sum_{i=1}^n t_i}{\ln \frac{1}{\alpha}} .

当失效数为1时，最多仅有 \(n\) 个结果不次于 \(Z_n\) ， 很容易可以计算 \(G(\theta)\) 的值。 在 \(n\) 较大而且 \(Z_n\) 中失效数较多时， \(G(\theta)\) 的求和(3)中项数很多，计算量很大。 当 \(n=10\) 时， \(\mathcal D\) 约有一千项， 当 \(n=20\) 时， \(\mathcal D\) 约有一百万项， 还在可以穷举计算的范围内。 但是，当 \(n \geq 30\) 时， \(\mathcal D\) 就有10亿项， 每计算一次 \(G(\theta)\) 都需要很长时间。