title: Order statistics (e.g., minimum) of infinite collection of chi-square variates?

text: This is my first time here, so please let me know if I can clarify my question in any way (incl. formatting, tags, etc.). (And hopefully I can edit later!) I tried to find references, and tried to solve myself using induction, but failed at both. I'm trying to simplify a distribution that seems to reduce to an order statistic of a countably infinite set of independent $\chi^2$ random variables with different degrees of freedom; specifically, what is the distribution of the $m$th smallest value among independent $\chi^2_2,\chi^2_4,\chi^2_6,\chi^2_8,\ldots$? I would be interested in the special case $m=1$: what is the distribution of the minimum of (independent) $\chi^2_2,\chi^2_4,\chi^2_6,\ldots$? For the case of the minimum, I was able to write the cumulative distribution function (CDF) as an infinite product, but can't simplify it further. I used the fact that the CDF of $\chi^2_{2m}$ is $$F_{2m}(x)=\gamma(m,x/2)/\Gamma(m)=\gamma(m,x/2)/(m-1)!=1-e^{-x/2}\sum_{k=0}^{m-1}x^k/(2^k k!).$$ (With $m=1$, this confirms the second comment below about equivalence with an exponential distribution with expectation 2.) The CDF of the minimum can then be written as $$F_{min}(x) = 1-(1-F_2(x))(1-F_4(x))\ldots = 1-\prod_{m=1}^\infty (1-F_{2m}(x)) $$ $$= 1- \prod_{m=1}^\infty \left(e^{-x/2}\sum_{k=0}^{m-1}\frac{x^k}{2^k k!}\right).$$ The first term in the product is just $e^{-x/2}$, and the "last" term is $e^{-x/2}\sum_{k=0}^\infty x^k/(2^k k!)=1$. But I don't know how (if possible?) to simplify it from there. Or maybe a totally different approach is better. Another potentially helpful reminder: $\chi^2_2$ is the same as an exponential distribution with expectation 2, and $\chi^2_4$ is the sum of two such exponentials, etc. If anyone is curious, I am trying to simplify Theorem 1 in this paper for the case of regression on a constant ($x_i=1$ for all $i$). (I have $\chi^2$ instead of $\Gamma$ distributions since I have multiplied by $2\kappa$.) Thanks for any help!   Dave


title: Linear regression with slope constraint

text: I want to perform a very simple linear regression in `R`. The formula is as simple as $y = ax + b$. However I would like the slope ($a$) to be inside an interval, let's say, between 1.4 and 1.6. How can this be done?


title: what's the challenge in big data analysis?

text: Could you guys give me a review on the challenges and common techniques in dealing with big data problem? For example, what I know is that in when sample size is limited and dimension is high, the estimation of covariance would be an issue, what we did is to suppressing the big eigenvalues and increasing the small eigenvalues.


title: How can I predict values from new inputs of a linear model in R?

text: I've created a linear model in R: `mod = lm(train_y ~ train_x)`. I want to pass it a list of X's and get its predicted/estimateed/forecasted Y. I looked at `predict()`, but I think that is for something else, or I just don't know how to use it. I'm guessing by taking the coefficients of my model, I could manually plugin the test_x variables one-by-one, and get a predicted Y, but I'm guessing there is a more efficient way to do this.


title: How to test the effect of variable on spatial distribution of points in R

text: I'm trying to find statistical test (or procedure) which is able to **discover patterns in spatial distribution of points.** I sketch the problem by giving the example about position of bird territories and personality of birds. **1.** Load this example data ( _they are completely fabricated_ ). Data consists of geographical coordinates of territories and the level of aggression for each bird. **Example data:**               x_coor<-rep(c(2,2.5,3,3.5,4,17,17.5,18,18.5,19,19.5,20,20.5,21,21.5,22,17.5,     18,18.5,19,31.5,32,32.5,33,33.5,1,1.5,2,3,4,5,16,17,18,19,20,21,22,23,24,25,     26,16,16.5,17,18,19,20,27.5,28.5,29.5,30.5,31,31.5,32,33,34,2,3,4,5,6,7,8,9,     12,13,14,15,16,18,20,21,22,23,24,25,26,28,29,30,31,32,33,34,1,7,8,10,14,22,     26,28,33),c(5,5,5,3,2,1,2,3,3,3,2,3,2,2,1,1,1,4,4,3,2,3,3,8,4,8,1,4,4,4,4,3,     2,2,2,2,2,2,3,3,3,1,1,1,3,5,5,7,1,2,3,3,3,1,5,5,5,4,1,2,2,3,2,2,3,1,2,4,3,3,     2,3,1,5,2,4,1,2,3,2,1,1,1,1,1,1,1,1,1,1,1,1,1,2))     y_coor<-c(15,15.5,16,16.5,17,16,16.5,17,17.5,18,16,17,17.5,18,18.5,17.5,18,     18.5,18,18.5,22,22,22.5,21.5,22,22.5,21.5,22,22.5,21.5,22,22.5,21.5,22,21.5,      22,22.5,22,22.5,22,22.5,22.5,22.5,7,6,6.5,7,7.5,5.5,6,6.5,7,5,5.5,6,16,16.5,     15,15.5,16,14,14.5,15,11,11.5,12,12.5,13,13.5,14,14.5,11.5,12,12.5,13,10,11,     12,13,14,15,16,17,18,13,14,19,20,14,15,20,21,15.5,16.5,20,21,17,18,19,20,21,     22,23,21,24,20,24,20,23.5,20,23.5,21,23.5,21,23.5,22,23,24,22,23,24,22,23,24,     23,7,6,4,5,8,2,3,4,5,8.5,1,2,3,4,8,1,2,3,4,5,6,7,16,16,17,16,17,18,16,17,18,     10,14,15,17.5,9,10,11,12,13,8,9,10,16,17,8,9,10,14,15,6,8,23,25,11,9,23,13,     25,11,15,22,18,20,14,23,16,18,20,6,3,9,1,6,22,24,4,9,19,2,12,26,11,18,10,18,     26,12,2,4,6,8,26,10,19,2,4,6,8,16,20,26,14,22,25,9,20,12,7,20,6,18,1,2,8,26,     15,15,10,4,1,25)     aggr_low<-1:35     aggr_mid<-25:65     aggr_big<-50:75     aggr_max<-70:100     aggression<-c(sample(aggr_low,75,replace=TRUE),sample(aggr_mid,100,     replace=TRUE),sample(aggr_big,62,replace=TRUE),sample(aggr_max,10,replace=TRUE))     my.data <- data.frame(x_coor, y_coor, aggression)     rbPal <- colorRampPalette(c("blue","red"))     my.data$Col <- rbPal(10)[as.numeric(cut(my.data$aggression,breaks = 10))]     par(mfrow=c(1,2))     plot(x_coor,y_coor,pch = 20,col = my.data$Col)      **2.** As you can see in the plot, the territories of low aggressive birds ( _blue_ ) have tendency to cumulate together, whereas aggressive birds are solitary ( _red_ ). I spend a long time by searching for statistical method which could confirm such pattern (in this case clearly visible). My first try was to correlate distance matrices of aggression of birds and spatial position of territories. Then applying `mantel.test` for similarity of these two matrices. This however did not work because `mantel.test` is **unable to absorb the information about spatial isolation of clusters** of low aggressive territories. It expect, that if low aggressive birds are together, they will be together all in one big cluster (not in four different clusters as visible from the plot). So `mantel.test` is unable to confirm or disprove the hypothesis about the influence of aggression on spatial distribution of territories. **3.** I have notified that four cluster do not know about each other, thanks to long distances among them, so I try different approach. I made again the distance matrix of all territories...               geo.dist<-dist(cbind(x_coor, y_coor))     geo.dist<-as.matrix(geo.dist)      ...but I keep **only distances of nearest neighbours by setting distance threshold** (manually to be the 3 distance units)               geo.dist[geo.dist > 3] <- NA      This creates the _imaginary circles_ around each territory and I was able to count the number of neighbours for each one of them.               no.neighbours<-numeric(ncol(geo.dist))     for (i in 1:ncol(geo.dist))     no.neighbours[i]<-sum(!is.na(geo.dist[ ,i]))-1      And as you can see form this plot...               plot(aggression,no.neighbours,col = my.data$Col, pch=20)      There is **clear negative correlation** between number of neighbours and aggression.               cor.test(aggression, no.neighbours)      **Is it correct to use this steps to confirm (or falsify) this hypothesis? ("Agression of birds influence the position of territories") Does anybody know some more appropriate solution or the name of commonly used test for such problems?**


title: $n$-th order Markov model with $m$ conditional probabilities

text: # Initial Problem I am working on a problem that determines the completion of projects. There are eight stages for the project to be in before completion, with the ninth stage being completion. There is a certain number of days between each project ranging from [0, $\infty$), i.e. the project can go from Stage 1 to Stage 2 on the same day. I have data for about 10000 projects with 6000 being completed and 4000 not- yet-completed and are currently at one of the given stages. If, out of the 6000 completed projects 5000 of them took 0 days to go from Stage 1 to Stage 2, that should be (with $S_{12}$ the route between Stages 1 and 2) as $P(S_{12}|t=0)=\frac{5}{6}$. If 100 projects took 1 day, then that should be $P(S_{12}|t=1)=\frac{1}{60}$ and so on and so forth. # Method I think the best way to go about this is by using an absorbing Markov chain, with the ninth stage as the absorbing node. (Please excuse any errors in assumptions or usage of terms with Markov chains as my experience is limited to reading Wikipedia and posts here.) So the way I think the chain should look: Each stage has a line with an arrow pointing to the next stage (except the last) along with an arrow pointing from each stage to itself. # Issues   1. "The past doesn't matter." Even though that's a rule for Markov chains, we do care if it took, say, 90 days to go from Stage 1 to Stage 2. Thus, when we get to Stage 6 that may imply that this project tends to lag between stages. This being the case, should it be an 8th order Markov Chain? 9th order? (Not sure if the fact that it's absorbing changes things.)   2. Conditional Probabilities: Above, I mentioned using probabilities such as $P(S_{12}|t=0)=\frac{5}{6}$. Before getting more into the Markov chain aspect, does it make sense how I computed the conditional probabilities from the data? Assuming they are correct, my understanding of computing Markov chains is that one eventually gets to $xP^n$ and then runs that for $n=${however long they want the process}. My $P$ matrix is now a 9x9 matrix, but it seems like I should have a different $P$ matrix for each time interval.      * Meaning: For 0-day projects the first element is $\frac{1}{6}$ since there is a $\frac{1}{6}$ chance of the project staying in Stage 1, and the element one row below should be $\frac{5}{6}$ since that's the probability of moving to Stage 2 in 0 days. This would be in my $P_0$ matrix. My $P_1$ would have $\frac{59}{60}$ as the first element and $\frac{1}{60}$ as the second. Let's say I do this all the way up to $P_{730}$ and the last matrix would be $P_{>730}$ (to make things easier).     * Now that I have 732 probability matrices, what do I do with them? As mentioned above, if I wanted to do 10 days with a regular Markov matrix I would just make $n=10$ for the eigenvalue/vector decomposition of $xP^n$. But instead do I just do $x P_0 P_1 P_2 P_3 P_4 P_5 P_6 P_7 P_8 P_9$? Eventually we would like to get to a point that for any project we can predict what date it should be finished. Any help would be much appreciated.


