title: How to analyze this data using rjags, or any other way?

text: There are three groups that received different treatments and then learn a task which they are scored on over the course of 60 sessions. 1) **How to choose a function to fit?** -The best choice in this case would be `s=asym/(1+exp((midpt-x)/rate))`, but say you didn't know that. 2) **How to compare the three groups to see if they are different or not?** -For example Group 1 is better/slower than Group 2. 3) **Say there is a third factor measured (e.g., weight) how to look for correlation between weight and the curve fit for each individual?** -For example use the mean estimates. ![enter image description here](http://i.stack.imgur.com/MiiLK.png) **R Code to Generate Data:**               x<-1:60     n<-5     dat=NULL     par(mfrow=c(4,4))     for(g in 1:3){     for(subj in 1:n){     asym<-rnorm(1,10,3)     midpt<-runif(1,5,40)     rate<-runif(1,0,5)     s=asym/(1+exp((midpt-x)/rate))     s<-rnorm(length(s),s,1)     s[which(s<0)]<-0     dat<-rbind(dat,cbind(g,subj,x,s))     plot(x,s, type="p", pch=16, col=rainbow(3)[g], xlab="Time",ylab="Score",     main=paste("Group",g,":","Subject#",subj)     )     }     }     colnames(dat)<-c("Group","Subject","Time","Score")      **The data for the plot:**               > dput(dat)     structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,      1, 1, 1, 1, 1, 1, 1, 1, 1, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 3, 3, 3, 3, 3,      3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,      3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,      3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,      3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,      3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,      3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,      3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,      3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,      3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,      3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,      3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,      3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,      3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,      3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,      3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 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, 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, 3, 3, 3, 3, 3,      3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,      3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,      3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4,      4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,      4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,      4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,      5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,      5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,      5, 5, 5, 5, 5, 5, 5, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,      1, 1, 1, 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, 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, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,      3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,      3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4,      4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,      4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,      4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5,      5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,      5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,      5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 1, 1, 1, 1, 1, 1, 1, 1,      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,      1, 1, 1, 1, 1, 1, 1, 1, 1, 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, 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, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,      3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,      3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,      3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,      4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,      4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,      4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,      5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,      5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 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, 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, 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, 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, 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, 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, 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, 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,      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, 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, 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, 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, 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, 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, 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, 0, 1.03158619583643, 0.0700928693799847,      0, 0, 0.533480242356568, 0, 0, 0.313246505626844, 0, 0.201551652833782,      0, 0, 8.58420236367213, 7.2151582207413, 6.76349698933596, 8.43232503189182,      7.84113269530773, 7.48036740442613, 7.22860126336753, 6.23889089933047,      7.5993822730141, 8.97711255338355, 8.51408365480236, 7.33679466179734,      8.07177922774387, 7.46193250657315, 6.69880073822405, 8.7239879783316,      8.63777777133611, 7.32811771015052, 8.21150317871178, 6.75210022817865,      9.19584687427251, 8.84904099669932, 8.12546344066828, 7.42498538242589,      7.43402031096714, 7.95799796160173, 6.75415984451499, 8.41849594339,      6.45648891349807, 8.36526948269262, 6.73507595185595, 9.59230196004733,      6.8481757733426, 7.5336436175937, 7.26145341491774, 6.84374696765792,      5.77076068283607, 9.84247282482329, 7.78925766764691, 6.50711775948628,      8.19902246920613, 8.01867704068551, 7.71330298589825, 8.50368409962827,      8.36500008512313, 8.37090214729723, 7.28633938219456, 0, 1.65008045947966,      0.575606733495696, 0, 0.786060958385394, 0, 1.54869682228415,      0.168946971854925, 0, 0, 0, 0, 0, 0.213101769560987, 0.440700799230945,      0, 0, 0.190648393014687, 0, 0.774590152871335, 0, 0, 0, 0, 0.178064415218869,      0, 0, 0.376138177422834, 0.458985077565505, 0.917398964689014,      1.78416461506729, 1.37112772722134, 2.63275499887878, 0.247691262888474,      2.71308742063024, 3.41901470336746, 4.11919704486841, 4.04233691698445,      4.97741716048949, 4.5916456213906, 3.74309103176392, 6.25015151641352,      5.01570149643103, 5.02144478232446, 6.44037624126331, 4.36299946803298,      4.81397614036387, 4.21616606222475, 3.94904834726585, 6.2725119824346,      4.96221630092668, 4.85059793371767, 6.98468854857044, 6.13751640244375,      5.31975082309877, 4.95974732120037, 4.96944203644392, 5.04695972379456,      5.55090724593748, 4.9072728753408, 0, 0, 1.27148334131862, 0,      0, 0, 0, 0.549249296594036, 0, 0, 0, 0, 0.797894117158937, 0.76474221921926,      0.384456323004612, 0.753724350710486, 0.210172266893594, 0.778029999658896,      0, 1.20840679374069, 2.27515993940812, 2.33021789487285, 2.48682959545432,      2.77439234964342, 2.57654149782979, 3.6899693083328, 2.84497438582907,      4.53576090191172, 5.31904443998464, 5.83484731859952, 4.16995891072675,      7.46357804432578, 8.0189267920945, 10.1186546998081, 9.17802346129456,      9.91280676028886, 8.6358637300104, 10.6171626191218, 10.60862571782,      8.90168548865334, 11.6613913866556, 11.494344312576, 9.83605794475883,      10.5086341235427, 8.57098156695857, 10.6578914329804, 11.3242140265126,      10.6128797597814, 8.93858406115906, 11.442477317725, 10.2773057203746,      11.422141318439, 10.4578041679639, 11.2805966086809, 10.9729630215866,      10.6400309657672, 10.5325561935863, 10.3806803061043, 12.2462712098331,      11.1739112019198, 0, 1.37519992104938, 0.00230181625672307, 0.498711700484401,      0.271234286784188, 0, 0.546777568628416, 1.02624761127694, 0,      0.802875563963303, 0, 0.75661766368099, 0.859889392170508, 0.771882758754063,      2.68818858778084, 0, 1.69857213256685, 0, 0, 0, 0.509209770872453,      0, 0.75119289550595, 0, 0.792697602543835, 0.99601355947128,      0.615579450984568, 2.25557535234532, 1.14413794471848, 0.111971817227507,      1.29934159296573, 0.109732809889403, 1.44834021140093, 4.61542498626519,      5.20534003980441, 7.58289763660421, 7.05378158533781, 6.48087001005784,      7.46752691991914, 7.06731997153742, 8.43297260855126, 7.26972090539512,      7.38136834636119, 7.40544432675767, 7.75642154862361, 6.76324320314802,      5.59456019542772, 7.40022245947996, 7.48835520637033, 7.65877703591718,      8.21148150997923, 6.6599237915552, 8.0160411529529, 6.60139606432752,      8.65128003881004, 5.79775256308877, 6.85878904841619, 7.84674726214019,      4.3467120746636, 7.55232892761264, 0, 0.499583639143117, 0, 0,      0.31485702801172, 0.3515816928799, 0.931603213795884, 0.206120594561129,      0, 0, 0.360932275745327, 0.433044022640932, 0, 2.39516067822733,      2.72351140421486, 0.484924178938609, 0, 0.224022539659828, 0.911042602837107,      1.31501526265012, 1.94337905383981, 0.595289303995586, 2.41167836489199,      4.88047662904194, 3.60477176392402, 4.70295731707599, 5.46521105212445,      5.91105387042816, 7.92754559701199, 6.30579526052163, 7.94051926136399,      8.93734485763105, 8.10797265204026, 7.47519696905344, 7.87790798680286,      8.37458434707145, 9.88966658479798, 9.3828993940551, 10.7502415658999,      8.59224045217057, 8.97257663583808, 9.81078210053123, 8.50519146337723,      9.80410794437355, 9.90656980648206, 8.95149994048142, 8.08214817120231,      9.15819742292325, 9.35229009595692, 8.66137979115357, 7.61365572420257,      10.1643174026938, 10.0183955498814, 7.4284240739572, 8.70163824293303,      8.97812942423644, 8.40596894398334, 8.9355272887446, 10.7662692661888,      10.0322226889031, 0, 1.17127841967856, 0, 0, 1.33105575939346,      1.26518317997976, 0.909776451243087, 0, 1.17988956763035, 0,      0.762002789890116, 0, 0.832083512904214, 0.45869046938725, 0.64237344364855,      0, 0, 1.70209058186332, 0.0156597425215317, 0.530568921018249,      1.50289220053205, 1.20254177632711, 1.90295146480201, 2.81271374240921,      3.68747746746046, 4.52706950129025, 3.93442281955956, 2.5515226646154,      5.01710413340776, 6.26260874497632, 5.49950756302542, 6.14255023572,      7.64342305367938, 9.89655193964118, 8.77398103235339, 7.64631421131576,      7.77582538877235, 8.30983971196275, 7.8411138994189, 8.23146449905929,      8.78847752071418, 7.98107466126393, 10.070698869224, 11.4707954191209,      10.7802260341204, 9.57883214145681, 9.45484978547647, 9.48982907075886,      9.1970407665382, 11.6749555300959, 9.63111689744601, 8.58499382409839,      11.1715368888387, 9.05676905125274, 9.62946264711658, 9.58942799488546,      9.85519469529828, 11.3423886432642, 10.2388217496706, 9.321088820638,      2.54991829591086, 1.46892808704928, 1.33731418668811, 2.61646582984654,      2.66441558379904, 1.94155905469566, 2.04033379550069, 2.48570760395896,      4.42042378981082, 7.39539169831503, 7.13750763644472, 7.56179540351745,      8.3182611498041, 9.08500888550381, 9.70808482796405, 7.73515608944475,      8.84736021599179, 9.91632970815538, 6.85549285375218, 9.28524174344136,      9.90944985145035, 9.80872079378178, 7.50153753628047, 9.43229031930222,      8.38966493676363, 9.09120122077017, 8.03664950919612, 9.3394983271658,      8.75907130568531, 9.48293724329285, 12.4992084970297, 8.65195448377916,      8.12867324807331, 9.35139039060214, 7.45991827498198, 7.99054707729177,      9.33353246911462, 10.0140575953154, 9.62351889572776, 8.24274614700205,      8.58318657774853, 10.7558360231622, 8.02327882376913, 9.71166452940178,      10.0019440376244, 9.45855647023341, 9.69119327483915, 10.7993477282601,      10.0843411976038, 9.11504744632914, 6.8274470318345, 8.99834673719741,      7.79205423949022, 9.36012575556776, 7.89124643127951, 9.30281848417072,      8.74860790084056, 9.5283366549812, 9.38344434187134, 10.2903751565271,      0.453745586179272, 0, 0, 0, 0.873987314878141, 0, 0.722662147827297,      0.403389511042975, 0.306843142235151, 1.02980154677912, 1.26162245184229,      0, 0.738903729490879, 1.09034307241339, 0, 1.58160726939009,      1.00967640268439, 0, 0.80185808962151, 2.32319612854146, 0.680322447550744,      3.45272401214232, 3.98795894572399, 3.07813212524925, 2.74649863082109,      6.22716830523846, 6.0055534600751, 8.61611296275224, 7.43196060556798,      10.8845246381959, 9.35652157761648, 13.535727140111, 11.8228054818836,      13.6576101040097, 14.5985384582905, 14.720679614782, 16.4464360616433,      15.650288548452, 17.1098996963077, 15.9900299666647, 15.1733513084352,      13.460940246974, 17.1560722762101, 16.3594343339547, 17.6334893606005,      15.4742543462313, 17.9252909798281, 16.0685717969627, 16.5922216823379,      17.4019406800853, 16.4827827623591, 16.2283964030323, 17.3342536279416,      18.0418680171866, 16.7384801226295, 16.6891351719263, 16.9196389730558,      17.4360442891675, 17.9082826953212, 17.9080493357515, 0, 2.12708907689137,      0.180057921261079, 1.10881706586506, 0, 0.0121433242856085, 0.297984364375141,      0, 0, 0, 2.44423336550259, 0.203063335897161, 0, 0, 1.6324104596689,      2.62199381849, 4.65386129828348, 8.1727461408628, 8.46401772635808,      9.53664285078371, 10.6167429148595, 14.1642559994959, 14.1737553430174,      11.7871799896131, 11.6842610420266, 13.3115074910252, 12.542611801357,      11.7928591701523, 13.3234935331678, 13.3703045475931, 11.1816912869155,      12.7919758124473, 14.1524237174757, 13.7740194661304, 13.4163111242968,      12.1519239667949, 12.3411763952143, 14.1146449302467, 10.9799552937331,      11.9519542581919, 13.5564141397704, 12.6603319807821, 13.4818352225032,      12.4199788441116, 12.8984348945313, 12.4901614053039, 12.6826270147467,      12.8972689922525, 12.8761297786996, 14.5082862520908, 13.3088273305523,      11.689683656794, 11.6905619013621, 12.9372150020277, 11.5355247094442,      13.2387324569227, 13.2672299010805, 13.8721895345878, 13.6468900741745,      10.2407254663892, 0, 0, 0, 0.647542861902527, 0, 0, 0.150952145577755,      0.957392963769722, 0.0609913610066249, 0.120528953161183, 0,      0, 0, 0.704879442804958, 1.78408096310525, 0, 1.65023390832557,      0, 0.0923652900702186, 0, 0, 0, 0.0450577278912485, 0.751091209975893,      1.43744886707744, 0, 0, 0.402523264684891, 0, 2.18009214813759,      0.771506994461835, 0.366208114547776, 2.94463744259694, 3.23887911393408,      1.32442373753528, 2.39016677599377, 2.11611119824864, 4.4636528481399,      6.31011316007659, 5.93705422891922, 6.86234582501057, 7.04191497749316,      8.65532314557916, 7.94975578786001, 8.32465895135107, 10.2544694107277,      10.612243226875, 11.0923080512734, 11.330145953257, 10.9264546212034,      11.2641173880016, 13.725651001296, 12.4933952691916, 11.6957249070898,      9.83082658829007, 10.1737819351518, 10.0278938736921, 10.2885620365588,      11.4721405494207, 10.4049091351876, 0, 0.749971887045118, 0.0935661865127177,      0, 0.276694847632314, 2.16959186220499, 1.03997487041489, 1.3157745452179,      2.55444161550264, 5.55564346433049, 6.5606015783716, 6.26361090505197,      8.21561991560817, 8.18272868467541, 10.3969243810899, 11.2467666988384,      11.9401702885542, 11.4225328143267, 12.0539674679589, 12.0500044044004,      11.9905170877494, 10.2276294132972, 13.8701081425486, 12.2345025406528,      12.1581776929614, 12.2557621066668, 14.8327853422983, 13.5706592392689,      11.670858423701, 12.1264565841421, 10.9882135554883, 11.910278334423,      12.2622517792821, 11.1803535948265, 12.3218944257386, 12.4022799379156,      13.3489881438595, 12.1519181213977, 11.8546766797642, 12.803931564254,      12.0702080416288, 11.7787770410977, 12.891356400286, 13.5826348158223,      13.7637998730129, 12.4586804497261, 13.5973025845184, 11.7437527504006,      11.9123359201718, 12.5433135233633, 11.6795905695007, 12.8596795881007,      12.2905970480821, 12.6113804745829, 11.2696221922899, 11.6150967135138,      12.5065320547622, 11.0547024118768, 12.4539473082803, 12.7320989917665,      0.828387327199212, 0, 0.678802434035732, 0.862374178714765, 1.2684255030382,      0, 1.83576620899006, 4.71180773591447, 7.15854703611923, 9.98309590659522,      12.238244011289, 10.1107587439567, 11.9735917806089, 10.4329498397988,      11.5227705911928, 11.0857295162145, 12.2887936081285, 11.5414959841271,      13.0408942545176, 12.1105317096253, 12.278626227706, 12.928520689136,      12.7511178109053, 12.8603069946746, 12.0171195930741, 10.8521389192071,      11.1200911135579, 12.6196298934609, 13.0379900612457, 11.9235143903274,      10.5663812910848, 9.58631186189439, 12.3408962856898, 11.8850627357227,      11.4640087047838, 11.6422475778865, 11.7972751414621, 12.4895822567521,      12.5922359409969, 11.3441135167686, 11.1062090592898, 13.5286995172176,      10.9736620416305, 9.97146752611849, 12.1805784426637, 10.4799209638073,      12.000342774528, 10.0789727398843, 11.8248401081463, 10.9804844946854,      11.835301869748, 12.6690753183117, 13.0419456703609, 10.5160699826278,      11.0290656985927, 11.2908514684923, 13.2547962666313, 11.716640532936,      10.8753781493541, 12.8659989966435, 0.233214052916719, 0, 0,      0, 0.666797775504629, 1.03113449232196, 0, 0, 0.235911824664366,      0.643988682300099, 2.35216791250754, 4.48251004203957, 6.27829837596023,      7.80685682111985, 8.85867907137818, 11.0234994200981, 10.659303622617,      8.83821767964211, 10.6285040270975, 11.9610855075861, 9.98850567573665,      9.38183764798606, 9.12964903103025, 9.79518026168435, 10.2486037764971,      10.7399919605559, 10.4000561384545, 10.1195650266721, 10.8172633161224,      11.7290460025207, 9.55500863048794, 9.51997696957129, 8.136245145294,      10.0898702564584, 8.83124284475855, 9.53702720995183, 10.234965256881,      11.8832880596309, 11.0260846401941, 10.3616281489605, 11.1595975463324,      10.6065362312525, 11.7789145668303, 10.7391985633112, 8.7345559786342,      9.83721824425054, 9.87969316669263, 11.0054313800222, 9.45630790909011,      10.8026069853639, 10.7823164169656, 8.78222450286674, 9.88401914578533,      9.49819408578393, 11.1070575565997, 9.80251606293417, 11.6800557898553,      10.6628812540799, 10.3481870385697, 8.99926635802634, 2.0658714541966,      0, 0, 0, 0.643696068708256, 1.95882560158578, 1.46771793144107,      0, 0, 0, 0, 1.57807910644584, 0, 3.15566724378049, 4.15202850585111,      4.49986707099665, 6.29955267584495, 8.11631515638246, 8.6413891325569,      8.63053991475974, 8.02359584421251, 9.22297649026951, 8.5220554417732,      7.83394559530608, 8.05150333895255, 8.81497999527858, 8.32251493454089,      8.68646050203423, 8.30964562140176, 6.24790293399363, 7.50865404587946,      7.92394725060917, 7.85279621262117, 9.79363556954982, 8.33078748739039,      11.0010688083092, 7.87335403086467, 8.04400683453104, 8.87349257925123,      5.34416893325697, 7.16402027526085, 9.18892534710063, 7.37307377041812,      9.12042575487411, 7.13803128469491, 9.53932820803311, 7.2471488352217,      8.13462753251326, 9.32994365735403, 8.9143655639326, 8.68129951409304,      8.06584003431567, 7.00267444564897, 10.8518954588216, 8.09952530895893,      9.52528768838161, 8.0767001453374, 7.54472269908506, 8.2017301034284,      8.93945165288032, 0.664837492834661, 1.42967883870151, 1.21461210510671,      1.52247428784209, 1.5857475270794, 4.21788588756755, 4.07472703605436,      4.52783139621065, 4.61304565650941, 4.97223959138618, 7.17290663409386,      7.55151237688616, 7.23137033191666, 5.6507989759957, 9.61894330860323,      8.44955835598088, 7.7741137106094, 7.07096361575731, 9.79144128185851,      10.5769775728822, 9.65529577781692, 9.13039573776741, 9.25919197691465,      9.45335765533184, 8.9662422782295, 9.18219296678602, 9.47685142345247,      10.0756859334663, 10.0817513792461, 11.0084378809858, 9.78118846695106,      9.76090117884639, 8.92844816510047, 8.16642032644479, 10.0468897494144,      9.04894852060757, 10.9828155120775, 9.93150745149563, 9.30769292314918,      11.0035921724632, 12.2009103535073, 11.6412800632603, 11.2191465078207,      10.2601299220248, 10.3150450691667, 9.14815709751079, 10.1162531475645,      8.50753892237624, 9.53758392124005, 10.8216774355975, 9.5098305251634,      10.4837560161139, 10.2617540070955, 10.038654859975, 9.00812936944552,      10.7869945837923, 9.91834237561623, 11.0188750966019, 11.2219462128563,      10.9428246709094), .Dim = c(900L, 4L), .Dimnames = list(NULL,          c("Group", "Subject", "Time", "Score")))


