siddharth's space

How to Remote Desktop to Ubuntu:
follow instructions at:
https://aws.amazon.com/premiumsupport/knowledge-center/connect-to-linux-desktop-from-windows/
one bit of instruction is missing in above link and that has to do with enabling RDP traffic to your server (allow a client machine to communicate with server over RDP). You will need to add the RDP rule to the security policies. Find it in AWS console.

If you see a blank gray screen when you RDP to AWS:
http://askubuntu.com/questions/603314/gray-screen-when-i-try-to-remote-desktop-to-ubuntu

How to Flatten a folder in unix (very useful – need to do this every once in a while):
http://unix.stackexchange.com/questions/52814/flattening-a-nested-directory?newreg=4fb80858719b4a8e9fae5eebac6ded3b

How to Resize images using imagemagick:
http://stackoverflow.com/questions/10802606/how-to-batch-resize-images-in-ubuntu-recursively-within-the-terminal

find . -name "*.jpg" | xargs mogrify -resize 50%

above can run for a long time if you have thousands of files in the directory. In that case use this script:

for f in `find . -name "*.jpg"`
do
    echo resizing $f ...
    mogrify -resize 50% $f
done

this will also run for same amount of time but will display messages to console as it is resizing each file

Authorizing Inbound Traffic for Your Linux Instances:
http://docs.aws.amazon.com/AWSEC2/latest/UserGuide/authorizing-access-to-an-instance.html

How do I make `ls` show file sizes in megabytes?
http://unix.stackexchange.com/questions/64148/how-do-i-make-ls-show-file-sizes-in-megabytes

To see the disk space consumed by files in a directory:
$du -shb .
http://stackoverflow.com/questions/1241801/total-size-of-the-contents-of-all-the-files-in-a-directory

To send file as email attachment:

cat errors.txt | mailx me@email.com

http://stackoverflow.com/questions/17359/how-do-i-send-a-file-as-an-email-attachment-using-linux-command-line

To copy file contents to clipboard:
http://stackoverflow.com/questions/1620018/copy-all-the-lines-to-clipboard

Logging out:
When you are done with your session, its always a good practice to end a ssh session by typing $logout rather than closing the window. Similarly if you are using RDP, don’t just close the window to end the session. Logout of RDP using the GUI and you might want to uncheck the box that says save this session for future logons.

Deep Learning:

http://neuralnetworksanddeeplearning.com/

https://www.kaggle.com/c/facial-keypoints-detection/details/deep-learning-tutorial

Most important lesson of software development if you are a manager:

in case the point is not clear: First of all, the productivity does not increase in proportion to the # of devs i.e., if you double the devs, it won’t double the productivity and deliverables. More importantly, the graph does not reach an asymptote. After a certain point there is a decrease in productivity as more devs are added to a team because of various reasons:

• people conflicts and politics
• impossible for one person to keep track and abreast of every changeset that is happening. as a result, left hand does not know what right hand is doing

generate data:
x=rnorm(1000,mean=0)
y=rnorm(1000,mean=10)
obs1=data.frame(rbind(cbind(rnorm(1000),rnorm(1000)),cbind(rnorm(1000,mean=10),rnorm(1000,mean=10))),as.factor("Class A"))
obs2=data.frame(rbind(cbind(x,y),cbind(y,x)),as.factor("Class B"))
colnames(obs1) = c("x", "y", "class")
colnames(obs2) = c("x", "y", "class")
df=rbind(obs1,obs2)
make scatter plot:
#scatter plot
#dev.new()
png(file="scatter.png")
plot(obs1$x,obs1$y,col=colors[[1]],xlab="x",ylab="y",main="scatter plot")
points(obs2$x,obs2$y,col=colors[[2]])
dev.off()

see histograms:
for(i in 1:2)
{
#dev.new()
png(file=paste("hist-",names[[i]],".png",sep=""))
hist(obs1[,i],col=colors[[1]])
hist(obs2[,i],col=colors[[2]],add=TRUE)
legend(2000,9.5, c("Class A","Class B"))
dev.off()
}

build classifier:
model=glm(class~.,family="binomial",data=df)
response=predict(model,type="response")
M=NULL
for (threshold in seq(0.01,0.99,0.01))
{
M=cbind(M,calc(response,df[,3],levels(df[,3]),threshold))
}
M = cbind(c(1,1),M,c(0,0))
fp = M[1,]
tp = M[2,]
auc = -auc(fp, tp)
#dev.new()
png(file="roc.png")
plot(fp,tp,type="l",xlab="false positive",ylab="true positive",main="Logistic Regression classifier")
#lines(c(0,0),c(1,1))
dev.off()

