高尔顿钉板R语言实验

基于Matlab的Galton钉板问题黄自力高鹏黄安康摘要在概率论的发展过程中，最早出现的研究对象是一种计算概率的数学模型，称为古典概型。

一般的说，若随机试验满足下列两个条件：（1）它的样本空间只有有限多个样本点;（2）每个样本点出现的可能性相同，称这种实验为有限等可能实验或古典概型，galton钉板实验就是其中之一。

关键词galton顶板二项分布 poisson分布正文在概率论的发展过程中，最早出现的研究对象是一种计算概率的数学模型，称为古典概型。

一般的说，若随机试验满足下列两个条件：（1）它的样本空间只有有限多个样本点;（2）每个样本点出现的可能性相同，称这种实验为有限等可能实验或古典概型，galton钉板实验就是其中之一。

Galton钉板试验是英国生物统计学家Galton设计的。

在一板上钉有n排钉子，如图示，其中n=5。

右图中15个圆点表示15颗钉子，在钉子的下方有n+1个各子，分别编号为0，1，2，…，n。

从Galton钉板的上方扔进一个小球任其自由下落，在下落的过程中当小球碰到钉子时，从左边落下与从右边落下的机会相等。

碰到下一排钉子时又是如此。

最后落入底板中的某一个格子，图中用一条折线显示小球下落的一条轨迹。

向Galton钉板扔进一个小球，显然不能预测小球回落到哪一个格子，如果不断重复扔进过程，将会发生什么结果呢?关于Galton“高尔顿等人关于回归分析的先驱性的工作，以及时间序列分析方面的一些工作，…是数理统计学发展史中的重要事件.”──摘自《中国大百科全书》（数学卷）高尔顿是英国人类学家、生物统计学家.1822年2月6日生于伯明翰,1911年1月17日卒于萨里郡黑斯尔米尔.高尔顿是生物学家达尔文的表弟.他早年在剑桥学习数学，后到伦敦攻读医学.1860年当选为皇家学会会员，1909年被封为爵士.1845—1852年深入到非洲腹地探险、考察.高尔顿是生物统计学派的奠基人，他的表哥达尔文的巨著《物种起源》问世以后，触动他用统计方法研究智力遗传进化问题，第一次将概率统计原理等数学方法用于生物科学，明确提出“生物统计学”的名词.现在统计学上的“相关”和“回归”的概念也是高尔顿第一次使用的，他是怎样产生这些概念的呢？1870年，高尔顿在研究人类身长的遗传时，发现下列关系：高个子父母的子女，其身高有低于其父母身高的趋势，而矮个子父母的子女，其身高有高于其父母的趋势，即有“回归”到平均数去的趋势，这就是统计学上最初出现“回归”时的涵义.高尔顿揭示了统计方法在生物学研究中是有用的，引进了回归直线、相关系数的概念，创始了回归分析.开创了生物统计学研究的先河.他于1889年在《自然遗传》中，应用百分位数法和四分位偏差法代替离差度量.在现在的随机过程中有以他的姓氏命名的高尔顿─沃森过程（简称G─W 过程）.高尔顿发表了200篇论文和出版了十几部专著，涉及人体测量学，实验心理学等领域，其中数学始终起着重要作用.他在统计学方面也有贡献，高尔顿在1877年发表关于种子的研究结果，指出回归到平均值（regression toward the mean ）现象的存在，这个概念与现代统计学中的“回归”并不相同，但是却是回归一词的起源。

r语言实验报告R语言实验报告介绍•本文旨在对R语言实验报告进行相关介绍和指导。

准备工作•在开始编写R语言实验报告之前，需要进行一些准备工作：–安装R语言环境–确保安装必要的R包–理解实验要求和相关数据集实验报告结构•一个完整的R语言实验报告通常包含以下几个部分：1. 标题•实验报告的标题应简明扼要地描述实验内容。

2. 引言•引言部分应包含以下内容：–实验的背景和目的–实验所采用的数据集和方法的简要介绍3. 数据分析•数据分析部分是实验报告的重点，应包含以下内容：–数据的读取和预处理–数据的可视化–统计分析方法的应用–结果的解释和讨论4. 结论•结论部分应总结实验的结果，并对实验的目的和方法进行评价。

5. 参考文献•参考文献部分应列举实验报告中所引用的相关文献。

编写要点•在编写R语言实验报告时，需要遵守以下要点：1. 语法规范•使用清晰、准确的语法表达实验过程和结果。