title: Prove the equivalence of the following two formulas for Spearman correlation

text: From wikipedia, Spearman's rank correlation is calculated by converting variables $X_i$ and $Y_i$ into ranked variables $x_i$ and $y_i$, and then calculating Pearson's correlation between the ranked variables: ![Calculate Spearman via wikipedia](http://i.stack.imgur.com/YXr3K.png) However, the article goes on to state that if there are no ties amongst the variables $X_i$ and $Y_i$, the above formula is equivalent to ![second formula to calculate Spearman](http://i.stack.imgur.com/5StVI.png) where $d_i = y_i - x_i$, the difference in ranks. Can someone give a proof of this please? I don't have access to the textbooks referenced by the wikipedia article.


title: Correct approach using t-test and Wilcoxon rank sum test?

text: I want to compare 2 datasets. I have a biological structure A with its atom coordinates which are a few thousands. I also have from these atom coordinates the coordinates over time of about 20 nanoseconds. From this structure A I made about 70 simulations. After the simulation I calculated the distance between atoms within this structure and afterwards the average distance of these 70 simulations for each time step which means that I finally have the average distance of 20 nanoseconds of the structures A. Now I also have a structure B with which I made 80 simulations and performed the same procedure as described above. Finally I want to compare the average distance over time of structure A and structure B. 1) I make a Shapiro-Wilk or Shapiro-Francia test to make sure that the data has a normal distribution 2) If the data is normally distributed I use a t-test, otherwise a Wilcoxon rank sum test Is this a correct approach?


title: Normally distributed errors - Why not use the observed residual histogram?

text: In the context of the classical linear regression model (with all standard assumptions), we know that when the error term is normally distributed, least squares is minimum variance among all linear unbiased estimators and we know the exact statistical distribution of the t statistics. However, if the residual histogram suggests that the error term is something other than normal, then the t statistics will not have the claimed t distribution (nor will the F statistic be F distributed). Surely, it does not make sense to assume normality if the residuals suggest otherwise. Instead, if we want to know the theoretical distribution of the t statistic, why not assume that the errors are distributed in population based on what we observe from the residual histogram?


title: Dealing with multicollinearity by separating common variance

text: On this page, I read about multiple ways to deal with multicollinearity. I determined that the last method suggested may be best for my purposes. I copied it below: > Treat the common variance as a separate variable and decontaminate each > covariate by regressing them on the others and using the residuals. That is, > analyze the common variance as a separate variable. I think that I get this conceptually but do not know how to do it in SPSS or R. If someone could explain how to do it in SPSS (preferably) or R, I'd appreciate it!


title: What machine learning technique I need to apply?

text: I have a unix shell data. Each session starts with `<begin>` and `<end>`. Between `<begin>` and `<end>`, the info contains users command in the shell session. From this data, I need to figure out how many users are on the system with probability. For example, the is something like this:               <begin>     pwd     cd ~     <end>     <begin>     ls     vi     <end>     <begin>     cat     ls     pwd     <end>      I have put this data into this type of data frame: head(result,20)                 user  dat          1  elm          2 gzip          2  <1>          2   ls          2   ci          2   -l          2  <1>         2   ci         2   -l         2  <1>         2   ci         2   -l         2  <1>         2   ci         2   -l         2  <1>         2   ls         2   ls         2  <1>        2   ls      I need this dat to be like this:               user ls cd vi cat      1   5   0  0  0      2   3   2  5   6      etc, any ideas how I would do this?


title: What's the fundamental difference between these two regression models?

text: Suppose I have a bivariate responses with significant correlation. I am trying to compare the two ways to model these outcomes. One way is to model the difference between the two outcomes: $$(y_{i2}-y_{i1}=\beta_0+X'\beta)$$ Another way is to use `gls` or `gee` to model them: $$(y_{ij}=\beta_0+\text{time}+X'\beta)$$ Here is a foo example:               #create foo data frame          require(mvtnorm)     require(reshape)     set.seed(123456)     sigma <- matrix(c(4,2,2,3), ncol=2)     y <- rmvnorm(n=500, mean=c(1,2), sigma=sigma)     cor(y)     x1<-rnorm(500)     x2<-rbinom(500,1,0.4)     df.wide<-data.frame(id=seq(1,500,1),y1=y[,1],y2=y[,2],x1,x2)     df.long<-reshape(df.wide,idvar="id",varying=list(2:3),v.names="y",direction="long")     df.long<-df.long[order(df.long$id),]         df.wide$diff_y<-df.wide$y2-df.wide$y1               #regressions     fit1<-lm(diff_y~x1+x2,data=df.wide)     fit2<-lm(y~time+x1+x2,data=df.long)     fit3<-gls(y~time+x1+x2,data=df.long, correlation = corAR1(form = ~ 1 | time))      What's the fundamental difference between `fit1` and `fit2`? And between `fit2` and `fit3`, given they are so close on the $p$ values and estimates?


title: Similarity measure between two text documents

text: Let the first document involve the words $\{x_1,x_2,\ldots,x_n\}$ and the second one be composed of $\{y_1,y_2,\ldots,y_m\}$ where $n$ is not necessarily equal to $m$. I have a similarity measure that works for elements of the sets; $s(x_i,y_j)$ for all $i$ and $j$. This similarity may not be symmetric, that is, $s(x_i,y_j)\neq s(y_j,x_i)$. But I'm not familiar about the similarity of two sets of instances. I thought that the average of the best similarities is reasonable: $$s\left(\{x_1,x_2,\ldots,x_n\},\{y_1,y_2,\ldots,y_m\}\right) = \frac{1}{n}\sum_{i=1}^n \max_{j=1}^m s(x_i,y_j)$$ What are the main methods used to compute $s\left(\{x_1,x_2,\ldots,x_n\},\{y_1,y_2,\ldots,y_m\}\right)$ in such a case? Particularly the similairty of two text documents?


title: Goodness of fit of a SUR model

text: I know that McElroy R^2 is a measure of goodness of fit for Seemingly Unrelated Regressions (SUR models), but how can one judge that the estimated equations are well fit by using the McElroy R^2?


title: R scripts accompanying the book of "Analysis of Financial Time Series" by Ruey S. Tsay

text: Are there any R scripts that accompany the book of `Analysis of Financial Time Series, Third Edition ` I went to the book's webpage, and just found the data set.


title: Sum of squares proof where $N=n_{UC}+n_{TX}$

text: $TX$ is variable that indicates treatment status ($TX=1$ if the patient gets the new treatment, and $0$ otherwise, and $UC = 1 - TX$ indicates they got the standard treatment). Of $N$ patients, $n_{UC}$ get usual care and $n_{TX}$ get new treatment. If $n_{UC} + n_{TX} = N$, show that $\sum_i(TX_i-\overline{TX} )(TX_i-\overline{TX}) = n_{UC}\cdot n_{TX}/N$.


title: improving user preferences based on behavior

text: I have a project to improve user experience in our software by gathering data on their actions while using the software and comparing that to preferences they made when initializing the software. The idea is that if they say they like Coke but they keep drinking Pepsi, we would (1) ask "Are you sure you like Coke, not Pepsi?", (2) just change their user settings to Pepsi, or (3) incorporate their usage history into whatever logic already uses their initialization preferences. Not all user preferences are exclusively categorical. Some may be continuous (30% Coke 70% Pepsi). Now, I could probably solve this programmatically but it seems like a statistical or machine learning approach might be better. Can someone tell me if this is a particular class of problem and/or point me to any literature that attempts to solve this from a statistical (or ML) point of view? Thank you


title: kernelized l1 norm and the representer theorem

text: I'm trying to derive a kernel-ized $l_1$ penalty for logistic regression. I have been looking at the slides Learning with Sparsity Inducing Norms along with the slides on Regularization and Variable Selection via the Elastic Net. On page 20 of reference 2 they add a $l_1$ penalty of the form $\lambda_1\sum_{i=1}^n|\alpha_i|$ with little explanation. Given we want a $l_1$ norm of the form $\sum_{i=1}^D|w_i|$, I think we can use the representer theorem to help things along. On page 36 of 1, "Regularization and representer theorem", they say that for the $l_2$ norm the minimized function w.r.t weights $w$ is of the form $w = \sum_{j=1}^n\alpha_j\Phi(x_j)$ and then say that this is also equal to $\sum_{j=1}^n\alpha_jk(.,x_j)$. I tried looking at the reference for this equality but found it totally incomprehensible. As best as I can figure, the $l_1$ norm should be something like $||\tilde{w}||_1 = \sum_i|\tilde{w}_i| \le |\sum_i\sum_{j=1}^N\alpha_j\Phi(x_j)|^{(i)}$ $|\sum_i\sum_{j=1}^N\alpha_j\Phi(x_j)|^{(i)} = |\sum_i\sum_{j=1}^N\alpha_jk(.,x_j)|^{(i)} = |\alpha^TK1|$ where $1$ is the column vector of ones. My question is two-fold:   1. Why is $f = \sum_{j=1}^n\alpha_j\Phi(x_j) = \sum_{j=1}^n\alpha_jk(.,x_j)$?   2. What should the $l_1$ norm look like?


title: Deep learning libraries/software with good tutorials or examples

text: I am looking for deep learning libraries that have very good tutorials or examples, since I would like to learn by examples or doing tutorials even if these are not the best deep learning libraries. If they also run on GPUs much better.


title: histogram normalised to area 1

text: i have a histogram with the y-axis showing the proportion in percentage. That makes sense to me but now i have read that histograms can be normalized with the result that the area of the rectangles is 1. Can someone explain me the idea,reason behind this? Many thanks!


title: Is a beta coefficient of 0.3441 strong?

text: I am wondering if there is a table or some reference that explains the strength of the beta coefficient.


title: Time series with an order of integration greater than 1

text: The issue of unit roots and what $ I(1) $ variables do to a regression estimation is covered clearly online. However, it is of course possible to have order of integration greater than 1. What does it actually mean to say a time series is $ I(n) $ where $n>1$ ? What would a series like this look like? Are there any examples of series we would expect to display this property? Does it present any further problems than an $ I(1) $ variable?