title: Limiting behavior of a martingale

text: This is a homework question: Suppose that $X_0=1$ and that for $n\geq 1$ $$X_n\sim \left\{ \begin{array}{l l} U(0,X_{n-1}) & \quad \text{with probability $1-X_{n-1}/2$}\\ U(X_{n-1},1) & \quad \text{with probability $X_{n-1}/2$} \end{array} \right.$$ Determine the limiting behavior of $X_n$. My attempt: 1) Check if this is a (sub/super)martingale 2) By the martingale convergence theorems, since it is bounded (below by 0, above by 1) this random variable must converge to something. 3) It can't possible converge to some value $x\in (0,1)$ since if $X_n=x$, then theres a high probability that $X_{n+1}$ is more than $\epsilon$ units away from $x$ for any $\epsilon>0$. So it should converge to either 0 or 1 (or both). 4) My guess is $X_n \overset{a.s.}{\longrightarrow} X\sim \text{bernoulli}(p=1/2)$. A quick simulation however suggests that it converges almost surely to just 0!               quiz5_number2 <- function(reps=1000, path=NULL){       x0<-1       path<-c(path,x0)       for(i in 1:reps){         y<-rbinom(1, size=1, prob= 1-x0/2)          x1<- runif(1, min=0,max=x0)*y + runif(1, min=x0,max=1)*(1-y)         x0<-x1         path<-c(path,x0)       }       path     }     path<-quiz5_number2(100)     plot(path)      ![enter image description here](http://i.stack.imgur.com/th4TL.png) How would you answer this problem? I realize the problem seems poorly worded. Blame my teacher :) EDIT: On further analysis it is quite simple to see that while 0 is an 'absorbing' state (if you get to zero you will stay there), 1 is not an absorbing state (if you are at 1 you have 0.5 chance to be drawn from uniform(0,1) and 0.5 chance to stay). So I guess it makes a lot of sense that it converges to 0. Thank you for reading.