2. 结果的解释•对于结果的解释，应该尽量采用简洁明了的语言，避免使用过于专业的术语或过于复杂的句子结构。

3. 图表的使用•图表是实验报告中常用的可视化工具，应合理使用图表来展示数据和结果，并配以简洁明了的图题和注解。

4. 逻辑性和连接性•实验报告应具有良好的逻辑性和连接性，各部分之间应有明确的联系和衔接，以确保整篇报告的连贯性。

结语•编写一份规范、完整的R语言实验报告需要系统的学习和实践，希望本文对您有所帮助。

参考文献•[参考文献1]•[参考文献2]继续编写一份更详细的R语言实验报告：R语言实验报告介绍•本文旨在对R语言实验报告进行相关介绍和指导。

准备工作•在开始编写R语言实验报告之前，需要进行一些准备工作：–安装R语言环境：确保在电脑上成功安装R语言的最新版本。

–确保安装必要的R包：根据实验需求，安装并加载所需的R包，例如ggplot2、dplyr等。

–理解实验要求和相关数据集：认真阅读实验要求，理解实验的目的和需求，并熟悉所使用的数据集。

R语言实验指导书（二）2016年10月27日实验三创建和使用R语言数据集一、实验目的：1.了解R语言中的数据结构。

2.熟练掌握他们的创建方法，和函数中一些参数的使用。

3.对创建的数据结构进行，排序、查找、删除等简单的操作。

二、实验容：1.向量的创建及因子的创建和查看有一份来自澳大利亚所有州和行政区的20个税务会计师的信息样本 1 以及他们各自所在地的州名。

州名为：tas, sa, qld, nsw, nsw, nt, wa, wa, qld, vic, nsw, vic, qld, qld, sa, tas, sa, nt, wa, vic。

1)将这些州名以字符串的形式保存在state当中。

2)创建一个为这个向量创建一个因子statef。

3)使用levels函数查看因子的水平。

2.矩阵与数组。

i.创建一个4*5的数组如图，创建一个索引矩阵如图，用这个索引矩阵访问数组，观察结果。

3.将之前的state，数组，矩阵合在一起创建一个长度为3的列表。

4.创建一个数据框如图。

5.将这个数据框按照mpg列进行排序。

6.访问数据框中drat列值为3.90的数据。

三、实验要求要求学生熟练掌握向量、矩阵、数据框、列表、因子的创建和使用。

实验四数据的导入导出一、实验目的1.熟练掌握从一些包中读取数据。

2.熟练掌握csv文件的导入。

3.创建一个数据框，并导出为csv格式。

二、实验容1.创建一个csv文件（容自定），并用readtable函数导入该文件。

2.查看R语言自带的数据集airquality（纽约1973年5-9月每日空气质量）。

3.列出airquality的前十列，并将这前十列保存到air中。

4.查看airquality中列的对象类型。

5.查看airquality数据集中各成分的名称6.将air这个数据框导出为csv格式文件。

（write.table (x, file ="", sep ="",s =TRUE, s =TRUE, quote =TRUE)）三、实验要求要求学生掌握从包中读取数据，导入csv文件的数据，并学会将文件导出。

Galton钉板实验一、实验内容某车间有200台车床互相独立的工作，由于经常需要检修、测量、调换刀具等种种原因需要停车，这使每台车床的开工率只有60％。

而每台车床在开动时需耗电1kW，显然向该车间供电200kW可以保证有足够电力供这些车床使用，但是在电力比较紧张的情况下，给这个车间供给电力太多将造成浪费，太少又影响生产。

如何解决这一矛盾？一种解决方案是保证有基本足够的电力供应该车间，比如要求在8小时的生产过程中允许有半分钟的电力不足，半分钟约占8小时的0.1％，用概率论的语言就是：应供应多少电力才能以99.9％的概率保证不会因为电力不足而影响生产？问题：（1）计算分布函数在某些点的取值F(m)，m＝0，1，2， (200)并将它绘于图上，辅助某些必要的计算，求出问题中所需要的供电功率数。

（2）将8小时按半分钟分成若干时间段，共有8*60*2＝960个时间段。

用二项分布模拟8小时车床运行的情况。

观察已算得的供电功率数是否能基本满足车间正常工作，写出你的结论。

