一种广义的The Lady Tasting Tea 实验的贝叶斯网络分析框架

谢钢（Gang Xie）

Office of Research Services and Graduate Studies Charles Sturt University Australia

写在前面

The Lady Tasting Tea 是统计学史上最著名的经典实验之一，通常被认为是现代统计假设检验（statistical hypothesis testing）的起点。近百年来，这一实验几乎出现在所有统计学教材中，并成为解释 p-值 和显著性检验（test of significance）的经典案例。

2024 年，我用贝叶斯网络（Bayesian Network，BN）首次建立了The Lady Tasting Tea 实验的推荐推断分析模型，提出了一种能够计算受试者能力后验概率（posterior probability）的方法。与传统的显著性检验相比，该模型不仅能够计算 p-值，还能够根据实验结果估计受试者具有不同能力水平的可能性。

本文是在前期工作的基础上进一步发展而来的。文章提出了一种基于模拟—EM（Simulation–EM）方法的建模框架，使 The Lady Tasting Tea 实验能够推广到任意偶数个茶杯（N）以及任意数量能力水平（k）的实验设计。随着实验规模不断扩大，模型参数数量将呈指数增长，传统依赖人工计算条件概率表的方法已难以实施，而本文提出的方法能够有效解决这一问题。

本中文译文根据我的英文预印本 （链接： https://www.researchsquare.com/article/rs-8634718/v1 ）整理而成（有ChatGPT的参与翻译）。这个英文预印本就是那个我在https://blog.sciencenet.cn/blog-3503579-1540150.html 里提到的被Nature 的子刊Scientific Reports拒稿的文本。有兴趣的读者可阅读原文看看其质量如何，是否应该被拒稿。正文采用中文学术写作习惯进行整理翻译，数学公式保持原样，R 程序保留原文，以方便读者与原论文相互对照。

摘要

Lady Tasting Tea 实验被广泛认为是现代统计假设检验的起源。在我们此前的研究中，曾利用贝叶斯网络（Bayesian Network，BN）对这一经典实验进行了建模，但当时的模型仅考虑了三个能力水平，且所有条件概率均需人工计算和输入，因此难以推广到更复杂的实验设计。

本文提出一种新的建模方法，使贝叶斯网络能够适用于更加一般化的 The Lady Tasting Tea 实验。本文的主要贡献在于提出了一种模拟样本数据生成方法，可支持任意偶数个茶杯（N）以及任意数量能力水平（k≥2）的实验设计。基于 R 模拟程序，本文建立了一个包含 10 杯茶、5 个能力水平的贝叶斯网络模型，并利用 Netica 内置的 Expectation–Maximisation（EM）算法完成模型参数学习。整个模型包含 21 个节点和 22,511 个条件概率参数，若采用人工方式建立几乎不可行，而 Simulation–EM 方法能够高效完成模型训练。

虽然随着实验规模扩大，模型参数空间仍将呈指数增长，但本文提出的方法为推广 The Lady Tasting Tea 实验提供了一种具有良好可扩展性的建模框架。除传统的 Fisher 显著性检验外，该框架还能够计算受试者能力水平的后验概率，从而为统计推断提供更加丰富和完整的信息。

关键词：

Lady Tasting Tea；Bayesian Network（BN）；Simulation-based inference；Fisher's test of significance；posterior probability

1 引言

The Lady Tasting Tea 实验是统计学发展史上的经典实验，被普遍认为是现代统计假设检验（statistical hypothesis testing）的起点。在 The Design of Experiments 一书中，R. A. Fisher 借助这一实验阐述了统计显著性检验（test of significance）的基本思想。

实验的背景十分简单。一位女士（Dr. Blanche Muriel Bristol）声称，仅凭品尝即可判断一杯奶茶究竟是先加牛奶（milk-first）还是先加茶（tea-first）。为了验证这一说法，Fisher 设计了一项实验：准备八杯奶茶，其中四杯采用 milk-first，另外四杯采用 tea-first，并以随机顺序呈现给受试者。受试者需要判断每一杯奶茶属于哪一种冲泡方式。

在 Fisher 的框架下，首先建立零假设（null hypothesis）：受试者并不具备辨别两种冲泡方式的能力，其判断完全来自随机猜测，即每一杯茶被正确判断的概率均为 50%。

