title: Can LSA be used for document similarity?

text: I have to find the similarity between two documents. The two documents are simple text documents and i have to report a score. I was using cosine similarity initially. But I was told that LSA is a better means. But when I got to read a few tutorials i always noticed that they used more than two documents. So is it effective to include LSA when i need to find similarity between two documents alone?


title: Multivariate models for ordinal variables

text: I would like to fit a multivariate regression model of an ordinal random vector using ordinal variables as covariates. I am wondering if there is anything implemented in R or other software that could do that. Or, alternatively using a Bayesian approach.


title: Do Random Forests exhibit prediction bias?

text: I think this is a straightforward question, although the reasoning behind why or why not may not be. The reason I ask is that I have recently written my own implementation of a RF and although it performs well it is not performing quite as well as I had expected (based on the Kaggle Photo Quality Prediction competition data set, the winning scores and some of the subsequent info that came available about what techniques were used). The first thing I do in such circumstances is plot prediction error for my model, so for each given prediction value I determine the mean bias (or deviation) away from the correct target value. For my RF I got this plot: ![Prediction value versus bias from correct target value](http://i.stack.imgur.com/k95AK.png) I'm wondering if this is a commonly observed bias pattern for RF (if not then it could perhaps be something specific to the data set and/or my implementation). I can of course use this plot to improve predictions by using it to compensate for the bias, but I'm wondering if there's a more fundamental error or shortcoming in the RF model itself that needs addressing. Thankyou. == ADDENDUM == My initial investigation is at this blog entry Random Forest Bias - Update


title: Selection of "best split" at a given node in RF

text: At a given node, random forests sample `k` (aka `mtry`) variables (predictors). Then they pick the variable for which they get the best split. How is this specifically determined? From what I understand RF do something like this:               for variable in variables_sampled:      for threshold in possible_thresholds:          check data purity (i.e. entropy) according to the divisions A and B where           A: (variable > threshold) and B: (variable <= threshold)      Is this understanding correct? How are `possible_thresholds` chosen?


title: Finding MVUE for a Poisson distribution

text: The number of breakdowns Y per day for a certain machine is a Poisson random variable with mean $\lambda$. The daily cost of repairing these break downs is given by $C=3Y^2$ If $Y_1, Y_2, ..., Y_n$ denote the observed number of breakdowns for $n$ independently selected days find an MVUE for $E(C)$. We can use the Rao-Blackwell Theorem. We know that $E(C) = E(3Y^2)=3[V(Y) + (E(Y))^2]$ and $E(Y)=\lambda=V(Y)$. With some calculations we see that $E(Y^2)= \lambda + \lambda^2$ $\sum_{i=1}^n Y_i=\bar Y$ is a sufficient statistics for $ \lambda$ So I am assuming we can replace $\lambda$ with $\bar {Y}$ I am just unsure how to go from here. Can you guys help me pull strings together?


title: Drawing from the multivariate Student's t-distribution

text: So I am trying to figure out if there is a nice decomposition for sampling from the multivariate Student's t-distribution like there is for sampling from the multivariate normal distribution: http://en.wikipedia.org/wiki/Multivariate_normal_distribution#Drawing_values_from_the_distribution There are quite a few `R` packages out there for doing this but none of them seem intuitive to me to implement. What I am trying to accomplish it to sample from the multivariate t-distribution given that I have a specified number of degrees of freedom, a known mean and covariance matrix. I need to be able to sample from the multivariate t-distribution because I have analytically integrated out parameters from my posterior distribution. I.e., $$p(\psi|y)=\int_{R^p}\int_{R^+}p(\beta,\sigma,\psi|y)d\sigma^2d\beta$$ which after the integration results in $\psi\sim\text{Multivariate Student-t}(\mu,\Sigma,\nu)$. Since I have a closed analytical form for $\psi$, I would like to be able to sample directly from the multivariate Student t-distribution.


title: PCA and Discriminant Analysis Under Multicollinearity

text: I'm trying to run discriminant function analysis on 27 acoustic, temporal and descriptive variables I obtained from animal calls. However, a number of these variables are collinear because they were calculated from other variables. For example: f0 range = f0 high - f0 low. Previous work using the same variables utilized 2 PCAs to create component scores (PCA 1: Frequency; PCA 2: Component Scores 1, Temporal and Descriptive, on which the discriminant analyses were run (individuality, context etc.) However, I think this was not the correct option because 1) the PCA correlation matrix is "not definite positive" and 2) the components aren't going to be correct anyway because the PCAs aren't detecting these collinear relationships. Accordingly, I was wondering if there is a method that allows me to reclassify the variables into scores before DA that takes collinearity into account? Or, is my best option to remove the derived variables or their constituent elements before running the analysis (Using the example above, to remove f0 range OR remove f0 high AND f0 low)?


title: Linear regression multiple output in python

text: Say I have predictor array `x:(n,px)` and a predicted array `y:(n, py)`. What would be the best way in Python to calculate all regression coefficients (linear) from `x` to each dimension of `y` (`1...py`)? The output of the whole thing would be a matrix `(py, px)` (for each output, `px` parameters). I could easy iterate over output dimensions, however that would be inefficient as I will recalculate the pseudo inverse matrix of `x`. Is there any efficient implementation out there?


title: LGOCV caret package R

text: i am learning data mining through book . During classification chapters about Neural Networks the authors have below code. I have below questions:               ## pre2008 <- 1:nrow(training) ## training is a dataset that has training data     ctrl <- trainControl(method = "LGOCV",                          summaryFunction = twoClassSummary,                          classProbs = TRUE,                          index = list(TrainSet = pre2008),                          savePredictions = TRUE)     nnetGrid <- expand.grid(.size = 1:10, .decay = c(0, .1, 1, 2))     maxSize <- max(nnetGrid$.size)         set.seed(476)         nnetFit <- train(x = training[,reducedSet],                           y = training$Class,                      method = "nnet",                      metric = "ROC",                      preProc = c("center", "scale"),                      tuneGrid = nnetGrid,                      trace = FALSE,                      maxit = 2000,                      MaxNWts = 1*(maxSize * (length(reducedSet) + 1) + maxSize + 1),                      trControl = ctrl)      LGOCV - when do we use it? I read the post, but still not clear. the post says that it is a variant of LOOCV for hierarchical data. but my Y variable is not hierarchical :( twoClassSummary - can it be used only when we have two classes? can i used it for say Iris data? * * * LGOCV is also known as Monte-Carlo Cross Validation. More details are available here.


title: Comparing two time period of same binary variables

text: We are a financial institution which falls under this ‘Big Bank’ category. As a part of our routine analysis on input variables that feed several of our risk models, we would like to setup an input monitoring system that can trigger any drastic changes in the trend/pattern in our input variables over different period of time spanning form 2007 to 2013. We have decent understanding on how to perform this analysis on continuous variables, but we don’t know how to do this on binary variables that has two outcome(yes=1, No=0).I would like to take the same variable at different time period ( for e.g Var1 from Jan2007 to Jan2008 Vs Var1 from Jan2008 to Jan2009) and see if there is a significant change in the trend/distribution. I would assume my question boil down to a point where we are comparing the binomial distribution of the same variable at two different time period. If I do that I am afraid my samples wont be independent to each other, that threw away the option of using Chi-Squared test, What’s the right strategy to perform this analysis, how do I carry out this test preferably in SAS, if SAS doesn’t have right procedure or function, then R will be the next choice.


title: Gibbs sampling in document clustering

text: I'm trying to understand the idea behind the Topic Models in document clustering. In Latent Dirichlet Allocation, it is necessary to approximate the posterior distribution of topics over the document. One of the methods is Gibbs sampling - I'm not sure if I understand how it works - I'm basing my understanding on this paper: http://faculty.washington.edu/jwilker/559/SteyversGriffiths.pdf. It starts with random word-topic assignment and iteratively, new samples are computed, conditioned on previous results, and after some number of iterations, it starts to converge. But I don't really get the idea, why multiple samples need to be kept. Isn't the result of the last iteration the best approximation of the desired distribution? I will be grateful for possibly very simple explanation - my knowledge in this area is really poor.


title: Understanding of GLM

text: I have scoured around, reading posts on Cross Validated (Difference between logit and probit models) and also looking at references including Dobson and McCullagh and Nelder, e.g. http://www.statsci.org/glm/books.html so I am aware that this topic is well trodden. Nevertheless, I am trying to articulate and formalize my understanding of GLM and while several posts have helped me with that, I am conscious of gaps and the possibility of an unsound foundation of my understanding. In simple linear regression we have some set of observations $(x_i, y_i)$ pairs and treat $y_i$ as a realization of a random variable, $Y_i$, distributed as $Y_i \sim N(\mu_i, \sigma^2)$. The means ($\mu_i$) depend on the predictor but the variance is constant. We model $\mu_i = \beta_0 + (\beta_1 \times x_i)$ (or matrix equivalent ) which I believe is the same as saying $Y_i = \mu_i + \epsilon_i$ with $\epsilon_i \sim N(0,\sigma^2)$. I am not sure but I think a correct way to state $\mu$ is $\mu_i = E[Y_i|X=x]$. Can someone confirm that? Anyway, we might transform the response to achieve linearity (maybe taking logs) and in that case we are modelling $log(Y_i) = \alpha_0 + (\alpha_1 \times x_i) + \epsilon_i$ where $\epsilon_i$ now is assumed to have a log normal distribution. We generalize by decomposing the model into:   1. A structural component ($\beta_0 + (\beta_1 \times x_i)$)   2. A link g(.)   3. A response distribution (or random component) (member of the exponential family - Guassian, bin,gamma etc) Let's say that the observations are assumed to come from a distribution in the exponential family and to keep things simple assume that it is the Gaussian distribution. Again, the expected value comes out as $E[Y_i|X=x] = \mu_i$, but does so from the first derivative of the b(theta) term in the exponential family form (http://www.amazon.com/Generalized-Edition-Monographs-Statistics- Probability/dp/0412317605) (pg 29). In the case where we assume a binomial distribution, the expected value comes out at np but it still comes from the first derivative of the b(theta) term when the distribution is expressed in the exponential family form. Note, I don't fully appreciate or understand the derivation of the mean or variance at this stage, perhaps someone could provide a layman's explanation? Instead of modelling the mean as was done for simple linear regression, we now model a transformation of the mean so instead of saying $\mu_i = \beta_0 + (\beta_1 \times x_i)$ we are saying $g(\mu_i) = \eta_i = \beta_0 + (\beta_1 \times x_i)$ where g is some link function (invertible and differentiable). I think this is a key (yet still somewhat confusing to me) distinction between SLR and GLM. In SLR we transform the response ($y_i$) and model that, in GLM we transform the expected value ($\mu_i$ in the case of the Gaussian example) and model that. Another way of saying it is that in the SLR case we model $E[g(Y_i|X=x)] = \beta_0 + (\beta_1 \times x_i)$ but in the GLM world we are modelling $\eta_i = g[E(Y_i|X=x)] = \beta_0 + (\beta_1 \times x_i)$. My question is around verifying my understanding and statement of the essence of the foundations of GLM and the differences between GLM and the traditional linear model. Thanks.


title: How does this trick for estimating bit error rate break down?

text: Imagine you are receiving a message over and over via a lossy data path. The path causes bit errors but does not affect the length of the message (or shift any bits). You don't know the _actual_ message but you will get many copies of it so you can reconstruct it (e.g. with a median filter). You want to estimate the error rate of the channel by observing the sequence of messages. I solved this problem a while ago by "learning" the expected message and counting errors as $\text{population}(\text{learned} \oplus \text{current})$. Today it struck me that the learning component may be unnecessary: For low error rates you could approximate the rate as $\frac{\text{population}(\text{previous} \oplus \text{current})}{2}$. That is effectively measuring the previous frame against the current _and vice versa_ simultaneously. At low error rates half of the errors can be expected to be due to $\text{previous}$ and half due to $\text{current}$. My assumption about the independence of the errors in $\text{previous}$ and $\text{current}$ breaks down as the error rate increases and repeated errors become more likely. How does the actual error rate impact the rate given by my estimation?


title: How to generalize Particle Filters (w.r.t. multiple states)

text: I'm using particle filters for inference in a hidden markov model with an infinite state-space. My current state-variable is multidimensional and there are interdependencies between some dimensions. I thought therefore, that it would be beneficial, if I split it up into sub-states and model ther dependences in a small bayesian network. Thus, I would only maintain particles for each of the smaller-dimensional sub- states and hopefully get better results with fewer particles. I thought, that representing N one-dimensional distributions needs fewer samples than representing the N-dimensional joint-distribution, represented by the Bayesian Network. Does this make sense? This leads to a Dynamic Bayesian Network, but I don't know how sampling, likelihood-weighting and resampling should be adapted for such a particle filter with multiple states. **Update:** Looks like `Nonparametric Belief Propagation` is an algorithm (class) that extends particle filtering techniques to arbitrarily structured graphical models.


title: Sum of two independent Student t variables with same dof is t distributed?

text: What is the distribution of the difference of two-t-distributions suggests that the sum of two t-distributions is never t distributed. With t distribution I mean the (non-standardized) t distribution with location and scale parameter. Now, let $Y_1,Y_2$ be independently t distributed with same dof $\nu$, location $\mu$ and scale $\sigma$. Then the sum is given by $X = Y_1+Y_2 = (\mu+\sqrt{\nu/V}\sigma Z_1) + (\mu+\sqrt{\nu/V}\sigma Z_2)$, where $V$ is $\chi^2$ distributed with $\nu$ degrees of freedom and $Z_1,Z_2\stackrel{iid}{\sim}N(0,1)$. Since the sum of two independent $N(0,1)$ is $N(0,2)$, it follows that $X = 2\mu+\sqrt{\nu/V}\sqrt{2}{\sigma}Z$, where $Z\sim N(0,1)$. Thus $X$ is t distributed with unchanged dof, twice the location and new scale $\sqrt{2}\sigma$. So am I missing something?


title: Evaluation of k-means output for >3D

text: I'm implementing the k-means algorithm (in R Map-Reduce) and I wanted to verify if the output I'm getting is close enough to the true centroids of the cluster. This is how I'm verifying with a 2D dataset currently: I plot both the dataset and the centroids I've got as output and see if the centroids are close to the centre of the clusters visually. I think I can do that for 3D data too. But I don't know how to verify this with higher dimensional data that cannot be plotted. It looks really stupid to plot the data and visually verify each time, right? So let me tell you why I'm doing this: > The centroids don't come in a particular order. The 1st centroid in this > trial might be in the 2nd position in the next, so I can't find the distance > between the matrix of my output and the matrix of, say, the R's default > kmeans output (if I'm verifying my output with R's **kmeans** ). Sorting > with respect to any one dimension to compare sounds stupid, since any > dimension can be a lot more sensitive to the data when compared to another. So, for now, I'm verifying 2D data visually. Do I have to use dimensionality reduction? Does someone have ideas about how I can verify higher dimension data?


title: Is there a formula or rule for determining the correct sampSize for a randomForest?

text: I'm playing with a randomForest and have found that generally increasing the sampSize leads to better performance. Is there a rule / formula / etc that suggests what the optimal sampSize should be or is it a trial and error thing? I guess another way of phrasing it; what are my risks of too small of a sampSize or too large (overfitting?)?


title: Neural network and logistic regression not yielding output for certain inputs

text: I am working on a project on gesture recognition(in c# language). In some part I had to use neural networks using the sigmoid activation function. I am using the aforge.net(http://www.aforgenet.com/) library for the same. The system developed gets trained when I use a certain part of the training data. However when I use a different portion of the training set I give ans as NaN(not a number). Initially I thought it would be due to some variable overflow. However the system supports large number of inputs from one section of the training data and doesn't support as much as one inputs from the other section. Now this makes me reach a conclusion that something is wrong with that input. However I am not able to what may be probably wrong. I have tried fitting a logistic regression function which yielded similar results. Are there any constraints while modelling the input. Any help would be appreciated.


title: Time variable in Longitudinal data set mixed model question

text: It has been a few years since I fit a mixed model, so I have gone on a massive review session on old notes, books (Pinheiro and Bates, Faraway, etc) and going through the posts on SO and CV about mixed models. It has been great, but I am left with a few questions that are likely so basic they are missed in most of these materials. The data set I am working with has a time column that starts at 1 and goes to 38 by period. I also have 3 variables, one for each, year, quarter and period. In my fixed effect, should I be modeling sales~Time+quarter+year+period or should I ignore the quarter, year and period variables and just use my aggregated time variable. Should I put ordered(time) inthere? On the one hand, I imagine having the additional quarter/period/year variables in there might be redundant, since they are already in the Time variable, but at the same time, it might be good to try and use them to control for seasonality? I am wondering if my fixed effect should be               Sales~Time+Quarter+year+period+Region+Policy      or               Sales~Time+Region+Policy      Also, if anyone has advice on fitting models with more than 1 or two random effects, it would be much appreciated. In all of my notes, its mostly just specific designs with 1 random effect (sometimes intercept and sometimes slope + intercept), but rarely multiple. I currently am trying to look at 5 variables I would like to be random effects (at least look at them). Is it as simple as (1|A) + (1|B) + (1|C) + etc in lmer? Any trouble I am walking myself into? Thank you as always CV. Edit: I am interested in a change in policy that has potentially harmed sales. I have a 3 year data set and in the middle of that data set we changed our policy and I am fitting LMM to see first if there is actually a negative effect and second to get a vague idea of how big it is, if its a decrease of .001%, we don't care, but if its a larger number it should be examined further (which naturally will be my next task). Edit2: For random effects (which I am still in the process of tinkering with) I have the product, who is paying and purchase quantity. What makes the cut here is still a work in progress, currently using LRT and AIC (mostly AIC, because there seems to be reservations about LRT with lmer sometimes). There are 38 months in the set and each month has sales values for each type of payment recieved (there are 5), each product sold (3), the region it is sold in (~30 here). The other variables I have are sales, number of things sold, the team that did the selling (in just because the change in policy relates to it) and a dummy variable for change in policy. The model I am working around right now is looking something like:               m1 = lmer(Sales~ Time+Policy+Team+(Product|Territory)+(salesqty|Territory)+ (payer|Territory), data=data )       Though like I said, this is by no means final, just what I have running in R at this second.


title: How to do community detection in a weighted social network/graph?

text: I'm wondering if someone could suggest what are good starting points when it comes to performing community detection/graph partitioning/clustering on a graph that has **weighted** , **undirected** edges. The graph in question has approximately 3 million edges and each edge expresses the degree of similarity between the two vertices it connects. In particular, in this dataset edges are individuals and vertices are a measure of the similarity of their observed behavior. In the past I followed a suggestion I got here on stats.stackexchange.com and used igraph's implementation of Newman's modularity clustering and was satisfied with the results, but that was on a unweighted dataset. Are there any specific algorithms I should be looking at? Many thanks for any advice ~l


title: How to detect intermittent time series?

text: I need to automatically identify if a time series is intermittent or not. Depending on the result I'll use one or another method for forecasting it. Is there any test to detect intermittent time series?


title: glmer() warning message help

text: I'm running a glmm with the lme4 package in R that looks like the following:               glmer(y ~ factor1 + factor2 + (1 | RE), data=data, poisson)      And I get the following warning message:               Warning message:     In mer_finalize(ans) : gr cannot be computed at initial par (65)      Both `glmer(y ~ factor1 + (1 | RE), data=data, poisson)` and `glmer(y ~ factor2 + (1 | RE), data=data, poisson)` work just fine without a warning message. Can someone explain to me what this warning means, why I'm getting it, and what, if anything, I might do about it?


title: Comparing MLE parameters using bootstrapped confidence intervals

text: I have MLE curve parameter estimates for 3 populations (2 per population), and am looking for a clever way to compare them. At the moment, I am non- parametrically bootstrapping my datasets in order to create _n_ datasets, which I then compare using a standard MANOVA procedure. I have a number of issues with this, one of them being that the variances of my bootstrapped parameters are heteroskedastic and non normally distributed. Ideally, I was thinking that I could compare the bootstrapped confidence intervals around my parameters, i.e. show that the bootstrapped confidence intervals around each of my popuation parameters "overlap". I was looking at bivariate comparisons, but many of these rely on the fact that it is a MV gaussian distribution, which according to multivariate shapiro tests, I do not have. To wrap things up, my questions come down to this:   1. What would people see as the most sound procedure?   2. I had a look through this post Comparing points in a bivariate space, however, am curious about the second answer and if anybody has any suggestions as to where I could find information on non normal mulitvariate analysis. Thanks for your help.


title: It it legitimate to apply a one-way ANOVA to data that have been normalised to the untreated controls?

text: I was wondering if it is OK to use one-way ANOVA after normalisation to the untreated controls? This is in an animal model of wound healing where there are 4 wounds per animal, one untreated and 3 treated with different drugs. I have noticed that the healing rate differs between animals so in order to reduce the standard deviation and account for the inherent healing rate in the animals I normalise the data to the untreated controls in each animal so that the controls = 1. This means that they don't have a SD. Is it still OK to use one-way ANOVA or would there be a better test to use?


title: Unbiased estimate of the semi-partial correlation

text: * Is the sample semi-partial correlation a biased estimate of the population semi-partial correlation?   * If it is biased, what is an unbiased estimator of the population semi-partial correlation?   * Are there are any references in the statistical literature that discuss this?


title: How to derive jackknife bias for variance and mean

text: I am having a hard time to understand how one derives the jackknife bias for the variance and mean. 1) Why do we need an inflation factor of $(n-1)$ when calculating the jackknife bias of the mean? 2) How does one derive the jackknife bias for the variance?