title: Adjusting any power analysis with FPC?

text: Is it correct (perhaps at least as a good approximation) to determine sample size via a power analysis that assumes an infinite population and then to adjust the required sample size by a finite population correction? For example, I am using the formula from Fleiss: ![enter image description here](http://i.stack.imgur.com/x3SOQ.jpg) To calculate the sample size required to test the difference between two independent proportions, with unequal sample sizes. If I get out of this formula that the two samples sizes are: n1= 300 n2= 600 (as r in the formula, the ratio n2/n1, was set to 2) can I then apply the FPC to the total sample size to correct for the finite population? ![enter image description here](http://i.stack.imgur.com/IYQ3i.jpg) For example, say the population size is really 1000. Then, can I apply the correction to the total sample size (300+600=900) and recalculate the sample sizes from the two populations as: new_total_sample = $\frac{900*1000}{900+(1000-1)}=474$ and then allocate between the two populations as new_n1 = 158 and new_n2 = 316? Thanks! **ADD After Dr. Lumley's response:** I. I am a little confused on the difference between a hypothesis test and a confidence interval for a finite population difference. Are you saying it is legitimate to estimate the confidence interval of say, the difference between the proportion of men and women answering 'Yes' but not a p-value to determine if it is zero (i.e. to reject the null)? In the language of your package, I am thinking of               design  <- svydesign(id=~1, fpc=~N,data=experimentData)     svyglm(respondedYES~as.factor(isMale), design=design, family=binomial)      where you would construct a CI around the coefficient for isMale II. Also, I understand I think, the issue with the example of how the proportion of men and women could never be the same given your lack of a common factor argument, but what about hypothesis testing with a finite population in the case of an experimental design (where the researcher controls the treatment assignment) - where this issue is not present? That is the realm I am trying to learn how (if) the fpc can be used.


title: What is the name of the survey method for determining the rank order of multiple items using a multiple question process?

text: I have seen surveys that use the following method to rank multiple items against each other. Let's say you have 25 items and you want to know their rank order (according to the survey-taker's preference). You pick 5 at random and ask the survey-taker to pick their favourite and least favourite of the group. You can then pick another 5 and repeat the process until the list is exhausted. You can then pick the favourites and rank them against each other and then do the least favourites and then the ones in the middle and through an iterative process, rank the entire list. I'm guessing this survey method has a name and I was hoping someone could tell me what it is. If you know of a survey site that can implement it, that would be very helpful. Thanks.


title: Does Indirect Least Squares satisfy OLS assumptions?

text: I thought about it and I think ILS satisfies almost all OLS assumptions but I don't know if it satisfies homoskedasticity. Please drop hints if possible, thanks!


title: Estimating probability of 2 dependent variables

text: I am working on a programming problem of 2 dependent variables X and Y. X is made of feature vector of size [128 x 1] and Y is the observation. My problem is that the 2 variables are dependent that is given X we should compute probability of Y i.e P(Y|X) or given Y we should compute P(X|Y). I am uncertain as which model can capture this property of the 2 variables. Hidden Markov model I think will only capture one way. Is there any probability distribution method to capture the distribution in both ways.


title: What does it mean to bootstrap a classifier?

text: I'm reading a paper on data mining, and the authors mentioned that "we bootstrapped a classifier". I understand that bootstrapping means sampling with replacement, but I didn't understand how that relates to classifiers.


title: Excluding one combination of two factors in a linear model

text: Imagine you are modelling an outcome based on two discrete factors, `factor1` and `factor2`, which can be yes or no. You could do something like:               model <- lm(outcome ~ factor1 + factor2)      Now suppose that out of the four combinations of values for the two factors there is one that can never happen. For instance if `factor1` is true `factor2` can be either true or false, but if `factor1` is false `factor2` can only be false. Does one need to account for that in the model? And, if yes, how so?


title: For data similar to audio, how to determine if there are 1 or 2 categories?

text: I have to compare pairs of audio strems as 1d time series. Looking at the aligned trajectories I need to either cluster them together or assume they arise from independent generators. I remember hearing of Dirichlet processes being applied. Could try to defined a hidden Markov model on the structure, not sure yet how the emission probabilities would work out, but it must be possible. any recommendations? Best EDITS- set of 1D data time series s in {1...N}, another set of trajectories q in {1...N} all over a discrete time domain t=0...t=T under a uniform interval. Task, determine whether under certain assumptions (eg. iid, normality etc), is the generator under these assumptions a valid null hypothesis to have generated the two categories of trajectories s and q? Or is the null hypothesis abandoned for a more complex model? So do we use an autoregressive model? a Gaussian process would be hard to interpret the inverse covariance. Some parametric model would seem to be valid. Don't know... any ideas?


title: What is the daily job routine of the machine learning scientist?

text: I'm a master CS student in a German university now writing my thesis. I will be done in two months I have to make the **very hard** decision if I should continue with a PhD or find a job in the industry. My reasons for doing a PhD:   * I'm a very curious person and I feel I still lack **too much** knowledge. I want to learn a lot and the PhD will help me for that, since I can do more good courses and read tons of papers and be an expert in data mining and machine learning. I love math but wasn't good at it in my undergrad (bad uni). Now in this German Uni I feel I developed a lot of great math skills and I want to improve that because I really really love math! (I was really really bad in math in my undergrad and during my lifetime but now I see I can do math well!)   * I will work with intellectually challenging stuff.   * I need to be honest and say that also I hate to see someone else with a _higher_ degree than me. So if I walk into the street and see someone with a PhD, I don't have to say "oh wow this guy is smarter than me." I prefer to be on the other side. ;) My reasons for NOT doing a PhD:   * I read on the internet about doing a PhD or not doing it. I found out that in the most and usual cases people with a PhD do the same kind of work of people with masters. (that was a general observation in computer science, not about ML/DM).   * I can start a career and make a lot of money in 1 or 2 years, then I can probably start my own company.  What is not clear yet: I still don't know what is my ultimate goal at the end. Is it to have a famous little company? Or is it to be a famous scientist? I still don't have an answer for this question yet. To help me make a decision I want to know two things:   * What is it like to work as a data scientist/machine learner with a master degree in the industry? What kind of work you do? Especially when I read those ads on Amazon as a machine learning scientist, I always wonder what they do.   * The same question as before, but with a PhD. Do you do something different or the same thing as with masters?   * Am I going to deal with challenging interesting problems? Or some boring stuff? As a slight note: I have seen a guy with a PhD in machine learning (in Germany) and is working in a company that promotes a machine learning software. As I understood most of his job is training people to use the methods and software (decision trees ..etc). It would be great if I can get some answers of experiences related to Germany/Switzerland in some famous good companies.


title: Comparison between MaxEnt, ML, Bayes and other kind of statistical inference methods

text: I'm in no way a statistician (I've had a course in mathematical statistics but nothing more than that), and recently, while studying information theory and statistical mechanics, I met this thing called "uncertainty measure"/"entropy". I read Khinchin derivation of it as a measure of uncertainty and it made sense to me. Another thing that made sense was Jaynes description of MaxEnt to get a statistic when you know the arithmetic mean of one or more function/s on the sample (assuming you accept $-\sum p_i\ln p_i$ as a measure of uncertainty of course). So I searched on the net to find the relationship with other methods of statistical inference, and God was I confused. For example this paper suggest, assuming that i got it right, that you just get a ML estimator under a suitable reformulation of the problem; MacKey, in his book, says that MaxEnt can give you weird things, and you should't use it even for a starting estimate in a Bayesian inference; etc.. I'm having trouble finding good comparisons. My question is, could you provide an explanation and/or good refences of weak and strong points of MaxEnt as a statistical inference method with quantitative comparisons to other methods (when applied to toy models for example)?


title: Linear regression effect sizes when using transformed variables

text: When performing linear regression, it is often useful to do a transformation such as log-transformation for the dependent variable to achieve better normal distribution conformation. Often it is also useful to inspect beta's from the regression to better assess the effect size/real relevance of the results. This raises the problem that when using e.g. log transformation, the effect sizes will be in log scale, and I've been told that because of non-linearity of the used scale, back-transforming these beta's will result in non- meaningful values that do not have any real world usage. This far we have usually performed linear regression with transformed variables to inspect the significance and then linear regression with the original non-transformed variables to determine the effect size. Is there a right/better way of doing this? For the most part we work with clinical data, so a real life example would be to determine how a certain exposure affects continues variables such as height, weight or some laboratory measurement, and we would like to conclude something like "exposure A had the effect of increasing weight by 2 kg".


title: How do I compare the repeatability of two sets of experiments with a low number of samples?

text: I'm investgating the repeatability of an experimental process (actually measuring lift in a wind tunnel), which can be measured using a single number, Q. However, I only have limited access to the tunnel so the number of experiments that I can run is restricted. If I have 5 experiments with the tunnel in one configuration: Q1 = [1.18, -0.41, -0.66, 0.98, 0.1] and five in another: Q2 = [-0.36 -0.73 -1.47 0.15 -0.31] How confident can I be that the repeatability of the tunnel in the second configuration is better than the first? i.e. the standard deviation of the data from the first configuration is 0.81 and from the second configuration is 0.60, but these values will clearly change if I do additional experiments. Can I somehow calculate the standard deviation of the standard deviation? Any advice on how I could calculate this would be appreciated and more generally how I can determine when I have enough experiments to say something like I am 95% confident that configuration 2 has better repeatability than configuration 1.


title: Is Area under curve a composite function

text: I have some data examples. If I split the data into three parts and the have some scores for each example of the three parts and then calculate individual AUCs for the three parts In the next case, I keep all the data examples in one cluster, use the same scores as before and calculate the AUC over this whole dataset. Is it possible that the AUC of the whole bulk could be much greater than the individual parts?


title: How to diagonalize a large sparse symmetric matrix, to get the eigen values and eigenvectors?

text: The matrix could be as large as $2500\times 2500$, what is the best algorithm to do that, is there some algorithm that is easy to write a program, is there any convenient packages for that?


