title: Find k of n items with least pairwise correlations

text: I have a matrix of pairwise correlations between n items. Now I want to find a subset of k items with the least correlation. Thus there are two questions:   1. Which is the appropriate measure for the correlation within that group?   2. How to find the group with the least correlation? This problem appears like a kind of inverse factor analysis to me and I'm pretty sure that there is a straight-forward solution. I think this problem actually equals the problem to remove (n-k) nodes from a complete graph so the remaining nodes are connected with minimum edge weights. What do you think? Thanks for your suggestions in advance!


title: Document classification on Twitter: LDA or something non-generative?

text: On Twitter, I would like to assign weight to users based on how similar their bios are to a topic (e.g. technology would be a topic, the training set would be the bios of users who are known to be interested in technology). I was thinking of using LDA, where I would get the theta value (document-topic association) of the new document (which would be the bio of the Twitter user I want to classify) and use it as the weight. Does this make sense? Should I rather use a discriminative algorithm since classification is supposed to be a discriminative problem?


title: Cross validation statistics

text: I have a question about the generation of statistics from cross validation schemes for performance estimation (two classes classification in my case). Before each experiment (K iterations) I randomly generate K folds from my pool of examples. Then I run the usual K-fold cross validation scheme to obtain performances for each iteration. Then I average these performances to obtain statistics about the experiment (K iterations). My problem is that sometimes, some test folds contains only examples from one class (let's say C), thus I can't compute test performances for the other class (N) for this iteration. My question is, how should I compute statistics for the whole experiment in this case (one or many iterations for which I have performances for only one class). Thank you


title: Aggregating votes from different users with some constraints and probabilities

text: Suppose we have $N$ users that can vote for a given question by "Yes" or "No". Each user have a probability which is considered as the confidence that we have in that user (e.g. a user may have a probability of 0.8 to answer the question correctly, and a probability of 0.2 to answer it wrongly). Each user requesst an amount of money as a reward for answering the question. I (as the asker of the question) have a limited budget for paying the users. Given this budget, which users should I ask to maximize the final answer to be correct (i.e. the answer obtained after aggregating the asked answers)? **Additional informations:** "what does it mean for a vote to be correct": if the question is "does 10/2 equals 5?" then "Yes" is a correct answer to this question; now if the probability associated to the user is 1, that means that the user will answer this question correctly because it have a very good knowledge of the domain of the question (mathematics); if the probability associated to the user is 0, this means that the user will answer this question wrongly; if the probability associated to the user is 0.5, then his answer may be correct ("Yes 10/2=5") with a probability of 0.5, or maybe wrong wt probability 1-0.5. What I want is to maximize the fact of getting the good final answer by aggregating the obtained answers (want to know that the answer to my question is "Yes" 10/2=5). For this, which users should I ask (given my budget, and given that each user request a given amount of money to accept to answer my question), and how to aggregate their answers in order to have the correct one. What is called "vote" is the classical aggregation of the answers, by taking the most common answer among the received ones (if there is more "Yes" answers then ok I will consider the final result as "Yes), but in my real example I have some constraints: probabilities associated to users, ammount of money requested by each user to answer my question, and the total amount of money that I have, so I need to also consider these constraints. If you want, forget about the word "vote" in the question, just tell me who to ask and how to aggregate the answers that I get.


title: Linearity in local sensitivity analysis

text: As it is know local sensitivity analysis attempt to quantify the local impact of input factors on the model, through partial derivatives: a derivative of the outputs accordingly to the inputs, when the last ones are slightly perturbed. Many methods are used to estimates this derivatives: Direct method, Green Function Method, etc. However, the input-output relationship is generally assumed to be linear (from Saltelli textbooks on the subject). Why it is assumed to be linear in a local neighborhood, even the whole system may not be linear. I guess it may have a relation with Taylor expansion expression. I am trying to understand the real reason behind this assumption. Thank you.


title: Understanding the intra-class correlation coefficient

text: Suppose we fit a simple mixed / multilevel model with no predictors, which some call a _variance components model_ : $$ y_{ij} = \beta_0 + u_j + e_{ij} $$ where $u_j$ are the level-2 residuals and $e_{ij}$ are the level-1 residuals, and we obtain estimates, $\hat{\sigma_u^2}$ and $\hat{\sigma_e^2}$ for the variance of $u_j$ and $e_{ij}$ respectively. I get it that $$ \nu = \frac{\hat{\sigma_u^2}}{\hat{\sigma_u^2} +\hat{\sigma_e^2}} $$ is the proportion of total variance in the response due to between-level-2 units, and so it makes sense that this is called the variance partition coefficient. What I don't get is why $\nu$ can also be interpreted as the correlation between response scores for 2 randomly selected observations from within the same level-2 unit (ie the intra-class correlation). It makes sense that $\nu$ should be _related_ to this correlation, but why should be be exactly equal ?


title: Definition of "optimal" instruments

text: The book I read (Davidson,MacKinnon - Econometric Theory and Methods) describes the definition of "optimal instrument variables" as the following: > Usually, and this is seen very often in other textbooks which cope with > regression models, we want to estimate ${\mathbb E}\left(y | X\right)$. > Another way of saying the same thing but in a more general sense is that we > want to estimate ${\mathbb E}\left(y | \Omega\right)$. $\Omega$ is a > information set containing all variables which determine the outcome of $y$. > > In the context of instrument variable estimation we consider the usual model > $y = X\beta + u$ with homoscedastic $u$'s. Now it is assumed that at least > one explanatory variable of X is endogenous. Hence we want to find a matrix > $W$ of instruments which minimize the asmyptotic covariance matrix of > $\beta_{IV}$. Here I get lost. It is assumed that in general $X$ can be assumed to satisfy the relation $X = \bar X + V$ where $E(V|\Omega) = 0$ and $\bar X = E(X|\Omega)$. Now $\bar X$ is supposed to be the optimal matrix of instruments but I don't get what $\bar X$ is supposed to be. Yes, it is the the expected outcome of $X$ given all the information which is in $\Omega$, but if one variable in $X$ is assumed to be endogenous and I replace this variable with a linear combination of the remaining variables which are assumed to be exogenous and that $\bar X$, I ask myself why $\bar X$ is a viable instrument. Why is this $\bar X$ not correlated with $u$ no more? I'm not looking for a technical solution nor for a mathematical proof of some sort just. I would like to understand the idea behind this approach.


title: Mixed effects count modeling

text: I am trying to predict the probability that someone's vehicle fails its emission inspection a discrete number of times using a Poisson count model (actually a negative binomial hurdle model, but the same logic applies), $$ P(y_i | X_i) = \dfrac{\exp(-e^{X\beta})(e^{X\beta})^{y_i} }{y_i!} $$ where the probability of vehicle $i$ failing emissions $y$ times is a function of different predictors. This is a straightforward problem. What is not straightforward is that households frequently own multiple vehicles, and I need to accommodate the within-household effects. In some ways I can think of this this as an unbalanced panel.               HH  |  VEH  | FAILURES |  X      A       1        0       x_a1      A       2        0       x_a2      B       1        2       x_b1      C       1        1       x_c1      C       1        0       x_c2      .......      It looks like a panel, but $y_t$ has no practical relationship to $y_{t-1}$, nor does the $X$ for different vehicles. Most panel methods I know involve taking the differences of one variable or another. Is there some other way to incorporate correlated observations that would help me with this problem? * * * Update: There are two packages for R that claim to do this; `MCMCglmm` and `glmmADMB`. The latter appears to have a simpler call syntax, as the former requires conversational ability in Bayesian inference. Unfortunately, I cannot get the example scripts for `glmmADMB` to run on my system (OSX 10.7, R 2.15.2), because they crash with a matrix type error.


title: Issues on computing Pearson correlation coefficient for two vectors

text: When computing the Pearson correlation coefficient between a vector $X = \{x_1,...,x_n\}$ and a vector $Y = \{y_1,...,y_n\}$, we need to compute it as $$\frac{1}{n}\sum_{i=1}^{i=n}(x_i-\bar{x})(y_i-\bar{y}).$$ Here $\bar{x} = (1/n)*\sum_{i=1}^{i=n}x_i$. My question concerns the case where $x_i$'s, $i=1\dots,n$, are in different order of magnitude. For instance, when $x_1 = 1000$, and $x_2 = 1$ is there any specific procedure to preprocess these different $x_i$'s? Or can we directly compute $\bar{x}$ as the above? Is there any article or book discussing this topic?


title: Smoothing 2D data

text: The data consist of optical spectra (light intensity against frequency) recorded at varying times. The points were acquired on a regular grid in x (time) , y (frequency). In order to analyse the time evolution at specific frequencies (a fast rise, followed by an exponential decay), I would like to remove some of the noise present in the data. This noise, for a fixed frequency, can probably be modelled as random with gaussian distribution. At a fixed time, however, the data shows a different kind of noise, with large spurious spikes and fast oscillations (+ random gaussian noise). As far as I can imagine the noise along the two axes should be uncorrelated as it has different physical origins. What would be a reasonable procedure to smooth the data? The goal is not to distort the data, but remove "obvious" noisy artefacts. (and can over- smoothing be tuned/quantified?) I don't know if smoothing along one direction independently of the other makes sense, or if it's better to smooth in 2D. I've read things about 2D kernel density estimate, 2D polynomial / spline interpolation, etc. but I'm not familiar with the jargon or the underlying statistical theory. I use R, for which I see many packages that seem related (MASS (kde2), fields (smooth.2d), etc.) but I cannot find much advice on which technique to apply here. I'm happy to learn more, if you have specific references to point me to (I hear MASS would be a good book, but perhaps too technical for a non- statistician). _Edit:_ Here's a dummy spectrogram representative of the data, with slices along the time and wavelength dimensions. ![image2d](http://i.stack.imgur.com/pk98t.png) The practical goal here is to evaluate the exponential decay rate in time for each wavelength (or bins, if too noisy).


title: How to test equal means when the observations are based on estimates?

text: We have two groups ($\mathbf{X}$ and $\mathbf{Y}$) to test the equality of the means. Suppose $$ \mathbf{X}_i \sim MVN(\boldsymbol{\mu}_\mathbf{X},\boldsymbol{\Sigma}_\mathbf{X}), \qquad \mathbf{Y}_j \sim MVN(\boldsymbol{\mu}_\mathbf{Y},\boldsymbol{\Sigma}_\mathbf{Y}) $$ where $1 \leq i \leq n_\mathbf{X}$ and $1 \leq j \leq n_\mathbf{Y}$. Here $MVN$ stands for _Multi-Variate Normal_. The means and covariance matrices are unknown, also $\boldsymbol{\Sigma}_\mathbf{X} \ne \boldsymbol{\Sigma}_\mathbf{Y}$. However, $\mathbf{X}_i$ and $\mathbf{Y}_j$ can only be estimated by $\mathbf{X}_i^*$ and $\mathbf{Y}_j^*$ respectively and $$ \mathbf{X}_i^* \sim MVN(\mathbf{X}_i, \boldsymbol{\Sigma}_{\mathbf{X},i}), \qquad \mathbf{Y}_j^* \sim MVN(\mathbf{Y}_j, \boldsymbol{\Sigma}_{\mathbf{Y},j}) $$ where $\boldsymbol{\Sigma}_{\mathbf{X},i}$ and $\boldsymbol{\Sigma}_{\mathbf{Y},j}$ are known. How can we test whether $\boldsymbol{\mu}_\mathbf{X}=\boldsymbol{\mu}_\mathbf{Y}$ or not?


title: Finding a rejection zone without an alpha being given

text: **Question:** In one faculty the distribution of the heights of the students is N(180,25) In the second faculty it's N(160,20) We take some students from both faculties and put them together in a class. There are 2 hypothesis regarding how they are distributed: $H_0$: $\frac 14$ of them are from faculty I, and $\frac 34$ are from faculty II. $H_1$: $\frac 14$ of them are in faculty II, and $\frac 34$ are from faculty I. If we look at a random student- If his height is bigger than 168 we reject $H_0$. What is the rejection zone? **What I did:** To try to formalize it I rephrased this in math terms: X is a Bernulli RV representing a single random student, s.t. p is the probability of a student to be in faculty I. Therefore $H_0: p=0.75$ I think the statistic S(X) is the height of a single random student, so I want to find out the prob. under $H_0$ that his height is bigger than 168. I did it this way : $0.75P(\frac {S(X)-180}{25}>\frac {168-180}{25})+0.25 P(\frac{S(X)-160}{20}<\frac{168-160}{20})=0.5994$ (using R) Is this result the alpha used to calculate the rejection zone $C_a$? (We've been told that $P_{H_0} (S(X)\in C_a)=a$ is how to do a statistic test) I would love some thorough explanation (maybe some good links too) for this(not only hints) as I feel my understanding of the whole subject is still lacking. Thanks a lot


title: How to solve $\mathrm dX(t)=X(t)^2 \mathrm dt+X(t)\mathrm dB(t)$ with condition $X(0)=1$?

text: I want to solve the stochastic differential equation $$\mathrm dX(t)=X(t)^2 \mathrm dt+X(t)\mathrm dB(t)$$ with condition $X(0)=1$.


title: Most common classifier used with PCA?

text: Which is the most used classifier applied once PCA is performed? I computed PCA on my data, and now I have the points projected on a line, I was wondering what is the preferred classifier. Is it a discriminant function for the normal density a good choice? Or is it better to use SVM?


title: Comparing 2 independent non-central t statistics

text: One estimate of the 'quality' of a portfolio of stocks is the _Sharpe ratio_ , which is defined as the mean of the returns divided by the standard deviation of the returns (modulo adjustments for risk free rate, etc). The sample Sharpe ratio is the sample mean divided by the sample standard deviation. Up to a constant factor ($\sqrt{n}$, where $n$ is the number of observations), this is distributed as a (possibly non-central) $t$-statistic. Are there known techniques for comparing ~~the mean of~~ independent variables distributed as non-central $t$-statistics? ~~Of course, there are non- parametric tests of mean, but is there something specific to the case of noncentral $t$?~~ (I'm not sure what I meant by that.) **edit** : the original question is somewhat ambiguous (well, it's actually not what I want). Is there a way to test the null hypothesis: population Sharpe ratio of $X$ equals population Sharpe ratio of $Y$, given independent collections of observations drawn from $X$ and $Y$? Here Sharpe ratio is mean divided by standard deviation. **edit** : given $n_x, n_y$ observations of $X, Y$, construct sample means, standard deviations, to get sample Sharpe ratios: $\hat{S}_x = \frac{\hat{\mu}_x}{\hat{\sigma}_x}, \hat{S}_y = \frac{\hat{\mu}_y}{\hat{\sigma}_y}$. Then $t_x = \sqrt{n_x}\hat{S}_x$, and $t_y = \sqrt{n_y}\hat{S}_y$ are distributed as non-central t-statistics with noncentrality parameters $\sqrt{n_x}S_x$ and $\sqrt{n_y}S_y$, where $S_x, S_y$ are the population Sharpe ratios of $X, Y$. Given these independent observations, I wish to test the null hypothesis $H_0: S_x = S_y$. In one form of the problem, one only has the summary statistics $n_x, n_y, \hat{\mu}_x, \hat{\mu}_y, \hat{\sigma}_x, \hat{\sigma}_y$. For large sample sizes, $t_x, t_y$ are approximately normal, I believe but the small sample size case is also of interest (funds often quote performance based on monthly returns).


title: Using information geometry to define distances and volumes…useful?

text: I came across a large body of literature which advocates using Fisher's Information metric as a natural local metric in the space of probability distributions and then integrating over it to define distances and volumes. But are these "integrated" quantities actually useful for anything? I found no theoretical justifications and very few practical applications. One is Guy Lebanon's  work where he uses "Fisher's distance" to classify documents and another one is Rodriguez' ABC of Model Selection… where "Fisher's volume" is used for model selection. Apparently, using "information volume" gives "orders of magnitude" improvement over AIC and BIC for model selection, but I haven't seen any follow up on that work. A theoretical justification might be to have a generalization bound which uses this measure of distance or volume and is better than bounds derived from MDL or asymptotic arguments, or a method relying on one of these quantities that's provably better in some reasonably practical situation, are there any results of this kind?


title: Estimating distribution parameters from few data points

text: Say I'm doing stats on the height of adults from various countries. I assume the heights of adults from one country are normally distributed, and ignore sex differences (I also ignore the fact that neighbouring countries tend to have similar populations). I have a bunch of data by country, but for some countries I have very few data points, which can lead to quite big errors on my estimate of standard distributions. Is there a way I can use the data from countries for which I have a lot of data to get better estimates - say if I notice standard deviation in those countries is always between 7 and 8.5 cm, but my dataset for Nepal (for which I have 9 samples) has a standard deviation of 9.5 cm, I should probably correct that downwards. But how? Is there a formula for this? When calculating the parameters for Nepal, shouldn't the data from the other countries allow me to have a "prior" distribution of expected means and deviations, which I would then update by taking the actual data from Nepal into account? How would I formalize this methodology? (I got this while looking for a simple reduction of the problem that prompted previous question, which didn't get an answer yet - I'm still mostly looking for good methodologies for thinking about this kind of problem).


title: Using R2WinBUGS, how to extract information from each chain?

text: I am new to Gibbs Sampling, and have been using WinBUGS, but I find that it is not well-suited towards storing/presenting results, so I have been calling it from R using the R2WinBUGS package. The data is apparently stored as a "bug" class. I converted it to coda to run diagnostics, and it displays each of the chains, but I am confused as to the $ extension for each individual chain. I cannot find any good documentation for the "coda" class (the cran instructions are not helpful). My code is below:               > bugs.sim <- bugs(data, inits, parameter, "gl.bug", n.chains = 5,                        codaPkg = TRUE, DIC = FALSE, n.iter = 5000)     > codaobject <- read.bugs(bugs.sim)      As you can see, I have 5 chains, and I would like to take the mean and standard deviation of each. How do I go about doing this? I can use the codaobject to take the Geweke diagnostic of each chain, it displays each chain as `[[i]]` ($i=1,\dots,5$). Thanks in advance. And any references to a detailed documentation for R2Winbugs, would also be greatly appreciated.


title: Removing factors from a 3-way ANOVA table

text: In a recent paper, I fitted a three-way fixed effects models. Since one of the factors wasn't significant (p > 0.1), I removed it and refitted the model with two fixed effects and an interaction. I've just had referees comments back, to quote: > That time was not a significant factor in the 3-way ANOVA is not of itself a > sufficient criterion for pooling the time factor: the standard text on this > issue, Underwood 1997, argues that the p-value for a non-significant effect > must be greater than 0.25 before treatment levels of a factor can be pooled. > The authors should give the relevant p-value here, and justify their pooling > with reference to Underwood 1997. My questions are:   1. I've never heard of the 0.25 rule. Has anyone else? I can understand not removing the factor if the p-value was close to the cut-off, but to have a "rule" seems a bit extreme.   2. This referee states that Underwood 1997 is the standard text. Is it really? I've never heard of it. What would be the standard text (does such a thing exist)? Unfortunately, I don't have access to this Underwood, 1997.   3. Any advice when responding to the referees. * * * Background: this paper was submitted to a non-statistical journal. When fitting the three-way model I checked for interaction effects.


title: Based on z-score, is it possible to compute confidence without looking at a z-table?

text: How are the values in a z-table computed? Is it possible to compute confidence without looking up a z-table?


title: Panel study is a quasi-experimental study? Quasi-experimental is the same as correlational?

text: I'm working with panel data. This panel dataset includes data from pupils of two kinds of schools:   * State schools (G1) where pupils change to secondary school when they are 12 years old.   * Private schools (G2) where pupils remain in the same school. They can follow studying in the same school when they are 12 years old.  I want to compare the change in a dependent variable between these two groups of pupils, controlling for some other independent variables. So, is this a quasi-experimental design? Maybe G1 can be seen as a treatment group and G2 as a control although there isn't a randomized assignment? Or is this an observational study? And finally, what exactly is a correlational study? Is it quasi-experimental or observational?


title: Getting started with time series in R

text: With some great help from this forum, I have been able to get up and running with some basic time series analysis in R. Right now, my needs are mostly univariate time series. Here is my question: I can read in daily data from database into a data frame. I have two columns, date which is understood by R as POSIXct and the second which is the value of interest and numeric. What is the best, most straightforward way to make this a ts object where R understands the start/end dates and represents daily observations dynamically? It seems to me that I shouldn't be required to (when coercing my object to a ts object) manually tell it the start and end dates when the data frame is already has it. For some context, I have been able to aggregate other data from daily to weekly, but find myself doing things in ways that just seem long and unnecessary considering R already understands my raw data as time. As you can tell, I am new to R and time series in R, but I figure that since R is so powerful, there probably is a pretty easy way around my issues. Any help you can provide will be very much appreciated.


title: Scaling data in machine learning algorithm

text: I am currently working on a linear classifier, which uses the statistical learning paradigm; that is, no knowledge about the distribution from which that training data are drawn from is available. The resulting optimal solution is a hyperplane of the form $\mathcal{H}: \mathbf{w}\cdot\mathbf{x}+b$, where $\mathbf{w} = (w_1,\cdots,w_n)^T$ is the normal vector to the hyperplane, $b$ is a real bias term, and $\mathbf{x}\in\mathbb{R}^n$ is a vector that belongs to the $n$-dimensional column space $\mathbb{R}^n$. The separating hyperplane, $\mathcal{H}$, is computed by solving an unconstrained, non-linear minimization problem. I have observed that the objective function (i.e., the function that has to be minimized) has different behaviour when training data take large values. For example, let $n=2$, i.e., we are in the 2D Euclidean space. Although my algorithm achieves to find a good solution when the 2D training samples lie on a neighborhood of zero (they are all inside a disk of radius, say 50), when I translate my training points by $(1000,1000)^T$, I find out that the objective function has not a single, but infinite minimizers that lie on a line. Apparrently, this is not desirable... I have thought that a possible way to address this issue would be the scaling of my data. What is your opinion? It is my intention to make a discussion on the above issue, if anyone else is interested in. It would be great to read some feedback, while I will come back soon. Thanks for any useful comment!


title: Correlation among repeated measures

text: Can someone explain what 'correlation among repeated measures' in G*power exactly means? I believe it is the test-retest reliability, but why is it important and what are the effects on the inferred sample size?


title: Probability of drawing a cursed coin not from a certain box

text: There are three boxes of gold coins. In the first box, 1 out of 3 coins are cursed, in the second 2 out of 3 are cursed, and in the third 1 out of 4 are cursed. Furthermore, the first box is ornate, and so has 3/5 chance of being picked, while the other two boxes have a 1/5 chance. The experiment consists of selecting a box, and then selecting a coin. What is the probability that someone drew a cursed coin and did not pick box 1? My approach: I recognize that there are non independent events. And to find the probability of not picking box 1 with a cursed coin I need to find the probability of choosing a cursed coin and the other two boxes then add them together. \begin{align} P(Curse \cap Box 2) &= P(Box 2)P(Curse | Box 2)\\ &= (1/5)(2/3)\\ &= 2/15 \end{align} \begin{align} P(Curse \cap Box 3) &= P(Box 3)P(Curse | Box 3)\\ &= (1/5)(1/4)\\ &= 1/20 \end{align} \begin{align} P(Curse \cap Box 2) + P (Curse \cap Box 3) &= 2/15 + 1/20\\ & = 0.1833 \end{align} Therefore, $P(Curse \cap Not Box 1) = 0.1833$. I feel like the logic is correct and would like a confirmation.


title: Feature Selection Methodology

text: I have a question regarding methodology in Machine Learning. Assume we have a classification problem, the size of training set is 100 every instance in training set is represented by 200 features. In order to eliminate meaningless features we evaluate Information Gain / Mutual Information of every feature and eliminate those with low mutual information. Then in order to evaluate the model we use 10-fold cross validation. The question is after all these manipulation whether the calculated evaluation can be used as real evaluation, in other words, are we allowed to perform all these manipulation and get reasonable model with reasonable evaluation. **Ideas** : Initially the entire process looked for me reasonable, but then I though if I am allowed to eliminate features based on the Infogain of train set or better to do it on validation set, on the other hand, in any case the trained model has access only to train set, so if features has low infogain it will have a minor weight in the model. I would appreciate if you could share your thoughts.


title: Estimated variance using linear factor models

text: This might sound like a trivial question, but I just can't get my head around it! I am trying to evaluate the performance of different risk models (linear factor models), by computing the **_bias statistic_**. It is simply a function of estimated variance of the asset you are testing for and its forward return. Let's consider the following model, $M_1$ $R$ = a + $b_1$$F_1$ + $b_2$$F_2$ + e = a + $B*F$ + e where _R_ is the return of the asset, _F_ are the factor premiums, _B_ are the factor exposures, _a_ is a constant and _e_ is the residual The estimated asset variance using $M_1$ is: $var(R)$ = $B^T$$\Sigma$$B$ + $\Delta$ where $\Sigma$ is the variance-covariance matrix of factors F, and $\Delta$ is the variance of the regression residuals. However, if I consider model $M_2$, which has completely different factors that $M_1$ $R$ = a + $b_3$$F_3$ + $b_4$$F_4$ + e and compute the estimated variance of the asset, I get the exact same value!! Here is the link to a snapshot of the data: data. The first column (after the date column) is the dependent variable, where the rest of the columns are factors $F_1$, $F_2$, $F_3$, $F_4$. Here is the code I used to compute the forecasted risk:               library(zoo)          data = read.csv(csv, stringsAsFactors = FALSE, na.strings= "NA")     data1.zoo = zoo(data[,c(2,3,4)], order.by = as.Date(data[,1], format="%Y-%m-%d"))     data2.zoo = zoo(data[,c(2,5,6)], order.by = as.Date(data[,1], format="%Y-%m-%d"))      ## Function to estimate forecasted variance               get_est_var <- function(z) {        df = as.data.frame(z)        reg = lm(sp500 ~., data = as.data.frame(z), na.action = na.omit)        beta = t(as.matrix(reg$coefficients[-1]))        res = as.matrix(reg$residuals)        sigma = as.matrix(cov(df[,-1]))             est.var = (beta %*% tcrossprod(sigma, beta)) + (deviance(reg)/df.residual(reg))           est.var     }          risk1 = sqrt(get_est_var(data1.zoo))     risk2 = sqrt(get_est_var(data2.zoo))      You will see that the values risk1 and risk2 are pretty similar although the factors used are completely different!!! Can anyone please enlighten me as to how do we compute **model specific** estimated asset variance? Thanks for your help! ## REFERENCE Here is a reference regarding the computation of variance in the context of factor models: link. Please take a look at the formula at **the bottom of page 7**. The covariance of the factors formula is **on top of page 4**. Please note that every authoritative source on portfolio management corroborates these formulas. Here is another reference about the implementation of risk models in R: link to implementation in R. Please see **page 10 and page 12**.


title: Sampling distribution of the radius of 2D normal distribution

text: The bivariate normal distribution with mean $\mu$ and covariance matrix $\Sigma$ can be re-written in polar coordinates with radius $r$ and angle $\theta$. My question is: **What is the sampling distribution of $\hat{r}$, that is, of the distance from a point $x$ to the estimated center $\bar{x}$ given the sample covariance matrix $S$?** Background: The true distance $r$ from a point $x$ to mean $\mu$ follows a Hoyt distribution. With eigenvalues $\lambda_{1}, \lambda_{2}$ of $\Sigma$, and $\lambda_{1} > \lambda_{2}$, its shape parameter is $q=\frac{1}{\sqrt{(\lambda_{1}+\lambda_{2})/\lambda_{2})-1}}$, and its scale parameter is $\omega = \lambda_{1} + \lambda_{2}$. The cumulative distribution function is known to be the symmetric difference between two Marcum Q-functions. Simulation suggests that plugging in estimates $\bar{x}$ and $S$ for $\mu$ and $\Sigma$ into the true cdf works for large samples, but not for small samples. The following diagram shows the results from 200 times   * simulating 20 2D normal vectors for each combination of given $q$ ($x$-axis), $\omega$ (rows), and quantile (columns)   * for each sample, calculating the given quantile of the observed radius $\hat{r}$ to $\bar{x}$   * for each sample, calculating the quantile from the theoretical Hoyt (2D normal) cdf, and from the theoretical Rayleigh cdf after plugging in the sample estimates $\bar{x}$ and $S$. ![enter image description here](http://i.stack.imgur.com/Mujgq.png) As $q$ approaches 1 (the distribution becomes circular), the estimated Hoyt quantiles approach the estimated Rayleigh quantiles which are unaffected by $q$. As $\omega$ grows, the difference between the empirical quantiles and the estimated ones increases, notably in the tail of the distribution.


title: Visualizing high dimensional binary data

text: What is a good way to visualize high dimensional (say n=10) **binary** data? I remember reading something about that a few years ago. Say for instance, you want to plot / cluster pizzas based on their topping, e.g. ham, chicken, mushrooms etc.


title: How to find the exact and likelihood ratio p-value for a hypothesis test that a population sex ratio is 1:1?

text: > Skulls were excavated from a dig. The sex of each skull was determined by > anatomical appreciation. For Sites A, the data was Male=162, Female=110. > > a) Find the exact p-value for a hypothesis test about whether the sex ratio > of the population buried at this sites was 1:1. > > b) Instead use a likelihood ratio test. Find the likelihood ratio test > statistic, and resulting p-value. Compare with (a). Any help would be appreciated.