title: GEE (or GLMM) in SPSS: Interpreting outputs and model selection

text: I am attempting to analyze my (experimental psych) data in SPSS, and I have a few questions regarding the kind of analysis I should be using (GEE or GLMM), how I should be interpreting the output, and how I should be selecting the best fitting model. (Disclaimer: Regression analyses in general are fairly new to me, as my data can usually be analyzed using simple non-parametric statistics or ANOVAs). Here is a little about my current data set: the DV is binary (a yes/no response) and I am interested in looking at the effect of several IVs (all categorical, 2 levels each) on that measure. Given that all my variables are experimentally manipulated, I believe they would be considered fixed rather than random effects. One of my IVs (condition: A, B) is a within-subjects factor, and the rest of my IVs are between-subject factors. These are my main concerns: 1) Should I be using GLMM or GEE? I get the same results for both, in terms of which effects are significant (which is good), but the parameter estimates are obviously different. I’m leaning towards GEE, given that I don’t have any random effects. Is this appropriate given the nature of my data? 2) If I only have dichotomous variables, does the correlation structure in GEE really matter? 2) I’m not quite sure how to interpret the parameter estimates for either GEE or GLMM, particularly in the case of interaction effects. Based on what I’ve read, I have been interpreting my main effects using the exponential coefficient (this is the odds ratio, yes?). So if I set my reference category to 0 (a “no” response), then I should be interpreting the exponential coefficient as the probability of participants responding “yes” in a given condition, relative to the other condition? Ex: If Exp(beta) for condition A is 3.267, then the probability that participants will say “yes” in condition A is 3.276 times the probability that participants will say “yes” in condition B, all other things being equal? (Is this how one would usually report this kinds of result?) Even if this is accurate, I’m still not sure how to figure out what is going on with the significant 2-way interactions. (Also, since the output only gives parameter estimates for one level of each of the effects. Does this mean I should to re-run the analysis with a different reference category to have a clear picture of the nature of these effects?) 3) I’m unsure about the process of selecting the appropriate model. Is there a standardized procedure for doing this? Just by trial and error I think I have found the best model (using GEE), but I’m not sure how to describe my model selection process when I report it, or whether mine was a valid means of doing so. Basically, I ran a simple model with just main effects, then dropped all the non-significant factors. Dropping the non-significant main effects actually seemed to reduce the fit of the model (why would that be?) so I kept all the main effects in the model, and then added all the 2- and 3-way interactions to see if any interactions were significant. Then I dropped everything except the main effects and the two 2-way interactions that were significant, and compared that final model to the original main-effects-only model. The model that included the two interactions was better (i.e. had a lower QICC and QIC) than the main-effects-only model, so I stopped there. However, one of the main effects that was significant in the main-effects-only model is no longer significant in the final model. I'm not sure why that would be, or what this means. Thank you in advance!