二、实验过程问题（1）编写程序如下：function bin() %200台车床正常工作的台数满足二项分布p=0.6; %正常工作概率n=200; %200次事件x=[0:5:n];f=binocdf(x,n,p);bar(x,f);axis([-1 201 0 1]); %坐标分配end运行结果：将上述程序的取样间隔改为一时，即x=[0:5:n]; 改为x=[0:1:n];结果如下：通过观察上面两幅结果，得出大约在m=140KW时电力才能以99.9％的概率保证不会因为电力不足而影响生产。

问题（2）模拟车床运行情况的函数代码为：function bin1n=200;p=0.6;m=960;rand('seed',3);R=binornd(n,p,1,m); %模拟服从二项分布的随机数,生成1*960的矩阵for i=1:n+1 %开始计数k=[];k=find(R==(i-1)); %找出R中等于(i-1)元素下标，并存于向量k中h(i)=length(k)/m; %计算落在编号i-1的格子的小球频率endx=[0:1:n];Bar(x,h);axis([-1 201 0 1]) ; %画频率图end运行后生成的分布图为：输入以下代码，计算服从n=200,p=0.6的二项分布的随机变量的分布列的理论值：function bin2n=200;p=0.6;x=[0:1:n];f=binopdf(x,n,p);bar(x,f);axis([-1 201 0 1]);end得到理论分布图为：通过对两图的对比可以看出，当进行大量次重复投球后，小球的堆积形状和理论上的分布情况（随机变量X的分布列）非常接近。

实验三数组的运算、求解方程（组）和函数极值、数值积分【实验类型】验证性【实验学时】2 学时【实验目的】1、掌握向量的四则运算和内积运算、矩阵的行列式和逆等相关运算；2、掌握线性和非线性方程（组）的求解方法，函数极值的求解方法；3、了解 R 中数值积分的求解方法。

【实验内容】1、向量与矩阵的常见运算；2、求解线性和非线性方程（组）；3、求函数的极值，计算函数的积分。

