用向量作逆变换

设向量x长度为n,
其中保存了1到n的正整数的一个排列。
把x看成是在集合{1,2,…,n}上的一个一一变换，
求向量y使得y能够表示上述变换的逆变换。
即任给长度为n的向量z, z[x]表示按照x的次序重新排列z的元素，
而z[x][y]则应该恢复为z。

斐波那契数列计算

设数列x0=0,x1=1，
后续值按如下公式递推计算：

xn=xn−2+xn−1, n=2,3,…

这样的数列叫做斐波那契数列。
希望编写R函数，
输入n，返回计算的xn的值。

穷举所有排列

设向量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,…,n。
从中有放回地抽取m个并记录每次的号码，
穷举m个号码中多少个1，多少个2，…,
多少个n这样的结果。

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

升降连计数

x <- c(
 9, 14, 17,     # (1) 
 3, 16, 18, 20, # (2) 
15, 13,         # 降连
 2,  4,  5, 12, # (3) 
 1,  6, 10,     # (4) 
 7, 19,         # (5) 
 8, 11)         # (6)

plot_runs_prob1 <- function(x){
  plot(x, type="b", lwd=1)
  segs <- matrix(c(
    1,3,  4,7,  10,13,  14,16,  17,18,  19,20),
    ncol=2, byrow=TRUE)
  for(i in seq(nrow(segs))) {
    lines(seq(segs[i,1], segs[i,2]), 
          x[seq(segs[i,1], segs[i,2])], 
          col="red", lwd=2)
  }
}
plot_runs_prob1(x)

n <- 10000
set.seed(102)
x <- sample(n)

x <- c(
  9, 14, 17, # (1)
  3, 3, 9,   # (2)
  3, 3,      # (3)
  2, 20,     # (4)
  17, 17,    # (5)
  12, 10,   # 下降
  1, 5, 5, 16, # (6)
  10, 11)    # (7)

     start end
[1,]     1   3
[2,]     4   6
[3,]     7   8
[4,]     9  10
[5,]    11  12
[6,]    15  18
[7,]    19  20

set.seed(101)
x <- sample(1:5, size = 50, replace=TRUE)

高斯八皇后问题

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

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

任务：

• 编程穷举所有的不同摆法，显示结果。
• 一个摆法通过旋转90度、180度、270度，
  左右对调、上下对调以及这些变换的组合，
  可以转换为另一种摆法。
  按是否可以互相转换将所有摆法分成不同的组，
  显示这些组。
• 用图形显示每一种摆法，并分组显示。
• 将问题推广到9×9的格子摆9个皇后的问题。

最小能量路径

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

下面的例子来自MIT计算思想公开课课件。

考虑一个m×n的矩阵，如：

M=⎛⎝⎜⎜⎜⎜⎜713597361297515936737316912866⎞⎠⎟⎟⎟⎟⎟

考虑从第一行某个元素出发到最后一行某个元素的路径，
每条路径从第一行的某个元素出发，
然后进入下一行的正下方或正下方的左侧元素或右侧元素，
直到进入最后一行的元素为止，
每条路径有m个元素，
计算路径上元素总和。
将M的(i,j)元素看作是该位置的“能量”，
试图找到从顶层（第一行）到底层（最后一行）的“最小能量”路径。
比如，上面的例子从第一行出发，
依次经过的列号为

6,5,4,3,2,

最小能量值为13。

智者千虑必有一失

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

圆桌夫妇座位问题

在一张圆桌上用餐时，
n对夫妇随机入座。
计算没有任何一位妻子和她的丈夫相邻的概率。
通过推导可得此概率为

pn=1+∑k=1n(−1)kCkn(2n−k−1)!2k(2n−1)!.(1)

例如p2=1/3, p3=4/15=0.2667,
p4=0.2952, p5=0.3101。

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

城市间最短路径

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

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

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

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

Daubechies小波函数计算

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

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

ϕ(x)=2‾√∑khkϕ(2x−k)

一种特殊的尺度函数是只在有限区间上非零的, 叫做紧支集的。
紧支集尺度函数可以在给定{hk}后用以下迭代公式生成：

η0(x)=ηn+1(x)=I[−0.5,0.5](x)2‾√∑k=02N−1hkηn(2x−k)

其中N是正整数, N=2时h0=0.482962913145,
h1=0.836516303738, h2=0.224143868042, h3=−0.129409522551。
已知ϕ(x)的支集(不为零的区间)为[0,2N−1],
ηn(x)的支集包含于[−0.5,2N−1]中。

