title: Correlation coefficient for data table

text: I would like to display correlation coefficients in a table (ideally - with p-value). However, my code produces exactly the same values for each period (so something is obviously wrong). Could you give me any advice:               #first of all, I read my data table from CSV file:     imported <- read.table (file="/home/someone/data_for_R.csv", header=TRUE, sep='\t', quote='"\'', dec=',', fill=FALSE, comment.char="#",  na.strings = "NA", nrows = -1, skip = 0, check.names = TRUE, strip.white = FALSE, blank.lines.skip = TRUE)          # Typing: class(imported[["Period"]]) produces:     # [1] "factor"          #Typing: levels(imported[["Period"]]) produces:     # [1] "Summer 2010" "Summer 2011" "Winter 2010" "Winter 2011" "Winter 2012"          xx <- imported[c("Period","Data1.MEAN","Data2.MEAN")]     result <- by(xx, xx$Period, function(x) {cor(xx$Data1.MEAN, xx$Data2.MEAN)})         result.dataframe <- as.data.frame(as.matrix(result))         result.dataframe$C <- rownames(result)      **EDIT:** Code which reads file from Github:               library(RCurl)     x <- getURL("https://raw.githubusercontent.com/kedziorm/testowe/master/data_for_R.csv")     imported <- read.csv (text=x, header=TRUE, sep='\t', quote='"\'', dec=',', fill=FALSE, comment.char="#",  na.strings = "NA", nrows = -1, skip = 0, check.names = TRUE, strip.white = FALSE, blank.lines.skip = TRUE)     xx <- imported[c("Period","Data1.MEAN","Data2.MEAN")]     result <- by(xx, xx$Period, function(x) {cor(xx$Data1.MEAN, xx$Data2.MEAN)})         result.dataframe <- as.data.frame(as.matrix(result))         result.dataframe$C <- rownames(result)      **EDIT:** This should finally work:               x <- "Period\tDate\tData1.MEAN\tData1.MEDIAN\tData2.MEAN\tData2.MEDIAN\tData3.MEAN\tData3.MEDIAN\nWinter 2010\t26-03-2010\t0,3580917\t0,307479\t0,551191\t0,612853\t0,3476462\t0,3996462\nWinter 2010\t26-04-2010\t0,3016958\t0,2643808\t0,417791\t0,393714\t0,2811050286\t0,3061050286\nSummer 2010\t03-07-2010\t0,1916181\t0,1816603\t0,390925\t0,37385\t0,2183438286\t0,2923438286\nSummer 2010\t04-07-2010\t0,2548711\t0,1738567\t0,4349834\t0,4957131\t0,2467746286\t0,3437746286\nWinter 2011\t01-11-2010\t0,3393042\t0,2870481\t0,497295\t0,538132\t0,3210420857\t0,3690420857\nSummer 2011\t04-06-2011\t0,222748\t0,2218226\t0,363823\t0,275725\t0,2309696\t0,2809696\nSummer 2011\t05-06-2011\t0,241889\t0,1918457\t0,373566\t0,292997\t0,2306573429\t0,2966573429\nWinter 2012\t07-11-2011\t0,2264874\t0,2601413\t0,373048\t0,274139\t0,2456219143\t0,2756219143\nWinter 2012\t08-11-2011\t0,2414665\t0,2662565\t0,314382\t0,279857\t0,2348871429\t0,2598871429\nWinter 2012\t09-11-2011\t0,2817838\t0,2325952\t0,376063\t0,468148\t0,254412\t0,287412\nWinter 2012\t10-11-2011\t0,2476841\t0,2667485\t0,406902\t0,476582\t0,2632384571\t0,3632384571\n"     imported <- read.csv (text=x, header=TRUE, sep='\t', quote='"\'', dec=',', fill=FALSE, comment.char="#",  na.strings = "NA", nrows = -1, skip = 0, check.names = TRUE, strip.white = FALSE, blank.lines.skip = TRUE)     xx <- imported[c("Period","Data1.MEAN","Data2.MEAN")]     result <- by(xx, xx$Period, function(x) {cor(xx$Data1.MEAN, xx$Data2.MEAN)})         result.dataframe <- as.data.frame(as.matrix(result))         result.dataframe$C <- rownames(result)


title: Performing multiple comparisons on linear mixed model

text: I have collected data on gas fluxes from plots of soil subjected to 5 different treatments ("D2", "K2", "M", "N", and "O2"), which also possessed variable clay contents. The experiment was laid out in a randomized complete block design, with 4 replications. Within each plot, two separate measurements of flux were performed. 4 blocks x 5 treatments x 2 repeated measurements for 40 total observations I have built a linear mixed model with nlme, and performed multiple comparisons with Tukey's test, as shown below. The data.frame cum_fluxes can be downloaded from: https://www.dropbox.com/s/e58k4vdnevw2fsm/cum_fluxes               require(multcomp)     require(nlme)       cum_fluxes <- dget("cum_fluxes")     n2o_flux.lme <- lme(cum_flux_n2o_ln ~ treatment * clay + block, random = ~1|subsampling, data = cum_fluxes)     aov.n2o <- anova(n2o_flux.lme)     comp.treat <- glht(n2o_flux.lme, linfct=mcp(treatment="Tukey"))     print(summary(comp.treat))      I have two questions. On the one hand, the resulting ANOVA table looks like this (significant effect of treatment with p = 0.0279):                              numDF denDF   F-value p-value     (Intercept)        1    20 14815.691  <.0001     treatment          4     7     5.281  0.0279     clay               1     7     1.326  0.2874     block              3     7     0.674  0.5951     treatment:clay     4     7     0.810  0.5568      However, there are apparently no significant differences found when considering pair-wise comparisons of the five treatments. Did I make a mistake somewhere, or am I interpreting the results incorrectly? Second, I get the following warning message when I run the glht() command:               Warning message:     In mcp2matrix(model, linfct = linfct) :       covariate interactions found -- default contrast might be inappropriate      This makes some sense to me: since I allow for interaction between a covariate (quantitative regressor – i.e., clay content) and a factor (so, non-parallel lines), the difference between two factor levels depends on the value of the covariate. However, I can't figure out how and to what value I should specify the covariate when I compare the factor levels?


title: A concise way to describe the distribution of words in a text

text: I have a text, and I count, for each distinct word, the number of times that word appears in the text. In some texts, there are many rare words, so the distribution will have a long "tail" of words that appear once. In other texts, there are common words that usually appear more than once, so the "tail" will be shorter. What is a standard measurement that describes this difference?


title: Unbalanced repeated measures

text: I'm new to using R and therefore I have some doubts. I have a data table with measures gland secretion of chickens at two ages different, sometimes we just need a measure therefore is not balanced. Also usually we caught two chickens per nest but sometimes we just have a chicken. My question is whether the secretion of the gland is affected by age, for it following the models that I found online, I have done the following and let me know if it is correct. First the data table I entered would be:               NEST        Chick          Secretion             Age     257             1             1.429               5     257             2             1.890               10     127             1             1.567               5     127             1             1.900               10     127             2             1.320               5     127             2             2.000               10          > attach (data)     > NEST = factor (NEST)     > AGE = factor (Age)     > library (nlme)     > lme (Secretion ~ AGE, random = 1 ~ | NEST)      Is my model correct?


title: Why can't I run simply asymmetric garch in the multivariate garch dcc analysis?

text: I have a question that I firstly ran the MGARCH DCC model with GARCH(1,1) specification for all equations, and it worked well. Then, I would like to run MGARCH DCC model with SAARCH(1) specification for all equations. But STATA says: > option saarch() not allowed I don't know why since earch worked well for univariate GARCH modelling for variable 1 and 2. Could anybody help me with this? Thank you very much! The stata output is:               mgarch dcc (var1 = , arch(1) saarch(1)) (var2 = , arch(1) saarch(1)), distribution(t)     option saarch() not allowed


title: K-means clustering for usage profiling