full code:

colors = list(rgb(0,0,1,1/4),rgb(1,0,0,1/4))
x=rnorm(1000,mean=0)
y=rnorm(1000,mean=10)
obs1=data.frame(rbind(cbind(rnorm(1000),rnorm(1000)),cbind(rnorm(1000,mean=10),rnorm(1000,mean=10))),as.factor("Class A"))
obs2=data.frame(rbind(cbind(x,y),cbind(y,x)),as.factor("Class B"))
colnames(obs1) = c("x", "y", "class")
colnames(obs2) = c("x", "y", "class")
df=rbind(obs1,obs2)
#scatter plot
#dev.new()
png(file="scatter.png")
plot(obs1$x,obs1$y,col=colors[[1]],xlab="x",ylab="y",main="scatter plot")
points(obs2$x,obs2$y,col=colors[[2]])
dev.off()

Dp=D(as.matrix(obs1[,1:2]),as.matrix(obs2[,1:2]))
names=c("x","y")
for(i in 1:2)
{
	#dev.new()
	png(file=paste("hist-",names[[i]],".png",sep=""))
	hist(obs1[,i],col=colors[[1]])
	hist(obs2[,i],col=colors[[2]],add=TRUE)	
	legend(2000,9.5, c("Class A","Class B"))
	dev.off()
}

model=glm(class~.,family="binomial",data=df)
response=predict(model,type="response")
M=NULL
for (threshold in seq(0.01,0.99,0.01))
{
	M=cbind(M,calc(response,df[,3],levels(df[,3]),threshold))
}
M = cbind(c(1,1),M,c(0,0))	
fp = M[1,]
tp = M[2,]
auc = -auc(fp, tp)
#dev.new()
png(file="roc.png")
plot(fp,tp,type="l",xlab="false positive",ylab="true positive",main="Logistic Regression classifier")
#lines(c(0,0),c(1,1))
dev.off()

create training data:

n=10000
x=runif(n)
y=runif(n)
df=data.frame(x,y,as.factor(x>y))

visualize:

colors = list(rgb(0,0,1,1/4),rgb(1,0,0,1/4))
names(df) = c("x", "y", "class")
classes=levels(df[[3]])
obs1=as.matrix(subset(df,class==classes[[1]])[,1:2])
obs2=as.matrix(subset(df,class==classes[[2]])[,1:2])
# make scatter plot
dev.new()
plot(obs1[,1],obs1[,2],xlab="x",ylab="y",main="scatter plot",pch=0,col=colors[[1]])
points(obs2[,1],obs2[,2],xlab="x",ylab="y",main="scatter plot",pch=1,col=colors[[2]])

result:

D = \frac{|\mu_1-\mu_2|}{\sigma_1+\sigma_2}:

> Dp
        x         y 
0.7121875 0.7208857 

histograms of x and y:

build classifier:

model=glm(class~.,family="binomial",data=df)
summary(model) # prints summary

Output:

Call:
glm(formula = class ~ ., family = "binomial", data = df)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-0.11832   0.00000   0.00000   0.00000   0.08847  

Coefficients:
              Estimate Std. Error z value Pr(>|z|)
(Intercept)  5.765e-01  1.923e+01   0.030    0.976
feature1     9.761e+04  8.981e+04   1.087    0.277
feature2    -9.761e+04  8.981e+04  -1.087    0.277

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1.3863e+04  on 9999  degrees of freedom
Residual deviance: 2.9418e-02  on 9997  degrees of freedom
AIC: 6.0294

Number of Fisher Scoring iterations: 25

algorithm fails to converge because of perfect separation:

Warning messages:
1: glm.fit: algorithm did not converge 
2: glm.fit: fitted probabilities numerically 0 or 1 occurred 

compute AUC:

response=predict(model,type="response")
groundTruth=df[,3]
M=NULL
for (threshold in seq(0.05,0.95,0.05))
{	
	M=cbind(M,calc(response,groundTruth,classes,threshold))	
}
M = cbind(c(1,1),M,c(0,0))	
fp = M[1,]
tp = M[2,]
auc = -auc(fp, tp)

where:

auc=function(x,y)
{
	d = diff(x, 1)
	e = as.numeric(lapply(1:(length(y)-1), function(i)y[i]+y[i+1]))	
	return(sum(0.5*d*e))
}