由于实验设计预先固定为 8 杯茶（4 杯 milk-first、4 杯 tea-first），所有可能的实验安排共有 70 种，且每一种出现的概率完全相同。因此，如果受试者能够正确识别全部 8 杯茶，则随机猜测得到这一结果的概率仅为

Fisher 将这一较小的概率解释为反对零假设的有力证据，并据此认为，实验结果支持受试者确实具有辨别两种冲泡方式的能力。

另一方面，如果受试者仅正确识别 8 杯中的 6 杯，则随机猜测产生这一结果的概率为

因此，在单侧检验中，随机猜测得到 6 杯或更多正确判断的概率（即 p-value）为

在这种情况下，Fisher 将认为实验结果不足以拒绝零假设。然而，从统计推断的角度来看，一个更值得回答的问题是：根据实验结果，受试者真实的能力水平究竟是多少？

在我们此前的研究中，曾利用贝叶斯网络（BN）对 The Lady Tasting Tea 实验进行了建模（见图1）。该模型能够根据受试者的表现，在给定能力水平先验分布（prior distribution）的条件下，计算各能力水平对应的后验概率（posterior probability）。

图1

与传统显著性检验不同，贝叶斯网络模型关注的是

而 Fisher 方法计算的 p-value 实际上对应于

两者在概率意义上并不相同，因此能够回答的问题也不同。前者直接对应研究者真正关心的问题，而后者只是衡量在零假设成立时观察到当前数据的可能性。

基于上述工作，本文进一步提出一种Simulation–EM 建模方法，使 The Lady Tasting Tea 实验能够推广到任意偶数个茶杯和任意数量能力水平的实验设计。本文首先利用 R 程序生成满足实验设计要求的大规模模拟样本数据，然后借助 Netica 的 EM 算法完成模型参数学习，从而避免人工构建条件概率表带来的巨大工作量。

本文余下内容安排如下。第 2 章介绍实验概率分布、模拟样本数据生成方法以及贝叶斯网络模型的构建过程；第 3 章介绍所建立的 BN 模型及其应用实例；第 4 章讨论本文方法与传统显著性检验之间的区别，以及该方法的意义与未来研究方向。

2 方法（Methods）

本文提出的一般性的推广 The Lady Tasting Tea 实验建模方法主要包括三个部分：首先分析实验中正确识别冲制方式的茶杯数量的概率分布；其次介绍模拟样本数据的生成方法；最后说明贝叶斯网络模型的构建过程。

2.1 正确识别茶杯数量的概率分布

在 Fisher 的显著性检验框架下，The Lady Tasting Tea 实验关注的问题是：受试者是否具有辨别奶茶冲泡方式的能力。对应的零假设（null hypothesis）认为，受试者完全依赖随机猜测，因此每杯茶被正确判断的概率均为 0.5。

然而，从统计建模的角度来看，我们真正关心的是受试者的能力水平（ability level）。设实验共提供 杯茶，受试者正确判断单杯茶的冲制方式的概率为 ，其中：

• ：完全随机猜测；

• ：始终判断正确；

• ：具有不同程度的辨别能力。

如果每一杯茶的测试结果彼此独立，则正确识别的茶杯数量的统计规律服从二项分布，其概率可表示为论文中的公式

（1）。

这里需要强调的是，Lady Tasting Tea 实验与普通二项实验并不相同。实验要求最终必须判断出相同数量的 milk-first 和 tea-first 茶杯，因此各次测试结果之间并非独立事件，而是存在约束关系。受试者前面的判断结果会影响后续测试结果的可能性，这也是该实验区别于标准二项分布的重要特点。

对于这种具有依赖关系的实验，目前只有在 （完全随机猜测）时能够得到解析表达式，其概率由论文中的公式

（2）给出。

利用该公式，可以方便地计算不同实验设计下随机猜测产生各种实验结果的概率。例如，在经典的 8 杯茶实验中：

• 全部判断正确的概率为 ；

• 正确判断 6 杯的概率为 。

同样，当实验扩展到 10 杯茶时，由于共有

种不同的冲泡顺序，因此各种实验结果的概率仍可依据公式（2）计算。

但是，当受试者能力水平 时，其测试结果的计算的原理公式           目前尚不存在对应的一般解析公式。这也是本文采用模拟方法建立推广模型的主要原因。

2.2 推广实验模拟样本数据的生成

由于 The Lady Tasting Tea 实验的设计规则完全确定，因此可以利用 R 程序模拟整个实验过程，并生成用于训练 BN 模型的大规模样本数据。