【实验方法或步骤】第一部分、课件例题：1.向量的运算x<-c(-1,0,2)y<-c(3,8,2)v<-2*x+y+1vx*yx/yy^xexp(x)sqrt(y)x1<-c(100,200); x2<-1:6; x1+x22.x<-1:5y<-2*1:5x%*%ycrossprod(x,y)x%o%ytcrossprod(x,y)outer(x,y)3.矩阵的运算A<-matrix(1:9,nrow=3,byrow=T);AA+1 #A的每个元素都加上1B<-matrix(1:9,nrow=3); BC<-matrix(c(1,2,2,3,3,4,4,6,8),nrow=3); C D<-2*C+A/B; D #对应元素进行四则运算x<-1:9A+x #矩阵按列与向量相加E<-A%*%B; E #矩阵的乘法y<-1:3A%*%y #矩阵与向量相乘crossprod(A,B) #A的转置乘以Btcrossprod(A,B) #A乘以B的转置4.矩阵的运算A<-matrix(c(1:8,0),nrow=3);At(A) #转置det(A) #求矩阵行列式的值diag(A) #提取对角线上的元素A[lower.tri(A)==T]<-0;A #构造A对应的上三角矩阵qr.A<-qr(A);qr.A #将矩阵A分解成正交阵Q与上三角阵R的乘积，该结果为一列表Q<-qr.Q(qr.A);Q;R<-qr.R(qr.A);R #显示分解后对应的正交阵Q与上三角阵Rdet(Q);det(R);Q%*%R #A=Q*Rqr.X(qr.A) #显示分解前的矩阵5.解线性方程组A<-matrix(c(1:8,0),nrow=3,byrow=TRUE)b<-c(1,1,1)x<-solve(A,b); x #解线性方程组Ax=bB<-solve(A); B #求矩阵A的逆矩阵BA%*%B #结果为单位阵6.非线性方程求根f<-function(x) x^3-x-1 #建立函数uniroot(f,c(1,2)) #输出列表中f.root为近似解处的函数值，iter为迭代次数，estim.prec为精度的估计值uniroot(f,lower=1,upper=2) #与上述结果相同polyroot(c(-1,-1,0,1)) #专门用来求多项式的根，其中c(-1,-1,0,1)表示对应多项式从零次幂项到高次幂项的系数7.求解非线性方程组（1）自编函数： (Newtons.R)Newtons<-function (funs, x, ep=1e-5, it_max=100){index<-0; k<-1while (k<=it_max){ #it_max 表示最大迭代次数x1 <- x; obj <- funs(x);x <- x - solve(obj$J, obj$f); #Newton 法的迭代公式norm <- sqrt((x-x1) %*% (x-x1))if (norm<ep){ index<-1; break #index=1 表示求解成功}; k<-k+1 }obj <- funs(x);list(root=x, it=k, index=index, FunVal= obj$f)} # 输出列表（2）调用求解非线性方程组的自编函数funs<-function(x){ f<-c(x[1]^2+x[2]^2-5, (x[1]+1)*x[2]-(3*x[1]+1)) # 定义函数组J<-matrix(c(2*x[1], 2*x[2], x[2]-3, x[1]+1), nrow=2,byrow=T) # 函数组的 Jacobi 矩阵list(f=f, J=J)} # 返回值为列表 : 函数值 f 和 Jacobi 矩阵 Jsource("F:/wenjian_daima/Newtons.R") # 调用求解非线性方程组的自编函数Newtons(funs, x=c(0,1))8.一元函数极值f<-function(x) x^3-2*x-5 # 定义函数optimize(f,lower=0,upper=2) # 返回值 : 极小值点和目标函数f<-function(x,a) (x-a)^2 # 定义含有参数的函数optimize(f,interval=c(0,1),a=1/3) # 在函数中输入附加参数9.多元函数极值（1）obj <-function (x){ # 定义函数F<-c(10*(x[2]-x[1]^2),1-x[1]) # 视为向量sum (F^2) } # 向量对应分量平方后求和nlm(obj,c(-1.2,1))（2）fn<-function(x){ # 定义目标函数F<-c(10*(x[2]-x[1]^2), 1-x[1])t(F)%*%F } # 向量的内积gr <- function(x){ # 定义梯度函数F<-c(10*(x[2]-x[1]^2), 1-x[1])J<-matrix(c(-20*x[1],10,-1,0),2,2,byrow=T) #Jacobi 矩阵2*t(J)%*%F } # 梯度optim(c(-1.2,1), fn, gr, method="BFGS")最优点 (par) 、最优函数值 (value)10.梯形求积分公式（1）求积分程序： (trape.R)trape<-function(fun, a, b, tol=1e-6){ # 精度为 10 -6N <- 1; h <- b-a ; T <- h/2 * (fun(a) + fun(b)) # 梯形面积 repeat{h <- h/2; x<-a+(2*1:N-1)*h; I <-T/2 + h*sum(fun(x)) if(abs(I-T) < tol) break; N <- 2 * N; T = I }; I}（2）source("F:/wenjian_daima/trape.R") # 调用函数f<-function(x) exp(-x^2)trape(f,-1,1)（3）常用求积分函数f<-function(x)exp(-x^2) # 定义函数integrate(f,0,1)integrate(f,0,10)integrate(f,0,100)integrate(f,0,10000) # 当积分上限很大时，结果出现问题integrate(f,0,Inf) # 积分上限为无穷大ft<-function(t) exp(-(t/(1-t))^2)/(1-t)^2 # 对上述积分的被积函数 e 2 作变量代换 t=x/(1+x) 后的函数integrate(ft,0,1) # 与上述计算结果相同，且精度较高第二部分、教材例题：1.随机抽样（1）等可能的不放回的随机抽样:> sample(x, n) 其中x为要抽取的向量, n为样本容量（2）等可能的有放回的随机抽样:> sample(x, n, replace=TRUE)其中选项replace=TRUE表示有放回的, 此选项省略或replace=FALSE表示抽样是不放回的sample(c("H", "T"), 10, replace=T)sample(1:6, 10, replace=T)（3）不等可能的随机抽样:> sample(x, n, replace=TRUE, prob=y)其中选项prob=y用于指定x中元素出现的概率, 向量y与x等长度sample(c("成功", "失败"), 10, replace=T, prob=c(0.9,0.1))sample(c(1,0), 10, replace=T, prob=c(0.9,0.1))2.排列组合与概率的计算1/prod(52:49)1/choose(52,4)3.概率分布qnorm(0.025) #显著性水平为5%的正态分布的双侧临界值qnorm(0.975)1 - pchisq(3.84, 1) #计算假设检验的p值2*pt(-2.43, df = 13) #容量为14的双边t检验的p值4.limite.central( )的定义limite.central <- function (r=runif, distpar=c(0,1), m=.5,s=1/sqrt(12),n=c(1,3,10,30), N=1000) {for (i in n) {if (length(distpar)==2){x <- matrix(r(i*N, distpar[1],distpar[2]),nc=i)}else {x <- matrix(r(i*N, distpar), nc=i)}x <- (apply(x, 1, sum) - i*m )/(sqrt(i)*s)hist(x,col="light blue",probability=T,main=paste("n=",i), ylim=c(0,max(.4, density(x)$y)))lines(density(x), col="red", lwd=3)curve(dnorm(x), col="blue", lwd=3, lty=3, add=T)if( N>100 ) {rug(sample(x,100))}else {rug(x)}}}5.直方图x=runif(100,min=0,max=1)hist(x)6.二项分布B（10，0.1）op <- par(mfrow=c(2,2))limite.central(rbinom,distpar=c(10,0.1),m=1,s=0.9)par(op)7.泊松分布: pios（1）op <- par(mfrow=c(2,2))limite.central(rpois, distpar=1, m=1, s=1, n=c(3, 10, 30 ,50)) par(op)8.均匀分布：unif(0,1)op <- par(mfrow=c(2,2))limite.central( )par(op)9.指数分布：exp(1)op <- par(mfrow=c(2,2))limite.central(rexp, distpar=1, m=1, s=1)par(op)10.混合正态分布的渐近正态性mixn <- function (n, a=-1, b=1){rnorm(n, sample(c(a,b),n,replace=T))}limite.central(r=mixn, distpar=c(-3,3),m=0, s=sqrt(10), n=c(1,2,3,10)) par(op)11.混合正态分布的渐近正态性op <- par(mfrow=c(2,2))mixn <- function (n, a=-1, b=1){rnorm(n, sample(c(a,b),n,replace=T))}limite.central(r=mixn, distpar=c(-3,3),m=0,s=sqrt(10),n=c(1,2,3,10)) par(op)第三部分、课后习题：3.1a=sample(1:100,5)asum(a)3.2（1）抽到10、J、Q、K、A的事件记为A，概率为P(A)=1(5220)其中在R中计算得：> 1/choose(52,20)[1] 7.936846e-15（2）抽到的是同花顺P(B)=(41)(91) (525)在R中计算得：> (choose(4,1)*choose(9,1))/choose(52,5) [1] 1.385e-053.3#(1)x<-rnorm(1000,mean=100,sd=100)hist(x)#(2)y<-sample(x,500)hist(y)#(3)mean(x)mean(y)var(x)var(y)3.4x<-rnorm(1000,mean=0,sd=1) y=cumsum(x)plot(y,type = "l")plot(y,type = "p")3.5x<-rnorm(100,mean=0,sd=1) qnorm(.025)qnorm(.975)t.test(x)由R结果知：理论值为[-1.96,1.96]，实际值为：[-0.07929,0.33001]3.6op <- par(mfrow=c(2,2))limite.central(rbeta, distpar=c(0.5 ,0.5),n=c(30,200,500,1000))par(op)3.7N=seq(-4,4,length=1000)f<-function(x){dnorm(x)/sum(dnorm(x))}n=f(N)result=sample(n,replace=T,size = 1000)standdata=rnorm(1000)op<-par(mfrow=c(1,2)) #1行2列数组按列（mfcol）或行（mfrow）各自绘图hist(result,probability = T)lines(density(result),col="red",lwd=3)hist(standdata,probability = T)lines(density(standdata),col="red",lwd=3) par(op)。