text: I am trying to use k-means clustering to profile mobile device usage behaviour for IT users. My data consists of different system and user level variable/readings like number of calls/sms, cpu/memory usage, number of users and system applications/services etc. The readings are taken every 5 minutes. The idea I have is to use say 1 month's data for training, i.e. clustering, and then use the future data to compare with existing clusters and try to find (dis)similarity between the two. The assumption is different users will have different usage; hence readings from USER B will not fit into clusters from USER A. Now two questions I have:   1. After training (clustering), how do I compare new data with existing clusters to determine (dis)similarity, i.e. new data belongs to same user or not? I am thinking of finding nearest cluster and then checking if the point lies within this cluster's boundary.    2. I am using Silhouettes plot to determine the clustering quality. I get some negative values e.g see![one of the plots here](http://i.stack.imgur.com/OvTiB.jpg). Should I be concerned? or is it normal to have some negative values?


title: Feature selection in the training set

text: I have a classifier, and I am using leave one out cross-validation to assess its performance. On each iteration, I divide the dataset into training and testing sets. The testing set is just the subject that I am going to evaluate (leave one out). Now, I **divide the training set into folds,** and I do feature selection like this: I run my **filter** feature selection algorithm on every fold. When I am done, I have a voting algorithm to obtain the final set with the features that were selected in each fold. I understand that this procedure is adequate when you have a small sample like in my case (subjects = 30, features = 960). **My question is why, if at all, would it be a bad idea to do feature selection on the whole training set instead of dividing it into folds?**


title: Finding the rejection region of a given case

text: > A casino dealer suspects that his die isn't a fair die. He thinks that the > die gives more chances to the even numbers (the probability for an even > number is larger than $0.5$). To test his suspects he decides to have an > experiment, he writes down the outcome of $60$ independent die rolls. The > dealer would say that the die is unfair if the number of even numbers in > those $60$ rolls will be larger than $35$. I'm asked to:   1. Formalize the dealer's assumption.   2. What is the rejection region ($C$)?   3. What is the probability for both mistakes ($\alpha$ and $\beta$)? I started by saying that $X$, the number of even results of a fair cube $X\sim B(\frac{1}{2}, 60)$, and thus $X \sim N (30, 15)$ $H_0: \mu = 30$   $H_1: \mu > 30$ So this is the dealer's assumption. The rejection area should be derived from the fact that the dealer will reject $H_0$ if and only if the even numbers on the experiment will be less or equal to $35$. I'm having a hard time finding $C$ and calculating $\alpha$ and $\beta$.


title: Introductory reading on Copulas

text: For some time now, I have been looking for a good introductory reading on Copulas for my seminar. I am finding lots of material that talk about theoretical aspects, which is good, but before I move onto them I am looking to build a good intuitive understanding on the topic. Could anyone suggest any good papers that provide a good foundation to a beginner (I have had 1-2 courses in statistics and understand marginals, multi-variate distributions, inverse transform, etc., to a reasonable extent)?


title: What is the most appropriate test for this situation?

text: I want to statistically test the mean difference of a continuous variable (let's call this variable 1) between four groups (nominal) in my sample. Variable 1 has a very skewed distribution. This test also needs to be corrected for three other continuous variables that are certainly confounding variables when it comes to explaining the "variable 1" differences between the four groups. Many thanks!


title: How to design and implement an asymmetric loss function for regression?

text: ## Problem In regression one usually computes the **mean squared error** (MSE) for a sample: $$ \text{MSE} = \frac{1}{n} \sum_{i=1}^n\left(g(x_i) - \widehat{g}(x_i)\right)^2 $$ to measure the quality of a predictor. Right now I'm working on a regression problem where the goal is to predict the price customers are willing to pay for a product given a number of numeric features. If the predicted price is too high no customer will buy the product, but the monetary loss is low because the price can simply be decremented. Of course it should not be too high as then the product may not be bought for a long time. On the other hand if the predicted price is too low, the product will be bought quickly without the chance to adjust the price. In other words the learning algorithm should predict slightly higher prices which can be decremented if necessary rather than underestimating the true price which will result in an immediate monetary loss. ## Question How would you design an error metric incorporating this cost asymmetry? * * * ## Possible Solution A way to define an asymmetric loss function would be to simply multiply by a weight: $$ \frac{1}{n} \sum_{i=1}^n \left| \alpha - \mathbb{1}_{(g(x_i) - \widehat{g}(x_i)) < 0} \right|\cdot \left(g(x_i) - \widehat{g}(x_i)\right)^2 $$ with $\alpha \in (0,1)$ being the parameter we can adjust to change the degree of asymmetry. I've found it here. This seems like the most straight forward thing to do, while maintaining the quadratic loss.


title: Restrictions on the type of HMM observations

text: I am working on developing an HMM (or DBN) to detect vigilance from time- series observations of eye-closure. Vigilance is defined as a binary variable (vigilant or non-vigilant). While I understand the general framework and assumptions of HMMs I am a bit confused about the notion of observations or emissions. Specifically what can be used as an observation? My confusion stems from the fact that each individual observation provides very little information about the vigilance. In light of this I would like to consider the following types of features:   1. Mean or SD of eye-closure over a short time window   2. Spectral components (from a DFT) over a short time window   3. The output of a static ML algorithm (e.g. a Decision Tree) trained to predict vigilance based on eye-closure observations I am wondering if each of these types of features would be acceptable? and How the use of these features changes the estimation of the HMM parameters? I would particularly appreciate pointers to any relevant literature. Thank you


title: How to replicate the results from Australian Institute of Sport?

text: I am taking a little data analysis course using SAS software. I was pretty much told to generate the results using the data and the description provided in a specific StatSci page which I am going to also copy- paste here: > **Australian Institute of Sport** > > Keywords: multiple regression, one or twoway analysis of variance, factors > > **Description** Data on 102 male and 100 female athletes collected at the > Australian Institute of Sport, courtesy of Richard Telford and Ross > Cunningham.               Variable    Description          Sport       Sport     Sex         male or female     Ht          Height in cm     Wt          Weight in kg     LBM         Lean body mass     RCC         Red cell count     WCC         White cell count     Hc          Hematocrit     Hg          Hemoglobin     Ferr        Plasma ferritin concentration     BMI         Body mass index = weight/height^2     SSF         Sum of skin folds     %Bfat       % body fat      > **Source** Richard Telford and Ross Cunningham, Australian National > University. Cook, R. D., and Weisberg, S. (1994). An Introduction to > Regression Graphics. Wiley, New York. > > **Analysis** This works well as an example of a one-way or two-way layout by > sport and sex, as well as an example of multiple regression. (I have used > log-Ferr.) But I am a little confused and cannot really ask anyone else since it's the beginning of the weekend. What does the analysis part mean? Do they try to predict Ferr using Sport and Sex? What meaningful results should I try to obtain from this data? I'd really appreciate the ideas and maybe the SAS software procedures involved which I'd try to further implement and tailor for this analysis.


title: Graphical comparison of regression models

text: Will a graph on predicted and measured values plotted for two models separately be helpful in comparing them?


title: Drawing data from "population" for regression analysis

text: We have a response variable $Y$ and predictor $X$, and we draw $n$ samples $(Y_1,X_1), \ldots, (Y_n, X_n)$ from the population of interest to do a regression analysis. Under the assumptions of a simple linear regression model, my question is a conceptual one: how do we really think about the response on the $i$th unit, $Y_i$? Do we say it's drawn from the level or subpopulation of individuals with $ X = x_i $, or from the aggregate population over all the values of $X$? Moreover, while we assume that the response $Y$ in every subpopulation defined by $X$ is normal with equal variances, how do we think about the aggregate population from which $Y_i$ is drawn?


title: Machine learning algorithms/approaches for class recommendations?

text: I am asking a theoretical question about machine learning in terms of clustering. Is it possible, given a set of data of classes that students have taken in a semester to recommend additional classes that students should take if they selected some classes? I am thinking along the line of forming clusters of classes and figuring out if a particular set of picked classes match with a pre-existing set of classes. Then, recommend the class that are in the set. But I am new to machine learning, and so welcome any other suggestions of algorithms. In addition, this is not particularly theoretical, so feel free to ignore: but does anyone know any particular software that can accomplish this? I know LensKit is a software to handle recommendations but it seems to need ratings (which I do not have). I welcome any mathematical manipulations that can turn clusters into "ratings." Thanks.


title: Comparing different classifiers (using ski-kit cross validation values)

text: Thanks for taking the time to read this I'm new to Machine Learning and so am going through a Kaggle competition to help me improve but I have a question. How can I compare different classifiers?? My Python isn't as neat as I would like it to be but i think it's correct. Please let me know if I'm doing something strange               train = pca(train)     cfr = KNeighborsClassifier(n_neighbors=neighbours, algorithm="kd_tree")     cv = cross_validation.KFold(len(train), n_folds=10, indices=True)          results = []     i = 0     count =0     for traincv, testcv in cv:         ClassPred = cfr.fit(train[traincv], target[traincv]).predict(train[testcv])          for j in range(0,(len(train)/10)):             labelindex = testcv[j]             if (ClassPred[j] == target[labelindex]):                 i = i+1         accuracy = (i/2000.0)*100         i=0          results.append(accuracy)         count = count + 1         print "accuracy for fold", count, " : ", accuracy,"%"         print "time after fold", count, " : ", elspasedfold,"%"          elapsed = (time.clock() - start)     #print out the mean of the cross-validated results     print "Results for RBF, c = 10.0, gamma=0.1 \n" + str(np.array(results).mean())     print "Time taken is %ds" % elapsed`.      As you can see I'm doing 10 fold cross validation on the training data and hopefully producing an 'accuracy' value out of the other end. My question here is how can I compare different classifiers?? If I choose cfr as a different classifier (e.g. Random Forest) and get a value for that how do i statistically compare the two? My initial thoughts are to use a two-tailed t-test on the values for each fold (so 10 values per classifier) instead of just taking using the one value it currently outputs, which I think is an average of all the folds combined. I would need to then get a p-value to see if the differences are significant. I am not sure how I would go about implementing this however or if this is the correct thing to do. My stats is a little patchy (Which I'm working on) but any help anyone can give me would be much appreciated. So to clarify, using python and sci-kit what would be the best way to compare the performance of different classifiers on my data? Thanks! EDIT: The tutorial assessment criteria is that 10 fold cross validation must be used.


title: Bayesian vs OLS . Degrees of freedom. do we need to bother about d.o.f for bayesian estimation?

text: Im new to bayesian and still trying to figure out the process in detail. I read somewhwere that one of the advantages of bayesian estimation is that it helps save on degrees of freedom and can be used when there are lesser number of data points and more number of parameters. Can someone through some light on this? i know its got something to do with MC MC but not very clear on how that helps.


title: Compute probability

text: Suppose that $X$, $Y$ and $Z$ are $\text{i.i.d.} \sim \text{Uniform}(0,1)$. Let $t > 0$ be a ﬁxed constant. (i) Compute $P(X/Y \leq t)$   (ii) Compute $ P(XY \leq t)$   (iii) Compute $ P(XY/Z \leq t)$ I found the solution for (i) part undermining the different values of $t$.   $ P(X/Y \leq t) = \int_0^1\int_0^{yt}dx dy = t/2 $ ; when $t\leq 1$ I am finding hard time to lake limits for different values of $t$ i.e when $t\leq 1 $ and $t>1$


title: Test for differences between (among) related, but not matched, samples

text: When two samples are related, or dependent, but the observations are not matched, are there any tests that will determine if the samples (means or otherwise) are different? I've searched extensively and have only found tests for matched samples, which is not what I need.


title: Why are Dunnett's p-values smaller than those from two-sample t-tests?

text: Since Dunnett corrects for multiple comparisons I would assume the p-value would always be higher than a direct two-sample t-test. I came across a situation where this was not the case. Possible reasons why? Example:               library(multcomp)     a1=data.frame(group=rep(c('a','b','control','d'),each=3),                   y=c(13.3,13.8,13.6,10.9,14.4,11.1,18.8,15.2,16.6,1.6,1.2,.9))     a1$group=relevel(a1$group,ref='control')     fit=aov(y~group,data=a1)     summary(fit)     test.dunnett=glht(fit,linfct=mcp(group="Dunnett"))     summary(test.dunnett)               t.test(a1$y[a1$group=='control'],a1$y[a1$group=='a'])$p.value     t.test(a1$y[a1$group=='control'],a1$y[a1$group=='b'])$p.value     t.test(a1$y[a1$group=='control'],a1$y[a1$group=='d'])$p.value      Output:                           Df Sum Sq Mean Sq F value   Pr(>F)         group        3  413.0  137.67   74.99 3.37e-06 ***     Residuals    8   14.7    1.84                         Simultaneous Tests for General Linear Hypotheses          Multiple Comparisons of Means: Dunnett Contrasts               Fit: aov(formula = y ~ group, data = a1)          Linear Hypotheses:                      Estimate Std. Error t value Pr(>|t|)         a - control == 0   -3.300      1.106  -2.983  0.04298 *       b - control == 0   -4.733      1.106  -4.279  0.00704 **      d - control == 0  -15.633      1.106 -14.131  < 0.001 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1     (Adjusted p values reported -- single-step method)          [1] 0.08504438     [1] 0.03780491     [1] 0.00342318      The Dunnett p-values (.04, .007, .001) are smaller compared to the t-test p-values (.085, .038, .003).


title: Saved glmer object in R does not show up correctly in R

text: I've run a mixed logit model with glmer and saved the resulting glmer object to disk as an .Rdata file. But when I load it and print it using print(myglmerobject,corr=0) I see a list of variables, correlations, etc. rather than the nice piece of results this command generates. Any ideas how I can fix this? It's a big dataset and it takes me hours to run the model again.


title: What is the probability of n people from a list of m people being in a random selection of x people from a list of y people?

text: If I am selecting 232 people from a pool of 363 people without replacement what is the probability of 2 of a list of 12 specific people being in that selection? This is a random draw for an ultra race where there were 363 entrants for 232 spots. There is an argument about whether the selection was biased against a certain group of 12 people. My initial attempt at calculating this was that there was 232 choose 363 possible selections. The number of combinations of any one person from the list of twelve is 1 choose 12 + 2 choose 12 + ... + 11 choose 12 + 12 choose 12. Thus 1 choose 12 + 2 choose 12 .... / 232 choose 363. Which ends up being a very low number which is clearly too low. How do I calculate this?


title: Probability of choosing N out of K correctly from a total population of X

text: Out of a total population of 50 items, there are 20 "correct" items, what is the probability of me selecting 11 (or more) of these correct items if I randomly choose 20? How can I do this in R? I was thinking of using `binom.test()` but that assumes replacement, and my questions involves 'no-replacement' when making your selections.


title: In Kneser-Ney smoothing, how are unseen words handled?

text: From what I have seen, the (second-order) Kneser-Ney smoothing formula is in some way or another given as $ \begin{align} P^2_{KN}(w_n|w_{n-1}) &= \frac{\max \left\\{ C\left(w_{n-1}, w_n\right) - D, 0\right\\}}{\sum_{w'} C\left(w_{n-1}, w'\right)} + \lambda(w_{n-1}) \times P_{cont}(w_n) \end{align} $ with the normalizing factor $\lambda(w_{n-1})$ given as $ \begin{align} \lambda(w_{n-1}) &= \frac{D}{\sum_{w'} C\left(w_{n-1}, w'\right)} \times N_{1+}\left(w_{n-1}\bullet\right) \end{align} $ and the continuation probability $P_{cont}(w_n)$ of a word $w_n$ $ \begin{align} P_{cont}(w_n) &= \frac{N_{1+}\left(\bullet w_{n}\right)}{\sum_{w'} N_{1+}\left(\bullet w'\right)} \end{align} $ where $N_{1+}\left(\bullet w\right)$ is the number of contexts $w$ was seen in or, simplier, the number of distinct words $\bullet$ that precede the given word $w$. From what I've understood, the formula can be applied recursively. Now this handles known words in unknown contexts nicely for different n-gram lengths, but what it doesn't explain is what to do when there are out-of- dictionary words. I tried following this example which states that in the recursion step for unigrams, $P_{cont}(/) = P^0_{KN}(/) = \frac{1}{V}$. The document then uses this - quoting Chen and Goodman - to justify the above formula as $P^1_{KN}(w) = P_{cont}(w)$. I fail to see how it works out in the presence of an unknown word $w = \text{unknown}$ though. In these cases $P_{cont}(\text{unknown}) = \frac{0}{\text{something}}$ since, obviously, the unknown word doesn't continue anything regarding the training set. Likewise the count of n-grams is going to be $C\left(w_{n-1}, \text{unknown}\right) = 0$. Furthermore, the whole $\sum_{w'} C\left(w_{n-1}, w'\right)$ term might be zero if a sequence of unknown words - say, a trigram of OOD words - is encountered. What am I missing?


title: Variable selection: Why certain categories are chosen but not others?

text: I'm doing variable selection using the Lasso. To explain my response variable I have several predictors, both categorical and numerical, but I have problems to explain the process that underlies when Lasso selects only a category from a variable with several categories. For example, one of my predictors is a categorial variable with four levels, and Lasso just selects one of them. So, Lasso is working with the whole variable (the four categories) but some may "enter" and some not. How can I explain this? It's something related with Analysis of Covariance? I hope my question makes sense and I would appreciate a not very mathematical answer. Thanks


title: Compute the cumulative hazard of a time interval

text: I'm a little confused in calculating the cumulative hazard within a time interval. We know that $H(t)=\int^{t}_{0}h(u)du$, if I have $\Delta t=t_1-t_0$ (1) $H(\Delta t)=\int^{\Delta t}_{0}h(u)du$ (2) $H(\Delta t)=\int^{t_1}_{0}h(u)du-\int^{t_0}_{0}h(u)du$. Which is the correct one? (1) or (2)? I felt the (1) may be the correct.


title: Extreme Value Theory and heavy (long) tailed distributions

text: I'm analyzing data about which I have a strong suspicion that it is self- similar (Hurst parameter ranging from 0.60 to 0.78 depending on estimation method and sample sequence). I also observe high realization values much more often (compared to experiments where no self-similarity is suspected) which suggests that the generating distribution can have a long (heavy) tail. It seems that Extreme Value Theory cannot be used in this case to reason about the tail behavior. In case of self-similarity long range dependence does obviously violate the i.i.d. requirement, but does the inability to use EVT in my case apply more broadly to all long tailed distributions?


title: Hidden Markov Model with MFCC coefficients

text: I don't know whether this is the correct forum for this but here goes: I'm trying to implement a Hidden Markov Model to be able to predict and find the best sequence/path for a training file. So far, I have the (MFCC) Coefficients of a signal, and I am looking to train this so I can then compare (two data sets) using the Viterbi algorithm to find the best path. But I have some confusions:   1. If I estimate $\Pi$, can I then use the backward-forward algorithm to find the probabilities at a given state?   2. Instead of using the backward-forward algorithm, can I use the Viterbi algorithm to find the probability?    3. If I'm trying to identify two probable things (two outcomes) and `N = 2` what would `M` be? I assume that `N` = states and `M` = observations? But, the vector I have contains: N = 13, M = 450, so these values cannot therefore be used for N->M for the HMM.    4. Do I therefore send the training data to the forward-backward algorithms, which predicts the probability for each state, then this will give me a final output... Which I can then compare with the Viterbi de-coder?


title: Logistic regression, Chi-square, and study design

text: I have a study in which I have developed a new predictor (binary) for a disease (also a binary variable). The study has two parts. In the first part, I want to test if my predictor is strongly associated with disease. I am planning to use a chi-square test on the predictor vs. disease contingency table (2x2) for this first part. In the second part, I want to test if my predictor is complementary to existing predictors (which are binary or continuous). To do this, I am planning to compare 2 nested logistic regression models to predict disease: model 1 with my predictor & existing predictors, and model 2 with only existing predictors. I will use likelihood ratio test or Akaike information criterion for the comparison. My main question is: should I be using logistic regression for the first part also, instead of a chi-square test? Or, is chi-square test more powerful than logistic regression to test association in 2x2 contingency tables? These questions are related to a previous thread, but my questions were not fully answered there. Also, is logistic regression the best way to test my hypothesis in the second part of the study? Finally, if the answer to the test in the second part of my study is yes, will the first part become too redundant to be included in the study? Thank you.


title: Plotting a histogram estimate together with an analytical solution of the mean and covariance parameters

text: How do I plot a histogram estimate of this marginal distribution p(x1):               import numpy as np     import matplotlib.pyplot as plt     from scipy.stats import norm     linalg = np.linalg          N = 100     mean = [1,1]     cov = [[0.3, 0.2],[0.2, 0.2]]     data = np.random.multivariate_normal(mean, cov, N)     L = linalg.cholesky(cov)     # print(L.shape)     # (2, 2)     uncorrelated = np.random.standard_normal((2,N))     data2 = np.dot(L,uncorrelated) + np.array(mean).reshape(2,1)     # print(data2.shape)     # (2, 100)     plt.figure()     plt.scatter(data2[0,:], data2[1,:], c='green')         plt.scatter(data[:,0], data[:,1], c='yellow')     plt.show()          # Plotting histograms and fitting normal distributions     plt.subplot(211)     plt.hist(data[:,0], bins=20, normed=1, alpha=0.5, color='green')     plt.hist(data2[0,:], bins=20, normed=1, alpha=0.5, color='yellow')     x = np.arange(-1, 3, 0.001)     plt.plot(x, norm.pdf(x, *norm.fit(data[:,0])), color='green')     plt.plot(x, norm.pdf(x, *norm.fit(data2[0,:])), color='yellow')     plt.title('Var 1')          plt.subplot(212)     plt.hist(data[:,1], bins=20, normed=1, alpha=0.5, color='green')     plt.hist(data2[1,:], bins=20, normed=1, alpha=0.5, color='yellow')     x = np.arange(-1, 3, 0.001)     plt.plot(x, norm.pdf(x, *norm.fit(data[:,1])), color='green')     plt.plot(x, norm.pdf(x, *norm.fit(data2[1,:])), color='yellow')     plt.title('Var 2')          plt.tight_layout()      together with an analytical solution given the mean and covariance parameters: N = 100 2-dimensional samples x = (x1,x2)T ∈ R2 drawn from a 2-dimensional Gaussian distribution, with mean µ = (1,1)T and covariance matrix Σ = (0.3 0.2 0.2 0.2) I'm told that I can use this from Bishops Pattern Recognition and Machine Learning (taken from page 88):               p(xa) = p(xa, xb) dxb      Can you write me the analytical expression for the marginal distribution p(x1)? I'm using Python.


title: Regarding topic candidate list for topic modeling analysis (LDA)

text: I am confused on one issue of topic modeling analysis. When we perform the topic modeling analysis over a collection of documents, do we need to provide a set of topic candidate list as the input for the analysis? Or the resulting topics are just extracted from the input documents themselves?


title: Fitting the inverse Weibull distribution to data in R

text: I am trying to fit the Weibull distribution to a dataset in R. I have a dataset, "loss", as a .csv file, and use `fitdistr(loss,"invexp")`. I received the message: > Failure infitdistr(loss, "invweibull") : unsupported distribution How can I fit the inverse Weibull Distribution to a dataset in R? * * * Thank you! I am wonder both how to get this done in R and why I have an error when using a function fitdistr to make maximum-likelihood estimation. My error message is about function "invweibull" in R. The same problem I have for the frechet distribution. My dataset is the following: > 336.893 2.468.001 572.651 426.581 501.593 433.167 708.431 736.761 1.405.323 > 2.161.790 492.029 1.429.564 585.885 563.030 408.345 4.840.762 308.791 > 1.201.796 1.463.367 849.336 510.141 282.088 648.475 699.729 849.677 606.299 > 380.988 I am not sure if I can simply fit a Weibull to the inverse of the observations and say then that parameters are the same for inversed weibull. I use paket "actuar" and function "invweibull" and "frechet" and become a failure unsupported distribution. Simple question is, why the estimations can not be done for this dataset and may be somebody know how is it possible to fit inverse weibull and frechet distributions in R? Thank you very much for your respond to my question!


title: "autocorrelation" of nominal variable

text: I have data on patients attendance of a series of sessions of counseling. At each session, a patient could be A, S or D. Each patient has data for 12 sessions, and there are ~100 patients. I am interested in how the variables are 'autocorrelate'. I could just make 3x3 tables of (e.g.) time 1 vs. time 2, time 1 vs. 3, time 2 vs. 3, etc. but I am curious if there is something more sophisticated. Googling this revealed some measures for spatial autocorrelation, e.g join- counts, but nothng specifically for temporal autocorrelation. I could also try to adapt the join counts to temporal data, but, again, am curious if someone has done this already. I have access to R and SAS packages


title: Significant difference between two groups in a right /wrong quiz

text: I did an exam, which the comment given was that my answers to two questions were too "short" or lacking… (I do not know what was lacking). I would like to learn and understand. As I remember I think I used something called an on line calculator for categoric values when I came to my answers. The exam questions were:   1. Students representing two different geographical ares had a quiz test.      * group 1 had 36 right and 13 wrong     * group 2 had 29 right and 19 wrong Is there a significant difference between these groups (p<0,05)?   2. What if the quiz result where like this:     * group 1 360 right, 130 wrong     * group 2 290 right, 190 wrong Is there a significant difference between this groups (p<0,05)? My answers were:   1. Yes there is a significant difference   2. No there is not a significant difference. My questions to you are:   * Are my answers correct?   * And what is the proper way of showing how you get to this/ or correct answer (without using a calculator)? Meaning, what are the steps I need to take to solve this kind of questions? I am most grateful if anyone have time to explain this to me.


title: Implementing EM clustering for a mixture markov model

text: I have a mixture Markov model (containing K clusters, or components) that I am trying to train, e.g perform clustering over a set of varying length sequences. Each component of the model is a first order Markov chain. The model is defined by the following: ![enter image description here](http://i.stack.imgur.com/O4TuQ.gif) where ![enter image description here](http://i.stack.imgur.com/HjzXm.gif) is the model describing behaviour of sequences in the k-th cluster. Each kth model component(cluster) has initial probability distribution vector (theta_KI) and a transition matrix (theta_tK) - see equation 2. Each sequence seq_i in the set of sequences has a different length L (so the L in the 2nd equation would be different for each sequence in the set). There is n possible visible states in the model and m sequences in each cluster, and the total amount of sequences is M. In this case, how would I define e- and m-steps for my EM algorithm? My understanding is that I have to maximize log likelihood function of the model, which would be a sum of likelihoods across all vectors v_i ( synonym with sequence seq_i)- something like: ![enter image description here](http://i.stack.imgur.com/05Mmi.gif) However, I get stuck not knowing if this is the right definition for what needs to be maximized and which arguments need to be found for that. I have looked for literature but it is very sparse, almost nonexistent. Another question: is there actually an R library that has functions that can be extended to include such a model for EM algorithm?


title: Interpreting interaction between dummy IV and continuous moderator with log DV

text: I urgently need your help to read the results of my regression analysis for my master thesis, which I need to hand in next week. My DV is the natural logarithm of R&D expenditures, IV is a dummy variable labelled **Decline** , where 1= Firm is in a decline and 0= Firm is not in decline. The moderator is managerial ownership measured in percentage (0-100% or 0-1.0), labelled **MOWN**. My thesis has two goals: 1) Find the relationship between IV (Decline) and DV (R&D Expenditures, log) 2) Find the moderating effect of MOWN on the relationship between IV-DV I regressed Decline, MOWN, and Decline X MOWN simultaneously and got following coefficients. $log(RD) = 6.984 - 0.852Decline - 4.703Mown + 3.030MownDecline$ The coefficient for Decline is **-0.852***** , which means that _firms invest 85% less in R &D expenditures when Decline=1_ (correct?) Coefficient for MOWN: **-4.703* _**Coefficient for MOWN x Decline: **3.030_** Contant: 6.984*** How do I interpret the results of the interaction now? When a firm is in decline (Decline=1), _firms will invest 303% more in R &D expenditures with a one-percentage point increase in MOWN?_. Secondly, can I treat the MOWN variable as additional IV as well? Saying that MOWN, the moderator, has also a direct impact in the DV? Saying that _"a one- percentage point increase in MOWN decreases DV by 470%"_ Thank you a lot for your help. I really appreciate it!


title: Weibull and Poisson proportional hazard test - survival data Stata

text: In a Cox model the proportional hazard model can be tested by for example:                stset end, failure(x==1) enter(start) origin(start) scale(365.25)      stcox a b c d e        estat phtest      stphtest, detail /*test for proportional hazards*/      stphplot, by(a) *plotting the different hazards by time       stphplot, by(b)      stphplot, by(c)      stphplot, by(d)      stphplot, by(e)      My question is if there is any similar test, except for the `stphplot` I could do for either a Weibull or an exponential regression in Stata to check if the proportional hazard assumption holds? I started like:                stset end, failure(x==1) enter(start) origin(start) scale(365.25)      streg a b c d e, dist(weibull)        stphplot, by(a) *plotting the different hazards by time       stphplot, by(b)      stphplot, by(c)      stphplot, by(d)      stphplot, by(e)      Thanks for your help! Linda


title: What is the standard procedure for evaluating a user-based CF algorithm with a dataset offline?

text: I have read some papers and other materials about the evaluation of recommender systems (RS). Most of them discuss the various properties of RS (e.g. accuracy, diversity, etc.), and metrics for different tasks (e.g. RMSE, precision, recall, etc.), and some protocols. But I am still not very clear about data partitioning and the detailed validation procedure. [Shani, Guy, and Asela Gunawardana. "Evaluating recommendation systems." Recommender systems handbook. Springer US, 2011. 257-297.] For example, I want to evaluate a basic user-based CF algorithm with rating feedback. The data format is user-item-rating. There're 100 users in the dataset. Now I want to use 5-fold cross validation method. This is my way:   1. According to users, I randomly separate the original dataset into 5 folds. Each fold contains 20 users.   2. For all users I set one same time instant, so that before and after this time instant, each user contains a number of ratings. The data of any user x before the instant is for creating user profile, the data after the instant is for validating rating prediction (user x as user being predicted) or just being hidden(user x in neighborhood selection).   3. For each user i in all 5 folds, I create a user profile represented by rating vector{} for user i.   4. For fold-1, for each user i in fold-1 I calculate user similarities between user i and the users in the fold-2 ~ fold-5, and set an appropriate value of k-NN to obtain the neighbourhood for each user.   5. For each user i in fold-1, I get predicted ratings from his neighbourhood users (items no t used by the tested user) by an appropriate rating averaging equation.   6. calculate prediction error for user i.   7. Repeat step 2~6 for each fold (5 times all) and get the average prediction error. This is one round.   8. Use different parameter combinations(e.g. k-NN neighborhood size), repeat step 2~7 with many times, so as to obtain the lowest prediction error (meaning best parameter combination).   9. With the best parameter combination, I get the lowest prediction error, and this is the final evaluation result. I have two questions:   * Is this procedure **correct** for evaluating a user-based CF algorithm?    * If so, what are the specific data in this procedure corresponding to **the concepts of "training set", "testing set" and "validation set"** from Statistics?


title: Interpreting discrimination values in IRT polytomous (grm) models. Are there any cut-offs or anchors?

text: I am working on a project that is looking at the discriminative properties of 6, 5-point Likert items in a scale. I am using the ltm package in R to examine this within a item response theory (irt) paradigm, using a graded response model (grm). For each item, I am given a "discrimination value" that corresponds with the slope of the ICC at each level given the ability level. For my six items, the values are (approximately): 2.1, 1.8, 1.3, 0.8, 2.7, 2.3. How would I interpret these? Is there a cut-off for "good" discriminative properties? Are there values related to "weak", "medium", or "strong" levels of discrimination? Thanks for your help!


title: How do I compute the density of this data set that is made up of two different 3D-distributions?

text: A sequel to this question. I have a dataset where:   * $\frac{4}{5}$ of the points are drawn from: $(x, y) \sim \mathcal{U}_{2}(0,30)$, $(z) \sim \mathcal{U}_{1}(14.5, 15.5)$.   * $\frac{1}{5}$ of the points are drawn from: $(x, y, z) \sim \mathcal{U}_{3}(0,30)$ Where $\mathcal{U}_{d}(x,y)$ is to be interpreted as a $d$ dimensional set of points which are in each dimension drawn from the range between $x$ and $y$. ## The implementation I have implemented this in matlab like this: General Init:               dim = 3;          uniP = 1/5;     wallP = 4/5;          uniformN = ceil(N * uniP);     wallN  = ceil(N * wallP);      First distribution (wall):               % parameters     lowerWall = [0,0,14.5];     upperWall = [30,30,15.5];          % values     [wallD] = blockUniformDist(lowerWall, upperWall, wallN, dim);      Second distribution (noise):               % parameters     lower = 0;     upper = 30;          % values     [uniformD] = uniformDist(lower, upper, uniformN, dim);      Combine data and compute the density:               % Data     data = [wallD; uniformD]          % Density         uniDensity = 1 / ((upper - lower) ^ dim);          wallDensity = 1;     for i=1 : dim        wallDensity = wallDensity/(upperWall(i)-lowerWall(i));      end          wallSpace = (data(:,3) < upperWall(3)) & (data(:,3) > lowerWall(3));          trueValues = wallP * wallDensity .* wallSpace  + ...         uniP  .* (ones(N, 1) * uniDensity);      The wallSpace is a boolean array that indicates for each observation in data whether or not it lies within the wall. Since the range of the wall is equal to the range of the wall is equal to the range of the uniform data in dimension one and two I only consider the third dimension. If a point with index `i` isn't part of the wall its density is `uniDensity`, since `wallSpace(i)` is zero for such walls `trueValues(i)` equals uniDensity. A point with index `j` whose z is between 14.5 and 15.5 is in the wall, its density should thus be $\frac{4}{5} \cdot$ `wallDensity` \+ $\frac{1}{5} \cdot$ `uniDensity`. Since `wallSpace[i]` is one for these points, this is the density that is placed in `trueValues[j]`. `blockUniformDist(lowerWall, upperWall, wallN, dim)` is defined as:               function [ data ] = blockUniformDist( lower, upper, N, dim )     %BLOCKUNIFORMDIST Samples N values from a uniform distribution with      % dim dimensions.     % INPUT     %   - lower:   The lowest value allowed (per dimension)     %   - upper:   The highest value allowed (per dimension)     %   - N:    Number of samples to be taken     %   - dim:  Dimension of the distribution     % OUTPUT     %   - data: A vector of samples form the distribution          % values     data = rand(N, dim);          for i= 1 : dim             data(:,i) = lower(i) + data(:,i).*(upper(i) - lower(i));         end      end      And `uniformDist` as:               function [ data ] = uniformDist( lower, upper, N, dim )         %UNIFORMDIST Samples N values from a normal distribution with mean mu         %and standard deviation sd in dimension dim.         % INPUT         %   - lower:   The lowest value allowed         %   - upper:   The highest value allowed         %   - N:    Number of samples to be taken         %   - dim:  Dimension of the distribution         % OUTPUT         %   - data: A vector of samples form the distribution              % values         data = lower + rand(N, dim) .* (upper - lower);     end      ## The result The result of this is that each observation has one of two densities either $8.962962963e-04$ or $7.407407407e-06$. Plotting the data set with the density dictating the colour (points with density $7.407407407e-06$ in red and points with density $8.962962963e-04$ in blue) results in: ![The dataset plotted, with all points with density $7.407407407e-06$ in red, and points with density $8.962962963e-04$ in blue.](http://i.stack.imgur.com/5I0RP.jpg) ## The actual question Shouldn't the points in the denser area of the plot all have the same density, and thus all the same colour?


title: Predictive Model from Counts Data

text: I have some data that is the number of times a person visited a doctors office over a course of $5$ years. I want to create a model that would be able to predict the most likely number of counts that a person might visit the doctors office. When I looked at the histogram of the counts distribution I noticed that it was skewed right so I am thinking about going with a poisson model. My questions are:   1. Is this a good model to predict the number of counts?   2. What other types of models/methods would it be good to look into?   3. Once I have my model how can I test the accuracy? Any help, comments or references would be greatly appreciated! Thank you!


title: feature selection for longitudinal data

text: I have a longitudinal data which looks like this. Number of time points are different for each ID. Y is the binary response variable (take values 0 & 1) and X1-X20 are either continuous or categorical features. ID Timepoint X1 X2 ...X20 Y 1 09/19/2009 1 09/27/2010 2 02/05/2009 2 02/27/2010 2 03/01/2011 There are missing values for X1-X20 indicated by NA. I want to apply Random Forest method for feature selection (X1-X20). Not sure how to proceed for such longitudinal data since the samples are not independent. Will appreciate any suggestions.


title: Asymptotic distribution of uniform order statistics

text: It can be shown that for an iid sample from a Uniform(0, 1) distribution, \begin{equation} n(1-U_{(n)}) \rightarrow exp(1) \\\ n(U_{(1)}) \rightarrow exp(1) \end{equation} To see this just try finding the distribution function of the left hand side and then take the limit to infinity. Now, it turns out we can show that they are actually asymptotically independent and it's possible to derive the joint asymptotic distribution of them. My question is how can we do this. I found a book that states this but omits the proof. I'm hoping it wasn't omitted because of it's simplicity because I can't see the answer. **Edit:** Here I attempt to be fancy with characteristic functions and use a trick similar to something I've seen in a book. However, something is still wrong in this calculation but I can't figure it out. Let \begin{equation} \Phi(t_1,t_2) = \mathbb{E}\exp{(it_1 n^2 U_{(1)} + i t_2 n^2 (1-U_{(n)})}) \end{equation} Then notice the characteristic function we desire will be $\Phi(\frac{t_1}{n},\frac{t_2}{n})$. Now we expand this about 0 by differentiating under the expectation (Taylor expansion). Noting that the Uniform order statistics are actually beta distributed we can get their expectations from the formula for the beta mean. We get \begin{equation} \Phi(\frac{t_1}{n},\frac{t_2}{n}) = 1 + (i \frac{n^2}{n+1} , i \frac{n^2}{n+1}) (\frac{t_1}{n}, \frac{t_2}{n})^T + o(\|(\frac{t_1}{n}, \frac{t_2}{n})\|) \end{equation} The problem is that this goes to $1 + it_1 + it_2$ which is wrong because it should factor into the characteristic functions of two independant exponentials. Can someone figure out where I went wrong in the calculation above?


title: Survey Weighted Regression Without FPC in R

text: I'm using the svydesign package in R to run survey weighted logit regressions as follows:               sdobj <- svydesign(id = ~0, weights = ~chweight, strata = ~strata, data = svdat)          model1 <- svyglm(formula=formula1,design=sdobj,family = quasibinomial)      However, the documentation states a caveat about regressions without specifying finite population corrections (FPC):                   If fpc is not specified then sampling is assumed to be     with replacement at the top level and only the first stage of      cluster is used in computing variances.      Unfortunately, I do not have sufficient information to specify my populations at each level (of which I sampling very little). Any information on how to specify survey weights without FPC information would be very helpful.


title: Bayesian network inference using pymc (Beginner's confusion)

text: I am currently taking the PGM course by Daphne Koller on Coursera. In that, we generally model a Bayesian Network as a cause and effect directed graph of the variables which are part of the observed data. But on PyMC tutorials and examples I generally see that it not quite modeled in the same way as the PGM or atleast I am confused. In PyMC the parents of any observed real world variable are often the parameters of the distribution that you use to model the variable. Now my question really is a practical one. Suppose I have 3 variables for which data is observed (A, B, C) (lets assume they are all continuous variables just for the sake of it). From some domain knowledge, one can say that A and B cause C. So we have a BN here - A, B are the parents and C is the children. now from the BN equation P(A, B, C) = P(C | A, B) * P(A) * P(B) I can say A and B are some normal distributions with some mu and sigma, but how do I model P(C | A, B) ? The general idea I want to learn, is how do I learn this BN using PyMC so that I can query the BN. Or do I have to augment the BN with parameters of the model in some fashion. Is this problem solvable using pymc? or have I got some fundamentals wrong? Any help would be appreciated!


title: Generating Vectors under Constraints on 1 and 2 norm

text: Update: I left out some important information in my previous description... I am actually dealing with a special problem, which is better described as follows: Given user-specified parameters $\alpha$ and $\beta$ (where $0 < \alpha < \frac{1}{n+1}$, $\beta > 0$), generate at least one $n$-vector $A=[a_1,a_2,....,a_N]$ such that: $[ (a_1)^2 + (a_2)^2 + .... + (a_N)^2 ] = (\beta - \alpha)*\alpha$, where $\alpha = 1 - (a_1 + a_2 + .... + a_N)$ and each component of $A$ is restricted to lie in the interval $(\alpha,1)$. *Note that I am using rounded parentheses to indicate that the interval does not include its boundary points. * * * While working on a broader problem in developing a computer program, I realized that I could save a lot of time if I could generate n-vectors (for fixed n), subject to the following constraints:   1. the 1-norm of each vector must be equal to some user-specified scalar $\lambda$ (note: $\lambda$ will always be between 0 and 1)   2. the square of the 2-norm of each vector must be equal to some user-specified scalar $\phi$ (note: $\phi$ will always be positive)   3. Each vector component must be between 0 and 1 Of course, there are also logical constraints on the choices of $\lambda$ and $\phi$ given the above constraints... but I suppose we can assume that these parameters are well chosen by the user. Ultimately, I am searching for a way to accomplish this (or come as close as possible to accomplishing this). I'm not sure how to 'get off the ground', however.


title: Moment generating function of multinomial distribution

text: How would one find the moment generating function of the multinomial distribution, $\underline{X} \sim \mathrm{multinomial}(n, \underline{p})$? I know that by definition we have $$M_X (\underline{\theta}) = \mathbb{E} \exp{(\underline{\theta}^T X)} = \mathbb{E} \exp{\sum_{i=1}^k \theta_i X_i }$$ but then I can't see how we can go from here? I can't really get my head around what to do with the expected value of a vector. I don't think the $X_i$'s are independend so we can't factor out the expected value. Any explanation would be appreciated!


title: PCA Model Selection Using AIC

text: I want to use the Akaike Information Criterion (AIC) to choose the appropriate number of factors to extract in a PCA. The only issue is that I'm not sure how to determine the number of parameters. Consider a $T\times N$ matrix $X$, where $N$ represents the number of variables and $T$ the number of observations, such that $X\sim N\left(0,\Sigma\right)$. Since the covariance matrix is symmetric, then a maximum likelihood estimate of $\Sigma$ could set the number of parameters in the AIC equal to $\frac{N\left(N+1\right)}{2}$. Alternately, in a PCA, you could extract the first $f$ eigenvectors and eigenvalues of $\Sigma$, call them $\beta_{f}$ and $\Lambda_{f}$ and then calculate $$\Sigma=\beta_{f}\Lambda_{f}\beta_{f}'+I\sigma_{r}^{2}$$ where $\sigma_{r}^{2}$ is the average residual variance. By my count, if you have $f$ factors, then you would $f$ parameters in $\Lambda_{f}$, $Nf$ parameters in $\beta_{f}$, and $1$ parameter in $\sigma_{r}^{2}$. Is this approach correct? It seems like it would lead to more parameters than the maximum likelihood approach as the number of factors increases to $N$.


title: Outlier detection using Medial Rule

text: I have a data set captured every minute. Using Tukey's test for outlier detection (i.e., if number points above Q2+2.3IQR is above 10% of total distribution) then I mark the one minute window as anomalous. It works fine for the first minute. In reality, the system contines the behavior for the next few minutes and the algorithm will not mark the next few windows as anomalies as the median is shifted way above and continuous. By domain knowledge, the system continues in anomalous state for next two to three minutes and then settles down. To my understanding, I need to add one more algorithm along with Tukey's test so that the shifts are also measured and marked as anomalies Any help on finding an algorithm is much appreciated. PS: I couldn't attach as I don't have enough reputation. Let me shed more light... My problem is something different from masking. Window 1: Distribution: Right tailed normal distribution Outlier based on Median Rule: Yes, as certain number of points lies above Q2+2.3IQR Median of the distribution = 23 Windows 2: Distribution: Normal distribution but little large spread, as the effect of window 1 continuous in window 2 Outlier based on median rule: No, the spread (variance) of distribution is more and there are no point lies outside Q2+2.3IQR Median of the distribution = 40 Windows 2: Distribution: Normal distribution but little large spread, as the effect of window 1 continuous in window 2 and 3 Outlier based on median rule: No, the spread (variance) of distribution is more and there are no point lies outside Q2+2.3IQR Median of the distribution = 60 Suddenly the the behavior drops to something like it was in window 1, From the domain knowledge the window 2, window 3 are anomalies. On the other hand we can argue that this cannot be treated as anomalies as the system behavior repeats such a way in learning phase


title: Why does sphericity diagnosed by Bartlett's Test mean a PCA is inappropriate?

text: I understand that Bartlett's Test is concerned with determining if your samples are from populations with equal variances. If the samples are from populations with equal variances, then we fail to reject the null hypothesis of the test, and therefore a principal component analysis is inappropriate. I'm not sure where the problem with this situation (having a homoskedastic data set) lies. What is the problem with having a data set where the underlying distribution of all your data is the same? I just don't see the big deal if this condition exists. Why would this make a PCA inappropriate? I can't seem to find any good information anywhere online. Does anyone have any experience with interpreting why this test is relevant to a PCA?


title: Transformation of dependent variable: how to interpret it?

text: I recently reread some statistics books and noted something weird: They all discuss the assumptions of linear regression and mention the need for a normal distributed dependent variable. In the next step, there is always a discussion of variables which are non-normal and potential "cures" like sqrt, log, inverse etc. But the books are never going to discuss how to interpret the transformed variable - which strikes me as odd, since a sqrt function changes the interpretation from "the GDP of a state increases by 1 when the social-freedom index increases by 1" to "the square root of GDP of a state increases by 1 when the social-freedom index increases by 1", which is quite a difference. Can somebody point me to some good explanation about this?


title: are there any nonparametric forecasting methods?

text: Are there any good statistical non-parametric forecasting methods besides machine learing methods like neural netwworks/decision trees etc. for time series analysis ? If so, are there any R packages that implement them ? Thanks


title: Relation between confidence interval and testing statistical hypothesis for t-test

text: It is well known that confidence intervals and testing statistical hypothesis are strongly related. My questions is focused on comparison of means for two groups based on a numerical variable. Let's assume that such hypothesis is tested using t-test. On the other side, one can compute confidence intervals for means of both groups. Is there any relation between overlapping of confidence intervals and the rejection of the null hypothesis that means are equal (in favor of the alternative that means differ - two-sided test)? For example, a test could reject the null hypothesis iff the confidence intervals do not overlap.


title: Logistic regression and maximum entropy

text: I have read (e.g. here) that a (multinomial) logistic regressor corresponds to a maximum entropy classifier. My question is, how does one end up with the formula for logistic regression starting with the maximum entropy principle?


title: What are the assumptions of ARIMA modeling for forecasting time series?

text: What are the assumptions of ARIMA / Box-Jenkins modeling for forecasting time series?


title: Formulas for probabilities in Bayes theorem

text: In continuation to this question $p(h|D) = \frac{p(D|h)p(h)}{p(D)}$ $p(h) = $prior probability of hypothesis $h$ $p(D)$ = prior probability of training data $D$ $p(h|D)$ = probability of $h$ given $D$ $p(D|h)$ = probability of $D$ given $h$ I am inputting to a machine a training data $D$ in form of a table with $n$ training instances. The machine read only $m$ training examples among $n$ training examples in order to get a hypothesis as output. I have the following understanding: . $p(h) = $prior probability of hypothesis $h$ = $\frac{\left\vert{Version- space}\right\vert}{\left\vert{Hypothesis-space}\right\vert}$ $p(D)$ = prior probability of training data $D$ = $\frac{\left\vert{Observed- training-data}\right\vert}{\left\vert{Training-data}\right\vert}$ = $\frac{{m}}{{n}}$ $p(D|h)$ = probability of $D$ given $h$ = $\frac{\left\vert{Training-examples- satisfied-by-hypothesis-h}\right\vert}{\left\vert{Training data}\right\vert}$ = $\frac{\left\vert{Training-examples-satisfied-by- hypothesis-h}\right\vert}{{n}}$ I do not know how to interpret $p(h|D)$ = probability of $h$ given $D$ in the above forms. Is my understanding about formulas for corresponding probabilities is correct and what is formula for $p(h|D)$?


title: Generalized linear models with continuous proportions

text: We are having problems running generalized linear models with proportional data. For example, we have data like this:               Species  Trait(Diet)   IndividualsinForest   TotalIndividuals   ProportionForest      X        Insectivore   300.5                 500.7              0.60                  Y        Frugivore      32.3                  47.6              0.67                    And we want to determine whether trait influences the proportion of individuals in forest. Note that the individual numbers are continuous (have decimals), because the original counts have been DISTANCE-adjusted. Most on the models that we have seen in R deal with count data. For example, see M.J. Crawley _Statistics: An Introduction using R_ , Chapter 10: Analyzing proportion data (pdf). These models use the `cbind` statement to bind together the number of successes and the number of failures and then relate that to the predictor:               model.1 <- glm (cbind(IndividualsinForest, (TotalIndividuals-IndividualsinForest))~Diet,                      family=binomial)      However, there seems to be two problems to using this approach for our situation:   1. We are using `glm` on a `cbind` that combined vectors with continuous (decimal point) values; it is not clear whether this is allowed, although it runs in R. We have played around with multiplying everything by 1000 to get rid of the decimals, but find the results change a lot with such a technique (very different results if multiply by 10, by 100, by 1000)   2. `cbind` runs a weighted regression, as explained by Crawley. This means that the model will weight the data of Species X, above, with 500.7 total individuals more heavily than the data of Species Y, with 47.6 individuals. But conceptually, we do not want a weighted regression, because we want to treat all species equally: a rare species’ data is as important to us as an abundant species’ data. Species is the unit of replication here. We tried to force a non-weighted regression using the same code, by adding columns to the data, for example,               Species   ProportionForest   PInForest   PTotal      X         .60                60          100      Y         .67                67          100       And then running:               model.3 <- glm(cbind(PInForest, (PTotal-PInForest))~Diet, family=binomial)      But the resulting analyses are very overdispersed, and running them with `family=quasibinomial` gives strange results (completely non-sensitive tests, with large p-values for fairly clear differences). Hence our current analytical strategy is to fall back on arcsine transforming the proportion and then running a general linear model. But we'd prefer though to run a generalized linear model because the arcsine transformation seems to be a relict of old pre-computer statistics. Any ideas?


title: Distribution of the hyperbolic distance between random points in the Poincaré disc

text: Let two points at polar coordinates $(r_i, \theta_i)$ and $(r_j, \theta_j)$ be hyperbolic points in $\mathbb{H}^2_\zeta$ with curvature $K=-\zeta^2$. The radial coordinates of these points are exponentially distributed with exponent $\alpha$: $p(r) = \alpha e^{\alpha(r - R)}$ with $R \approx \ln(N)$ (the radius of the Poincaré disc), $N$ being the number of points in this disc, exponent $\alpha > 0$ and $r \in [0, R]$. The angular coordinates are uniformly distributed: $p(\theta) = \frac{1}{2\pi}$ with $\theta \in [0, 2\pi]$. For simplicity, let $\zeta = 1$. As a result, the hyperbolic distance between points in $\mathbb{H}^2_1$ is defined as follows: $d_{i, j} = \text{acosh}[\cosh(r_i)\cosh(r_j) - \sinh(r_i)\sinh(r_j)\cos(\theta_{i, j})]$ with $\theta_{i, j} = \pi - |\pi - |\theta_i - \theta_j||$. What's the distribution of distances between random points with the above characteristics in the hyperbolic disc of radius $R$? **Edit** I've been playing with the generation of random points using $p(r)$ and $p(\theta)$ for different values of $\alpha, N=1000$ and $\zeta=1$. These are the resulting point densities and hyperbolic distance distributions: ![Hyperbolic Densities](http://i.stack.imgur.com/RUotZ.png) By eye, the hyperbolic distances seem to be $Beta$ distributed. Of course a proper analysis has to be carried out to confirm this hypothesis...


title: Differentiation and dimensionality reduction

text: I have a matrix $x=N\times M$ of $N$ data points, where each one has $M$ features. Also, $y$ is the binary labels vector. In my case, $N$ is much smaller than $M$, so before running a classifier like SVM I used PCA to reduce the dimensionality. I picked the first 100 components which explained 85% of the variance. SVM with those 100 features gave me a poor discrimination of AUC ~0.6. Recently I thought on another method for dimensionality reduction, which I'm not sure about its validity: Each feature is a vector of $N$ values, which can be separated into 2 different vectors according to the labels. So for each feature, I can calculate how well it can differentiate the two classes by its own, using AUC, or even just the p-value of a t-test (given the histograms are normal). The problem with this method, is the dependency of the features vectors. Does anyone know how can I pick the best $k$ features vectors, which together differentiate the classes in the best way?


title: filtering discrete data

text: I looking for good ways to estimate a true state of the system from a set of observations, that take several integer values. E.g. I have an sequence of observations X_i of the unknown variable Y. Both the observations and the true value of Y belong to a set {1,2,3,4,5}. Are there any specific techniques for this situations beyond moving averages? E.g. is there a Kalman filter set up for this setting, or anything better?


title: Analyzing categorical data where the last category is unbounded

text: This is a bit of a long-winded data vis question. We're gathering data from a water quality test measuring fecal contamination. The test measures fluorescence in 11 discrete chambers. The water is known to be more contaminated as more of the chambers fluoresce and there's a non- linear way to quantify the water's contamination based on the chamber count. That is, one chamber corresponds to a contamination level of 2, five chambers is 46, ten chambers is 210 and all eleven chambers is.. an unknown level greater than 210. I'm having a hard time figuring out how to present that last case. I can plot contamination levels but how should I visualize or analyze "greater than 210?" Is there a term for this? Thanks!


title: SVM options in scikit-learn

text: Just curious about two options in the scikits SVM class. Does anyone know what Scale C and shrinking do? Not much in the documentation unfortunately. Thanks


> **Possible Duplicate:**   >  How to do a pretty scatter plot in R? How to use ggplot to do the binned quantile plots(one type of scatter plot)? Hi all, I have done scatter plot: plot(x, y). Now I wanted to do binned quantile plots... can ggplot2 help me? For example, we bin x data into 10 bins. For each bin, we draw the 10 deciles of the corresponding y data in that bin as points/dots. And then accross all bins, we would like to connect the corresponding decile points/dots together(something like equi-quantile or equi-decile curves)... How do I do that in R or ggplot2? Is there an existing function/command that can do this? Thanks a lot!


title: A multitrait-multimethod matrix and data set

text: I am working my way through Multitrait-Multimethod Matrix in a psychometrics class. We're only required to be able to analyze them but I'd really like to be able to construct them. I think I am able to do that (basically it's a rearangement of the correlation matrix with coefficient alpha for the diagonals). I am pretty sure I got it but want a **data set** and a **finished MTMM example** (something that looks like the format below):![enter image description here](http://i.stack.imgur.com/ckBrr.png) That way i can work my way through and make sure my final product looks like it's supposed to. I've attempted for a number of hours to find something but am frustrated and hoping someone has such a resource. I use R so i don't care what format the data set is in. Thank you in advance for your help.


title: The approach of labelling a collection of documents using a shared topic

text: I have a collection of documents, and know they may share a single topic. Is there a way to identify this topic? I know LDA (Latent Dirichlet Allocation) is an approach. But LDA result is to associate each document with a set of topics. However, I am more interested to find a common topic that can be used to label this collection.


title: Probability of events

text: Here is a probability problem: you observe .5 cars on average passing in front of you every 5 minutes on a road. What is the probability of seeing at least 1 car in 10 minutes? I'm trying to solve this in 2 ways. The first way is to say: P(no car in 5 minutes) = 1 - .5 = .5. P(no car in first 5 minutes and no car in second 5 minutes) = P(no car in first 5 minutes) * P(no car in second 5 minutes) by independence. Therefore P(at least 1 car in 10 minutes) = 1 - .5*.5 = .75. However, if I try the same, with a Poisson distribution with rate lambda = .5 per unit of time, for 2 units of time, I get: P(at least 1 car in 2 units of time) = 1 - exp(-2*lambda) = .63. Am I doing something wrong? If not, what explains the discrepancy?


title: Multivariate normal posterior

text: This is a very simple question but I can't find the derivation anywhere on the internet or in a book. I would like to see the derivation of how one Bayesian updates a multivariate normal distribution. For example: imagine that $$ P({\bf x}|{\bf μ},{\bf Σ}) = N({\bf \mu}, {\bf \Sigma}) \\\ P({\bf \mu}) = N({\bf \mu_0}, {\bf \Sigma_0}) $$ After observing a set of ${\bf x_1 ... x_n}$, I would like to compute $P({\bf \mu | x_1 ... x_n})$. I know that the answer is $P({\bf \mu | x_1 ... x_n}) = N({\bf \mu_n}, {\bf \Sigma_n})$ where $$ \bf \mu_n = \Sigma_0 (\Sigma_0 + \frac{1}{n}\Sigma)^{-1}(\frac{1}{n}\sum_{i=1}^{n}{\bf x_i}) + \frac{1}{n}\Sigma(\Sigma_0+\frac{1}{n}\Sigma)^{-1}\mu_0 \\\ \bf \Sigma_n = \Sigma_0(\Sigma_0 + \frac{1}{n}\Sigma)^{-1}\frac{1}{n}\Sigma $$ I am looking for the derivation of this result with all the intermediate matrix algebra. Any help is much appreciated.


title: ROC vs. Accuracy

text: If you want to compare two learning algorithms, which metric is better to use in general: ROC or accuracy? I understand that in ROC, you get both the sensitivity and specificity?


title: Sequential testing

text: I have been exposed to R before but after being away from the language for 10 years I find my self having to start from the start. I have a large data set with over 1 million rows and 20 variables. I know there are quicker and faster ways of mining that data but my customer insists on using hypothesis testing to identify significant segments/clusters. I am looking for some R code that would be able to perform a sequential t-test across each of the 20 variables one at a time. Only returning significant variables and samples in table format. Would really appreciate the help.


title: How can I get residuals output in a variable or to act like a data frame in R?

text: I have an object:               noise.lm      it's just a simple linear model with an X and Y. when I type in               resid(noise.lm)      it produces the residuals like this:                        1           2           3           4           5           6      -0.40501681 -1.16960373 -0.66636219 -0.41798296  1.37942381  0.84336242       I would like to access the residual values without the observation numbers, maybe by making a data frame like this                  residual     1  -0.40501681     2  -1.16960373     3  -0.66636219     .       .     .       .     .       .      Or is there an easier way to store just the residual value in a variable like:               residual_variable          -0.40501681, -1.16960373, -0.66636219, -0.41798296,  1.37942381,  0.84336242      An example of what I would like to do would be something like testing for autocorrelation               plot(length(residual_variable),residual_variable)      or when trying to check normality with a histogram... Thanks in advance


title: Least angle regression for a set of vectors?

text: As far as I know, LARS solves the following problem (using the same notation as Efron et al. Least angle regression): > Given a vector y, and a matrix X. Pick some column vectors from X, and > express y as a linear combination of these column vectors. In my problem, I have a set of vectors, which can form a matrix Y. I would like to express each column vector of Y as a linear combination of some column vectors from X. I hope that the fitting quality could be good (L2), and the number of column vectors I use from X is small (L1). Therefore, I am wondering if there is any systematic study (any papers?) on applying LARS to fit a set of vectors, rather than just one vector. Below is my solution. Comments are highly appreciated. Suppose Y = (y1 y2 ...) and X = (x1 x2 ...), where yi, xi are vectors. Treat Y as a high-dimensional 1D vector (flatten it) like $Y' = \left(\begin{array}{c} y1\\\ y2\\\ ...\end{array}\right)$ Then we define the new basis matrix X' as (the blank part are all zeros) $X' = \left(\begin{array}{ccc} X & & \\\ & X &\\\ & & ...\end{array}\right)$ Then we just apply LARS directly to fit Y' using X'. Once we have picked columns from X', we merge them into a set of columns from X, which will be the final "basis vector" set to fit Y. One thing I feel uncomfortable about the above solution is that, it does not consider the fact that the same set of basis vectors from X will be used to express all column vectors of Y. Only the final step merges the result from LARS. Your help is appreciated. Thanks.


title: How to approximate measurement uncertainty?

text: At the moment I use standard deviation of the mean to estimate uncertainty: ![definition = stdev over sq.root of N](http://upload.wikimedia.org/math/8/a/e/8ae657bf495238caf69d82d754e70ed8.png) where _N_ is in hundreds and mean is a time series (monthly) mean. I present it then like this: ![mean plus-minus stdev about the mean](http://rogercortesi.com/eqn/tempimagedir/eqn1898.png) for each element (month) in the (annual) time series. Is this valid? Is this appropriate for time series?


title: State-of-the-art in smoothing splines

text: What is the state-of-the-art in the efficient computation of smoothing splines? The algorithm I see mentioned most often is that of Reinsch, dating back to 1967. As I understand it, the most expensive part is performing a Cholesky decomposition. What other algorithms for smoothing splines have there been since then?


title: How are CP (Cost Complexity) values calculated in RPART (or decision trees in general)

text: From what I understand, the cp argument to the `rpart` function helps pre- prune the tree in the same way as the minsplit or minbucket arguments. What I don't understand is how CP values are computed. For example               df<-data.frame(x=c(1,2,3,3,3,4), y=as.factor(c(TRUE, TRUE, FALSE, TRUE, FALSE, FALSE)), method="class")     mytree<-rpart(y ~ x, data = df, minbucket = 1, minsplit=1)      Resulting tree...               mytree     n= 6           node), split, n, loss, yval, (yprob)           * denotes terminal node          1) root 6 3 FALSE (0.5000000 0.5000000)         2) x>=2.5 4 1 FALSE (0.7500000 0.2500000) *       3) x< 2.5 2 0 TRUE (0.0000000 1.0000000) *      Summary...               summary(mytree)          Call:     rpart(formula = y ~ x, data = df, minbucket = 1, minsplit = 1)       n= 6                    CP nsplit rel error    xerror      xstd     1 0.6666667      0 1.0000000 2.0000000 0.0000000     2 0.0100000      1 0.3333333 0.6666667 0.3849002      Where's the .666 and .01 coming from?


title: How to retrieve the prediction matrix according to the formula of a regression model in R?

text: Suppose I have a logistic regression model such like this:               set.seed(123)     df<-data.frame(     y=rbinom(100,1,0.5),     x1=rnorm(100,10,2),     x2=rbinom(100,20,0.6))          fit<-glm(y~x1*x2,data=df,family="binomial")     coef(summary(fit))                    Estimate Std. Error    z value  Pr(>|z|)     (Intercept)  5.08314564 6.43692399  0.7896855 0.4297115     x1          -0.66691041 0.64071095 -1.0408912 0.2979260     x2          -0.28338654 0.51254819 -0.5528974 0.5803337     x1:x2        0.04037126 0.05100223  0.7915588 0.4286180      Does somebody know how to get the prediction matrix in a format like this:                   intercept x1       x2  x1:x2     1   1        10.506637 10  105.06637     2   1        9.942906  17  169.02941     3   1        9.914259  10  99.14259     4   1        12.737205 11  140.10925


title: Question on Matrix Operation Notation

text: I have confused myself in attempt to express a 2-step operation in matrix notation-- The data consists of 2 "chunks" of information--quantity of identical items sold in different stores and their prices at each store. Suppose I have the data in the following format-- ![enter image description here](http://i.stack.imgur.com/o2Jgu.png) * * * Suppose I make the quantity and price information into separate matrices (so that the "quantity matrix" is a 3-by-3 matrix consisted of columns 2,3,4 data from the above matrix and the "price matrix" is a 3-by-3 matrix consisted of columns 5,6,7 data from the above matrix. Suppose I was given a "quote" vector for each store's budget estimation. The "quote" vector consists of 3 elements [quote for item A, quote for item B, quote for item C] How could I express the following operations in matrix notations? For example, for store 1, I want to multiply item A quantity by its price (call it "overall cost of A"), item B quantity by its price (overall cost of B) and item C quantity by its price (overall cost of C), then I want to calculate the sum of the overall cost of all 3 items. I want to do this operation for all the stores.


title: How to perform parameters tuning for machine learning?

text: I have a very basic question regarding parameter tuning using grid search. Typically some machine learning methods have parameters that need to be tuned using grid search. For example, in the following standard formulation of SVMs: $$ \min_w\frac{1}{2}\|w\|^2 + C\sum_{i=0}^{N}{\rm{hinge\ loss}}(x_i,y_i,w) $$ We have to tune the model parameter $C$. If I am given the _training_ and _test_ sets, then I would first split the _training_ set into `validation_train` and `validation_test`. Tune the parameters using the `validation_train` and `validation_test` and then using the best parameters, retrain on the complete training set and finally perform testing on the test set. So my question is which of the following is better on an average?   1. Would it be better to use the final model from validation, which was trained on `validation_train` for final testing, because the parameters were optimized on this training set?   2. Or would it better to use the entire training set and retrain the model with the best parameters from grid search? Although the parameters were not optimized for this set, we have more final training data in this case.


title: PCA and random forests

text: For a recent Kaggle competition, I (manually) defined 10 additional features for my training set, which would then be used to train a random forests classifier. I decided to run PCA on the dataset with the new features, to see how they compared to each other. I found that ~98% of the variance was carried by the first component (the first eigenvector). I then trained the classifier multiple times, adding one feature at a time, and used cross-validation and RMS error to compare the quality of the classification. I found that the classifications improved with each additional feature, and that the final result (with all 10 new features) was far better than the first run with (say) 2 features.   * Given that PCA claimed ~98% of the variance was in the first component of my dataset, why did the quality of the classifications improve so much?    * Would this hold true for other classifiers? RF scales across multiple cores, so it's much faster to train than (say) SVM.   * What if I had transformed the dataset into the "PCA" space, and run the classifier on the transformed space. How would my results change?


title: How to minimize residual sum of squares of an exponential fit?

text: I have the following data and would like to fit a negative exponential growth model to it:               Days <- c( 1,5,12,16,22,27,36,43)     Emissions <- c( 936.76, 1458.68, 1787.23, 1840.04, 1928.97, 1963.63, 1965.37, 1985.71)     plot(Days, Emissions)     fit <- nls(Emissions ~ a* (1-exp(-b*Days)), start = list(a = 2000, b = 0.55))     curve((y = 1882 * (1 - exp(-0.5108*x))), from = 0, to =45, add = T, col = "green", lwd = 4)      The code is working and a fitting line is plotted. However, the fit is visually not ideal, and the residual sum of squares seems to be quite huge (147073). How can we improve our fit? Does the data allow a better fit at all? We could not find a solution to this challenge on the net. Any direct help or linkage to other websites/posts is greatly appreciated.


title: How should I analyze data with repeated measures variables and integer between-subject factors?

text: I ran a study, where a single variable was measured repeatedly for each participant, under different conditions. Analyzing this data with SPSS has been rather straightforward, by using the General Linear Model -> Repeated Measures ANOVA, and specifying each measure in the within-subject factors. My problem arises when I try to include some post-experiment data I collected, like a test that all participants completed after the experiment, with integral scores ranging from 0-24. I want to include this to see how these scores relate to the within-subjects measures, and adding this data as a between-subject factor makes sense, but due to the nature of the test score variable I have up to 25 categories that make the data impossible to work with. It seems like this should be a common issue, and my Googling has helped me a bit, but most of the examples I have seen use categorical variables like gender for the between-subject factors. What do I do if the between-subject factor isn't categorical? I will be very grateful for any insight, comments and tips you have. Thanks!


title: How to calculate precision and recall when some of the test data remains unclassified

text: Consider a situation where we are running a classifier (the actual classification algorithm doesn't matter here), and the class labels are given based on a score. If score > 0, the data point is labeled A, if score < 0, the data point is labeled B. All the training data contains data points with positive or negative scores. However, in my test data, there are a few points that return score = 0. How should I measure the precision and recall in this scenario where some points can't be classified into any class? * * * [ _additional information from a comment below_ ] I faced this in a sentiment classification task. The normalized scores are in the range [-1,1], with 0 being the score for documents with no sentiment. It so happened that I had no neutral documents in my training data, but in the test data, some documents returned a score of 0.


title: chi-square test to check relationshio between season and win/loss

text: Given records for a team in 2 semesters Fall and Spring. We want to check if there is a relationship between the semester "Fall or Spring" and result "win or loss". Can chi-square test be used? Example of data               Fall , Spring      win , loss       win , loss     win , loss     win , win      win, win     loss, loss     loss, win     loss, loss      Many thanks


title: How to create a "male" variable in stata

text: My stata regression currently has the female variable, but I want to test whether gender matters on another variable or not. How would I go about creating the male variable? Or do I even need it? Will dropping the female variable work as well?


title: Which model should I use for Cox proportional hazards with paired data?

text: I am hoping someone can help me with which model (frailty, strata or cluster) I should use for my data. I have paired data so I need to take that into account when modelling the Cox PH and am unsure which model will give me a more accurate result. My study was looking at the time it took for a person to become calm after being subjected to a particular stimulus. Each person was subjected to two different stimuli, on separate days. They were randomly assigned which stimulus was first. I have modeled this with survival analysis (time-to-event) but I now need to take into account that the data is paired. Any help in regards to when you would use frailty, strata or cluster models would be great.


title: How to deal with a skewed class in binary classification having many features?

text: I am doing data analysis in the mobile ad targeting domain. I have around **18 features** and for a combination of these features, the result is either True or False **(1/0)** depending on whether the impression was clicked or not. The problem here is that the **output class is highly skewed**. Click though rate is around **0.4%**. (i.e value is 1 only 4 out 1000 times). I have a data set of **2 million rows** and I am using 90% as train set and 10% as test set. I have used **logistic regression** from **sckit-learn** package in **python**. Now after training my model I get all values for test set as 0. Please tell me what the problem could be and what should I do to solve it? P.S. : I have tried increasing my data set size and also reducing the number of features(even to just one feature). If I see the probability of each class(0/1) in the test set, I get around 0.002 - 0.005 for a 1. Thanks


title: Issues with calculating gradient descent operation

text: I have this issue when using gradient ascent. I have some synthetic data and after my first iteration the objective function decreases and from the second iteration it keeps on increasing. Is it possible for a convex function. I also have the learning rate very low 0.001. Any suggestions.


title: Prediction model on online game economy

text: I want to study the economy of an online game. In specific I want to examine if there is a possibility to create a prediction model. I would try to describe the whole concept and I am asking for papers or models that could fit in my occasion. I have made a research for papers, but I haven't found anything. Users have many different ways for incomes. Daily free money, selling, quests etc. A user can create items like potions and then sell them to other users through an auction. Items don't have a standard price and users define the price by biding on it. Also, items have a few properties of healing etc. I want to create a model to predict if the price of items will increase or decrease due to other changes in the game. For example, if I increase the daily free money to users, I will put more money on the economy and want to see if the prices will change somehow. Or if I will change the properties of the item how this would affect the whole economy. I have plenty of real historic data to use, but I haven't found anything related with this.


title: How to test whether there is a significant difference in mean squared error between two datasets?

text: I want to test whether Alzheimer's disease causes a change in brain aging compared to healthy patients. Therefore I have constructed a linear regression model of spectral parameters of brain recordings(the features, or dependent variables) versus age (age is the independent variable). Now I wish to fit the model on the healthy patients, and then use the coefficients to calculate the expected age of the Alzheimer's patients - comparing the mean squared error of the healthy dataset and the Alzheimer's dataset should help show whether or not there is a difference to the aging due to the disease. (i.e. if the model that works well for healthy patients fails miserably for Alzheimer's patients then there is probably a difference) I guess I will fit the linear regression model on 80% of the healthy patients (the training set), holding back 20% (as a test set) to calculate the MSE with. I would use cross-validation but then I will end up with as many different sets of coefficients as I have folds and how would I know which to fit the Alzheimer's patients with? The mean of the coefficients perhaps? An advantage to cross validation though is that I would have a mean and standard deviation of the MSE estimates from the sets of healthy patients so I could use that to determine the significance of deviation between the healthy and diseased MSE's, which is handy. I guess I could also sample subsets many times from the Alzheimer's patients and create a set of MSE estimates which I could then calculate the standard deviation and mean of to get some idea of the variance there as well so I have an idea about how sensitive it is to that particular dataset. (Should I do this with replacement i.e. bootstrapping? Why or why not?) Any advice is greatly appreciated.


title: Does k-NN with k=1 always implies overfitting?

text: I found somewhere such statement, but on the other hand in some sources I found, that it is ok. What about risk of overfitting while using 1-NN in binary classification problem where explanatory variables are TF-IDF values (cosine measure)?


title: Bayesian MMSE estimators from a transformation of the observations

text: Consider a random variable X whose value we want to estimate using a Bayesian MMSE estimator. Let $O_1(X)$ be a set of observations which depend on $X$ in some complex way (captured by $P(O_1|X)$) then the MMSE estimator is the conditional mean $\hat{X}_1=\mathbb{E}[X|O_1]$. Now consider another set of observation which is a possibly complex transformation of the first $O_2(O_1)$ and the corresponding MMSE estimator $\hat{X}_2=\mathbb{E}[X|O_2]$. Now it is obvious that if the mapping $O_2(O_1)$ is deterministic then $MSE[\hat{X}_1]\le MSE[\hat{X}_2]$ from the minimality of the MSE of $\hat{X}_1$. My question is if this is true also when the mapping $O_2(O_1)$ is probabilistic, that is defined by a conditional distribution $P(O_2|O_1)$. Intuitively, it should be as any other stochasticity in the mapping just seems to introduce additional noise as it does not depend on $X$. But I wonder if one can show this explicitly. EDITED: as the MMSE estimator is unbiased, indid $MSE[\hat{X}_i]=V[\hat{X}_i]$


title: Algorithm to evaluate whether you should buy now or wait

text: I'm currently in search of an algorithm that can determine whether or not it's time to buy something (an item, a stock, a service, etc.) given an history of prices (30, 50, 100, ...). My idea is something like this : > Given the history of prices, you should buy now because the price is likely > to be going up. Or, > Given the history of prices, the price is likely to drop even lower, hence > you should wait to buy. I've been Googling to find this, but i think I might not have entered the right keywords, because I couldn't find anything. That, or this does not exist because it's not reliable. Statistics are not my main strength, but algorithms are. If you have knowledge of a mathematical study on this topic, I'll take it too.


title: Why variances from nested random effects appear identical?

text: Could anyone shed some light on a slightly confusing result from a lmer4 model in R? I have a nested random effect, but the variance explained by each level seems identical, and I can't see why. I have a data set of bee observations in controlled arrays of snapdragons, with different relative frequencies of two colour morphs. There are five frequency treatments, and each is replicated three times, and each of these 15 experiments has a separate bumblebee hive. For each array I have observations on multiple bees, and I want to know how often they move between the same colour morph or switch. I'm modelling this as a binary mixed model in lme4, using a nested random effect to account for the variation due to different bees, different hives, and different replicates:               glmer(next_yellow ~ freqY + nflowers + height +                        (1 | bee/hive/replicate),                     data=ftrans,family=binomial)      The fixed part accounts for whether the frequency of yellow plants, flower number or flower height predict whether the next plant visited will be yellow, and is all clear. What's confusing me is the output for random effects:               Random effects:      Groups              Name        Variance Std.Dev.      rep:(hive:bee) (Intercept) 0.24964  0.49964       hive:bee       (Intercept) 0.24964  0.49964       bee            (Intercept) 0.24964  0.49964      Number of obs: 1652, groups: rep:(hive:bee), 286; hive:bee, 286; bee, 286      Inspecting the data suggests that there's tons of variation between bees, and very little for replicates, but they're still all estimated as 0.24964. Is there a statistical for this? Or is it likely I've mucked something up with the coding of variables or something (I've looked at this, and I can't see anything wrong, and everything is coded as factor levels). Any insight much appreciated!


title: Finding uncertainty in coefficients from polyfit in Matlab

text: I have a data set of x- and y-values, that I want make a linear fit on. Using `polyfit(x,y,1)` I get the coefficients `a` and `b` for a linear fit `ax = b` for this data, but I would also like to find the uncertainty or standard deviation for these coefficients. Does anyone know an easy way of doing this? My Google-fu only gave me this result, and seeing as the last answer in that thread is a correction to first answer, I don't know if I should trust any of the answers. The answer from that thread is:               [z,s]=polyfit(x,y,1);     ste = sqrt(diag(inv(s.R)*inv(s.R')).*s.normr.^2./s.df);


title: Recreate logistic regression equation from table of odds data

text: I'm reading the technical manual for a linking study between two assessments. It's pretty clear that the table is model output from a fitted logistic regression equation. Here's what pass odds look like on test 2 as a function of score on test 1 (RIT score): ![enter image description here](http://i.stack.imgur.com/8lwbx.png) It seems silly to use a lookup table that rounds to 5 when the model that made that table could give a better estimate. But how do I recreate that equation from this output? I have a good sense of how I would fit this model if I had the raw data, but I'm not sure what to do here. Not `glm(family=binomial)` because the data I have is are odds ratios, not pass / no pass (i.e., 1s and 0s), right? Here's the data:               PASS <- c(0, 0, 0, 0.01, 0.01, 0.01, 0.02, 0.04, 0.06, 0.1, 0.15, 0.23,      0.33, 0.45, 0.57, 0.69, 0.79, 0.86, 0.91, 0.94, 0.96, 0.98, 0.99,      0.99, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)          RIT <- c(120L, 125L, 130L, 135L, 140L, 145L, 150L, 155L, 160L, 165L,      170L, 175L, 180L, 185L, 190L, 195L, 200L, 205L, 210L, 215L, 220L,      225L, 230L, 235L, 240L, 245L, 250L, 255L, 260L, 265L, 270L, 275L,      280L, 285L, 290L, 295L, 300L)