title: In non-negative matrix factorization, are the coefficients of features comparable?

text: I'm using Alternating Nonnegative Least Squares Matrix Factorization Using Projected Gradient. The result (I use 2 as rank) is like this:               Original (Matrix V)     [[ 10  20   1   9]      [  1   2   1  10]      [  7   6   1   5]      [  5   4   1   5]      [100 200   1  90]]          Basis matrix (Matrix W)     [[ 0.0211552   0.09923254]      [ 0.032219    0.00224676]      [ 0.0259072   0.02522768]      [ 0.02491072  0.01388348]      [ 0.19794036  1.00108155]]          Coef (Matrix H)     [[ 144.00101426   66.1643417    62.06074785  170.96878852]      [  69.14663274  187.0342459     0.           54.62548645]]      The Coef seems comparable. For example, in the first column of Coef, the 144 is greater than 69. Does that mean in the first column of original data, the Feature1 is more important than Feature2? I'm not sure about the previous question. Because the result might also be like this:               Basis matrix (Matrix W)     [[ 0.211552   0.09923254]      [ 0.32219    0.00224676]      [ 0.259072   0.02522768]      [ 0.2491072  0.01388348]      [ 1.9794036  1.00108155]]     Coef (Matrix H)     [[ 14.400101426   6.61643417    6.206074785  17.096878852]      [  69.14663274  187.0342459     0.           54.62548645]]      The error function (something like V-WH) might be the same. However, in this situation, in the first column of Coef, Feature2 becomes more important. I think there might be some constraints for `Matrix W` and `Matrix H` to make one of them in the same order of magnitudes so that the `Coef` is comparable. I checked the paper about Alternating Nonnegative Least Squares Matrix Factorization Using Projected Gradient, but I can't say I understand the algorithm. Does anyone have ideas about this?


title: Why arrange variables by causality in bivariate regression?

text: Suppose we have variables $(X,Y)$ and we have theory tell us that $X$ $\overset{\text{cause}}{\implies} Y$. Perhaps they're time-series variables and it would be common to see something like this: $$\boxed{Y_{t+1} = a + b X_{t} + e_{t+1}.\,\,\,\,\,(1)}$$ My question is why it's not just as legitimate to specify: $$\boxed{X_t = \alpha + \beta Y_{t+1}+\varepsilon_{t+1}\,\,\,\,\,(2)}$$ where $\alpha = \frac{-a}b$, $\varepsilon_t = \frac{-e_t}b$, $\beta = \frac1b$ ? I understand that in finite samples these equations usually don't hold (for the estimates). However this doesn't answer the question of why the former is more legitimate. Also, the application isn't forecasting, so I don't need the ease of being able to plug in $X_t$ into $(1)$ and arriving at $\text{forecast}(Y_{t+1})$ with no effort. This question occurred to me when I was reading an empirical paper. They had variables $(A,B,Y)$ and they thought that $A,B \overset{\text{cause}}{\implies} Y$. Specifically they thought a good approximation of the relationship was: $$\boxed{Y_{t+1} = a + b(A_t - B_t) + e_{t+1}.\,\,\,\,\,(3)}$$ However they had theoretical reasons to think that there is a break in $b$ when $(A_t - B_t)$ was "large" versus a "small" $(A_t - B_t)$. They then went on a big spree of testing all sorts of different dummy variable combos. They even went on to specify a non-linear smooth transition model. However, in this spot it seems far more sensible, interpretable, easier (and probably more statistically powerful) to just do: $$\boxed{Q_\tau(A_t-B_t|Y_{t+1}) = \alpha(\tau) + \beta(\tau)Y_{t+1} + \varepsilon_{t+1}(\tau)\,\,\,\,\,(4)}$$ where $Q$ is the quantile function. Is this wrong?


title: Probability distribution of $\sigma$

text: Suppose $n = 5$, $s = 0.0771$. The question is: what's the $P(\sigma < 7.63\%)$?


title: Recommendation for a book about recommender systems

text: Can you recommend a book with good information that can be applied to developing a recommender system?


title: Optimal polynomial order in equality-constrained linear regression

text: In fitting a polynomial $p_N = \sum\limits_i^N a_i x^i$ to my dataset $(x,y)$ using a least squares criteria, I'm interested in determining the optimal polynomial order $N$. The caveat is that, owing to problem specific restrictions, the polynomial must satisfy the following boundary conditions: $p_N(0) = 0$, and $\frac{d}{dx}p_N(L)=0$. ($L$ is known). I gather that if the problem was unconstrained, the order could be determined using 2 successive ANOVA F-tests: the first to determine whether or not order $N$ adequately explains the variance in the observed $Y$, and if it passes, a second test to determine whether the added explained variance by order $N$ over order $N-1$ is statistically significant. Is there a similar test that I could apply that accounts for the presence of boundary conditions? Thank you in advance for your time!


title: How to sequence z-score creation and log transformation in an experiment with multiple groups and conditions?

text: In a psycholinguistic task, participants listened to and viewed stimuli, and were asked to make acceptability judgements on them:   * 4 conditions    * 4 groups    * Rating scale from 1-5 I have been advised to use z scores and log transformation (for R) on the ratings scores: ### Questions:   * Should the ratings be computed into Z scores before log transformation?   * Should computation (whether z or log first) be done by group or for the whole data set?   * Should computation (whether z or log first) be done by condition or for the whole data set?


title: ARIMA loss of effectiveness due to non-Gaussianity

text: I am looking for references regarding to loss of effectiveness of ARIMA due to different misspecification. For instance the following article:   * Chung Chen and George C. Tiao Journal of Business & Economic Statistics Vol. 8, No. 1 (Jan., 1990), pp. 83-97 tells what happens if there is random level shifts in data, but they are ignored and approximated by original ARIMA model. The following:   * Johannes Ledolter, International Journal of Forecasting, Vol. 5, No. 2, 1989, pp. 231-240 tells what happens if there is random outlier ignored. What I am looking for (and failed to find) is what happens when the data is not Gaussian. Imagine, that Box-Cox transformation has been performed not carefully and we use not optimal arguments. Or, even with proper Box-Cox transformation, it is not guaranteed that our initial data were distributed by Box-Cox power law exactly, then our transformed data are not normal exactly. What will be estimators of forecast variance increase (decrease?) and bias increase in this case?


title: Linear regression with non-linear constrained coefficients

text: I have a simple 3-factor linear model: $y = a_1x_1 + a_2x_2 + a_3x_3$, and I want the coefficients to satisfy $a_3 = a_1 \times a_2$. How should I code it in R?


title: holdback validation - test and validation equal but training much better - acceptable?

text: I have a dataset that, after modelling (with bootstrap aggregation of trees), I find I have a very high and thus over-fitted training result, but my held- back validation and test portions are practically identical, and are of a level I am happy with. I am not quite sure if this is a valid model in this case. I'm using a $60:20:20$ split. I was hoping that due to the test and validation portion being similar, the model can be assumed reliable for the portion of the test data explained. Please correct me if I am wrong.


title: Paneldata regression where some variables have common unit roots

text: I've got paneldata on Y, A, B, C. And I am trying to model dependent variable Y on independent variables A, B, C. First I checked for unit roots (in Eviews). Y, A, B, C have no individual unit roots, but Y has a common unit root, as well as A. How should I continue? I already regressed Y on A, B, C but now I am hesitant since I'm not sure if I may do this, since Y and A have common unit roots. Any help would be appreciated. T=1900 (trading days) N=90 (French stocks) Y=liquidity measure A=stock price. B=no. Shares C=some value/growth stock score


title: Multiple weights in multilevel models

text: I am currently working on a random intercept multilevel model using the European Social Survey round 6 dataset. It is a 2-level model with individuals (level 1) nested within countries (level 2). To simplify things, imagine the following regression: $Y_{ij}=\beta_{0j}+\beta_{1j}X_{ij}+e_{ij}$ where the dependent variable is **trust in the European Parliament** on a scale from 0-10, and the level-1 predictor is **gender**. In the data set, there exists two weights:   * _Design weight_ : The design weights are inclusion probabilities for individuals $i$ in countries $j$. The design weight corrects for slightly different probabilities of selection, thereby making the sample more representative of a ‘true’ sample of individuals from each country.   * _Population size weights_ : The population size weight makes an adjustment to ensure that each country is represented in proportion to its population size. The population size weight is calculated as PWEIGHT= [Population size]/[(Net sample size in data file)*10 000] My question is: do I need to specify the population size weights when I run the multilevel model? I tend to get different results. Below is the regression with design weights apllied (I am using Stata):               . xtmixed trstep gndr [pw = dweight]|| land:, mle var           Obtaining starting values by EM:           Performing gradient-based optimization:           Iteration 0:   log pseudolikelihood =  -92442,22       Iteration 1:   log pseudolikelihood =  -92442,22  (backed up)          Computing standard errors:          Mixed-effects regression                        Number of obs      =     39923     Group variable: land                            Number of groups   =        24                                                          Obs per group: min =       579                                                                    avg =    1663,5                                                                    max =      2711                                                               Wald chi2(1)       =      5,91     Log pseudolikelihood =  -92442,22               Prob > chi2        =    0,0151                                            (Std. Err. adjusted for 24 clusters in land)     ------------------------------------------------------------------------------                  |               Robust           trstep |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]     -------------+----------------------------------------------------------------             gndr |   ,1147821   ,0472334     2,43   0,015     ,0222063    ,2073578            _cons |   4,144926    ,117911    35,15   0,000     3,913825    4,376027     ------------------------------------------------------------------------------          ------------------------------------------------------------------------------                                  |               Robust                  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]     -----------------------------+------------------------------------------------     land: Identity               |                       var(_cons) |   ,3184852   ,0689119      ,2084065    ,4867066     -----------------------------+------------------------------------------------                    var(Residual) |    5,93535   ,2514202      5,462477    6,449158     ------------------------------------------------------------------------------      And here is the regression with both design weights, population size weights and scaling applied:               . xtmixed trstep gndr [pw = dweight]|| land:, mle var pweight(pweight) pwscale(size)          Obtaining starting values by EM:           Performing gradient-based optimization:           Iteration 0:   log pseudolikelihood = -81334,099       Iteration 1:   log pseudolikelihood =  -81333,24       Iteration 2:   log pseudolikelihood =  -81333,24            Computing standard errors:          Mixed-effects regression                        Number of obs      =     39923     Group variable: land                            Number of groups   =        24                                                          Obs per group: min =       579                                                                    avg =    1663,5                                                                    max =      2711                                                               Wald chi2(1)       =     10,73     Log pseudolikelihood =  -81333,24               Prob > chi2        =    0,0011                                            (Std. Err. adjusted for 24 clusters in land)     ------------------------------------------------------------------------------                  |               Robust           trstep |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]     -------------+----------------------------------------------------------------             gndr |   ,1680609   ,0513105     3,28   0,001     ,0674942    ,2686276            _cons |   3,745146   ,1854299    20,20   0,000      3,38171    4,108582     ------------------------------------------------------------------------------          ------------------------------------------------------------------------------                                  |               Robust                  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]     -----------------------------+------------------------------------------------     land: Identity               |                       var(_cons) |   ,2645594   ,0598105      ,1698583     ,412059     -----------------------------+------------------------------------------------                    var(Residual) |   6,070198   ,3532338      5,415894    6,803549     ------------------------------------------------------------------------------      I can't figure out how the population weights influence the ML estimates. The official ESS documentation says: > When comparing data from two or more countries but without reference to the > average (or combined total) of those countries, only the design weight need > be applied. When comparing data of two or more countries and with reference > to the average (or combined total) of those countries, both design and > population size weights should be applied. (ESS Documentation) The question is, whether I actually just compare countries, or compare them to an overall mean. The latter seems to be correct, as the parameter estimates actually relate to the overall mean $\mu_{00}$. But I may be wrong. The Stata manual has an extensive section on weighting here but it did not clear things up for me. Any suggestions?


