color cube

Just having some fun with 3-D plotting

setup.colors=function(detail=50){
    out=matrix(NA,detail^3,3)
    row=0
    for (red in seq(0,1,1/(detail-1))){
        for (blue in seq(0,1,1/(detail-1))){
            for (green in seq(0,1,1/(detail-1))){
                row=row+1
                out[row,1]=red
                out[row,2]=blue
                out[row,3]=green
                }
            }
        }
    return(out)
    }

color.plot=function(palette,fixed.limits=TRUE){
    require(rgl)
    n=dim(palette)[1]
    if (fixed.limits){
        plot3d(palette[1:n,],col=rgb(palette[1:n,1],
            palette[1:n,2],
            palette[1:n,3]),
            xlim=c(-0.01,1.01),ylim=c(-0.01,1.01),
            zlim=c(-0.01,1.01))
        } else {
            plot3d(palette[1:n,],col=rgb(palette[1:n,1],
                palette[1:n,
                palette[1:n,
                }
    }

colors=setup.colors(100)
color.plot(colors)

It's cool because you can even spin it around and look at it from other angles.



This uses the regular plot function to look at the color intensity growing from 0 to 1.

### changing Green
for (i in seq(0,1,1/(100-1))){
    plot(colors[colors[,2]==i,c(1,3)],col=rgb(
        colors[colors[,2]==i,
        colors[colors[,2]==i,2],
        colors[colors[,2]==i,3]),
        xlim=c(-0.01,1.01),ylim=c(-0.01,1.01),lwd=2,pch=20,
        xlab="Red",ylab="Blue",main=paste0("Green: ",i))

    }

factorization

Here's a recursive function that returns the factors of an integer:

factors=function(n){
n=abs(n)
out=NULL
new=0
if (floor(n) != ceiling(n))
return (factors(floor(n)))
for (i in 2:n){
if (n/i == floor(n/i)){
new=n/i
break
}
}
if (new==1){
return (i)
} else {
return (c(i, factors(new)))
} 
}

The n=abs(n) and return (factors(floor(n))) statements protect against negative and non-integer inputs. The idea is take a number n and divide it by the numbers 2 up to n until that division results in another whole number. If it does and that number is n itself then the function returns i=n, else there is another factorization to be made and so returns i with the next factor by the statement return (c(i, factors(new))).

> factors(223092870)
[1]  2  3  5  7 11 13 17 19 23
>
> factors(1073741824)
[1] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
>
> factors(1765469)
[1] 1765469
>
> system.time(factors(223092870))
   user  system elapsed 
   0.09    0.39    0.48 
>
> system.time(factors(1073741824))
   user  system elapsed 
   1.13    1.53    2.66 
>
> system.time(factors(1765469))
   user  system elapsed 
   2.51    0.00    2.51 


Improvements can be made in at least one place. After n/2 we know there aren't going to be any new factors, but the code currently iterates through all the values from n/2 to n. Making some slight changes can decrease the computation time significantly, especially for prime numbers.

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.