title: Fisher's exact test for results of Wilcoxon-Mann-Whitney

text: I have two groups of subjects (9 and 13 people). And data from two condition for this groups. In each condition for each subjects I have some measurements. But number of measurements are different for conditions (for second condition more then for first). First I use Wilcoxon-Mann-Whitney test for every subject to find subjects with significant difference between conditions. Second I use Fisher's Exact Test for understand is there significant difference in numbers of subjects with difference between condition for my two groups. Say, in first group 4 subjects show significant difference, and 3 in second. So my fisher's table is [[4,5],[3,10]]. Is this make any sense? and if so, Shall I use any correction for multiple comparisons?


title: How to test and examine interaction effect in 2 by 3 by 6 mixed ANOVA design?

text: ### Background: I have two groups of children - one group with normal hearing and one with permanent hearing losses - call them NH and PHL. Each child has had a hearing test, consisting of 6 different sounds (.25, .5, 1, 2, 4, and 8KHz) assessed under three different test conditions (classic, remote, and hybrid) Each child therefore has 18 data points - three conditions by six sounds. ### Questions: How can I test the hypothesis that, for each of the six sounds, the three test conditions returned the same hearing thresholds and that there is no difference between the groups in this respect? That is, both normal hearing and hearing impaired children generated the same hearing test results in the three test conditions. I am not interested in the actual differences in hearing levels between the groups as I already know how well they hear - I just need to know that their hearing levels are the same no matter how I test them. If there are differences between the groups or between the conditions, how can I examine the interaction effects to see whether the differences were greatest for particular sounds or in particular conditions?


