R语言学习记录-4

#估计二重积分
f = function(x,y)
{
	temp = exp( - (x + y)^2)
	return (temp)
}
mc2 = function(n,f)
{
	x = runif(n)
	y = runif(n)
	sum(f(x,y))/n
}
#取1000个均匀随机数进行验证,可得
mc2(1000,f)

#计算二重积分
(exp(4) - exp(3)) * (exp(2) - exp(1))
#下面用Monte-Carlo方法计算,定义函数
f = function(x,y)
{	exp(x+y)}
mc3 = function(n,a1,b1,a2,b2,f)
{	x = runif(n)
	y = runif(n)
	x1 = (b1 - a1) * x + a1
	y1 = (b2 - a2) * y + a2
	sum((b2-a2)*(b1-a1)*f(x1,y1))/n
}
#可见,随着随机点达到10000时,所得到的结果与真值非常接近.
mc3(10000,1,2,3,4,f)

#程序控制结构
if结构
f = function(x)
{	if (x > 0)
	return (sqrt(x))
	else 
	return (0)
}
#或者
f=function(x){ifelse(x>0,sqrt(x),0)}
f(5)
#需要注意的是,上述程序中如果将f(x)中变量定义为向量时,
程序会运行出错,因为if判断条件是标量的真值或假值,
而不能判断向量.另外,有多个if语句时,else与最近的一个if匹配.
可使用if…else if …else if …else的多重判断结构表示多分支,
多分支也可以使用swith()函数.例如
f = function(x)
{	if (x < -5) 0
	else if (x < 5) 1
	else if (x < 8) 2
	else 3}

f(-8)

#例4.3.1
#Pareto 分布是意大利经济学家Pareto 将其作为一种收入分布而介
绍的,它不仅在经济收入模型中得到应用,在其它领域中也得到广泛的应用.
#定义密度函数
dpareto = function(x,a,theta)
{	if (x > a) z = t/a *(a / x)^ (thera + 1) 
	else z = 0
	z}
#分布函数为
ppareto = function(x,a,thera)
{	if (x > a) z = 1 - (a/x)^thera
	else z = 0
	z
}
ppareto(1,2,3)
ppareto(5,2,3)
#所以定义的分位数函数为:
qpareto = function(p,a,thera)
{	if (p < 1 & p >= 0) z = a * exp( - log(1 - p) / thera)
	else z = 'FALSE'
	return (z)
}
qpareto(11,2,3)
qpareto(0.7,2,3)

###
for循环
for循环是对一个向量和列表逐次处理,格式为:
for(name in value){表达式1
表达式2
……}
###
t = numeric(15)
for (i in 1 : 15){t[i] = (i - 1) * i}
t

#以例4.1.4中讨论的艾拉姆咖分布随机数的产生方法为例,可以运行以下程序
library(rootSolve)
Fa = function(x,t)
{	1 - (1 + 2 * t * x) * exp(-2 * t * x)}

u = runif(10)
t = numeric(10)
for(i in 1 : 10)
{
	m = function(y)
	{	
		Fa(y,0.25) - u[i]
	}
	t[i] = uniroot.all(m,c(0,12))
}
t

Fa = function(x,t)
{	1 - (1 + 2 * t * x) * exp(-2 * t * x)}
rfa = function(n)
	{	
		u = runif(n)
		t = numeric(n)
		for(i in 1 : n){
		m = function(y){Fa(y,0.25) - u[i]}
		t[i] = uniroot.all(m,c(0,15))
	}
t
}
rfa(100)


#例4.3.3
p = numeric(365)
for(i in 1 : 365)
{p[i] = 1 - prod((365 : (365 - i + 1)) /365)}
p[10]

###
另外,for循环也可用于构造矩阵,只是需要for循环嵌套.
在矩阵给定的初值用array(0,c(n,n))形式定义,
array()函数的完全使用为array(x, dim=length(x), 
dimnames=NULL),其中x是第一自变量,应该是一个向量,
表示数组的元素值组成的向量.dim参数可省,
省略时作为一维数组(但不同于向量).
在for循环中类似于一维情形下numeric()函数的作用.

n = 8
x = array(0,c(n,n))
for(i in 1:n){
	for(j in 1:n){
		x[i,j] = i * j}
		}
x

#while循环
while循环是在开始处判断条件的当型循环,
它与for循环的不同之处在于:若知道终止条件(循环次数)就用for循环,
若无法知道运行次数,则用while循环或repeat循环.

#用while循环编写一个计算1000以内的Fibonacci数列
f = 1
f[2] = 1
i = 1
while(f[i] + f[i + 1] < 10000000){
f[i + 2] = f[i] + f[i + 1]
i = i + 1}
f