R语言学习记录-3

#最大值和最小值分布
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)