calc=function(response, groundTruth, classes, threshold)
{
	type1 = classes[1]
	type2 = classes[2]
	n=length(response)
	predicted = as.factor(ifelse(response<threshold,type1,type2))
	I = which(groundTruth==type1)
	fp = length(which(predicted[I] != type1)) / length(I)
	I = which(groundTruth==type2)
	tp = length(which(predicted[I] == type2)) / length(I)
	return(c(fp,tp))
}

plot AUC:

plot(fp,tp,type="l",xlab="false positive",ylab="true positive",main="Logistic Regression classifier")

result:

LDA:

lda_classifier=lda(class~., data=df)
Call:
lda(class ~ ., data = df)

Prior probabilities of groups:
 FALSE   TRUE 
0.5007 0.4993 

Group means:
       feature1  feature2
FALSE 0.3346288 0.6676169
TRUE  0.6710111 0.3380432

Coefficients of linear discriminants:
               LD1
**feature1  4.280490
feature2 -4.196388**

finally you add add some noise to make logistic regression converge:

df=data.frame(x+rnorm(n,sd=0.1),y+rnorm(n,sd=0.1),as.factor(x>y))

scatter plot:

ROC:

Here is a function to compute (false positive, true positive) pair given response, ground truth, classes and threshold:

# response is e.g. predict(model,type="response") 
calc2=function(response, groundTruth, classes, threshold)
{	
	type1 = classes[1]
	type2 = classes[2]
	n=length(response)
	tpden = 0
	fpden = 0
	fp = 0
	tp = 0
	for(i in 1:n)
	{	
		predicted = ifelse(response[i]<threshold,type1,type2)
		actual = groundTruth[i]
		if (actual == type1)
		{
			fpden = fpden + 1
			if (predicted != actual)
			{
				fp = fp + 1
			}
		}
		else
		{
			tpden = tpden + 1
			if (predicted == actual)
			{
				tp = tp + 1
			}
		}		
	}
	fp = fp/fpden
	tp = tp/tpden		
	return(c(fp,tp))
}

This function is 100x slower than below which does the same thing:

calc=function(response, groundTruth, classes, threshold)
{
	type1 = classes[1]
	type2 = classes[2]
	n=length(response)
	predicted = as.factor(ifelse(response<threshold,type1,type2))
	I = which(groundTruth==type1)
	fp = length(which(predicted[I] != type1)) / length(I)
	I = which(groundTruth==type2)
	tp = length(which(predicted[I] == type2)) / length(I)
	return(c(fp,tp))
}

To benchmark, install rbenchmark package and use it like below:

> benchmark(calc(response,groundTruth,classes,0.5), calc2(response,groundTruth,classes,0.5),replications=10)
                                              test replications elapsed relative user.self sys.self user.child sys.child
1  calc(response, groundTruth, classes, 0.5)           10    0.14    1.000      0.14        0         NA        NA
2 calc2(response, groundTruth, classes, 0.5)           10   13.06   <strong>93.286</strong>     13.05        0         NA        NA

getwd to get current working directory (analogous to pwd)
setwd to set current working directory
options(digits=3) will print numbers with 3 decimal digits
runif(20) – create vector of 20 uniformly distributed random numbers
dir – to list files in a directory
dir.create – to create a directory
source(‘myscript.R’) – to execute a script
install.packages(Rwave) – to install Rwave package
update.packages() to update packages
installed.packages() shows details of installed packages
help(package=’Rwave’) will list functions in package Rwave
library(Rwave) – to load Rwave package
library() – shows what packages you have saved in your library
search() shows which packages are loaded and ready for use
.libPaths() shows where your library is located e.g. C:/Program Files/R/R-3.1.1/library
sink(‘filename’) will direct output to filename
sink() will restore output to screen
pdf(figure.pdf) save figure as pdf
png(figure.png) save figure as png
length(x) returns length of x
dim(x) returns dimensions of x (should be matrix or array)
str(x) tells structure of x
class(x) tells what class is x type of
mode(x) tells how x is stored
names(x) gives names of components (e.g., columns) of x
5 main types of collections in R – vectors, matrices, arrays, data frames and lists
c to create vectors
matrix to create matrix
array to create array (matrix with more than 2 dimensions) e.g., array(1:24,c(2,3,4)) creates 2x3x4 array
data frame is a matrix whose columns can have different data types df=data.frame(col1,col2,col3)
df[1:2] will return first 2 columns of data frame
the names argument of data.frame ctor can be used to give names to the various columns. We can then index by column names instead of having to use integers to index.
use table function to cross tabulate 2 columns of a data frame e.g., table(df$col1,df$col2)
use attach and detach to avoid having to prefix name of data frame again and again
also see with which provides similar functionality
Nominal and ordinal variables in R are called factors (think of enums in other languages). Both are created using factor function. Using ordered=TRUE in the ctor tells R if you are creating ordinal variables (whose values can be ordered). The order can also be specified using the levels param e.g.,