title: Really a knifes edge?

text: Something nice and topical. I just read these two items on the news:   * Obama is in the lead by 50.4% to 48%, with 61% of votes counted. (Ohio)   * With 86% of the vote counted, Virginia is still sitting on a knife edge. Romney is hanging on to a lead of 49.9%, but Obama is snapping at his heels on 48.7%. Intuitively they sound like "dead certs" not "knife edge", because the populations involved must be huge. Then I realized I didn't know exactly what to do to prove my thesis. Formalizing my thesis: I'm 99% sure that the current leader won't change once all votes are counted. But the question I have is: how many voters do there have to be in each state to be able to claim this? NOTE: to be clear on definition, if there are N voters in Ohio, and Obama is leading 50.4:48 based on counting 0.61 x N votes, what is N to have 99% confidence that the Obama ratio will be >=50? **UPDATE:** The comments have explained to me that the 86% was of precincts, not voters, with precincts supposed to represent around 2500 voters. But more importantly there may be an urban bias in the precincts that get counted last. If you wish to post an answer please assume my original assumptions: that it is 86% of actual voters, and that there is no bias in the votes counted so far. P.S. Apparently, Virginia ended up with Obama 1,868,191 (50.57%) to Romney 1,767,692 (47.85%).


title: Calculating the K-fold cross validation confusion matrix

text: I know this is a dumb question, but I can not find an answer (I guess everyone else knows this). I want to verify: When I do k-fold cross validation of a decision tree, the software produces a confusion matrix for the test data. Is each instance assigned the most frequent class in the cross validation results to create the confusion matrix?


title: Discrete variable, dummies with only one occurance; how do work with this in a regression?

text: I am roughly running an OLS regression of prices and the number of network providers in cities in the UK. Let $P$ stand for prices and $net$ stand for the number of network providers. The frequency of $net$ is:               net   frequency     0        100     1        140     2        100     3        40     4        10     5        10     6        1     7        1      After generating dummies d0, ..., d7, for each state of $net$ (for example, d0 is equal to 0 when there are at least one company on the market and zero otherwise), I ran a regression of $P$ over d0,...,d6. My question is not about this model per se, but about if it is correct in general to keep the dummies d6 and d7 given that they are equal to 1 only once. Should I group d5, d6 and d7? (that is consider d5+, a dummy for the case where there are 5 or more companies on the market), or the correct thing to do is to keep the given dummies as it is? On the general model I am considering, the outcomes change a bit when I work with d5+ instead of d5,d6,d7. Would you know a reference dealing with this situation?


title: How do we create a confidence interval for the parameter of a permutation test?

text: Permutation tests are significance tests based on permutation resamples drawn at random from the original data. Permutation resamples are drawn without replacement, in contrast to bootstrap samples, which are drawn with replacement. Here is an example I did in R of a simple permutation test. (Your comments are welcome) Permutation tests have great advantages. They do not require specific population shapes such as normality. They apply to a variety of statistics, not just to statistics that have a simple distribution under the null hypothesis. They can give very accurate p-values, regardless of the shape and size of the population (if enough permutations are used). I have also read that it is often useful to give a confidence interval along with a test, which is created using bootstrap resampling rather than permutation resampling. Could you explain (or just give the R code) how a confidence interval is constructed (i.e. for the difference between the means of the two samples in the above example) ? **EDIT** After some googling I found this interesting reading.


title: Eliciting priors from experts

text: How should I elicit prior distributions from experts when fitting a Bayesian model?


title: MANOVA and correlations between dependent variables: how strong is too strong?

text: The dependent variables in a MANOVA should not be "too strongly correlated." But how strong a correlation is too strong? It would be interesting to get people's opinions on this issue. For instance, would you proceed with MANOVA in the following situations? a) Y1 and Y2 are correlated with r=0.3 and p<0.005 b) Y1 and Y2 are correlated with r=0.7 and p=0.049 Some representative quotes in response to onestop:   _"MANOVA works well in situations where there are moderate correlations between DVs" (course notes from San Francisco State Uni)   "The dependent variables are correlated which is appropriate for Manova" (United States EPA Stats Primer)   "The dependent variables should be related conceptually, and they should be correlated with one another at a low to moderate level." (Course notes from Northern Arizona University)   "DVs correlated from about .3 to about.7 are eligible" (Maxwell 2001, Journal of Consumer Psychology)_ n.b. I'm not referring to the assumption that the intercorrelation between Y1 and Y2 should be the same across all levels of independent variables, simply to this apparent grey area about the actual magnitude of the intercorrelation.


title: Estimating the dimension of a data set

text: A colleague in applied statistics sent me this: > "I was wondering if you know any way to find out the true dimension of the > domain of a function. For example, a circle is a one dimensional function in > a two dimensional space. If I do not know how to draw, is there a statistic > that I can compute that tells me that it is a one dimensional object in a > two dimensional space? I have to do this in high dimensional situations so > cannot draw pictures. Any help will be greatly appreciated." The notion of dimension here is obviously ill-defined. I mean, I could run a curve through any finite collection of points in high dimensional space, and say that my data is one-dimensional. But, depending on the configuration, there may be an easier or more efficient way to describe the data as a higher dimensional set. Such issues must have been considered in the statistics literature, but I'm not familiar with it. Any pointers or suggestions? Thanks!


title: Can someone explain the concept of nontrivial joint sufficient statistic for a pair of parameters?

text: How to find joint sufficient statistic for a pair of unknown parameters from a bivariate random sample? Also, what exactly does nontrivial mean in this context? Please elaborate and be simple in your illustration, as much as you can.


title: How do I compute the density of this distribution

text: Ferdosi et al define six artificial datasets to test density estimation methods. Part of the fourth dataset is defined as: $$ Uniform(x,y) = [0,100], Gaussian(z) = [M = 50, var = 5] $$ Where $M$ is defined as the mean of the distribution and $var$ as the variance. To test my own density estimator using this dataset I'd like to do two things:   1. Sample $N$ items from this dataset.    2. Compute the 'true' density of each of these samples drawn. ## Sampling For step (1) I tried the following in Matlab:   1. Draw $N$ samples from a uniform distribution (this gives a $N \times 1$ matrix).              uniform = makedist('uniform', 'lower', 0, 'upper', 100);     x = uniform.random(N, 1);        2. Draw $N$ samples from a uniform distribution (this gives a $N \times 1$ matrix).               y = uniform.random(N,1);        3. Draw $N$ samples from a Gaussian distribution (resulting in a $N \times 1$ matrix).              z = mvnrnd(M,var,N);        4. Combine these two matrices to a $N \times 3$ matrix.              data = [x, y, z];      This seems to give the correct results, that is to say the resulting plot is the same as the plot shown in the article. Is this the correct way to sample this distribution? If not how should I go about it? Is there a difference between sampling three different vectors or sampling the $x$ and $y$ vector from a bivariate distribution? ## Computing the density This wikipedia artikel defines ''the joint probability density function $fX,Y(x, y)$ for continuous random variables as: $$ f_{X,Y}(x,y) = f_{Y\mid X}(y|x)f_X(x) = f_{X\mid Y}(x\mid y)f_Y(y)\ $$ where fY|X(y|x) and fX|Y(x|y) give the conditional distributions of Y given X = x and of X given Y = y respectively, and fX(x) and fY(y) give the marginal distributions for X and Y respectively.'' But I don't know how to use it, if it is the correct way, to determine the density of the defined distribution.


title: peak detection - points / sigma ratio

text: I would like to ask you (statistics experts) if the following approach is kosher or a nonsense. ### Problem My EM based detection algorithm ends up (for some data) with a result that looks like attached picture. With second peak I am satisfied, but first one is only a noise. ![enter image description here](http://i.stack.imgur.com/oukSY.png) ### Approach Following procedure proved to be useful for me...   1. Take standard deviation of cluster ($sigma$) and number of points in cluster ($N$)   2. Calculate $ratio = N / sigma$   3. If $ratio < threshold$ then remove the cluster from a resulting list of clusters ### Questions   1. I suspect that this ratio will tell nothing about if data in cluster are normally distributed or not. However, since I want only nice bell shaped peaks that fit to data histogram, this should be a solution for me. Is it correct approach, or am I doing something which is "nasty" for statistics?   2. My threshold value is based on my observations I made on data from a device, which I have been developing software for. Is there any way how to calculate optimal value for threshold?


title: Choice of loss function in correlation matrix prediction

text: There is a random vector $X=(X_1,\ldots,X_p)$, with $p$ large, $E[X]=0$ and $V[X_j]=1\ \forall \, j=1,\ldots,p$, but the correlations are different from zero. We cannot assume multivariate normality for $X$ (the marginal distributions are slightly skewed and heavy-tailed, as common in financial data). Say we have a training set of $n$ realizations of the vector $X$, and we have some method of estimating the $p\times p$ correlation matrix. The goal is to compare the performance of several methods for predicting the correlation matrix in new data (hold-out/validation set). _What is a good loss function between the estimated covariance matrix and the observed covariance matrix in the hold-out set?_ So far, to my mind come 1. a predictive quasi-normal log-likelihood, 2. sum of absolute/squared errors between the off-diagonal elements of the estimated correlation matrix and the observed correlation matrix of the held-out data. Any better suggestions?


title: Trying to understand NMF

text: I'm trying to understand how NMF derived, and I got the basic idea of NMF, that is, it tries to appoximate original matrix $V$ with $WH$, where $V$ are non-negative, and $W,H$ are constrained to be non-negative. **Question 1** When I approach this problem, it's pretty intuitive for me that, I should try to minimize the reconstruction error, $$Err(W,H) = \min_{W,H}\|V-WH\|$$ , if say $Err(W,H)$ is the objective function, then I have a problem to define the constraints, how to constrain $W,H$ to be non-negative? **Question 2** The paper I'm reading, it doesn't take the reconstruction error as the objective function, instead it takes the perspetive that NMF tries to reconstruct the original $V$ from $WH$ adding some **Poisson** noise (last 2 paragraphs of page 3). This confuses me badly, why **Poisson**? What about Gaussian or gamma or beta? **Question 3** No matter what the objective function is, I think I can always solve the problem by **gradient descent** or **newton's method** , can I? And the reason why the **multiplicative update** method proposed is because it has better speed than gradient descent?


title: How to determine the accuracy of a multiple linear regression model?

text: I searched a method to determine the accuracy of a linear regression model. I found that I should calculate r-squared. Is this the only method or are there other methods?


title: Variance of residual when error variance is known