每一组模拟数据包含三部分信息：

• 受试者预设能力水平（ability level）；

• 实验中茶杯的实际排列顺序；

• 受试者对应的判断结果。

模拟过程可概括为以下四个步骤。

（1）初始化实验设计

设实验共有

杯茶，其中一半采用 milk-first，另一半采用 tea-first。随后，对所有茶杯进行随机排列，生成一次实验所对应的真实冲泡顺序。

（2）生成受试者判断结果

根据能力水平的不同，模拟过程分为两种情况。

第一种情况：随机猜测（）

受试者对第一杯茶完全随机判断。之后每一步判断都需要考虑前面已经给出的结果，使最终判断中 milk-first 与 tea-first 的数量始终满足实验设计要求，不会出现某一类别数量超过允许范围的情况。

第二种情况：具有辨别能力（）

每一次判断均按照能力水平 决定：

• 以概率 作出正确判断；

• 以概率 作出错误判断。

与此同时，模拟程序仍然保证最终判断结果满足实验设计约束，即两类茶杯数量保持一致。

（3）生成输出数据

模拟结束后，将程序中的数值编码转换为对应的文字标签：

• 1 → milk

• 0 → tea

随后组合形成一组完整模拟样本，包括：

• 能力水平；

• 茶杯排列顺序；

• 全部判断结果。

论文附录给出了 onesample(pp, NN) 函数的完整 R 程序代码。

（4）生成训练数据集

为了训练推广后的 BN 模型，需要重复运行上述模拟过程大量次数。

例如，本文采用如下五个能力水平：0.5，0.6，0.75，0.9，1.0。

利用 R 随机抽取 100,000 个能力水平，并分别运行 onesample() 函数，即可生成包含 100,000 次实验结果的模拟数据集。生成的数据最终保存为    simuData_Tea_Milk.csv 的数据文件。该数据集随后用于 Netica 中 BN 模型参数的学习。

表 1(A) 和 表1(B) 给出了模拟数据集前六行的示例。其中：第一列表示能力水平；后续 10 列记录实验中茶杯的真实排列顺序；最后 10 列记录受试者对各杯茶的判断结果。

2.3 BN 模型的构建

本文所有贝叶斯网络（BN）模型均采用 Netica 软件构建。

BN 本质上是一种联合概率分布模型，其网络结构采用有向无环图（Directed Acyclic Graph，DAG）表示。图中的每个节点对应一个随机变量，节点之间的连线表示变量之间的条件依赖关系。

一个完整的 BN 模型包含两个组成部分：

• 网络结构（Structure）：确定各变量之间的依赖关系；

• 模型参数（Parameters）：利用条件概率表（Conditional Probability Table，CPT）量化各变量之间的概率关系。

在前期研究中，我们依据 Lady Tasting Tea 实验的设计，人工建立了 BN 网络结构，并逐项计算所有 CPT 参数。

对于 图 1 所示模型，共需要填写 3,573 个 CPT 参数。

随着能力水平由 3 个扩展到 5 个，图 2 的 CPT 参数数量进一步增加至 5,615 个。

可以看出，随着实验规模不断扩大，人工建立 CPT 已变得十分困难，不仅工作量巨大，而且容易出现人为错误。

为此，本文提出 Simulation–EM 方法。该方法仍由研究者人工建立 BN 网络结构，但模型参数不再人工计算，而是利用大量模拟样本数据，通过 Netica 内置的 Expectation–Maximisation（EM）算法自动学习完成。

对于推广后的 10 杯茶、5 个能力水平模型，共包含 21 个节点和 22,511 个 CPT 参数。借助 Simulation–EM 方法，整个参数学习过程均可自动完成，大幅提高了建模效率，并使 Lady Tasting Tea 实验能够推广到更复杂的实验设计。

训练 图 3 所用的模拟数据作为补充材料提供。建模时，应确保 CSV 文件中的变量名称与 Netica 模型中的节点名称完全一致，以保证模型能够正确读取数据并完成参数学习。

图2

图3

3 结果（Results）

本文建立了两个 BN 模型。第一个模型对应经典的 8 杯茶实验，但将受试者能力水平由前期研究中的 3 个水平扩展为 5 个水平；第二个模型则采用本文提出的 Simulation–EM 方法，建立了适用于 10 杯茶、5 个能力水平的推广模型。下面分别介绍这两个模型及其应用。

