public
Created

mc_int

  • Download Gist
mc_int.r
R
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75
mc.int<-function(mod.in,int=NULL,n=10000,int.base=0,plot=T,plot.pts=T,
plot.poly=T,cols=list('black','green','blue'),round.val=2,val=F){
require(scales)
#variable names from model
if(is.null(mod.in$call$formula)){
x.names<-colnames(eval(mod.in$call$x))
y.names<-colnames(eval(mod.in$call$y))
x.range<-range(eval(mod.in$call$x))
}
else{
forms<-eval(mod.in$call$formula)
dat.names<-model.frame(forms,data=eval(mod.in$call$data))
y.names<-as.character(forms)[2]
x.names<-names(dat.names)[!names(dat.names) %in% y.names]
x.range<-range(eval(mod.in$call$data)[,x.names])
}
#integration range if null
if(is.null(int)) int<-x.range
#create data for model fit
x.fit<-seq(x.range[1],x.range[2],length=n)
pred.dat<-data.frame(x.fit)
names(pred.dat)<-x.names
y.fit<-predict(mod.in, newdata=pred.dat)
dat.fit<-data.frame(x.fit,y.fit)
#create random integration data
x.rand<-runif(n,int[1],int[2])
y.rand<-runif(n,int.base,max(y.fit))
dat.rand<-data.frame(x.rand,y.rand,y.under=rep(cols[[1]],n),stringsAsFactors=F)
names(dat.rand)[1]<-x.names
dat.rand$y.under[y.rand<predict(mod.in,newdata=dat.rand[,1,drop=F])]<-cols[[2]]
y.under<-dat.rand$y.under
#get integration
area.all<-diff(range(x.rand))*diff(range(y.rand))
area.und<-area.all*sum(y.under==cols[[2]])/length(y.under)
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(int.base,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.x<-data.frame(cord.x)
names(cord.x)<-x.names
cord.y<-c(int.base,predict(model,newdata=cord.x),int.base)
cord.x<-rbind(int[1],cord.x,int[2])
dat.poly<-data.frame(cord.x,cord.y)
names(dat.poly)[1]<-x.names
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
 
}

Please sign in to comment on this gist.

Something went wrong with that request. Please try again.