R语言学习记录-3

超、凢脫俗 2023-05-30 04:58 23阅读 0赞

#最大值和最小值分布
    f1 = function(x)#串联系统
    { 
    pg = pgamma(x,2,0.01)
    dg = dgamma(x,2,0.01)
    5 * dg * (1 - pg)^4
    }
    f2 = function(x)#并联系统
    { 
    pg = pgamma(x,2,0.01)
    dg = dgamma(x,2,0.01)
    5 * dg * (pg)^4
    }
    
    G = function(y,f)
    { 
    J = integrate(f,lower = 0,upper = y)$value
    return (1 - J)
    }
    
    G(100,f2)
    G(100,f1)
    G(200,f2)
    G(200,f1)
    G(300,f2)
    G(300,f1)
    
    
    #混合分布
    #以正态分布为例,构造一个两样本混合正态分布
    mixN2 = function(x,p)
    { 
    p * dnorm(x) + (1 - p) * dnorm(x,10,1)
    }
    z = seq(-4,20,0.1)
    op = par(mfrow = c(1,3))
    plot(z,mixN2(z,0.5),type="l",ylim=c(0,0.4))
    plot(z,mixN2(z,0.2),type="l",ylim=c(0,0.4),col=2)
    plot(z,mixN2(z,0.8),type="l",ylim=c(0,0.4),col=3)
    par(op)
    
    
    mixN3 = function(x,p,mu,sigma)
    { 
    mixf = p[1] * dnorm(x,mu[1],sigma[1]) + p[2] * dnorm(x,mu[2],sigma[2])+ p[3] * dnorm(x,mu[3],sigma[3])
    }
    
    
    z = seq(-4,20,0.1)
    op = par(mfrow = c(1,3))
    plot(z,mixN3(z,rep(1/3,3),c(0,3,6),c(1,1,1)),type="l",ylim=c(0,0.4))
    plot(z,mixN3(z,0.2),type="l",ylim=c(0,0.4),col=2)
    plot(z,mixN3(z,0.8),type="l",ylim=c(0,0.4),col=3)
    par(op)
    
    #R中虽然自带许多常见的分布,而有的分布却不存在现成的程序包,这就需要用函数的方
    #式来定义它.以艾拉姆咖分布为例.
    fa=function(x,t){ 4*t^2*x*exp(-2*t*x)}
    Fa=function(x,t){ 1-(1+2*t*x)*exp(-2*t*x)}
    z=seq(0,15,0.1)
    plot(z,fa(z,1),type="l",main="艾拉姆咖分布密度函数")
    lines(z,fa(z,1/2),col=2)
    lines(z,fa(z,1/3),col=3)
    lines(z,fa(z,1/4),col=4)
    
    #随机数的产生
    rt(10,6)
    #在R中,uniroot.all(F,c(a,b))函数可以解出(a,b)上关于函数F=0的所有根,需要加载root-
    #Solve包,运行下列程序:
    library(rootSolve)
    u = runif(1)
    Fa = function(x,t){ 1-(1+2*t*x)*exp(-2*t*x)}
    m = function(y)
    { 	J = Fa(y,0.25) - u
    	return (J)
    }
    uniroot.all(m,c(0,15))	
    #这样,我们就可以得到一系列参数为0.25的艾拉姆咖分布随机数.
    这种取随机数的方法虽然原则上可行,但每次产生一个随机数都需要运行后三条语句才得到一个结果,
    效率较低.后面章节中将介绍for循环语句可以有效的提高工作效率.
    4.2 随机数的产生及应.
    #拓展问题：如何同时生成多个随机数
    #随机数的应用(Monte- Carlo方法)
    #均匀分布的随机数一个重要的应用就是Monte- Carlo方法
    
    #一元函数的定积分值估计
    mc = function(n,g)
    	{ 
    	x = runif(n)
    	J = mean(g(x))
    	return (J)
    	}
    
    g1 = function(x)
    { 	exp(-x^2)
    }
    
    mc(10000000,g1)
    #准确值
    integrate(g1,0,1)
    
    z=seq(0,1,length=100)
    plot(z,runif(1000),ylim=c(0,1))
    lines(z,g1(z),lwd=2,col=2)