text: Suppose the true model is $y = 2x + \text{error}$, where $\text{error}$ is normally distributed with variance 1 and mean 0. See the R code below. Variance for any residual from the fitted model (under OLS) is given by $\text{Var}(\text{error}) \times (1 - \text{hatvalue})$. Ex: var(residual) for the first observation is 0.99. Question: How can I verify that var(for residual for the first observation) is indeed equal to 0.99 using the data? Do I have to do sampling a good number of times but with different y values which is due to uncertainity in error and then fit a model, next compute var(of all the residuals obtained at observation 1 using all samples) and then check if it is atleast close to 0.99? Appreciate if anyone can please clarify.                x <- rchisq(300,1)      error <- rnorm(300,0,1)      #true model      y <- 2*x + error      mod <- lm(y~x)      h <- hatvalues(mod)      #var_ei = var(epsilon) * (1-h(i))      #for example - for the first residual i = 1      #since variance(epsilon) = 1      var_e1 = 1 * (1 - h[[1]])      #var_e1 = 0.99


title: Hausman test for Independence of Irrelevant Alternatives (IIA), negative value

text: I get a negative value for the Hausman test for the independence of irrelevant alternatives (IIA) assumption. Now I find contradicting conclusions in the literature: some people say that I should take the absolute value, others say that the value should be rounded to $0$. What should I do?


title: How is a Sparse RBM different from a Gaussian-Bernoulli RBM

