title: Does SPSS give factor scores for mean subtracted variables?

text: From the Johnson's book, I can see the orthogonal factor model is as follows: $X-\mu=LF+\epsilon$ So if I use SPSS for extracting factor scores from a factor analysis of several variables, then am I going to get factor scores for mean subtracted variables? Actually my intention is to use the factor scores as IVs along with some IVs that were not used in the factor analysis, which is, as far as I know is a valid procedure. But if I get estimated values of $F$ (factor scores) from $(X-\mu)$, then what should be done to the IVs which were not included in the factor analysis during regression? What should I do with the DV? Should they also be mean subtracted?


title: Thesaurus for statistics and machine learning terms

text: Does there exist any reference thesaurus for statistics and machine learning terms? I know that Wikipedia articles often contain synonyms, but I would like to have a mere thesaurus that I could go through easily (vs. a full encyclopedia) to make sure I know all the jargon.


title: On what basis do we differentiate a standard normal distribution from normal distribution?

text: The normal distribution is a prerequisite for ANOVA. It is not clear what we are supposed to understand by "standard normal distribution". Does ANOVA need a _standard_ normal distribution? The answer can be given in the context of meta-analysis of sample correlation.


title: Test is data is derived from mixture or just noisy normal

text: Is there a way to quantify \ test \ describe how likely it is that high dimensional data come from a single Gaussian or a mixture of Gaussians (with different means \ SDs) or not? What I am thinking of is similar to the idea of a t-test, but when the separation into groups is unknown.


title: Conditional Probability With Naive Bayes

text: I have recently begun exploring R I have downloaded a data set that contains flight times for when a plane takes off at Origin and Lands at destination The data has up to three months’ worth of information and looks something like this Plan_ ID, Origin_Airport, Destin_Airport, Sched_Dep, Sched_Arriv, Act_Dep, Actual_Arriv I changed the data so it looks somewhat like the below Dep_Day_Month: FRI-3 Movement: AAAAAA (ORKORK - DUBDUB) Num_Flights 1 OT_Dep 1 OT_Arr 0 OT_Dep_Per 1.0000000 OT_Arr_Per 0.00000000 So for three months’ worth of data I split it into Days so Dep_Day_Month can consist of Mon-3 to Sun-3 which is Monday to Friday of the month of March. This Monday to Friday breakdown exists for 3 months. I also added an additional Monday to Friday for the previous week The Movement is a concatenation of the flight number and the origin Destination. Num of flights is the number of flights on the route/month/day. OT Dep is the number of flights that left on time (so any time before or equal to the on time departure). OT_Arr is the count of the number of on time flights. The last two columns are a percentage of this The final aim of the experiment is to take all the Mondays and place them into a table, all the Tuesdays into a table, all the Wednesdays…., then predict if the flight will arrive on time based on its past performance of the lane maybe using a Bayesian model. Basically it’s a very simple form of routing analysis (for one leg at the moment..always ORKORK to somewhere else) I have two questions I’m hoping to get help with J 1) Bayesian Models require that the input all be independent in order to work effectively (I have read research papers where it has good results on non- independent variables but that’s for another time). This I think can be achieved if I combine just the OT_Dep_Per and OT_Arr_Per into a single column. I think this can be done using Dependant probability. So for example what is the probability flight A will arrive on time given Flight A Departed on time. What im not clear on is how to implement this if for example Flight A arrives on time if it departed late because pilots can generally catch up in the air if the wind is good and they floor it) 2) Do you think it’s a good way to approach the problem? Splitting the flight performance up by day by Month for the variables adding the previous week to get the most recent trends Thank you very much for your help


title: Poisson arrival

text: Users arrive according to Poisson process with rate $\lambda$. If every third user is removed, then do the remaining users form a Poisson process with rate $2\lambda/3$? If every other user is removed, then do the remaining users form a Poisson process with rate $\lambda/2$? and so on..?


title: How can cost decrease by more than 100 percentage points?

text: I saw the following graph today in the Atlantic (link), but I can't understand how the cost of televisions can go decrease by more than 100 percentage points in the last ten years. How is that mathematically possible? Or did the authors make a mistake? ![enter image description here](http://i.stack.imgur.com/VlG8C.png)


title: Autocorrelation of convolution integral

text: Work out the autocorrelation $r_Y(\tau) = E[Y(t)Y(t+\tau)]$ with $Y(t) = \int_{-\infty}^{\infty} h(t-u) x(u)$ and $X$ a WSS, ergodic process I always get: $h(t)* h(t+\tau) * r_X(\tau)$ (with $*$ convolution) ## My approach (after multiplication) $\int h(t-u) ( \int h(t+\tau -u')r_X(u'-u) \text{d}u' ) \text{d}u$ The inner integral is the definition of a convolution hence $= \int h(t-u) ( h(t+\tau)*r_x(t+\tau-u) ) \text{d}u$ This gives again a convolution integral with time variable t. so, $= h(t)* h(t+\tau) * r_X(\tau) $ ## Correct solution: $ r_X(\tau) * h(\tau) * h(-\tau) $ ## Question What am I doing wrong?


title: Methods for measuring snowball effects in a "complete" longitudinal dataset

text: I'm looking for ways to test for "cumulative advantage" effects in a longitudinal dataset (see image) ![Data set example](http://i.stack.imgur.com/W3LQ2.jpg) I guess the data set is principally similar to this: http://www.caldercenter.org/whatis.cfm , but to be precise, the figures are weekly marked shares of a fixed selection of products. The full data set contains nearly 25,000 products tracked over 76 weeks. What I want to investigate is:   1. Whether the rate of growth in marked share is a function of current values of marked share.   2. Whether small "advantages" at an early stage grow larger over time. I have basic command of SPSS and methods like correlation and regression, and have some rough ideas like e.g. for each week "manually" measure the correlation between the (accumulated) marked share and the marked share of the following week, but I'm quite sure better suited approaches exsist. Could someone please be so kind as to point me to possibly fruitful directions?


title: General reference to Rényi entropy

text: Is there a book or a review article serving as a good general reference to Rényi entropy, its applications and related concepts?


title: (When) Does simulating a larger sample from a small sample yield better results?

text: I am trying to classify a set of $p$ predictors into 5 classes. But my sample size $n$ is rather small, so I fear I won't get a very robust estimate. Now an idea would be to subset my data for each of the 5 response classes, and simulate more data in each class. I would e.g. assume a multivariate normal distribution for the predictors, and then (I think some people call this a _parametric_ bootstrap) estimate $\mu$ and $\Sigma$ from the multivariate normal, and with those parameters then simulate 10000 new observations in this class. I would do this for every class and then have 50000 additional observations. Given that I simulate the data from only the knowledge I have from my small sample, is there anything I can gain with this procedure? Robustness yes, but I will certainly not get more information, right? Does this approach make sense at all? Is there maybe already a better way to solve this problem?


title: difference in training and testing procedure of model

text: Can anyone please tell me the difference in training and testing of a model. I have developed 5/6 different single pass online learning algorithm (ets, ets+, evolving fuzzy modelling, SOFNN, Gustafson-Kessel, participatory learning). I have collected 3000 pts of data for training and 500 points of data for testing from MG time series in *.mat file in simulink. How do the training and testing procedures differ?


title: Contrasts in mixed model

text: Suppose I have a data frame like this (four groups, for each n = 12):               > head(my_data,n=8)          Subject Group Condition    Value     1       1     A    Cond_1 12.40407     2       1     A    Cond_2 14.89856     3       2     B    Cond_1 13.87088     4       2     B    Cond_2 14.31305     5       3     C    Cond_1 13.40773     6       3     C    Cond_2 13.48016     7       4     D    Cond_1 13.76183     8       4     D    Cond_2 12.60769     ...     ...      I want to examine if there are differences in `Value` between the four groups A, B, C, D, considering also the effect of `Condition`. So I have to do a mixed model. Before, I set specific contrasts:               > AvsB <- c(1,-1,0,0)     > AvsC <- c(1,0,-1,0)     > A vsD <- c(1,0,0,-1)     > AvsAvg <- c(1,0,0,0)     > contrasts(my_data$Group) <- cbind(AvsB,AvsC,AvsC,AvsAvg)      And fit the full model (suppose I have already compared the various models with the intercept model):               > model <- lme(Value~Group+Condition+Group:Condition, random=~1|Subject/Condition,                    data=my_data, method="ML")      However, sometimes the following error appears:               Error in MEEM(object, conLin, control$niterEM) :     Singularity in backsolve at level 0, block 1      Other times, the error doesn't appear but the output doesn't show all the contrasts. For example, here `AvsAvg` is disappeared:               Fixed effects: Value ~ Group + Condition + Group:Condition                                    Value Std.Error DF   t-value p-value     (Intercept)               21.793124 1.2435711 44 17.524631  0.0000     GroupAvsB                 -0.255563 2.1539283 44 -0.118650  0.9061     GroupAvsC                  0.176721 2.1539283 44  0.082046  0.9350     GroupAvsD                  0.492208 2.1539283 44  0.228516  0.8203     ConditionCond_2           -0.016164 0.2029039 44 -0.079661  0.9369     GroupAvsB:ConditionCond_2 -0.270565 0.3514399 44 -0.769874  0.4455     GroupAvsC:ConditionCond_2  0.217103 0.3514399 44  0.617754  0.5399     GroupAvsD:ConditionCond_2 -0.277043 0.3514399 44 -0.788309  0.4347      Generally, it appears to me as if a maximum of three contrasts are shown. **Am I mistaking something with the contrasts or with the model?** Thanks in advance.


title: Outlier detection across multiple conditions?

text: \-----EDIT: Sorry, I was somewhat in a hurry when I wrote my question. Let me rephrase it. I intend to conduct a Mixed ANOVA analysis consisting of a Repeated Measures DV and a between-subjects IV (a categorical variable that has three conditions). In order to draw the right conclusion I would first like to determine whether there are outliers in the RM scores by considering the standardized scores for each variable that makes up the RM variable and see whether I find +/-3 standard deviations. However, according to my understanding these standardized values are dependent on the mean, which I expect may differ per condition. Thus, how should I determine my standardized scores? Should I use filters to determine the scores per condition or can I just compute them across all three conditions? I have tried both and there is a difference, some values are edging on 3SD when analyzed across all conditions but they become outliers when I compute the scores per condition (or the other way around). Regarding themy other sources, I don’t have an example available right now but when I looked into this yesterday I mainly found descriptions that stated that the SD scores should be calculated for the DV and any variable above |+/-3| is an outlier. This, of course, I already know and it doesn't help me answer my question. Hope it becomes a bit clearer now. * * * I am sorry if this question has been asked and answered before but I tried searching for it but I could not find a satisfying answer. I plan to run a Mixed ANOVA with a repeated measures and a categorical variable and I would like to check for outliers in the dependent variables (the repeated measurements). However, I am unsure on whether I should check for outliers in each of the RM variables (determine the standardized scores for all of the RM variables and look at outliers within each of them across all conditions), or whether I should consider my data on a per condition basis (set up and use a filter per condition, and then check for outliers in the variables per condition). In general explanations about outliers I only found examples in which the entire variable is considered at the same time and no conditions were discussed so I hope someone can help me out. Thanks a lot.


