mc_int_ex

mc_int_ex.r

mc.int.ex<-function(x.range=c(-3,3),int=c(-1.96,1.96),n=10000,plot=T,
  plot.pts=T,plot.poly=T,cols=list('black','green','blue'),round.val=2,val=F){
	
	require(scales)
	
	#create data for model fit
	x.fit<-seq(x.range[1],x.range[2],length=n)
	y.fit<-dnorm(x.fit)
	dat.fit<-data.frame(x.fit,y.fit)
	
	#create random integration data
	x.rand<-runif(n,int[1],int[2])
	y.rand<-runif(n,0,max(y.fit))
	dat.rand<-data.frame(x.rand,y.rand,y.under=rep(cols[[1]],n),stringsAsFactors=F)
	dat.rand$y.under[y.rand<dnorm(dat.rand[,1])]<-cols[[2]]
	
	#get integration
	area.all<-diff(range(x.rand))*diff(range(y.rand))
	area.und<-area.all*sum(y.rand<dnorm(x.rand))/n
	out.txt<-paste('Integration from',int[1],'to',int[2],'=',round(area.und,round.val),
		sep=' ')
	
	if(plot){
		
		plot(x.rand,y.rand,xlim=x.range,ylim=c(0,max(y.fit)),type='n',main=out.txt)
		points(dat.fit$x.fit,dat.fit$y.fit,type='l')

		if(plot.poly){
			
			#create data for integration polygon
			cord.x<-seq(int[1],int[2],length=(n*diff(int)/diff(x.range)))
			cord.y<-c(0,dnorm(cord.x),0)
			cord.x<-c(int[1],cord.x,int[2])
			dat.poly<-data.frame(cord.x,cord.y)
			
			polygon(dat.poly[,1],dat.poly[,2],col=alpha(cols[[3]],0.5),border=NA)
			
			}
		
		if(plot.pts){
			
			points(dat.rand[,1],dat.rand[,2],col=alpha(dat.rand$y.under,0.5),pch=19)
			
			}

		}

	if(!val) return(cat(out.txt))
	area.und
	
	}