3.1 所建立的 BN 模型

本文首先建立了 图 2 所示的 BN 模型。

该模型共包含 17 个节点。其中：

• 底部 8 个节点表示实验中依次提供的 8 杯茶，每个节点只有 milk-first 和 tea-first 两种状态，两种状态的先验概率均为 50%。

• 中间 8 个节点表示受试者对每一杯茶的判断结果。

• 顶部节点表示受试者能力水平（ability level）的先验分布。

与前期研究不同，本文将能力水平扩展为 5 个等级：0.5（完全随机猜测），0.6，0.75，0.9，1.0（始终判断正确）。

图 2 中顶部节点默认采用均匀先验分布，即五个能力水平具有相同的先验概率。

整个模型共包含 5,615 个条件概率表（CPT）参数。

在 BN 中，每个节点表示一个随机变量，节点之间的连线表示变量之间的条件依赖关系。节点中显示的百分比表示各状态对应的概率，右侧条形图则直观展示了这些概率分布。

图 3 展示了本文提出的一般性的推广 BN 模型。

该模型对应 10 杯茶、5 个能力水平的实验设计。

与 图 2 不同，该模型并不是通过人工计算所有条件概率建立，而是首先利用 R 程序生成 100,000 组模拟实验数据，再利用 Netica 的 Expectation–Maximisation（EM）算法完成参数学习。

整个模型共包含：

• 21 个节点

• 22,511 个 CPT 参数

需要指出的是，由于 图 3 属于基于模拟样本学习得到的模型，因此它与 图 1 和 图 2 的理论模型存在一个细微区别。

图 1 和 图 2 中所有概率均依据实验设计直接计算，因此属于理论上的精确模型；而 图 3 则属于统计学习模型，其参数来源于有限数量的模拟样本，因此能力水平节点的先验分布以及各茶杯节点的初始概率只能近似满足理论值，而不可能完全相同。不过，当模拟样本数量足够大时，这种差异可以忽略，对模型推断结果几乎没有影响。

3.2 BN 模型的应用示例

上述 BN 模型均建立在标准贝叶斯统计框架之上，因此不仅能够计算传统显著性检验中的 p-值，还能够计算受试者能力水平的后验概率。这也是本文模型相较于传统 Fisher 方法最大的优势。

图 2 所示 8 杯茶、5 个能力水平模型的使用方法与前期研究中的三能力水平模型基本一致。

本文首先利用该模型说明 BN 推断过程，然后再介绍推广模型（Figure 3）的应用。

由于 BN 本质上是一个联合概率分布模型，因此既可以计算

也可以计算

前者对应于传统显著性检验所关心的问题，即在给定能力水平条件下观察到实验结果的概率；后者则是贝叶斯统计真正关心的问题，即在已经观察到实验结果之后，不同能力水平成立的后验概率。

图 4 给出了一个典型示例。

图4

实验中 8 杯茶的排列顺序为：TTMMMTTM

即：Tea-first、Tea-first、Milk-first、Milk-first、Milk-first、Tea-first、Tea-first、Milk-first。受试者成功判断了全部 8 杯茶。

在这一条件下，BN 模型根据默认均匀先验分布计算得到受试者能力水平的后验分布，其结果显示在 Figure 4 顶部粉色节点中。

如果固定某一能力水平，例如假定受试者能力水平为 0.75，则模型还可以进一步计算对应实验结果出现的概率，即

由此即可得到传统意义上的 p-值。

表 2 汇总了不同能力水平对应的计算结果。

需要特别指出的是，在 The Lady Tasting Tea 实验中，实验结果不仅与正确判断茶杯数量有关，还与茶杯实际提供顺序有关。

例如，当实验结果均为 全部 8 杯判断正确时，如果茶杯排列顺序由

TTMMMTTM改为TTTTMMMM，则对应的 p-值 将发生变化（能力水平为 0.5 时除外）。这是因为 The Lady Tasting Tea 实验中各判断之间存在依赖关系，而不是普通二项分布中的独立试验。

目前，对于能力水平大于 0.5 的情况尚不存在相应解析公式，因此本文提出的 BN 模型实际上提供了一种数值计算方法，可适用于任意满足实验设计要求的实验情景。

Netica 的 Findings Calibration 功能允许研究者自行设定能力水平的先验分布。