title: Latent variable model and graphical models

text: At first I was under the impression that latent variable models and graphical models were different things altogether, but after reading some papers about the former class of models, it seems they have more in common than I thought. For instance, in both cases there's a specification of the conditional probabilies. Is there any distinction between the two? How can I view these in relation to each other?


title: Duda, Hart, Stork No Free Lunch Discussion

text: Please see this question regarding Duda, Hart, and Stork's No Free Lunch Thm Discussion Hi all, I was having trouble understanding the description of the NFL theorem in Duda, Hart, and Stork. My issues are the same as brought up in the reference above. I read the answer too, but I'm still having the following problems. Equation (1) in the book, makes sense to me. It says that the expected error rate for a learning algorithm given the training set is: \begin{equation} \mathcal{E}[E|\mathcal{D}] = \sum_{h,F}\sum_{x \notin \mathcal{D}} P(x)[1-\delta(F(x),h(x)]P(h|\mathcal{D})P(F|\mathcal{D}) \end{equation} (Actually, I believe the $P(x)$ should be $P(x|x \notin \mathcal{D})$, but that won't make a difference in the proof.) Now, as stated in the text, $P(h|\mathcal{D})$ is the probability that the algorithm will produce the hypothesis function $h$. So with just a slight change of notation, we could write that the expected error for algorithm $k$ is: \begin{equation} \mathcal{E}_k[E|\mathcal{D}] = \sum_{h,F} \sum_{x \notin \mathcal{D}} P(x)[1-\delta(F(x),h(x)]P_k(h|\mathcal{D})P(F|\mathcal{D}) \end{equation} So the conditional error for algorithm $k$, given target function $F$ is \begin{equation} \mathcal{E}_k[E|F, \mathcal{D}] = \sum_{h}\sum_{x \notin \mathcal{D}} P(x)[1-\delta(F(x),h(x)]P_k(h|\mathcal{D}) \end{equation} Now, the text states "the expected off-training-set classification error when the true function is $F(x)$ and the probability for the kth candidate learning algorithm is $P_k(h(x)|\mathcal{D})$ is given by \begin{equation} \mathcal{E}_k[E|F,n] = \sum_{x \notin \mathcal{D}} P(x)[1-\delta(F(x),h(x)]P_k(h(x)|\mathcal{D}) \end{equation} What I don't understand is this: In the case of a stochastic algorithm, there are many hypothesis functions $h$ that could be produced by algorithm $k$. But in the above equation, only one $h(x)$ appears, and there is no summation over $h \in \mathcal{H}$. The penultimate equation above seems a better way to describe the classification error for algorithm $k$ given $F$. That is why I agree with the author of the above referenced question (especially his #1 and #2). If the equation in the text is correct, then, which $h$ is it? Anybody have any ideas? I must be missing something. Thanks!


title: Determining customer price sensitivies by mining purchase histories

text: Currently, I am looking for methods to calculate some kind of price sensitivity index for customers of an online store based on the purchase histories. I've found a lot of literature in which panel data was used to calculate the price elasticities for a certain product (category) first, and then a regression model was applied to determine the influencing attributes (e.g. high income, ...). My question is, if there any approaches to derive price sensitivities (or a similar measure) on a customer/product (category) level based on historical purchase data.


title: Which test to use to measure if target has been met?

text: The photo below is an example of the result i would like to obtain using the test which i cannot find in the article itself, i was hoping maybe any of you would know. I have tried to carry out chi square test for goodness of fit but what i really want is to able to compare BP between two groups (men vs. women). Would anyone shine a light and tell me what the test name is called? ![enter image description here](http://i.stack.imgur.com/CkB3O.jpg)


title: need some help in data analysis in SPSS

text: I Conducted a study in some patients it included:   1. age   2. gender   3. duration   4. some medical tests   5. pre. and post questionnaire   6. post test results I divided the patients into two groups acoording to duration of complain I need to correlate these two groups to assess wether this difference in durations plays a role in patients test results and performance. How do i corelate them please?


title: Parameters, Estimates

text: I lack some knowledge in the concepts of parameters, estimates and moment (math and stats). I can't find an online easy-to-understand source of information about these concepts. Would you help me with a link or even better with a speech ? Here is one thing that I'd aimed to understand: Source: a book called _Principles in Population Genetics_ > $\sigma_x^2 = \text{E}[x-\text{E}(x)^2]$ > > $\text{Var}(x) = \text{mean}(x^2) - \text{mean}(x)^2$ > > But it is a biased estimator because... > > $\text{E}[\text{Var}(x)] = (n-1)\sigma^2/n$ > > A non-biased estimator for $\sigma^2$ is therefore > > $[n/(n-1)]\text{Var}(x)$ **What is the meaning of all this ?**


title: Plotting a Probit Time Series in R

text: I've searched around a bit, but am having trouble finding an answer to my specific question that isn't too technical. I want to predict the probability that we are at a turning point – arbitrarily defined – for a given time series. I have historical data from which I have determined some times that are "turning points" (1 for yes, 0 for no) and I have a list of regressors for which I actually have data. From this, I would like to a) use a probit or logit model to estimate coefficients on the regressors I have and b) use future (or out-of-sample) values for these regressors to forecast/predict the chance that we are at a turning point now. What I'm having trouble with is putting time series and probit analysis together. I call:               glm(formula = turning ~ . - Date, family = binomial(link = "probit"), data = d1)      I would like to use the `predict` functionality to get a time series of probabilities. Two final thoughts   1. I haven't done much statistics before, so I want to avoid certain pitfalls. My definition of a "turning point" is rather specific in that for over 123 monthly observations (historical) only 4 meet the criterion. Does this hinder the ability for the probit model to work effectively, should I "cushion" each of these points in a range (i.e. every time plus or minus 2 of my current turning point would be switched to 1).   2. Are there better methods for binary forecasting? I notice that in baseball batting average forecasting a beta distribution is used to model the priors and updated. Is there a logical way to bring that into this?


title: Cross-Validation: Conditional test error vs expected test error

text: My textbook on cross-validation is Element of Statistical Learning (Hastie et. al., 2nd ed.) In sections 7.10.1 and 7.12, they talk about the difference between conditional test error $E_{(X^*,Y^*)}[L(Y, \hat{f}(X))|\tau]$ and expected test error $E_\tau [E_{(X^*,Y^*)}[L(Y, \hat{f}(X))|\tau]]$. Here $\tau$ is the training data set, $L$ is the loss function, $\hat{f}$ is the model. $E$ is the expectation. They explained that CV only estimates expected test error well. My question is, is there any reason why we would care about the conditional test error? The only reason I could think of is that we want to answer the question 'If God puts n data set on the table, but only let us take 1 home to fit our model, which one should we choose?'


title: How to use 1.7159 * tanh(2/3 * x) as activation function?

text: I have a simple neural network and it works with the logistic function as activation function. Now I want to avoid the saturation problem by substituting the logistic function by the hyperbolic tangent:               #define SIGMOID(x) (1.7159*tanh(0.66666667*x))      #define DSIGMOID(S) (0.666666667/1.7159*(1.7159-(S))*(1.7159+(S)))      But the network never converges, the MSE stays the same throughout the training. Here's my training samples:               double training_data[][4]={                 {0, 0,  0,  -1},                 {0, 0,  1,  1},                 {0, 1,  0,  1},                 {0, 1,  1,  -1},                 {1, 0,  0,  1},                 {1, 0,  1,  -1},                 {1, 1,  0,  -1},                 {1, 1,  1,  1}};      The network does converge if I use the original (non-scaled) hyperbolic tangent function, that is:               #define SIGMOID(x) (tanh(x))     #define DSIGMOID(S) (1-((S)*(S)))      Do I miss something? E.g. Scaling the output to match the range (-1.7159, 1.7159) or anything?


title: Choosing an appropriate part of an unreliable dataset

text: I have a dataset with $\approx2000$ entries (for example, model of car). For each car model I know the weight of the car, and the power output. I don't know the price or age, which is likely to affect the power to weight ratio of the car. I have a reasonable amount of faith in the data, but some manufacturers might fib a bit (maybe 5% of the data is unreliable). **How can I determine, in a defensible way, the best power output for a given weight of car?** I'm mathematically literate to an engineering undergraduate level, but don't have any statistics experience beyond high school. I've thought of the following:   * Sort cars by power/weight ratio and then take the top $n\%$. Problem is the determination of $n$.   * Graph the data and plot a trendline through it. This will underestimate the best power output by a large factor.   * Take the best power to weight ratio and just use that. Problem is that power to weight ratios do not necessarily scale linearly with car weight (it may be easier to produce a small car with a high power to weight ratio for example).   * Bin the data into weight categories, then do the same as above. Problem is that taking a single datapoint from an unreliable dataset doesn't seem right.   * As above, but trying to cross-check the data against outside information (for example, magazine performance tests). This information isn't available for all cars, so probably unfeasible. For my curiosity, what is the name of these class of problems? It seems like something that would occur commonly, but I don't know how to search for it. Is there an English to Statistics translation tool somewhere?


title: Looking for step by step example of sampling from DAG in Bayesian model

text: I am looking for a tutorial type example that shows the step by step process sampling from a simple hierarchical model. For example, I am trying to study the distribution of `p` in a Bernoulli experiment where I have a set of 10 data / observations (`h[i]`).               model {       p ~ dunif( 0, 1 )       for( i in 1 : 10)  {         h[i] ~ dbern( p )       }     }      I am not clear as to how the WinBUGS (or similar samplers) would generate the correct samples for `p` values given my `h1-10`. Is there a paper or article that explain this?


title: Random forest with a mix of categorical and lexical features

text: I'm working on a binary classification problem for which the dataset is mostly composed of categorical features, but also a few lexical ones (i.e. article titles and abstracts). I'm experimenting with Random Forests, and my current strategy is to build the training set by appending the $k$ best lexical features (chosen with univariate feature selection, and weighted with $tf- idf$) to the full set of categorical features. This works reasonably well, but as I cannot find explicit references to such a strategy of using hybrid features for RF, I have doubts about my approach: does it make sense? Am I "diluting" the power of the RF by doing so, and should I rather try to combine two classifiers specializing on both types of features?


title: Analyzing presence/absence data for art history

text: How can I best analyze the presence/absence data that I have gathered about the illuminations in a medieval manuscript? I have identified individual motifs in chapter headings and created a spreadsheet that has 1 or 0 for each (there are 37 motifs in total). There were at least 4 people who worked on the manuscript, which is clear from simple motif that appears on all of the pages analyzed. I did PCA and color- coded the data points with this control variable, and it appears that the presence/absence of these other motifs does, to some extent, correlate with it. I was thinking that co-occurence analysis would be the next direction to explore, although I would be eager to learn of any other statistical techniques that might be useful for this type of data. Are there any pitfalls that I would need to be aware of, or any R packages that people would recommend? Specifically, I want to figure out which motifs correlate with each of the 4 artists. If possible, I'd also like to be able to take into consideration the possibility that similarities could also result form shared models that the artists were looking at. More specifically, I have run a correlation analysis using pearson's method on the pages created by a one of the four artists. How can I find closest correlations from the resulting dataframe, and how can I determine if they are statistically significant?


title: Proof of asymptotic variance

text: How do you prove that $X_n - E[X_n] = O_p(\sqrt{Var(X_n)})$ It's used in my textbook and I don't know where they get it from.


title: When using a Neural Network to classify more than two classes, is it better to have multiple output nodes (one for each class) or one output node?

text: Currently, I am using a neural network to classify data into one of three groups (a logistic activation function is used on all but the output nodes). I can train the neural network in two ways: 1) For each observation $X_{i}$ I can have one output variable $O_{i}$ that can take three values: $1,2,$ or $3$. Or, 2) For each observation $X_{i}$ I can have three output variables $O_{i1},O_{i2}$, and $O_{i3}$ that would take on the values $0$ or $1$ based on what the class label for that observation is. So for example, if observation $i$ belongs to class 2, $O_{i2}$ would be $1$, whereas $O_{i1}$ and $O_{i3}$ would be $0$. Are the above two approaches equivalent? Is one better than the other? Are there any theoretical reasons to believe one should produce better results than the other?