title: How to calculate ROC(Ag) in this sound perception study?

text: I am studying individual differences in sound perception where I ask participants to listen to a pair of sounds and judge whether they are the same or different. The answer choice and scoring are as below: (number) is the score               [Correct answer = Same]      Definitely Same (1) Probably Same (0.5)      I Don't Know (0) Probably Different (0) Definitely Different (0)            [Correct answer = Different]     Definitely Same (0) Probably Same (0)      I Don't Know (0) Probably Different (0.5) Definitely Different (1)      I then convert the score to the standard $d'$, $z(\text{h})-z(\text{f})$. Therefore the $d'$ coefficient is their "score". Can someone tell me whether this is the right way to do? I have been told the best measure of sensitivity is the area under the ROC (Ag). But I don't understand how is this applicable to this design if I want to calculate sensitivity index for each person. (Detection Theory Book, p.64)               Ag= 1/2 {SUM symbol} (Fi+1 + Fi) (Hi+1 + Hi).      The index $i$ tracks the ROC points so that (F1, H1) equals (0,0), (F2, H2) is the first point to the right, and the last point is (1, 1). Can someone correct me what I have done wrong or how to calculate based on A(g) formula (i.e. what is Fi+1)


title: What is the difference between AdaBoost.M2, AdaBoost.M1, Gentle AdaBoost, RealAdaboost and the Original Adaboost?

