d <- read.table('inp_reg1.txt',skip=1) vec_x <- d[,1] vec_y <- log10(d[,2]) nn <- length(d) xm<-mean(vec_x) sx<-sd(vec_x) ym<-mean(vec_y) sy<-sd(vec_y) n<-360 vec_deg<-seq(0,2*pi,length=n) rho <- cor(vec_x,vec_y) pp=0.95 vec_rr<-sqrt(-2*(1-rho^2)*log(1-pp)/(1-2*rho*sin(vec_deg)*cos(vec_deg))) pltx<-xm+sx*vec_rr*cos(vec_deg) plty<-ym+sy*vec_rr*sin(vec_deg) plot(vec_x,vec_y,xlim=c(45,50),ylim=c(1,4),type="p",ann=F) par(new=TRUE) plot(pltx,plty,xlim=c(45,50),ylim=c(1,4),type="l",lty=2,xlab="Elevation (m)",ylab="Observed value") library(ellipse) polygon(ellipse(rho,scale=c(sx,sy),centre=c(xm,ym),level=0.95,npoints=360),xlim=c(45,50),ylim=c(1,4))