图 5 展示了三种不同的先验分布，包括两种均匀分布和一种非均匀分布。

在相同实验结果下，模型分别计算得到三种不同的后验分布。

这一结果说明：后验概率不仅受到实验数据的影响，同时也受到能力水平先验分布的影响。这正是贝叶斯统计推断的重要特点。

图5

图3 所示推广模型同样可以完成上述各种推断分析。

不过，由于该模型采用模拟数据进行参数学习，因此理论上仍可能出现一种特殊情况，即受试者判断结果违反实验设计约束，例如某一类别茶杯数量超过允许范围。论文将这种情况称为 overguess。

图6(A) 给出了一个 overguess 示例，该情景不符合 Lady Tasting Tea 实验设计，因此不能作为有效实验结果。

图 6(B) 则给出了对应的合法实验情景，可作为实际分析时的替代方案。

在实际应用中，只要保证输入数据满足实验设计要求，推广模型即可完成与经典模型相同的贝叶斯推断分析，同时能够处理更加复杂的实验设计。

4 讨论（Discussion）

本文提出的 BN 模型为什么能够比传统显著性检验更好地回答 The Lady Tasting Tea 实验所关注的问题？回答这一问题，首先需要明确实验真正希望解决的科学问题。

对于 The Lady Tasting Tea 实验而言，我们真正关心的并不是：

这位女士是否真的具有辨别奶茶冲泡方式的能力？

而是：

根据实验结果，这位女士真实的能力水平究竟是多少？

这两个问题看似相近，但在统计推断中的含义却完全不同。

Fisher 的显著性检验关注的是：如果零假设成立，当前实验结果出现的概率有多大；而研究者真正希望了解的是：在已经观察到实验结果之后，各种能力水平成立的可能性分别有多大。

换句话说，传统显著性检验回答的是

而贝叶斯方法关注的是

两者并不能相互替代。

这一差别并不仅仅是数学公式上的不同，而是统计推断思想上的根本区别。

在 Fisher 的方法中，较小的 p-值 通常被解释为反对零假设的有力证据。然而，从逻辑上讲，反对零假设并不等同于支持某一个具体的备择假设。

事实上，Neyman 曾指出，如果在没有考虑统计功效（statistical power）以及明确备择假设的情况下，仅依赖显著性检验作出科学结论，可能导致错误甚至具有误导性的推断。

例如，假设受试者真正的能力水平只有 0.55，即仅略高于随机猜测。如果实验样本量较小，则统计功效可能非常低，即使受试者确实具有一定辨别能力，显著性检验仍然很难拒绝零假设。此时，仅依据 p-值 作出结论，并不能准确反映受试者真实的能力水平。

本文提出的 BN 模型并不试图回答"是否拒绝零假设"这一问题，而是直接根据实验结果估计不同能力水平的后验概率，从而为研究者提供更加完整的推断信息。

从方法论角度来看，BN 模型的另一项优势在于它本质上是一个联合概率分布模型（joint probability distribution model）。

模型中的每个节点对应一个随机变量，节点之间的依赖关系利用 Bayes' Theorem 进行描述，再通过概率链式法则（probability chain rule）组合成为整个系统的联合概率分布。

因此，在同一个模型中，任意变量都可以作为分析目标，其余变量则作为观测条件（findings）。

这种统一的概率框架，使 BN 模型既能够计算传统显著性检验中的

也能够直接计算研究者更加关心的

这也是本文模型能够同时提供 p-value 和后验概率分析结果的理论基础。

有人认为，Bayesian 方法需要人为指定先验分布，因此具有一定的主观性（subjectivity）。本文认为，这种观点值得进一步讨论。

无论采用 Frequentist 方法还是 Bayesian 方法，统计建模都不可避免地包含研究者对模型、数据和假设的选择，因此一定程度的主观判断始终存在。关键并不在于是否完全消除主观性，而在于这些假设是否合理、是否能够根据新的证据不断修正。

从这一意义上说，本文更加倾向于把统计推断理解为一种假设情景分析（what-if analysis）工具，而不是一种能够最终"证明"科学结论的方法。

BN 模型允许研究者在不同先验假设下重复分析实验结果，并比较不同假设对推断结果的影响。这种分析方式更符合科学研究不断修正、不断更新认识的过程。