title: NN: Should we apply weight decay to the bias?

text: In CS294A lecture notes, Andrew Ng writes (about autoencoders): "Usually weight decay is not applied to the bias terms... Applying weight decay to the bias units usually makes only a small different to the final network, however". Is there any particular reason for which we shouldn't apply weight decay to the bias terms? Does it reduce the performance of the network?


title: ks.test and ks.boot - exact p-values and ties

text: I am confused by the behaviour of ks.test (package stat) a) in the presence of ties and b) if one-sided while doing a two-sample test. Documentation: "Exact p-values are not available for the two-sample case if one-sided or in the presence of ties." I ask if black (experiment) and red (control) follow the same distribution function without knowing the underlying distribution function. In my hands exact p-values are computed if one-sided and in the presence of ties (according to the warning message). But two-sided the p-value is just < 2.2e-16 but not an "exactly" reported. If interested you may download the data as .Rda (length of the vector ~ 9000): https://www.dropbox.com/s/xl29jvpurkbwqpm/black.Rda?dl=0 https://www.dropbox.com/s/5biptm1xet36v3v/red.Rda?dl=0 Example:               ks.test (black, red)              Two-sample Kolmogorov-Smirnov test          data:  black and red     D = 0.0731, p-value < 2.2e-16     alternative hypothesis: two-sided          ks.test (black, red)$p.value           [1] 0     Warnmeldung: # means warning message     In ks.test(black, red) :       im Falle von Bindungen sind die p-Werte approximativ # "Bindungen" means ties          ks.test (black, red, alternative="g")$p.value # not as expected          [1] 1.235537e-23     Warnmeldung:     In ks.test(black, red, alternative = "g") :       im Falle von Bindungen sind die p-Werte approximativ          ks.test (black, red, alternative="l")$p.value          [1] 0.0005651143     Warnmeldung:     In ks.test(black, red, alternative = "l") :       im Falle von Bindungen sind die p-Werte approximativ      I tried `ks.boot` (package "Matching") that claims to work for two.sample tests with ties and "provides correct coverage even when the distributions being compared are not entirely continuous." Same story. I get exact p-values for one-sided conditions only. For instanche:               ks.boot (black, red, alternative="l")          $ks.boot.pvalue     [1] 0.001          $ks          Two-sample Kolmogorov-Smirnov test          data:  Tr and Co     D^- = 0.0275, p-value = 0.0005651     alternative hypothesis: the CDF of x lies below that of y          $nboots     [1] 1000          attr(,"class")     [1] "ks.boot"      Did I missunderstand the sentence "exact p-values are not available for the two-sample case if one-sided or in the presence of ties?" I thought the sense was: No exact p-value if one-sided or ... Are the p-values of ks.test (two.sample, one sided) "correct"? In terms of delivering exact p-values ks.boot was not superior. Can anybody please comment on this? Thanks Hermann


title: Whether to include time as a covariate in repeated measures analysis

text: Suppose we want to model the temperature ($1$ = hot, $0$ = cold) of 5 regions over the following time periods: 1 month, 6, months, and 4 years. Suppose a predictor variable is $X$ (amount of deforestation, 1 = a lot, 0 = very little). In modeling the main effect (deforestation), should we include time as a predictor?


title: What is the best training method for 15-30k records with 5-12 features to classify data into 2 groups?

text: I don't know if this is a suitable question here. I have data from 2 groups, neutral and positive. Here how the data looks like               A   02.840  00.960  00.950  00.995  00.995  00.602  1     B   04.730  01.000  01.000  00.999  01.000  01.000  1     C   NULL    00.970  00.005  00.078  00.987  00.000  0     D   NULL    01.000  01.000  00.010  00.675  00.451  1     E   00.530  01.000  00.022  00.231  00.443  00.000  0      The first column is the record name/ID and the last column is its class/group. The six columns in the middle are the features which some of them might be NULL. With the knowledge that I studied years ago, I'm gonna try ANN and SVM first. But I'm not sure if I am on the right track. So I wonder if anyone here can guide me.


title: When two events $A$ and $B$ have no result in common

text: I and my friends just had a little discussion whether the events are independent or dependent if they have no outcome in common. I thought that they have to be independent. When two events are independent, then $P(A)=P(A\mid B)$. Is the information given in the question enough to establish this? If you look at it as a Venn diagram, then if there is no overlap between A and B, then they are independent. But my friend objected and said that this depends on the sample space. So let's say we have two events: $P(A)=6/12$, $P(A\mid B)=2/4$, and $P(B)=4/12$, then obviously $P(A\mid B)$ equals $P(A)$. But for me, something smells fishy here. I know this is more a stochastic question that a statistic question, but maybe someone can help.


title: Optimal parameter selection by repeated  k-fold

text: I am working on Lasso problem and the selection of the optimal tuning parameter with $k$-fold procedure, say $k=10$. Since this procedure relies on random subsampling, value of the optimal parameter will change each time I repeat the procedure. As an example, it can be 0.32, then 0.41, then 0.29, etc. Two questions:   1. Can I use repeated $k$-fold and average the results?   2. How do I compute the standard error in order to use one standard rule?


title: Problem with Developing a Model for an Exploratory Study (PCA EFA)

text: * ** _Background:_ *** I am doing an academic study that is more exploratory in nature (my advisor did not want me to develop any hypothesis). So I read books and conducted interviews to try and gather a large list of items that people feel are important and influential to the topic being studied. From that list, I developed them into questions to ask in a survey which I got a good number of responses (n>500). The 50 questions survey include a list of 35 items that people think are important for the topic and 15 demographic questions about the respondent. * ** _My Problem:_ *** Now I have a problem on choosing the correct analyses. What I want to achieve is through the analyses of the data, come up with a model that consists of all the items that people find as important for the topic being studied. But at the same time, my advisor doesn't want me to reduce the items from the survey too much because he feels all the items in the survey are important, but rather, based on the items respondents indicated as important (have high mean score?), group them together to form factors in a model (without reducing too much data) and look for any relationships. I am not sure which type of analyses would be suitable for this.  So far I think I chose the wrong type of analysis so I appreciate some help. After reading up on potential analyses that are suitable for studies that do not have hypotheses, I conducted a Principal Confirmatory Analysis on the non- demographics questions and formed them into 4 factors. But in doing so, there were about 10 non-demographic items that were deleted. My advisor is now asking me if PCA is the best analysis to choose because he saw that most of the mean scores for all the questions were high (at around 3-4 out of 5) so he assumed those are all important and should all be included into a model. I believe what he wants is all the questions with a high mean score which indicate that the respondent felt was important should really be kept in the model and should not be deleted from the model. _Ultimately, I am looking for suggestions as to what analyses would be suitable for the outcome I am trying to achieve...to develop a model that consists of all the items that people find as important from the survey (and if the mean score is any indication, a very high number of the items are score 3-4 out of 5 in the likert)_ I would appreciate any suggestions.


title: First quick glance at a dataset

text: Please pardon my ignorance, but... I keep finding myself in a situation, where I'm faced with a bunch of new data I managed to find. This data usually looks something like this:               Date     Number1  Number2  Category1  Category2     20120125      11      101        Dog      Brown     20120126      21       90        Cat      Black     20120126      31      134        Cat      Brown     (...)      Usually at first glance I can't really tell if there are any trends here. The correlations between the various columns may not be very significant, but I would be delighted if I didn't have to manually create a plot for every possible combination of columns/categories. Is there a tool out there which would accept a table of data along with information which columns should be treated as numbers, dates and categories and then proceed to plot:   * correlations between each two numerical columns   * correlations between each two numerical columns, with separate trend lines for each category   * each number column as a time series,   * each number column as a time series, separated by category,   * etc. In the end this would generate a large number of plots, most of which would show only noise. Ideally, the tool could score the plots by correlation and in the end display a slideshow starting with the highest scoring plots. This would be a very imperfect, but useful first glance at the dataset. So? Is there a tool everyone uses for this and I just don't know about it, or is this something we need to make?


title: Updating the sample variance, sample by sample

text: > **Possible Duplicate:**   >  Computation of new standard deviation using old standard deviation after > change in dataset I know the sample variance of n samples and its mean and n itself. Now I would like to add one more sample so I have n + 1. Is there a way how to add the sample without recalculation the whole term. (computational efficient). Approximation would be fine too.