title: Compare classifiers based on AUROC or accuracy?

text: I have a binary classification problem and I experiment different classifiers on it: I want to compare the classifiers. which one is a better measure AUC or accuracy? And why?               Raondom Forest: AUC: 0.828  Accuracy: 79.6667 %                SVM: AUC: 0.542  Accuracy: 85.6667 %


title: What to do with a variable that loads equally on two factors in a factor analysis?

text: After performing a factor analysis on a set of variables, I have one variable that loads equally on two factors.   * What should I do with this variable that loads equally on two factors?   * Should I remove this variable from the factor analysis, and rerun the factor analysis?


title: Correlated error term residual in logit regression: what are my options?

text: I have estimated a model with many interactions of both continuous and factor explanatory variables. The model is to be used for prediction. My model has performed reasonably in out-of-sample testing. However, I have found that my errors are correlated. I took my fitted values of y and ran several regressions like this: $$y\sim \text{fitted}(y)\cdot \text{explanatory variable}$$ There were many instances where an explanatory variable already included in the model was significant. What are some good techniques to resolve this issue? Side note:   I have no strong a priori views on what the "true" model should look like.


title: How exactly is sparse PCA better than PCA?

text: I learnt about PCA a few lectures ago in class and by digging more about this fascinating concept, I got to know about sparse PCA. I wanted to ask, if I'm not wrong this is what sparse PCA is: In PCA, if you have n data points with p variables, you can represent each data point in p dimensional space before applying PCA. After applying PCA, you can again represent it in the same dimensional space, but, this time, the first principal component will contains the most variance, the second will contain the second most variance direction and so on.. So you can eliminate the last few principal components, as they will not cause a lot of loss of data, and you can compress the data. Right? Sparse PCA is selecting principal components such that these components contain lesser non-zero values in their vector coefficients. How is this supposed to help you interpret data better? Can anyone give an example?


title: Statistical Problem

text: A psychiatrist hires you to study whether her clients self-disclosed more while sitting in an easy chair or lying down on a couch. All clients had previously agreed to allow the sessions to be videotaped for research purposes. You randomly assigned 10 clients to each condition. The third session for each client was videotaped and an independent observer counted the clientsʹ disclosures. You reported that ʺclients made more disclosures when sitting in easy chairs $(M = 18.20)$ than when lying down on a couch $[M = 14.31, t(18) = 2.84, p < .05]$, two-tailed.ʺ Explain these results to a person who understands the t test for a single sample but knows nothing about the t test for independent means.


title: Calculating Second-Order Tikhonov Regularization Parameter in Mathematica