text: Sparse RBM are described in this paper Gaussian-Bernoulli RBMs are describe (kinda poorly) in this paper as Gaussian Units, and more clearly in this masters thesis. Sparse RBM: (Quoting direct from the paper, Page 3) $$P(v_{i}|\mathbf{h})=\mathcal{N}(c_{i}+\sum_{j}w_{ij}h_{j},\,\sigma^{2})$$ $$P(h_{j}|\mathbf{v})=logistic\left(\frac{1}{\sigma^{2}}(b_{j}+\sum_{i}w_{ij}v_{i}\right)$$ Gaussian-Bernoulli RBMs (Quoting less directly from the masters thesis., pages 38 and 36. Some substitutions/expansions made) $$v_{i}|\mathbf{h}\sim\mathcal{N}(c_{i}+\sum_{j}w_{ij}h_{i},\,\sigma_{i}^{2})$$ $$P(h_{j}=1|\mathbf{v})=logistic\left(b_{j}+\sum_{i}w_{ij}\frac{v_{i}}{\sigma_{i}^{2}}\right)$$ So the difference would be a scaling factor of $\frac{1}{\sigma}$ on $b_i$. and $b_i$ is a learned factor controlling the threshold, so will learn away a scaling factor. They seem to have the same PDFs, though to formulation of the energy function is superficially different, I think it simplifies to be the same. I the difference is $\sigma$ in a sparse RBM is a definable parameter -- that controls sparseness. But $\sigma_i$ in a Gaussian-Bernoulli is the variance of the the features $i$. They are both trained the same way with Contrastive Divergence. is this right? are there more connections between them? Have i got it all wrong?


title: Using fitted ARIMA model to forecast new time series

text: I'm new to ARIMA analysis and I'm trying to understand how fitted ARIMA model can be used to forecast new time series, given starting point only. Should the model be refitted or the estimated parameters are sufficient to produce forecasts? Thanks


title: Confusion matrix calculation in random forest classifier in R

text: After training using random forests on the iris dataset, I get an OOB error and a confusion matrix.               data(iris)     (model <- randomForest(Species~., iris, ntree=500,importance=T,do.trace = 100) )     model$oob.times      The help mentions that the confusion matrix is based off the OOB data. But I see the confusion matrix reports values using the same number of samples as used in the training input. Can anyone explain as to how exactly the confusion matrix reports its error values? Does it use the OOB error value and scale it to the training data-set size, or does it pick samples of OOB data and run the RF again?


title: Forecast mean and variance for group data

text: Apologies if this is a bit of a simple question, but I haven't been able to find any answer to this over the past week and it's driving me crazy. Background Info: I have a dataset that tracks the weight of 5 individuals over 5 years. Each year, I have a distribution for the weight of individuals in the group, from which I calculate the mean and standard deviation. Data is as follows:                 Year = [2002,2003,2004,2005,2006]       Weights_2002 = [12, 14, 16, 18, 20]       Weights_2003 = [14, 16, 18, 20,20]       Weights_2004 = [16, 18, 20, 22, 18]       Weights_2005 = [18, 21, 22, 22, 20]       Weights_2006 = [2, 21, 19, 20, 20]      The Question: How do I project annual distributions of weight for the group the next 10 years? Ideally, I would like the uncertainty about the mean to increase as time goes on. Likewise, I would like the uncertainty about the standard deviation to increase too. Phrased another way, I would like to project the distributions of weight going forward, accounting for both:   1. Natural Variance in the Data   2. Increasing uncertainty. Any help would be greatly, greatly appreciated. Is this even possible? If anyone can suggest how to do this in R, that would be even better. I've tried using the Holt-Winters method but it seems that I can only input the mean. Thanks guys!


title: Measuring distribution of growth rates

text: I have monthly sales data for 50,000 companies. Unfortunately, I have very limited statistics training. I suspect/hope that by some measure of growth these companies would fit into a standard normal distribution, but I can't figure out what that measure would be. By default we measure growth as a ratio, e.g. monthly growth is               (this_months_sales / last_months_sales) - 1      But since sales are never negative, the range of this number is -1 to infinity. I can't take the log of this number, either, since it goes through zero. I want to get a mean and standard deviation of growth so I can characterize relative performance consistently. I mean, how much harder is it to achieve 100x growth compared to 50x growth? Is that as hard as going from 1x growth to 2x growth? How should I proceed?


title: Posterior distribution for LDA and Newdata

text: I am using the 'topicmodels' package in R. I tested the posteriori probability for newdata over jss_LDA result by this code :               post <- posterior(jss_LDA, newdata = dtmNew)      For each words i have the probability over all topics. What is the signifiance of each value probability in terms of newdata?


title: How to determine the factors correlated with observed data?

text: I have box-office collection data on a number of movies. I also have the production budget, director name, lead actor, actress, language and other meta data related to the movie. I want to know which factors determine the box office collection. If a certain factor (say the director) is correlated with the collection, then I would like to determine its effect. I was initially thinking of linear regression with dummy variables. However, I realised that I do not want to know whether a certain director, certain actor or certain language play a role, but rather whether direction, actor, etc. have any role in determining the revenue of the movie collectively. I would appreciate any help about the possible directions I should pursue.


title: Recommended effect size measure for repeated measures t-tests and repeated measures contrasts on experimental data

text: ### Context This came up recently in a consulting context. A researcher was performing repeated measures t-tests based on experimental data. Some of the analyses involved comparing one condition with another. Other analyses involved performing contrasts comparing one or more conditions with one or more other conditions. ### Question   * What effect size measure would you recommend using in relation to (a) repeated measures t-tests (b) repeated measures contrasts comparing one or more group means with one or more other group means?   * If you were reporting a d-based effect size measure, which measure of the within group standard deviation would you use?   * Are there any references that you would recommend? I have a few thoughts, but I'm keen to get your suggestions.


title: Visualizing multi-dimensional data (LSI) in 2D

text: I'm using latent semantic indexing to find similarities between documents (thanks, JMS!) After dimension reduction, I've tried k-means clustering to group the documents into clusters, which works very well. But I'd like to go a bit further, and visualize the documents as a set of nodes, where the distance between any two nodes is inversely proportional to their similarity (nodes that are highly similar are close together). It strikes me that I can't accurately reduce a similarity matrix to a 2-dimensional graph since my data is > 2 dimensions. So my first question: is there a standard way to do this? Could I just reduce my data to two dimensions and then plot them as the X and Y axis, and would that suffice for a group of ~100-200 documents? If this is the solution, is it better to reduce my data to 2 dimensions from the start, or is there any way to pick the two "best" dimensions from my multi- dimensional data? I am using Python and the gensim library if that makes a difference. Thanks!


title: Using latent Dirichlet allocation for information retrieval

text: I am working on understanding various document ranking algorithms like (TF- IDF, LSI, language models, etc) by actually implementing them. I want to understand LDA and using various resources to understand the algorithm. What I don't understand is how we come up with the latent (hidden) variables/topics. Can someone please explain it to me using examples like:               Doc1: "shipment of gold damaged in a fire",     Doc2: "delivery of silver arrived in a silver truck",      Doc3: "shipment of gold arrived in a truck"};      Query: "gold silver truck"      I will really appreciate any help in this regard. Thanks!


title: Ordinal data, different scales, comparison between variables

text: I would appreciate any help you could give me I have ordinal data which I want to standardise in order to compare averages and data spread between variables. I need to standardise since the variables have different scales (one is out of 7, the other out of 15). I have transformed the original values into percentages and then found the mean and standard deviations for each variable in order to allow for comparison. Is this ok? Somehow converting into percentage values and then taking the mean and SD of this doesn't seem right... Any help would me much appreciated :)


title: ARMA regression

text: I am studying stock price. I first define that the ln(price) has an obvious unit root of 1 lag. I was wondering if someone could explain me how to proceed to construct a correct ARMA model...? Thanks is advance.


title: What can cause PCA to worsen results of a classifier?

text: I have a classifier that I'm doing cross-validation on, along with a hundred or so features that I'm doing forward selection on to find optimal combinations of features. I also compare this against running the same experiments with PCA, where I take the potential features, apply SVD, transform the original signals onto the new coordinate space, and use the top $k$ features in my forward selection process. My intuition was that PCA would improve the results, as the signals would be more "informative" than the original features. Is my naive understanding of PCA leading me into trouble? Can anyone suggest some of the common reasons why PCA may improve results in some situations, but worsen them in others? Thanks


title: Are very large log likelihood and delta AIC values problematic for model selection?

text: I am using AICc for small sample sizes to compare 8 _a priori_ models (including null model). I fitted my models using a GLMM due to the nested nature of my data and defined the family as 'poisson' based on a visual inspection of the error structure from my residual vs fitted values plot. After running my analysis using lmer from the lme4 package in R, I obtain very large log likelihood values, e.g., LL = -120995125, which makes the delta AICc values very large as well, e.g., difference between 'best' fit model and second ranking model is 12422654. This makes comparison within the model set rather dubious. I was wondering if anyone has encountered a similar problem or if there is something I am missing when using family=poisson for lmer in my model fitting.


title: What does taking the logarithm of a variable mean?

text: > **Possible Duplicate:**   >  When are Log scales appropriate? Say you have a bar graph displaying data, for example, "Cost of Computer Orders by Population", and you are trying to analyze the data and find a distribution. The information does not indicate anything so you take the logarithm of the variable and the graph then resembles a normal distribution. What does taking the logarithm of the information indicate?


title: Bayesian approach to trend detection in non-parametric data

text: I have a series of data points, and I want to see how much evidence there is that the points are getting bigger over time. The data themselves are counts, but for various reasons I don't want to build a model of this as a Poisson process with potentially changing lambda. I'd like to only look at order statistics. I know there are various ad-hock approaches I could take (like randomly sampling pairs of points, and seeing how often the most recent of the two is the larger). However, what I'd really like are recommendations for the best approach to take from a Bayesian viewpoint. In particular, I'd like an approach that can continue to evaluate new data points as they come in, continually updating prior beliefs about the strength of the trend. Implicit in this, obviously, is some kind of non-parametric measure of the strength of the trend, so it seems I'm asking for a suggestion for that as well. Hopefully you can excuse the two-fer on account of the clear dependency between the metric and the approach.


title: Data entry tool for sparse table

text: Please suggest a data entry tool (if one exists) for a single large sparse table. Something like... row x, checkbox list of 300 columns, check 5 that pertain, enter data for those 5. row x+1, same checkbox, check different 14 that pertain, enter data for those 14. row x+2, checkbox does not contain column, add column title and data type, check new column created and 33 others that pertain, enter data for those 34. Although possible, it's inconvenient to do in a worksheet - too much tabbing/arrowing around.


title: How to demonstrate failure of CLT in R?

text: I was assigned to demonstrate the Central Limit Theorem (CLT) in R in my statistics class. I already made some progress with simulation using simple.sim in R. I want to prepare 3 examples of the CLT with continuous distributions ( _rnorm_ , _rt_ and _rf_ ) and 3 examples with discrete distributions ( _rbinom_ , _rpois_ and _rtriangle_ for triangle distribution). I also want to demonstrate the failure of the CLT. I realized that this could happen when the variance of the distribution does not converge. I can't find any discrete distribution that has non-converging variance and is implemented in R.


title: Stepwise versus L2 regularized logistic regression: dataset-specific performance

text: I have two data sets from different collections. The second data set is smaller. They were both analyzed with the same methods in order to derive feature sets of 10-30 features each. Each feature set was produced the same way for both data sets. Then, I run many Logistic Regressions to fit both data sets with all feature sets. Additionally, all of the experiments were repeated with both **L2 regularized** and **Stepwise** Logistic Regression. The observation is that the best fitting for the first data set was done with Stepwise Logistic Regression, while for the second one was done with L2 regularized Logistic Regression. This is quite consistent, i.e. for all 15 experiments on each dataset. _Why did each method perform better on either data-set? May this have to do with particular characteristics of each data-set?_ For example, I know that L2 deals better with multicollinearity and lower ratio of observations/variables. Can I assume that L2 performed better on 2nd data set because it had multicollinearity? Then, L2 does not zero out any coefficients, which Stepwise in fact does. Can I say that Stepwise did better on 1st data set because it may have had some too noisy features that needed to be zeroed out?


title: Are random effect models same as Bayesian versions of fixed effect models?

text: Are random effect models same as Bayesian version of fixed effect models, in the sense that the parameters or coefficients are random variables? For example, is the Bayesian linear regression model viewed as a random effect model? Thanks!


title: Showing that an ordinal logistic model with a proportional odds assumptions is the binary logistic regression with a different prior

text: Can some one show that a Bayesian ordinal logistic model with a proportional odds assumptions is the binary logistic regression with a different prior. That is you are willing to take a multivariate normal prior for the regression parameters and an independent prior of the cutoffs, which can be truncated normal distribution. Truncation to impose order on the cutoffs.


title: R: glmulti for mixed models returns several best models, Automated Model Selection (multilevel analysis, hierarchical model, nested data)

text: I searched the entire web including this forum on some help on how to use the glmulti package in order to identify the "optimal" fixed part of a mixed model with a given random part. However, I could not find anything to solve my problem. PLEASE HELP! I followed the introduction provided in the glmulti.pdf, which you can find in the R packages folder after installing the glmulti package. The first step, I'm trying to do, is to identify the "optimal" fixed part of a mixed model with a given random part. However, the output of glmulti contains a number of different models (more than 50), which are all evaluated same as good by the information criteria (IC). This does not make sense to me! I also used the f-test and t-test obtained by REML estimation and compared nested models using ML estimation to compare the models provided by glmulti. In terms of AIC, BIC and logLik, all models score differently, even though glmulti evaluated them the same. I also used a subset of my predictors in order to compared my results from glmulti with results I obtained by the MiMUn package, which also allows to identify the "best" mixed model across a number of candidate models. The MiMUn package identified one single "optimal" model, which is one of the models provided by the glmulit package. Using f-tests and t-tests obtained by REML estimation and by comparing nested models using ML estimation, I could verify that MiMUn selects the "optimal" model. The disadvantage of using MiMUn is, that it does not allow for as many predictors in the fixed part as glmulti. If I use MiMUn to evaluate my entire set of predictors, it process for some 24 hours before it crashes. If I use glmulti on my entire set of predictors, it comes up with more 50 models, which are evaluated same as good. Therefore, I really would like to get the glmulti package running. See bellow how I apply my code on a subset of my predictors:                 ## ## the wrapper       lmer.gmulti <- function (formula, data, random = "", ...) {       lmer (paste (deparse (formula), random), data = data, REML = F, ...)       }            ## the model       test = glmulti(log(y) ~                   a + b + c + d + e + f + g,                 level=1,                 fitfunc = lmer.gmulti,                 random = "+(1| x)",                data= my_data)      All predictors are continuous variables, in spite of the random intercept which is categorical. Alternatively, I tried a different wrapper function, which I found in a documentation of an older version of the glmulti package. I tried it on the most recent version of glmulti and on the former version 0.7. Both did not produce the correct result. See the code below:                    lmer.gmulti_2011 <- function (formula, data, random = "", ...) {          lmer (paste (deparse (formula), random), data = data, ...)          }


title: What is the best test to compare pattern similarity in different condition?

text: I have several subjects, for each subject I have two measure (pre and post treatment) of several parameters, I have calculated a multivariate similarity measure (Mahalanobis distance) between each dataset and a particular condition (target). My similarity measure varies across time and I need to measure if there is an effect of the treatment, in other words if the parameters post treatment are more similar to the target parameters than the measures taken before treatment. Could I compare the distribution of the similarity using a Kolmogorov-Smirnov test and state that the two distributions are statistically different? Should I test each subject separately or measure the distribution using similarity from all subjects?


title: Testing whether or not a result is significantly different from another?

text: If I have a table with results (columns being different feature subsets and rows being different machine learning algorithms) how can I test whether or not some of the results (0-1 loss) is not significantly different from the best result? I have found a wide range of different statistical significance tests, but non of them seem to fit the purpose entirely. A solution in Excel or Python will be preferable. ![An example...](http://i.stack.imgur.com/JRkmm.png) I'm especially thrown off about which information is needed? Is the table with results enough in itself?


title: Statistical error in Bayesian framework

text: Residuals and errors are related but not exchangeable. In Wikipedia I read: > In statistics and optimization, statistical errors and residuals are two > closely related and easily confused measures of the deviation of an observed > value of an element of a statistical sample from its "theoretical value". > The error of an observed value is the deviation of the observed value from > the (unobservable) true function value, while the residual of an observed > value is the difference between the observed value and the estimated > function value. I would agree with this - but only if we look at it from a frequentist perspective. That is, we have to know the one true population value in order to be able to talk about statistical errors. Now, as far as I understand, in the Bayesian framework there are no such things as true values. What does that mean for the use of the term "error"? In principle one could probably still talk about "errors" and assume that they, themselves have another underlying distribution but somehow it seems not correct to me to use the above definition in a framework without true values? Any thoughts on this?


title: What is the relevance of bootstrapped confidence intervals on markov chain transition matrices?

text: Am I right in thinking that bootstrapping is just replicating your original data x amount of times? So I have a series of $n=20$ sequences of varying length containing the letters ACGTE. The probability of having C after to A is $$ P(X_C=x_C|X_A=x_A), $$ and which can be calculated from the set of sequences, thus $$ \widehat{P_{ij}}=\displaystyle \frac{x_{ij}}{\sum^mx_{ik}}\quad \text{where m=ACGT or E}. $$ Assume this is $\widehat{P}\rightarrow P$ as $n\rightarrow \infty$. However I don't have enough sequences to produce a full transition matrix $P$, instead I have a sparse matrix: \begin{align} \widehat{P}=&\begin{pmatrix} 0& 0 & 0& 1.00& 0\\ 0& 0 & 0.22& 0.44& 0.33\\ 0& 0.19& 0.57& 0.10& 0.14\\ 0.03& 0.19& 0.09& 0.50& 0.19\\ 0& 0.06& 0.13& 0.56& 0.25\\ \end{pmatrix} . \end{align} This source suggests that I can use bootstrapping to increase my sample size and get a better approximationg of $P$: arxiv.org/pdf/0905.4131v1.pdf However, what does this really tell me? Nothing more useful if I don't have more real data right? How do confidence intervals (http://hal.inria.fr/docs/00/46/88/04/PDF/estimvar_final.pdf) make it credible all of a sudden?


title: Two methods of bootstrap significance tests

text: Using bootstrap I calculate p values of significance tests using two methods:   1. resampling under the null hypothesis and counting the outcomes at least as extreme as the outcome coming from the original data   2. resampling under the alternative hypothesis and counting the outcomes at least as distant from the original outcome as the value corresponding to the null hypothesis I believe that the 1st approach is entirely correct as it follows the definition of a p value. I'm less sure about the second, but it usually gives very similar results and reminds me a Wald test. Am I right? Are both methods correct? Are they identical (for large samples)? * * * Examples for the two methods (edits after DWin's questions and Erik's answer): > Example 1. Let's construct a bootstrap test similar to the two sample T > test. Method 1 will resample from one sample (obtained by pooling the > original two). Method 2 will resample from both samples independently. > > Example 2. Let's construct a bootstrap test of correlation between x₁…xₐ and > y₁…yₐ. Method 1 will assume no correlation and resample allowing for (xₑ,yₔ) > pairs where e≠ə. Method 2 will compile a bootstrap sample of the original > (x,y) pairs. > > Example 3. Let's construct a bootstrap test to check if a coin is fair. > Method 1 will create random samples setting Pr(head)=Pr(tail)=½. Method 2 > will resample the sample of experimental head/tail values and compare the > proportions to ½.


title: One sample t-test, unsure of test value

text: I would like to use a one sample t-test where accuracy is 50-50.. The data is rated on a three point likert scale of 1-performed task incorrectly, 2 - performed task with minor mistakes, 3-performed task accurately I have calculated the mean and standard deviation of my data and that is fine, I am just confused about what test value to use.. (newbie) Do I use .5 as in 50% or do I use 1.5 as this is the center of the possible high score of 3 that can be achieved?


title: Is a post-hoc test supposed to be performed on transformed or original data?

text: This might be a silly question, but I could not find an equivalent one on CrossValidated. I was wondering: if I transform my data in order to meet the assumption of homoscedasticity before running a two-way ANOVA, should I also use the transformed data to perform Tukey's post-hoc test and find out which treatments are significantly different from each other? I can sense I will get an answer that tells me it depends on my data, the kind of transformation and the resulting data, but there might an obvious explanation or a rule of thumb somewhere.


title: Standard error of a weighted mean when observations are not independent

text: I have a situation where I have a bunch of noisy observations with known (and normal) sampling error and a non-trivial co-variance structure between the observations. In an uncorrelated setting, I'd generate weights by taking the inverse variance of the observations and calculate the standard error of the mean by taking the reciprocal of the sum of weights. How does that generalize in the face of an arbitrary correlation matrix?


title: KL divergence between two multivariate Gaussians

text: I'm having trouble deriving the KL divergence formula assuming two multivariate normal distributions. I've done the univariate case fairly easily. However, it's been quite a while since I took math stats, so I'm having some trouble extending it to the multivariate case. I'm sure I'm just missing something simple. Here's what I have... Suppose both $p$ and $q$ are the pdfs of normal distributions with means $\mu_1$ and $\mu_2$ and variances $\Sigma_1$ and $\Sigma_2$, respectively. The Kullback-Leibler distance from $q$ to $p$ is: $\int \left[\log( p(x)) - \log( q(x)) \right]\ p(x)\ dx$, which for two multivariate normals is: $\frac{1}{2}\left[\log\frac{|\Sigma_2|}{|\Sigma_1|} - d + Tr(\Sigma_2^{-1}\Sigma_1) + (\mu_2 - \mu_1)^T \Sigma_2^{-1}(\mu_2 - \mu_1)\right]$ Following the same logic as this proof, I get to about here before I get stuck: $=\int \left[ \frac{d}{2} \log\frac{|\Sigma_2|}{|\Sigma_1|} + \frac{1}{2} \left((x-\mu_2)^T\Sigma_2^{-1}(x-\mu_2) - (x-\mu_1)^T\Sigma_2^{-1}(x-\mu_1) \right) \right] \times p(x) dx$ $=\mathbb{E} \left[ \frac{d}{2} \log\frac{|\Sigma_2|}{|\Sigma_1|} + \frac{1}{2} \left((x-\mu_2)^T\Sigma_2^{-1}(x-\mu_2) - (x-\mu_1)^T\Sigma_2^{-1}(x-\mu_1) \right) \right]$ I think I have to implement the trace trick, but I'm just not sure what to do after that. Any helpful hints to put me back on the right track would be appreciated!


title: Regression with non-normally distributed residuals

text: There are several posts on this site talking about the need of normality when interpreting the meaning of the p.value of a linear regression. But not much I think is said about how to deal with non-normal data set. On this post, they give some solutions when the distribution is long-tailed. I'm dealing with a case where my residuals (and my dependent variable) has a multimodal distribution (as it can be seen on the following Kernel density graphs) and takes discrete values (as it can be seen in the other graphs). My model takes "FP" as dependent variable and "Design Complexity" and "Sample size". On the following graph, each point and whiskers represent 20 points. I strongly expected to have an effect of "Design Complexity" but didn't know wether "Sample size" may have an effect. As I couldn't reject the possibility that an interaction might exist between "Sample size" and "Design Complexity" I observed both. ![enter image description here](http://i.stack.imgur.com/mF8v3.png) Here is the R code for ths model.               summary(aov(FP~Obs.size*Design.Complexity, data=data.and.factors))      But my residuals are definitely not normally distributed: note: All the following graphs are presented in both kernel density and regular x-y plot               plot(density(residuals(aov(data.and.factors$FP~data.and.factors$Design.Complexity*data.and.factors$Obs.size))))      ![enter image description here](http://i.stack.imgur.com/J5XkS.png) ![enter image description here](http://i.stack.imgur.com/Hhf56.png) My dependent variable is not normally distributed either and its distribution change from one "Design Complexity" value to another (see below) Distribution of FP ![enter image description here](http://i.stack.imgur.com/y4YOg.png) ![enter image description here](http://i.stack.imgur.com/NGtID.png) Distributions of FP for Design Complexity equal 1, 2, 3, 4, 5 and 6. ![enter image description here](http://i.stack.imgur.com/pyx0i.png) ![enter image description here](http://i.stack.imgur.com/fiFNA.png) **How can I get trustfull p.values ?**


title: Regression analysis with variables of different lengths

text: I would like to perform multiple regression with two variables that have the same unit but have very different gradient lengths. (one has a range from 1-20, the other has a range from 1-40) How do I cope with variables of different lengths in regression? Is centering of the data enough to avoid that the variable with the longest gradient has the biggest impact? cheers Leke


title: Time to Event modeling, fixed but different durations

text: I am looking at probability of an event ($E$) for a number of customers. Each customer qualifies for the analysis through a qualifying Action ($A$), and has a finite Duration ($D$) to complete Event. There are a number of interim actions that should have an effect on reaching the Event. Unlike the time-to- event modeling I've done before, the Duration is known ahead of time, but is different for every customer. That is, customer 1 may have a duration of 3 months from their qualifying Action; whereas customer 2 may have a duration of 3 weeks. The duration is always known at the time of the qualifying Action. Consider an advance car rental booking. Making the reservation is the qualifying Action. The Event of interest is whether the customer completes a full online profile. The Duration is the time period between contract & pickup. Interim actions might include receiving an email from the agency, calling into a help desk, starting the profile. I think that incorporating duration is important because a customer booking 1 day in advance might (should) have a different propensity to complete an online profile than one who books 6 months in advance. My instinct is to transform the duration to a proportion of the finite duration. Concretely, qualifying Action is time 0 and everyone has a duration of 100. But I've not seen that actually done in any of my literature review, and it feels like this may lose important information.


title: Comparing 8 categories (degree) against % who did a gap year within that degree - which test?

text: Want to see if there is a statistically significant difference between the % of students who did a gap year across 8 different degrees. $$ \begin{array}{l|r} \text{Degree} & \text{Gap } [\%] \\\hline\\\hline \text{Architecture}& 8.82 \\\hline \text{Biology}& 22.41 \\\hline \text{Languages}& 22.22 \\\hline \text{Mathematical Sciences}& 4.55 \\\hline \text{Mechanical Engineering}& 12.43 \\\hline \text{Pharmacy}& 30.00 \\\hline \text{Politics}& 37.04 \\\hline \text{Psychology}& 32.65 \end{array} $$ Which stats test is used? Been reading conflicting things online and in my stats book. Thought a one-way ANOVA would be it, so tried in SPSS but only it recognised one variable (gap %) so wouldn't run the test. ![Screenshot of SPSS](http://i.stack.imgur.com/wu19e.png) Perhaps the layout is wrong? Is the cohort categorical and the gap % scale?


title: Multiple imputation for outcome variables

text: I've got a dataset on agricultural trials. My response variable is a response ratio: log(treatment/control). I'm interested in what mediates the difference, so I'm running RE meta-regressions (unweighted, because is seems pretty clear that effect size is uncorrelated with variance of estimates). Each study reports grain yield, biomass yield, or both. I can't impute grain yield from studies that report biomass yield alone, because not all of the plants studied were useful for grain (sugar cane is included, for instance). But each plant that produced grain also had biomass. For missing covariates, I've been using iterative regression imputation (following Andrew Gelman's textbook chapter). It seems to give reasonable results, and the whole process is generally intuitive. Basically I predict missing values, and use those predicted values to predict missing values, and loop through each variable until each variable approximately converges (in distribution). Is there any reason why I can't use the same process to impute missing outcome data? I can probably form a relatively informative imputation model for biomass response ratio given grain response ratio, crop type, and other covariates that I have. I'd then average the coefficients and VCV's, and add the MI correction as per standard practice. But what do these coefficients measure when the outcomes themselves are imputed? Is the interpretation of the coefficients any different than standard MI for covariates? Thinking about it, I can't convince myself that this doesn't work, but I'm not really sure. Thoughts and suggestions for reading material are welcome.


title: CI for logistic regression

text: What does it means if no CI was given for binary logistic regression analysis in SPSS output? ![enter image description here](http://i.stack.imgur.com/rAGyE.jpg)


title: How to interpret type-II error probability while doing A/B testing?

text: I've been trying to familiarize myself with some of the statistics behind conducting A/B testing on a website and trying to pick a relevant sample size. I've come across several articles that talk about power analysis (this is the best one that I have come across), and I'm trying to dive a little more into the theory.   1. The 37 signals post talks about doing a power analysis by using R's power.prop.test with a type-I error probability of 5% and a type-II error probability of 20%. This means that if after I run an A/B test with the given number of samples and the metric I am looking for is higher than the value I plugged into to the power analysis, it is very probable that the second variation had a positive effect. If it is not, it is very probable that the variation had no effect or a negative effect. Is this correct?   2. Assume that I want to iterate very quickly and therefore run a power analysis with a type-II error probability of 50% or even higher. My type-I error probability is still 5%. If I run the same A/B test with the lower sample number and the metric is again higher, can I still say that it is very probable that the B variation had a positive effect? However, if it is not higher, because my type-II error possibility is so high, I can make no conclusions about whether the second variation was ineffective or not?


title: Showing $X^{-1}(\mathcal F^\prime)$ is $\sigma$-algebra

text: Let $(\Omega,\mathcal F)$ and $(\Omega^\prime,\mathcal F^\prime)$ be two measurable spaces and $X:\Omega \to \Omega^\prime $, show that $X^{-1}(\mathcal F^\prime)$ is a $\sigma$-algebra. I know that I must show three things, so I think first I write the definition of $X^{-1}(\mathcal F^\prime)$ and then I check the properties.


title: How to make sure explaining variable actually explains explained variable, and not vice versa?

text: I'm somehow stuck about the way I should verify regression models. For instance, an extremely over-simplified version of the Solow Growth Model might look like this: $$ growthrate = \beta_0 + \beta_1 savingsrate + u $$ So how can I make sure now that there is not only correlation, but causation in exactly that direction? How can I make sure that the real relation is **not** $$ savingsrate = \beta_0 + \beta_1 growthrate + u $$ I understand that _endogoneity_ deals with the case where u is related to the explaining variable; _simultanity_ deals with the case when explaining and explained variable are determined at the same time; but how about cases where explained variable is determined _before_ explaining variable? I somehow feel that I miss a very basic point here, but cannot figure out where exactly.


title: Fast linear regression robust to outliers

text: I am dealing with linear data with outliers, some of which are at more the 5 standard deviations away from the estimated regression line. I'm looking for a linear regression technique that reduces the influence of these points. So far what I did is to estimate the regression line with all the data, then discard the data point with very large squared residuals (say the top 10%) and repeated the regression without those points. In the literature there are lots of possible approaches: least trimmed squares, quantile regression , m-estimators, etc. I really don't know which approach I should try, so I'm looking for suggestions. The important for me is that the chosen method should be fast because the robust regression will be computed at each step of an optimization routine. Thanks a lot!


title: Likelihood ratio test in R

text: I am in desperate need for help. I am trying to calculate the likehood ratio test in R, but I don't have allot of experience using R. For example, to calculate the following Suppose $X_1, X_2,\ldots, X_n$ is a random sample from a normal population with mean $\mu$ and variance $16$. Find the test with the best critical region, that is, find the most powerful test, with a sample size of $n = 16$ and a significance level $\alpha = 0.05$ to test the simple null hypothesis $H_0: \mu = 10$ against the simple alternative hypothesis $H_A: \mu = 15$. This is the code I have up until now               #X values     #Added the mean=10 thanks to the comments here     X=rnorm(16,10,sd=4)          mlog1=function(media,x,sdev){     sum(-dnorm(x,mean=media, sd=sdev, log=T))     }          #L(10)     #also added sdev     out1=nlm(mlog1,10,x=X,sdev=4)     #L(15)     #also added sdev     out2=nlm(mlog1,15,x=X,sdev=4)          k1=out1$estimate         k2=out2$estimate          #L(10)/L(15)         print( k1/k2)      This is the very basic axample and the solution, done by hand, can be found here. I know you also have the LRT function, but I don't know how to apply this function. Please help, the code looks correct but it always returns 1 or 0.9999999 for any value in the parameters.


title: How to find principal components from variance-covariance matrix?

text: If I have the following var-cov matrix: $\Sigma_{A,B,C} = \left(\begin{array}{ccc} 1 & 2 & 3 \\ 2 & 4 & 6 \\ 3 & 6 & 9 \\ \end{array}\right)$ (1) How can I find a constant a, so $avv'$ becomes a projection matrix of rank 1? (2) How can I use this result to find the direction of the PCA from (A, B, C)'? (3) Is it possible to find the variance of the first principal component?


title: Phi MOA, Chi-Square and Null Hypothesis, What can be learned?

text: In short, is there anything to be learned from performing a Phi measurement of association (MOA) on a Chi-Square test where it failed to reject the null hypothesis? $$\phi = \sqrt{\chi^2/N}$$ where $\chi^2$ is the value returned by Chi-Square and $N$ is the number of observations. Or is the phi coefficient only useful on significant tests? I haven't touched stats in a few years but I couldn't find a explicit remark on the usefulness of a measurement on non-significant tests. Thanks.


title: Variational inference for nested Chinese restaurant process

text: I recently read paper by Chong Wang and David M. Blei "Variational Inference for the Nested Chinese Restaurant Process". And I couldn't understand the next part (from p.5): > The variational update functions for **_W_** and **_x_** depend on the > actual distributions we use, and deriving them is straightforward. If they > include an infinite sum then we apply similar techniques as we did for > _q(vi)_. Could anyone explain (and give an example) how to implement variational update for **_W_** and **_x_**? And if there are any other good papers/tutorials/code to understand variational inference for hLDA? Thanks for any suggestions.


title: Tools to detect jumps in a linear time series

text: I have a financial time series that has a linear down trend, but sometimes a jump happens (see image below). What statistical methods can I use to detect these jumps as early as possible? ![linear trend with jumps](http://i.stack.imgur.com/QgU9E.png)


title: How to analyze differences of two groups with replicates?

text: I am trying to see if there are differences between two groups (background genetic mutants) and some values associated to them. Every group has 4 different gene types with two replicates each and its respective value.               > data     bg              type     value     a            3motif 0.7371446     a            3motif 0.7519591     a     captured_exon 0.8030102     a     captured_exon 0.8110690     a  helicase_protein 0.7331316     a  helicase_protein 0.7479729     a replicase_protein 0.6897468     a replicase_protein 0.6967064     b            3motif 0.6934017     b            3motif 0.7076932     b     captured_exon 0.7725613     b     captured_exon 0.7816896     b  helicase_protein 0.6938436     b  helicase_protein 0.7153722     b replicase_protein 0.6630786     b replicase_protein 0.6831855      I'd first like to see if there are differences between the different types of each group (background). So, while looking for a possible answer I found that the following anova test would fit on my needs:               >aov(value~bg+Error(type), data = data)          Call:         aov(formula = value ~ bg + Error(type), data = data)          Grand Mean: 0.7300979           Stratum 1: type          Terms:                      Residuals     Sum of Squares  0.02462776     Deg. of Freedom          3          Residual standard error: 0.09060494           Stratum 2: Within          Terms:                              bg   Residuals     Sum of Squares  0.004222238 0.001159024     Deg. of Freedom           1          11          Residual standard error: 0.01026479      Estimated effects are balanced      But I don't know how to iterpret the results, would you please help me?


title: Calculate the error variance in a linear regression model

text: I am trying to calculate the error variance for the following question but I don't have clue where to start. Could anyone please help? ![enter image description here](http://i.stack.imgur.com/yaxsG.jpg) ![enter image description here](http://i.stack.imgur.com/gC1ha.jpg)


title: Mixture of probits: understanding truncated-based likelihoods

text: I am trying to implement a mixture model of probits to infer the best decision boundary for every latent subpopulation. When doing Gibbs sampling, we eventually have to compute $P(y^* | w_c)$ where $y^*$ is the auxilary variable in the augmented version and $w_c$ are the weights or the $\beta$ coefficients in a regression. (for a quick reference to the simple version see Wikipedia. The mixture model does not affect the way $y^*$ and $w$ are sampled.) Debugging my code, I found that the likelihood of $p(y^* | w_c)$ gives odd results since sometimes they prefer the wrong decision boundary (Test 1). However, the posterior samples pf $w | y^*$ seem to be ok (Test 2) **Test 1 (fail)** : I sample all the $y^*$ given the true $w_{true}$. Then, given the sampled $y^*$ I compare their likelihood given $w_{true}$, $P(y^* | w_{true})$, with their likelihood given some other $w$, $P(y^* | w_{wrong}) $ (in the mixture model, this comparison is done at every gibbs iteration to decide whether this subpopulation should jump to another cluster). What happens is that $P(y^* | w_{wrong}) > P(y^* | w_{true})$. The original is like this, where the colors represent the value of y* for those with $y>0$ (upper points) and $y<0$ (lower points). The likelihood of this $w$ is -172.0354 ![Data and true boundary](http://i.stack.imgur.com/Pzz39.png) Now the likelihood of a clearly wrong $w_c$ is... higher!: ![Data and and false boundary](http://i.stack.imgur.com/j0chG.png) This is another example taken from the real Gibbs sampler. According to this, the third line is more likely than the second one (the true one): ![gibbs](http://i.stack.imgur.com/fzeqQ.png) **Test 2 (ok): Draws from posterior $w | y^*$**. To be sure that given those $y^*$ the most likely $w$ are the ones around the true $w$, I generated samples from the posterior $w | y^*, b, v, y$ where b and v are the prior mean and variance of the normal distribution that generates the w. The result is that no matter my priors b and v, the y* make the posterior to very centered around the true $w$: ![w posteriore](http://i.stack.imgur.com/1nZmU.png) My likelihood function (python) is:               logsum = 0     for n in range(len(y_user)):         # if use R function, do not normalize limits a and b         if y_user[n]==1:             low = 0.000001 - dot(X_user[n], w_cluster) # 0 is NOT included in the range             upp = np.inf         else:             low = -np.inf             upp = 0 - dot(X_user[n], w_cluster) # 0 IS included in the range                  likelihood = truncnorm.pdf(y_star_user[n], low, upp, loc=dot(X_user[n], w_cluster), scale=1)         logsum += np.log(likelihood)     return logsum        How is it possible that some clearly wrong lines are more likely than the true ones? Is my likelihood function wrong ? **Update** : I normalize the _w_ all the time to avoid numerical issues in the truncated normal function. I think this is a mistake since you avoid the classifier to be very confident about points that are far from the line. But I can't see whether this might be the cause of the problem. Maybe they both the good classifier and the bad are so "shy" that the data is, "by chance" (a missfortunate set of $y^*$) making the bad one win?


title: Log-Likelihood Ratio Test: Difference of Equations

text: I have noticed that the sign at the front of the log-likelihood ratio formula changes depending on what source I have been looking at. Below are the two formula I have seen, exactly as they were written (in case the notation is offering something I am missing): $LLR = -2 ( LL_{\rm Simple} - LL_{\rm General})$ (Source) or $LLR = 2(LL_{\rm new} - LL_{\rm baseline}$) (Source: _Discovering Statistics Using SPSS_ , By Andy Field) I was wondering if someone could tell me the significance of the sign at the front of the model? My suspicion is that it is something to do with whether you are introducing or removing parameters. Thanks.


title: Ways to Estimate Date of Conception

text: Why don't epidemiologists ask people when they had sex? Why do they estimate time of conception from the birth date? This seems to introduce unnecessary error. Wouldn't it be better to record the time you had sex in some sort of database (assuming it was the time it resulted in pregnancy)? Note that fertilization of the egg by sperm can happen without sex.


title: Updating ARIMA model

text: My question is about updating the parameters of a regression with ARIMA errors model as new (monthly) data becomes available each month. Similar question were asked here before:   * Updating ARIMA models at frequent intervals   * How to update forecasts in an ARIMA model?   * What do you consider a new model versus an updated model (time series)? My concern is, when do I know (or how can I test) that my model structure (ARMA terms, exogenous variables) is still suitable? Say in the simplest case my model is y = const + b*x + e, where e ~ ARIMA I estimated in my fitting sample b=-3. A new data point (of the next month) is available and I re-fit the model and obtain b = -3.5 and so on. **When do I know, that something is wrong or how can I continuously (with each update) test my model for consistency?** I know, that the Chow-Test a useful test for model stability. But I always have only one data-point to update. My idea is:   * Check each month the one-step-ahead forecast with the actual value. If the model gives still reasonable results -> keep it.    * Collect 6 new monthly data points (so half yearly) and check with a Chow-test the model stability between the initial fit-sample and the new sample (6 monthly data points).


title: How to do model selection for Cox Proportional Hazards Model when stratifying?

text: I am performing backward selection for the Cox Proportional Hazards model with several variables (it is required that I do backward selection). One of these variables (let's call it "A") does not meet the proportionality assumption in the full model. So I would like to consider the possibility of stratifying on variable A. How should I perform backward selection in this scenario? I am considering performing backward selection on the full model while stratifying on A (call the model selected "model B"). I will then compare fit statistics between model B with stratification and model B without stratification. If I see that the stratification is providing only a marginally better fit, then I will perform another backward selection on model B without stratification. Two questions:   1. Does this seem reasonable or is there a better way?   2. How do I formally test if model B without stratification is better than model B with stratification? A likelihood ratio test would not work because the models are not nested and comparing AIC or BIC doesn't make sense because there are the same number of variables in each model. Edit: any general advice on model selection in COX PH is appreciated as well!


title: Feature selection before SVM

text: I have a simple but difficult question. Does feature selection before SVM help? I have a data set that has ~1100 features but a lot of these are redundant data / uncorrelated data. Can someone give me a reason why it should / should not help? Thanks a lot.


title: How to calculate purity?

text: In cluster analysis how do we calculate purity? What's the equation? I'm not looking for a code to do it for me. ![enter image description here](http://i.stack.imgur.com/5YcWc.png) Let $\omega_k$ be cluster k, and $c_j$ be class j. So is purity practically accuracy? it looks like were summing the amount of truly classified class per cluster over the sample size. equation source The question is what is the relationship between the output and the input? If there's Truly Positive(TP), Truly Negative (TN), Falsely Positive(FP), Falsely Negative (FN). Is it $Purity = \frac{TP_K}{(TP+TN+FP+FN)}$?


title: Multiple Regression: Variation Explained

text: I'm stuck on a homework problem and I'm hoping you all might be able to help me out. I'd be very appreciative! We were given this data set: Data. I ran a regression in Minitab with GPA as the response and IQ and Self-Concept as the predictors. I've answered a number of questions correctly with this regression, but I can't figure out how to interpret the data to answer the following two questions:   1. What percent of the variation in GPA is explained by IQ?     2. Additional amount of variation explained due to the addition of Self-Concept to the model? This is what I got from Minitab: ![enter image description here](http://i.stack.imgur.com/KWb2A.png) Is this enough information to discern the answer to these questions? If so, how do I interpret this information to find the percent of variation explained by each predictor?


title: svm classification visualization

text: I wish to make a visualization of the separating hyperplane in the SVM algorithm along with my training features. As my feature vectors have 8 dimensions, we cannot directly visualize this. The only potential solution I can think of is to get a projection down to 2 dimensional space.   1. Will this solution even make sense? I know that SVMs will often make the data linearly separable in higher-dimensions and this may not occur in 2 dimensional space. Would this at least be a decent approximation?   2. How do people visualize their classification and their training features in practice? Thanks for all the help!


title: How to construct parametric bootstrap confidence intervals in R without using boot package?

text: How can I construct parametric bootstrap confidence intervals in R without using the `boot` package?


title: Name for a classifier which builds models on disjoint input space regions

text: I wrote a meta classifier which instead of aggregating various models on the whole data (like bagging, boosting, random forest), it splits the input space into some disjoint regions and builds a base model for each region. The reasons why one would want to try this approach would be:   * remove all cross-regional interferences, for example build one model for each US state, thus letting the model to fit only on data from a specific state   * for computational reasons, when too much data is available and there is such an input variable which can partition the whole space in smaller regions The input space is split according to a disjoint predicate set. Because the input space regions are disjoint, at prediction time, for each input data point, a proper base model would be selected for fitting the results. The question is how could be named this procedure. My expectation is that there is some theory and this kind of meta algorithm has already a name. If not, what would be a proper name for it?


title: How to report ANOVA F-values when they are zero?

text: I conducted an ANOVA and one of the F-values for a main effect is .000 on SPSS, and the other is .003 on SPSS. As I am rounding to two decimal places, how do I report these? Does F = 0.00, or is F < 0.01, or do I make an exception and report a third decimal place for the .003?


title: How to perform a spectral density analysis in R using the multitaper package?

text: To Burr and other experts on this topic: I tested, well rather ad-hocly the R-package `multitaper`, and I succeeded getting graphical outputs. I got graphs for (a) spectrum, (b) spectrum with confidence interval and a graph for the F-test.   I have a couple of questions related to the package and the output which I hope will be commented.   1. What is the F-test testing – what is the exact null-hypothesis (I suppose the null is the “white-noise” hypothesis – a flat spectral density)?    2. Why is not the critical value for the F-test integrated in the F-test plot (as a horizontal line)? The package does not produce any numerical output from the test so we have no idea how “close” to reject the null-hypothesis.   3. Is the significance level in the F-test optional, for example 5 or 10%?   4. When you evaluate the spectrum, how do you transform frequencies to years when you have yearly data?   5. It is also possible to plot the spectrum with confidence interval. What is the big deal of showing the confidence interval without relating it to the null-hypothesis?    6. You recommended plotting the spectrum from spec.mtm on a log scale. Is it possible to integrate the log-command in the code:                resSpec <- spec.mtm(willamette, k=10, nw=5.0, nFFT = "default",                 centreWithSlepians = TRUE, Ftest = TRUE,                 jackknife = FALSE, maxAdaptiveIterations =100,                  plot = TRUE, na.action = na.fail)      and how do we do it?


title: Is this a repeated measures experiment?

text: We conducted an experiment on hens with seven treatments (medicines) replicated three times. Every treatment was implemented on a flock of fifteen young hens and three of these treated hens were slaughtered after 7, 14, 21, 28, and 35 days and measurements were observed. I wonder whether   * these data should be analyzed as repeated measure (split plot nomenclature)   * if Treatment$\times$Time interaction is significant then orthogonal polynomial contrast (i.e., linear, quadratic, and cubic etc.) should be used. Or there is another better approach to analyze this kind of data? I'd appreciate if someone can share/tell a similar worked example. Thanks


title: Online bachelor's degree in Statistics

text: A person who is confined to a nursing home and is very talented at statistics/mathematics desires to earn an online bachelor's degree in statistics (or related discipline) from a good university. Do you have any suggestions for this person? Thank you!


title: Empirical AUC in validation set when no TRUE zeroes

text: In a cross-validation setting (LASSO penalized logistic regression), I'm calculating AUC. However, I'm interested in the variability of these estimates over the folds (this will give me _an_ indication of the stability of my model selection over the folds). As such, I want to find the empirical AUC in each of the 10 validation sets, and then calculate the variance over them. This poses a problem, as sometimes a validation set only holds only observations that have true outcome 1 or only observations that have true outcome 0. I don't know of a way to calculate the AUC in such a setting. What would be the sensible approach here?   * Ignore this 'fold' in the calculations regarding AUC   * Give it some value anyway, like 0.5   * Perhaps you can suggest a way of approximating the AUC in such cases (adding 1 fake observation of the other kind and assume its predicted probability is either 0, 0.5 or 1?)   * Don't try this variance over the folds idea.


title: quantile threshold of t distribution and normal distribution

text: I want to calculate the quantiles of a t distribution and a normal distribution fitted to my data (can be found here: http://uploadeasy.net/upload/nk0f.rar). The t distribution is leptokurtic, therefore there will be an area, where the quantile of the t distribution is larger than the quantile of the normal distribution for the same probability level. Then there will be a threshold, at which they are the same and then an area, in which the quantile is smaller. Consider this picture: ![Quantile](http://i.stack.imgur.com/YDDOG.png) Ignore the "VaR" and read it as a quantile: At 0,05, the quantile of the leptokurtic density is smaller thant the quantile of the equivalent normal distribution. At 0,01 the quantile of the leptokurtic density is larger. I want to find this threshold. I am looking at the right side of the distribution. My code for calculating the quantiles refers to a special probability level. My questions are:   1. Are my calculations correct? Especially the plot?:   2. How can I get the threshold value? I am wondering about the following: The threshold value is NOT at the intersection of both densities at the tail right? My code gives the following picture: ![DensityTail](http://i.stack.imgur.com/FDLvY.png) I mean, is it correct, that the quantile of the t distribution is larger than the quantile of the normal distribution if a probability level is considered, which leads to quantiles before the intersection? Yes, I think so?! But the intersection is not the threshold right?! In this case, the threshold where the quantile of the t distribution is larger than the quantile of the normal distribution lies before the intersection point, does this always hold?


title: partial eta squared

text: Would it please be possible to help me with a statistics query. I have asked for the partial eta squared on the output of a one-way within subjects ANOVA from SPSS and it has given me the value .137. According to Cohen (1992) the classifications for effect sizes should be; r=.10 small, r=.30 medium and r=.50 large (I am unsure whether these classifications can be attributed to partial eta squared?) I appreciate this is a very stupid question but please may I clarify my value falls in the small effect category? Additionally, does this suggest that the proportion of variance explained by the IV on the DV and not explained by other variables in the analysis is 13.7%? Any help would be greatly appreciated!! Thank you!


title: Automating determination of number of clusters from a kmeans cluster analysis

text: I use kmeans for clustering a set of data. However, I have to specify the number of clusters. The problem is that sometimes I need 2 and other times I need 3 clusters.   * Is there a clustering algorithm that could incorporate that feature in it?


title: Linear mixed effect models: how to construct a suitable null model

text: Lets take as an example a repeated measures design with 10 subjects that are all reading the same letter strings and pressing a button as soon as they determine whether the string is valid English word, producing reaction times (RT). I wish to determine whether word length has a significant effect on the produced RT (it should), using a linear mixed effect model. Using R and the lme4 package I construct the following model:               m = lmer(RT ~ 1 + word.length + (1 + word.length|subject), data=rt.data)      As you can see, I allow both the intercept and the slope to vary randomly across subjects, as I suspect that the effect of word length might be larger for slow readers than fast readers. Understanding that p-values are not trivial in these types of models, my approach is to construct a NULL model, containing only the random effects. But I'm not sure whether this should be:               m.null = lmer(RT ~ (1 + word.length|subject), data=rt.data)      or:               m.null = lmer(RT ~ (1 | subject), data=rt.data)      In the end, I wish to perform an anova between the model with word length and the NULL model like so:               anova(m, m.null)      which should give me a p-value whether the addition of word length actually makes the model fit better and thus whether word length actually influences the RT.


title: Applying logistic regression with low event rate

text: I have a dataset in which the event rate is very low ( 40,000 out of $12\cdot10^5$). I am applying logistic regression on this. I have had a discussion with someone where it came out that logistic regression would not give good confusion matrix on such low event rate data. But because of the business problem and the way it has been defined, I can't increase the number of events from 40,000 to any larger number though I agree that I can delete some nonevent population. Please tell me your views on this, specifically:   1. Does accuracy of logistic regression depend on event rate or is there any minimum event rate which is recommended ?   2. Is there any special technique for low event rate data ?   3. Would deleting my nonevent population would be good for the accuracy of my model ? I am new to statistical modeling so forgive my ignorance and please address any associated issues that I could think about. Thanks,