title: What do they analyze

text: Please, this is maybe a meta-statistics question but let me ask this. What would you tell a mathematician to answer the question: "Mention few examples what kind of data analysis the insurance companies or banks do?" Something a la "They have the set of claims and compute the mean."... (obviously a stupid example).


title: Identifying modes in floating point data

text: My new computer clock runs at a rate that changes in odd steps over time, even after I tuned it via the Linux adjtimex software. Here is a plot of the change in the cumulative clock drift for each of about 1700 samples taken every 10000 seconds, with a bunch of missed points and outliers when I was off the network and ntpdate wouldn't work. E.g. early on, the clock gained about 0.09 seconds every 10000 seconds (9 ppm). ![clock-rate.png](http://bcn.boulder.co.us/~neal/media/images/clock-rate.png) I'm looking for some clever library functions that can automatically identify what I'll call the various statistical "modes" of this data set - i.e. there is one mode in the middle where y is about -0.02 for a long time, then another one near 0.09 early on and one at 0.202 at the end. Most code for finding the mode of data that I've seen deals with integers and discrete data, but as you see this has pretty messy floating point values. At any rate, ideally I'd like a summary that automatically finds the modes I identified above, and also gives me a standard deviation for each. Start/stop points for each mode for extra credit. Python code preferred.


title: Test identicality of discrete distribution

text: Is there a standard way to test whether two vectors were drawn with the same **discrete** distribution in R? Something like a Kolmogorov-Smirnov test, but for discrete distributions. I think two-sample chi-squared test would be appropriate. Does any package provide it? I can't get `chisq.test` to work for me.


title: When should one _avoid_ Bootstrapping?

text: > **Possible Duplicate:**   >  What are examples where a "naive bootstrap" fails? Perhaps this question will turn out to be a bit soft, but I think it may have a hard answer. Feel free to move it to discussion or CW if it turns more philosophical. I am a computer scientist by training, and perhaps because of the old proverb about how things appear when you're holding a hammer, bootstrapping seems like the most natural solution to just about every problem I see in hypothesis testing. For example, today, while answering a question here in CV, I even suggested bootstrapping for what was really just a straightforward t-test situation (without noticing what it reduced to!). This lead me to a couple of questions: First, is there a reason to avoid boostrapping other than the tendency to inadvertently conceal assumptions about the underlying data or distributions? That is, **are there reasons to avoid bootstrapping other than "It's easier to shoot yourself in the foot."?** Second, perhaps a little more philosophically, in many fields, my own included, one sees a lot of bad statistics. I'm curious about whether, for a person operating with little statistical education, it is easier to break things when using conventional frequentest methods of hypothesis testing or using bootstrapping. Obviously the conventional methods are harder to break in the sense that assumptions are clearly stated, but in my experience, people often do not read the fine print, and ignore the assumptions anyway. **In the case that a person blindly applies each method to the data, do you think they are more likely to draw bad conclusions using bootstrapping than they would when applying conventional methods of inference without heeding the underlying assumptions?**


title: Distribution of extremal values

text: If an item follows normal distribution, average also follows normal distribution. What about minimum and maximum?


title: How to assess differential risk of disease across three groups after adjusting for other risk factors?

text: I have 3 different groups: A, B, and C:   * A (has a medical condition) has 30 entries    * B (another medical condition) has 31 entries    * C (control group) has 55 entries All participants have information on the same set of variables. I want to assess whether there is a statistically significant difference in risk of a specific disease between:   * Group A and C    * Group B and C    * Group B and A Also, all groups need to be adjusted for other known risk factors.


title: Drop highly correlated items?