任务：
编写计算ϕ(x)的R程序，
通过20次迭代计算，
输出ϕ(x)在[0,2N−1]区间的256个等间隔点上的函数值并作图。
在迭代过程中，
应不重新结算函数格子点的值，仅计算新加入的格子点的值。
不要利用递归函数，
递归函数需要每次重新计算已经计算过的函数值。

房间加热温度变化

某个房屋带有天花板，天花板与屋顶之间有一定的空间。
用壁炉保持房间温度。
为了研究房间内与天花板上方的温度变化，
建立了如下的微分方程组：

dx1dt=dx2dt=0.35(−9.7sin(t+3)π12+8.3−x1(t))+0.46(x2(t)−x1(t))+11.10.28(−9.7sin(t+3)π12+8.3−x2(t))+0.46(x1(t)−x2(t))

其中x1(t)是房间在t时刻的温度（单位：摄氏度），
x2(t)是天花板上方在t时刻的温度，
t是单位为小时的时间。
设t=0时x1(t)=x2(t)=4。

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

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

线性回归实例

1. 用readr包的read_csv函数读入class.csv为tibble类型的数据框，
   命名为da。
2. 将其中的weight列取出，保存为变量y。
3. 将数据中的age和height两列取出，并在左侧添加一列全是1的列，
   构成一个3列的矩阵，保存为变量X。
4. 用如下的公式计算以体重为因变量，以年龄和身高为自变量的线性回归模型系数：
   β̂ =(XTX)−1XTy

   其中y是包含了因变量观测值的列向量。

5. 用如下公式计算残差的方差估计：
   e=y−Xβ̂ ,σ̂ 2=1n−3∑i=1nê 2i.

   其中n是样本量，ei是向量e的第i个分量。

6. 显示回归系数和σ̂ 2的值。
7. 运行如下的R回归程序，与自己计算的结果比较：

lm(weight ~ age + height, data=da) |> summary()

核回归与核密度估计

考虑核回归问题。核回归是非参数回归的一种, 假设变量Y与变量X之间的关系为：

Y=f(X)+ε

其中函数f未知。
观测到X和Y的一组样本Xi,Yi, i=1,…,n后, 对f的一种估计为：

f̂ (x)=∑i=1nK(x−Xih)Yi∑i=1nK(x−Xih)

其中h>0称为窗宽，
窗宽越大，
得到的密度估计越平滑。
K叫做核函数, 一般是一个非负的偶函数, 原点处的函数值最大, 在两侧迅速趋于零。
例如正态密度函数, 或所谓双三次函数核：