text: I am trying to map slowness of underwater sound velocity in a river using some tomographic device. The location of each acoustic receiver/transmitter is shown in the picture below ![enter image description here](http://i.stack.imgur.com/dM7wr.png) To find the slowness m={S11,S12,S13, ...S112,...,S71,S72,S73, ... , S712} from the observations d={t_R1L1,t_R1L2,...,t_R4L4} I need to solve the regularized least-square equation min||Gm-d||^2+a^2||Lm||. I, therefore, established matrices                 G={{0, 12.1524, 2.7609, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 15.6297, 0, 0, 0,        0, 0, 0, 0, 0, 0, 0, 0, 15.6297, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,        5.65726, 9.97241, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 15.6297, 0, 0, 0, 0,       0, 0, 0, 0, 0, 0, 0, 15.6297, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,        14.3434, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 14.912, 0, 0, 0,        0, 0, 0, 0, 0, 0, 0, 0, 15.6283, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,        15.6283, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.50649, 11.1218, 0, 0, 0,       0, 0, 0, 0, 0, 0, 0, 0, 15.6283, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,        15.6283, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 11.0604, 2.88671, 0, 0, 0,       0, 0, 0, 0}, {0, 0, 19.1877, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,        2.16882, 17.9406, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.58072, 15.5287,       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6.99261, 13.1168, 0, 0, 0, 0, 0,        0, 0, 0, 0, 0, 0, 9.40451, 10.7049, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,       11.8164, 8.293, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 14.2283, 0.877225,       0, 0, 0}, {0, 0, 17.3355, 7.70369, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,       10.5772, 15.6647, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2.61622,        18.2809, 5.34477, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 12.9361, 13.3058,       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.97514, 18.2809, 2.98585, 0, 0,        0, 0, 0, 0, 0, 0, 0, 0, 0, 15.295, 10.9468, 0, 0, 0, 0, 0, 0, 0, 0,        0, 0, 0, 7.33406, 14.2199}, {0, 0, 0, 0, 14.4946, 0, 0, 0, 0, 0, 0,        0, 0, 0, 0, 9.24233, 9.34103, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2.41548,        16.1679, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 18.5834, 0, 0, 0, 0, 0, 0, 0,       0, 0, 0, 14.172, 4.41137, 0, 0, 0, 0, 0, 0, 0, 0, 0, 7.34514,        11.2382, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 17.054, 0, 0, 0, 0, 0, 0, 0,        0, 0, 0, 0}, {0, 0, 0, 0, 11.8019, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,        15.1311, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 15.1311, 0, 0, 0, 0, 0, 0,       0, 0, 0, 0, 0, 15.1311, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 15.1311,        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 15.1311, 0, 0, 0, 0, 0, 0, 0, 0, 0,       0, 0, 13.5034, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 2.1671, 10.7843,        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 16.6048, 0, 0, 0, 0, 0, 0, 0, 0, 0,       0, 0, 7.58393, 9.02091, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 16.6048,        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 9.34735, 7.25749, 0, 0, 0, 0, 0, 0,       0, 0, 0, 0, 0, 16.6048, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 11.1108,        1.36223, 0, 0, 0}, {0, 0, 0, 0, 1.29176, 15.5556, 0, 0, 0, 0, 0, 0,        0, 0, 0, 0, 0, 5.29103, 16.3088, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,        4.53791, 17.0619, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3.78479, 17.815,        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3.03167, 18.5681, 0, 0, 0, 0, 0, 0,       0, 0, 0, 0, 0, 2.27856, 19.3212, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,        1.52544, 16.2157}, {0, 0, 0, 0, 0, 0, 18.5392, 0, 0, 0, 0, 0, 0, 0,        0, 0, 0, 19.4035, 2.06375, 0, 0, 0, 0, 0, 0, 0, 0, 0, 19.9029,        1.56439, 0, 0, 0, 0, 0, 0, 0, 0, 0, 20.4023, 1.06503, 0, 0, 0, 0, 0,       0, 0, 0, 0, 20.9016, 0.565678, 0, 0, 0, 0, 0, 0, 0, 0, 0.433035,        20.9679, 0.0663213, 0, 0, 0, 0, 0, 0, 0, 0, 0, 19.7006, 0, 0, 0, 0,        0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 14.1187, 0, 0, 0, 0, 0, 0,        0, 0, 0, 0, 0, 16.3486, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 9.75723,        6.5914, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 16.3486, 0, 0, 0, 0, 0, 0, 0,        0, 0, 0, 4.73928, 11.6094, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 16.3486, 0,       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 14.59, 0, 0, 0, 0, 0, 0, 0}, {0, 0,        0, 0, 0, 0, 1.65553, 11.4627, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,        15.1901, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 15.1901, 0, 0, 0, 0, 0, 0,       0, 0, 0, 0, 0, 15.1901, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 15.1901,        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 15.1901, 0, 0, 0, 0, 0, 0, 0, 0, 0,       0, 0, 7.70533, 3.70495, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0.447913,        15.4943, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 10.2405, 8.21959, 0, 0, 0,       0, 0, 0, 0, 0, 0, 0, 0, 17.5151, 0.944913, 0, 0, 0, 0, 0, 0, 0, 0,        0, 0, 0, 18.4601, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6.32976, 12.1303,       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 13.6044, 4.85561, 0, 0, 0, 0, 0,        0, 0, 0, 0, 0, 0, 15.1623}, {0, 0, 0, 0, 0, 0, 0, 0, 6.77971,        17.4671, 0, 0, 0, 0, 0, 0, 0, 0, 0, 14.9997, 11.4365, 0, 0, 0, 0, 0,       0, 0, 0, 5.00342, 18.2163, 3.21656, 0, 0, 0, 0, 0, 0, 0, 0,        13.2234, 13.2128, 0, 0, 0, 0, 0, 0, 0, 0, 3.22712, 18.2163, 4.99286,       0, 0, 0, 0, 0, 0, 0, 0, 11.4471, 14.9891, 0, 0, 0, 0, 0, 0, 0, 0,        0, 17.4914, 6.76916, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0,        0, 0, 0, 0, 18.2098, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 16.1093, 3.74475,       0, 0, 0, 0, 0, 0, 0, 0, 0, 13.0672, 6.78682, 0, 0, 0, 0, 0, 0, 0,        0, 0, 10.0251, 9.8289, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6.98307, 12.871,        0, 0, 0, 0, 0, 0, 0, 0, 0, 3.941, 15.913, 0, 0, 0, 0, 0, 0, 0, 0, 0,       0, 17.7183, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0,        14.3166, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 15.6093, 0, 0, 0, 0, 0, 0,       0, 0, 0, 0, 0, 15.6093, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 9.15931,        6.45003, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 15.6093, 0, 0, 0, 0, 0, 0, 0,       0, 0, 0, 0, 15.6093, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 11.7252, 0,        0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 2.36708, 11.8831, 0, 0, 0, 0, 0,        0, 0, 0, 0, 0, 0, 15.5369, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 15.5369,       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 14.6006, 0.936286, 0, 0, 0, 0, 0,        0, 0, 0, 0, 0, 0, 15.5369, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 15.5369,       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 12.7613}}      and               L=Flatten[       　　　Table[             M = ConstantArray[0, {7, 12}];,             M[[i, j]] = -4;             M[[i, j + 1]] = 1;             M[[i, j - 1]] = 1;             M[[i + 1, j]] = 1;             M[[i - 1, j]] = 1;             Flatten[M],             {j, 2, 11},             {i, 2, 6}],             1];      Given that `{{U, V}, {Lamb,M}, X} = SingularValueDecomposition[{G, L}]` , where `G=U.Lamb.Conjugate[Transpose[X]]` and `L=V.M.Conjugate[Transpose[X]]`, **Why the generalized singular values of G and L, gama, become indeterminate or zero? And, how can I fix this issue?**               gama=lambda/mu      where               lambda=Sqrt[Diagonal[Transpose[Lamb].Lamb]]           ={0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,          0., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,          0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,          0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,          0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}      and               mu=Sqrt[Diagonal[Transpose[M].M]]       ={0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,      0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,      1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,      1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,     1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.}


title: Predicting data based on two data sets

text: I'm brand new to SAS programming, and I'm having difficulty figuring out this issue: I have two sets of data 1st set has 4 columns (income is a continous range of numbers, sex is 0 for male - 1 for female, and employed is 0 for yes - 1 for no.               id | income | sex | employed      2nd set has 3 columns               id | income | sex       The ids across the two sets do NOT match. What I'm trying to do is create a predictive model of what `employeed` might appear for each `id` in the 2nd set using the 1st set of data as a guideline. After reading through some documentation, I _think_ scoring the data is the way to go, but I have no idea how I would set that up.


title: Posterior distribution of the parameter knowing the prior

text: I have exponentially distributed probability of event $E$ $$P(E|a) = a \exp(-aE),$$ where $a$ is the rate parameter of the exponential distribution. Now the probability distribution for $a$ is a Gamma distribution with parameters $d$ and $b$ $$P(a|d,b) = b^d a^{d-1} \exp(-ba) / \Gamma(d)$$ Knowing the prior distribution over $a$ and the likelihood over $E$, $P(E|a)$, how do I derive the posterior distribution for $a$? I think (not sure!) that it should of the form $P(a|E,d,b)$. I also think that the result should be a standard probability distribution.


title: checking for autocorrelation with many observations and few time periods

text: How would I go about checking for autocorrelation if I had a few thousand observations for each time period and had about 15 different time periods? The data set I am working with has a lag variable which is obviously very influential but I am concerned about the bias in my other predictors when I use this lag variable.


title: Fitting distributions

text: I have some data on serving sizes of a particular food (raw oysters), and I am trying to determine what probability distribution I should select to model these data for a risk assessment (RA). The R package I am planning to use for the RA (mc2d) uses a gamma distribution for serving size in its example, although the package allows you to pick pretty much any distribution you need. I am used to selecting distributions largely based upon the nature of the data. For instance, since I have worked often with ecological community data, I have used the Poisson distribution for rare occurrence events (finding an animal). However, in this case I can't seem to find what people generally use for serving sizes. Does it matter what distribution you use as long as it fits your data? Eyeballing it, the gamma distribution looks like it will fit pretty well if a = 1. But I don't really understand the gamma distribution despite reading about it, and can't piece together in what situations it is used. It actually seems as if the only thing that matters about the gamma distribution is that you select the two parameters to model the shape correctly. Furthermore, I imagine there are other distributions with other set parameters that would fit these data similarly, and it makes me uncomfortable not knowing the difference between them. Edit: for instance, the Erlang distribution also looks very similar to what I need at k=1. How does one decide between two such as these?


title: How to find sufficient statistic?

text: I have trouble finding the following sufficient statistics.   How do you do this?   $$X\sim \Gamma(\alpha, \beta)$$ $$f(x;\alpha, \beta)=\frac{e^{-x/\beta}x^{\alpha-1}}{\Gamma(\alpha)\beta^\alpha}$$ **Question** : Is $log(X_1+X_2)$ a sufficient statistic for beta? I am using the first way, conditional probability, but how to you plug it in and reduce? But there's a log in it. I understand there are two ways to find it:   1. Sufficiency principle: $$P(X_1=x_1, X_2=x_2,...X_n=x_n|T(X_1,...,X_n))$$ does not depend on $\theta$.   2. Factorization theorem: $$P_\theta(x_1, x_2,...,x_n)=h(x_1, x_2,...x_n)g_\theta(T(x_1,x_2,...,x_n))$$ where $h$ only depends on the observation vector $(x_1,...x_n)$ and not $\theta$.


title: Constrained Random Walk Prior BUGS/JAGS

text: I'm currently trying to implement a model along the lines of Owen (2009) and Knorr-Held (2000) in JAGS. In particular, I am struggling to model the constrained random walk priors. We have that $a_{ t }\quad =\quad { a }_{ t-1 }+{ u }_{ t },$ ${ u }_{ t }\quad$ ~ $\quad N(0,Q)$ and also $Q\quad =\quad { \sigma }^{ 2 }(I-\frac { 1 }{ n } { 11 }^{ T })$ I understand that we are centring ${ u }_{ t }\quad$ in order to make the elements of $a_{ t }\quad$ sum to zero but i'm unsure how to implement this in jags. This is what i have come up with so far                   att[i[t]] <- att[i[t-1]] + u[i[t]]         u[i[t]] ~ dnorm(0,Q[i])      Any help would be greatly appreciated!


title: Estimating the mutual information for two signal samples

text: I have two signals _s1_ and _s2_ , sampled 170 times each.               x = 0.2:0.2:34;     s1 = rand(size(x));     s2 = randn(size(x));      The calculation of the MI (mutual information) between two discrete variables requires knowledge of their marginal probability distribution functions and their joint probability distribution. I am estimating each signal's marginal distribution using this Kernel Density Estimator.               [~,pdf1,xmesh1,~]=kde(s1);     [~,pdf2,xmesh2,~]=kde(s2);      I am estimating the joint pdf using this 2D Kernel Density Estimator.               [~,pdf_joint,X,Y]=kde2d(data);      I created a function which takes as input the original signals, their marginal pdfs, and their joint pdf, and computes the Mutual Information. Unfortunately, I seem to have some hug bug in this function, which I can't figure out. The MI should always be a positive number, but I am getting complex and/or negative numbers! The function I wrote for computing the MI is shown below.               function mi = computeMI(s1, s2, pdf1, xmesh1, pdf2, xmesh2, pdf_joint, X, Y)          N = size(s1, 2);     p_i = zeros(1, N);     p_j = zeros(1, N);          for i=1:N         p_i(i) = interp1(xmesh1, pdf1, s1(i));         p_j(i) = interp1(xmesh2, pdf2, s2(i));     end;          mi = 0;          p_ij = zeros(N, N);          for j=1:N         for i=1:N             p_ij(i, j) = interp2(X, Y, pdf_joint, s1(i), s2(j));             delta_mi = p_ij(i,j) * (log2(p_ij(i,j) / (p_i(i) * p_j(j))));             mi = mi + delta_mi;         end;     end;      Thank you very much for your help.


title: How can I test for a difference between ordered groups?

text: Let's say I have a questionnaire about a toy and I ask:   1. The person to rate the toy on a scale of "bad, medium, good"   2. How much the person thinks it cost to make the toy I want to test whether there's a (monotonic) relationship between the two variables. What's the best way of doing this?   1. One way is to run two different one-sided t-tests, comparing (bad, medium) and then (medium, good), but I guess this loses some information.   2. Another way is to use an ordinal logistic regression, where cost is the explanatory variable and the toy rating is the response.   3. Another way is to maybe use something like Spearman rank correlation? This seems kind of funny to me, though, so I'm not sure.   4. Something with (ordered) ANOVA? (I've never really understood ANOVA that well, so I'm not sure.) What's a good approach to try? Bonus 1: What are the pros/cons of the different approaches above? Bonus 2: points for pointing to an R function I could use.


title: Optimizing Principal Component factor weightings over time

text: I was given the returns of a cross-asset class portfolio of ETFs and I conducted PCA to obtain factors on dates from T-n, T-3, T-2,..., T. What I would like to do is decompose the market moves from T+1, T+2, ... onwards into combinations of the PCs. My questions is, what sort of algorithm or optimization method can I use to obtain a set of factor weightings that explain the maximal amount of variance in the market moves going forward. Furthermore, this optimization method should be capable of outputting a list of maxima , not just one, so that the stability of the optimal factor weightings can be assessed over time, since the optimal factor weightings for one day/time period may not be the same as the optimal factor weightings over a larger period.


title: Variance of the sum of correlated random variables

text: i'm trying to compute the variance of the random variable $$X = \frac{1}{N}\sum_{i=1}^N x_i$$ where $x_i$ are correlated identical random variables (mean and variance defined) obtained from a stationary stochastic process. I've reached the point where, with algebra only, I obtain $$\left\langle \left(X-\langle X\rangle\right)^{2}\right\rangle = \frac{\left\langle x^{2}\right\rangle -\left\langle x\right\rangle ^{2}}{N}+2\frac{1}{N^{2}}\sum_{\Delta=1}^{N}\sum_{i=1}^{N-\Delta-1}\left\langle \left(x_{i}-\left\langle x\right\rangle \right)\left(x_{i+\Delta}-\left\langle x\right\rangle \right)\right\rangle $$ I identify $$R\left(\Delta\right)=\frac{1}{\left\langle x^{2}\right\rangle -\left\langle x\right\rangle ^{2}}\frac{1}{N-\Delta-1}\sum_{i=1}^{N-\Delta-1}\left\langle \left(x_{i}-\left\langle x\right\rangle \right)\left(x_{i+\Delta}-\left\langle x\right\rangle \right)\right\rangle$$ to be the autocorrelation function and thus I can write $$\left\langle \left(X-\langle X\rangle\right)^{2}\right\rangle = \frac{\left\langle x^{2}\right\rangle -\left\langle x\right\rangle ^{2}}{N}+\frac{\left\langle x^{2}\right\rangle -\left\langle x\right\rangle ^{2}}{N}2\sum_{\Delta=1}^{N}\frac{N-\Delta-1}{N}R\left(\Delta\right)$$ Can I further simplify this? I mean, does the quantity $$\tau=\sum_{\Delta=1}^{N}\frac{N-\Delta-1}{N}R\left(\Delta\right)$$ means anything in the context of statistics that makes sense to write $$\left\langle \left(X-\overline{X}\right)^{2}\right\rangle = \frac{\left\langle x^{2}\right\rangle -\left\langle x\right\rangle ^{2}}{N}\left(1+2\tau\right)$$ In the (likely) case this is a standard text book calculation, can you provide a reference on a (good) book that presents it?


title: What is geometric autoregressive process?

text: Can anyone give a definition for Geometric Autoregressive Process? Any specific properties? And, in what fields is this mostly applied? To add some context to the question, here is a section of the original Mathematica code that mentions the process:                n = 1000;      \[Lambda] = n \[Gamma];      \[CapitalDelta] = \[Sigma]/Sqrt[\[Lambda]];           p = Table[(1 - z) Exp[t \[CapitalDelta]] + z, {t, -n, n}];            (* (p - z)/(1 - z) follows a geometric AR (1); *)


title: High Ljung-Box p-values at large lags

text: I am trying fit an ARIMA model to stock returns. I have reached a decent model using the AIC criterion. However, the ljung-box p value under a diagnostic plots are pretty weird. The null hypothesis get rejected at higher lags. I tried modifying the parameters, but L-B p value betters only marginally, with a loss in AIC. Any help how I can balance the two ? Also any reasons why the p-value is so low for higher lags. Ihave attached the diagnostic's image: ![tsdiag](http://i.stack.imgur.com/AL5p5.png)


title: Prediction with randomForest (R) and missing values

text: I have a fine `randomForest` classification model which I would like to use in an application that predicts the class of a new case. The new case has inevitably missing values. Predict won't work as such for NAs. How should I do this then?               data(iris)     # create first the new case with missing values     na.row<-45     na.col<-c(3,5)     case.na<-iris[na.row,]     case.na[,na.col]<-NA          iris.rf <- randomForest(Species ~ ., data=iris[-na.row,])     # print(iris.rf)          myrf.pred <- predict(iris.rf, case.na[-5], type="response")     myrf.pred     [1] <NA>      I tried `missForest`. I combined the original data and the new case, shaked it with `missForest`, and got imputed values for NAs in my new case. Too heavy computing though.               data.imp <- missForest(data.with.na)      But there must be a way to use rf-model to predict a new case with missing values, right?


title: When and how to use weights for sequence analysis in social science?

text: ### Weighting in sequence analysis So far, I have scarcely found papers that address the issue of weighting for sequence analysis (using for example the optimal matching algorithm). Sequence analysis normally involves several steps:   1. setting or calculation of substitution and insertion/deletion costs,   2. computation of distance matrices and    3. following cluster analyses or discrepancy analyses[1]. At least, the R-package TraMineR (see Gabadinho et al. 2010 and Gabadinho et al. 2011, p. 11) and the Stata-ado SEQCOMP by Laurent Lesnard makes it possible to include weights at step 1 and 3.   Furthermore, Lesnard explicitly recommends the usage of sample weights for steps 1 and 3: > "Sample weights should only be used to calculate transition matrices, and > consequently substitution costs. Instead of counting the number of > transitions, it is simply the weighted number of transitions that should be > taken into account. The matching procedure in itself, namely, the comparison > of pair of sequences, does not require any weights; it is by definition a > one to one procedure. However, sample weights should be turned on to > interpret results, for instance, if cluster analysis is used, the size of > the clusters obtained must be weighted."   >  Lesnard (2010: 415, endnote 12) ### Open questions Nonetheless, there does not seem to be a consensus in the literature when and which weights are needed or useful.   * What do you think is the best rationale for applying weights in sequence analysis?    * When should sequences be weighted?    * Do you use cross-sectional sampling weights or longitudinal weights accounting for sampling probabilities as well as panel attrition?    * How do you apply weights if you have unbalanced panel data?    * The usage of weights in TraMineR is well documented; but do you have examples for the usage of weights with a Stata-ado?  ### References   * Gabadinho, Alexis, Gilbert Ritschard, Matthias Studer and Nicolas S. Müller (2010): Mining sequence data in R with the TraMineR package: A user's guide, University of Geneva.    * Gabadinho, Alexis, Gilbert Ritschard, Nicolas S. Müller and Matthias Studer(2011): Analyzing and visualizing state sequences in R with TraMineR, in: Journal of Statistical Software, Vol. 40, No. 4, pp. 1-37.    * Lesnard, Laurent (2010): Setting Cost in Optimal Matching to Uncover Contemporaneous Socio-Temporal Patterns, in: Sociological Methods and Research, Vol. 38, No. 3, pp. 389-419.    * Studer, Matthias, Gilbert Ritschard, Alexis Gabadinho and Nicolas S. Müller (2011): Discrepancy Analysis of State Sequences, in: Sociological Methods and Research. Vol. 40, No. 3, pp. 471-510.  [1] See Studer et al. (2011) for a presentation of discrepancy analysis that is an ANOVA like approach for distance matrices.


title: Not all subjects undergo conjunction of two treatments, two-way within-subject (repeated-measures) ANOVA correct?

text: I am wondering about the correct statistical analysis for an experiment involving 2 independent variables (treatment A, 4 levels and treatment B, 2 levels) and one dependent variable (C, score). Because subjects undergo all levels of both treatments, I am assuming a two-way within-subject (repeated- measures) ANOVA (4x2) should be used. However, there is a complicating factor: we used a setup in which we measured C for 2 subjects together, according to the following scheme: \begin{array}{| l | c | r |} \hline & \text{B level 1} & \text{B level 2}\\ \hline \text{A level 1} & \text{Subject 1} & \text{Subject 2}\\ \text{A level 2} & \text{Subject 2} & \text{Subject 1}\\ \text{A level 3} & \text{Subject 1} & \text{Subject 2}\\ \text{A level 4} & \text{Subject 2} & \text{Subject 1}\\ \hline \end{array} When looking at treatment A, all subjects undergo levels 1 to 4, so within- subject design. Also, when looking at treatment B, all subjects undergo levels 1 and 2, so again within-subject design. However, it is not the case that all subjects undergo all levels of A in conjunction with all levels of B, e.g. subject 1 has undergone A level 1 in conjunction with B level1, but not with B level 2. And vice versa for subject 2. Hence, I am inclined to think that regarding interaction between treatment A and B, a within-subject design might not be correct. Any suggestions would be appreciated.


title: Sample size for a very skewed A/B Test

text: I would like to perform an A/B Test on my website. I have basic knowledge on know how to do a basic test statistics, but I'm not sure on how to choose the sample size. In particular, if I have an event with a conversion rate (before the campaign) very very low, let's say p = 1/10^6. To apply the C.L.T. on p for a given interval of confidence, I guess I need a very much larger sample size (for both the control sample and the treatment sample). How can I find the best size of the sample, given that a so skewed distribution of p?


title: Obtaining the outcome hazard ratio? [Survival analysis]

text: I am reading this paper: https://dipot.ulb.ac.be/dspace/bitstream/2013/96837/1/2011-025-ARIAS_DEHON- theroads.pdf On page 23 (of the paper, not the pdf) they have a table showing the "outcome hazard". They say they calculated it using the formula for on page 21 (Equation 3). I don't think I am understanding the formula correctly, so could someone help me how they obtained those values? Using time period 1 as an example, I thought the formula meant the following: e^-1.16/(1+e^-1.16) = 0.2386, which is close to their answer (0.2383), but not quite.


title: Two sided Fisher's exact test with uneven number of samples

text: I'm not a statistician so maybe this question doesn't make sense, but here goes... I'm trying to write some code that does a 2-sided Fisher's exact test on some data that looks like this:               | GROUP | SAMPLE_NUMBER | TOTAL_RESULT | BAD_RESULT |     -----------------------------------------------------     | 1     | 1             | 100          | 0          |     | 1     | 2             | 100          | 2          |        | 1     | 3             | 100          | 5          |     | 1     | 4             | 100          | 1          |     -----------------------------------------------------     | 2     | 1             | 100          | 2          |     | 2     | 2             | 100          | 1          |     -----------------------------------------------------     | 3     | 1             | 100          | 6          |     | 3     | 2             | 100          | 3          |     -----------------------------------------------------     | 4     | 1             | 100          | 7          |     | 4     | 2             | 100          | 9          |      The book I'm reading says that I need to sort the results into duplicate pairs and get the Fisher's P-value for each pair, like so:               | GROUP | PROPORTIONS COMPARED | P-VALUE |     | 1     | 0/100 vs. 2/100      | 0.2487  |     | 1     | 5/100 vs. 1/100      | 0.1203  |     | 1     | 2/200 vs. 6/200      | 0.1758  |     | 2     | 2/100 vs. 1/100      | 0.6231  |     | 3     | 6/100 vs. 3/100      | 0.3337  |     | 4     | 7/100 vs. 9/100      | 0.6153  |      Because group 1 has 4 sets of data, the book says I need to compare the first and second, the third and fourth and then the first two with the last two, hence there are three rows in the above table for group 1. I kind of get this but what would I do if there were an uneven number of samples (e.g. 3) from group 1?


title: Asymptotic vs. bootstrap test statistic, size and power properties

text: I am running a VAR based Granger non-causality tests. I've obtained asymptotic and bootstrap $p$ values for Wald joint test of 0 restriction on a set of lagged variables. It appears that bootstrap based models fail to reject Granger non-causality more often than the models based on asymptotic theory. My question is: How do I evaluate the size and power properties of these approaches/tests? I know how to evaluate the size properties of the bootstrap based tests, but I'm not sure how to approach the asymptotic theory based tests.


title: Standard errors of the MLEs

text: Can anybody tell me how to find numerical values for standard errors of the MLEs of Weibull distribution using the uncensored real data set on the breaking stress of carbon fibres (in Gba) reported by Cordeiro GM, Lemonte AJ. The β-Birnbaum–Saunders distribution: an improved distribution for fatigue life modeling. Comput. Stat. Data Anal. 2011; 55: 1445–1461. The data are (n = 66):               3.70, 2.74, 2.73, 2.50, 3.60, 3.11, 3.27, 2.87, 1.47, 3.11, 3.56, 4.42, 2.41, 3.19,      3.22, 1.69, 3.28, 3.09, 1.87, 3.15, 4.90, 1.57, 2.67, 2.93, 3.22, 3.39, 2.81, 4.20,      3.33, 2.55, 3.31, 3.31, 2.85, 1.25, 4.38, 1.84, 0.39, 3.68, 2.48, 0.85, 1.61, 2.79,      4.70, 2.03, 1.89, 2.88, 2.82, 2.05, 3.65, 3.75, 2.43, 2.95, 2.97, 3.39, 2.96, 2.35,      2.55, 2.59, 2.03, 1.61, 2.12, 3.15, 1.08, 2.56, 1.80, 2.53      I use Mathematica so kindly give answer for Mathematica users.


title: Different fitting models using auto arima and tbats

text: I have one year of daily data for forecasting. while using auto.arima to find the best fit model, it gives me ARIMA(3,1,3). However, when I used tbats to find the best fit model, it gave me the following: BATS(1,{2,2},-,-) I wonder why they provided different models ( If I get it right, the BATS model is an ARMA(2,2)). (ARIMA predictions makes more sense as BATS predictions give negative numbers which does not make sense) I also believe that my data have weekly seasonality but non of these models include any seasonality index. Attached is my timeseries plot. ![enter image description here](http://i.stack.imgur.com/9sC9v.jpg)


title: Goodness of fit index in HMM

text: I am working on a dataset which I am training using the Baum Welch algorithm in order to produce an HMM. Because of the fact that my dataset is rather small I am using Leave-one-out cross validation (LOO). Now, I need to decide an appropriate index for the goodness of fit. The thing is that I cannot decide which is the optimal choice and moreover how can I apply it. For instance, can I use the squared R2, or the sum of log likelihoods, the BIC, or AIC. I've read for all these but I am a bit lost. Should I apply it on the training data? And if yes, how? Perhaps by using the forward algorithm in order to get the predicted emission probability? But if yes, what kind of sequence should I use as input? Which is the real value and which the predicted in my case? Thank you.


title: Multivariate Mixed model in nlme

text: I have two correlated response variables ($y_1,y_2$) explained by the same covariate set ($x_1,x_2$), each with mean = 0. I have a random grouping factor ('group') which is heteroskedastic ($group \sim N(0,\sigma^2_{u1})$ for $y_1$ and $group \sim N(0,\sigma^2_{u2})$ for $y_2$). I would like to model the correlation between the responses through the random error which is also heteroskedastic. $\epsilon_{1} \sim N (0,\sigma^2_{e1})$, $\epsilon_2 \sim N(0,\sigma^2_{e2}), $ and $Cov(\epsilon_{1},\epsilon_2)=\sigma_{e1,e2}$. So, the equations for the $j^{\text{th}}$ observation in the $i^{\text{th}}$ group look like this: $$ y_{1ij} = \beta_{11} x_{1ij} + \beta_{12}x_{2ij} + group_{1i} + \epsilon_{1ij} \\ y_{2ij} = \beta_{21} x_{1ij} + \beta_{22}x_{2ij} + group_{2i} + \epsilon_{2ij} \\ $$ I am using the following code to generate an example dataset. It is followed by the nlme code that I attempted (which gave me an error). I would appreciate any guidance on how to do this right. Thanks!               library(MASS)     library(matrixcalc)     require(reshape2)     library(nlme)               ni <- c(rep(15,6),rep(20,6))          smpl <- sum(ni)          des <- factor(x=rep(x=1:length(ni),times=ni))          Z   <- model.matrix(~des-1,data=des)               mydata <- data.frame(des)          colnames(mydata) <- c('group')               means <- c(0,0,0,0)                                               # Covariance parameter values of Responses and Covariates           y1_sig <- 7           y2_sig <- 80           x1_sig <- 100           x2_sig <- 150                 cov_fe <- matrix(c(                              7*y1_sig     ,  -0.4*y1_sig*y2_sig,    0.75*y1_sig*x1_sig,   -0.65*y1_sig*x2_sig,                              -0.4*y2_sig*y1_sig,       y2_sig*y2_sig,   -0.6 *y2_sig*x1_sig,    0.8*y2_sig*x2_sig,                              0.75*x1_sig*y1_sig,  -0.6*x1_sig*y2_sig,         x1_sig*x1_sig,   -0.5*x1_sig*x2_sig,                             -0.65*x2_sig*y1_sig,   0.8*x2_sig*y2_sig,   -0.5 *x2_sig*x1_sig,        x2_sig*x2_sig                              ),4,4)               set.seed(101)          fe <- mvrnorm( n = smpl, mu = means, Sigma = cov_fe, tol = 1e-6, empirical = FALSE)           mydata$X1 <- fe[,3]               mydata$X2 <- fe[,4]          # random values     u1_sig <- 10     e1_sig <- 8     u2_sig <- 150     e2_sig <- 100          cov_ref <- diag(c(u1_sig*u1_sig ,u2_sig*u2_sig   ))     cov_rer <- diag(c(e1_sig*e1_sig , e2_sig*e2_sig  ))          ref <- mvrnorm(n = length(ni), mu = rep(0,2), Sigma = cov_ref, tol = 1e-6, empirical = FALSE)     reff <- Z %*% ref                        rer <- mvrnorm( n = smpl,  mu = rep(0,2), Sigma = cov_rer,  tol = 1e-6, empirical = FALSE)          mydata$Y1 <- fe[,1] + reff[,1] + rer[,1]         mydata$Y2 <- fe[,2] + reff[,2] + rer[,2]          # Stacking data into a univariate form     resp1 <- data.frame(mydata$Y1, mydata$Y2)     mydata_resp = melt(resp1,  variable_name='Y')          mat <- cbind(mydata$X1, mydata$X2, mydata$group)     # New dataset     mydata1 <- data.frame(direct.prod(diag(2),mat))     colnames(mydata1) <- c('X11','X21','group1','X12','X22','group2')          mydata1$variable <- mydata_resp$variable     mydata1$value <- mydata_resp$value     mydata1$group1 <- factor(mydata1$group1)     # Create binary variable to differentiate between stacked responses     mydata1$D       <- as.integer(mydata1$variable == "mydata.Y1")          # nlme call     mdl.nlme <- lme(fixed = value ~ 0 + X11 + X21 + X12 + X22 ,                     random = ~ -1| D /group1,                     correlation = corSymm(form=~1|D),                     data=mydata1)


title: Using random sampling, is it better for the samples to be disjoint?

text: Background I would like to measure the performance of a model trained on 3k samples, because this number of samples might be feasible to obtain in practice. I have a larger set of samples to choose from (about 25k), but generating the model is expensive in terms of time and computational effort, so I can only do about 5 runs. In each run, I will pull 3k samples for training, and test on the remaining samples. **Question** Should the five sets of 3k samples used for training be disjoint? In other words, will this produce a better estimate of this model's performance than five sets of 3k samples selected independently at random? Note that in neither case is it possible to use all samples for training at least once.


title: How do you call this type of "regression" and how to do it?

text: I have got data that looks like               HAIR COLOUR  SKIN COLOR  NUMBER     "black"      "brown"     1101     "red"        "brown"     102      "yellow"     "brown"     103     "black"      "yellow"    1104     "red"        "yellow"    105      "yellow"     "yellow"    106     "black"      "white"     1107     "red"        "white"     108      "yellow"     "white"     109      (in the reality, there are more factors, more values of each factor and it looks more messy). I need to explain the numbers based on the factors, e.g. all large numbers are when HAIR COLOR = "black", and, if the impact of HAIR COLOR is subtracted, it is greated for "white" and smaller for "brown". Or, in any way, explain the numbers in a summarized way. I wonder if there is a name for this type of analysis, and how to carry it out (maybe there are some standard techniques)? N.B. The numbers mean nothing, I have made up this table just illustrate my question.


title: Expected value and variance of trace function

text: For random variables $X \in \mathbb{R}^h$, and a positive semi-definite matrix $A$: Is there a simplified expression for the expected value, $\mathop {\mathbb E}[Tr(X^TAX)]$ and variance, $Var[Tr(X^TAX)]$? Please note that $A$ is not a random variable.


title: Confused about which variables should be dependent and independent in this MANOVA

text: I am currently doing research on social media and how it can effect depression, anxiety and stress scale (DASS) scores. I have two age groups (18-25, 26-35) and also took scores on narcissistic personality inventory (NPI). I'm struggling to know what statistical test to do. I was thinking MANOVA with age as my independent variable and time spent on social networks, NPI scores and DASS scores as my three dependent variables. Is this correct? My research question is: Does social network use affect DASS scores and NPI scores in two age groups?


title: Credibility Intervals

text: I'm trying to understand when credibility intervals are useful? Are there examples of real world situations where credibility intervals are the better thing to use compared to confidence intervals? Note that by "useful", I mean maximizing some concrete real-world objective (so not for instance trying to get posterior intervals for one's subjective beliefs which can be useful but not something I'm looking for) Thanks


title: How to avoid multicolinearity in SVM input data?

text: Do you know of any techniques that allows one to avoid and get rid of multicolinearity in SVM input data? We all know that if multicolinearity exists, explanatory variables have a high degree of correlation between themselves which is problematic in all regression models (the data matrix is not invertible and so on). My question is actually a bit more subtle as relevant data has to be selected. In other words, multicolinearity must be avoided while keeping relevant input data. So in a way, I guess that my question is also about reducing the input data matrix dimension. So, I've thought about PCA, which permits both reducing dimension and getting uncorrelated vectors (the PCs), but then I don't know how to deal with the eigenvector's signs.


title: Should I start with correlations first before doing any testing?

text: I'm using open source data to analyze innovation in a business context for my thesis. I've got the questions from the survey in an Excel file and although there's only 13 questions there's over 600+ variables and it is confusing for a beginning SPSS user like me. I want to see if there's a relationship between "innovation optimism" and the "environment" of a firm's home country (Japan). Is it best to start with looking for correlations? I'm not sure how many variables I can select at one time to test this. (i.e., if bivariate correlations only let me test 2 variables at a time how do I test multiple ones to see if there's a correlation?) Also, how do I select only Japan and not all the 22 countries surveyed?


title: What's the difference between Maximizing Conditional (Log) Likelihood or Joint (Log) Likelihood while estimating parameters of a model?

text: Consider a response **y** and data matrix **X**. Suppose I'm creating a model of the form - **y ~ g(X,$\theta$)** (g() could be any function of X and $\theta$) Now, for estimating $\theta$ using Maximum Likelihood (ML) method, I could go ahead either with Conditional ML (assuming I know the form of conditional density **f(y|X)** ) or with Joint ML (assuming I know the form of joint density **f(y,X)** or equivalently, **f(X|y) * f(y)** ) I was wondering if there are any considerations in going ahead with either of the above two methods apart from the assumption about the densities. Also, are there any instances (specific types of data) where one method overpowers other most of the time?


title: Possible causes for the state noise variance to become negative in a Kalman Filter?

text: I am having some trouble debugging an application of a linear discreet Kalman Filter. From time to time, I find that there are diagonal elements of the covariance matrix that become negative. This is clearly erroneous and causes the whole filter to blow up. What are the possible causes and remedies for this problem? I know there is a known issue with rounding errors if the error measurement update is expressed as $P^t = (I - K_tH_t)P^{t-1} $ and so I have implemented this, as suggested in several sources, as $P^t = (I - K_tH_t)P^{t-1}(I-K_tH_t)^T + K_tR_tK_t^T $ which is supposed to be the sum of two positive definite matrices, ensuring that $P^t$ is positive definite. In my case, this isn't so. The elements go negative not by small, rounding errorish values, but by comparable values to the positive elements. My Kalman Filter implementation has been pretty battle hardened by this point so I don't think there are any glaring bugs. I haven't encountered this problem before. I am playing around with a new dataset which is causing these errors. The main suspect at this point that I can come up with is that fact that at times, entire columns of the $H_t$ matrix can be zero. Recall that $H$ is defined as $z_t = H_tx_{t-1}$ where $z$ are the measurements and $X$ the state vector. As I say, sometimes whole columns of $H$ are zero, meaning that there are no measurements that contain any information about some particular state in the measurements for the given timestep. This may be an issue, but on the other hand I would hope that the filter would simply leave that state unchanged, and add a bit of process noise to its variance, without having any major dramas. I feel like this situation must have happened before with some other datasets without any ill effect. Where else could I look to see where the problem is emerging? The $H$ matrix contains negative values, but that should not be an issue. $K$ also is negative at times, however following through the logic of the Filter construction, the negative elements in $H$ that led to negatives in $K$ should cancel back out in the error measurement update step. Clearly that isn't happening though. Are there any standard approaches to ensuring the robustness of a Kalman Filter that I need to employ here, or any other suggestions on where to look to find the cause?


title: Approximate distribution of normal squared

text: I am studying for a test, one section of which will cover the delta method. This problem came from that section: Let $X\sim N(\mu,n^{-1})$. Find an approximate distribution of $X^2$. (It also asks for exact distribution, but I can do that). Can you help me get started? I can't see how I can apply delta method or central limit theorem etc. to find some approximation. I know that $\sqrt{n}(X-\mu) \sim N(0,1)$, but I don't see how that can help either.


title: How to extract factor scores for each observation in IRT?

text: library(ltm)     fit <- rasch(LSAT)     factor.scores(fit)      This code can generate latent trait scores for response patterns. How can I extract latent trait scores for each observation?


title: How is the support in the Apriori algorithm calculated in the case of duplicates?

text: I am trying to implement Apriori algorithm. However I have a small doubt what to do when the same item appears more than once in the one basket,I have 2 transactions say               T1 = {A,A,C}     T2 = {A,X}      What is the support of A ? Is it 3 or 2 ?


title: Likert Scales multiple answers and items

text: I am trying to introduce my data into SPSS. At the moment i am very confused. My questionnaire comprises 8 statements, where participants have ranked from strongly disagree to strongly agree their opinions. so far so simple. the only problem is, that for each of the 8 statements there are 8 characters that the participant has to rank. Could you please let me know how can i introduce such multiple items Likert scales?


title: A stochastically increasing exponential family for which $\lim_{\theta\rightarrow\inf\Theta}\mbox{P}_\theta(X\leq x)\neq 1$

text: ## Question A little something that I've been wondering about for a while: Let $P_\theta$ be a stochastically increasing (one-parameter) exponential family on the sample space $\mathcal{X}$ with $\Theta\subset\mathbb{R}$ being its natural parameter space, i.e. $\Theta$ being the set of values for which the cdf $F_\theta$ defines a probability measure. Is it always true that $$F_\theta(x)\nearrow 1\qquad\mbox{as}\qquad\theta\searrow \inf\Theta,$$ and $$F_\theta(x)\searrow 0\qquad\mbox{as}\qquad\theta\nearrow \sup\Theta$$ for all $x\in\mathcal{X}$ which are not in the boundary of $\mathcal{X}$? * * * ## Definitions, an example, speculation A distribution $P_\theta$ over $\mathcal{X}$ parametrized by $\theta\in\Theta$ is stochastically increasing if, for fixed but arbitrary $x\in\mathcal{X}$, $F_\theta(x)$ is decreasing in $\theta$, where $F_\theta(x)=\mbox{P}_\theta(X\leq x)$. An example is the binomial distribution, where $\mathcal{X}=\{0,1,2,\ldots,n\}$ and the natural parameter space is $\Theta=(0,1)$. $$F_\theta(x)=\sum_{k\leq x}{\binom{n}{k}}\theta^k(1-\theta)^{n-k}$$ is decreasing in $\theta$. Example: ![](http://i.stack.imgur.com/15FoB.png) In this setting we have $$F_\theta(x)\nearrow 1\qquad\mbox{as}\qquad\theta\searrow \inf\Theta=0,$$ and $$F_\theta(x)\searrow 0\qquad\mbox{as}\qquad\theta\nearrow \sup\Theta=1$$ for all $x$ not on the boundary of $\mathcal{X}$. At the boundary, i.e. when $x\in \{0,n\}$, only one of the limits is attained. Note that if we restrict the parameter space to be a proper subset of $(0,1)$, such as $(0.2,0.8)$, this is no longer true: $F_\theta(5)\nearrow F_{0.2}(5)=0.174\ldots<1$ as $\theta\searrow \inf\Theta=0.2$. This property seems to hold for all commonly used exponential families, so my guess would be that if a counterexample exists it has to be somewhat pathological. It could perhaps involve functions that cause $\exp(<\theta,T(x)>)$ to blow up for some finite $(\theta,x)$, making its integral infinite. In some sense this would make the parameter space "restricted". A "trivial" example is perhaps the Bernoulli(p) distribution, since its sample space equals its own boundary, so that only one of the limits can be attained for each point in the sample space. But that's a wee bit boring, and I'd much rather have example where at least one point of the sample space is not in the boundary.


title: Checking normality of a Likert scaling question

text: I have likert scale responses (1-5 rating) of 20 respondents for a set of 25 questions. So that means that each respondent has answered 25 questions (with 1-5 rating). If I need to check the normality distribution of all the 25 questions together for a particular respondent, how do i do that?


title: Monte Carlo integration help needed

text: I'm trying to simulate these two integrals using Monte Carlo simulation: $$ \int_{-\infty}^\infty \exp(-x^2) dx, \quad \mbox{and } \int_{-\infty}^\infty \exp(-|x|) dx . $$ When I use `runif(n,-Inf,Inf)`, I get `NaN` (number very close to zero) (using R). I also tried converting that to a polar coordinates integral ($drdθ$ integral) but then I don't really know how to do the double integration needed. Any ideas?


title: Finding the probability of normality via Anderson-Darling, Shapiro-Wilk , and Kolomogrov-Smirov

text: I have quite a few distributions that were generated by some system, and I am trying to find the probability of accordance to normality for these distributions. So, the probability of accordance to normality could be found as: # P(p>0.05) = ( #p(AD test)>0.05 )/ total_no_of_distributions. where p(AD test) refers to the AD test p-value resulted from fitting from fitting one distribution to normal, ( #p(AD test)>0.05 ) is the sum of all p's greater than 0.05, and P(p>0.05) is the probability of accordance to nomrmality. My issue is that I intend to randomly try several methods, e.g., AD test, KS test, SW test, or LF test, to test the normality of the distributions that I have. A colleague has told me that this method might not be accepted because it is based on several methods and not one, which is quite orthodox. My idea is to propose a method that makes use of all the available normality testing methods, since it is difficult decide which one is more powerful. I appreciate any comment and feedback to this issue. Regards, Rau


title: Meaning of a neural network as a black-box?

text: I often hear people talking about neural networks as something as a black-box that you don't understand what it does or what they mean. I actually I can't understand what they mean by that! If you understand how back-propagation works, then how is it a black-box? Do they mean that we don't understand how the weights that were computed or what?


title: Is paired t-test valid for half-normal data?

text: ok i am not a statistician. I am doing paired comparison for several outcomes of cross validated interpolation results. I have to consider only the absolute value in my analysis and so although original data follows a normal distribution the resultant absolute values are only half normal. So is the t-test still valid and sound and in either case (yes/no) what is the reasoning behind it ?


title: Fit and predict logistic regression model for censored data?

text: Suppose I have a data frame such that:               set.seed(2014)     df<-data.frame(y=rbinom(100,1,0.3),futime=as.integer(rnorm(100,100,10)),          age=rnorm(100,50,5),sex=rbinom(100,1,0.5))     df[df>100]<-100      where `y` is the outcomes and `futime<100` are right censored of follow up time. I would fit a logistic regression model to predict the outcomes as if they are not censored. I saw there are some options such as `tobit` model, `survival` model and `censReg` package, but seems not very similar situation. Would somebody know the best way for this using `age` and `sex` as predictors?


title: How can any statistics be calculated from a infinite population.

text: My understanding of some statistics: For a given experiment a finite number of samples can be taken, defining the sample size. However this experiment may have an infinity large population size. Given that in an infinitely large population size any possible outcome will occur an infinite number of times, how can we truly say one outcome is more likely than another. Mathematically lets say from our sample we've found that there is a 0.1 chance of drawing a blue ball, and a 0.0000001 chance of drawing a red ball. If the balls are always replaced these odds do no change. However since the ball can be drawn any number of times this will happen: blue ball count: 0.1 * x red ball count: 0.0000001 *x where x is any number. Clearly the odds do not change with increasing x. But if x is infinitely large then: blue ball count: 0.1 * infinity = infinity red ball count: 0.0000001 8 infinity = infinity So how can we confidently state any odds when the population size is infinity? P.S. This is just out of interest so I would appreciate as much detail, links etc. as possible


title: Two normal distributions

text: I will mark this as homework, even though it's just general interest. The scenario is a little silly, because it's just a real world case for a theoretical problem I've been thinking about. A man enters a hot dog eating contest every year and in the contest he has to eat 10 hotdogs in 5 minutes (yes it is different to normal hotdog contests): As soon as the man eats his 5th hotdog (half way to finish them) the time is recorded, and when he eats his final hotdog the time is recorded again. SO we have two sets of data that can both be approximated with normal distributions, i.e. how long it take sot eat 5 hot dogs and how long to eat 10. Call the first mean and std dev $\mu_{h}, \sigma_{h}$ (h for halfway) and the second distribution's mean and std dev $\mu_{f}, \sigma_{f}$. Both means are in minutes. This year the man enters the contest and it is seen that he takes $\mu_{h} +\epsilon$ minutes to finish the fifth hotdog, for some $\epsilon > 0$, i.e. he is slower than usual. My question is, what time can we expect the man two finish the contest, and how to calculate this? Thinking about it, I thought it would be best to consider the two distributions as a bivariate normal distribution. So if we have the first random variable $X$, as the time taken to get halfway, and $Y$ as the time taken to finish, I want $E(Y|X =\mu_{h} +\epsilon )$. I don't know if this is correct thinking, please correct me if I am wrong. ALso I don't know how to go about calculating the above expected value. Any help is welcome.


title: What is the reason why we use natural logarithm (ln) rather than log to base 10 in specifying function in econometrics?

text: What is the reason why we use natural logarithm (ln) rather than log to base 10 in specifying functions in econometrics?


title: Would it make sense to do a 2SLS regression for Nominal GDP and Money Supply?

text: Would it make sense to do a 2SLS regression for Nominal GDP and Money Supply (using predetermined variables of government expenditure and investment)? If it does make theoretical sense, is this because GDP and Money Supply are endogenous? I am unsure why GDP is endogenous... For example if there is an increase in money supply, then I know there is an increased in GDP (since AD increases); however, I fail to see how an increase in GDP has a feedback effect and causes an increase in money supply.


title: Setting PyMC model with two different time series data

text: I've been working with PyMC for a bit, and am stuck on this one. I see example on fitting time series, in the tutorial and others like: http://lighthouseinthesky.blogspot.com/2011/10/curve-fitting-part-5-pymc.html and http://healthyalgorithms.com/2010/10/19/mcmc-in-python-how-to-stick-a- statistical-model-on-a-system-dynamics-model-in-pymc/ However, in each case, we either have a single variable, or data for two or more variables for the same time points (so the vectors are all of the same length). What happens if my data is like:               t=[1,2,3]     x=[.5,.8,.2]      and               t=[1.5,2.5]     y=[.1,.5]      with some simple model, like               x=a*t+b     y=c*t+d      the model isn't important. the point is that the "observed" data have different time points. I haven't been able to find any example like this, without perfectly aligned data. Are there any? Is there a simple example that someone could point out? thanks!                Brian Blais


title: clustered-stratified random sampling

text: i've got questions about clustered-stratified random sampling. let's say I want to do a research in narcotics. the population that I will be working with is a combination of academics, students, and government employees (law enforcers, DEA, etc). I consider the population's unknown because I couldn't get the exact number of the population. from that population, I'm thinking of clustering the population into two groups, those from academics and the practitioners (the government employees). From that clusters, I developed two levels, the academic cluster contains lecturer group and students group. the practitioners group contains high-level group and low-level group. Have I already clustered and stratified my sample in right manner? How do I decide the exact sample for the cluster without knowing specific population of each group and without knowing total population? and how do I decide the exact sample for each level of the cluster? thank you very much, I appreciate your time and suggestions :)


title: Regression model where errors are not normally distributed

text: From a Physics equation I have the following model: $$W=\beta_0+\beta_1Z$$ $\beta_0$, $\beta_1$ are fixed values for which I want to find a $1-\alpha$ confidence region. I have $(z_1,w_1),(z_2,w_2),\ldots,(z_n,w_n)$ observations. I get these observations from evaluating the function $f: \mathbf{R}^{l+k+2}\rightarrow \mathbf{R}^2$. $$f\left(X_1,\ldots,X_l,Y_{1i},\ldots, Y_{ki},P_{1i},P_{2i}\right)=\left(f_1\left(\ldots\right),f_2\left(\ldots\right)\right)=(Z_i,W_i),$$ where $X_j\sim\mathcal{N}(\mu_j,\sigma^2)$ $Y_{ij}\sim\mathcal{N}(\theta_{ij},\sigma^2)$ $P_{ik}\sim\mathcal{N}(p_{ik},\frac{p_{ik}}{1000})$; $\sigma^2$ is known, but not $\mu_j,\theta_{ji},p_{ji}.$ $X,Y,P$ are measures from an electronic instrument. The equation $w_i=\beta_0+\beta_1z_i$ holds when $\left(z_i,w_i\right)=f\left(\mu_1,\ldots,\mu_l,\theta_{1i},\ldots,\theta{ki},p_{1i},p_{2i}\right)$ ($\bf{not}$ when $\left(z_i,w_i\right)=\mathbf{E}(f(\ldots))$ Note that there are both common $\left(X_j\right)$ and not common variables $\left(Y_{ji} , P_{ki}\right)$ for all observations. Still It happens that $Y_{ji}=Y_{st}$ and $P_{1i}=P_{2t}$ for some $j$ $i$ $s$ and $t$. (Meaning that variables $Y$ and $P$ are used for more than one observation, but not all) $f_1$ and $f_2$ are rational functions. Both numerators and denominators are quadratic functions when restricted to one variable. They are well defined. So, **what can i do to get confidence regions for $(\beta_0,\beta_1)$?** A regression model and Maximum Likelihood techniques are tempting. \begin{equation} \begin{pmatrix} w_1 \\ \vdots \\ w_n \end{pmatrix} = \begin{pmatrix} 1 & z_1 \\ \vdots \\ 1 & z_n \end{pmatrix} \begin{pmatrix} \beta_0 \\ \beta_1 \end{pmatrix} + \begin{pmatrix} \epsilon_1 \\ \vdots \\ \epsilon_n \end{pmatrix} \end{equation} $\bf{\epsilon}=\bf{W-\beta_1Z-\beta_0}$ (Is this right?) As $(Z,W)$ has a complicated joint distribution, so does $\bf{\epsilon}$. I can figure out what to do when each $(Z_i,W_i)\sim\mathcal{N_2}((z_i,w_i),\sum)$. But this is not the case. I understand that there may not be an easy answer for these questions, but if anyone can give some general tips or point me to the right literature i can work from there. Thanks for reading.


title: 4 point and 5 point Likert scales in a new questionnaire

text: I'm looking to confirm whether having a 4 point Likert Scale for some constructs, and a 5 point Likert scale for other constructs will be suitable when developing a questionnaire? I have 7 constructs with a 4-point scale, then 2 constructs with a 5 point scale and one construct with a different 4-point scale Will this affect validation process when testing for construct validity? And which analysis would be best - e.g. polychoric correlations? Thank you.


title: Find minimal set of variables for regression (PCA?)

text: Given a large matrix with 10,000 rows (variables) and 20 columns (sampling timepoints), I am trying to build a linear model to find the sampling timepoint given some of the variables. So in essence, I want to find a small subset of variables that explains the data well enough that I can build a linear model like `timepoint ~ v1 + v2 + v3 + v4`. I hope that this will allow me to only take measurements of the few variables that actually make a difference instead of measuring all 10,000 of them every time. I have already tried using PCA for this, by doing a PCA on the whole matrix and then using the PCs as my variables. This works to some degree, but it does not actually solve my problem, since I would still have to measure all the variables for this to work. When I look at the loadings for the different PCs I usually get a large number of variables for each of them, so I cannot use this to determine the important variables either. Is there a way I can do this?


title: Is it possible to predict/forecast results of pairwise comparison matrices based on one or two initial samples?

text: I'm curious if it's possible to predict the results of a pairwise comparison matrix that will be sent out to 100 people based on the results of one or two filled out surveys? **The matrix has the following format:**               Indicators    Var_1    Var_2    .....    Var_n     Var_1         NA     Var_2                  NA     .....                           NA     Var_N                                    NA      There are three ways to to fill this survey:   1. column more important   2. Column and row are equally important   3. row more important I have received 2 filled out surveys so far and just wondering if from those 2 results I can predict how the results will look like for the 100? Thanks, D


title: How to include level-2 variables in a random intercept model when fitting HLM with lme4?

text: I want to include 2 covariates that are group characteristics as fixed effects in predicting a variable that is at the individual level. How to include them in the lmer formula? At this point my random effects are given just by the intercept and my fixed effects are only individual level variables.               a <- lmer(DV ~ 1 + v1 + v2 + (1 | Team), data = dat)      Should they they be included as random effects?


title: How does simulation work

text: I am trying to use the Zelig package to create a linear model and perform 'what if' analysis and run a simulation to determine what the quantiles will be if give an value. My data frame looks like this:               dput(head(mydata, 10))     structure(list(DATE = structure(c(1367409600, 1367413200, 1367416800,      1367420400, 1367424000, 1367496000, 1367499600, 1367503200, 1367506800,      1367510400), class = c("POSIXct", "POSIXt"), tzone = ""), Tot = c(41.19,      48.04, 58.76, 58.92, 63.96, 39.2, 47.9, 52.83, 54.89, 56.12),          VOLUME = c(17447212L, 22345495L, 24855225L, 24756301L, 23329007L,          15601614L, 20119332L, 21975461L, 22538957L, 21344801L)), .Names = c("DATE",      "Tot", "VOLUME"), row.names = c(NA, 10L), class = "data.frame")      I ran these steps:               z.out <- zelig(Tot~VOLUME, data=mydata, model="ls")     x.out <- setx(z.out, VOLUME=20000000)     s.out <- sim(z.out, x.out)      `summary(s.out)`, reveals the following:               Predicted Values: Y|X          mean    sd    50%   2.5%  97.5%       64.872 0.489 64.886 63.902 65.801      Is 97.5% a percentile? 97.5% of the time predicted value will be within 65.801? I understand the linear model and what if analysis. But I really don't get how simulation works. Can somebody explain to me in layman's terms how simulation works? Looks like sim function estimates the parameters based on "maximum likelihood models using asymptotic normal approximation to the log-likelihood", anybody knows what that is and any good examples? given the model and independent variable value, how is it that using MLM, we get different result every time, any ideas or great examples? This is what docs say:" Simulate quantities of interest from the estimated model output from zelig() given specified values of explanatory variables established in setx(). For classical maximum likelihood models, sim() uses asymptotic normal approximation to the log-likelihood." Any ideas what does this mean?


title: The gamma distribution and Poisson processes

text: I know that the gamma distribution with parameters $k$ and $\theta$ can be used as a model for the occurrence of events. The requirement on the events would be that their occurrence is random and the mean time between them equals $\theta$ and that they belong to a Poisson process. Now I have a question regarding the applicability of these: I want to model a biological growth process where we think of the resulting length as random. What we do know is that the underlying processes (biochemical etc.) are Poisson processes. Is it in the context of the statement above appropriate to model the length with a gamma random variable?