title: What is the difference between a particle filter (sequential Monte Carlo) and a Kalman filter?

text: A particle filter and Kalman filter are both recursive Bayesian estimators. I often encounter Kalman filters in my field, but very rarely see the usage of a particle filter. When would one be used over the other?


title: Regression on wavy angular data?

text: I'm a newbie at stats/machine learning so please bare with me. This plot is a result of an experiment that attempts to find the perceived angle of a stimulus. ![enter image description here](http://i.stack.imgur.com/QRVeW.png) The stimulus is placed at a position (see the X axis of the plot) and participants reply with the perceived angle (see the Y axis)--in other words, the question is "what direction do you think it's coming from?" Perceived angles therefore range from 0 to 360. The plot shows this strange sinusoidal pattern about the best fit line (in blue). How would I go about capturing that wavy pattern in a regression? I would like to obtain an equation for perceived angle as a function of position. A linear regression is not a bad approximation, but I would like to see if there's anything better. I also tried to do a regression with a few added polynomial terms, but I couldn't get a reasonable fit (perhaps I just messed up). Thoughts/advice appreciated.


title: R² (squared) from a generalized linear mixed-effects models (GLMM) using a negative binomial distribution

text: I try to compute the marginal and conditional R² for a GLMM using a negative binomial distribution by following the procedure recommended by Nakagawa & Schielzeth (2013) . Unfortunately, the supplementary material of their article does not include an example of a negative binomial distribution (see the online version of the article stated below, I also added their code below). I fitted my model using the glmmPQL function from the MASS package.               full_model  <- glmmPQL ( Y~ a + b + c,  random = ~ 1 +  A | location       , family = negative.binomial (1.4 ) ,data= mydata       In particular, I do have the following problems:   1. First, I need to extract the fixed effect design matrix of my model. However, full_model @X or model.matrix(full_model) does not work. I also tried to set the argument x=TRUE before extracting the matrix. Well, this should not be too tricky, but the following problems are.    2. Second, I need to specify the distribution-speciﬁc variance of my model. Examples in the article (see table 2 & and the supplementary R code of the online article) specify this for a binomial and a Poisson distribution. With some deeper statistical knowledge, it should not be difficult to specify this for a negative binomial distribution.    3. Finally, I would need to know if glmmPQL uses additive dispersion or to multiplicative dispersion. In the paper, they state: "we only consider additive dispersion implementation of GLMMs although the formulae that we present below can be easily modiﬁed for the use with GLMMs that apply to multiplicative dispersion. " Thus, in case glmmPQL uses multiplicative dispersion, I would need further help to adjust the formula. Can anybody help? Thanks, best Philipp P.S. R-code is welcome. Nakagawa & Schielzeth (2013) A general and simple method for obtaining R 2 from generalized linear mixed-effects models. Methods in Ecology and Evolution 2013, 4, 133–142. doi: 10.1111/j.2041-210x.2012.00261.x Their R script:                 #A general and simple method for obtaining R2 from generalized linear mixed-effects models       #Shinichi Nakagawa1,2 and Holger Schielzeth3       #1 National Centre of Growth and Development, Department of Zoology, University of    Otago, Dunedin, New Zealand       #2 Department of Behavioral Ecology and Evolutionary Genetics, Max Planck Institute for Ornithology, Seewiesen, Germany       #3 Department of Evolutionary Biology, Bielefeld University, Bielefeld, Germany       #Running head: Variance explained by GLMMs       #Correspondence:       #S. Nakagawa; Department of Zoology, University of Otago, 340 Great King Street,    Dunedin, 9054, New Zealand       #Tel:  +64 (0)3 479 5046       #Fax: +64 (0)3 479 7584       #e-mail: shinichi.nakagawa@otago.ac.nz                  ####################################################       # A. Preparation       ####################################################       # Note that data generation appears below the analysis section.       # You can use the simulated data table from the supplementary files to reproduce exactly the same results as presented in the paper.            # Set the work directy that is used for rading/saving data tables       # setwd("/Users/R2")            # load R required packages       # If this is done for the first time, it might need to first download and install the package       # install.package("arm")       library(arm)       # install.package("lme4")       library(lme4)                 ####################################################       # B. Analysis       ####################################################            # 1. Analysis of body size (Gaussian mixed models)       #---------------------------------------------------            # Clear memory       rm(list = ls())            # Read body length data (Gaussian, available for both sexes)       Data <- read.csv("BeetlesBody.csv")            # Fit null model without fixed effects (but including all random effects)       m0 <- lmer(BodyL ~ 1 + (1 | Population) + (1 | Container), data = Data)            # Fit alternative model including fixed and all random effects       mF <- lmer(BodyL ~ Sex + Treatment + Condition + (1 | Population) + (1 | Container), data = Data)            # View model fits for both models       summary(m0)       summary(mF)            # Extraction of fitted value for the alternative model       # fixef() extracts coefficents for fixed effects       # mF@X returns fixed effect design matrix       Fixed <- fixef(mF)[2] * mF@X[, 2] + fixef(mF)[3] * mF@X[, 3] + fixef(mF)[4] * mF@X[, 4]            # Calculation of the variance in fitted values       VarF <- var(Fixed)            # An alternative way for getting the same result       VarF <- var(as.vector(fixef(mF) %*% t(mF@X)))            # R2GLMM(m) - marginal R2GLMM       # Equ. 26, 29 and 30       # VarCorr() extracts variance components       # attr(VarCorr(lmer.model),'sc')^2 extracts the residual variance       VarF/(VarF + VarCorr(mF)$Container[1] + VarCorr(mF)$Population[1] + attr(VarCorr(mF), "sc")^2)            # R2GLMM(c) - conditional R2GLMM for full model       # Equ. XXX, XXX       (VarF + VarCorr(mF)$Container[1] + VarCorr(mF)$Population[1])/(VarF +    VarCorr(mF)$Container[1] + VarCorr(mF)$Population[1] + (attr(VarCorr(mF), "sc")^2))            # AIC and BIC needs to be calcualted with ML not REML in body size models       m0ML <- lmer(BodyL ~ 1 + (1 | Population) + (1 | Container), data = Data, REML = FALSE)       mFML <- lmer(BodyL ~ Sex + Treatment + Condition + (1 | Population) + (1 | Container), data = Data, REML = FALSE)            # View model fits for both models fitted by ML       summary(m0ML)       summary(mFML)                 # 2. Analysis of colour morphs (Binomial mixed models)       #---------------------------------------------------            # Clear memory       rm(list = ls())       # Read colour morph data (Binary, available for males only)       Data <- read.csv("BeetlesMale.csv")            # Fit null model without fixed effects (but including all random effects)       m0 <- lmer(Colour ~ 1 + (1 | Population) + (1 | Container), family = "binomial", data = Data)            # Fit alternative model including fixed and all random effects       mF <- lmer(Colour ~ Treatment + Condition + (1 | Population) + (1 | Container), family = "binomial", data = Data)            # View model fits for both models       summary(m0)       summary(mF)            # Extraction of fitted value for the alternative model        # fixef() extracts coefficents for fixed effects        # mF@X returns fixed effect design matrix       Fixed <- fixef(mF)[2] * mF@X[, 2] + fixef(mF)[3] * mF@X[, 3]            # Calculation of the variance in fitted values       VarF <- var(Fixed)            # An alternative way for getting the same result       VarF <- var(as.vector(fixef(mF) %*% t(mF@X)))            # R2GLMM(m) - marginal R2GLMM       # see Equ. 29 and 30 and Table 2       VarF/(VarF + VarCorr(mF)$Container[1] + VarCorr(mF)$Population[1] + pi^2/3)            # R2GLMM(c) - conditional R2GLMM for full model       # Equ. XXX, XXX       (VarF + VarCorr(mF)$Container[1] + VarCorr(mF)$Population[1])/(VarF +     VarCorr(mF)$Container[1] + VarCorr(mF)$Population[1] + pi^2/3)                 # 3. Analysis of fecundity (Poisson mixed models)       #---------------------------------------------------            # Clear memory       rm(list = ls())            # Read fecundity data (Poisson, available for females only)       Data <- read.csv("BeetlesFemale.csv")            # Creating a dummy variable that allows estimating additive dispersion in lmer        # This triggers a warning message when fitting the model       Unit <- factor(1:length(Data$Egg))            # Fit null model without fixed effects (but including all random effects)       m0 <- lmer(Egg ~ 1 + (1 | Population) + (1 | Container) + (1 | Unit), family = "poisson", data = Data)            # Fit alternative model including fixed and all random effects       mF <- lmer(Egg ~ Treatment + Condition + (1 | Population) + (1 | Container) + (1 | Unit), family = "poisson", data = Data)            # View model fits for both models       summary(m0)       summary(mF)            # Extraction of fitted value for the alternative model        # fixef() extracts coefficents for fixed effects        # mF@X returns fixed effect design matrix       Fixed <- fixef(mF)[2] * mF@X[, 2] + fixef(mF)[3] * mF@X[, 3]            # Calculation of the variance in fitted values       VarF <- var(Fixed)            # An alternative way for getting the same result       VarF <- var(as.vector(fixef(mF) %*% t(mF@X)))            # R2GLMM(m) - marginal R2GLMM        # see Equ. 29 and 30 and Table 2        # fixef(m0) returns the estimate for the intercept of null model       VarF/(VarF + VarCorr(mF)$Container[1] + VarCorr(mF)$Population[1] + VarCorr(mF)$Unit[1] + log(1 + 1/exp(as.numeric(fixef(m0)))))            # R2GLMM(c) - conditional R2GLMM for full model       # Equ. XXX, XXX       (VarF + VarCorr(mF)$Container[1] + VarCorr(mF)$Population[1])/(VarF +    VarCorr(mF)$Container[1] + VarCorr(mF)$Population[1] + VarCorr(mF)$Unit[1] + log(1 +                                                                             1/exp(as.numeric(fixef(m0)))))                 ####################################################       # C. Data generation       ####################################################            # 1. Design matrices        #---------------------------------------------------            # Clear memory       rm(list = ls())            # 12 different populations n = 960       Population <- gl(12, 80, 960)            # 120 containers (8 individuals in each container)       Container <- gl(120, 8, 960)            # Sex of the individuals. Uni-sex within each container (individuals are sorted at the pupa stage)       Sex <- factor(rep(rep(c("Female", "Male"), each = 8), 60))            # Condition at the collection site: dry or wet soil (four indiviudal from each condition in each container)       Condition <- factor(rep(rep(c("dry", "wet"), each = 4), 120))            # Food treatment at the larval stage: special food ('Exp') or standard food ('Cont')       Treatment <- factor(rep(c("Cont", "Exp"), 480))            # Data combined in a dataframe       Data <- data.frame(Population = Population, Container = Container, Sex = Sex, Condition = Condition, Treatment = Treatment)                 # 2. Gaussian response: body length (both sexes)       #---------------------------------------------------            # simulation of the underlying random effects (Population and Container with variance of 1.3 and 0.3, respectively)       PopulationE <- rnorm(12, 0, sqrt(1.3))       ContainerE <- rnorm(120, 0, sqrt(0.3))            # data generation based on fixed effects, random effects and random residuals errors       Data$BodyL <- 15 - 3 * (as.numeric(Sex) - 1) + 0.4 * (as.numeric(Treatment) - 1) + 0.15 * (as.numeric(Condition) - 1) + PopulationE[Population] + ContainerE[Container] +          rnorm(960, 0, sqrt(1.2))            # save data (to current work directory)       write.csv(Data, file = "BeetlesBody.csv", row.names = F)                 # 3. Binomial response: colour morph (males only)       #---------------------------------------------------            # Subset the design matrix (only males express colour morphs)       DataM <- subset(Data, Sex == "Male")            # simulation of the underlying random effects (Population and Container with variance of 1.2 and 0.2, respectively)       PopulationE <- rnorm(12, 0, sqrt(1.2))       ContainerE <- rnorm(120, 0, sqrt(0.2))            # generation of response values on link scale (!) based on fixed effects and random effects       ColourLink <- with(DataM, 0.8 * (-1) + 0.8 * (as.numeric(Treatment) - 1) + 0.5 *    (as.numeric(Condition) - 1) + PopulationE[Population] + ContainerE[Container])            # data generation (on data scale!) based on negative binomial distribution       DataM$Colour <- rbinom(length(ColourLink), 1, invlogit(ColourLink))            # save data (to current work directory)       write.csv(DataM, file = "BeetlesMale.csv", row.names = F)                 # 4. Poisson response: fecundity (females only)       #---------------------------------------------------            # Subset the design matrix (only females express colour morphs)       DataF <- Data[Data$Sex == "Female", ]            # random effects       PopulationE <- rnorm(12, 0, sqrt(0.4))       ContainerE <- rnorm(120, 0, sqrt(0.05))            # generation of response values on link scale (!) based on fixed effects, random effects and residual errors       EggLink <- with(DataF, 1.1 + 0.5 * (as.numeric(Treatment) - 1) + 0.1 *   (as.numeric(Condition) - 1) + PopulationE[Population] + ContainerE[Container] +   rnorm(480,                                                                                0, sqrt(0.1)))            # data generation (on data scale!) based on Poisson distribution       DataF$Egg <- rpois(length(EggLink), exp(EggLink))            # save data (to current work directory)       write.csv(DataF, file = "BeetlesFemale.csv", row.names = F)


title: Regularization and Mean Estimation

text: Suppose I have some i.i.d. data $x_1, \ldots, x_n \sim N(\mu, \sigma^2)$, where $\sigma^2$ is fixed and $\mu$ is unknown, and I want to estimate $\mu$. Instead of simply giving the MLE of $\mu = \bar{x}$, one could estimate (1) $\mu = \lambda \mu_0 + (1 - \lambda) \bar{x},$ for some "prior best guess" $\mu_0$. This also has a nice Bayesian interpretation: we place a prior $\mu \sim N(\mu_0, \sigma^2_0)$ on $\mu$, and $\lambda$ is the weighted precision. I seem to recall that this also has an explicit L2 regularization- interpretation (i.e., we choose some penalty and minimize the squared loss to get the above estimate), similar to things like lasso and ridge regression, but I can't remember how it goes. Can anyone explain what the L2 regularization-interpretation of (1) is? [More general answers, where the data isn't necessarily normal distributed, are welcome as well.]