factor(data,levels=c(1,2),labels=c("male","female"))

is basically saying that data is an enum where {male=1,female=2}
Lists are well lists. they are indexed using double square brackets [[ and ]]
to read data from a flat file use read.table

df=read.table(file,header=TRUE,sep=",",row.names="name")

header indicates whether the first record in file contains the column headings
sep is separator used to separate columns
row.names specifies the columns that can be used to ID a record (sort of like PK in SQL)
dev.new() similar to figure in matlab (create new figure)
more commends: dev.next, dev.prev, dev.set, dev.off
par(mfrow=c(m,n)) to create a figure with m by n plots in it
graphics.off() to close all graphs
plot commands starts a new graph
lines command adds to current plot
abline can be used to draw vertical and horizontal reference lines on a plaot
x%%y same as x%y in C# 5%%2 = 1
x%/%y returns quotient 5%/%2=2
!= not equal to
== equal to
!x not x
x|y x OR y
x&y x AND y
isTRUE(x) test if x is TRUE
cut allows you to divide a numeric range into intervals returning a factor. These intervals can then be used with tapply to do binning.
Sys.Date() returns today’s date
date() returns current date and time
order() to sort data
which gives indices of a vector that are TRUE
sample() function can be used to do sampling
diff(x,lag=1) will return difference of successive elements
to apply a function on a particular column, you can use transform() e.g.

newData=transform(myData,myVar=scale(myVar))

seq(to,from,by) e.g. seq(0,10,2) returns 0,2,4,6,8,10
lapply and sapply apply a function over a list
to measure cross correlation between signals use

ccf(x, y, lag.max = NULL, type = c("correlation", "covariance"), plot = TRUE, na.action = na.fail, ...)

to measure just the correlation, use cor
see this answer to find the time at which maxima occurs in cross correlation

Find_Max_CCF<- function(a,b)
{
 d <- ccf(a, b, plot = FALSE)
 cor = d$acf[,,1] // gives cross correlation function
 lag = d$lag[,,1] // gives vector of time lags
 res = data.frame(cor,lag)
 res_max = res[which.max(res$cor),] // gives time lag at which max cc occurs
 return(res_max)
} 

to apply a circular shift to a vector, i came across this code:

<pre>circshift <- function(a, sz) {
    if (is.null(a)) return(a)
    
    if (is.vector(a) && length(sz) == 1) {
        n <- length(a)
        s <- sz %% n
        if (s > 0) a <- a[c((s+1):n, 1:s)]

    } else if (is.matrix(a) && length(sz) == 2) {
        n <- nrow(a); m <- ncol(a)
	    s1 <- sz[1] %% n
	    s2 <- sz[2] %% m
        i1 <- if (s1 > 0) c((s1+1):n, 1:s1) else c(1:n)
        i2 <- if (s2 > 0) c((s2+1):m, 1:s2) else c(1:m)
	    a <- a[i1, i2]

    } else
        stop("Length of 'sz' must be equal to the no. of dimensions of 'a'.")

	return(a)
}</pre>

to convert a list of vectors into a matrix use this:

M=do.call(rbind,lapply(data$M,function(x) extract_features_internal(x,fs)))

to extract filename from full path use basename(fullPath)
to compute average correlation between columns of a matrix use:
C = cor(M) # computes correlation between columns of M
correlation = mean(C[lower.tri(C)]) # get average correlation

use proc.time() to get time
See this link http://www.ltrr.arizona.edu/~dmeko/notes_6.pdf to really understand what goes behind the scenes when you call spec.pgram or spec.ar or spectrum
Use drop to Delete the dimensions of an array which have only one level.
To find indices of top N elements in a list or vector use this:

order(R, decreasing=TRUE)[1:N]

To filter a data frame on values of its column use this:

obj = subset(table,column_name==class)

above will return those rows of table where column_name equals class

http://jsfiddle.net/siddjain/wwLCD/
http://jsfiddle.net/siddjain/Hq63y/1/