text: I know some questions have been made for example between gentle and ada, but I was wondering about gentle and the others. For instance, I saw two different works talking about AdaBoost.M1 and AdaBoost.M2 and another work only addressing Gentle, Real and the original Adaboost. It seemed to me that the Gentle was the latest and best one among gentle, real and original ; and that M2 was superior over M1. I am not sure whoever what is the difference between (.M1 and .M2) versus the Gentle, Real and Original. There is a package in R called ada which builds on top of Real, Gentle and the Original Adaboost, but since it does not address neither .M1 nor .M2 I couldn't compare which is one is better.


title: When should one use multiple regression with dummy coding vs. ANCOVA?

text: I recently analyzed an experiment that manipulated 2 categorical variables and one continuous variable using ANCOVA. However, a reviewer suggested that multiple regression with the categorical variable coded as dummy variables is a more appropriate test for experiments with both categorical and continuous variables. When is it appropriate to use ANCOVA vs. multiple regression with dummy variables and what factors should I consider in selecting between the two tests? Thank you.


title: Understanding text Clustering

text: I have not had much luck explaining my question so please do make a comment if you don't understand what I am asking. **Core Question:** How can I approach the problem of clustering text[not text document] based on the similarity of each word. **Example :** If I have a document having a lot of text in csv such as               column1        , column2 , column3 , column4     c:\Home\Users ,data1 , data2 , data3     c:\trials\Element ,dat1 , dat2 , dat3      **What I want to do** : I want to cluster this text based on similarity of each column of the csv . How would you approach this kind of problem. which algorithms can I try. **What I have tried** : I have tried generating a proximity matrix for each of the columns using the edit distance as a distance measure . I then added the proximity matrix of each column together . I then feed into the K-means but My result was not too good.Am I approaching this wrongly? Most of the examples on clustering I have found are using numbers in their attributes I have not yet found one that uses text . I will appreciate any link to clarify this . Please I am a new to the field and some of my wording might not be direct. Please ask me any question so I can clarify the problem.


title: Analysis hierarchical circular mixture data

text: I have circular data such that multiple human participants were, each shown a color from a color wheel, asked to remember it for a "retention interval", then report it back by clicking a color wheel. Each participant was tested many times in each of two conditions testing conditions (short retention interval vs long retention interval). I have written and evaluated code that, for a given participant and condition, will model their data (represented as degrees of error between stimulus and report) as a mixture of a Von Mises distribution and a uniform distribution, estimating the proportion of Von Mises as well as the concentration of the Von Mises (the code also optionally permits estimating the location of the Von Mises, if one is unwilling to assume unbiased error, but for now I am willing to make this assumption for this data). However, the code was really hacked together from help that others provided me (see author's note in the previously linked paper) and I really don't understand the EM algorithm enough to extend it to achieve either of the following: 1) Permit estimation of data from multiple conditions from the same participant whilst enforcing a shared value of one or more parameters. This would allow me to evaluate, within each participant, the degree to which various hypotheses about the differences between conditions (no difference, difference on all parameters, difference on just the mixing proportion, difference on just the Von Mises concentration) are supported by the data. I'd likely want to use AIC-corrected (or cross-validated, if AIC is unavailable) likelihoods for comparing these hypotheses. 2) Implement (1) plus a hierarchical strategy that permits simultaneous estimation of parameters across the entire set of observed participants. I don't fully understand the concept of shrinkage yet, but I gather that it might usefully apply here, making (2) a more powerful approach than simply applying (1) to each participant individually and aggregating the results. Can anyone provide assistance with (1) or (1)&(2)? -- Following a request to clarify the model, here is my attempt. I'm more familiar with representing models algorithmically than symbolically, so I'll use R code. The data from the ith participant in the jth condition is assumed to be generated by a process like:               library(CircStats)     data_ij = c(          rvm( n_ij*rho_ij , pi , k_ij )         , runif( n_ij*(1-rho_ij) , 0 , 2*pi )     )     data_ij = data_ij[ sample(n_ij) ] #shuffle to make mixture identity latent      where `n_ij` is the number of observations made (determined by the experiment, so not a to-be-estimated parameter), `rho_ij` is the proportion of observations sampled from the Von Mises, and `k_ij` is the concentration of the Von Mises. With 2 or more conditions measured within each participant,               rho_ij = rho_i + conditionEffectOnRho_i*condition_j     k_ij = k_i + conditionEffectOnK_i*condition_j      where `conditionEffectOnRho` and `conditionEffectOnK` are vectors supplying the effect on `rho` and `k`, respectively, in each condition. I don't have strong a priori expectations for the distribution of `rho_i` and `k_i`, nor am I sure that it is appropriate to restrict the variance of condition effects across participants to zero (indeed, the covariance of the effects might be of theoretical interest), hence the `_i` associated with `conditionEffectOnRho` and `conditionEffectOnK`.


title: Compare two samples