基于Matlab的Galton钉板问题基于Matlab的Galton钉板问题黄自力高鹏黄安康摘要在概率论的发展过程中，最早出现的研究对象是一种计算概率的数学模型，称为古典概型。

一般的说，若随机试验满足下列两个条件：（1）它的样本空间只有有限多个样本点;（2）每个样本点出现的可能性相同，称这种实验为有限等可能实验或古典概型，galton钉板实验就是其中之一。

关键词galton顶板二项分布 poisson分布正文在概率论的发展过程中，最早出现的研究对象是一种计算概率的数学模型，称为古典概型。

一般的说，若随机试验满足下列两个条件：（1）它的样本空间只有有限多个样本点;（2）每个样本点出现的可能性相同，称这种实验为有限等可能实验或古典概型，galton钉板实验就是其中之一。

Galton钉板试验是英国生物统计学家Galton设计的。

在一板上钉有n排钉子，如图示，其中n=5。

右图中15个圆点表示15颗钉子，在钉子的下方有n+1个各子，分别编号为0，1，2，…，n。

从Galton钉板的上方扔进一个小球任其自由下落，在下落的过程中当小球碰到钉子时，从左边落下与从右边落下的机会相等。

碰到下一排钉子时又是如此。

最后落入底板中的某一个格子，图中用一条折线显示小球下落的一条轨迹。

