R语言学习记录-4

痛定思痛。 2023-06-13 06:10 17阅读 0赞

#估计二重积分
    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