text: I have two IVs that are highly correlated with each other at 0.979 (Pearson) & 0.919 (Kendall's). IV1: Quality of response   IV2: Quality of Technical Advice Sample Size: 252 Considering the similarity in the phrasing of the IVs, could it mean that the respondents do not perceive a difference between the two measurements? Would it be appropriate to drop one of the items?


title: Coerce contingency table to report unrepresented classes

text: I have a contingency table computed from the out of sample observations. In some cases not all classes are represented in the out of sample and hence do not appear. I want to coerce the contingency table to report zeros for classes not in the out of sample. How can I do it efficiently. For example my contingency table is                 x z     x 1 2     y 1 3     z 2 2      but i want it to look like                 x y z     x 1 0 2     y 1 0 3     z 2 0 2


title: Neural Network: Spikes using conjugate gradient

text: I'm using OpenANN to train a neural network with one hidden layer and a softmax output layer with cross entropy as the error function. For my application, the conjugate gradient algorithm seems to be a good performing algorithm compared to for example LBFGS. The image below shows the training error and validation error for one configuration of the neural network (20 inputs, 4 hidden neurons, 12 outputs). Is it 'normal' for CG to have spikes in the performance during training, and if so, what is the explanation for the spikes? The spikes only appear for CG and not for LBFGS and LMA. ![enter image description here](http://i.stack.imgur.com/EOY6Z.png)


title: How do I detect the number of distributions in a set of data?

text: I'm trying to detect subpopulations of data in a set of data. The problem is that the data may have more than one subpopulation and the subpopulations may not be normally distributed. For example, the data could look like this.... ![enter image description here](http://i.stack.imgur.com/98yw9.png) The Y axis is the count and the X axis is the value. What I want is to determine the number of subpopulations and their peak values. The programming language I'm most familiar with is Python, however, my methods have been less than successful. It becomes problematic because in the example I showed there is a subpopulation with two peaks. My method looks for peaks and valleys so in this example it would identify three subpopulations even though there are actually two. I'm new to statistics so I was wondering if anyone here could suggest a way to detect subpopulations of data.


title: How to test my hypothesis?

text: With a paired sample t-test I found out that participants communicated more pride about a team achievement than about an individual achievement. What I would like to examine now, is if this difference between communicated pride about a team achievement and an individual achievement is larger for Indians compared to Americans. I have a sample consisting of 148 Indians and 268 Americans. All participants had to indicate their communicated pride over both a team achievement and an individual achievement (so it is a within subjects design). Which test should I run? I thought about computing a new variable which is computed as communicated pride over a team achievement - comunnicated pride over an individual achievement. Then I have to determine if these differ from each other, but I don't know how!?


title: Non-constant standard deviation in residuals

text: I am fitting a model in the frequency domain, and my fit looks as follows: ![Model fit](http://s3.postimg.org/8aop8zepv/fit.jpg) As you can see, the model function does not fit the data perfectly, especially in the higher frequencies. So, I examined the residuals and found that a polynomial of the first degree is the missing term in the model function. What concerns me more is that variance of the data decreases with frequency. What implications does that have on my analysis? (It seems that looking at the autocorrelation of the residuals also doesn't prove very useful.) Any ideas would be very helpful.


title: Multiple Discriminant Analysis, Linear Discriminant Analysis, and Multidimensional scaling - how are they related?

text: Some time ago when I took a Pattern Classification class, the "concept" was introduced as Multiple Discriminant Analysis: You want to project your data onto a subspace (if you are interested in dimensionality reduction) that yields a better separation between different classes (e.g., in a supervised class. task with multiple labeled classes). I think I understood the technique and it is kind of related to a principal component analysis - with the difference that you have 2 scatter matrices: the scatter matrix within and the scatter matrix for in between classes ... Anyway, later on when I was trying to dig into so machine related topics, I stumbled across the terms "Linear Discriminant Analysis (LDA)" and "Multidimensional Scaling (MS)" and I am wondering how those are connected/related. If I understood correctly, you can basically start with 2 categories   \- Linear Discriminant Analysis   \- Quadratic Discriminant Analysis which are describing/categorizing the "complexity" of the problem. And if I speak of "Multiple Discriminant Analysis (MDA_" it basically revers to how many dimensions each class has!? Basically, that MDA is a subcategory of LDA and QDA!? And Multidimensional Scaling would be more like a super-category, where it just refers to the process of scaling a multidimensional feature space, e.g., via MDA or PCA etc.!?


title: Measures of goodness-of-fit using multiply imputed data in Zelig

text: I am running a logistic regression model in R using multiply imputed data created using Amelia II, which I am then analyzing using Zelig. I would like to be able to report some measures of goodness-of-fit (e.g. likelihood ratio, pseudo R-squared, Hosmer-Lemeshow), however none are provided in the default Zelig output and I haven't been able to figure out a way to extract any from the `zelig()` object. Do measures of goodness-of-fit need to be calculated differently when using multiply imputed datasets? Are there any R packages that are able to do this? I have looked into several packages that provide measures of goodness-of-fit, such as pscl, however they only work on glm objects, not MI objects created when using Amelia and Zelig. Thanks in advance for your help!


title: Estimating ARMA equation using lm() in R

text: Is there a way to estimate an ARMA equation using the `lm()` function in R without using `arima()`?


title: Feature Normalization/Standardization before or after Feature Selection?

text: Should the process of feature normalization/standardization be done before or after the feature selection process?


title: Removing Time-Series Variance from Panel Data

text: We are working with panel data. But we want to study only the cross-section part of the panel data. So can anybody please tell me how to do any kind of data transformation, so that I can remove the time-series variance from the panel data? Thank you.


title: Identifying What Causes a Variable to Increase

text: Say I have a dataset with several continuous and categorical variables, and I want to identify what variables (values or properties of these variables) may cause one of the continuous variables to increase. How can I model this problem? Several approaches come to mind (e.g. **logistic regression** , **ANOVA** , **decision trees** ) but most (all?) of them require that I define a **condition** on my target variable and thus make it a **categorical** variable and/or treat the problem as a **classification problem**. But what if my **dependent** variable is a continuous variable and I am looking for properties of the independent variables (continuous and categorical) that make my dependent variable to increase? Is **linear regression** the best I can do? What if my goal is to find **ranges** of the predictors that correlate statistically with increases of my dependent variable? What **models** do I have available for this problem, and most importantly, what **questions** should I ask myself about the problem to identify the proper one?


title: Back Transforming Rates in Poisson GLM with Box and Cox Transformation

text: Suppose I have fitted a Poisson GLM to model rates as follows:               > fit.1=glm(response~X1+X2+ offset(log(population)),family=poisson,data=...)      I can get the estimated rates by using two new values for X1 and X2 and population=1 as follows:               >new.data=data.frame(X1=new.X1,X2=new.X2,population=1)     >estimated.rates=predict(fit.1,newdata=new.data,type="response")      Following Subsection 13.4.5 in Introduction to Linear Regression Analysis, 5th Edition, I then double check the Deviance Residuals using `qqnorm` and found that they are not even approximately normally distributed. So to revise my model, I used a box and cox transformation (lambda) to improve the model and fitted a new quasi-Poisson model:               > fit.2=glm(I(response^lambda)~X1+X2+ offset(log(population)),family=quasipoisson,data=...)      I can again estimate the rates similarly as before:               > estimated.rates.2=predict(fit.2,newdata=new.data,type="response")          The `estimated.rates.2` is on the scale of the transformed response (using the box and cox transformation). Is it possible to back transform the rates `estimated.rates.2` into the original scale (i.e. similar to `estimated.rates`)?


title: AR models on non stationary data

text: i am currently reading Diebold and Li's 2006 paper: Forecasting the term structure of government yields where the authors fit, albeit simple, AR(1) models on clearly non stationary data. Why is this allowed in certain situations? and what are the consequences/implications of doing so? Thanks


title: computing EP messages in factor graph

text: I'm using infer.net for implementing a dynamic Bayesian network model. every node in each layer is dependent on all nodes in previous layer (we want to train wji, the weights of edges between each node j from the previous level (time t-1) to node i in current level (time t)). node values are shown by yi (value of node i). The dependency is a complicated function like: yi(t) = y(t-1).(1-alpha.deltaT) + beta . 1 / (1 + exp(-WeightedSigma)) WeightedSigma = Sigma (Wij * yj(t)) Now I want to put a factor node in my factor graph implementation in infer.NET. I don't know how should I implement such factor node and its message passing. I'm so confused. Can anyone please help me in this case? I appreciate any help. Thanks


title: what are the alternative open source tools for PredictionIO?

text: PredictionIO is good to be enough for content discovery and recommendation but it seems it does not support classification. Then I should use a different tool then Prediction IO for my prediction server. Do you know good alternatives? I know Google Prediction Server, Azure ML, BigML and looking for more alternatives to compare them. Best :)


title: Jaynes' Derivation of Herschel-Maxwell for Normal Distribution

text: I am reading the following paper: http://www-biba.inrialpes.fr/Jaynes/cc07s.pdf and cannot seem to figure out how Jaynes is deriving (P2) and below (specifically the log arithmetic log[f(x)/f(0)]....) Can someone help me out?


title: Question about probabilistic principal component analysis

text: I have a rather basic question about Probabilistic Principal Component Analysis, which I am now trying to apply to a real-world problem. In PPCA, the crucial assumption is that the generating process of the observations in $R^n$ is $t=Wx +\sigma^2\epsilon$, where x are iid standard gaussians in $R^q$ (with $q \le n$) and $\epsilon$ are iid standard gaussian vectors in $R^n$. The authors find that the MLE solution is $\hat W=U_q (\Lambda_q - \sigma^2 I)^{1/2}$, where $\Lambda_q$ is the diagonal matrix of the $q$ largest eigenvalues of the empirical covariance matrix, and $U_q$ is the corresponding eigenvectors submatrix, and $\sigma^2$ is the average of the remaining smaller eigenvalues (see sec. 3.2 of the linked paper). My question is simple. The covariance matrix is $C=WW'+\sigma^2 I$. Directly plugging in the above result, we obtain $\hat C=U_q\Lambda_q U_q$ (the $\sigma^2$ terms cancel out). Can this be correct? The covariance matrix would be rank-deficient. Am I missing something?


title: Why is using cross-sectional data to infer / predict longitudinal changes a Bad Thing?

text: I'm looking for a paper which I hope exists, but don't know if it does. It could be a set of case studies, and / or an argument from probability theory, about why using cross-sectional data to infer / predict longitudinal changes may be a Bad Thing (i.e. is not necessarily so, but can be). I've seen the mistake made in a couple of big ways: inferences were made that because richer people in Britain travel more, then as society gets richer, the population as a whole will travel more. That inference turned out to be untrue for an extended period - more than a decade. And a similar pattern with domestic electricity use: cross-sectional data implies big increases with income, that don't manifest over time. There are several things going on, including cohort effects and supply-side constraints. It would be very useful to have a single reference that compiled case studies like that; and / or used probability theory to illustrate why it is that using cross-sectional data to infer / predict longitudinal changes can be very very misleading. Does such a paper exist, and if so, what is it?


title: Difference in frequencies in two groups of unequal size

text: I want to compare the frequency of retweets between 2 groups of unequal size. I've got the total tweet count and the retweet count for each user. I was hoping the R Cookbook webpage would help but the tests presented there seem to be more suited for repeated measures and equal sample sizes for the two groups: http://www.cookbook-r.com/Statistical_analysis/Frequency_tests/ I will be grateful for advice on what test to use here. My dataset looks a bit like this:               userid tweetCount ReTweetcount Group     1      45        3             A     2      100       25            A     3      23        0             B


title: Dependency between two independent experiments pointless?

text: Let's say I have two experiments of the same type, where I measure something x times in one location and y times in another location, is there any point in looking for a dependency in values between these two experiments? In case it matters, there's a difference in sample size and both samples have a near normal distribution. For example, if I would want to measure the length of people in Asia and the length of people in Europe, I could investigate a null hypothesis test H0 'The average length of Asians is smaller than the average length of Europeans.'. This is a valid study. However, it seems awkward to wonder if the length of Europeans could be influenced by the length of Asians. After all, the length can be influenced by other potential factors, like food, weather, genes, ... If the world population were suddenly to take in growth hormones, then the average length of both Asians and Europeans are likely to increase. It seems silly to assume Europeans increased in length because Asians did. Europeans did increase in size along with Asians, but not because of a direct result. Is it only worth considering whether two variables have any sort of dependency if they have paired values (like age and weight)? If not, what would be the logical next step to take in order to test for dependency? I've also noticed that I cannot perform any kind of correlation test since I do not have any paired values. Can I assume as an immediate result that my two variables have no linear dependency? I tried to find the answer myself, but there are quite a lot of tests and things to be aware of; it gets a little overwhelming for a novice.


title: How to quantify sampling intensity for Repeated Measures ANOVA

text: I have a dataset that I have been running one-way repeated measures ANOVA on. The people I am working with are interested in answering how much of an actual change in a measurement would have to happen in order for that change to be considered statistically significant. So the way to measure the smallest detectable difference that I have been using is from Zar's Biostatistical Analysis, 5th edition 2010: $$ \delta = \sqrt{\frac{2ks^{2}\phi^{2}}{n}} $$ (Equation (10.34) in Zar's book). So if we can specify a significance level, power level, sample size, and have a estimate for $\sigma^{2}$ we can solve for $\delta$. ($k= p-1$, $n=$sample size, $s^{2}$= estimate of $\sigma^{2}$, and $\phi^{2}$ is the critical value for a specific significance and power level). The nice thing about the minimum detectable difference (MDD) is that it is in the units of what we are studying. So I have two questions. One: Is the MDD the only option to determine how much of a change had to have occurred to detect a difference? I have only ran into the MDD in Zar's book and with the people I am working with. The main question. Sampling has not been consistent for this data set. Here is a plot of a subset of the data set. (Sorry for using a URL link for the picture. It is huge and I couldn't figure out how to scale it down without re- creating it for this question.) What is a good way to measure the effect of sampling intensity? For example, if we say our total sample size remains the same but had a consistent sample of every 7 years instead of having one that fluctuates like it does in reality, what kind of change in the MDD would that create? Only way I have thought of testing this is to systematically drop one year's samples and see how that influences the MDD. Although that would not be keeping a consistent sample size... This question is a lot more complicated then I had hoped! I am looking for any advice or input into these two questions. Any help would be greatly appreciated. If there's anything else I am missing or that is confusing, please let me know. I usually do the ANOVA in SAS and calculate the MDD using R.


title: MLE estimator - dividing the log likelihood by n gives different result

text: The log likelihood is as follows: $nln\beta + (-\beta-1)\sum ln(x_i)$ Dividing the log likelihood by n gives, $ln\beta + \frac{(-\beta-1)}{n}\sum ln(x_i)$ Using these two log likelihoods, I got the same MLE estimator, which is $\frac{n}{\sum ln (x_i)}$. However, I was confused by the different asymptotic properties. For the first log-likelihood, the information matrix is $\beta^2/n$, for the second, it is simply $\beta^2$. I understand that simply because of algebra, the results will be different. But I want to understand the statistical property: why do these two MLEs for two log-likelihood have different asymptotic variance? How is the difference related to the fact that the log-likelihood was divided by $n$?


title: Fitting different models to my data in R

text: I have number of seeds fallen at different distances. I need to know which model fits these data, if any model does. There are different models that have been used to fit these data: negative exponential, Weibull, lognormal. I need to know how to fit all of them in R, starting from 0, generating the data. I only managed to fit the model 2-parameter asymptotic exponential, but I didn't know which parameter to look at, or generate, to know if it was a good fit. The data directly from my Excel is this:               Distancia (distance)      5.4 9.69 11.32 12.84 12.98 14.92 17.76 17.94 17.96 22.17          61.37 63.18 87.88 94.24 112.02 116.91 45.23 141.2          Abundancia (number)     1 2 9 49 24 0 18 3 5 3 10 0 3 0 7 9 6 3


title: Whitening a mixture of Gaussians

text: Suppose I have a dataset that was sampled from a mixture of Gaussians: $$ X \sim \sum_i w_i \mathcal{N}(\mu_i, \Sigma_i)$$ Technically, I can center and then whiten $X$ so that it has zero mean and unit covariance. But something about that process seems unsatisfactorily oversimplified: for example, what if the means of the mixture were arranged along a semicircle, such that no values in $X$ actually occurred at the overall sample mean ? It seems like this is a case where whitening wouldn't really have much meaning in a global sense. However, given only a dataset, one could theoretically fit a mixture of Gaussians model to it, and then use the estimated covariances to whiten points in local regions of the data space. My question is, has this process already been formalized somehow under some name that I have yet to come across? Does it even make sense? What are some good references I ought to consult to learn more about "local whitening"?


title: Are pairwise-convex cost functions also globally-convex?

text: I would like to get a better intuition for a cost function of an optimisation problem. The cost function depends on a number of real-valued arguments. I would like to find out if the cost function is convex. Obviously I can visualise the cost function for 2 arguments in a 3-dimensional plot. I did this for all pairs of arguments and it turned out that the cost function is convex for all these pairs. Now I would like to know: does this imply that the cost function is also globally convex?


title: Constraint optimisation