title: Generate uniformly distributed weights that sum to unity?

text: It is common to use weights in applications like mixture modeling and to linearly combine basis functions. Weights $w_i$ must often obey $w_i ≥$ 0 and $\sum_{i} w_i=1$. I'd like to randomly choose a weight vector $\mathbf{w} = (w_1, w_2, …)$ from a uniform distribution of such vectors. It may be tempting to use $w_i = \frac{\omega_i}{\sum_{j} \omega_j}$ where $\omega_i \sim$ U(0, 1), however as discussed in the comments below, the distribution of $\mathbf{w}$ is not uniform. However, given the constraint $\sum_{i} w_i=1$, it seems that the underlying dimensionality of the problem is $n-1$, and that it should be possible to choose a $\mathbf{w}$ by choosing $n-1$ parameters according to some distribution and then computing the corresponding $\mathbf{w}$ from those parameters (because once $n-1$ of the weights are specified, the remaining weight is fully determined). The problem appears to be similar to the sphere point picking problem (but, rather than picking 3-vectors whose $ℓ_2$ norm is unity, I want to pick $n$-vectors whose $ℓ_1$ norm is unity). Thanks!


title: Probability of Game Outcomes (with ties!)

text: Let's say I'm playing 10 "games". For each game, I know the probability of winning, the probability of tying, and the probability of losing. From these values, I can calculate the probability of winning X games, the probability of losing X games, and the probability of tying X games (for X = 0 to 10). I'm just trying to figure out the probability of each outcome after playing all 10 games. For example, what is the probability of winning 7, losing 2, and tying 1? Any ideas? Thanks! * * * Here's some example data for 2 games: Game 1:   * win: 23.3%   * tie: 1.1%   * lose: 75.6% Game 2:   * win: 29.5%   * tie: 3.2%   * lose: 67.3% Based on this, we can calculate the probability after playing 2 games of: * * *   * 0 wins: 54.0%   * 1 win: 39.1%   * 2 wins: 6.9% * * *   * 0 ties: 95.8%   * 1 tie: 4.2%   * 2 ties: 0.0% * * *   * 0 losses: 8.0%   * 1 loss: 41.1%   * 2 losses: 50.9% * * * Based on these numbers, what is the probability of each outcome? The possible outcomes (W-L-T) would be:   * 2-0-0   * 1-1-0   * 1-0-1   * 0-1-1   * 0-2-0   * 0-0-2 Obviously this is pretty trivial with 2 games, but is there a generic formula for finding the probability of _W_ wins, _T_ ties, and _L_ losses after X games?


title: Clustering genes in a time course experiment

text: I have seen a few queries on clustering in time series and specifically on clustering, but I don't think they answer my question. **Background:** I want to cluster genes in a time course experiment in yeast. There are four time points say: _t1_ _t2_ _t3_ and _t4_ and total number of genes _G_. I have the data in form a matrix _M_ in which the columns represent the treatments (or time points) _t1_ _t2_ _t3_ and _t4_ and the rows represent the genes. Therefore, _M_ is a Gx4 matrix. **Problem:** I want to cluster the genes which behave the same across all time points _t1_ _t2_ _t3_ and _t4_ as well as within a particular time point _ti_ , where i is in {1, 2, 3, 4} (In case we cannot do both the clusterings together, the clustering within a time point is more important than clustering across time points). In addition to this, I also want to draw a heatmap. **My Solution:** I use the R code below to obtain a heatmap as well as the clusters using `hclust` function in R (performs hierarchical clustering with euclidean distance)                   row.scaled.expr <- (expr.diff - rowMeans(expr.diff)) / rowSds(expr.diff)              breaks.expr <- c(quantile(row.scaled.expr[row.scaled.expr < 0],                                    seq(0,1,length=10)[-9]), 0,                                    quantile(row.scaled.expr[row.scaled.expr > 0],                                    seq(0,1,length=10))[-1] )                   blue.red.expr <- maPalette(low = "blue", high = "red", mid = "white",                          k=length(breaks.expr) - 1)              pdf("images/clust.pdf",              height=30,width=20,pointsize=20)         ht1 <- heatmap.2(row.scaled.expr, col = blue.red.expr, Colv = FALSE, key = FALSE,            dendrogram = "row", scale = "none", trace = "none",           cex=1.5, cexRow=1, cexCol=2,           density.info = "none", breaks = breaks.expr,            labCol = colnames(row.scaled.expr),           labRow="",           lmat=rbind( c(0, 3), c(2,1), c(0,4) ), lhei=c(0.25, 4, 0.25 ),           main=expression("Heat Map"),           ylab="Genes in the Microarray",           xlab="Treatments"           )         dev.off()      I recently discovered `hopach` package in _Bioconductor_ which can be used to estimate the number of clusters. Previously, I was randomly assigning the number of bins for the heatmap and cutting the tree at an appropriate height to get a pre-specified number of clusters. **Possible Problems in my solution:**   1. I may be not clustering the genes within a particular treatment and clustering genes only across treatments or vice versa.   2. There may be better ways of obtaining a heatmap for the pattern I want to see (similar genes within a treatment and across treatments).   3. There may be better visualization methods which I am not aware of. **Note:**   1. _csgillespie_ (moderator) has a more general document on his website in which he discusses all the aspects of time course analysis (including heatmaps and clustering). I would appreciate if you can point me to an articles which describe heatmaps and clustering in detail.   2. I have tried the `pvclust` package, but it complains that _M_ is singular and then it crashes.


title: Calculating and interpreting effect sizes for interaction terms

text: Can an effect size be calculated for an interaction effect in general and more specifically using the F-statistic and its associated degrees of freedom? If yes, should this be done and what is the appropriate interpretation of the effect size given the nature of interaction effects? Initially, I assumed the answer is "no." However, a quick Google search turned up statements about computing effect sizes for interaction terms. Any assistance in clarifying this issue will be greatly appreciated!


title: Clustering with some cluster centers fixed/known

text: Thanks for reading my question. I have several thousand data points scattered on an (x,y) grid that I am trying to cluster. The data points are not uniformly distributed across the grid, but are concentrated in certain areas. I am most interested in identifying the centers of the clusters as representing starting points that minimize the average (Euclidean) distance from a point to the nearest cluster center. Depending on the specific model and data set, there are between 3 and 7 clusters. The number of clusters is known beforehand in each instance, and does not need to be determined by an algorithm. In each situation, some of the centers of the clusters are known (0 to 4 known starting points), but the rest are unknown. The goal is to identify the centers of the unknown clusters that minimize the average distance to a center across all the data points. How can I run a clustering algorithm where I can specify a certain number of cluster centers, and solve for the others? I am using R, and have looked primarily at package mclust. My thought was that specifying priors for mean and scale to Mclust with very small scale variances for the known centers and very large scale variances (uninformative prior) for the unknowns would be a good approach, but I am having trouble coding it. The available examples for specifying priors in the package documentation aren't terribly helpful (to me), and might be used for a completely different purpose than what I'm trying to do. My attempt to code this in R looks something like:               # create data matrix of points     x <- rnorm(100, 50, 25)     y <- rnorm(100, 50, 25)     my.data <- cbind(x,y)          # Two known centers, rest are unknown so provide mean of x,y as default starting point     #     # known centers are (25, 10) and (90, 65), assume midpoints of grid for others     x.prior.mean <- c(25, 90, 50, 50, 50)     y.prior.mean <- c(10, 65, 50, 50, 50)          # Provide small scale (variance) for known centers, large scale for unknown centers (uninformative prior)     x.prior.scale <- c( 0.1, 0.1, 100, 100, 100)     y.prior.scale <- c( 0.1, 0.1, 100, 100, 100)          # Create a cluster model with no prior specified     my.clust.noprior <- Mclust(data=my.data, G=5)          # Now add a prior for mean     my.clust.prior <-   Mclust(data=my.data, G=5, prior=priorControl(mean=cbind(x=x.prior.mean, y=y.prior.mean)))          # Compare what I think are the centers of the clusters (mean of parameters).     # The centers in the prior-specified case don't seem to reflect the known centers     my.clust.noprior$parameters$mean     my.clust.prior$parameters$mean          # Commented out, but attempting the following statement that adds scale parameter yields an error:     #    Error in chol.default(priorParams$scale) : non-square matrix in 'chol'     #     # my.clust <-   Mclust(data=my.data, G=5, prior=priorControl(mean=cbind(x=x.prior.mean, y=y.prior.mean), scale=cbind(x.prior.scale, y.prior.scale)))      Is there a way to accomplish what I'm trying to do? I am open to using other R packages besides mclust if there's one better suited for this problem. Thank you