向Galton钉板扔进一个小球，显然不能预测小球回落到哪一个格子，如果不断重复扔进过程，将会发生什么结果呢?关于Galton“高尔顿等人关于回归分析的先驱性的工作，以及时间序列分析方面的一些工作，…是数理统计学发展史中的重要事件.”──摘自《中国大百科全书》（数学卷）高尔顿是英国人类学家、生物统计学家.1822年2月6日生于伯明翰,1911年1月17日卒于萨里郡黑斯尔米尔.高尔顿是生物学家达尔文的表弟.他早年在剑桥学习数学，后到伦敦攻读医学.1860年当选为皇家学会会员，1909年被封为爵士.1845—1852年深入到非洲腹地探险、考察.高尔顿是生物统计学派的奠基人，他的表哥达尔文的巨著《物种起源》问世以后，触动他用统计方法研究智力遗传进化问题，第一次将概率统计原理等数学方法用于生物科学，明确提出“生物统计学”的名词.现在统计学上的“相关”和“回归”的概念也是高尔顿第一次使用的，他是怎样产生这些概念的呢？1870年，高尔顿在研究人类身长的遗传时，发现下列关系：高个子父母的子女，其身高有低于其父母身高的趋势，而矮个子父母的子女，其身高有高于其父母的趋势，即有“回归”到平均数去的趋势，这就是统计学上最初出现“回归”时的涵义.高尔顿揭示了统计方法在生物学研究中是有用的，引进了回归直线、相关系数的概念，创始了回归分析.开创了生物统计学研究的先河.他于1889年在《自然遗传》中，应用百分位数法和四分位偏差法代替离差度量.在现在的随机过程中有以他的姓氏命名的高尔顿─沃森过程（简称G─W 过程）.高尔顿发表了200篇论文和出版了十几部专著，涉及人体测量学，实验心理学等领域，其中数学始终起着重要作用.他在统计学方面也有贡献，高尔顿在1877年发表关于种子的研究结果，指出回归到平均值（regression toward the mean ）现象的存在，这个概念与现代统计学中的“回归”并不相同，但是却是回归一词的起源。

R语⾔中进⾏期权定价的Heston模型原⽂链接：/?p=12111在本⽂中，我将向您展⽰如何模拟股票价格的Heston随机波动率模型。

Heston模型是针对具有随机波动性的期权，并于1993年申请了债券的货币期权。

对于固定的⽆风险利率，其描述为：通过使⽤这种模型，可以得出欧洲看涨期权的价格。

这是函数的描述。