title: Curve based clustering of multivariate data (time series like data)

text: > **Possible Duplicate:**   >  Reducing no of variables subsetted based on depth for PCA I have a question, I am trying to apply a method to my research area, which has not been aplied yet, based on publications and literature I could find. The case is the following: I have a dataset with x variables and these variables are changing with depth for every observation. I will show a simple dataset: ![sampledata](http://i.stack.imgur.com/328vV.jpg) Sorry I had no idea how to enter a table so I eneterd a pic. It is important how the curve looks like, so it would be nice to cluster the observations based on their curves(characteristic), which reoresents the x variables. My questions are: What would be the most appropriate method? Can I reduce the number of curves with PCA? Can I combine the curves of the x variables to one? Or is there a method, which can cluster my observations based on their x curves? I'd prefer an R or Python solution, because I have no acccess to matlab. Thank You in advance -v


title: How would one calculate the percent error on a vote?

text: Maybe I'm not asking for the right measure. Basically, there are 9million votes already in, and 4.8million no, and 4.2million yes. There are a million votes left to count. What are the chances that the true outcome of the vote was really a yes? (800,000 would have to be yes votes of the remaining million).


title: How do I apply the method of moments for estimating parameters in a sum-relationship?

text: We have a model relationship between three random variables like this: $$ U = C + S $$ I have a ton of measurements of realizations of $U$, as well as a ton of realizations of $C$. But the interesting thing is the density of $S$. The measurements are not paired, so I can not take the differences $U-C$. One solution we tried was called _deconvolution_ , which does exactly what we want, i.e. extract the density of a "noise" variable $S$. This approach brought us fine results. Alternatively, we wanted to try a parametric approach, specifically the method of moments. **My first question:** Is the following algorithm correct, assuming gamma distributed $U$, $C$ and $S$? This should exactly reflect the method of moments, if I understood it correctly:   1. Fit a gamma distribution to $U$ and $C$, and store the parameters $k_U$, $k_C$, $\theta_U$ and $\theta_C$, using the wikipedia notation.   2. Compute the expected value and variance of $U$ and $C$, using $\mu = k\theta$ and $\sigma^2 = k\theta^2$. Denote them as $\mu_U$, $\sigma^2_U$, $\mu_C$ and $\sigma^2_C$.   3. Assuming independence, compute $\mu_S = \mu_U - \mu_C$ and $\sigma^2_S = \sigma^2_U - \sigma^2_C$.   4. Solve the equation system $\mu_S=k_S\theta_S$ and $\sigma_S^2 = k_S\theta_S^2$ for $\mu_S$ and $\theta_S$. Then I have the two parameters, and thus the gamma density of $S$. But the kernel density estimates of $U$ and $C$ are very different from the fitted gamma densities. So, I'd like to try this approach with Weibull distribution assumptions. This leads to **my second question:** How would I proceed when I assume Weibull distributed $U$, $C$ and $S$? The formula for the mean and variance look very ugly over at Wikipedia. Are there analytical solutions or do I have to go R and do some numerical optimizing there? Finally, if you are still there, **my third question:** Deconvolution fails as soon as $C$ and $S$ are correlated. With the method of moments, I think I can just see how the distribution of $S$ would change by assuming a correlation of, say, 0.5. I would do this by altering only the step 3. in my first question, and computing $\sigma^2_S = \sigma^2_U - \sigma^2_C - 2 Cov(U,C)$. Is that a valid approach?


title: Which statistical method to use in survey-based quantitative research?

text: I hope you can help me because at the moment I am confused and helpless. I am final year PhD student, and my supervisors asked me to find out statistical methods other than exploratory (percentage of findings) to present my data. Here are my research data where we want to present consumer attitude and reasons towards disposal, with individual responses (N=70) collected via a questionnaire with the following items:   1. Product is functional? yes/no/not answered   2. If it is working can you pinpoint the failure point? yes/no/not answered give the failure point   3. If it is working, did you buy new product? yes/no/not answered   4. If it is not working did you buy new products? yes/no/not answered?   5. Are you aware of other options? yes/no/not answered I want to know how can I measure reliability, and what other statistics analysis I can present: GLM, ANOVA , MANOVA, etc.


title: How do you create a multivariate distribution with both continuous and discrete data?

text: I understand how to create a multivariate normal distribution to handle multiple sources of continuous data, and I understand how to create a multivariate categorical distribution to handle the discrete data. But I can't find any information on distributions over both continuous and discrete data at the same time. Can someone point me to a book/paper/blog that describes how to do this?


title: Reliabilities or factorial analysis first?

text: I've computed a scale with several theoretically defined factors. Now I am doing a factorial extraction to see if EFA discovers the same dimensions as in theory. My question is:   1. Do I need to check internal reliabilities of the theoretical scales before doing the factorial analysis, that is take the items that are supposed to belong to a given scale and check the internal consistency of this scale? Then based on reliability coefficients, do I eliminate certain items that lower scale reliability and with the left items conduct the factorial analysis?    2. Do I first extract the factors, try to determine the scales, whether they are consistent with the theoretical defined ones or not, and then take the items found to comprise a scale with factor analysis and try to determine their internal consistency? In other words, do I first check scales reliability and then conduct EFA, or first conduct EFA then check scale reliablity, after having eliminated some items?


title: Mixed model specification in R

text: I have nested structure design with the following model:               Yijklt = m + ai + bij + rijk + sikl + (bs)ijkl + eijkl           where i = 1, 2, ..................., s; number of sets      j = 1, 2, ........................, r; number of replications      k = 1, 2, ........................., m; males      l = 1,2,............................, f; females       t = 1, 2, ................, n; progeny           where m is the general mean           ai is the effect of the ith set           bij is effect of jth replication in the ith set           rik is the effect of kth male in the ith set           sikl is the the effect of the lth female to kth male in the ith set           (bs)ijkl is the interaction term      The "non-mixed" model formula in R is as follows:                model <- lm(y ~ set + replication %in% set + male %in%                  set + female %in% male %in% set + replication %in%                  female %in% male %in% set)      model <- lm(y ~ set + replication:set + male:set + female:male:set + replication:                 female:male:set)      What would be its corresponding model in mixed model using lme4 or nlme package in R? I want to do it right. Please note that all terms in the new mixed model are random, I just want to estimate variance components associated. Edits: work out example based on the answers below by @Aaron               library(lme4)     set.seed(5)     d <- expand.grid(set=factor(1:3), male=factor(1:3),female=factor(1:3),      progeny = factor(1:3), replication=factor(1:4))     d$y <- rnorm(nrow(d))


title: Wilcoxon rank sum test for significant differences between two time points

text: I have data like so …(Col1:Companyname Col2:Data at timepointA; Col3:replicate Data at timepointA; Col4:Data at timepointB; Col5:replicate Data at time point B               GeneA      0    0    0    0     GeneB   2666 2254 1833 1849     GeneC      0    0    0    0     GeneD      0    0    0    0     GeneE      0    0    0    0     GeneF      0    0    0    0     GeneI      0    0    0    0     GeneH   2007 2062 1751 1767     ….(3000 more rows) ….      I am trying to implement a non-parametric method the Wilcoxon rank sum test (WRST) to see which Genes are significantly different between two broad groups(timepointA and timepointB). end up with another column(Pval) after using this approach. From the theory I have seen WRST is used to determine broader differences and returns a single p-value is it appropriate for more focused questions like are individual genes significantly different (from two broad groups)? would the equation look like this?   Wg = Rg1 − (n1(n1 + 1))/2 where g is the gene of interest (row from data) n1 is the sample size of that group and Rg1 is the rank sums of that group. I would end up with the same Wg all the time for every gene! help/clarification is appreciated.


title: How to classify country names given possible alternate spellings or abbreviations?

text: Let's say I have a list of users who have specified the country they reside in by typing in something. I want to find the total number of users who came from the US, the UK, and everywhere else. But, since users typed in their country manually my data looks like this:               USER   COUNTRY     1      USA     2      United States     3      US     4      UK     5      England     ...       Are there any public datasets or good algorithms to solve this type of problem?


title: Maximum likelihood estimation using kalman filter for an economical two factor model, whats the measurement noise?

text: I'm trying to implement a Kalman Filter for the parameter estimation of a linear gaussian two factor model in Matlab. (Schwartz Smith model for commodity prices) In other words: I try to compute the log-likelihood of the parameters. **My model:** $X_t = A X_{t-1} + \epsilon_X$ , with X beeing two-dimensional. $Y_t = \begin{pmatrix} 1 \\\ 1 \end{pmatrix}^T X_t + \epsilon_Y$ , with Y beeing one-dimensional. $A$ is time invariant and only depends on the parameters $\theta$, which I would like to determine. **My question:** Whats the variance of $\epsilon_Y$ ? I know, that it normally represents the noise of the measurement process, but I don't know what would be the equivalent in an economic context? I managed to implement the filter and the results for experimenting with the variance are (mean = 0):   * For high variance: High log-likelihood in unreasonable areas far away from the real values.    * For very (!) small variance: Reasonable results close to real values, but extremly small _log_ -likelihoods around $-10^{12}$. I fear, that this will cause numerical issues in the upcoming maximum likelihood estimation. (I plan to use Metropolis Hastings) Any help would be appreciated! Thanks!


title: Error while using dlmModReg

text: I have started coding a Kalman filter in R. I am using dlmModreg to build an object of class dlm, which I am planning to use as my input to dlmFilter. But, I am stuck in the dlmModreg step itself. I am providing the parameters for dlmModreg as described in 1. This is my code               linear.reg     <- lm(y.learn$hospitalizations ~ y.learn$google);     flu.model      <- dlmModReg(y$google,dV=25.66,dW=0.001016);      I get the following error               Error in dlmModReg(y.learn$google, dV = 25.66, dW = 0.001016) :      Inconsistent dimensions of arguments      For dV I used the residual variance obtained from the linear model(linear.reg) However, this error vanishes when I remove dW. Also if I change my code as               dlmModReg(y.learn$google, dV = 25.66, m0 = 0.023711, c0 = 0.001016)      I get the error below               Error in dlmModReg(y.learn$google, dV = 25.66, m0 = 0.023711, c0 = 0.001016) :      unused argument(s) (c0 = 0.001016)      Here m0 is the Beta(google) estimate and c0 is the Beta(google) variance. I am not sure what is the right value to be used for parameters dW, dV, m0, c0?


title: Putting limits on estimated coefficient values

text: I recently came across an example of using eviews to put non-linear/convex limits on estimated coeficients on OLS. I could track down a link to a forum describing how this is done (here)... But I don't actually know anything more about how these models are estimated in practice. Has anyone a link to a paper/book/source code/explanation on how this is _implemented_ in practice? # Question update: Using cardinal's answer to my orginal question i have been able to solve this problem using $\verb+optim()+$ --just replacing $\beta$ by $1-\exp(\beta)$ in the minus log likelihood equation and minimizing itteratively. Now i would like to solve it using $\verb+nls()+$. So for example, in the following case:               library(MASS)     #coefficients bounded by 1     n<-100     p<-2     x<-mvrnorm(n,rep(0,p),diag(p))     t1<-4     t2<-1.5     t3<--2     y<-cbind(1,x)%*%c(t1,t2,t3)     y<-y+rnorm(n,0,2)      i get $1-\exp(\hat{t}_2)=1$ and $1-\exp(\hat{t}_3)=-2$. So my question now is can this be done with $\verb+nls()+$?


title: Finding x percentile in a graph

text: I have the following set of data:               Time   Counts   Total Counts   % Complete     1      2              2           2/10     2      0              2            2/10     3      4              6            6/10     4      3             9            9/10     5      1             10          10/10      In other words, after 1 minute went by, there were 2 counts; after 2 minutes, there were 0 counts, etc. I can easily bar graph this in Excel or R. But I'm having trouble figuring out how to easily find the point at which the "Counts" are X% complete without summing up the total counts after each additional row and dividing by the total count in the fourth column. Is there a way to get quickly get the row in which the data is X% complete in Excel or R?


title: Proper use and interpretation of zero-inflated gamma models

text: Hello StackExchange users, **Background:** I am a biostatistician presently wrestling with a dataset of cellular expression rates. The study exposed a host of cells, collected in _groups_ from various donors, to certain peptides. Cells either express certain biomarkers in response, or they don't. The response rates are then recorded for each donor-group. Response rates (expressed as percentages) are the outcome of interest, and peptide exposure is the predictor. Note that observations are clustered within donors. Since I only have the summary data, I am treating the donor-wise response rates as continuous data (at least for now). The complication arises from the fact that I have many zeroes in my data. Far too many to be ignored. I am considering a zero-inflated gamma model to deal with the fact that I have skewed continuous data coupled with an overabundance of zeroes. I have also considered the Tobit model, but this seems inferior since it assumes censoring at a lower bound, as opposed to genuine zeroes (econometricians might say the distinction is moot). **Question:** Generally speaking, when is it appropriate to use a zero- inflated gamma model? That is, what are the assumptions? And how does one interpret its inferences? I would be grateful for links to papers that discuss this, if you have any. I've found a link on SAS-L in which Dale McLerran provides NLMIXED code for a zero-inflated gamma model, so it appears to be possible. Nonetheless, I would hate to charge forth blindly. Very much appreciated,   Brenden


title: Validating a model for a set of DNA sequences

text: Consider a set of DNA sequences i.e. strings composed from the four letters A, C, G, T. We model these sequences as a random vector, $\mathbf{X}=(X_1,\dots, X_n)$. Each member of the set of sequences is a realization of this random vector. In this case the $X_i$ only take the values A, C, G, T. More generally let $r$ be the number of letters. We define a family of probability distributions for X. Our aim is to find the probability distribution that is the best fit for the data. Our assumption is that the variables $X_i$ can be divided into disjoint subsets, such that variables in each subset are dependent among themselves, and the subsets are mutually independent. It suffices to define the distribution on one of these subsets. This distribution is defined as follows. Let $\mathbf{S}=(X_{j_1},\dots, X_{j_k})$ be a sub-vector of $\mathbf{X}$ for some subset $\\{j_1,\dots, j_k\\} \in \\{1,2\dots, n\\}$. Then $\mathbf{S}$ takes $r^k$ possible values. Then we define a categorical probability distribution on $\mathbf{S}$. $P(S=(s_{1}, s_{2}, \dots, s_{k})) = p_{s_{1}s_{2}\dots s_{k}}$ So the distribution of $\mathbf{S}$ is defined by $r^k$ probabilities. Observe that the collection of random variables $\\{X_{j_1},\dots, X_{j_k}\\}$ with this distribution is not independent in general. We use MDL to find a suitable model for the data, in combination with a genetic search algorithm. This gives good results when applied to simulated data, with the caveat that it requires a sufficiently large data set to converge to the correct model. However, the important question is how well this model validates against real data. I believe this can be divided into the question of how much explanatory power this model has when applied to real data, as well as what predictive power it has. I thought about the explanatory power issue, and I'm not sure what can be usefully done. Maybe this model just has too many parameters from the perspective of model fitting. However, if anyone has suggestions, please let me know. As regards prediction, I think I can use cross-validation as well as permutation resampling. For cross-validation, I can divide the data into K portions, and set aside in turn each of the portions to use for prediction, and select the model with the other K-1 portions. So, the question is then how one should do prediction. Since the aim here is to check whether the predicted model for the validation set is the correct one. Usually, one needs a decision rule, which (I think) in this case would be of the form likelihood > value in this case, since the higher the likelihood, the more likely the data comes from the model. I'm not sure how to formalize this, or if one can justify it from some theory. I've found precisely zero references of anyone doing this, which suggests I may be doing something wrong. I've looked at some decision theory stuff but most results seem to assume conditions are met that are not true in this case. I'm also not sure how to set the cutuff. One possible approach would be to use the empirical distribution of the training set and set some cutoff based on that. A related alternative approach would be to generate ROC curves, and then generate some summary statistics from that, like Area Under Curves. For permulation resampling, I can take the data set and randomly permute it across the strings. So the different values of $X_i$ in the different strings will be exchanged among themselves. This will produce a data set which has the same marginal distributions as the original data set. but where the individual variables are uncorrelated. One can then check to see how well the model distinguishes this permuted data set with the original data set. It seems to me this if there is clear separation, i.e. the likelihoods of the permuted data set are well below the likelihoods of the original data set, then that is evidence that the original data does have correlation structure. I'm not sure how to quantify this separation.


title: Is it ever okay to ignore heteroskedastic residuals and continue with analysis?

text: My data is misbehaving and I can't seem to get residuals with constant variance despite doing more transformations than Optimus Prime. Is it ever okay to just continue with analysis in and just make a note that some assumptions were violated as the data is difficult or impossible to fit with a linear model? ![enter image description here](http://i.stack.imgur.com/48BCJ.png)


title: Fisher's Exact Test in contingency tables larger than 2x2

text: I was taught to only apply Fisher's Exact Test in contingency tables that were 2x2. Questions:   1. Did Fisher himself ever envision this test to be used in tables larger than 2x2 (I am aware of the tale of him devising the test while trying to guess whether an old woman could tell if milk was added to tea or tea was added to milk)   2. Stata allows me to use Fisher's exact test to any contingency table. Is this valid?   3. Is it preferable to use FET when expected cell counts in a contingency table are < 5?


title: Estimated standard errors for weak instrumental variable

text: I would like an explanation on the statement in bold below. At first glance, I'd think that a weak instrumental variable would yield a even bigger standard error estimate. _"When instruments are weak, however, two serious problems emerge for two- stage least squares. First is a problem of bias. Even though two-stage least squares coefficient estimates are consistent — so that they almost certainly approach the true value as the sample size approaches infinity — the estimates are always biased in finite samples. When the instrumental variable is weak, this bias can be large, even in very large samples. Second, **when an instrumental variable is weak, two-stage least squares’ estimated standard errors become far too small**. Thus, when instruments are weak, confidence intervals computed for two-stage least squares estimates can be very misleading because their mid-point is biased and their width is too narrow, which undermines hypothesis tests based on two-stage least squares."_ Murray, Michael P.. Avoiding Invalid Instruments and Coping with Weak Instruments. Journal of Economic Perspectives, v. 20, pp. 111-132.


title: How to analyze residuals of Poisson log-linear model?

text: I have bird count data and use classical poisson loglinear model, i.e. we have counts `obs(i,j)` \- observed count for site `i` and year `j`, and the model is:               ln(model(i,j)) = site_effect(i) + year_effect(j)      where `model(i,j)` is expected count for site `i` and year `j`, and it is assumed that counts are poisson distributed, i.e. `obs(i,j) ~ Poiss(lambda = model(i,j))`. Serial correlation and overdispersion was taken into account. Now I need to **analyse residuals of this model as a dependent variable in another regression with some explaining variables** (Please don't tell me I should put these variables in the above model - to be brief, there are some technical issues behind that.) **The problem now is:**   1) how to define the residuals   2) what is their distribution   3) which regression use to explain them.   **Possible solutions of residual definition and analysis:**   A) obs/model with another poisson loglinear regression   B) log((obs+1)/model), normal distribution -> normal linear regression   C) log((obs+1)/(model+1)), normal distribution -> normal linear regression   D) ... any other??   **Ad Solution A** : use `obs/model`,i.e. divide observed value by the expected value from model (for each i,j). But what is their distribution? I would say Poisson, but the numbers are not whole numbers! **So I multiplied by 20 and rounded to whole numbers to be able to fit poisson. Is this a good idea, or is it nonsense? Or how should I fit poisson to real number continuous random variable?** Here is the distribution: ![distribution of 20*obs/model, rounded to whole numbers, with Poisson fit](http://i.stack.imgur.com/HzsJG.gif) Note that the Poisson distribution doesn't fit (used goodfit() test in R). So there is a problem that obs/model are real numbers - not whole natural, as the poisson distribution. Moreover, I'd like to avoid another poisson loglinear regression if possible, as I'm quite infamiliar with how to analyze the explaining variables, the explained/unexplained variability etc. as in the normal regression. **Ad Solution B/C** : use something like `log((obs+1)/model)` or `log((obs+1)/(model+1))`. There is a problem of zero observed counts (there are no zero model counts), so I can't just use `log(obs/model)` as I would like to. How to solve this problem?? It is better to use solution B or C? Here are the distributions. ![Solution B](http://i.stack.imgur.com/hjhPr.gif) ![Solution C](http://i.stack.imgur.com/pzQTX.jpg) Note that also in this case, the normal distributions don't fit. Nothing really fits, but this is quite normal with real count data. :-) **So, summary** :   \- which solution you think fits best?   \- how to overcome the problems each solution has?   \- is it a big mistake to use B/C and go for normal linear regression, if I'd like to avoid the poisson? Should I use B or C?   Thanks a lot!


title: Implementing MMAC algorithm

text: I was reading a paper by Fadi A. Thabtah and et al., "MMAC: A New Multi-class, Multi-label Associative Classification Approach". I was trying to implement their MMAC algorithm. In MMAC algorithm that on page 3 of the article. I do not understand the algorithm, can some one please help me to convert it to a pseudo code or program code, it will be really helpful, thanks! Link article: https://docs.google.com/document/d/1YuXgPrHtH33xhtqWg2T3hSS290fnWK3JBCSnQWvsboE/edit?usp=sharing


title: Ways to overcome small number of survey responses

text: For one of my academic project I ran a sample survey. But unfortunately the respondents were not very proactive & I get just 25 responses for my study. With this small response I guess I can't do any statistical studies (e.g. anova, regressions, etc). So, I was wondering if based on the sample's mean & standard deviation I can simulate a larger dataset & then do my statistical analysis. I'll also be grateful if somebody can tell me the method to do simulation with given mean & sd in R software. Regards,   Ari


title: Community detection in mobile social network

text: Can anyone guide me regarding any good survey papers and implementations of recently published papers on community detection in mobile social network. I am willing to implement any existing algorithm and then R&D my own algorithm. I would be using ONE simulator for the purpose of community detection based on real data traces from a college campus.


title: How should standard errors for mixed effects model estimates be calculated?

text: In particular, how should the standard errors of the fixed effects in a linear mixed effects model be calculated (in a frequentist sense)? I have been lead to believe that the typical estimates (${\rm Var}(\hat\beta)=(X'VX)^{-1}$), such as those presented in Laird and Ware [1982] will give SE's that are underestimated in size because the estimated variance components are treated as though they are the true values. I have noticed that the SE's produced by the `lme` and `summary` functions in the `nlme` package for R are not simply equal to the square root of the diagonals of the variance-covariance matrix given above. How are they calculated? I am also under the impression that Bayesians use inverse gamma priors for the estimation of variance components. Do these give the same results (in the right setting) as `lme`?


title: Larger coefficients (economically) or flipped sign when using fixed effects instead of OLS?

text: Is it possible to get a larger coefficient (either a larger negative or a larger positive) when moving from OLS to a fixed effects regression? Furthermore, is it possible/likely for a coefficient sign to flip when moving from OLS to fixed effect regression? Could this be completely 'normal' or is it a sign of problems with the model / data?


title: Beta distribution from mean and quantile

text: Suppose I'm given the mean and one quantile (e.g. the 20% quantile) of a random variable $x$, and I want to find the parameters $\alpha$ and $\beta$ of a Beta distribution that has the same mean and quantile. Is there an efficient way to do it? Using just the mean, I know that since $\bar{x} = \frac{\alpha}{\alpha+\beta}$, we have $\beta = \frac{\alpha}{\bar{x}} - \alpha$. So we only really have one parameter to estimate. But I'm unsure how to use the quantile information to take the next step. Maybe there's something I can do with the Incomplete Beta when I know the ratio $\frac{\beta}{\alpha} = \frac{1-\bar{x}}{\bar{x}}$? I have access to `R` myself, so I could use a numerical optimizer for this, but ideally I need a method that can be carried out in Excel in someone else's environment. Excel does have `BETA.DIST()` and `BETA.INV()` functions available. A look-up table would be fine, but a closed-form formula would be better if it's possible.


title: Prioritization based on three factors

text: **Background:** Sales reps visit doctor and detail about a product/drug. One visit is termed as one call. In return he writes the prescriptions to doctors prescribing that particular drug. **Problem Statement:** There are three variables to be considered based on which an algorithm needs to be developed which will help prioritize the doctors.   1. **Potential Untapped:** Prescriptions written by the doctor about products which are analogs/competitive to the client's product.   2. **Potential Tapped:** Prescriptions written by the doctor about the client product.   3. **Call Value:** It is the ratio of total calls given to a doctor to the total prescriptions he has written. E.g.: A doctor who was given 10 calls, and who in turn wrote 10 prescriptions is less important than a doctor who was given 5 calls, and who in turn wrote 10 prescriptions... The aim is to use these 3 variable and come up with an algorithm that will help to prioritize the doctors available.


title: Analyzing the behavior of a non-linear model

text: I am working with a random forest model and would like to better understand how the variables impact the forest output. Using partial plots is useful, but this is basically a linear treatment of the model, because it assumes there is no interaction between the variables. Are there known methods for evaluating the non-linear-ness of the a random forest model or for analyzing how variables interact with each other within the model?


title: What is the expected dot product of two evolving vectors?

text: Suppose we have the following random process. We start with two vectors $a_1=(0)$ and $b_1=(0)$. In going from $i$ to $i+1$, we will a perturbation to $a_i$ and $b_i$. With probability $p$, we perform Case 1, otherwise perform Case 2.   * Case 1: We pick an element of $a_i$, say at coordinate $x$ (chosen uniformly at random from all coordinates $1,2,\ldots,|a_i|$, where $|a_i|$ is the length of $a_i$), and append it to the end of $a_i$ to form $a_{i+1}$. Similarly, to form $b_{i+1}$ append the element at coordinate $x$ from $b_i$ to the end of $b_i$. The remainder of $a_{i+1}$ and $b_{i+1}$ is the same as $a_i$ and $b_i$, respectively. [The copied element is at coordinate $x$ in both $a_i$ and $b_i$.]   * Case 2: Choose $x$ and $y$ to be two coordinates chosen uniformly at random from $1,2,\ldots,|a_i|$. If $x \neq y$ then set $a_{i+1}[x]=a_i[x]+1$ and $b_{i+1}[y]=b_i[y]+1$ (if $x=y$ then don't do anything [it's not clear at this point whether or not I want to enforce $x \neq y$]). Again, the remainder of $a_{i+1}$ and $b_{i+1}$ is the same as $a_i$ and $b_i$, respectively. > **Question** : What is the expected dot product of $a_i$ and $b_i$? That is, > what is $\mathrm{E}(\sum_x a_i[x] \cdot b_i[x])$? So here's an example of what these vectors could look like:               i a_i b_i     1 (0) (0)     2 (0,0) (0,0) [case 1 x=1]     3 (0,0,0) (0,0,0) [case 1 x=1]     4 (1,0,0) (0,1,0) [case 2 x=1 y=2]     5 (1,0,0,1) (0,1,0,0) [case 1 x=1]     6 (1,0,0,1) (0,1,0,0) [case 2 x=4 y=4]     7 (1,0,0,2) (0,1,1,0) [case 2 x=4 y=3]      I'm looking at a related problem studying evolving random graphs. In the random graph problem, the two vectors represent the in-degrees and out-degrees of an evolving network over time. In the graph, we can duplicate vertices or add or delete an edge (which is related to cases 1 and 2 above). In this case, the dot product therefore counts the number of directed paths of length 3. However, the above problem differs from the one I'm considering since (a) it's not necessarily the case that $\sum_x a_i[x]=\sum_x b_i[x]$ and (b) there's nothing to stop parallel edges arising here. Although, I'm hoping that this simplified question will be more answerable and the techniques can be re-used.


title: Multivariate time series forecasting using Support Vector Regression in R

text: I have data like this:               pce        pop        psavert    uempmed    unemploy     507.8      198712     9.8        4.5        2944     510.9      198911     9.8        4.7        2945     516.7      199113     9.8        4.6        2958     513.3      199311     9.8        4.9        3143     518.5      199498     9.8        4.7        3066      I am trying to use SVM - regression to fit the data like this               svmRbftune <- train(unemploy ~ pce + pop + psavert + uempmed,                         data = EconomicsTrain, method = "svmRadial",                         tunelength = 14, trControl = trainControl(method = "cv"))          svmRbfPredict <- predict(svmRbftune, EconomicsTest)      How can I forecast two to three steps ahead?


title: What is a good way of estimating the dependence of an output variable on the input parameters?

text: I am evaluating a scenario's output parameter's dependence on three parameters: A, B and C. For this, I am conducting the following experiments:   * Fix A+B, Vary C - Total four sets of (A+B) each having 4 variations of C   * Fix B+C, Vary A - Total four sets of (B+C) each having 3 variations of C   * Fix C+A, Vary B - Total four sets of (C+A) each having 6 variations of C The output of any simulation is the value of a variable over time. For instance, A could be the area, B could be the velocity and C could be the number of vehicles. The output variable I am observing is the number of car crashes over time. I am trying to determine which parameter(s) dominate the outcome of the experiment. By dominate, I mean that sometimes, the outcomes just does not change when one of the parameters change but when some other parameter is changed even by a small amount, a large change in the output is observed. I need to capture this effect and output some analysis from which I can understand the dependence of the output on the input parameters. A friend suggested Sensitivity Analysis but am not sure if there are simpler ways of doing it. Can someone please help me with a good (possibly easy because I don't have a Stats background) technique? It would be great if all this can be done in R. **Update:** I used linear regression to obtain the following:               lm(formula = T ~ A + S + V)          Residuals:          Min       1Q   Median       3Q      Max      -0.35928 -0.06842 -0.00698  0.05591  0.42844           Coefficients:                 Estimate Std. Error t value Pr(>|t|)         (Intercept) -0.01606    0.16437  -0.098 0.923391         A            0.80199    0.15792   5.078 0.000112 ***     S           -0.27440    0.13160  -2.085 0.053441 .       V           -0.31898    0.14889  -2.142 0.047892 *       ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 0.1665 on 16 degrees of freedom     Multiple R-squared: 0.6563, Adjusted R-squared: 0.5919      F-statistic: 10.18 on 3 and 16 DF,  p-value: 0.0005416       Does this mean that the output depends mostly on A and less on V?


title: Linear discriminant analysis and the "kernel trick"?

text: This is problem 12.10 in "The Elements of Statistical Learning": > Suppose you wish to carry out a linear discriminant analysis (two classes) > using a vector of transformations of the input variables $h(x)$. Since > $h(x)$ is high-dimensional, you will use a regularized within-class > covariance matrix $W_h + \gamma I$. Show that the model can be estimated > using only the inner products $K(x_i, x_{i'}) = \left < h(x_i), h(x_{i'}) > \right >$. How can I go about showing that regularized linear discriminant analysis can be estimated using only inner products, as in the "kernel trick" that is often used with SVM's?


title: What is "the shortest half of the data"?

text: Here is a histogram (realized with JMP) displaying two types of box plot called outlier box plot and quantile box plot. ![enter image description here](http://i.stack.imgur.com/I7MVb.png) Right below, there are a bunch of explanations of the meaning of the different features of the outlier box plot. ![enter image description here](http://i.stack.imgur.com/0DKbI.png) The red bracket is "the shortest half of the data (the densest region)". What does this expression exactly mean? How different is it from the quartiles? Source


title: How to efficiently repeat a function on a data set in R?

text: Let's say I have a function `funToRepeat`:               funToRepeat <- function(df){         rows <- nrow(df)         cols <- ncol(df)         err <- runif( rows * cols )         dfNew <- df + as.data.frame( matrix(err, nrow=rows, ncol=cols) )         dfNew     }      I want to operate `funToRepeat` on `dfTest` (a toy data set), multiple times, say `n`.               dfTest <-       structure(list(V1 = 1:2, V2 = 3:4), .Names = c("V1", "V2"),        row.names = c(NA,-2L), class = "data.frame")                                        How would I do this and average the `n` outputs of `funToRepeat` in `R`, efficiently? A very bad example with `n=10` may be:               dfAvg <- as.data.frame(list(V1=c(0, 0), V2=c(0, 0)))     for(i in 1:10){         dfAvg <- dfAvg + funToRepeat(dfTest)     }     dfAvg <- dfAvg/10      Please notice, I am not trying to bootstrap. I want to operate `funToRepeat` on my data set and take the average of the `n` outputs. I want to do this as I am adding `err`, a random term, to my data.frame and want to provide a reasonable answer. (I am not doing exactly this in real life, but something very similar)


title: Iteratively reweighted least squares : asymmetric weights

text: For robust m-estimation, all the convergence results I'm aware of assume symmetric weights (eg: Huber function) in their formulation of the iterative reweighted least squares algorithm. Does the method apply to asymmetric weights, and does the corresponding convergence results hold? Thanks


title: Multiple two-way interactions in simple slope analysis

text: Currently, I am working on a multiple regression model that includes a total of 9 variables: Two education dummy code variables (categorical) as a control, one target independent variable (continuous), three moderator variables (continuous), and three interaction terms (that combine moderators with primary independent variable). The model appears as so: Block 1: Dummy Code Education 1, Dummy Code Education 2   Block 2: Target independent variable, three moderator variables   Block 3: Three interaction terms Some quick notes:   \- The three moderators are subscale scores from a single measure, but show no issue with collinearity   \- One of the education dummy codes is significant throughout each model   \- Out of the three interaction terms, only one of them is significant. What I am trying to do is to analyze the significant interaction using simple slope analysis, but trying to account for the other two interaction terms. The equation would look something like this: y = b0 + b1(IV) + b2(Moderator) + b3(Sig.Interaction) + b4(Moderator 2) + b5(Moderator 3) + b6(Interaction2) + b7(Interaction 3) + error It is easy enough for me to create a simple slope analysis with just holding everything constant except the IV, moderator, and sig interaction using a -SD/Mean/+SD values. However, the lines on the plot would changed due to the influences of the other factors. My question is what values should I enter into the other factors in order to account for the changes in the lines? Because if I just enter a 1 or 0, it just simply me holding everything constant again.


title: Regression for categorical y and proportional and discrete xs

text: I have a model with a categorical dependent variable and independent variables that contain a proportional variable also. It appears that multinomial would be an ideal fit, but am not certain on how to model an x variable containing proportional data in mlogit. Is there any other regression type that would be applicable in this case ? Thanks in advance.


title: How should you express a negative binomial distribution in an exponential family form?

text: When $X$ is $X_1$,...$X_n$, how do you express the following negative binomial distribution in an exponential family form? $$ f(k;r,p)\equiv\text{Pr}(X=k)=\binom{k+r-1}{k}(1-p)^{r}p^k~~~~~\text{for}~k=0,1,2,\ldots $$


title: SVD on a transposed matrix

text: I struggle with the interpretation of the singular values obtained by the SVD, also in relation to what the principal components of the PCA tell me. The mathematical background has been described here. **1\. Question** I wondered if it matters when I transpose my matrix and then apply SVD/PCA. Therefore I chose a 2x2 matrix to keep it simple and plot-able (columns are plotted against each other using R):               x = matrix(c(1,0,2,3),nrow=2)     x          [,1] [,2]     [1,]    1    2     [2,]    0    3     plot(x)      ![plot\(x\)](http://i60.tinypic.com/jrsoqv.jpg) Applying SVD and PCA:               svd(x)$d     [1] 3.6502815 0.8218544          prcomp(x)     Standard deviations:     [1] 1 0      PCA does perfectly what I expect, one dimension exhibits no variance (since we have only two data points) while the other explains all variance. I thought the SVs of the SVD represent something similar, so why is the second SV > 0? **2\. Question** Now transposing the matrix (or: plotting the rows of x against each other)               y = t(x)     y          [,1] [,2]     [1,]    1    0     [2,]    2    3     plot(y)      ![plot\(y\)](http://i57.tinypic.com/2n1axcx.jpg)               svd(y)$d     [1] 3.6502815 0.8218544          prcomp(y)     Standard deviations:     [1] 2.236068 0.000000      While PCA does what I expect (one dimension exhibits all the variance and this time more since the data points are further away from each other in that dimension), the SV from the SVD are the same as above. First question still holds here, but additionally: Why didn't the singular values change? (The mathematical answer is that X*X and XX* have the same eigenvalues.) This also holds for datasets that do not have a square matrix and also higher dimensional data sets:               > m = cbind(1:10+rnorm(10),1:10+rnorm(10))     > prcomp(m)     Standard deviations:     [1] 4.5043248 0.7969853          > svd(m)$d     [1] 27.196875  3.064609          > prcomp(t(m))     Standard deviations:     [1] 3.163586e+00 1.618444e-15          > svd(t(m))$d     [1] 27.196875  3.064609      (I am aware that I didn't scale the data as it is suggested.)


title: multi-channel time series

text: I'm trying to classify data, which is represented by multi-channel time series: one observation contains 300 time series each of length 100. For my algorithm I'd like to get less time series of the same length per observation, so it seems that I need a special kind of PCA. I found something about empirical orthogonal functions and multi-channel SSA, but I'm not sure if it's what I need (it looks quite complicated for me). Could you please recommend some papers on this topic or at least the names of methods, which can do this. Thank you very much!


title: Detecting parts of a song

text: Hopefully this is not too subjective... I'm looking for some direction in efforts to detect the different "parts" of a song, regardless of musical style. I have no idea where to look, but trusting in the power of the other StackOverflow sites, I figured someone here could help point the direction. In most basic terms, one could detect different parts of a song by just grouping consecutive repeating patterns and calling those a "part". That's maybe not so hard -- computers are pretty good at detecting repetition in a signal, even when there's some small variation. But it's hard when the "parts" overlap, as they do in most music. It's hard to say what kinds of music would be most well-suited to this kind of system. I would guess that most classical-style symphonic music would be easiest to process. > Any ideas of where to look for research in this area?


title: Clustering probability distributions - methods & metrics?

text: I have some data points, each containing 5 vectors of agglomerated discrete results, each vector's results generated by a different distribution, (the specific kind of which I am not sure, my best guess is Weibull, with shape parameter varying somewhere around exponential to power law (1 to 0, roughly).) I am looking to use a clustering algorithm like K-Means to put each data point into groups based on the attributes of its 5 component distributions. I was wondering if there are any established distance metrics that would be elegant for these purposes. I have had three ideas so far, but I'm not a seasoned statistician (more of a beginning data-mining computer scientist) so I have little idea how far I am off track.   1. Since I don't know exactly what kind of distributions I'm dealing with, my brute-force approach to the problem was to chop each of the distributions (I have 5 per point) into each of its respective discrete data values (I pad each corresponding one to the same length with zeros at the end) and use each of these values as a separate attribute for the data point itself. I tried using both Manhattan distance and Euclidean distance as metrics based on these attributes, for both the PDF and CDF.   2. Again, since I don't know what kinds of distributions I have, I figured that if I was going to measure the distance between the overall distributions I could use some sort of non-parametric test pairwise between distributions, such as the KS-test, to find the likelihood that the given distributions were generated by different PDFs. I thought that my first option (above) using the Manhattan distance would be a sort of upper bound on what I might get using this approach (since the KS statistic is the max absolute value of the difference of the CDFs, where Manhattan distance is the sum of the absolute values of the differences in the PDFs). I then considered combining the different KS-Statistics or P-values within each data point, probably using Euclidean distance, but possibly just taking the max of all of these values.   3. Lastly, in an effort to use what little I can interpret about the shape of the distributions, I thought I might try estimating the parameters of the distributions as fit into a Weibull curve. I could then cluster the distributions based on differences in the two parameters of the Weibull distribution, lambda and k (scale and shape), probably normalized according to the variance of these parameters or something of the sort. This is the only case where I thought I might have an idea of how to normalize the parameters. So my question is, what measure/methods would you recommend for clustering of distributions? Am I even on the right track with any of these? Is K-Means even a good algorithm to use? Edit: Clarification of data. Each data point (each object `Obj` that I want to cluster) actually literally contains `5 vectors` of data. I know there are exactly 5 phases that these objects can be in. We'll say (for the purposes of simplification) that each vector is of `length N`. Each one of these vectors (call it `vector i`) is a probability distribution with integer `x-values` of 1 through N, where each corresponding y-value represents the probability of measuring `value x` in `phase i` of the object `Obj`. N is then the maximum x-value I expect to measure in any phase of the object (this is not actually a fixed number in my analysis). I determine these probabilities in the following manner:   1. I take a single `Obj` and put it in `phase i` for `k trials`, taking a measurement at each trial. Each measurement is a single whole number. I do this for each of 5 phases of a single object, and in turn for each object. My raw measurement data for a single object might look like: Vector 1. [90, 42, 30, 9, 3, 4, 0, 1, 0, 0, 1] Vector 2. [150, 16, 5, 0, 1, 0, 0, 0, 0, 0, 0] ... Vector 5. [16, ... ..., 0]   2. Then I normalize each of the vectors on its own, with respect to the total number of measurements in that given vector. This gives me a probability distribution in that vector, where each corresponding y-value represents the probability of measuring `value x` in `phase i`.


title: The $z$-test vs the $\chi^2$-test for comparing the odds of catching a cold in 2 groups

text: I just read in a rather well-respected (popular) science magazine (the German PM, 02/2013, p.36) about an interesting experiment (without a source, unfortunately). It caught my attention because intuitively I doubted the significance of the result, but the information provided was sufficient for reproducing the statistical testing. The researchers wondered whether getting cold in cold weather increases the odds of catching a cold. So they randomly split a group of 180 students into two groups. One group had to hold their feet into cold water for 20 minutes. The other kept their shoes on. Kind of a funny manipulation, I think, but on the other hand I am not a doctor and maybe doctors think funny. Ethical issues aside. Anyways, after 5 days, 13 of the students in the treatment group had a cold, but only 5 in the group that kept their shoes on. The odds ratio of this experiment thus is 2.87. Given the rather small sample size, I started wondering if this difference may be significant. So I conducted two tests. First a simple test of equality of proportions using the normal approximation. This test has $z=1.988$ with $p=0.0468$. My guess is that this is what the researchers tested. This is truely just significant. However this z-test is only valid in large samples, if I am not mistaken, due to the normal approximation. Furthermore, the prevalence rates are rather small and I wonder whether this may not affect the coverage rate of the confidence interval of the effect. So my second try was a chi-square test of independence, both with Monte-Carlo simulation and standard Pearson Chi-square. Here I find p-values both about $p=.082$. Now that's all not so reassuring about the results. I wondered if there are more options to test this data and what your thoughts on the two tests are (in particular the assumptions of the first, significant, test)


title: How to evaluate/select cross validation method?

text: How can I decide which cross validation method is appropriate for my problem and data type? For instance, choosing among leave-one-out, or K-fold (and which K is appropriate?). Most of my searches end up on selection of model rather than selection of cross validation method.


title: Looking at residuals vs. residual percentages

text: Suppose I fit a linear regression to some data (say, weight vs. height), and all the standard linear regression assumptions are satisfied (in particular, the data is homoscedastic). For example, here's a random figure pulled from amstat.org that looks like it satisfies what I'm thinking of: ![enter image description here](http://i.stack.imgur.com/kvDMa.jpg) Now I'm doing some exploratory data analysis, so I want to look at examples where the linear regression is particularly off; that is, I want to sort all the individuals by how badly the regression predicts their weight from their height (so that, say, I can look at further details, like whether people whose weight the model underpredicts tend to eat a lot of junk food). My question is:   * Should I sort all the individuals by their raw residuals?   * Or should I sort all the individuals by the residuals as a percentage of the prediction, i.e., by `residual weight / predicted weight`? On the one hand, it seems like sorting by raw residuals might be the way to go, since standard linear regression errors are based off the squared residuals, and not the residual percentages. On the other hand, someone who weighs 70kg when their predicted weight is 50kg seems much more of an outlier than someone who weighs 120kg when their predicted weight is 100kg. Is it just a matter of preference or the particular model at hand?


title: Time series feature selection: 20 independent variable, 1 dependent variable, 13 observations

text: I understand the challenges associated with the small # of observations. Nevertheless, some of the predictors have very powerful signals. The question is, how does one use R to most efficiently code for building a model with the best combination of predictors. By best, my criteria is lowest out of sample MAPE for the last 3-4 observations in my time series. Is there an efficient way to build this procedure/test in R? Beyond code, any other insights to tackling this particular problem? Thank you


title: Where can one get real (health) actuarial/insurance probability data/distributions?

text: I'm researching some problems in the health/medical insurance sector and would like to use actual probability distributions (or current data) for doing some sort of risk analysis using simple simulations and the like, to see how different types of policies stack up against different people's perception of risk. I guess the only type of data I'd need IMO is the probability of fatality for various types of injuries (assuming the 'types' are quite granular). For example: 0.05% chance of a fatal injury, 25% chance of fracture etc., Basically some real data on which I can base my analysis to gauge the risk attitudes of buyers of the policy. Any ideas/pointers would be greatly appreciated. If the data is country specific I'm preferably targeting North America, South America and Europe.


title: How can one show a Kmeans solution is unique?

text: Suppose we are given a distribution P and a constant K. We wish to minimize the kmeans objective w.r.t centers ${C1,..Ck}$: ![Kmeans Objective](http://i.stack.imgur.com/Sw1lL.png) What constraints on $P$ are known to imply that the optimal solution is unique (up to reordering of centers)? Examples of ad hoc methods would also be appreciated. In particular intereset to me are continuous distributions. An example for one-dimensional $P(x)$ and $k=2$ with a unique optimal solution: $P(X=1)=0.5$ $P(X=-1)=0.5$ The only optimal centers are $C1=1, C2=-1$


title: A hair dresser's conundrum

text: My hairdresser Stacey always puts on a happy face, but is often stressed about managing her time. Today Stacey was overdue for my appointment and very apologetic. While getting my haircut I wondered: How long should her standard appointments be? (if the customer's preference for clean round numbers could be ignored, for a moment). Something to consider is a certain 'ripple effect' where one very late customer can lead to a string of delayed appointments. In reality, hair- dressers intuitively learn to space appointments longer and longer as they fear of these stressful days. But an optimum, elegant solution must be achievable by some statistical genius out there.. (if we dumb down reality a little) Let's assume a) hair cutting times are normally distributed and b) there's only one hair dresser. The cost of setting appointments _too long_ is obviously the hairdresser's time wasted waiting for the next appointment. Let's cost this wasted time $1 per minute. But if the appointment's not long enough, the next customer is kept waiting, which is a heavier cost of $3 per minute to customer-loving Stacey.   * Stacey works up to 8 hours per day, and has enough demand that she can fill as many appointments as she can fit in   * The mean hair cut takes her 30 minutes, with a std. dev of 10 minutes. (let's also assume men's cuts and women's cuts are the same!)  EDIT - some have rightly pointed out that Stacey could attend to EARLY customers ahead of their appointed time. This adds another layer of complexity, but if we treat this as quite a realistic problem we need to include it. Let's forget my 90/10 assumption and try for an assumption perhaps a little closer to reality.   * Some customers are late and some are early. The mean of customers are 2 minutes late with a standard deviation of 2 minutes (sounds reasonably near reality no?)  Exactly how long should her appointments be? * * * @alexplanation sorry I've moved the goal posts on you! I'm sure R readers appreciate your answer.


title: How to obtain simultaneous confidence intervals for correlated frequencies?

text: I observed 400 episodes of nursing care in a hospital. I tracked the movement of the nurses between 5 rooms $A,B,C,D$ and $E$. The maximum likelihood estimates of the probabilities of moving from room $i\rightarrow j$ are given by: \begin{equation} P_{ij}=\displaystyle \dfrac{\text{# of times from room $i\rightarrow j$}}{\displaystyle \text{Total # of transitions to any room}}\end{equation}   * Is there a way of defining a confidence interval on this maximum likelihood estimate $P_{ij}$?    * And for all maximum likelihood estimates of all possible room combinations ($n=25$ combinations)? This is confusing me because the $P_{ij}$ where $i=A..E$ and $j=A..E$ are dependent on each other. How can I account for this? **Reference** : I have come across this reference: http://arxiv.org/pdf/0905.4131v1.pdf This suggests that for n observations $X_i$, the empirical maximum likelihood estimate $\hat P_{ij}$ minus the actual probability $P_{ij}$ would tend to a normal distribution with mean 0 and standard deviation $\epsilon$. $$\sqrt{n}|\hat{P_{ij}}-P_{ij}|\sim N(0,\epsilon)\quad \text{as}\quad n\rightarrow \infty$$ How do I calculate $\epsilon$?


title: Vector multiplication in BUGS and JAGS

text: In R, c(3,1,0) * c(2,0,1) == c(6,0,0). This is not dot product and it's not cross product. First, what is the name for this product, and second, does it work in WinBUGS, OpenBUGS and/or JAGS?


title: Multidimensional Scaling: Interpreting output of different distance matrices (Euclidean or correlation)

text: I would like to understand the difference between using a Euclidean distance matrix or correlation matrix as input to a nMDS algorithm. I have completed MDS plots of both, and while similar, the distances between points are different (in some cases radically so). I understand that both measure different concepts of similarity, however what insights into the data do I get from using one method over the other? Thanks


title: Is median better than mean for reporting on website speed

text: Whilst I understand why `median` can be a better metric than `mean` in certain circumstances ( _e.g. you're creating/marketing a product, you want to use the`median` as it represents the biggest age group within your population and therefore real potential buyers, whereas `mean` age might be a number which doesn't correspond to anybody, or a very few people_) I am wondering whether `median` is a good metric when it comes to measuring website performance. Let's look at the below data representing `page load time` (in sec, the smaller the better): ![enter image description here](http://i.stack.imgur.com/fOZKX.png) `Sample 1` and `Sample 2` are identical, except for `C4`, meaning the website is only faster for 1 person. Yet, the `median` changes significantly between the 2 samples (nearly 40% improvement) whereas the `mean` is much more stable (less than 10% difference). Accordingly my take on this is that `mean` is more accurate are representing the page load time. **Do you agree? / Is my observation biased because of the particular samples I have chosen or would`mean` be more representative than `median` no matter the samples?**


title: Difficulty to intrepret pvclust results - only many low level tree clusters appear as significant

text: I'm trying to assess the uncertainty in hierarchical cluster analysis. It is a dataset composed of 409 observations and 27 variables (with a value ranging form 0 to 100). The dataset represents immunohistochemical scores in a gastrointestinal cancers. A meaninful clustering of observations and markers is observed with Pearson uncentered distance and average linkage.               hc <- hclust(Dist(t(imputedMatrix), method="pearson"), method="average")     hr <- hclust(Dist(imputedMatrix, method="pearson"), method="average")     heatmap.2(imputedMatrix, Rowv=as.dendrogram(hr), Colv=as.dendrogram(hc),   col=greenred(100), scale="none", ColSideColors=patientcolors, density.info="none", trace="none")          pv2 <- pvclust(imputedMatrix, method.dist="uncentered", method.hclust="average", nboot=10000)     plot(pv2, hang=-1)     pvrect(pv2, alpha=0.95)          clsig <- unlist(pvpick(pv2, alpha=0.90, pv="au", type="geq", max.only=TRUE)$clusters)          dend_colored <- dendrapply(as.dendrogram(pv2$hclust), dendroCol, keys=clsig,     xPar="edgePar", bgr="black", fgr="red", pch=20)     heatmap.2(imputedMatrix, Rowv=as.dendrogram(hr), Colv=dend_colored, col=greenred(1      However, when using pvclust to assess their uncertainty, many small low-levels subclusters are highlighted as significant, but not any higher level one. Also, a group of tumours (to the right in the plots) is indeed a control group that should be clearly distinguished at the highest level from the other ones. pvclust shows even there the same kind of pattern. Indeed the clusters of interest are the magenta, green and red in the color bar. Does this pvclust results support their existance (versus a clustering artefact by chance)? How could these pvclust results be interpreted? Maybe I am using the tool in a wrong way? Or the wrong tool for this kind of data? Thank you very much in advance. ![pvclust dendrogram](http://i.stack.imgur.com/nedQr.png) ![hclust and pvclust, the latter in red in the dendrogram](http://i.stack.imgur.com/KoaHW.png)


title: How to implement hierarchical clustering in $O(N^2)$ instead on $O(N^3)$

text: First a theoretical question. I know that natively, an hierarchical clustering algorithm is of complexity on the cube of number of samples N. This is due to the fact that in each iteration, one has to go over the entire distance matrix to find the smallest value. But, it is possible to implement it in lower complexity (of N sqr). How it is done and can you reffer me to an article or description of the implementation? Finally, I used scikit library to do HClustering. It was fine but it is limited to 15,000 samples or so. I want to cluster much more. can someone refer me to an implementation that fits larger numbers (and preferably runs at sqr N)?


title: Categorize statistical tests into univariate and multivariate methods

text: I am not sure about the following tests/methods whether they belong to the category of univariate or multivariate tests. Univariate tests/methods: t-test,ANOVA,ANCOVA, univariate linear regression (Y=a+bX) Multivariate tests/methods: MANOVA, multivariate linear regression (Y=a+bX+cX+…), cluster analysis, partial least squares discriminant analysis (PLS-DA), principal component analysis Can anyone confirm whether these tests/methods are assigned to the right categories?


title: Limits to tree-based ensemble methods in small n, large p problems?

text: Long time grazer, first time poster. I'm hoping to gather people's opinion on the following: Tree-based ensemble methods such as Random Forest, and subsequent derivatives (e.g., conditional forest), all purport to be useful in so-called "small _n_ , large _p_ " problems, for identifying relative variable importance. Indeed, this appears to be the case, but my question is how far can this ability be taken? Can one have, say 30 observations and 100 variables? What is the breaking point to such an approach, and are there any decent rules of thumb that exist? Preferably ones backed by evidence, using either simulated or real data sets? I haven't found much on the latter (here and here), so your thoughts/advice/(on topic) reference suggestions are most welcome! Cheers!


title: Markov chain prediction with unknown states

text: I want to make a prediction of the next state based on a training and test dataset. So I split my data into a train and test set and calculate the MLE on the train set and want to predict the next state on the test dataset. The problem now is, that there can be states in the test set, that never have been learned via the train set. How would I cope with that? A similar problem occurs, if I have a similar train/test setup, and want to calculate the log-likelihood on the test dataset, with states never occuring in the train set. Hope someone can help me out with that! Regards, Philipp


title: Spatial autocorrelation versus spatial stationarity

text: Let's assume we have points in two-dimensional space, and we wish to measure the effects of attributes $X$ on attribute $y$. The typical linear regression model is of course $$y= X\beta + \epsilon$$ There are two problems here: the first is that the $\epsilon$ terms may be spatially correlated (violating the independent and identical errors assumption), and the second is that the regression slope may vary throughout space. The first issue can be dealt with by incorporating spatial lag terms into the model, as in $$y=\rho W y + X\beta + \epsilon$$ We can even incorporate spatially autoregressive omitted variables (spatial fixed effects) with the spatial Durbin model described in the text by LeSage and Pace $$y=\rho W y + X\beta + WX\lambda + \epsilon$$ where $\rho$ is the strength of spatial correlation controlled by the weights matrix $W$. Clearly the form of the spatial lag will depend on assumptions about the form of the spatial correlation. The second problem has been addressed using "geographically-weighted regression" (GWR), a technique I am not as familiar with, but which is explained by Brunsdon, et al. (1998). As far as I can tell, it involves fitting an array of regression models to weighted sub-regions, thus getting an estimate of each $\beta_i$ that changes based on its space, $$\hat{\beta}_i = (X^TW_iX)^{-1}X^T W_i y$$ where $W$ is another spatial weights matrix, not necessarily different from the one above. **My question** : Is the first method (spatial autoregression) not sufficient to yield an unbiased estimate of the average marginal effect of $X$ on $y$? GWR seems like it is over-fitting: of course the $\beta$ changes in space, but if we want to know the average expected effect of a treatment without regard for its spatial position, what could GWR contribute? Here's my attempt at an initial answer:   1. If I want to know the premium for an additional bedroom _in a specific neighborhood_ , it seems that GWR would be my best option.    2. If I want to know the unbiased _global_ average premium for an additional bedroom, I should use spatial autoregressive techniques. Would love to hear other perspectives.


title: How to compare the means of two dependent groups when the groups are not paired?

text: I have performed a behavioral experiment on bees from one single hive. I let them fly under two different conditions and now I want to compare the average flight speed in the two conditions. The number of bees are different in the two conditions. In condition 1, I have 17 bees, and in condition 2, I have 19 bees. I don't know if the same individuals were flying in both conditions, because I didn't mark the individuals. However, it is most likely that at least a few of the bees were flying in both conditions. Moreover, the same bee might have been flying more than once in each condition. Because of this, I guess I have to treat my data as dependent (repeated measures), right? Now, my question is what test I should use to see if the two group means differ from each other? Since the two groups are dependent, I can't perform an independent samples t-test, but since my data is not paired, I can't perform a paired samples t-test either... Any kind of advice would be much appreciated! Thank you in advance, Therese Edit 1: The hives are estimated to consist of about 200 bees, but I can't say for sure how many bees this hive consisted of. Edit 2: First of all, thank you for your comments, Placidia. The two conditions were performed randomly over a couple of days, so I assume that the bees don't learn anything from the different conditions. However, I am wondering about the finite sample correction that you mentioned. In another similar experiment, I sampled 40 bees for each condition (but for each condition, I sampled the bees during two time periods, in the morning and in the afternoon, and did this over 2 days. In other words, I sampled 10 bees at a time). Do I have to make an adjustment for finite sample in that case? In yet another experiment, I sampled 17 bees from a wild population. Does that mean that the sample is infinite, and that I don't have to adjust for a finite sample? Thank you again, I really appreciate your help!


title: Two stage GMM estimator in Matlab

text: I am trying to create a simple GMM estimator for the mean of a normally distributed random variable using the first three odd central moments of a normal distribution (all of which should be zero theoretically). In two stage GMM, normally the first step is to minimise a least squares cost function of the errors of each individual moment condition within the sample; arriving at an initial estimate of the mean mu. The vector of errors after this first stage (evaluated at the first parameter estimate) is then used to create a covariance matrix of the moment conditions from which a weighting matrix is derived; then used in the second stage to arrive at an efficient estimate of the mean. I have coded this up in Matlab, and to the best of my knowledge, this has been done correctly. The issue I am having is that the covariance matrix is very close to singular - meaning that it doesn't have an inverse. The 1st stage appears to work ok - the issue is with the second stage not working as required due to this issue. Any help with diagnosing why I am running into this issue would be most appreciated! Here is the Matlab code.               clear; close all; clc;          %% Generate normal data mu = 10; var = 4; N = 10;          X = mu + sqrt(var)*randn(N,1);          %% Function definitions for basic GMM f_cost_c = @(mu_hat,degree) (1/N)*sum((X-mu_hat).^degree);          % Specify the number of odd moments to include in cost fn N_moments = 3;          f_cost_v = @(mu_hat)f_cost_c(mu_hat,1); for i = 1:N_moments-1 f_cost_v = @(mu_hat)[f_cost_v(mu_hat); f_cost_c(mu_hat,2*i+1)]; end          f_L_cost_c = @(mu_hat) f_cost_v(mu_hat)'*f_cost_v(mu_hat);          %% Find the value of mu_hat which minimises the basic squared cost (the 1st stage of 2S GMM) mu_hat_0 = 5; mu_hat = fminunc(f_L_cost_c,mu_hat_0); mu_hat_1 = mu_hat;          %% Use the value of the cost vector at the value found in 1st stage to generate weighting matrix v_cost = f_cost_v(mu_hat);          S_hat = (1/N_moments)*(v_cost*v_cost'); W_hat = inv(S_hat);          %% 2nd stage of 2SLS GMM          % Define new cost using estimated weight matrix f_L2_cost_c = @(mu_hat) f_cost_v(mu_hat)'*W_hat*f_cost_v(mu_hat);          mu_hat_0 = mu_hat; mu_hat = fminunc(f_L2_cost_c,mu_hat_0);


title: Confusion related to kriging

text: I was going through the wiki article related to kriging http://en.wikipedia.org/wiki/Kriging. However, I couldn't follow some derivations. ![enter image description here](http://i.stack.imgur.com/hIx4y.png) In the first figure for simple kriging, how come the variance of the error became the expression like that ![enter image description here](http://i.stack.imgur.com/48T0x.png) In the third figure, how come the variance of the error came to be expression like that. Can anyone please kindly give the derivation. I am really confused. In the above expression how come the


title: Hidden markov models: output observations defined by a (non-hidden) markov model?

text: Let me explain what my goal is: I would like to define a hidden markov model with two hidden states and say, five possible observations. As I understand (I'm quite new to HMMs), in each state HMM will output one of the observations, based on the given output probabilites. I would like to extend this behaviour in such a way that the output is not dependant only on that probabilities, but also on the last observation that was output by the HMM. You could think of it as if each hidden state had another (non-hidden) markov model embedded in itself, and use it to define the output. Both hidden states would have exactly the same markov models embedded, but with different transition probabilities. I hope it is clear what I'm asking. Any hint would be greatly apreciated.


title: Is it a good idea to use log scale on scree plots for PCA/ICA/FA?

text: I always found the concept of determining the "ideal" number of components/factors for an ICA/PCA/FA via a scree plot useful and quick, but also a bit shaky. In an effort to try to make the scree plot "elbow" more obvious (and without trying to fit a function to the scree plot data points). I was thinking whether plotting on log-scale coordinates might help shed some light. ![enter image description here](http://i.stack.imgur.com/sLXWk.png) The non-log plot (first from the left) is pretty bunched up, but it's as "obvious" as it gets with scree plots that the third factor is the elbow (and thus my ideal number of factors should be 3). Now, it also seems that whether my log scale is in base 2 or 10 makes little difference for the plot appearance - which I guess is good. Lastly, I find that the X/Y log plots show the most quaintness for the 3rd factor, making it the point where 2 types of "curves" in the plot meet. I also think this feature stands out better here than in the first (linear) plot. Is it safe to assume this is always the case? Is it worth recommending to use X/Y log coordinates for scree plots?


title: RNG, R, mclapply and cluster of computers

text: I'm running a simulation on R and a cluster of computers and have the following problem. On each of X computers I run:               fxT2 <- function(i) runif(10)     nessay <- 100     c(mclapply(1:nessay, fxT2), recursive=TRUE)      There are 32 computers, each with 16 cores. However, around 2% of the random numbers are identical. What strategies would you adopt to avoid this? I've been able to avoid this problem for fxT2 by setting a latency (i.e. delaying by a sec the time at which each job is send to each of the X computers). But it's seems very ad-hoc to fxt2. The problem is that in reality fxT2 is a long task involving pseudo random numbers. At the end of the process, I expect to get X*nessay reproduction of the same statistical experiment, **not** nessay reproductions. How to make sure that this is indeed the case and is there a way to check this?.


title: Rotation matrices and prior invariance for arbitrary dimensions

text: I have a question about a rotation matrix, which can be represented in 2 dimensions as: $$R_{2}(\theta)=\begin{bmatrix} \cos\theta & \sin\theta \\\ -\sin\theta & \cos\theta \end{bmatrix}$$ For some arbitrary angle $\theta$. This can be extended to an arbitrary number of dimensions by adding an identity matrix: $$R_{n}(\theta)=\begin{bmatrix} R_{2}(\theta) & 0 \\\ 0 & I_{n-2}\end{bmatrix}$$ I have found some "invariance" properties of a n-dimensional prior distribution when rotated in 2 arbitrary dimensions. My question is: can any rotation in arbitrary dimensions be represented by a sequence of 2-D rotations? It doesn't matter if the sequence is unique or not for my purposes. Or perhaps a better question is: if a prior distribution is invariant when rotated about 2 arbitrary dimensions, is it invariant when rotated about an arbitrary number of dimensions?


title: How to summarize multiple adjusted pvalue into a single measure?

text: Suppose I have a test with multiple hypotheses. For example, if I have X as a gene set (a set of symbols), and A, B, C pathways (sets of symbols), I am testing if A is enriched in X, B in X, C in X and so on. Then I am correcting the resulting pvalues for multiple hypothesis testing. Now I have several corrected pvalues, and I want to combine them to a single pvalue. Is this possible? I will add more info if needed.


title: self-organizing maps and order of samples

text: I've just started learning about self-organizing maps without much of a background in neural networks, so forgive me if this question seems trivial, but it seems that a SOM depends on which order the input training data points are fed into the SOM? So it's possible to obtain different visualizations for the same training data?


title: How are categorical variables used when fitting a decision tree in scikit-learn?

text: I am used to R, in which you can use factor(variable) to indicate a categorical variable. However, in scikit-learn, trying to pass a variable of strings causes the DecisionTreeClassifier to give an error. How are you supposed to use a categorical variable in scikit-learn? Is the only option for encoding a multi-class variable to use dummy variables?


title: detecting circadian rhythm in a time series

text: I have a sensor that can detect minute changes in distance. It produces a time series. I would like to point it at people and detect things like their sleeping pattern. How would one build a system that figures out what 'normal' means for this user and can say things like 'over the past 3 months your sleeping pattern has shifted 3 hours, this is bad'. Here is a picture of what the series looks like over the course of one day ![http://i.stack.imgur.com/9DbJy.jpg](http://i.stack.imgur.com/GJTF1.jpg) One transformation of it we call position, it tries to illustrate the position the subject holds (laying down, sitting up, etc) ![http://i.stack.imgur.com/h3G3j.png](http://i.stack.imgur.com/JvwoB.png) Any and all ideas are welcome. I've skimmed "Time Series Analysis and Its Applications: With R Examples" by Shumway and Stoffer. I think I should explore Spectral Analysis and Filtering due to the cyclical nature of the data. I've been using python/pandas and am looking into doing work with statsmodels and scikit-learn. EDIT1: Activity is the 'raw' data from the sensor. It is distance (higher values are closer to the sensor) but it's value lies in it's sensitivity. I want to detect breathing patterns if possible. As far as position, I'm not actually sure what the formula for that transformation was, I'll ask the author.


title: Basic statistics help

text: I know you guys are not allowed to give out any answers (I dont want any) I just need help in how to do it so I could practice on my own. Say I was given mean and the S.D. of x and y ONLY and also i am given the correlation coeffecient $r$. How would I find out the standard error of estimate for Y? I know the formula is $\sqrt{\frac{\sum{(y-y')^2}}{N}}$. and I am really confused on this. and also I am to find if given x how would i predict her y?


title: Constructs and multiple regression

text: I am using the decomposed theory of planned behavior to analyze predictor variables affecting teacher factors affecting the adoption of Web 2.0 technologies in the classroom. There are several questions that measure each construct including the dependent variable. Can I add each individual construct under transform data and then perform a multiple regression for each predictor variable against the dependent? I know I also need to check the scatter plots for outliers. Thanks, Berta


title: Metropolis-Hastings algorithm, using a proposal distribution other than a Gaussian in Matlab

text: I am currently working on my final year project for my mathematics degree which is based on giving an overview of the Metropolis-Hastings algorithm and some numerical examples. So far I have got some great results by using my proposal distribution as a Gaussian, and sampling from a few other distributions, however I am trying to go one step further by using a different proposal distribution. So far I have got this code (I am using Matlab), however with limited resources online about using different proposals it is hard to tell if I am close at all, as in reality I am not too sure how to attempt this, (especially as this gives no useful data output so far). It would be fantastic if someone could lend a hand if they know or forward me to some easily accessible information (I understand that I am not just asking coding advice, but Mathematics as well). So, I want to sample from a Gaussian using a proposal distribution of a Laplace, this my code so far:               n = 1000;       %%%%number of iterations          x(1) = -3;      %%%%Generate a starting point          %%%%Target distribution: Gaussian:          strg = '1/(sqrt(2*pi*(sig)))*exp(-0.5*((x - mu)/sqrt(sig)).^2)';     tnorm = inline(strg, 'x', 'mu', 'sig');          mu = 1;    %%%%Gaussian Parameters (I will be estimating these from my markov chain x)     sig = 3;               %%%%Proposal distribution: Laplace:          strg = '(1/(2*b))*exp((-1)*abs(x - mu)/b)';     laplace = inline(strg, 'x', 'b', 'mu');          b = 2;       %%%%Laplace parameter, I will be using my values for y and x(i-1) for mu               %%%%Generate markov chain by acceptance-rejection          for i = 2:n              %%%%Generate a candidate from the proposal distribution         y = laplace(randn(1), b, x(i-1));              %%%%Generate a uniform for comparison         u = rand(1);              alpha = min([1, (tnorm(y, mu, sig)*laplace(x(i-1), b, y))/(tnorm(x(i-1), mu, sig)*laplace(y, b, x(i-1)))]);                   if u <= alpha             x(i) = y;         else             x(i) = x(i-1);         end      end      If anyone can tell me if the above is completely wrong/going about it in the wrong way, or there are just a few mistakes (I am very wary about my generation of 'y' in for the for loop being completely wrong) that would be fantastic. Thanks, Tom


title: Is dependence of one classification on another a violation of the assumption of independence for ANOVA?

text: My subjects are mice, swimming in a pool for 60 s. The DV is how much time they spend in each quadrant of the pool. The performance of each subject is independent of every other subject. However, the software classifies each 0.5-s frame as being in one of the four quadrants. If the subject is in one quadrant he is by definition not in the other three quadrants. The typical way to analyze these data is with repeated-measures ANOVA, with each quadrant a measure, but this doesn't seem right to me if each of the observations in each quadrant is dependent on what happened in the other three quadrants. Am I thinking about this right or is repeated-measures ANVOVA OK in this case?


title: Random forest assumptions

text: I am kind of new to random forest so I am still struggling with some basic concepts.   In linear regression, we assume independent observations, constant variance…   * What are the basic assumptions/hypothesis we make, when we use random forest?    * What are the key differences between random forest and naive bayes in terms of model assumptions?