title: What can you do with 'crazy' data?

text: This question is more about an approach to a complicated data situation rather than particular statistical methods. I'm modeling our organization's electricity bills, and I have monthly billing data from 2008 to the present. I have a fairly accurate model for electricity usage, based on the month's average temperature. And I have a fairly good model for the basics of cost, based on usage, demand, and peak-season surcharges. (The cost model is what I really need, for budgeting purposes.) BUT the problem is the data. To give an example, in 2010 we received two different sets of credits, based on two different sets of overcharges in 2009, while at the same time we were evidently experiencing the hottest summer and coldest winter in a century, which means that for the next two years (2011-2013) we will be paying a surcharge to make up for that. It's a messy situation and I'm not sure how to handle it. My first attempt was to modify the data by shifting the refunds to the appropriate period in the past. Unfortunately, one of the refunds is calculated based on data I can't get, so it's just a guess on my part. And who knows how much the two-year surcharge will actually be. So I'm wondering what the proper approach would be, if any?   1. Is it even worth trying to adjust the data? Should I just accept that the bills are the bills and not care about how money is shifting around from year to year? I would think it would result in much higher variance and would mess up trends.   2. Could I try to use indicator variables to indicate when various surcharges and refunds hit? I can't believe this would work with only 3-4 years of data.   3. Should I try to model the various parts of the bill separately? Some (like the usage and demand charges) might be fairly stable, while others (like refunds or fuel charges) would be appropriately volatile.   4. Is it a mistake to try to model monthly bills at all? Should I model at the year level to hopefully smooth out things a bit? I'd appreciate any ideas or suggestions.


title: Whether to transform non-normal pre-test when running linear regression on transformed post-test?

text: I'm running linear regression model on a post-intervention test score controlling for pre-intervention test score. I used Box-Cox transformation on the post-intervention test score to normalize it. Since there no normality assumption imposed on independent variable, I plan not to transform the pre- intervention test score; but since the pre- and post-test are on a same scale, should I transform the pre-test as well just to be consistent?


title: What is the simplest way to classify airplan manuvers?

text: Suppose we have declared four motion types for air-plane. If we represent each maneuver with a trajectory line, what is the best classification method to retrieve the trajectory pattern with a similar trajectory shape to the new trajectory path? The proposed method should be able to classify partial paths so that we can predict maneuver type before it becomes complete. Thank for your help.


title: Creating classification features from wavelet transformed time series

text: I'm interested in using a wavelet transform, Haar for example, to create classification variables from time series data to use in logistic regression. Simple example. Let's say I'm trying to predict payment defaults and I have a person's monthly expense data and I know that someone with consistent expenses is better than someone with increasing expenses in the the most recent 4 months. If I have two sample borrowers: Borrower A - Good - expensesA = c(100,110,95,105), default = 0 Borrower B - Bad - expensesB = c(75,100,150,200), default = 1 If I am using logistic regression to create a classification model, using the R `wavelets` package `dwt()` function for a "haar" transform of the time series I get the features below. `W`s are wavelet coefficients and the `V`s the scaling coefficients. Do I need to use all four `W1` and `V1` values as variables to properly model this or is it okay to try just the `W1`s without `V1`s (or vice versa)? Is it worthwhile to try only the single `W2` and `V2`s as variables? Or is it better to try to use a clustering algorithm and label them based on clusters created from these coefficients? I know it of course also depends on the data, but I'm looking for a starting point regarding best practices.               tr = dwt(expensesA, filter = "haar")     tr          An object of class "dwt"     Slot "W":     $W1     [,1]     [1,] 7.071068     [2,] 7.071068          $W2     [,1]     [1,]   -5               Slot "V":     $V1     [,1]     [1,] 148.4924     [2,] 141.4214          $V2     [,1]     [1,]  205               Slot "filter":     Filter Class: Daubechies     Name: HAAR     Length: 2     Level: 1     Wavelet Coefficients: 7.0711e-01 -7.0711e-01     Scaling Coefficients: 7.0711e-01 7.0711e


title: Knnclassify() matlab

text: I am currently working on matlab and I am new to it. I would like to know how do i average the results from the folds (or otherwise combined) to produce a single estimation. My dataset has 2 class labels +1 and -1. and k=5. I use the knnclassify() method in matlab in order to perform cross validation. Kindly help!


title: Arima residual calculation and comparison with R

text: I have simulated an ARIMA(1,0,1) process using R. Below is the code.               library(stats)     library(forecast)     set.seed(101)     y <- arima.sim(n=1000,list(ar=0.2,ma=0.6))     plot(y)     qa <- arima(y,order=c(1,0,1))      Below is the output:               ARIMA(1,0,1) with non-zero mean           Coefficients:              ar1     ma1  intercept           0.2022  0.5895    -0.0713     s.e.  0.0430  0.0349     0.0605          sigma^2 estimated as 0.9236:  log likelihood=-1379.57     AIC=2767.14   AICc=2767.18   BIC=2786.77      I'm trying to see if the estimated residuals matches the calculated residuals for example, an ARMA (1,0,1) residual can be written as: $$ e_{t} = y_{t} - \mu -\phi y_{t-1}+\theta e_{t-1} $$ $$ \mu = intercept $$ $$ \phi=ar1 - coefficient $$ $$ \theta = ma1 - coefficient $$ Code for calculating residuals                     n = length(y)           e = rep(1,n)           e[1] = 0 ##since there is no residual at 1, e1 = 0           for (t in (2 : n)){             e[t] = y[t] +0.0713 - 0.2022 *y[t-1] + 0.5895*e[t-1]           }      Inputting all the numbers into equation say for residual 10 $$ e_{10} = y_{10} - (-0.0713) -0.2022y_{9} +0.5895e_{9} $$ calculated residual is:               > e[10]     [1] -0.8278073      residual from running arima in R is               > qa$residuals[10]     [1] -0.1424167      The estimated and calculated numbers don't match. My questions are   1. Am I doing something incorrect here ?   2. Why is there a difference between calculated and estimated numbers?


title: How to determine what method is being used?