callHestoncf(S, X, tau, r, v0, vT, rho, k, sigma){# S = Spot, X = Strike, tau = time to maturity# r = risk-free rate, q = dividend yield# v0 = initial variance, vT = long run variance (theta)# rho = correlation, k = speed of mean reversion (kappa)# sigma = volatility of volatility}现在，进⾏蒙特卡洛定价。

我们将为3个欧洲看涨期权定价，具有3种不同的执⾏价格。

我们在15年中使⽤100000个模拟，每个⽉进⾏⼀次。

以下是对仿真有⽤的参数：#Initial stock priceS0 <- 100# Number of simulations (feel free to reduce this)n <- 100000# Sampling frequencyfreq <- "monthly"# volatility mean-reversion speedkappa <- 0.003# volatility of volatilityvolvol <- 0.009# Correlation between stoch. vol and spot pricesrho <- -0.5# Initial varianceV0 <- 0.04# long-term variancetheta <- 0.04#Initial short rater0 <- 0.015 # Options maturitieshorizon <- 15# Options' exercise pricesstrikes <- c(140, 100, 60)为了使⽤模拟Heston模型，我们⾸先需要定义如何进⾏模拟。

R语言因子实验设计和解释案例分析报告示例1：两组比较示例2：多个组实例3：两个条件，两个基因型，一个交互项o野生型治疗效果（主效应）。

o突变体治疗的效果o没有治疗的突变型和野生型之间有什么区别？o通过治疗，突变型和野生型有什么区别？o基因型的不同反应（相互作用项）实例4：两个条件，三个基因型o基因型I的条件效应（主效应）o基因型III的条件效应。

o基因型II的条件效应。

o在条件A下III与II的影响o基因型III与基因型I的条件效应的相互作用项○基因型III与基因型II的条件效应的相互作用项。

为了允许iDEP中的复杂模型（/idep/），我尝试了解如何构建事实模型，并从DESeq2中提取期望的结果。

以下是基于DESeq2中resutls（）函数的帮助文档，以及Mike Love对用户提问的回答。

我想要做的一个重点是，当研究设计涉及多个因素时（参见上面关于基因型+治疗实例的图），结果的解释是棘手的。

与R中的回归分析类似，分类因素的参考水平构成了我们的分歧的基础。

然而，默认情况下，它们是按字母顺序确定的。

选择每个因素的参考水平是至关重要的。

否则你的系数可能会有所不同，这取决于你如何进入DESeq2的实验设计。

这可以通过R中的relevel（）函数完成。

参考级别是构成有意义比较基础的因素的基线级别。

在野生型与突变型实验中，“野生型”是参考水平。

在治疗与未治疗，参考水平显然是未经处理的。

例3中的更多细节。

例1：两组比较首先制作一些示例数据。

这是一个非常简单的实验设计，有两个条件。

这显示了可用的结果。

请注意，默认情况下，R会根据字母顺序为因素选择一个参考级别。

这里A是参考水平。

折叠变化定义为B与A比较。

要更改参考级别，请尝试使用“同一个”（）函数。

baseMean log2FoldChange pvalue padj gene9056 360.168909 -2.045379 0.0000000 0.0001366 gene3087 43.897516 -2.203303 0.0000173 0.0858143 gene3763 72.409877 -1.834787 0.0000434 0.1434712 gene2054 322.494963 1.537408 0.0000681 0.1689463 gene4617 6.227415 6.125238 0.0002019 0.4008408 如果我们想用B作为控制，并用B作为基线定义倍数变化。

高尔顿（G a l t o n ）钉板实验一、问题描述Galton 钉板试验是英国生物统计学家Galton 设计的。

在一板上钉有n 排钉子，如图示，其中n=5。

右图中15个圆点表示15颗钉子，在钉子的下方有n+1个各子，分别编号为0，1，2，…，n 。

从Galton 钉板的上方扔进一个小球任其自由下落，在下落的过程中当小球碰到钉子时，从左边落下与从右边落下的机会相等。

碰到下一排钉子时又是如此。

最后落入底板中的某一个格子，图中用一条折线显示小球下落的一条轨迹。

二、高尔顿钉板试验中的相关问题1、小球落入各个格子中的概率与频数做一个小球的高尔顿钉板试验，其落入第i 个格子的概率正好满足二项分布。

设高尔顿钉板有n 行钉，第n 行铁钉共有n 个，有（n+1）个空。

把这（n+1）个空由左到右依次编号为i=0，1，2，…，n 共（n+1）个空。

观察i=0这个空，小球从这个空落下的条件是：小球从第一次与铁钉碰撞后必须连续向左落下，即连续n 次选择向左落下，所以落入第i=0个空的概率为P （i=0）=C 0n （21）n (21)0。

观察i=1这个空，小球从这个空落下的条件是：小球从第一次与铁钉碰撞后连续n 次碰撞落下过程中，有且只有一次选择向右落下，其余都只能是向左落下，所以落入第i=1个空的概率为P （i=1）=C 1n （21）n-1(21)1。

小球从第一次与铁钉碰撞后连续n 次碰撞落下过程中，有i 次选择向右落下，其余都选择向左落下，所以落入第i 个空的概率为P （i ）= C i n （12）n-i (21)i （i=0，1，2，…，n ）。

故，当一个一个从顶部放入k 个小球，低槽中各格的理论频数为：h(i)=k ×P(i),(i=0,1,2,…，n).2、程序运行2.1 基本功能①输入小球数k、概率p;②计算高尔顿钉板n=4时，放入k个小球后，落入底槽各格中的实验小球数；③计算高尔顿钉板n=4时，放入k个小球后，落入底槽各格中的理论小球数；④动画演示每个小球下落路径及底槽各格小球数频率增长情况；④画出落入底槽各格中的实验小球数频率的柱状图；⑤画出落入底槽各格中的实验小球数、落入底槽各格中的理论实验小球数的频率曲线图；⑥关闭。