K(x)={(1−|x|3)30|x|≤1其它

与核回归类似，可以用核平滑方法估计总体分布密度。
设样本为Y1,Y2,…,Yn, 密度估计公式为

f̂ (x)=1n∑i=1n1hK(x−Yih)

其中K(x)满足∫∞−∞K(x)dx=1。
h是窗宽，窗宽越大，
估计的曲线越光滑。
h的一种建议公式为h=1.06Sn−1/5,
S为样本标准差。

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

二维随机模拟积分

设二元函数f(x,y)定义如下

=f(x,y)exp{−45(x+0.4)2−60(y−0.5)2}+0.5exp{−90(x−0.5)2−45(y+0.1)4}

(注意其中有一个4次方。)
求如下二重定积分

I=∫1−1∫1−1f(x,y)dxdy

f(x,y)有两个分别以(−0.4,0.5)和(0.5,−0.1)为中心的峰，
对积分有贡献的区域主要集中在(−0.4,0.5)和(0.5,−0.1)附近，
在其他地方函数值很小，对积分贡献很小。
可以用随机模拟方法估计I的值。

第一种估计方法是平均值法。
设Xi,i=1,2,…,N是均匀分布U(−1,1)的随机数，
Yi,i=1,2,…,N也是均匀分布U(−1,1)的随机数，
两者独立，可估计I为

Π1=4N∑i=1Nf(Xi,Yi).

第二种估计方法是重要抽样法。
设正态分布N(μ,σ2)的密度记为p(x;μ,σ2)，
令

=g(x,y)0.5358984p(x;−0.4,90−1)p(y;0.5,120−1)+0.4641016p(x;0.5,180−1)p(y;−0.1,20−1),−∞<x<∞,−∞<y<∞,

这是一个二元随机向量的密度，
是两个二元正态密度的混合分布。
设Ki,i=1,2,…,N是取值于{1,2}的随机数，
P(Ki=1)=0.5358984,
P(Ki=2)=1−P(Ki=1)。
当Ki=1时，取Xi为N(−0.4,90−1)随机数，
Yi为N(0.5,120−1)随机数；
当Ki=2时，取Xi为N(0.5,180−1)随机数，
Yi为N(−0.1,20−1)随机数，
这样得到的(Xi,Yi),i=1,2,…,N是g(x,y)的随机数。
令

Π2=1N∑i=1nf(Xi,Yi)g(Xi,Yi),

称为I的重要抽样法估计。

潜周期估计

设时间序列{yt}有如下模型:

yt=∑k=1mAkcos(λkt+ϕk)+xt, t=1,2,…

其中xt为线性平稳时间序列，
λk∈(0,π), k=1,2,…,m。
这样的模型称为潜周期模型。
如果有{yt}的一组样本y1,y2,…,yn,
可以定义周期图函数

P(ω)=12πn∣∣∣∣∑t=1nyte−itω∣∣∣∣2, ω∈[0,π].

这里ω是角频率。
对于潜周期数据，
在λj的对应位置P(ω)会有尖峰，
而且当n→∞时尖峰高度趋于无穷。
下面是一个样例图形。

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

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

set.seed(1); n <- 500; tt <- seq(n)
m <- 2; lam <- 2*pi/c(4, 7); A <- c(1, 1.2)
y <- A[1]*cos(lam[1]*tt) + A[2]*cos(lam[2]*tt) + rnorm(n)

编写R程序:

1. 编写计算P(ω)的函数，
   输入y=(y1,y2,…,yn)T
   和ω=(ω1,ω2,…,ωs)T,
   输出(P(ω1),P(ω2),…,P(ωs))。
2. 对输入的时间序列样本y1,y2,…,yn，
   编写函数用以上描述的算法估计m和{λj,j=1,2,…,m}。
3. 用上述模拟数据测试编写的算法程序。
4. 进一步地，用R函数fft()计算hj=P(πj/n),j=1,2,…,n。
5. 把整个算法用Rcpp和C++程序实现。

ARMA(1,1)模型估计

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

Xt=aXt−1+et+bet−1,t=1,2,…,T

其中0<|a|<1, 0<|b|<1，
−a≠b,
{et}独立同N(0,σ2e)分布。

1. 给定观测x1,x2,…,xT，
   写出在x0=0且e0=0条件下的条件对数似然函数。
2. 编写R程序，
   用条件最大似然估计方法估计参数a,b,σ2e。
   可使用数值优化程序如nlm()或optim()。
3. 对a=0.5, b=0.7, σe=1，
   T=100, 模拟生成xt的样本并重复模拟B=1000次，
   据此评估â ,b̂ ,σ̂ e的估计精度。
   模拟可使用R的arima.sim()函数。

若x0,e0已知，从x1,…,xT可递推地计算

et=xt−axt−1−bet−1, t=1,2,…,T

记θ=(a,b,σe,x0,e0)T，
有Xt|xt−1,…,x1,θ服从
N(axt−1+bet−1,σ2e)条件分布。
于是用联合密度的乘积公式可得

f(x1,…,xT|θ)====f(x1|θ)f(x2|x1,θ)…f(xT|XT−1,…,x1,θ)∏t=1Tdnorm(xt,axt−1+bet−1,σe)∏t=1T(2π)−1/2σ−1eexp{−121σ2e(xt−axt−1−bet−1)2}(2π)−T/2σ−Teexp{−121σ2e∑t=1Te2t}

其中dnorm(x,μ,σ)表示N(μ,σ2)的分布密度。

取x0=0,e0=0，给定a,b后递推计算{et}序列，
对数似然函数为

l(a,b,σe)=−Tlnσe−121σ2e∑t=1Te2t

注意其中的{et}也与待估参数a,b有关。

可以用数值优化算法估计求如上的问题的最大值点。
显然a,b的最大值点不依赖于σe的取值，
所以可以先求∑Tt=1e2t的最小值点。
这基本是一个最小二乘估计问题，
但是et−1也依赖于a,b的值所以不能直接用线性最小二乘求解。

VAR模型平稳性

称k元时间序列rt服从一个VAR(p)模型，
如果

rt=ϕ0+Φ1rt−1+⋯+Φprt−p+at(58.1)

其中ϕ0和{at}同VAR(1)的规定，
Φj是k阶方阵（k=1,2,…,p）。
利用向后推移算子（滞后算子）B可以将模型写成

(I−Φ1B−⋯−ΦpBp)rt=ϕ0+at

记

P(z)=I−Φ1z−⋯−Φpzp(58.2)

这是一个从复数z到k阶方阵P(z)的变换，
P(z)的每个元素是关于z的阶数不超过p的多项式。
称一元多项式函数det(P(z))（或记为|P(z)|）为模型的（逆序）特征多项式。

如果|P(z)|≠0, ∀|z|<1（z为复数），
则模型(58.1)是平稳的。

编写R函数，
输入三维（三个下标）的系数矩阵数组arrcoef,
返回特征多项式的所有复根。
arrcoef[,,1]保存Φ1,
arrcoef[,,2]保存Φ2,
等等。

作为例子，
考虑如下的三个模型，对每个计算特征多项式的复根并判断是否平稳：

(r1tr2t)=(0.2−0.60.31.1)(r1,t−1r2,t−1)+(a1ta2t)

(x1tx2t)=(0.5−0.25−1.00.5)(x1,t−1x2,t−1)+(a1ta2t)

rt=⎛⎝⎜⎜0.390.350.490.100.340.240.050.470.24⎞⎠⎟⎟rt−1+⎛⎝⎜⎜0.06−0.19−0.310.11−0.18−0.130.02−0.010.09⎞⎠⎟⎟rt−2+at

贮存可靠性评估

设某种设备的贮存寿命为随机变量X,
服从指数分布Exp(θ),
EX=θ>0。
假设有n台此种设备分别在时间t1,t2,…,tn进行了试验，
第i台试验成功用δi=1表示，
试验失效用δi=0表示。把n次试验的结果写成

Zn=⎛⎝⎜⎜⎜⎜t1t2⋮tnδ1δ2⋮δn⎞⎠⎟⎟⎟⎟(1)

这里δ1,δ2,…,δn认为是随机变量，
t1,t2,…,tn是固定的时间。
δi服从两点分布B(1,exp(−ti/θ))。
Zn取某个特定组合的概率为

P(Zn)=∏i=1n{[exp(−ti/θ)]δi[1−exp(−ti/θ)]1−δi}.(2)

要利用观测数据Zn估计θ的置信度为1−α的置信下限，
可以用陈家鼎《生存分析与可靠性》中的样本空间排序法。
注意到Zn中每个δi可以取1或0，
所以Zn的所有不同取值共有2n个，
记这些所有不同取值为，
在中定义如下的序：
设Z′n与Zn都是中的试验结果，
称Z′n不次于Zn,
并记作Z′n≳Zn，
如果如下两个条件之一成立：

1. ∑ni=1δ′i>∑ni=1δi;
2. ∑ni=1δ′i=∑ni=1δi,
   但∑ni=1δ′iti≥∑ni=1δiti.

记

G(θ)=∑Z′n≳ZnPθ(Z′n)=∑Z′n≳Zn∏i=1n{[exp(−ti/θ)]δ′i[1−exp(−ti/θ)]1−δ′i}.(3)

这是θ的严格单调增函数，
求解如下的方程

G(θ)−α=0(4)

得到解θ⎯⎯是θ
的置信度为1−α的置信下限。

当所有n次试验都失效时，
恒有G(θ)=1,
方程(4)无解，
取θ⎯⎯=0。

当n次试验都没有失效，
即失效数n−∑ni=1δi=0时,
方程为

exp(−1θ∑i=1nti)−α=0,

解得

θ⎯⎯=∑ni=1tiln1α.

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

针对n不太大的情形，
可以用精确的公式(3)计算G(θ)
并用(4)求解θ⎯⎯，
编写这个问题的纯R程序版本，
输入一组Zn值后（包括试验时间和试验结果），
输出用(3)和(4)求解θ⎯⎯得到的置信下限θ⎯⎯的值。
在文件store-reliab-data.csv中已经生成了10组模拟数据，
对这10组模拟数据计算相应的θ⎯⎯的值。
在此数据文件中，testid相同的行属于同一次试验的不同设备的时间和结果。

这个程序中用到比较多的循环，
考虑把主要计算部分用Rcpp包和C++代码实现，
看能够提高效率多少倍。

另一种计算G(θ)的方法是用随机模拟方法估计
G(θ)的值然后求解(4)。
随机模拟方法计算简单而且不受失效个数多少的影响，
但是有随机误差，在求解(4)要考虑到随机误差的影响。
随机模拟方法如下。
取模拟次数N,
对给定θ，
为了计算G(θ)，
模拟生成N组独立的寿命
(X(i)1,X(i)2,…,X(i)n),
i=1,2,…,N，
其中每个X(i)j服从Exp(θ)分布，
各分量相互独立。
计算δ(i)j=I{X(i)j>tj},
记

Z(i)n=⎛⎝⎜⎜⎜⎜⎜t1t2⋮tnδ(i)1δ(i)2⋮δ(i)n⎞⎠⎟⎟⎟⎟⎟

用

1N∑i=1NI{Z(i)n≳Zn}

来估计G(θ)。

小题分题型分数汇总

考虑中学某科的一次考试，
考卷各小题小题情况汇总在如下的用逗号符分隔的文本文件subscore-subtype.csv中:

序号,题型
1,选择题
2,选择题
3,选择题
4,选择题
5,选择题
6,选择题
7,选择题
8,选择题
9,选择题
10,选择题
11,简答题
12,简答题
13,填空题
14,简答题
15,简答题
16,简答题
17,简答题
18,写作

文件的中文编码为GB18030。
读入此数据为R数据框dm。
两列都读为字符型。

结果显示如下:

设部分学生的小题分录入到如下的用逗号分隔的文本文件
subscore-subscore.csv中:

学号,Y1,Y2,Y3,Y4,Y5,Y6,Y7,Y8,Y9,Y10,Y11,Y12,Y13,Y14,Y15,Y16,Y17,Y18
1138010104,3,3,3,3,0,3,7.5,2.5,3.5,6,5,4,2.5,5.5,0,4,2,45.5
1138010108,3,0,3,0,3,0,6,3,4,4,2,5,2,5,6,4,2,48.5
1138010114,3,3,3,3,0,3,5.5,3.5,4,6,2,5,2,5.5,6,1,2,44.5
1138010128,3,3,3,0,0,0,8.5,3.5,2.5,4,5,5,0,4.5,6,4,2,45
1138010126,3,3,3,3,3,3,7,0,4,6,5,5,2.5,4.5,6,4.5,2,44

文件的中文编码为GB18030。
读入小题分数据为R数据框ds，
分数读为数值型，
学号读为字符型。

结果显示如下:

在数据框dm中有每个小题的题型信息，
在数据框ds中有每个学生的每个小题的分数。
从这两个数据框，
汇总计算每个学生的题型分，
即每个学生选择题共考多少分，简答题共考多少分，
等等。

要注意的是，
最终的程序应该写成一个函数，
其中的计算不依赖于具体的题型名称、小题个数、小题与题型如何对应，
只要输入dm和ds两个数据框就可以进行统计。

如果没有这样的通用性要求，
这个问题就可以这样简单解决：

resm <- tibble(
  '学号'=ds[['学号']],
  '选择题'= rowSums(ds[, paste('Y', 1:10, sep='')]),
  '简答题'=rowSums(ds[,paste('Y', c(11,12,14:17), sep='' )]),
  '填空题'=ds[['Y13']],
  '作文'=ds[['Y18']]
)
resm |>
  arrange(`学号`) |>
  knitr::kable()

但是，
我们要求最后的函数在允许修改小题类型、试卷中小题个数、各小题的类型、学生分数数据的条件下也能正常工作。

类别编号重排

设聚类分析对若干个样本点聚成了3类，
结果如下：

dc0 <- tibble::tibble(
  obs = 1:10,
  g = c(3, 1, 3, 1, 2, 2, 1, 3, 3, 1)
)
knitr::kable(dc0)

这些类别的次序并没有意义，
重新聚类可能会将类别标签打乱重排但类别组成不变。
所以，
需要人为地指定合适的类别标签。

设对每个类进行概括统计得到了有代表性的统计量，结果如下：

dcstat <- tibble::tibble(
  g = 1:3,
  stat=c(2.3, 1.1, 3.0)
)
knitr::kable(dcstat)

按照统计量从大到小的次序重排3个组，
得到新的分组：

将原始观测添加新的分组标签ng，
并按ng排列，
结果为：

输入dc0和dcstat数据框，
如何用程序解决上述问题？

用R语言下载处理《红楼梦》htm文件

网上许多资源是html格式的文本文件。
比如，
《红楼梦》在许多网站可以浏览，
是在浏览器中按章节浏览。
我们希望将其下载到本地，
并转换为txt格式，
在手机或者电纸书阅读器中阅读。

这个任务涉及到R的文件访问，
字符型连接，
正则表达式，
中文编码问题。

设下载网站主页是
http://www.xiexingcun.net/hongloumeng/index.html，
各个章节的文件名是01.htm到99.htm，
100.htm到120.htm。

已下载的文件在如下链接中：hongloumeng.zip

将这些文件转换成txt格式，
然后合并成一个txt文件。

要求：
每一段仅用一行，中间不换行；
不同段落之间用换行分开；
每一回目开始有回目标题，
标题前后空行，标题格式为“第xxX回 标题内容”。