text: I was wondering what time series model would forecast future values in such a static but decreasing manner? ( I'm talking about the values around where the circle appears, as apparently, the data set that I have does that, judging from the graph, but I'm a beginner and not sure what model would do that? though I think it should be a simple one). I just need reasonable ideas :) Thanks! ![data](http://i.stack.imgur.com/vWStD.png)


title: Non-parametric z-score for comparison of outliers/extreme values across set of non-normal distributions

text: I have a network and for each node in the network I am calculating a certain property, essentially its 'connectedness' to a set of source nodes. In order to identify nodes which are better connected than nodes _at a similar geometric distance from the source nodes_ , I have grouped nodes into shells of given radii from the source nodes. So for instance, all nodes which lie within 0-2,2-4,4-6,... units of the source nodes. The following figure shows the distributions of 'connectedness' within the shells. ![enter image description here](http://i.stack.imgur.com/J71PY.png) I think the best way to describe most of the distributions is that they look roughly straight on a log-log plot (as you can see in the boxed example) i.e. they have wide/long-tails. Clearly, they are not Gaussian. Ideally, I would like to calculate a single score for each node, which gives a measure of how 'unexpectedly' well connected it is and is comparable across shells. I started by using a z-score $\frac{x_i - \bar x}{\sigma}$, but given the distributions are far from normal this doesn't seem like the best measure. I have switched to a z-score based on the median $\frac{x_i - \text{median}(x)}{\sigma}$, but this still feels unsatisfactory. Is there some kind of non-parametric measure which could be used to compare the significance of the outliers across shells and obtain a single ranking for all nodes in the network?


title: Visualizing overlapping categories?

text: I currently have a sample of 100+ job descriptions, which I am categorizing with the intention of identifying the most common areas of job activity within a company - including, and this is important - overlaps. What I am looking for is a software solution that would help me display the ensuing relationships visually. Ideally, I would like to have a visual "web" of sorts (think of a universe map in a 4X space strategy game) that can be rotated, zoomed in or out, and that highlights stronger overlaps. In the end, I would be able to display visually that most positions that do project manager also do financial control (for example), but very few project management positions do training delivery. Each job description is categorized by its responsibilities - e.g., project management, data research, recruitment, etc. It is important to emphasize that each document can be long to two or more categories (on average, 6-8) - e.g., the same position can be responsible for project management, recruitment and selection, fundraising, new project development, process analysis, etc. I am currently using SPSS Text Analytics for the categorization effort, which means that each document has a unique ID in the system and can be exported as such, along with the free text of the job description and the categories it has been assigned to. I currently have about 15 such categories in total, but I anticipate this number growing rapidly as the sample is expanded. Is there such software on the market, or can you think of a more effective way of presenting this information?


title: How to interpret augmented Dickey-Fuller unit root test in R

text: I am looking at a time series which has no obvious trend, but seems to have an intercept a little above zero. The results I get for `ur.df` function in R is the following:               logprice_df <- ur.df(test3, lags = 1, type= 'trend')     summary(logprice_df)          ###############################################      # Augmented Dickey-Fuller Test Unit Root Test #      ###############################################           Test regression trend                Call:     lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)          Residuals:          Min       1Q   Median       3Q      Max      -0.50614 -0.04394  0.00134  0.03859  0.64408           Coefficients:                   Estimate Std. Error t value Pr(>|t|)         (Intercept)  1.841e-02  8.268e-03   2.226  0.02626 *       z.lag.1     -1.573e-02  5.635e-03  -2.791  0.00537 **      tt           9.234e-06  1.080e-05   0.855  0.39272         z.diff.lag   1.411e-01  3.364e-02   4.195 3.01e-05 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 0.07512 on 865 degrees of freedom     Multiple R-squared: 0.02651,    Adjusted R-squared: 0.02314      F-statistic: 7.852 on 3 and 865 DF,  p-value: 3.572e-05                Value of test-statistic is: -2.791 2.6012 3.8997           Critical values for test statistics:            1pct  5pct 10pct     tau3 -3.96 -3.41 -3.12     phi2  6.09  4.68  4.03     phi3  8.27  6.25  5.34      The problem is that I do not understand how one shall interpret these values. What is `F-statistic`? And what is the probabilities below `Pr(>|t|)` in the first table, and also `Value of test-statistic`? I really appreciate any help.


title: What is the difference between online and batch Learning?

text: I currently read the paper Efficient Online and Batch Learning using Forward- Backward Splitting by John Duchi and Yoram Singer. I am very confused about the usage of the terms 'Online' and 'Batch'. I thought 'Online' means we update the weight parameters after processing one unit of the training data. Then we use the new weight parameters to process the next unit of the training data. However, in the paper above, the usage is not that clear.


title: Custom Testings Results Dashboard

text: I am creating a testing dashboard for displaying the results of our website testing. We run A/B tests and Multivariate tests. My testing code is sound (the gathering of the raw conversion data), but now I am working on finishing out our results reporting admin dashboard. (before we used a very old version of JMP to analyze our results but that was getting tedious so we/I decided to code up an admin dashboard with R) Here is an example set of our data. (for an A/B test that we ran)               >a=c(1.3571428571429,1.3941176470588,1.275,1.2408759124088,1.5842696629213,          1.5537634408602,1.4590909090909,1.4018691588785,1.495145631068,1.1794871794872,          1.5114503816794,1.2569832402235,1.5336538461538)     >b=c(1.4891304347826,1.4656084656085,1.44,1.5190839694656,1.6084656084656,          1.4166666666667,1.5665024630542,1.2533333333333,1.4081632653061,1.6173913043478,          1.468253968254,1.5392670157068,1.5081081081081)      And with that data I am running the following R code to calculate the p-value / statistical significance.               >t.test(a, b, paired=TRUE)      Now after a week or two we would like to know how long we would need to continue to run the test before we would get to a given level of statistical significance, given the current data that we have collected. I tried `n.ttest` but I have since been told that this is more meant for power analysis. So my question is what R function should I use to get this information? Also on a related note I would be interested to know if there is any other calculations I should be running on our test data that we need that I could also add to our testing dashboard.


title: Average over two variables: Why do standard error of mean and error propagation differ and what does that mean?

text: I'm doing an experiment with a cryostat to determine the critical temperature for lead. To avoid asymmetries, I determine the critical temperature both through heating (going from 2 K to 10 K) and cooling (10 K -> 2 K). Now I have two values, that differ slighty and I average them. So a measurement of (6.942 $\pm$ 0.020) K and (6.959 $\pm$ 0.019) K gives me an average of 6.951 K. **Now the question is: what is the error of that average?** One way to do it would be to calculate the variance of this sample (containing two points), take the square root and divide by $\sqrt{2}$. This gives me an SEM of 0.0085 K, which is too low for my opinion (where does this precision come from?) The other way is to say the the mean is a function of two variables, $\bar{T} = \frac{T_1 + T_2}{2}$, therefore by error propagation the error is $\Delta T = \frac12\sqrt{(\Delta T_1)^2+(\Delta T_2)^2}$, and that gives me a much more rational value of 0.014. I see how those values differ in terms of numbers, but **which one is correct when talking about the correct estimate for the standard deviation?**


title: What is a "good fit" Brier score and Harrell's C Index

text: This is a question I originally posted on r-help but it is more suited here. I will post the question and the answer I received from Dr. Winsemius and would be most grateful for any additional answers you can provide. I am evaluating survival models using Brier score (“peperr”) and Harrell’s C-Index (“Hmisc”). I am wondering:   1. What would be considered a “good fit” according to these scores (like the heuristic levels we have for R square in linear regressions) ?   2. Are there any papers to cite on the matter (I couldn’t find any) ?   3. Is there any paper to cite that discusses the limitation of using traditional reporting for model fit in survival analysis as opposed to these measures ? Dr. David Winsemius replied:   1. Frank Harrell's excellent text "Regression Modeling Strategies" has an extensive discussion of "goodness of fit" and the principles of model comparison. It's both too involved as well as off-topic for Rhelp. The other text to consult is Steyerberg's "Clinical Prediction Models".   2. I predict that the RMS bibliography would be an excellent place to start your search. Despite getting his name attached to what he calls the 'c-index', I don't think one could call Frank Harrell a proponent of that measure or any of the "competitors". It's really just a dressed up/transformed AUC. The message I have taken from reading his book and listening to presentations is that one should apply biologic tests of sensibility as well as careful investigation of the functional relationships between candidate predictors and the outcomes of interest. He speaks very disparagingly about automatic procedures.


title: Specify a Zero-inflated (Hurdle) Gamma Model in JAGS/BUGS

text: I'm trying to use a zero-inflated gamma model (or a gamma 'hurdle' model). The model is a mixture of logistic regression and generalized linear modeling. I can do this analysis in two steps: 1) do a logistic regression against presence/absence data and then 2) Use a generalized linear model with a gamma distribution on the positive values. I'm trying to set the model up where, if y = 0, then E(y) = p, but if y > 0, then E(y) is gamma distributed. I'm trying to set this up in BUGS/JAGS, because I've seen these models worked before for poisson-distributions. I tried to adapt the code from that link into a gamma distribution, but I can't quite seem to get it to work. I know the code won't work because my data have zeroes, and the likelihood function can't be evaluated with zeroes. However, even when restricted to positive values, I get errors for invalid parent values. Even so, I'm not even sure the model is specified correctly. Here is the model:               # For the ones trick     C <- 10000          # for every observation     for(i in 1:N){         # log-likelihood of the observation from the gamma likelihood         LogPos[i] <- (shape[i] - 1)*log(y[i]) - y[i]/scale[i] - N*shape[i]*log(scale[i]) - N*loggam(scale[i])         #likelihood         Lpos[i] <- exp(LogPos[i])              # redefine the shape and rate parameters as a function of the mean and sd         shape[i] <- pow(mu[i], 2) / pow(sd, 2)         scale[i] <- mu[i] / pow(sd, 2)              # mu is a function of MTD: use the inverse link         #mu[i] <- 1/eta[i]         mu[i] <- beta0 + beta1*MTD[i]                   # zero-inflated part, where w[i] is the probability of being zero         w[i] <- exp(zeta[i]) / (1 + exp(zeta[i]))         zeta[i] <- gamma0 + gamma1*MPD[i]              # ones trick         p[i] <- Lpos[i] / C         ones[i] ~ dbern(p[i])              # Full likelihood         Lik[i] <- Lpos[i]*(1 - w[i]) + equals(y[i], 0)*w[i]       }           # PRIORS     beta0 ~ dnorm(0, 0.001)     beta1 ~ dnorm(0, 0.001)          gamma0 ~ dnorm(0, 0.001)     gamma1 ~ dnorm(0, 0.001)          sd ~ dunif(0, 100)      Has anyone set a model up like this or have any advice on how to set it up correctly? **UPDATE** I've tried a new set of code that's similar, but slight different. I still have not gotten it to work               model{            # For the ones trick       C <- 10000            # for every observation       for(i in 1:N){              # make a dummy variable that is 0 if y is < 0.0001 and 1 if y > 0.0001. This is essentially a presence         # absence dummy variable         z[i] <- step(y[i] - 0.0001)              # define the logistic regression model, where w is the probability of occurance.         # use the logistic transformation exp(z)/(1 + exp(z)), where z is a linear function         w[i] <- exp(zeta[i]) / (1 + exp(zeta[i]))         zeta[i] <- gamma0 + gamma1*MPD[i]              # define the gamma regression model for the mean. use the log link the ensure positive, non-zero mu         mu[i] <- exp(eta[i])         eta[i] <- beta0 + beta1*MTD[i]              # redefine the mu and sd of the continuous part into the shape and scale parameters         shape[i] <- pow(mu[i], 2) / pow(sd, 2)         scale[i] <- mu[i] / pow(sd, 2)              # for readability, define the log-likelihood of the gamma here         logGamma[i] <- (shape[i] - 1)*log(y[i]) - y[i]/scale[i] - shape[i]*log(scale[i]) - loggam(scale[i])              # define the total likelihood, where the likelihood is (1 - w) if y < 0.0001 (z = 0) or         # the likelihood is w * gammalik if y >= 0.0001 (z = 1). So if z = 1, then the first bit must be         # 0 and the second bit 1. Use 1 - z, which is 0 if y > 0.0001 and 1 if y < 0.0001         logLik[i] <- (1 - z[i]) * log(1 - w[i]) + z[i] * ( log(w[i]) + logGamma[i] )              # Use the ones trick         p[i] <- logLik[i] / C         ones[i] ~ dbern(p[i])       }             # PRIORS       beta0 ~ dnorm(0, 0.001)       beta1 ~ dnorm(0, 0.001)            gamma0 ~ dnorm(0, 0.001)       gamma1 ~ dnorm(0, 0.001)            sd ~ dgamma(1, 2)          }      **UPDATE 2** I've gotten it to run by deleting the definition of logGamma[i] and putting it directly into the likelihood function, which now reads:               logLik[i] <- (1 - z[i]) * log(1 - w[i]) + z[i] * ( log(w[i]) + (shape[i] - 1)*log(y[i]) - y[i]/scale[i] - shape[i]*log(scale[i]) - loggam(scale[i]) )      The problem was that JAGS was trying to evaluate the log likelihood at all observations first, resulting in NA's for the 0's. In the new way, it only evaluates it if z = 1 (I think). I'm still having trouble getting the parameter estimates to line up. For example, the gamma's are almost identical to a logistic regression of the same form (hooray!). But the betas are pretty far off from a gamma GLM of the positive values. I don't know if this is normal or not, but I suspect there are still problems with my model specification.


title: Discovering subgroups based on unsupervised analysis of one variable

text: I am looking in a dataset that consists of frequencies of mutations of a certain type and want to identify distinct populations within that dataset so I can then use supervised analysis to pick out what is different between those groups from other data I've got. What approaches would you recommend? (Hierarchical clustering doesn't work because at least two variables are needed).


title: Best way to visually present relationships from a multiple linear model

text: I have a linear model with about 6 predictors and I'm going to be presenting the estimates, F values, p values, etc. However, I was wondering what would be the best visual plot to represent the individual effect of a single predictor on the response variable? Scatterplot? Conditional Plot? Effects plot? etc? How would I interpret that plot? I'll be doing this in R so feel free to provide examples if you can. EDIT: I'm primarily concerned with presenting the relationship between any given predictor and the response variable.


title: Given two dependent sequences, predict next output for one of them

text: I have a problem which can be generalized as follows. There is some black-box which produces actions {0, 1}. Having two sequences: **_X_** - series of previous actions, **_T_** - series of time intervals between actions, how can I predict next action? **_T_** in this problem is just additional data, I am not interested in predicting time of next action. Which methods can be used for this problem? My only idea right now is HMM, but maybe there are other options? **UPDATE:** After a mini investigation via google.com I've found that best results in problem of binary sequence prediction currently is obtained by evolving state transition rules of finite automata (see this paper). This method allows up to 70% accurate prediction of standard random generator algorithm used in most programming languages. So the question now is how to modify initial method in order to use additional information from **_T_** sequence in my question. The only way I can think of right now is: discretize T and build one more finite state machine which works on **_T_** sequence, then combine both FSMs into one system and try to evolve best set of transition rules for it. But, this idea looks too complicated, and unclear in details right now. There should be something simpler to try first.


title: Conditional Expectation of Multivariate Distributions

text: My book states that $\mathrm{E}[X_2] = \mathrm{E}[\mathrm{E}[X_2\mid X_1]]$ and gives a weird proof and a short paragraph about how awesome this result is and how intuitive this is. I neither understand the proof nor do I understand the "intuition". Can someone give me both? The textbook proof uses some non-standard notation and weird transformations but I am hoping for a clear lucid proof.


title: Skewed variables in factor analysis

text: I want to do principal component analysis (factor analysis) on SPSS based on 22 variables. However, some of my variables are very skewed (skewness calculated from SPSS ranges from 2-80!). So here are my questions:   1. Should I keep the skewed variables like that or could I transform the variables on principal component analysis? If yes, how would I interpret factor scores?    2. What type of transformation should I do? log10 or ln?   3. Originally, my KMO (Kaiser-Meyer-Olkin) is 0.413. Much literature recommends a minimum of 0.5. Can I still do factor analysis, or do I need to remove variables to raise my KMO to 0.5?


title: How to best display forecast deviation?

text: I have to represent some forecasting deviations into a chart. The problem is that I can’t find the proper way to do it. I don't know which approach to use in order fit with the data type I'm working with. The biggest is the traffic the more accurate the forecasts is. The following table represents the forecast deviation and traffic consumption for 5 users. The biggest deviation was made by the first user and not by the second user. How to represent that in a chart?               Deviation (%)      2.6     22.8       2.5     15.9      2.3                  Traffic (GB)      17.7     0.08      15.2      0.1     13.3      I can do the calculus of deviation*traffic/100 and to put those in a chart but maybe there is a better solution. Maybe I have to use other logic?


title: How to present results of a Lasso using glmnet?

text: I would like to find predictors for a continuous dependent variable out of a set of 30 independent variables. I am using Lasso regression as implemented in the glmnet package in R. Here is some dummy code:               # generate a dummy dataset with 30 predictors (10 useful & 20 useless)      y=rnorm(100)     x1=matrix(rnorm(100*20),100,20)     x2=matrix(y+rnorm(100*10),100,10)     x=cbind(x1,x2)          # use crossvalidation to find the best lambda     library(glmnet)     cv <- cv.glmnet(x,y,alpha=1,nfolds=10)     l <- cv$lambda.min     alpha=1          # fit the model     fits <- glmnet( x, y, family="gaussian", alpha=alpha, nlambda=100)     res <- predict(fits, s=l, type="coefficients")     res       My questions is how to interpret the output:   * Is it correct to say that in the final output all predictors that show a coefficient different from zero are related to the dependent variable?    * Would that be a sufficient report in the context of a journal publication? Or is it expected to provide test-statistics for the significance of the coefficients? (The context is human genetics)   * Is it reasonable to calculate p-values or other test-statistic to claim significance? How would that be possible? Is a procedure implemented in R?    * Would a simple regression plot (data points plotted with a linear fit) for every predictor be a suitable way to visualize this data?   * Maybe someone can provide some easy examples of published articles showing the use of Lasso in the context of some real data & how to report this in a journal?


title: (Non-linear) Transformation of confidence interval for multinomial parameters

text: I have a certain computational biology problem I wish to model. Say I have a vector $\vec{f}$ that yields $\vec{p}$, of which the explicit form amounts to picking $\vec{p}$ as an eigenvector from an eigenvalue-problem with the largest associated eigenvalue. $\vec{p}$ is a probability vector, i.e. its components sum to 1. I now have a multinomial distribution with $\vec{p}$ as the parameter. The data I have consists of $N$ samples from this multinomial distribution. As is the case in statistics, I do not know $\vec{p}$ (or $\vec{f}$ for that matter). Say I can construct some confidence interval for $\vec{p}$ of the multinomial distribution (there are a number of methods in literature), how do I transform this confidence interval to $\vec{f}$? The function is defined as $\vec{f} = \vec{h}(\vec{p}) = \left( Q^T diag\left( \vec{p} \right) \right)^{-1} \vec{p}$, where $Q$ is a stochastic transition matrix, i.e. its rows sum to 1. I know that if the transformation is monotonic, the boundaries maybe simply be transformed to yield a confidence interval for $\vec{f}$. So far I have linearised $\vec{h}(\vec{p})$ around $\hat{p}$ and used this with okayish results, but there remain some substantial departures from 'optimal' confidence intervals as observed from non-parametric bootstrap samples. Ultimately, I would like to work in this framework, and if possible not go to Bayesian (or non-parametric bootstrapping), as the dimensionality of the problem will scale up (realistically $\vec{f}$ might be 1000-dimensional or so). Thanks for any suggestions.


title: How to calculate hyperplane for SVM?

text: I am using SVM for document categorization. For that I have training sets i.e. set of training vector with their category. For example I have following sample data(assumption):                  -1 1:1 2:1        +1 2:1 4:1        -1 3:1 5:1        +1 1:1 6:1 Can any one say:      1) How to plot the data points in vector space(Sample diagram for the given test data will help me best)?      2) How to calculate **hyperplane** using the given sample?.      3) How to classify the new document using hyperlane for following data?                        2:1 4:1 4)Whether the kernel function are used for generating hypherlane efficiently? * * *               My sample data:     -1 3:1 11:1 14:1 19:1 39:1 42:1 55:1 64:1 67:1 73:1 75:1 76:1 80:1 83:1      -1 3:1 6:1 17:1 27:1 35:1 40:1 57:1 63:1 69:1 73:1 74:1 76:1 81:1 103:1      -1 4:1 6:1 15:1 21:1 35:1 40:1 57:1 63:1 67:1 73:1 74:1 77:1 80:1 83:1      -1 5:1 6:1 15:1 22:1 36:1 41:1 47:1 66:1 67:1 72:1 74:1 76:1 80:1 83:1      -1 2:1 6:1 16:1 22:1 36:1 40:1 54:1 63:1 67:1 73:1 75:1 76:1 80:1 83:1      -1 2:1 6:1 14:1 20:1 37:1 41:1 47:1 64:1 67:1 73:1 74:1 76:1 82:1 83:1      -1 1:1 6:1 14:1 22:1 36:1 42:1 49:1 64:1 67:1 72:1 74:1 77:1 80:1 83:1      -1 1:1 6:1 17:1 19:1 39:1 42:1 53:1 64:1 67:1 73:1 74:1 76:1 80:1 83:1      -1 2:1 6:1 18:1 20:1 37:1 42:1 48:1 64:1 71:1 73:1 74:1 76:1 81:1 83:1      +1 5:1 11:1 15:1 32:1 39:1 40:1 52:1 63:1 67:1 73:1 74:1 76:1 78:1 83:1      -1 5:1 16:1 30:1 35:1 41:1 64:1 67:1 73:1 74:1 76:1 80:1 83:1      -1 5:1 6:1 15:1 20:1 37:1 40:1 50:1 63:1 67:1 73:1 75:1 76:1 80:1 83:1      -1 5:1 7:1 16:1 29:1 39:1 40:1 48:1 63:1 67:1 73:1 74:1 76:1 78:1 83:1      -1 1:1 11:1 18:1 20:1 37:1 42:1 59:1 62:1 71:1 72:1 74:1 76:1 80:1 83:1      +1 5:1 18:1 19:1 39:1 40:1 63:1 67:1 73:1 74:1 76:1 80:1 83:1                           My Sample model file content:     1.0 3:1.0 11:1.0 14:1.0 19:1.0 39:1.0 42:1.0 55:1.0 64:1.0 67:1.0 73:1.0 75:1.0 76:1.0 80:1.0 83:1.0      1.0 4:1.0 6:1.0 15:1.0 21:1.0 35:1.0 40:1.0 57:1.0 63:1.0 67:1.0 73:1.0 74:1.0 77:1.0 80:1.0 83:1.0      1.0 2:1.0 6:1.0 16:1.0 22:1.0 36:1.0 40:1.0 54:1.0 63:1.0 67:1.0 73:1.0 75:1.0 76:1.0 80:1.0 83:1.0      1.0 5:1.0 6:1.0 15:1.0 20:1.0 37:1.0 40:1.0 50:1.0 63:1.0 67:1.0 73:1.0 75:1.0 76:1.0 80:1.0 83:1.0      1.0 5:1.0 7:1.0 16:1.0 29:1.0 39:1.0 40:1.0 48:1.0 63:1.0 67:1.0 73:1.0 74:1.0 76:1.0 78:1.0 83:1.0      0.47854266714671057 2:1.0 11:1.0 18:1.0 20:1.0 37:1.0 40:1.0 54:1.0 63:1.0 67:1.0 73:1.0 74:1.0 76:1.0 80:1.0 83:1.0      1.0 2:1.0 6:1.0 15:1.0 19:1.0 39:1.0 40:1.0 52:1.0 63:1.0 67:1.0 73:1.0 74:1.0 76:1.0 82:1.0 83:1.0      1.0 3:1.0 6:1.0 14:1.0 22:1.0 36:1.0 40:1.0 54:1.0 63:1.0 67:1.0 73:1.0 74:1.0 76:1.0 80:1.0 83:1.0      1.0 5:1.0 8:1.0 14:1.0 19:1.0 39:1.0 40:1.0 50:1.0 63:1.0 67:1.0 73:1.0 74:1.0 76:1.0 82:1.0 83:1.0      0.26128743031116053 3:1.0 7:1.0 18:1.0 22:1.0 36:1.0 40:1.0 50:1.0 63:1.0 67:1.0 73:1.0 74:1.0 76:1.0 82:1.0 83:1.0      1.0 5:1.0 15:1.0 22:1.0 36:1.0 40:1.0 63:1.0 67:1.0 73:1.0 74:1.0 76:1.0 82:1.0 83:1.0      1.0 1:1.0 15:1.0 21:1.0 35:1.0 45:1.0 62:1.0 68:1.0 72:1.0 74:1.0 77:1.0 80:1.0 93:1.0