text: Consider the constraint optimization $\text{argmin}_{\beta}(f(\beta)+\lambda g(\beta))$ can someone define $\beta(\lambda)$. That is, what is the relationship between $\lambda$ and $\beta$?.


title: Why are my p-values different from Stata's?

text: Simple linear regression in MATLAB. Something is wrong with my p-values. Stata reports that the p-values should be 0.000 and 0.207, whereas my code reports that they should be 0.0001 and 0.2988. All other statistics are the same as Stata's output. Please help me see the error of my ways.               %number of observations     n = 10;          x =    [1.0000, 0.4617; ...             1.0000, 0.7532; ...             1.0000, 0.2145; ...             1.0000, 0.5573; ...             1.0000, 0.5586; ...             1.0000, 0.4195; ...             1.0000, 0.4890; ...             1.0000, 0.9292; ...             1.0000, 0.2540; ...             1.0000, 0.7025];                    y =    [108.3156; ...             103.0389; ...             170.0043; ...             167.9905; ...             185.0679; ...             190.1774; ...             135.8128; ...             125.5962; ...             146.6757; ...             143.1218];                 %calculate coefficients     beta = (x'*x)\x'*y          %calculate residuals     u = y - x*beta;          %calcuatate sum of squares residuals     s_2 = (n-2)\u'*u;          %calculate t-statistics     design = s_2*inv(x'*x);          %calculate standard errors     stn_err = [sqrt(design(1,1));sqrt(design(2,2))]          %calculate t-statistics     t_stat(1,1) = sqrt(design(1,1))\(beta(1,1) - 0);     t_stat(2,1) = sqrt(design(2,2))\(beta(2,1) - 0);          t_stat          %calculate p-statistics     p_stat(1,1) = tpdf(t_stat(1,1),n-2) + tpdf(-t_stat(1,1),n-2);     p_stat(2,1) = tpdf(t_stat(2,1),n-2) + tpdf(-t_stat(2,1),n-2);          p_stat


title: Method to find pattern

text: Assume I tracked the behavior of users when they visit a Website. In detail I track the movement from which site to site they click within the Website. So I have a lot of tuples with (userID, site, time). Now I want to visualize if there are any patterns or clusters for the movement. Let's say most of the users click step by step through the Website, and another group visits site 1 then site 2, then they go back to 1 and then to 2 and then to site 3. What method can I use to classify the behavior?               a<-as.POSIXlt("2013-07-01 00:30:00")     b<-as.POSIXlt("2013-07-29 00:30:00")          aI<-as.numeric(a)     bI<-as.numeric(b)          times<-sample(seq(aI,bI,by=2),10000)     t<-sort(times)     class(t)<-c("POSIXt","POSIXct")     id<-seq(1,10050,20)     userID<-1     for(i in 1:200){        userID<-c(userID,sample(id[i]:id[i+1],50,replace=T))     }     userID<-userID[1:10000]          movement<-list(LETTERS[1:20],c("A","B","A","B","C","D","E","F","G","H","I")                   ,c("A","B","C","B","C","D","E","D","C","D","E")                   ,c("C","B","C","D","E","F","G","H","I","I"))          site<-character(10000)     for(i in unique(userID)){        p<-sample(movement,1)[[1]]        site[userID==i]<-p[1:length(userID[userID==i])]     }          table<-data.frame(userID=userID,site=site,time=t)


title: Confusion matrix of random Forest doesnot match predicted probabilities on train data

text: Based on an earlier question I balanced the classes such that the numbers in both classes is about similar. The random Forest gives next result:               > print(rFresult)          Call:      randomForest(formula = finresfh ~ ., data = rFdatasubset, importance = TRUE)                     Type of random forest: classification                          Number of trees: 500     No. of variables tried at each split: 14                  OOB estimate of  error rate: 35.53%     Confusion matrix:          1    2 class.error     1 1852  627   0.2529246     2 1022 1140   0.4727105      Prediction on the train set shows perfect separation in contrast to the confusion matrix:               > tab <- table(probability=round(predict(rFresult, newdata=rFdatasubset, type="prob")[,2],1), TRUE_status=rFdatasubset$finresfh)     > tab                TRUE_status     probability    1    2             0.1  978    0             0.2 1447    0             0.3   54    0             0.7    0   65             0.8    0 1551             0.9    0  543             1      0    3      The probability is estimated for the subjects to be in class 2. The "probability" table means the number of subjects with predicted probability level having a certain TRUE status. Can anyone explain why the estimated probabilities show a perfect separation but a totally different result in the confusion table?


title: how to partition the variance explained at group level and individual level

text: Based on results of Null model, some scholars suggest variance attributed to group level and individual level is computed using variance component output, whereas other suggest going the standard deviation scores. Is there consensus how to calculate ICC?


title: Pearson correlation of 3 variables

text: I am interested in whether or not a "correlation" of three variables is something, and if what, what would this be? _Pearson product moment correlation coefficient_ $$\frac{\mathrm{E}\\{(X-\mu_X)(Y-\mu_Y)\\}}{\sqrt{\mathrm{Var}(X)\mathrm{Var}(Y)}}$$ Now the question for 3 variables: Is $$\frac{\mathrm{E}\\{(X-\mu_X)(Y-\mu_Y)(Z-\mu_Z)\\}} {\sqrt{\mathrm{Var}(X)\mathrm{Var}(Y)\mathrm{Var}(Z)}}$$ anything? In R it seems like something interpretable:               > a <- rnorm(100); b <- rnorm(100); c <- rnorm(100)     > mean((a-mean(a)) * (b-mean(b)) * (c-mean(c))) / (sd(a) * sd(b) * sd(c))     [1] -0.3476942      We normally look at the correlation between 2 variables given a fixed third variable's value. Could someone clarify?


title: Confidence Intervals for the Mean of the Poisson Distribution 9-19 methods coded as an R package

text: I've trawled the internet looking for R packages to compute the confidence intervals for the mean of a Poisson distribution using the different methods cited by Patil and Kulkarni in their papaer "comparison of confidence intervals for the Poisson mean: some new aspects", REVSTAT vol. 10, no. 2, June 2012, pp 211-227. I have found in the R package, 'exactci', the routine poisson.exact, which gives a few methods. In the epitools R package there is pois.exact, pois.daly, pois.yar, and pois.approx. In the 'stats' R package there is poisson.test which uses the Clopper Pearson method for the confience interval of the mean. But in addition to these I was looking for those methods where the confidence interval is narrower, i.e. in a similar vein to the methods outlined for interval estimation of binomial proportions, as laid out in the paper by Agresti and Coull, 'Approxiamte is better than "Exact" for Interval Estimation of binomial proportions' the Amer. Stat. May 1998, vol. 52, no.2 pp 119-126. And on reading similar papers to the above on the Poisson distribution it was clear that for a good survey of confidence intervals, the score method and its variants should be included. It is these methods that do not appear to have been included in any R packages yet ? The closest I found was an R Package on Tolerance Intervals, ('Tolerance'). I don't claim to be totally clear on the distinction between CI and TI, but I think the latter guarantees something (e.g. measured length of a part) will lie in that interval to a given level of probability, but with the CI we are making a statement that the mean will lie in an interval to a given level of probablity. Is there an R package that does the CI calculation for the Poisson Distribution using these approximate methods (but not the normal dist. approx.) ? Thanks


title: Transformation of Martingale Differencce

text: We know non-zero measurable transformation of a dependent sequence is dependent. I want to ask if a non-zero measurable transformation of a sequence which is not martingale difference is not martingale difference?


title: Determining an optimum and distribution range from species abundance vs environmental gradient data

text: Hi and thanks for any help re the below This is a fairly complex question so I will attempt to ask it in a fairly basic manor. I have data on the abundance of 99 different species of estuarine macroinvertebrate species and the sediment mud content (0 - 100 %) in which each observation was obtained. I have a total of 1402 observations for each species (i.e. a massive dataset). Here is a subset of the raw data for one species to give you an idea of the data I'm working with (if I had 10 reputation points i'd upload a plot of real raw data:               Abundance: 10,14,10,3,3,3,3,4,5,5,0,0,0,0,0,0,0,0,0,0,0,0,0,6,6,6,0,0,0,0,12,0,0,0,34,0,0     Mud %:     0.9,4,2,10,13,14,6,5,5,7,22,27,34,37,47,58,54,70,54,80,90,65,56,7,8,34,67,54,32,1,57,45,49,4,78,65,45,35      The primary aim of my research is to determine an "optimum mud % range" (e.g. 15 - 45 %) and "distribution mud % range" (e.g. 0 - 80 %) for each of the 99 invertebrate species. As you can see the abundance data for the above species contains a significant number of zero values. Although this significantly skews any sort of model that I run on the data (i.e. GLM, GAM), even if I model the non-zero data only, the model for certain species does not fit the data at all well. So, my question is: what would be the best, most robust way to determine an "optimum" and "distribution" mud range for a given species, given that responses vary significantly between species?


title: How to interpret Box-Cox Transformation

text: I have used a Box-Behnken experimental design. I have a full quadratic model. However, I had to transform the response, $Y$ for the model to fit; I did this using a Box-Cox transformation with $\lambda=0.5$. For example, one of the regressions is like this: $$Y = 1.28 - 0.008X_1-0.025X_2-0.05X_3-0.13X_1^2-0.01X_2^2-0.006X_3^2\\\ +0.02X_1X_2-0.05X_1X_3+0.09X_2X_3$$ How do I interpret the terms?


title: Interpreting transformed dependent and independent variables

text: How would I interpret a transformed dependent variable (4th root) with some of its predictor variables transformed as well? In our study, we transformed our dependent variable to 4th root, $Y^\frac{1}{4}$, and one of our predictor variables to $\ln(X)$. How would I interpret that?


title: Observation vs model in a N-dimensional space

text: I have a set of observations which translates into $N_1(x_1,x_2,x_3,x_4)$ points where $x_1,...,x_4$ represent 4 parameters associated with each point. For example, the observed data would look like:               #x1 x2 x3 x4    # <-- Observ      1.2 0.5 46 10  # <-- 1      0.2 8.9 20 13  # <-- 2      0.8 4.5 18 9   # <-- 3      ...      0.2 4.3 55 7   # <-- N1      I also have 2 models, A & B, which attempt to reproduce that set of observations. Doing so each model gives me $NA(x_1,x_2,x_3,x_4)$ & $NB(x_1,x_2,x_3,x_4)$ points respectively:               #x1 x2 x3 x4   # <-- Model A      3.2 45 65 6   # <-- 1      0.2 8.9 20 14 # <-- 2      0.1 3.3 14 11 # <-- 3      ...      0.6 3.9 34 3  # <-- NA      and               #x1 x2 x3 x4    # <-- Model B      1.3 0.7 23 5   # <-- 1      1.0 7.4 17 9   # <-- 2      0.3 5.2 8 3    # <-- 3      ...      0.4 6.1 43 4   # <-- NB      **Important** : N1, NA and NB are **not** necessarily equal and there is **no** 1:1 correspondence between observed and modeled points. Think of it as a _cloud_ of points in a 4-dimensional space. Which statistic should I apply to try to quantify **which model is reproducing my observed set the best**?


title: Strategy to set the GBM parameters

text: I would like to know how you would set the parameters for the GBM model? 1) I could optimize the parameters sequentially. At first, using a large value for shrinkage and a small number of iterations (for faster computation), I would try to optimize n.minobsinnode, interaction.depth, bag.fraction, etc. Then using these parameters, I would scale the skrinkage to lower value and find the optimal number of trees. 2) All at the same time, using tons of scenario simulations.   For example, is it possible that the optimal interaction.depth to have different values for shrinkage 0.01 and 0.001?   If yes, it would not be fine to fix some parameters to optimize the others...