本文的应用实例表明，所建立的 BN 模型不仅能够完成 Fisher 显著性检验能够完成的分析，而且能够进一步给出受试者能力水平的后验概率，从而为 The Lady Tasting Tea 实验提供更加全面、更加符合逻辑的统计解释。

对于经典的 8 杯茶实验，BN 模型能够完整再现实验设计，并计算各种实验结果对应的概率；对于推广后的实验设计，本文提出的 Simulation–EM 方法则解决了模型参数数量迅速增长所带来的建模困难，使 BN 建模能够扩展到任意偶数个茶杯以及任意数量能力水平的实验。

虽然本文采用的是基于模拟样本和 EM 算法的参数学习方法，但这一思路具有良好的可扩展性。未来还可以进一步利用标准 Bayesian 建模方法（例如基于 MCMC 的 BUGS）开展相关研究，对推广后的 Lady Tasting Tea 实验进行更加深入的分析。

总体而言，本文提出的 Simulation–EM 方法不仅解决了推广 The Lady Tasting Tea 实验中的建模问题，也为利用贝叶斯网络研究类似实验设计提供了一种具有普遍意义的方法框架。

受资助信息披露Funding

本研究未获得任何外部经费资助。

致谢及声明Acknowledgment

本文的早期版本曾在 BNMA 16th Annual Bayesian Network Modelling Association Conference（Melbourne, Australia，October 2025）以会议报告的形式进行交流。作者感谢会议期间与会者提出的宝贵意见和建议。

参考文献（略）

附录Appendix:

# Generation of random samples of the data that follow the design of the Lady Tasting Tea experiment

#

   onesample <- function(pp, NN) {

#  Step 1: generation of the 2*NN cups with which NN cups of each kind

#  (milk or tea adding into the cup first) and the 2N cups will be served in a random order.

#  denote cups[i] = 1 if milk first, cups[i] = 0 if tea first.

  PrepareCups = c(rep(1,NN), rep(0,NN))

 cups = NULL

  cups = sample(PrepareCups, 2*NN, replace=F)  # this is a set of 2*NN cups arranged in a random order

# Step 2: generation of the corresponding the tester's assessment results

   testcups = 1:(2*NN)

  for (i in 1: (2*NN)) {

   if( pp == 0.5) {

    if (i == 1)  testcups[i] =  sample(c(0,1), 1)               

#

    if (i > 1 && i <= NN) {

      for ( j in 2:NN) {

        m = sum(testcups[1:(j-1)])

       p1 = (NN-m)/(2*NN-j+1); p0 = 1-p1

       testcups[j] = sample(c(0,1),1, prob=c(p0,p1))          }

                    }                   # end of if (i > 1 && i <= NN)

    if (i > NN) {

    for (k in (NN+1):(2*NN))  {

          m1 = sum(testcups[1:(k-1)])

          m0 = (k-1) - m1

     if(m1 >= NN) testcups[k] = 0

     if(m0 >= NN) testcups[k] = 1

     if(m1 < NN && m0 < NN) {

              m = sum(testcups[1:(k-1)])

       p1 = (NN-m)/(2*NN-k+1); p0 = 1-p1

       testcups[k] = sample(c(0,1),1, prob=c(p0,p1))    } 

         }                                 # end of  for (k in (NN+1):(2*NN))

             }                             # end of if (i > NN)

                         }               # end of if (pp == 0.5)

  else   {

        if (i <= NN) testcups[i] =  sample(c(cups[i],(1-cups[i])), 1, prob=c(pp,(1-pp)))

#

    if (i > NN) {

        for (k in (NN+1):(2*NN))  {

          m1 = sum(testcups[1:(k-1)])

          m0 = (k-1) - m1

     if(m1 >= NN) testcups[k] = 0

     if(m0 >= NN) testcups[k] = 1

     if(m1 < NN && m0 < NN) testcups[k] =  sample(c(cups[k],(1-cups[k])), 1, prob=c(pp,(1-pp)))

          }                   # end of  for (k in (NN+1):(2*NN))

             }                # end of if (i > NN)

# 

                  }             #  end of 'else'

                             } # end of the for loops

#

   cupsC = ifelse(cups==0,"tea", "milk")

   testcupsC = ifelse(testcups == 0, "tea", "milk")

   sample_1 = c(pp, cupsC, testcupsC )

   sample_1           # this is the resulting sample of one set of 2*NN cups 

               }                # end of the    onesample() function

############################################################################

支持单位：