text: similar questions have been asked but have not managed to get a conclusion from them. I am comparing two sets of samples, where ratios have been obtained for several analytes per sample. So the values are restricted to be in the [0,1] interval. One set of samples has been used as reference to set the limits for each analyte. In fact, the analytes are generally restricted to be in a range included within this interval (say between 0 and 0.27, or between 0.34 and 0.63, depending on analyte) I show an example where I have plotted the density plots for both sets of samples (reference in green) and also the QQplot. You can see that the possible ratios for this analyte are restricted to be less than 0.17. The limit was empirically set by the green distribution which is the reference. Using R commands and the sm package to graphically compare distributions               length((X[[44]][,8])      [1] 332     length((X[[45]][,8])     [1] 210     sm.density.compare(c(X[[44]][,8],X[[45]][,8]),group=c(rep(1,length(X[[44]][,8])),rep(2,length(X[[45]][,8]))),model="equal",band=FALSE)      ![Density plots. Green is reference sample](http://i.stack.imgur.com/nhBpf.png) Distributions seem normal but a normality test says otherwise:               shapiro.test(X[[44]][,8])                  Shapiro-Wilk normality test          data:  X[[44]][, 8]      W = 0.9693, p-value = 1.748e-06          shapiro.test(X[[45]][, 8])                 Shapiro-Wilk normality test          data:  X[[45]][, 8]      W = 0.9782, p-value = 0.002443      From the density plot, visually I would say that the two distributions are not different, since they are both within my limits (between 0 and 0.17) and nearly overlapping. In fact the red distribution seems to follow pretty well the reference distribution, which is what I expected. But if I test this formally, I get that there is a significant shift in the distributions.               t.test(X[[44]][,8],X[[45]][,8]) #regardless of non-normality since sample size is rather large                  Welch Two Sample t-test          data:  X[[44]][, 8] and X[[45]][,8]      t = 4.7662, df = 470.761, p-value = 2.506e-06     alternative hypothesis: true difference in means is not equal to 0      95 percent confidence interval:      0.003906504 0.009387357      sample estimates:      mean of x  mean of y      0.06604217 0.05939524          wilcox.test(X[[44]][,8],X[[45]][,8])                  Wilcoxon rank sum test with continuity correction          data:  X[[44]][, 8] and X[[45]][,8]      W = 44455.5, p-value = 6.545e-08     alternative hypothesis: true location shift is not equal to 0          ks.test(X[[44]][,8],X[[45]][,8])                  Two-sample Kolmogorov-Smirnov test          data:  X[[44]][, 8] and X[[45]][,8]      D = 0.2584, p-value = 6.951e-08     alternative hypothesis: two-sided      I also used the method proposed by Harrell-Davis in case comparing quantiles would show that in fact they are similar, but also got statistical differences (see http://www.r-bloggers.com/comparing-all-quantiles-of-two-distributions- simultaneously/).               qcomhd(X[[44]][,8],X[[45]][,8], q = seq(.1, .9, by=.1))     q  n1  n2      est.1      est.2 est.1_minus_est.2        ci.low      ci.up      p_crit p.value signif     1 0.1 332 210 0.04670665 0.04048838       0.006218272  0.0026910365 0.00928384 0.016666667   0.001    YES     2 0.2 332 210 0.05333198 0.04579936       0.007532620  0.0036450303 0.01091243 0.012500000   0.000    YES     3 0.3 332 210 0.05827195 0.05148210       0.006789846  0.0038079921 0.01072252 0.010000000   0.000    YES     4 0.4 332 210 0.06297303 0.05443425       0.008538777  0.0057801177 0.01134807 0.008333333   0.000    YES     5 0.5 332 210 0.06700634 0.05719850       0.009807839  0.0071635533 0.01236572 0.007142857   0.000    YES     6 0.6 332 210 0.07026849 0.06059111       0.009677378  0.0062359534 0.01242078 0.006250000   0.000    YES     7 0.7 332 210 0.07382629 0.06592428       0.007902005  0.0037463285 0.01219320 0.005555556   0.000    YES     8 0.8 332 210 0.07889714 0.07283984       0.006057305  0.0013907057 0.01080278 0.025000000   0.009    YES     9 0.9 332 210 0.08528227 0.08124071       0.004041560 -0.0004749558 0.00909416 0.050000000   0.079     NO      It seems that I am kind of looking for a method that would say what I want to say, but maybe I am biasing my analysis. Is there a reasonable justification (maybe statistically sound) that I can say to back my idea of both samples being comparable? * * * Thanks ever so much for your contributions. I did believe that I was going the wrong way. It seems indeed that could be an issue of not asking the right question. Effect sizes can be of great help in this context, asking "how difference are they?". Given that the distributions are not normal, I should probably not use differences in means and move to something like Cliff's delta or non-parametric approaches. Could another question to answer be "how probable is that the second set of data comes from the same population as the first -reference- sample?" Would that be an issue of goodnes of fit test? But I think I am trying to address that with the Kolmogorov-Smirnoff test, aren't I?


title: Discrete time analysis - laying out the data

text: I have two time-varying covariates (Residence Status [1 = In Residence, 0 = not in residence], and cumulative GPA. I would like to incorporate these variables into my Excel sheet before importing it into Stata and performing some survival analysis. However I am unsure about how to lay the data out. I currently have single records for each UserID and a 'time' variable so I can expand the dataset into a person-period format in Stata. However, I am unsure how to incorporate my Residence or cumulative GPA variables into my datasheet. I know that, ultimately, my data will look like the second part. However what I want to know is can I input the data like that **before** expanding the dataset? Or should I rather create a variable for each year (as illustrated in the first part of the worksheet)? ![enter image description here](http://i.stack.imgur.com/KXXNq.png)


title: Factors with only two variables in factor analysis

text: I am running a factor analysis and have a couple of questions. I have 10 variables, all of them come from a survey, with each answer is in the scale of 1 to 7. I have calculated a correlation matrix between the 10 variables, and looked for correlations of 0.5 or higher. These made me notice of 4 groups of correlated variables, two groups of 3 variables and two groups of 2 variables. Then I ran the factor analysis. According to the scree plot, I should take 4 factors. According to the eigenvalues, the 4th eigenvalue is 1.006, so it's slightly on the line, it is larger than 1, but not by much...In addition, after 4 factor the cumulative variance is 82%. After 3 it is only 72%. My question is, should I choose 4 factors or 3 ? Is it OK to have a factor with only two variables that construct it ? I know that a single variable factor is something we don't want, just like we don't want a sample-specific factors. But what is the rule regarding two variables in a factor ? I looked on the net, and found one document saying the minimum is 3. I didn't see any mathematical explanation and I am not convinced. What can you recommend me ? Thank you in advance !


title: Methods for analysing effects of percentage mortality on population data with zero/low abundances

text: I want to analyse the effects of percentage mortality from two sources (a predator and a disease) on the population abundance of a host measured at 12 sites for 8 years, with the main aim being to look for any delayed density- dependent effects of the predator or disease on the host population. I plan to use GLMMs/GAMMs. The population abundance and % mortality due to predation/disease are all measured from samples of the host (a moth larva) at each site/time period, but about 45% of samples failed to collect any hosts. Assuming these are genuine zeros, clearly there is no corresponding explanatory data (% predation/disease mortality) for these zero-abundance cases. Therefore, are there any methods to deal with this problem? Clearly, these cases could just be discarded, but would this bias the analysis etc? Indeed, prevalence estimates are also obviously very variable/unreliable when the sample size is very small, so it may also be best to discard all cases below a threshold sample size, but again does this have serious implications for biasing the analysis etc?


title: What happens to a $\chi^2_n / n$ random variable as $n \to \infty$?

text: I have read that $\chi_n^2 / n \to 1$ as $n \to \infty$. Is this true? If so how could I go about proving it?


title: Model selection and cross validation?

text: This question is a follow up for this one: Feature selection and cross- validation If I do cross-validation of a process that includes model (features) selection (so I re-do model selection for each fold), can I use the results of this cross-validation to inform the final model selection on the whole training set? For example, I may be more reluctant to include features that only were included in a few folds, or more lenient in including features that are marginal on the full model but were consistently included in most of the folds. I have the feeling that this would be wrong, because cross validation evaluated the process, not the particular features of the model. Then again, I want to have a good model, so I feel inclined to use all the information available to select the final model. Which approach is correct, if any? Would regularization be a better approach to feature selection? I thought about using bootstrapping as well, but I am using LOO cross-validation because I want/need to generate a prediction accuracy for each of my observations (which LOO does). Thanks. _Follow-up:_ I ended up ignoring the information gained from cross-validation. I chose not to include a "marginal" component based on better residuals distribution, ignoring the fact that it did get picked up quite regularly during cross-val. Not sure it was the right move, but I am happy with my model behavior. Besides, the process included PCA: the exact composition of each feature probably varied a bit for each fold, so it would have been risky or even wrong to use the information gained from cross-validation anyway.


title: Probability spread over number of days with finite wins

text: I am trying to figure out how to get odds based on a few factors and running into a few problems since I am a developer that is terrible at statistics/probability. My problem seems to be much like how the McDonalds Monopoly game functions with a few other variables. Here is what I am working with:   * 32 day contest   * 50 available cards   * 13 of 50 cards are "rare"   * Unknown number of entries   * 1,172 total prizes   * One entry each day draws 5 cards My largest problem is to be sure that of the 32 days, not all the prizes are given out in the beginning nor end but evenly distributed based off of the odds of getting a single "rare" card. Can anyone help with the math involved in figuring out what the odds should be for a single rare card based on number of prizes and days of the contest? Is there some other factor that I am missing as well? Do I need to estimate the number of people playing to figure out the probability? Help, of course, is very much appreciated.


title: R Package GBM - Bernoulli Deviance

text: All, I am _trying_ to study the GBM package in R. I. I wanted to try and figure out where the deviance, initial value, gradient and terminal node estimates came from. Please see this snippet: ![alt text](http://i.stack.imgur.com/xeCwF.jpg) To start out, I was not sure where the Bernoulli deviance came from. I expected that the deviance, would be -2 multiplied by the log-likelihood ![alt text](http://i.stack.imgur.com/GNSRh.jpg) , ignoring the weight variable w. Can anyone suggest what I am doing (or missing) wrong to dervive what is shown in GBM? II. I was also confused on where the terminal node estimates came from? Any help is very appreciated! Brian


title: Why does continuous Bayesian analysis seem to give this contradictory result?

text: Let's say you have a process that generates data according to r = sin(t) + epsilon, where epsilon ~ N(0,V) is Gaussian noise. The unconditional variance of r is 0.5 + V. Let's say we're forecasting r with a model m, and that our forecast is "perfect" in that m = sin(t). Construct v = r - m, which is the forecast error, and will be ~ N(0,V). According to Bayes, we then have p(r|v) ~ p(v|r)*p(r), which is the the product of two Gaussian PDFs, one with variance V, the other (0.5+V). This product will itself be a Gaussian with total variance T = 1/(1/V + 1/(0.5+V)). The funny thing is that T < V guaranteed, in fact T/V = (V+0.5)/(2V+0.5)! In other words, according to Bayes, the variance of your forecast error is less than the inherent noise in the data generation process itself!?! Isn't that impossible? Can anyone help me sort through this? Thank you in advance, -Jesse


title: How can I apply a Pareto tail to a truncated distribution?

text: Many income surveys (especially older ones) truncate key variables, such as household income, at some arbitrary point, to protect confidentiality. This point changes over time. This reduces inequality measures associated with the variable. I am interested in fitting a Pareto tail to the truncated distribution, replacing truncated values with imputed values to mimic the actual distribution. What's the best way to do this?


title: What is it that a statistician does?

text: When I tell my non-statistical friends I am a grad student pursuing a PhD in statistics they naturally say "oh so you want to be a professor?". I tell them no I actually plan to work in industry. Then they respond with, "and do what?". I have not found a good response to this question. I want to give them a range of interesting problems that statisticians work on, but my response is usually muddled or is too technical. I think most people picture us working in excel and calculating means and standard deviations. What is a good response to this question that is brief and interesting? Thanks! (not sure what a good existing tag is for this question)


title: What is the best Global Optimization algorithm to deal with integer values (Binary)

text: I want to use a global optimization algorithm to find the optimum set of values for n variables (n>100) and all variables take binary values (0 or 1). I am thinking in Genetic Algorithm, Particle Swarm Optimization, Differential Evolution (others?). But, which of them is suitable for integer (binary) variables? Note: I do not want to write customized function script to convert continuous values to binary x(x>=0.5)=1 and x(x<0.5)=0. This may bias the optimization process. Is there any guide on which algorithm to be used with which variable type? Regards,


title: rank data in subsets

text: I have a data asset where one "unique"ID corresponds to multiple cases, shown below. I am looking for a way to create a column ranking each ID by another column but keeping the rank unique for each ID.               ID    number     Rank(column I'm trying to add)           1       55         3     1       40         2     1       35         1     2       42         2     2       35         1     3       16         1     3       17         2      Any help would be greatly appreciated.


title: Is 'fair statistics' a thing?

text: Given that statistics can often be abused to deliberately present 'facts' to support a pre-existing viewpoint. (Lies, damned lies and statistics). And given confirmation bias. Is there an established method of choosing _which_ statistics to use, which ensures (or at least endeavours) to be as fair as possible? And does this area of enquiry have a name? **An example** If you want to compare the relative safety of flying vs. driving, you could choose from:   * fatalities per km   * fatalities per journey   * fatalities per time   * injuries per...  etc. Each one of these would yield different results which would 'pull the blanket' one way or the other. So how do you choose?


title: Intuition for uniform integrability as asymptotic light tails

text: I am trying to get an intuition to the concept of uniform integrability: $$ \lim_m \sup\lim_n \mathbb{E}[|X_n| I_{\{X_n \geq m \}}] \to 0 $$ Can this be seen as a form of "asymptotic light tails"? Can this intuition be formalized? And also- is this the same as the Lindeberg condition for the CLT?


title: Error when using 'mclust' package for mixture model estimation

text: When trying to perform a univariate _mixture model_ estimation, `R` package `mclust` produces the following _output_ with an **error** :               ----------------------------------------------------     Gaussian finite mixture model fitted by EM algorithm      ----------------------------------------------------          Mclust V (univariate, unequal variance) model with 9 components:              log.likelihood     n               df               BIC              ICL      -25576.3596860173 53940 26.0000000000001 -51436.0056895487 -84761.892259584          Clustering table:         1     2     3     4     5     6     7     8     9       3081  8923  8177  3017  8742  6568 10378  3650  1404      Error in object[as.character(G), modelNames, drop = FALSE] :        incorrect number of dimensions      **Reproducible example:**               library(mclust)          set.seed(12345) # for reproducibility          data(diamonds, package='ggplot2')  # use built-in data     myData <- log10(diamonds$price)          # determine number of components     mc <- Mclust(myData)     print(summary(mc))      **Additional diagnostic information:**               > str(mc)     List of 15      $ call          : language Mclust(data = myData)      $ data          : num [1:53940, 1] 2.51 2.51 2.51 2.52 2.53 ...      $ modelName     : chr "V"      $ n             : int 53940      $ d             : num 1      $ G             : int 9      $ BIC           : num [1:9, 1:2] -64689 -57330 -54802 -52807 -52721 ...      ..- attr(*, "dimnames")=List of 2      .. ..$ : chr [1:9] "1" "2" "3" "4" ...      .. ..$ : chr [1:2] "E" "V"      ..- attr(*, "G")= num [1:9] 1 2 3 4 5 6 7 8 9      ..- attr(*, "modelNames")= chr [1:2] "E" "V"      ..- attr(*, "oneD")= logi TRUE      $ bic           : num -51436      $ loglik        : num -25576      $ df            : num 26      $ hypvol        : num NA      $ parameters    :List of 4      ..$ Vinv    : NULL      ..$ pro     : num [1:9] 0.0606 0.1558 0.1628 0.0411 0.1655 ...      ..$ mean    : Named num [1:9] 2.7 2.86 3.03 3.24 3.39 ...      .. ..- attr(*, "names")= chr [1:9] "1" "2" "3" "4" ...      ..$ variance:List of 5      .. ..$ modelName: chr "V"      .. ..$ d        : num 1      .. ..$ G        : int 9      .. ..$ sigmasq  : num [1:9] 0.004915 0.006932 0.009595 0.000833 0.008896 ...      .. ..$ scale    : num [1:9] 0.004915 0.006932 0.009595 0.000833 0.008896 ...      $ classification: num [1:53940] 1 1 1 1 1 1 1 1 1 1 ...      $ uncertainty   : num [1:53940] 0.0132 0.0132 0.0134 0.015 0.0152 ...      $ z             : num [1:53940, 1:9] 0.987 0.987 0.987 0.985 0.985 ...      ..- attr(*, "dimnames")=List of 2      .. ..$ : NULL      .. ..$ : NULL      - attr(*, "class")= chr "Mclust"          > sessionInfo()     R version 3.1.1 (2014-07-10)     Platform: x86_64-pc-linux-gnu (64-bit)          locale:      [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C         LC_TIME=C                 [4] LC_COLLATE=C         LC_MONETARY=C        LC_MESSAGES=C             [7] LC_PAPER=C           LC_NAME=C            LC_ADDRESS=C             [10] LC_TELEPHONE=C       LC_MEASUREMENT=C     LC_IDENTIFICATION=C           attached base packages:     [1] stats     graphics  grDevices utils     datasets  methods   base               other attached packages:      [1] rebmix_2.6.1       ddst_1.03          evd_2.3-0          orthopolynom_1.0-5      [5] polynom_1.3-8      mixtools_1.0.2     segmented_0.4-0.0  boot_1.3-11             [9] mclust_4.3         MASS_7.3-34                 loaded via a namespace (and not attached):      [1] Rcpp_0.11.1      colorspace_1.2-4 digest_0.6.4     ggplot2_1.0.0         [5] grid_3.1.1       gtable_0.1.2     munsell_0.4.2    plyr_1.8.1            [9] proto_0.3-10     reshape2_1.4     scales_0.2.4     stringr_0.6.2        [13] tools_3.1.1


title: What is best formula style to write of statistical results in a scientific or medical paper?

text: I have written a dozen scientific papers and posters about medical trials and psychological trials with a lot of statistical analysis, and a common question "perturbed my certainties": > How do I write statistical significance so that it is easy to read, and > avoids repetition and panegyric? Generally speaking the most accepted formula is $$F(n,m) =f,\ p\lt 0.01$$ or $$t(n) =t,\ p\lt 0.001$$ or (Holm-Sidak) $$\alpha(n,m) =t,\ p\lt 0.001$$ or $$ F(n,m) =f,\ p=0.69$$ when there are 2 or 3 groups, I used $$F_{g1-g2}(n,m) =f,\ p_{g1-g2}=0.69,\ F_{g1-g3}(n,m) =f,\ p_{g1-g3}=0.69 \ldots\ .$$ My questions are   1. How to write formally correct and compact formula?   2. How to write formally correct repeated analysis vs group analysis?   3. How to highlight in formula the covariate variable?   4. How to write correlation?


title: Does log likelihood in GLM have guaranteed convergence to global maxima?

text: My questions are:   1. Are generalized linear models (GLMs) guaranteed to converge to a global maximum? If so, why?    2. Furthermore, what constraints are there on the link function to insure convexity? My understanding of GLMs is that they maximize a highly nonlinear likelihood function. Thus, I would imagine that there are several local maxima and the parameter set you converge to depends on the initial conditions for the optimization algorithm. However, after doing some research I have not found a single source which indicates that there are multiple local maxima. Furthermore, I am not so familiar with optimization techniques, but I know the Newton-Raphson method and IRLS algorithm are highly prone to local maxima. Please explain if possible both on an intuitive and mathematical basis! EDIT: dksahuji answered my original question, but I want to add the followup question [ **2** ] above. ("What constraints are there on the link function to insure convexity?")


title: Insignificant slope with statistically significant r value

text: I collected data of 60 countries to identify whether there is a relationship between average income per person and average life expectancy. However, I am having some trouble interpreting my lin-reg and r value.               y = 0.000437x + 67.68     r = 0.814      Does this inconsistency have to do with the nature of my data? Considering that that the data for income is in the 10-30 thousands, while the data for life expectancy is all below 100, would that influence the reliability of the persons r correlation test? If so, how? Also, I have checked and double checked both calculations a number of times, on the calculator, manually, on excel and through an online calculator and I still get the same answers. If someone could help me interpret this, that would be very helpful.


title: When is the inverse of the Fisher Information exact? (MLE)

text: I just have a quick question about MLE. Sometimes, when doing MLE problems, I see that the variance expression gotten from the inverse of the Fisher Information is exactly like what it should be and sometimes it isn't. Is there a reason that for some distributions this method is exact? Also, why is it that for some distributions, for example like $geom0(\pi)$, is it difficult to find an exact formula for $var(\tilde\pi)$? Thanks so much for your help! EDIT: $Geom0(\pi)$ is the geometric distribution where the PMF is $(1-\pi)^y\pi$ where $y=0,1,2,...\infty$ as opposed to $Geom1(\pi)$ where the PMF is $(1-\pi)^{(y-1)}\pi$ where $y=1,2,...\infty$. Sorry for the confusion!


title: Calculating the joint MGF of two random variables

text: This is an end-of-chapter question from a Korean textbook, and unfortunately it only has solutions to the even-numbered q's, so I'm seeking for some hints or tips to work out this particular joint moment generating function question. or I better state that I cannot find a way to solve this particular double integral) The question (roughly translated) states that Given the joint pdf of random variables X and Y, $$ f(x,y) = \frac{1}{\sqrt(2\pi)}e^{-x}e^{\frac{-(y-x)^2}{2}}, x \ge 0, -\infty \le y \le \infty $$ Find the joint mgf M(s,t), and for what values of s and t does the mgf exist? The following is how I approached this problem. $ M(s,t)=\int_{-\infty}^\infty\int_{0}^\infty e^{sx+ty}\frac{1}{\sqrt(2\pi)}e^{-x}e^{\frac{{-(y-x)}^2}{2}}dxdy $ $ =\int_{0}^\infty \int_{-\infty}^{\infty}e^{sx+ty}\frac{1}{\sqrt(2\pi)}e^{-x}e^{\frac{{-(y-x)}^2}{2}}dydx $ $ =\int_{0}^\infty \int_{-\infty}^{\infty} \frac{1}{\sqrt(2\pi)} e^{\frac{2(s-1)x-x^2}{2}}e^{\frac{ty-y^2}{2}}e^{xy} dydx $ $ =\int_{0}^\infty \int_{-\infty}^{\infty} \frac{1}{\sqrt(2\pi)} e^{\frac{-(x-(s-1))^2}{2}}e^{\frac{-(y-t)^2}{2}}e^{xy}e^{\frac{(s-1)^2+t^2}{2}} dydx $ I see $\frac{1}{\sqrt{2\pi}}$ in the intergral and y is to be integrated from $\infty $ to $-\infty$, so I believe that I need to make the expression in the form $\frac{1}{\sqrt(2\pi)}*e^-{\frac{(y-something)^2}{2}}*f(x)$ (pdf of a Normal distribution, and integral will be evaluated to 1 whatever the value something is) and then integrating the expression w.r.t x from 0 to $\infty$ from there. But I cannot find any useful integration technique to evaluate this double integral. I tried polar coordinate; making substitution $ x = r cos{\theta} $ and $ y = r sin{\theta} $ which will turn the integral to $ =\int_{\frac{-\pi}{2}}^{\frac{\pi}{2}} \int_0^{\infty} \frac{1}{\sqrt(2\pi)} r*exp(r*((s-1)cos{\theta}+tsin{\theta})-\frac{{r^2}(1-\frac{sin{2\theta}}{2})}{4}) drd{\theta} $ Evaluating this I got $\sqrt{\frac{2\pi}{3}}$ which I feel is incorrect as it's neither of a function of x nor y. So can anyone suggest me with any hints as to approach this question? p.s. Sorry that my English isn't good. I only spent my junior year of high school in Australia, and rest in Korea.


title: Comparing performance of two programs - problem with missing data

text: I'm measuring running times of two programs on $n$ inputs and I want to know whether one of the programs is faster. Each program has 30 seconds and 8 GB of memory for processing one input. If the program successfully solves the problem within the limits I record its running time, otherwise I record `t` when the time limit was exceeded or `m` when the memory limit was exceeded. Now I want to test hypothesis that the first program is faster than the second. But I don't know how to deal with `t` and `m` in the results. If the result were complete I would do paired t-test.


title: Large scale 1-Dimensional Gaussian Process Classification

text: I have a dataset with a small number of input points (e.g. +- 300), but millions of boolean outcomes.               Input:    Observations: #True:  #False     10          100           51    49     12          120           55    65     ...     300         1000          430   570      I would like to use a Gaussian Process to interpolate between (with all the nice things that a GP provides like error bars etc.). However the only code I have been able to find requires that the boolean outcomes are represented explicitly (e.g. +1/-1 for the GPML toolbox). This representation would be infeasibly large for a GP though. Can anyone point me towards code that can handle such a dataset (preferably Matlab code, but happy with any other suggestions too)?


title: Brisk algorithm

text: I am newbie in the field of object detection and object description algorithms. I need tutorials for Brisk Algorithm or can someone please explain in simple terms how the algorithm works? Thank you


title: What does deviance represent: a number of points?

text: I read very carefully the answer to the question "what is deviance?" already asked here, and I understand how to calculate and use deviance for model comparison, for example. However, a question remains as to _what deviance actually is_. As I was fitting a multinomial model to my data with either one or several predictors, the deviance for the most complex model was smaller by about 10000 and I was told the two models had "almost 10000 points difference in log-likelihood", which raised the question, does this number represent a number of points / a difference in number of points? I originally thought about some sort of difference in the number of well- classified outcomes between the two models. But, taking the example of another dataset I had with a total of 15707 points, 10 response categories, and depending on how many predictors I used the deviance values varied of course, but they were all around 15000, 16000 or 20000 - so way more "points" than I am analysing in my original sample. I just dont understand what these _points_ are. Any additional explanations on this would be great... thank you.


title: Image classification using histogram

text: My apologies for another stats beginner question. This feels like reinventing the wheel, unfortunately I just can't find the answer myself. Googling "average densities classification" does not help me. "classifying using histogram" brings up interesting pages, but too obscure for me. I have to classify x, and find a probability that x is of Class 0 or Class 1. I noticed that Class 0 and 1 have slightly different average densities: here are some plots showing the **average histograms for Class 0 (blue) and Class 1 (red)** : ![enter image description here](http://i.stack.imgur.com/Aq19v.png) Is there a way to use the histogram to classify new images? Intuitively, the closer the histogram is to the average histogram of the class, the higher the probability to belong to that class. I have tried using breaks of the histogram as feature, and this gives me some results, but I have a feeling I am taking the long way. Is there any ready-made R/Python functions to classify data/images based on histogram?


title: What is the purpose of precision variables?

text: Why do we need to include precision variables in a regression model (i.e., a variable that is associated with the outcome but not the predictor of interest)?


title: Inverse CDF of Mixture of Gaussians Without Sampling

text: It is fairly straightforward to compute the inverse cdf a gaussian mixture model by doing sampling (first selecting a component of the mixture according to the coefficient of the component, then sampling q ~ U(0,1) and then returning the inverse cdf of that gaussian at q). Of course, this seems pretty computationally intensive. It's not obvious to me that there is an analytical solution. One possibility that I'm considering is using some kind of regression model (either a random forest or a neural network) to approximate a function mapping the mixture coefficients, means, and variances of the gaussians to the inverse cdf. I would then train the model to reproduce the inverse cdf function. Then I would be able to get fast approximations for further queries to the inverse cdf.


title: least trimmed squared for regression

text: i'm new in statistics. hope you can help me on the following: i want to use least trimmed squared (LTS) for regression. below is the **coding in R** :               lts2_M1<-function(failure)     {     library(MASS)     y_log<-failure[,1]     x11<-failure[,4]     x2<-failure[,3]     fit0<-lqs(y_log ~ x11+x2,method = "lts", nsamp = "exact")     list(fit0)     }      the **result** is:               lqs.formula(formula = y_log ~ x11 + x2, nsamp = "exact", method = "lts")     Coefficients:     (Intercept)          x11           x2         -5.234269    -0.002685     0.110067       Scale estimates 0.2065 0.2301       **Question:** 1) Is that the correct way to perform LTS for regression? 2) In some of R example, i saw that they are using 'ltsReg'. what is the difference with the above method then? 3) From the result, residual standard error=0.2065. what is the formula used to calculate that? Would appreciate your help. thank you


title: Supervised learning based on phase space representation

text: Phase space learning Paper1 and Paper2 in neural network represents the input in higher dimension in auto associative learning. So, the network functions as an auto-associative memory where dynamical attractors are fixed points, each corresponding to one of the patterns that we want to store in the network. If the network has 3 input nodes, the the number of outputs will also be 3. Thus, for a time instant, $t$, Input is a vector of 3 features or 3 nodes. The Authors mention that a trajectory is traced out. In my case, for every time instant I have a feature vector of dimension 3 i.e there are three features. In my opinion, for every time step, t =1,2,...N, where N = number of examples, a learning algorithm is used to update the weights. Also. the Authors mention delay embedding for time series reconstruction which is an essential point in any phase space representation (Page 5 of Paper1). Based on these requirements, the following concepts are unclear and I will be thankful for an intuitive explanation. Q1. How come a feature vector is delay embedded and a trajectory is formed since we have only one example at every time instant? Q2. It is difficult to follow the training and learning as the details are not mentioned in the paper. I will appreciate an illustration so as to make the concept clear in how to apply phase space learning using any supervised algorithm where the objective function is the mean square error between the actual output of the network and the target. In my application, I only have the final target vector not a target vector for every time step.


title: Why precision and recall?

text: There are several concepts for classification performance : precision, negative predictive value, recall, and specificity. Why do people often choose only precision and recall together ? Why not negative predictive value or specificity? Why do they choose F score in addition to precision and recall? Can precision and recall alone deal with the imbalance problem between sample sizes of different classes? Thanks!


title: Statistical Significance of a learning Model

text: I built a learning model (for classification) based on a Random Forest classifier and i am asked to assess the statistical significance of its performances. Up to now, i trained and tested it on two different datasets A and B, respectively. What kind of test can i use?


title: Sample-adjusted meta-analytic deviancy macro/syntax for Excel, SPSS, MPLUS

text: I need to compute the sample-adjusted meta-analytic deviancy (SAMD) as part of an outlier search with meta-analytic data. Because we have a very large amount of meta-analytic data, computing this statistic by hand would take a little longer than forever. Does anyone know of an SPSS, Excel, or MPLUS syntax/macro available to compute SAMD? I do not have access to other statistical packages (e.g., SAS).


title: Statistical analysis and prediction of variations in anatomical features

text: I am doing a statistical analysis of variations in anatomical features, for example the shape of a pelvis bone. I found out that morphometric method can be used to quantify differences in geometry (e.g. shape and texture). My question is, what kind of method can I use to predict the shape variation of an anatomical feature based on a limited observation? Reference: 1. Predictive Properties of Statistical Shape Models 2. Prediction of biomechanical parameters of the proximal femur using statistical appearance models and support vector regression.


title: Can repeated measures be done on non-normally distributed variables using the Kruskal-Wallis test?

text: What do you recommend? I planned to do repeated-measures ANOVAs on 3 sets of data collected at 3 different points in time (initial n=1,500). The first data set appears to be non-normally distributed, and so I`ve gone on to do a series of Kruskal-Wallis H tests in SPSS (the independent variable has 4 groups, the dependent variable has 5 groups). What implications does this have on analyzing the 3 data collections over time? Should I do Kruskal-Wallis tests on all 3 data sets? Can they be used like repeated-meaures ANOVAs to compare results over time? Is there a better test? Thank you! pc


title: Thinning samples obtained using population MCMC methods

text: I am trying to use a Population Markov Chain Monte Carlo algorithm for parameter estimation in an ordinary differential equation model of gene regulation. The paper here gives a good summary of population MCMC. In most papers I have read, in the mutation step of an iteration, you only mutate one temperature. Thus if you have 10 temperatures, on average you will only propose a mutation of a chain once every 10 iterations. If you have a low crossover acceptance rate and exchange acceptance rate, like 10%, you still will likely only accept a proposal of a given temperature less than 20% of the time. This will result in high autocorrelation of the samples. Is the best way to resolve this to thin the samples by $\sim$ number of temperatures? Is it advisable to mutate every chain every iteration, and then perform a crossover or exchange move?


title: Error Bars for Peaks in Noisy Data

text: I'm doing an experiment where peaks in amplitude $A$ (the dependent variable) are expected as one varies the frequency $f$ (the independent variable). Based on our theoretical model,   * Away from the peaks, the value for the amplitude is approximately zero.   * Close to the peaks, the amplitude as a function of frequency approximates a Gaussian distribution. Currently, my method for identifying the peak is to take the mean and standard deviation of the amplitude, and defining peaks as data points where the amplitude is more than three standard deviations (above) the mean. This is possible as the data consists mostly of data where the amplitude is at zero. However, what I'm interested in is a good way to define a meaningful error bar for the value of the peak frequency at each peak. Close to the $i^{th}$ peak, what I'm currently doing (inspired by an idea of the full-width half maximum of a signal) is to take half the value of amplitude $A_{i,max}/2$ and find the range of values of $f$ for which the amplitude is higher than this value, taking the minimum and maximum of this range as the error bar. This error bar would then indicate the range of values where there is a 98% level of confidence that the peak frequency $f_i$ is indeed in that range. Is there any other meaningful means of handling the data, in particular in defining a meaningful error bar?


title: Acceptable R square Value

text: We conducted a parallel experiment with the conventional wet chemistry analysis for Sugar Cane sample for parameters POL, BRIX and Fibre content for payment purposes with Near Infrared Spectroscopy. We have developed a calibration equation for three parameters and the result for $R^2$ are as follows: 0.9 POL, 0.9 for Brix and 0.7 for Fibre. We would like to know if a $R^2$ of 0.7 for fibre can be use for prediction and if it is acceptable. Thanks


title: If a tennis match was a single large set, how many games would give the same accuracy?

text: Tennis has a peculiar three tier scoring system, and I wonder if this has any statistical benefit, from the point of view of a match as an experiment to determine the better player. For those unfamiliar, in normal rules a game is won by the first to 4 points, so long as you have a 2 point lead (i.e. if it is 4-2 you win, but 4-3 you need 1 more point, and keep going until one player is 2 ahead). A set is then a collection of games, and a set is won by the first to 6, again having to win by 2, except this time a special tie-breaker game is played instead of carrying on (except final set of Wimbledon etc...) The match is won by first to 2 or 3 sets depending on the competition. Now, tennis is also odd in that games are unfair. For any given point the server has a huge advantage, therefore each game the server alternates. In a tie-breaker game the serve alternates after every point, and it is the first to 7 points, again with a 2 point lead. Lets assume that player A has a probability of winning the point on their serve of $p_s$ and when receiving $p_r$. The question is this, suppose we A) just played tennis as a big "best of N games" match, how many games would give the same accuracy as normal best of 5 sets tennis B) just played tennis as a big tiebreaker game, how many points would give the same accuracy as normal best of 5 sets tennis? Obviously these answers will depend upon the $p_s$ and $p_r$ values themselves, so it would also be good to know C) What is the expected number of games & points played in normal tennis, assuming constant $p_s$, $p_r$ * * * **Defining "Accuracy"** If we assume that the skill of both players stays constant, then if they played for an infinite length of time, then one or other player would win almost surely, regardless of the format of play. This player is the "correct" winner. I'm pretty sure that the correct winner is the player for whom $p_r+p_s > 1$. A better format of play, is one that produces the correct winner more often, for the same number of points played, or conversely produces the correct winner with equal probability in few points played.


title: If not a Poisson, then what distribution is this?

text:
I have a data set containing the number of actions performed by individuals over the course of 7 days. The specific action shouldn't be relevant for this question. Here are some descriptive statistics for the data set: $$ \begin{array}{|c|c|} \hline \text{Range} & 0 - 772 \\ \hline \text{Mean} & 18.2 \\ \hline \text{Variance} & 2791 \\ \hline \text{Number of observations} & 696 \\ \hline \end{array} $$ Here is a histogram of the data: ![action histogram](http://i.stack.imgur.com/x0eIa.png) Judging from the source of the data, I figured it would fit a Poisson distribution. However, the mean ≠ variance, and the histogram is heavily weighted to the left. Additionally, I ran the `goodfit` test in R and got:               > gf <- goodfit(actions,type="poisson", method = "MinChisq")
     > summary(gf)
     Goodness-of-fit test for poisson distribution
     X^2                   df         P(> X^2)
     Pearson 2.937599e+248 771        0        The Maximum Likelihood method also yielded p-value = 0. Assuming the null hypothesis is: the data matches a Poisson distribution (the documentation doesn't specify this), then the `goodfit` test says we should reject the null hypothesis, therefore the data does not match a Poisson distribution. Is that analysis correct? If so, what distribution do you think will fit this data? My ultimate goal is to compare the mean number of actions between 2 samples to see if the means are different; is checking the distribution even necessary? My understanding is the typical tests (z-,t-,$\chi^2$ tests) don't work for Poisson distributions. What test should I use if the data is indeed Poisson- distributed?


title: ARIMA-GARCH model for exchange rates

text: I am currently working on a model for exchange rates, and I want to use an ARIMA-GARCH specification. More precisely, I work on the log-returns series. First of all, I perform multiple KPSS and ADF tests to get the conclusion that the series is stationnary up to a trend. I then use the auto.arima function to select the best model based on AIC or BIC series. My problem arises when I want to check the model adequacy by looking at the residuals. Since I work with data from January 1999 to December 2012, I have a time series of 3428 observations. When I look at the p-values of the Ljung-Box test, I get the conclusion that the different lags are not significant up to lag 23, and become significant from this lag onwards. Which conclusion should I make ? If I take the usual $\log(n)$ rule, I have to look around lag 8 and I thus have no problem. But sometimes in the literature, it is asked to look around $\sqrt(n)$ which is around lag 59, and I can not accept the model. Moreover, I want next to simulate trajectories for one year, that is, 250 points, and we can find advices saying it would be better to look around the lag corresponding to the projection horizon. Thus here, I would have a big problem, since I won't find a model which Ljung-Box p-values significative up to lag 250. Would you have any advice for me ? Thanks in advance. Ludovic


title: How to check which model is better in state space time series analysis?

text: I am doing time series data analysis by state space methods. With my data the stochastic local level model totally outperformed the deterministic one. But the deterministic level and slope model gives better results than with stochastic level and stochastic/deterministic slope. Is this something usual? All methods in R require initial values, and I read somewhere that fitting an ARIMA model first and taking values from there as initial values for state space analysis is one way; possible? or any other proposition? I should confess here that I'm totally new to state space analysis.


title: Understanding the Pareto distribution as applied to wealth

text: The Pareto distribution can be used to give a pdf for the wealth of a person chosen randomly from a population. (In fact, this was its origin. See, for instance, http://en.wikipedia.org/wiki/Pareto_principle ). I would like to explore the reciprocal question: Given the total amount of wealth in a population, what is the pdf of the portion that a randomly chosen person has. I conjecture that this is simply a constant times the Pareto distribution. More interestingly: What is the shape of the distribution curve, if the richest person would be at the 0 point on the x axis, the next richest person to the right, and so on - we would see a monotonically decreasing curve. But what is its shape? What is its derivative? It's quite likely that I'm not phrasing that question properly. Let me ask a more basic question: What is the appropriate terminology to explore the question? Give a probability distribution applied many times over, what is the shape of the resultant allocation curve?


title: Output layer of artificial neural networks when learning non-linear functions with limited value range

text: I want to approximate a non-linear function with a limited value range by an artificial neural network (feed forward, back propagation). Most tools and literature availabe suggest linear functions for the output neurons when doing regressions. However, I know a priori that my goal function is of limited range, therefore is it reasonable to use a function for the output neurons with limited value range, too? To be more explicit: My target function's values are in the range between 0 and 1, but the neural net does predict occasionally values that exceed this range (e.g. -1.3). Can i prevent the net from doing so, and is it reasonable? Thank you very much, Julian


title: What to do with negative Chow test?

text: I have just done a Chow test on a regression in order to see whether there is a structural break. I am a bit stumped however as my Chow test returns a negative number. Now what do I do? More specifically, this expression (from the Wikipedia entry for the Chow test): $\frac{(S_c-(S_1+S_2))/(k)}{(S_1+S_2)/(N_1+N_2-2k)}$ turns out negative for me. EDIT: It turns out it was a simple programming error.


title: Calculating average variance extracted (AVE) in R for checking discriminant validity (Fornell-Larcker criterion)

text: We are using `lavaan` in `R` to calculate CFAs (confirmatory factor analyses) and SEMs (structural equation models). We now want to test whether two latent constructs can be assumed to be unrelated (discriminant validity). Refering to > Fornell, C. & Larcker, D. F. (1981). Evaluating structural equation models > with unobservable variables and measurement error. Journal of Marketing > Research, 18, 39–50. we need to calculate the so called average variance extracted (AVE). I found the following R-help-mailing list post. However, I don't think the output is the AVE I need for the Fornell-Larcker criterion. Values are calculated for each item and not factors. Does anybody know how to obtain this value within `lavaan` or any other package? If it is not possible to automatically obtain the AVE for each construct, can somebody point out how to calculate it based on the `lavaan` output? Help would be much appreciated! Thanks.


title: Performance Evaluation for multiple classifiers in machine learning with Wilcoxon Rank Sum Test

text: I use four feature descriptors for five independent runs with each run randomly selecting the training set and testing set from the total labeled samples. SVM is the classifier. After acquiring the classification accuracy of each run for each descriptor, I get a accuracy matrix whose size is $5 \times 4$. Then I use the `ranksum` function in `MATLAB` to do Wilcoxon run sum test from pair-wise. Then I get the `p-value` matrix shown as follows: ![enter image description here](http://i.stack.imgur.com/wDyT3.jpg) I select $\alpha = 0.05$. If the `p-value` is less than $0.05$, then the difference is significant, and if the `p-value` is larger than $0.05$, then the difference between the two features is not significant. Is it right? Is the total procedure of statistical tests right? Thank you.


title: How to estimate distribution from DataBase?

text: I am currently working on a programme to compute the mutual information between a variable and a set of variables. Since I am working with discrete data, I am computing $\sum_{a\in A,b\in B,c\in C}P(A=a,B=b,C=c)*\log_{2}({P(C=c)P(A=a,B=b,C=c)\over P(A=a,C=c)*P(B=b,C=c)})$ where I look up in my DataBase, how often there is an event "A=a,B=b,C=c" and so on. From these counts I want to approximate the probability of these events. My supervisor asked me, if I really believe my approximation to be reasonable, since my DataBase is approximately 40000 entries long, each variable has 200 discrete values, and I use conditional independence tests with up to 6 variables in the conditional set. My argument is, that I estimate the distribution Maximum Likelihood-style, but on the other hand, I have to admit that I have too few data to back up my Maximum Likelihood estimation. So: what would you do? Stick to the Maximum Likelihood estimation? In WEKA, for example, they also do this. However: can you give me any hints what I should do?


title: How do I introduce features and their confidence values into classifiers

text: I would appreciate some advice on how best to weight or give more importance to a percentage with a larger denominator. Eg   * A. 1 out of 2 = 50%   * B. 5 out of 10 = 50%   * C. 500 out of 1000 = 50% Some of the data is sparse and continually emerging so   * A) could just be a blip and could change,    * B) is emerging and may vary slightly    * C) is becoming fairly certain Also, I'm trying to minimize the number of parameters so is there a way to weight the relative 50%'s for each scenario instead of creating a separate parameter?


title: How can I partition a distribution into two sub-populations with fixed bias? (simulation)

text: I am trying to simulate a selection model for a variable $Y$ dependent on covariate vector $X$, so that two groups/sub-sets $S=(0,1)$ of observations on $Y$ are created, which have a fixed difference in expectations (means) $\theta$ (i.e., selection bias) $$E(Y|S=1)-E(Y|S=0)= \theta$$ I am looking for solutions to approach this problem, where in principle I could also simulate two variables $Y_1$ and $Y_2$ and overlay these as a mixture distribution $Y$. However, note that the selection mechanism $S$ should be conditionally independent of $Y$ given $X$ (ignorable selection) $$P(S|Y,X) = P(S|X),$$ but not be marginally independent $$P(S|Y) \ne P(S).$$ Whereas it is easy to parameterize a model for $P(S|X)$, and also relate $E(Y|X)$, it seems very difficult to me to specify this model to achieve a fixed $\theta$ and keep variances of $Y$ across sub-sets approximately constant (the latter is not absolutely necessary, but desirable). If the mixture approach is chosen, I do not know how to overlay distributions in dependence of $X$. Thanks for suggestions!


title: Subjects with different order of magnitude values in panel data

text: I have quarterly sales data for a variety of stores and would like to estimate the effect of a regulation on sales. A panel type model would appear to be appropriate in this case, with the regulation as a dummy variable. However, the quantity of sales varies by two-orders of magnitude between stores. How should I account for this? Apologies for any problems with the question, I am something of a beginner.


title: k-way MinMax spectral clustering

text: Is there a k-way MinMax Cut Spectral Clustering which can be easily implemented? In the spectral clustering tutorial I found only 2-way MinMax cut.


title: What are the key statistical concepts that relate to data mining?

text: What topics in statistics are most useful/relevant to data mining?


What are pivot tables, and how can they be helpful in analyzing data?


title: divergence between nls (simple power equation) on non-transformed data and lm on log transformed data in R

text: I'm testing different ways to fit a simple power model to a dataset of 5 points. I don't get same results if I use nls on non-transformed data or lm on log10 transformed data in R whereas Excel gives exactly the same results (it's power fitting might be a hidden linear fitting on log transformed data). Here is an example of the code I'm using:               # data     r <- c(0.0778,0.0429,0.0320,0.0216,0.0148)     t <- c(4,10,40,100,400)     df <- data.frame(r,t)     df$log.r <- log10(df$r)     df$log.t <- log10(df$t)          # models     m1 <- nls(r ~ m * t^(-n), data=df, start=list(m=1, n=0.5))     m2 <- lm(log.r ~ log.t, data=df)          # simulate data     t.range <- seq(1, 1000, length.out=30000)     r.pred1 <- predict(m1, newdata=data.frame(t=t.range))     p1 <- data.frame(t=t.range, r=r.pred1)          log.t.range <- log10(t.range)     r.pred2 <- predict(m2, newdata=data.frame(log.t=log.t.range))     p2 <- data.frame(log.t=log.t.range, log.r=r.pred2)          # plot results     p1$log.r <- log10(p1$r)     p1$log.t <- log10(p1$t)     library(ggplot2)     ggplot(df, aes(x=log.t, y=log.r))+       geom_line(data=p1, color="red")+       geom_line(data=p2, color="blue", linetype="dashed")+       geom_point()+       ggtitle("red: nls on non-transformed data\ndashed blue: lm on log10 transformed data")+       theme_bw()      Thus, some questions:   * Which method (nls vs lm) would be the "best" in this case using R?   * Am I missing something in the use of nls?   * Could this be due to the small number of points (only 5) I'm trying to fit?


title: Big Data? Have we solved Small Data yet?

text: There has been a lot of attention on Big Data recently, where the problems are often more logistical (how to deal with large volumes of data) rather than statistical. At the other end of the spectrum we have "small data", where there is no logistical problem, but instead we have the problem of how to estimate the parameters of the model from a limited number of observations. This is often rather tricky as the true nature of the physical system may be non-linear and hence require many parameters to capture the "true" model. My question is, have we solved the problem of small data sufficiently well that we should all focus our energy on developing methods for big data?