title: saturated structural equation model

text: I have a problem with SEM claiming to be saturated where it actually should test for missing covariances I am running a SEM composed in the following way. The ones with “2” behind are quardratic terms SEM1 = " d_birds ~ Q1 + MAT + MAP + TS + VelT + SamplingEffort + NetworkSize + plantp+ MAP2 + TS2 + PS2 Q1 ~ MAT + MAP + TS + VelT + SamplingEffort + NetworkSize + plantp+ MAP2 + TS2 + PS2 " When I run it in R using the code: SEM1fit = sem(SEM1,data = td, fixed.x = T) summary(SEM1fit) it tells me that the model is saturated with df=0 and therefore no chi square test are performed as well. This I don’t understand since I include no covarianses in the model (i.e. regarded as a fixed term = 0) My question in: • Do you think it is still informative if I base the model on R square and estimates (which are given by the SEM output) or is there a better wey to analyse the data Best regards


title: Why Use the Cornish-Fisher Expansion Instead of Sample Quantile?

text: The Cornish-Fisher Expansion provides a way to estimate the quantiles of a distribution based on moments. (In this sense, I see it as a complement to the Edgeworth Expansion, which gives an estimate of the cumulative distribution based on moments.) I would like to know in which situations would one prefer the Cornish-Fisher expansion for empirical work over the sample quantile, or vice-versa. A few guesses:   1. Computationally, sample moments can be computed online, whereas online estimation of sample quantiles is difficult. In this case, the C-F 'wins'.   2. If one had the ability to forecast moments, the C-F would allow one to leverage these forecasts for quantile estimation.   3. The C-F Expansion can possibly give estimates of quantiles outside the range of observed values, whereas the sample quantile probably should not.   4. I am not aware of how to compute a confidence interval around the quantile estimates given by C-F. In this case, sample quantile 'wins'.   5. It seems like the C-F Expansion requires one to estimate multiple higher moments of a distribution. The errors in these estimates probably compound in such a way that the C-F Expansion has a higher standard error than the sample quantile. Any others? Does anybody have experience using both of these methods?


title: linear regression using heteroskedasticity robust standard errors in R

text: I want to perform an OLS regression on time series data using heteroskedasticity robust standard erros. So far i can come up with this:               model <- lm(I(y[2:T] - y[1:T-1]) ~ y[1:T-1])     regress <- coeftest(model, vcov=vcovHC(model, type = "HC0"))      Biscally this seems to work but when excluding the intercept term in the regression formula , i.e.               model <- lm(I(y[2:T] - y[1:T-1]) ~ 0 + y[1:T-1])      it yields NAs for all values except the Estimate so i am unable to compute distinct t statistics. I would be grateful for any advice on solving this.


title: Excel as a statistics workbench

text: It seems that lots of people (including me) like to do exploratory data analysis in Excel. Some limitations, such as the number of rows allowed in a spreadsheet, are a pain but in most cases don't make it impossible to use Excel to play around with data. A paper by McCullough and Heiser, however, practically screams that you will get your results all wrong -- and probably burn in hell as well -- if you try to use Excel. Is this paper correct or is it biased? The authors do sound like they hate Microsoft.


title: Smoothness of a surface

text: I am currently working on a model which takes two parameters and produces a measurement statistic. Think of it as Z = f(X,Y). Z is a matrix of my statistics and I am creating a surface plot of it in matlab. Basically, I am looking for a mathematical/analytical way of determining if the surface is smooth, or if it is jagged. Do large values tend to be clustered together or are they dispersed throughout the matrix? - that is my question. Basically, how mixed up are the values of my matrix? I need to run the model over different parameter sets and I want to be able to analytically determine which one of my surfaces is the smoothest, has the greatest clustering of large values, and ideally, has no negative values. Any help will be greatly appreciated and please let me know if you need any further information. Cheers


title: What are the assumptions for MANOVA and how I should apply it in R?

text: I would like to analyze effect of a treatment on my patient samples by MANOVA in R. I have measurement of 3 different protein level for treated and untreated patients. Now I would like to analyze effect of my treatment on the whole different protein measurements and not each individually (sort of summing them up together). What are the assumptions for MANOVA and how I should apply it on my dataset in R? An example line of my dataset:                       CD13    CD68     AT1     treatment     sample1 0.065   0.82     0.488   +     sample2 0       1.58     0.47    -      I have around 30 replicates for each treated and untreated patient group. Thanks


title: What is the difference between a loss function and decision function

text: I see that both functions are part of data mining methods such as Gradient Boosting Regressors. I see that those are seperate objects. But how is the relationship between both in general?


title: Regression with dummy for only one specific date

text: I am running a regression where I regress return of a stock on the market return. There is a specific date in my sample on which some event occurred. I believe the effect of this event to last only on that date. I want my model coefficients to be free of this effect because I want to predict what the return would be, had this event not occurred. I am thinking of including a dummy for that date (i.e. it is one for that date and zero for the rest of the sample). Is this an appropriate way to handle this issue? Thank you.


title: posterior covariance for linear inversion

text: I am performing a physical experiment where the unavoidably inaccurate setup likely creates errors in the inversion of the data for the model. Let's assume a linear relationship between data $\mathbf d$ and model $\mathbf m$ as follows: $$\mathbf d = \mathbf W \mathbf m,$$ where $\mathbf W$ is the operator (encoding the setup of my physical experiment). Standard linear inversion yields the model $$\mathbf m = \mathbf {C_m}^{-1} \mathbf W^\prime \mathbf {C_d}^{-1} \mathbf d$$ with the data covariance $\mathbf {C_d}$ and the model covariance $$\mathbf {C_m} = \left( \mathbf W^\prime {\mathbf {C_d}^{-1}} \mathbf W \right)^{-1}.$$ So far, so good. The operator $\mathbf W = \mathbf W (\mathbf p)$, however, depends in a non- linear way on a parameter $\mathbf p$, but a linearization, say $$\mathbf V = \dfrac {\partial \mathbf W} {\partial \mathbf p} \mathbf {p_0},$$ where $\mathbf {p_0}$ is the a-priori known (measured) parameter, will do. Now, how do I continue from here? That is how do I account for any (small) parameter uncertainty in the inversion for the model and model covariance?


title: Which the best programming language to handle time series?

text: I work with CSV files containing 10-minute data of temperatures of some years of a station, and I usually perform some calculations: mean, average, correlations between two stations, correlations in different seasons, correlations day/night, ... I program this calculations by using plain Python or Microsoft Excel VBA, by using loops and simple calculations. It's not very hard to code, but I wonder if there is a Python library or any other programming environment which helps me to save time and lines of code. I've heard about ScyPy and R, and I would be very grateful if you could share any experience with these tools. Thank you very much in advance.


title: Test for convergence within Gibbs sampler

text: I am running a Gibbs sampler for Multivariate Normal times Inverse Wishart posterior distribution with missing data imputation step. I am trying to check if my step of simulating covariance matrices from Inverse-Wishart converges. How would I go about it? In MCMC I usually run Gelman-Rubin diagnostic test. Thank you!


title: Specifying a difference in differences model with multiple time periods

text: When I estimate a difference in differences model with two time periods, the equivalent regression model would be a. $Y_{ist} = \alpha +\gamma_s*Treatment + \lambda d_t + \delta*(Treatment*d_t)+ \epsilon_{ist}$   * where $Treatment$ is a dummy which is equal to 1 if the observation is from the treatment group   * and $d$ is a dummy which is equal to 1 in the time period after the treatment occured Thus the equation takes the following values.   * Control group, before treatment: $\alpha$   * Control group, after treatment: $\alpha +\lambda$   * Treatment group, before treatment: $\alpha +\gamma$   * Treatment group, after treatment: $\alpha+ \gamma+ \lambda+ \delta$ Hence, in a two period model the difference in differences estimate is $\delta$. But what happens concerning $d_t$ if I have more than one pre and post treatment period? Do I still use a dummy that indicates whether a year is before or after the treatment? Or do I add year dummies instead without specifying whether each year belong to the pre or post treatment period? Like this: b. $Y_{ist} = \alpha +\gamma_s*Treatment + yeardummy + \delta*(Treatment*d_t)+ \epsilon_{ist}$ Or can I include both (i.e $yeardummy +\lambda d_t$)? c. $Y_{ist} = \alpha +\gamma_s*Treatment + yeardummy + \lambda d_t + \delta*(Treatment*d_t)+ \epsilon_{ist}$ In conclusion, how do I specify a difference in differences model with multiple time periods (a,b or c)?


title: MODMED in PROCESS generates insignificant interaction term yet significant conditional effect: how to interpret?

