breathing forecast

I found an opportune time to collect some data while my roommate slept (creepy, but whatever, it's for science). I timed his inhaling and exhaling (in seconds) on my phone and ran some R code.

With this data I look at time-series plots, density plots overlaid with a theoretical distribution for the data, and forecasting using an SARIMA model (library: astsa). For no reason other than fun, I forecast his next 10 breaths (five inhales and five exhales, which happened long ago).

The purpose of this post is to give examples of some of the things that can be done with the plot and sarima functions.

data:
The odd entries of breathe.dat were my roommate's inhales and the even entries were his exhales. The observations will undoubtedly include measurement error (I just used my phone after all and my reaction time isn't perfect) and also some error involved with his breathing. In this example I won't bother trying to model my roommate's breathing all that well. Why? Because who cares.

breathe.dat=c(
  1.2, 2.1, 1.1, 3.0, 1.2, 2.6, 1.8, 2.5,
1.3, 2.4, 1.8, 2.4, 1.5, 2.6, 1.4, 2.5,
1.2, 2.9, 1.5, 2.2, 1.4, 2.2, 1.3, 2.1,
1.4, 2.5, 1.6, 2.1, 1.3, 2.2, 1.3, 2.7,
1.4, 2.3, 1.5, 2.3, 1.5, 2.0, 1.3, 2.4,
1.3, 2.6, 1.5, 2.6, 1.1, 1.9, 1.0, 1.9,
1.1, 2.7, 1.1, 2.1, 1.4, 2.3, 1.3, 2.5,
1.3, 2.6, 1.4, 2.8, 1.3, 2.4, 1.4, 2.1,
1.2, 1.9, 1.4, 3.0, 1.2, 2.3, 1.4, 2.4,
1.1, 2.2, 1.4, 2.3, 1.3, 2.3, 1.4, 2.6,
1.5, 2.2, 1.3, 2.3, 1.5, 2.4, 1.3, 2.3,
1.3, 2.4, 1.4, 2.7, 1.2, 2.3, 1.3, 2.6,
1.6, 2.2, 1.5, 2.5, 1.4, 2.0)

n=length(breathe.dat) #n=102, 51 inhale, 51 exhale
ins=breathe.dat[seq(1,n,2)]  #gets subset entries 1,3,5,...,101
outs=breathe.dat[seq(2,n,2)] #gets subset entries 2,4,6,...,102

Here I split up the data in the inhales (ins) and exhales (outs).

plots:
# TIME-SERIES
path="c:/users/mickey/desktop/downloads/school/r play/plots/"
name=c("in","out","both")
dats=list(ins,outs,breathe.dat)
for (i in 1:3){
  png(paste0(path,"breathe ",name[i],".png"))
  plot(dats[[i]],type='b',ylab="Time (s)",
  main=paste0("Breathe ",name[i]),axes=FALSE)
    axis(1);axis(2);grid()
  dev.off()    #closes active device (the plot)
  }

# DENSITY

png(paste0(path,"hist in.png"))

hist(ins,main="Breathe In Histogram",freq=F,

    xlim=c(0.8,2.0),col='gray',xlab="Time")

curve(dnorm(x,mean(ins),sd(ins)),add=T,col='red')

legend(1.55,2.4,legend="Theoretical Normal",col="red",lty=1,

cex=0.75)

text(1.75,1.9,expression(paste(mu,"=1.350")))

text(1.75,1.75,expression(paste(sigma,"=0.164")))

dev.off()

png(paste0(path,"hist out.png"))

hist(outs,main="Breathe Out Histogram",freq=F,

xlim=c(1.5,3.5),col='gray',xlab="Time")

curve(dnorm(x,mean(outs),sd(outs)),add=T,col='red')

legend(2.7,1.45,legend="Theoretical Normal",col="red",lty=1,

cex=0.75)

text(3.05,1.18,expression(paste(mu,"=2.380")))

text(3.05,1.10,expression(paste(sigma,"=0.266")))

dev.off()

The data seem to follow a normal distribution. For this example, it isn't all that necessary or useful, but I find it interesting anyway.

forecast:

#forecasting

library(astsa)

acf(breathe.dat,100)

q=15

pacf(breathe.dat,30)

p=4

fore=sarima.for(breathe.dat,10,p,1,q)

# next 10 breaths

png(paste0(out,"forecast.png"))

upbound=fore$pred+fore$se*1.96

lowbound=fore$pred-fore$se*1.96

plot(breathe.dat[(n-20):n],type='b',xlim=c(1,21+length(fore$pred)),

ylim=c(min(lowbound)-0.1,max(upbound)+0.1),

main="10 Breath Forecast")

polygon(x=c(22:(21+length(upbound)),(21+length(lowbound)):22),

y=c(upbound,lowbound[length(lowbound):1]),col='blue')

points(21:22,c(breathe.dat[n],fore$pred[1]),type='b',col='red')

points(22:(21+length(fore$pred)),fore$pred,type='b',col='red',lwd=2)

dev.off()

All that work for something very useless, but there you go.