text: I have a moderated mediation model with a strong mediation effect, but with confusing results for the moderating effect (using PROCESS). Although the interaction effect in the regression model is insignificant, the conditional indirect effect does show that the mediated effect is conditional on the moderator (values above the 10th percentile are significant). To get a better understanding of whats going on I also used the MODMED macro which enables me to use the Johnson Neyman Technique. What this shows is that the moderator has a significant effect at values of the moderator higher than 3,38 (i'm using a 5 point scale). Yet again, the interaction in the regression model is insignificant. I've been going through the literature but am unable to find an answer. Would you be able to help me?


title: Aggregating samples for clustering time series

text: I want to cluster a set of 512 time series. The time series have sampling intervals of 1 day over a time period of 5 years. Thus, each time series consists of about 1800 samples. However, many of the samples are zero. So, my question is: Is it ok to aggregate the samples to larger sampling intervals, e.g. days or even years, despite the information loss that comes with aggregation? Isn't aggregation some form of data tampering? If not, how do I chose the appropriate temporal unit for aggregation?


title: Model selection using a standard p-value approach

text: I need to perform model selection using a standard p-value approach. Using logistic regression we would like to compare the following models:               Y =       A + B + C + D     Y = A*B + A + B + C + D     Y = A*C + A + B + C + D     Y = A*D + A + B + C + D      I performed the following analysis in R:               first  <- glm(Y~G+N+E+C,     family="binomial")     second <- glm(Y~G*N+G+N+E+C, family="binomial")     third  <- glm(Y~G*E+G+N+E+C, family="binomial")     fourth <- glm(Y~G*C+G+N+E+C, family="binomial")          summary(first)     summary(second)     summary(third)     summary(fourth)          anova(first, second, third, fourth, test="Chisq")      However, I think I do not have the right output here for a model selection based on p-values?               anova(first, second, third, fourth, test="Chisq")     # Analysis of Deviance Table     # Model 1: Y ~ G + N + E + C     # Model 2: Y ~ G * N + G + N + E + C     # Model 3: Y ~ G * E + G + N + E + C     # Model 4: Y ~ G * C + G + N + E + C     # Resid. Df | Resid. Dev | Df | Deviance | Pr(>Chi)     # 1   595   |  609.90    |                          # 2   594   |  609.90    | 1 | 0.000169 |   0.9896     # 3   594   |  609.81    | 0 | 0.087775 |             # 4   594   |  609.90    | 0 | -0.085001|        So, how to perform a model selection here, using a standard p-value approach?


title: Should I use an average ECDF?

text: This relates to a previous question of mine which didn't gain many responses, perhaps because it wasn't very clear nd well written. I hope this time I will be more accurate and get your much appreciated assistance. I am analyzing results of a biological experiment. The results given as a single value ( non-negative integer) per genomic position. I am interested in valleys, or local minimam over this series of values. I wish to control the false positives rate and get the significance for each local minima. I can shuffle the raw data which was used to produce the data. So what I do is to shuffle the raw data, create the new series of values,search for all local minima and keep their values. Now, I have something like this:               data_set  local_minima_values     =============================     true_data 4 9 1 27 12 0 0 2 5 32 0 1 5 70 2     sim_1 14 25 94 59 32     sim_2 52 0 14 74 82 12 54     ...      Note the number of local minima naturally varies between simulations. So, my idea was to calculate an ECDF for each simulation and then combine those ECDFs into a single "average ECDF" which represents the null hypothesis. Then, I can assign a p-value for each local minima from the true data, and see how significant ('surprising') it is. My questions are:   1. Does this make sense?   2. How do I create an average ECDF? I can't just merge the values from all simulation together and get and ECDF for this merged set, since the number of minima found in each simulation differs, and I think all simulations should have the same contribution to the average ECDF, or am I wrong?   3. How should I take the number of simulations (shuffles) into account? Thanks, Dave p.s. I'm working with R.


title: Graphing real-time data from a text file

text: I have a process which writes statistics from a server system to a file each second in this format:               label1 label2 label3     344 666 787     344 849 344     939 994 344      There are a number of different values which I need graphs for, and each value is added to the bottom of the file each second. I am looking for a nice way to show these numbers graphically, preferably a program which auto-updates each second and shows the different graphs. Can anyone recommend such a program for mac? Preferably a free one :)


title: Essential papers on matrix decompositions

text: I recently read Skillicorn's book on matrix decompositions, and was a bit disappointed, as it was targeted to an undergraduate audience. I would like to compile (for myself and others) a short bibliography of essential papers (surveys, but also breakthrough papers) on matrix decompositions. What I have in mind primarily is something on SVD/PCA (and robust/sparse variants), and NNMF, since those are by far the most used. Do you all have any recommendation/suggestion? I am holding off mine not to bias the answers. I would ask to limit each answer to 2-3 papers. P.S.: I refer to these two decompositions as the most used _in data analysis_. Of course QR, Cholesky, LU and polar are very important in numerical analysis. That is not the focus of my question though.


title: SPSS: Pearson's r not significant but confidence intervals do not include 0

text: I have calculated the value of Pearson's _r_ between two variables in SPSS, two tailed. Sample size is 81. The _r_ value is -.21. The _p_ value is .06. When I compute the 95% confidence intervals using bootstrapping with 2000 samples and the Bias corrected accelerated (BCa) option checked, I get confidence intervals [-0.36, -0.04] (does not include 0), indicating that in fact the result is significant at _p_ < .05. How do I explain the discrepancy between the _p_ value and the confidence intervals? I calculated the confidence intervals again with "percentile" instead of BCa checked and the confidence intervals were [-.378, .009] (does include 0), which makes more sense to me. How can I explain the apparent incongruence between the _p_ value and confidence intervals I attained when using the BCa method? A secondary issue is how I would report this (BCa) result in APA style.


title: Confidence interval for N, based on binomial

text: I've searched the web and these forums but couldn't find a subject that resembled my problem. If I missed something then I'd like to apologize upfront. My problem is this. I have a given probability p, given by an outside agency and the number of successes n from an owned datasource. I'm interested in the population total N and the confidence interval of N. Background: I know how many passengers there are for a specific bus agency using 2013 data (as in: I have the (micro)data from that specific company). Assumption: this is my n. I also know from transportation figures (from an outside agency) that this specific bus company had a market share (as measured in passengers) of 30% in 2013. Assumption: this is my p. This assumption brings with it certain other assumptions, that is fine. So now the question is: how many passengers were there in total (assuming this can not be taken from the transportation figures)? ie. what is N and what is the 95% confidence interval surrounding N? So my questions are: 1) Can I just use the relation: p=n/N, therefore the point estimat for N is: N=n/p? I'd assume so. 2) Can I just use the Clopper-Pearson interval to calculate an interval for p, then using that to get to an interval for N? 3) Suppose the 95% lowerbound for p=plb. Is the 95% lower bound for N then n/plb? Same logic for the 95% upper bound. 4) When calculating the Clopper-Pearson interval, can I just use: n=n, p=p, N=n/p? 5) If any of the answers to my questions 2-4 is 'No', then what is the correct approach? I realize that using this approach, N is not necessarily an integer. Is that problematic?


title: power of poisson test with large sample size

text: I would like to compute the power of two distributions which are poisson with $n_1$ being large and $n_2$ being medium sized. I have used the formula given by Krishnamoorthy and Thomson for the C-test (see below). $$\sum^\infty_{k_1=0}\sum^\infty_{k_2=0}\frac{e^{-n_1\lambda_1}(n_1\lambda_1)^{k_1}}{k_1!}\frac{e^{-n_2\lambda_2}(n_2\lambda_2)^{k_2}}{k_2!}I[P(X_1\geq k_1|k_1+k_2, p(\lambda_1/ \lambda_2))\leq\alpha]$$ However, as $n_1$ or $n_2$ get large the exponential function goes to $0$ which makes taking the sum of the probability rather pointless. Is there anyway I can compute the power with having an large uneven sample size? Thanks a lot for your help


title: Log approximation

text: Can somebody help me with this approximation: > as $m$ gets really large, $\dfrac{r}{m}$ gets really small and hence > $\log\left(1+\dfrac{r}{m}\right) \sim \dfrac{r}{m}$. I don't understand how did $\log\left(1+\dfrac{r}{m}\right)$ approximately equals to $\dfrac{r}{m}$? Please clarify.


title: Averaging z scores when doing a "group by"

text: I have a dataset where each row is an hourly measurement of certain fields (columns). For each column I then add another column that is its respective _z_ score relative to the entire population. If I then group by month by calculating the aggregate mean for all merging rows, are my _z_ scores still valid, or do I need to recalculate the _z_ score for each month? **EDIT** I realize I wasn't clear enough in my question. I have a dataset like this of hourly weather measurements going back several years: ![enter image description here](http://i.stack.imgur.com/Z8tIR.png) In panda I groupby and average each respective hour (So I average and group all the Jan 1st 2pm rows over all years... etc.)               df['hour'] = df['Date'].apply(lambda t: "%d-%d-%d" % (t.month, t.day,t.hour))     g = np.round(df.groupby('hour').mean(),decimals=4).reset_index()      which returns this: ![enter image description here](http://i.stack.imgur.com/70Il0.png) I then merge the grouped dataframe (g) onto the original df adding a column for each avg, reset the index to be of Date, and drop the 'hour' column from the grouped df:               new_df = pd.merge(left=df, right=g, on='hour', suffixes=('','_avg') )     new_df.set_index('Date').drop('hour',axis=1,inplace=False).head()      This returns: ![enter image description here](http://i.stack.imgur.com/zmdYV.png) What I would like however is to also add a "z score" column for each respective month/day/hour using the standard deviation of the entire population. I'm actually now sure how to do this but I think it would be similar to how the other "_avg" columns were computed. Assuming that is all possible and each row-hour has additional z-score columns for each metric this finally leads me to my original question... Later on outside the context of python/panda let's say I take only one year of the data, groupby month, and aggregate by mean. Will those averaged z scores by month still be representative of the original entire population? For example let's say I took the averages of June 2011 and found the temperature z-score to be 1.2. Can I say that the average temperature of June 2011 was 1.2 standard deviations away from the normal (total population) mean?


title: Time Series Data Mining Library?

text: Can anyone recommend a library for time series data mining tasks other than predictive modeling and statistical analysis? There seem to be a number for these purposes (e.g., Gretl), but nothing for the tasks of:   * Classification   * Clustering   * Subsequence search   * Anomaly detection   * Motif discovery   * Indexing I can't be the first person to want to use, e.g., iSAX, but Googling hasn't turned up anything other than a single Java library. Links to anything that one could interface with C/C++, Python, or Julia are especially welcome, but any help would be appreciated. Assurance that such a library does not exist would be helpful as well, since I would then know it wouldn't be a waste of time to write one...


title: How to avoid NaN in using ReLU + Cross-Entropy?

text: ReLU has a range of [0, +Inf). So, when it comes an activation value z=0/1 produced by ReLU or softplus, the loss value computed by cross-entropy : loss = -(x*ln(z)+(1-x)*ln(1-z)) will turn to NaN. As i know, my variables are run in theano.tensor type which cannot be modified after defined. So, I cannot simply change the values in z which are 0/1 into some approximate (such as 0.001/0.999). How to avoid this annoying situation without replacing the cross-entropy by MSE?