title: Restricted permutations

text: I have some multivariate data and want to investigate the effect of some environmental gradient. I want to use capscale but I don´t know how to deal with the permutation scheme. I have made up some artificial data, with 20 sites along a gradient ("env"):               ######### create some species data along a gradient     df <- structure(list(site = 1:20,                           sp1 = c(7L, 4L, 2L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L),                           sp2 = c(1L, 2L, 4L, 7L, 8L, 7L, 4L, 2L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L),                           sp3 = c(0L, 0L, 0L, 0L, 0L, 1L, 2L, 4L, 7L, 8L, 7L, 4L, 2L, 1L, 0L, 0L, 0L, 0L, 0L, 0L),                           sp4 = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 2L, 4L, 7L, 8L, 7L, 4L, 2L, 1L, 0L),                           sp5 = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 2L, 4L, 7L, 8L),                           sp6 = c(0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L),                           sp7 = c(0L, 0L, 0L, 0L, 0L, 0L, 2L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L),                           sp8 = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 3L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L),                           sp9 = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 4L, 1L, 0L, 0L)),                      .Names = c("site", "sp1", "sp2", "sp3", "sp4", "sp5", "sp6", "sp7", "sp8", "sp9"),                      class = "data.frame", row.names = c(NA, -20L))          # add some linear responses     df$sp6 <- round(seq(1, 8, 7/19), digits = 2) # linear response     df$sp7 <- round(seq(1, 4, 3/19), digits = 2) # no so strong linear response     df$sp9 <- round(seq(1, 6, 5/19), digits = 2)          # gradient     df$env <- 1:20      If I sampled only once I would do something like this:               # db-RDA sampled at one time     require(vegan)     mod <- capscale(df[ ,-c(1, 10)] ~ env, data = df, distance = "bray")     anova(mod, by = "terms", step = 999)  # assess the "significance" of contraining variable     plot(mod)      Now imagine I sampled the same data trice, but in three different months:               # now we replate exactly the same data 2 more times      repdf <- rbind(df, df, df)       repdf$time <- rep(1:3, each = nrow(df))     repdf$site <- factor(repdf$site)      If I would use unrestricted permutations, then this won't capture the repeated measures and the p-values would be to low. I could restrict the permutations within each sites (using strata = site in vegan), but this destroys only the temporal effect and yields to a p of 1 (because every permuation is the same):               repmod <- capscale(repdf[ ,-c(1, 11, 12)] ~ env, data = repdf, distance = "bray") # db-RDA     anova(repmod, by = "terms", strata = repdf$site, step = 999)        **My question:** How should I restrict the permutations assessing the effect of the gradient taking this temporal correlation into account? What permutationscheme should I use? **Some ideas:** a) Permute the strata ( = sites), but not within the strata. This will destroy the env gradient, so the p-value is only determined by this. b) Include "time" into the model (with interactions) and the run for time- effect a different permutation-scheme (permute within sites) than for env (permute sites, but not within sites). I know about the permute-package and can incorporate it into permutest.cca, so my question is more theoretical.


title: Gower's dissimilarity index

text: I would like to ask a question about Gower similarity/dissimilarity index. Is it ok to use the Gower dissimilarity measure with Ward linkage clustering? I was reading that the Gower similarity index should not be used with Ward linkage because the index is not metric. I was wondering if this is only the case for the similarity and not for the dissimilarity index that can also handle odrinal variables?!


title: Correlating two distributions to a third

text: My dataset consists of 4-million records, each a combination of age, race, and a discrete observed value. The goal is to find a correlation between the combination of age/race and the observation, so that given a distribution of ages and races, a distribution of expected observations can be calculated. I am looking into covariance matrices, but I am unsure how to handle the categorical data in race. What is the correct direction to begin investigating?               Age     Race     X     24      Black    32     26      White    37            ...      edit for newcomers: once I have obtained a linear model, I need to use distributions as inputs to receive a distribution as an output. I will not be predicting X for a single individual, but rather I need to determine the distribution of X for a given distribution of age/race. are there methods of transforming a linear model into something that can accomplish this? edit again: in the simplest terms... give the data that I have (Age, Race, Observed X) I need to come up with a model whose inputs are a distribution of ages and a distribution of races (i.e. not a single age/race combination) and whose output is the expected distribution of X. I have already formed a linear model using R that can give me the expected X given a single value of Age and Race, but the model needs to work on population-level inputs/outputs, not individual-level.


title: Correct order of performing imputation and variable selection

text: This is a general question about performing data analysis. I have a data set with ~1000 sample size and 200 features. Some of features have more than 50% missing or even higher. The missing pattern is unclear to me. I plan to do regression. Any one can provide some reference where should I start and what is the correct procedure to do? Suppose the variables needs no transformation. What might be done before regression include: A. looking at those with very high missing rate(>90% maybe?) and discard them B. impute the remaining variables C. perform variable selection Is it reasonable to do step A? If yes, how to determine the cut-off rate? Between imputation and variable selectin, which should be done first? Thank you for any reply.


title: AB Testing other factors besides conversion rate

text: So I've been learning about AB Testing and have used it to examine the form conversion rate of two different forms. However, I'm curious about testing whether a form with ads will generate more revenue than a site with just a form and no ads. So the two things I'm comparing:               - form with ads     - form with no ads      I'm trying to see which produces more revenue. So traditionally, a test of conversion rates would involve using a chi-square or fishers exact test to examine the values in a contingency table. I'm just not sure how to approach this question when it comes to revenue.                          Revenue     Conversions     Site 1 -    $500         100      Site 2 -    $400         70      Is this no different that a tradition ab test of conversion rates? Or could I just test wether a form made money or not.                         made money       made no money     form 1       50                200     form 2       5                 250


title: Logistic regression linear hyperplane

text: I'm performing logistic regression for binary classification in SAS and it outputs all of the coefficients to the logit model. I was wondering how you would take those coefficients and use them to convert the model into a linear hyperplane in the original space (i.e., hyperplane that splits the space into two classes, positive and negative).


title: How can I determine statistical significance in an A/B test in which the KPI is dependent upon two variables - one bernouli and one continuous?

text: In my work in online marketing, we frequently run A/B or multivariate web page tests. The Key Performance Indicator for these tests is overall Revenue. The treatment that nets the most revenue, either by influencing more conversions to sale, or by influencing a user to spend more per transaction, or a combination of these two, is the winner. However, how can we determine we have a statistically valid Revenue impact in a particular test? The standard-issue significance test does is not useful here, since the Revenue impact is obviously due to a combination of probability to convert + dollar value of conversion. We have played with a version of Overall Evaluation Criterion from the Taguchi playbook, but honestly I'm not statistically sophisticated enough to know if this approach is valid. I appreciate any suggestions.


title: Example Of Strict von Neumann Inequality

text: Let $r(\pi, \delta)$ denote the Bayes risk of an estimator $\delta$ with respect to a prior $\pi$, let $\Pi$ denote the set of all priors on the parameter space $\Theta$, and let $\Delta$ denote the set of all (possibly randomized) decision rules. The statistical interpretation of John von Neumann's minimax inequality states that $$ \sup_{\pi\in\Pi} \inf_{\delta\in\Delta} r(\pi, \delta) \leq \inf_{\delta\in\Delta}\sup_{\pi\in\Pi} r(\pi, \delta), $$ with strict equality guaranteed for some $\delta'$ and $\pi'$ when $\Theta$ and $\Delta$ are both finite. Can someone provide a concrete example where the _inequality_ is strict?


title: R language how to access wts of a class objects nnet

text: i have built a neural net class object using http://www.inside-r.org/r-doc/nnet/nnet and i am trying to access values returned by that object. Values section of page http://www.inside-r.org/r-doc/nnet/nnet says that "object of class "nnet" or "nnet.formula". Mostly internal structure, but has components wts the best set of weights found". But how to get those weights? i tried below commands but didnt get required results :( any help?                   neuralnetobject$wts              neuralnetobject$finalModel$wts


title: Is there a non-parametric test for whether a series' unconditional variance has changed over time?

text: I'm not sure how this would look. It seems you've have to specify regime breakpoints within the data, right? Or is there some rolling method that could compare window sizes against each other? Also key, is it possible to do this without specifying a model?


title: interpreting ancova and extracting estimated means

text: Imagine the following data:               ds <- data.frame(x=1:10,y=1:10,z=rep(c("A","B"),each=5))      the means for groups in `z` are:               library(plyr)     ddply(ds, "z", function(x) mean(x$y))     #  z V1     #1 A  3     #2 B  8      do a couple of models:               m1 <- glm(y ~ x , data = ds)     m2 <- glm(y ~ z , data = ds)      In `m1` the intercept is zero and the estimate for x is 1 as we would expect form the created data. Intercept meaning the value of `y` for `x=0`. In `m2` the intercept is `3` and the estimate for `B` is `5`. This is because of the contrasts between the intercept with the estimate i.e. the estimate for `B` is 8. These results make sense in terms of the means of my data as calculated separately above. so for `m1` I interpret intercept as `y where x=0` and `slope` and for `m2` as this is an anova, `intercept` as mean of first level and the rest are differences between that level and the first. However in the following models `m3` and `m4`               m3 <- glm(y ~ x + z , data = ds)     m4 <- glm(y ~ x * z , data = ds)          # m3     #Coefficients:     #              Estimate Std. Error    t value Pr(>|t|)         #(Intercept)  1.123e-15  6.201e-16  1.812e+00    0.113         #x            1.000e+00  1.720e-16  5.814e+15   <2e-16 ***     #zB          -7.628e-17  9.880e-16 -7.700e-02    0.941              # m4     #Coefficients:     #              Estimate Std. Error    t value Pr(>|t|)         #(Intercept)  1.123e-15  8.714e-16  1.289e+00    0.245         #x            1.000e+00  2.627e-16  3.806e+15   <2e-16 ***     #zB          -9.349e-17  2.305e-15 -4.100e-02    0.969         #x:zB         3.130e-18  3.716e-16  8.000e-03    0.994        The coefficients for intercept and `x` are the same as `m1` but how do I interpret the estimate for `B` and `A` as the estimate for `B` is `-7.628e-17`. My question is how do I interpret the categorical estimates in `m3` and from the output extract the estimated means, which I would report? and for `m4` is the slope for `x:zB` `1.000e+00 + 3.130e-18` and for `A` `1.000e+00 - 3.130e-18`?


title: Stacking/ensembling models with caret

text: I often find myself training several different predictive models using `caret` in R. I'll train them all on the same cross validation folds, using `caret::: createFolds`, then choose the best model based on cross-validated error. However, the median prediction from several models often outperforms the best single model on an independent test set. I'm thinking of writing some functions for stacking/ensembling caret models that were trained with the same cross-validation folds, for example by taking median predictions from each model on each fold, or by training a "meta-model." Of course, this might require an outer cross-validation loop. Does anyone know of any existing packages/open source code for ensembling caret models (and possibly cross-validating those ensembles)?


title: Gradient in GMM estimation

text: I have a question that might be trivial but I have not much knowledge on that method: I want to estimate a structural model with GMM and my model works in the sense that it estimated the right coefficients of simulated data. This works fine without adding a gradient (I'm using "gmm"-package in R), just the vcov-matrix of the parameters differs. So I'm a little bit confused which matrix is "the right" one and if I have to add the gradient or not. Thanks for your help!


title: Poisson and Exponential Distribution- Is the following question correct?

text: Accidents occur with a Poisson distribution at an average of 4 per week. i.e., $\lambda= 4$.   1. Calculate the probability of more than 5 accidents in any one week.   2. What is the probability that at least two weeks will elapse between accident? Query: Is is necessary to be clear that Part 2 is based on exponential distribution or time is exponentially distributed or something like this? Note: The question is based on integrating knowledge of the Poisson and exponential distributions.


title: What is exactly code vector and quantization vector of self organizing map?

text: I am trying to understand code vector in self organizing map. Could anybody explain me intuitively what it is exactly?


title: How to assess whether experimental measurements obtained from different technicians are biased?

text: Suppose I have a list of measurements from an experiment; for example, > 34 31 55 18 19 22 44 48 23 . . . But I then learn that these experiments were conducted by two different technicians, so I embolden those conducted by tech #1, and I leave untouched the measurements of tech #2: > **34** 31 55 **18** **19** 22 44 48 **23** . . . I want to know if the bold numbers are somehow biased relative to the unbolded numbers. How should I measure this? My intent is to compare the mean and variance of the overall sequence to the means and variances of each of the subsequences. But I'm a novice at statistics, so perhaps this naive approach is too shallow?


title: Good PCA examples for teaching

text: I'm teaching linear algebra to a class of engineers, social scientists and computer programmers. We just did singular value decomposition, and we have an extra day, so I thought I'd talk about the relation between singular value decomposition and principal component analysis. I have the theory part of the lecture written just fine, but am finding a hard time finding good examples to use. Here are the constraints:   * I want to show pictures. Ideally, the graphics should work well by themselves: Axes and data points in the scatter-plot should be labelled. English words are better than Latin species names.   * The question being studied should be interesting. Morphology of Nigerian fish, while important, is not a good way to grab a class's attention.   * In contrast to the preceding bullet point: Nothing on human racial differences; nothing on intelligence testing. That would lead to a lively discussion which would have nothing to do with the mathematical techniques.   * The mathematical analysis method should be basically pure PCA. The DW-NOMINATE project, while awesome, uses PCA as a starting point which is followed by a much more complicated hill climbing algorithm. I would think this would be easy. I can easily think of a dozen fun analysis projects I could do if I had the time to gather the data: Take the Pew Research polls and see whether PCA recovers the social policy/fiscal policy axis beloved of libertarians. Take a dozen measurements of typical dog breed physical characteristics and see if PCA can find the "sheep dog" cluster. Etcetera, etcetera... I'm looking for someone else who has already done the work so I can show it off. I'm afraid the timeline here is rather tight: I lecture tomorrow afternoon (Monday). I've spent most of the weekend googling PCA studies in different interesting fields and finding over and over that they aren't good fits.


title: Residual plot for simple linear regression

text: I'm trying to analyze a part of European Social Survey data.   The outcome has to be treated as continuous, however it can take only discrete values from 1 to 6. The predictor variable is age. I fitted a simple linear regression model and now I need to check model's assumptions.I have this residuals vs. predictor variable plot and a fit plot. Could you please help me to interpret them? As I understand SD of the error terms is not constant, but what else can we read from these plots? Thank you very much! ![this](http://i.stack.imgur.com/d4d3s.jpg) ![enter image description here](http://i.stack.imgur.com/xfKWo.jpg)


title: How should I choose SNR value that the object is recognisable by the human eye from background noise?

text: I have some images that are made of an object and some amount of background noise. I have sampled the object for example at 40%. So the object is becoming pixelated, then I inject some amount of noise to the object. I know that formula for SNR is:               SNR = 20log10 (Signal/Noise)      According to that formula the amount of my signal should be greater than the amount of noise in order for the signal to be recognisable. So when I’m sampling the object at 40% how much should be the background noise in order to recognise the signal? Or how should I relate my sampling rate of the object to the background noise?


title: Missing data due to absent parent

text: I am using the following regression: $$\text{Test score} = \beta_0+\beta_1\text{Mother's employment}+\beta_2\text{Mother's education}$$ where "Mother's employment" is a set of dummy variables indicating whether the mother works more than 35 hours a week, is unemployed or **is absent** , and "Mother's education" is also a set of dummy variables indicating if the mother has a high school diploma, a college degree or a PhD. If the mother is absent, then "Mother's education" is not applicable, i.e. there is no answer. How do I deal with this in **Stata**? Mean imputation? How do I do that with dummy variables?


title: Comparing two binary variables of unequal sizes

text: I measured a binary variable from two different populations, and now I'm trying to find out whether the different populations differ with regards to this variable. I could use a Chi-Square test, but that would necessitate that both populations have the same length. Is there an appropriate test for these circumstances? Thank you.


title: Are rejection regions always open/closed?

text: Does there exist any agreement on what must rejection regions look like topologically? If we identify the region of "acceptance" with the corresponding confidence interval (or confidence region in dimensions greater than 1) then, according to this question, it seems that rejection regions must be open. On the other hand, the rejection region given by the Neyman-Pearson lemma is usually defined in terms of a non-strict inequality, so at least in the continuous case it will be closed. (I know that this is a somewhat bizantine question, at least from a practical point of view.)


title: Initial steps with logistic regression

text: I have a dichotomous variable which I am trying to predict. In order to do this I have about 2000 variables out of which 1100 are continuous variables (real values between 0.0 and 1.0) and the rest are categorical variable (0 or 1). I am a total newbie and would like to know how to go about it. I am not aware of the relationships between these variables. Roughly the following are the steps that I am planning to take:   1. Should I try to use correlation analysis between each variable and the target variable?   2. Should I try some kind of clustering and remove some outliers and try to come up with a model.   3. What other steps should I take to identify relationships between variables and come up with a model? Thanks


title: Feature selection for disease classification based on tests

text: I have a dataset of around 100 different subjects Some of them have a disease, some do not (roughly 60:40 disease:no disease) They are subjected to a battery of 15 tests, to see if they are outside "normal" ranges. Just plotting the values for the different tests for disease vs. non-disease as different colours (using `matplot()` in R), I can see that the different groups follow distinct patterns across the different features. I then cluster the different groups (using `hclust()` in R) and if I cut the tree to make two clusters, the two groups separate fairly well into different clusters. My aim is to devise a set of rules from these tests, so if we test a new patient, we can decide whether or not they have a disease. So I need a classifier, to decide these rules, i.e. to work out which features to use, and what score cutoffs. What do people recommend?


title: What is the probability P(X > Y) given X ~ Be(a1, b1), and Y ~ Be(a2, b2), and X and Y are independent?

text: Given two independent variables with Beta distribution, $X \sim \text{Be}(a_1, b_1)$ and $Y \sim \text{Be}(a_2, b_2)$, how do you find the probability that the value of X is greater than the value of Y for a given observation? Does this probability have a name that I'm just blanking on?


title: Where can I find ADF test library or source code from c#

text: I would like to test for stationarity in cointegration. I intend to use an augmented dickey fuller test. However, I need one for c# - either a library or the source code. Or is your have source in a similar language c, java I can implement. Thanks


title: Is LOESS the appropriate way of visualizing my RT data?

text: I am conducting a Lexical Decision Task where my dependent random variable is Response Time (RT). My experimental design consists of 5 blocks of a 100 trials each. In each block, 50% of trials correspond to one experimental condition, and the other 50% to the other experimental condition. Order of presentation or trials in randomized. I am using the R programming language to do my analysis When I plotted the average RT per block, per condition, I saw that the difference between my conditions was quite prominent, whereas I expected differences to be small at first, and then become larger as blocks advanced. In order to see more in detail how the RT were changing across trials, I plotted RT across trials instead of blocks. Now, since there were a lot of data points, I had to use a smoothing technique called LOESS. This indeed smoothed the data pretty nicely but: 1) Apparently it made me loose data. When using the stat_smooth() function in R, it warned me of missing values:               1: Removed 31 rows containing missing values (stat_smooth).      2: Removed 3 rows containing missing values (stat_smooth).       2) The results are again surprising, to a point where I thing I might be doing something wrong. I wanted to attach one of my plots, but I lack the reputation points to attach images. In any case, the faster condition is _always_ faster, from the first trial where they appear. This doesn't make any psychological sense, and so I am wondering whether this smoothing technique is the appropriate one or not. Thanks in advance.


title: How to use stochastic variable without linking to data in PyMC?

text: When creating some observed data, I want to use a stochastic parameter but I do not want the stochastic parameter to be affected by this data - I don't want its likelihood to be affected by this data set. Is there a way to do that (I apologise if I have made any errors of terminology!) To be more concrete I am trying to do something like the following:               nu = mc.Gamma('nu', 3, 2)     # want to have nu fit to this data set     fit_obs = mc.NoncentralT('fit', fitmean, fitprec, nu, value=fit_data, observed=True)     # want to use nu in this model but do not want this data to      # effect the likelihood of nu     target_obs = mc.NoncentralT('target', targetmean, targetprec, nu, value=target_data, observed=True)      I understand in the above as it is the conditional likelihood of nu will be affected by both data sets (so it will be fit to both of them). EDIT: I have tried the following:               @mc.deterministic     def fixed_nu(nu=nu):         return nu      but as shown in Cam's answer this doesn't work.


title: How does Soboljansen deal with categorical input against continuous output?

text: I have two data.tables which correspond to each other:               ACTid<-read.table("ACTid.txt") %# Categorical input factor (1x1000)     Y<-read.table("Y.txt")     %# Continuous output (1x1000)     >ACTid     ...     992   1     993   1     994   1     995   6     996   3     997   6     998   3     999   1     1000  1               >Y     ...     992  1.074105e+01     993  8.108669e+00     994  6.361371e-01     995  2.352872e+00     996  7.873868e-02     997  1.582426e+02     998  1.435530e+01     999  1.913109e+01     1000 8.050628e+00          x <- soboljansen(model = sobol.fun, ACTid, Y, nboot = 100)     print(x)     plot(x)      How does `soboljansen` know what the relationship is between ACTid and cfu? And also how does it deal with categorical data? **UPDATE** My model is of following multiplicative form: > Let us simplify notation by allowing the number of germs (cfu) units picked > up by nurses' (i=1:1000) hands to be called $Y$. > > \begin{equation} Y_{i}=\displaystyle \lambda_i V_i A_i-\beta_iY_{i-1} > \end{equation} > > Where values are sampled via Monte-Carlo sampling for each run for: > $h=lognormal~(1.5,0.1)$, $\lambda=\Gamma(15,3)$, $\beta=$empirical non- > negative. $A=lognormal~(7,1.9)$ and $V=$empirical non-negative. The choice > of these parameters is further dependent on the type of care $ACTid=1:6$; sobol.fun= today is my first day with R. Matlab I'm competent with. Do you know of an example I can work through?


title: Moderation in repeated-measures design?

text: **Context** : Both dependent $(Y_1,~Y_2,~Y_3)$ and independent $(X_1,~X_2,~X_3)$ variables were measured repeatedly at three time points, $\text{Time}_1$, $\text{Time}_2$, and $\text{Time}_3$. Moderator $M$ is a continuous variable measured on one occasion, and is hypothesized to be unchanged throughout all three time points. **Preliminary analyses** : I ran moderated regression cross-sectionally for each time point, and interaction between $X$ and $Z$ is only significant at $\text{Time}_2$. This is because theoretically $\text{Time}_2$ has a unique characteristic that makes it different from $\text{Time}_1$ and $\text{Time}_3$. We believe that this characteristic made a difference in whether $Z$ played a moderating role in the relationship between $X$ and $Y$. **Question** : It appears that the relationship between $X$ and $Y$ not only depends on $M$, but also the measurement occasion/condition. I'm guessing there might be another level of interaction, but how do I test that when this additional interaction is occurring on a repeated-measures design?


title: How to deal with a fake scientist?

text: What makes you to believe that statistics is something is worth to think about. Isn't that a fake science? Do you think that statistics is a real science? Do you think statisticians should be paid?


title: Correct use of cross validation in LibsSVM

text: I am classifying data points from two different groups using `LibSVM`. I do the following:   1. Creating the input file for `LibSVM`. In the input file, I put all the data I have.   2. Scaling it (using `svm-scale`).   3. Using `grid.py` of `libSVM` for choosing `gamma` and `c` parameters.    4. Using `svm-train` for training. I use the entire dataset. I also use the `-v 10` option for 10-fold cross validation (`svm-train flag`). **My questions are** : a. Is the `-v 10` option of cross validation can replace the testing step? b. The result given by the steps above is suspiciously high (96%), and so I'm wondering if I am doing something wrong? c. Could the use of `grid.py` for parameter selection before the train + cross validation damage the results (as if I were testing on data I've already trained)?


title: Cleaning corpus

text: I have a classifier that I want to use on short documents and I need to improve quality of the corpus used for training. I have 250,000 documents and manual inspection is expensive. The collection of short documents has been pre-labelled by a collection of rules as good or bad. The accuracy of this pre-processing step has been assessed by sampling to be around 98% or higher. When I train the classifier on 80% of the corpus and test its accuracy on the remaining 20%, it has a 98% accuracy. What I am wondering is this: if I filter the corpus and remove the documents for which there is disagreement and retrain and repeat, can I expect to arrive at a corpus for which the classifier agrees with the original labeling process within some small margin of error (less than 1%)? EDIT: After finding a corpus that appears to be clean, I would manually inspect the documents that were rejected. I expect there would be around 5000 (2% of 250,000) Thanks!


title: How to translate a histogram to natural language?

text: Assume we have a set of histograms, let's say they describe age distribution of some people. We want to translate this:                 x       x x     x x x x x x      to: _mostly young_ ; and this:                     x         x x x     x x x x x x      to: _mostly middle-age_ ; and this:                 x       x x   x     x x x x x x      to, probably: _mostly young with some old_. There could be more complicated cases. Which fields can help with this? Fuzzy logic sounds like to be helpful, but I cannot see how it can help with distributions.


title: Bonferroni correction in generalised linear model

text: I am running a generalised linear model in R. I have a single response variable and a maximum of 4 possible explanatory variables. I am adding each explanatory variable to the model sequentially, based on whether the coefficient is statistically significant. If the coefficient for an explanatory variable is statistically significant at 0.05, the explanatory variable remains in the model. If the coefficient for an explanatory variable is NOT statistically significant at 0.05, the explanatory variable is removed from the model. I am wondering if instead of using 0.05, I should be using a Bonferroni corrected P value? Should I use a threshold of 0.05/4 = 0.0125?


title: OLS result with LASSO

text: I'm using the lasso function in MATLAB to perform regularized regression for parameter selection. I'm new to LASSO, so I generated some toy data and tested LASSO with lambda = 0 to make sure I would get the OLS result (since there would be no penalty for model size). Unfortunately, I'm getting a value pretty close to the OLS result, but not equal to it, which puzzles me. Here is my code:               x0 = linspace(1, 5, 100)';     X = [x0, x0.^2, log(x0)];     B = [1 3 5]';     Y = X*B;     [B_lasso, FitInfo] = lasso(X, Y, 'alpha', 1, 'Lambda', 0, 'Standardize', 1)     B_ols = (X'*X)^-1*(X'*Y)      Here, B_ols = [1, 3, 5]' as expected but B_lasso = [1.14, 2.99, 4.83]'. I'd really appreciate some input as to why B_ols and B_lasso both do not equal B in the case where lambda is 0 and there is no noise added to X*B when generating Y. Thank you!


title: Truncating error bars in publication graph?

text: I am preparing a graph for publication: it has 3 panels, and two groups (line graphs) in each panel, with error bars at each time point. For one of the panels, the last time point (with fewer observations) has quite wide error bars. **Is there a hard and fast rule about the range of the y-axis?** Anything to cite to say one way or the other? I can extend the range of the y axis to make sure the upper limit is included, but then it squashes the bulk of the data points together and doesn't really display the main trend of the bulk of the graph well. It's not about manipulating peoples views; it's about accurately portraying the (vast majority of the) section of interest as opposed to ensuring all points are included.


title: Transforming two sets of test scores to be "comparable"?

text: Given 2 sets of test scores (the same students take both tests), I'd like to transform them such that they can be added to give a combined score. The 2 tests have different variances, so adding them directly wouldn't work. Dividing each score by the standard deviation of that test then adding is better but this results in transformed tests that have different values for full marks. For example, if both tests were out of 50, then after transformation one is out of 40 and the other out of 100, and this is found disconcerting as it doesn't feel like we're weighting the tests equally. I'm wondering if there is a way to make both tests have the same variance while the score representing full marks is still equal in the two tests? For example, one way might be to convert both to percentiles and then add the percentiles, but I'm not sure if would be considered reasonable/acceptable practice, or if there is a better (e.g. nonlinear) way? Thanks a lot!


title: Understanding MCMC and the Metropolis-Hastings algorithm

text: Over the past few days I have been trying to understand how Markov Chain Monte Carlo (MCMC) works. In particular I have been trying to understand and implement the Metropolis-Hastings algorithm. So far I think I have an overall understanding of the algorithm but there are a couple of things that are not clear to me yet. I want to use MCMC to fit some models to data. Because of this I will describe my understanding of Metropolis-Hastings algorithm for fitting a straight line $f(x)=ax$ to some observed data $D$: 1) Make an initial guess for $a$. Set this $a$ as our current $a$ ($a_0$). Also add $a$ at the end of the Markov Chain ($C$). 2) Repeat the steps bellow a number of times. 3) Evaluate current likelihood (${\cal L_0}$) given $a_0$ and $D$. 4) Propose a new $a$ ($a_1$) by sampling from a normal distribution with $\mu=a_0$ and $\sigma=stepsize$. For now, $stepsize$ is constant. 5) Evaluate new likelihood (${\cal L_1}$) given $a_1$ and $D$. 6) If ${\cal L_1}$ is bigger than ${\cal L_0}$, accept $a_1$ as the new $a_0$, append it at the end of $C$ and go to step 2. 7) If ${\cal L_1}$ is smaller than ${\cal L_0}$ generate a number ($U$) in range [0,1] from a uniform distribution 8) If $U$ is smaller than the difference between the two likelihoods (${\cal L_1}$ - ${\cal L_0}$), accept $a_1$ as the new $a_0$, append it at the end of $C$ and go to step 2. 9) If $U$ is bigger than the difference between the two likelihoods (${\cal L_1}$ - ${\cal L_0}$), append the $a_0$ at the end of $C$, keep using the same $a_0$, go to step 2. 10) End of Repeat. 11) Remove some elements from the start of $C$ (burn-in phase). 12) Now take the average of the values in $C$. This average is the estimated $a$. Now I have some questions regarding the above steps:   * How do I construct the likelihood function for $f(x)=ax$ but also for any arbitrary function?   * Is this a correct implementation of Metropolis-Hastings algorithm?   * How the selection of the random number generation method at Step 7 can change the results?   * How is this algorithm going to change if I have multiple model parameters? For example, if I had the model $f(x)=ax+b$. Notes/Credits: The main structure of the algorithm described above is based on the code from an MPIA Python Workshop.


title: Confidence bounds for PDF

text: I build confidence bounds for estimating PDF of the empirical sample using bootstrapping:               data <- rnorm(1000)     d <- density(data)          boot <- replicate(100,                         { x <- sample(data, replace=TRUE);                           density(x, from=min(d$x), to=max(d$x))$y})           CI <- apply(boot, 1, quantile, c(0.025,0.975) )                hist(data, freq=F, ylim=c(0,max(CI[2,])))     lines(d, lwd = 2)     lines(d$x, CI[1,], lty=2)     lines(d$x, CI[2,], lty=2)      I have questions:   1. From what reasons I need to choose the number bootstrapping repeats?   2. How can I use confidence bounds for the determination of the minimum required sample size?


title: 3D space learning and prediction

text: I want suggestions about learning and predicting some object's position before hitting one out of four sides of a wall. I have some priority according to side of wall, and of course all the scenarios would be in 3D space having `x,y,z` coordinates of the object's trajectory coordinates. The simple scenario is shown in this picture.   ![Walls](http://i.stack.imgur.com/NIPzs.jpg) I want suggestions as I have data in 3D space having `x,y,z` coordinate sets, and I want to train and predict that four bouncing balls hit which side of the wall and which ball first in the case of many-to-one hitting. I have bouncing balls `set = {B1,B2,B3,B4}`: _Red, Green, Blue,_ and _Yellow_ respectively, and four sides of a wall having some priority. Which learning machine should I go for, and what optimization algorithm is best for my scenario, and on what basis?


title: Weighted SVM and Up-sampling

text: I'm performing classification with the libSVM package in R and am wondering about the correct procedure for weighting or up-sampling. I have a data set that is 19,396 observations of which only 81 are the positive class. Since I'm dealing with such an unbalanced data set I would like to use weighting or up- sampling to ensure the positive class is adequately represented in model building. In utilizing weighted SVM I'm just trying to figure out what the correct weight for the positive class would be. I took 2/3 of my data as a training set and ended up with the following number from each class:                  Y     N         53 12878      Should I then use 12878/52 = 247.6538 as the weight, or 19396/53, or 19396/81? Is there a rule of thumb for this? Also, if using up-sampling should I up-sample my training set only until it is balanced, resulting in a training set with 25,756 records and leave the test set alone? Thanks.


title: Explanation of order of integration in the context of cointegration

text: Please explain the meaning of 'order of integration' when talking about cointegration.   An explanation with some examples would be great.


title: Issues with solving large sparse linear equations

text: I have some issues solving sparse linear equations $Ax = b$ My matrix $A$ is sparse with dimension of 5 million by 5 million. Actually, it is a combination of two matrices. One is tridiagonal and the other matrix has a bigger band. For example the $29450^{th}$ row has non zero diagonal entry. Plus non zero element at $14725^{th}$, $44175^{th}$, $29451^{th}$, $29449^{th}$,$29605^{th}$,$29295^{th}$,column. So every row has almost 6 non zero elements plus non zero diagonal entry. I tried to use a couple of iterative solvers in matlab, but they don't seem to give good accuracy. Also if I run the iteration for a long time, I will have round off errors and all. What can I do to get the best estimate?


title: Correct approach for performing retrospective A/B testing on existing data

text: I have large amounts of data about events and historical behavior (clicks, purchases etc) of users on a site. I would like to be able to correlate the inputs (like clicks, purchases, other events) with the outputs of interest (lifetime order value of customer etc). Can I use the exact same math and process that I would use for a typical A/B test, where I don't have the data ahead of time? Or, is there a different approach that should be taken in this case?


title: Test to determine preferential distribution of samples in two groups

text: I have data structured hierarchically in this way, where single observations ("Member") are grouped in larger groups:               Major_group    Member       Sample1    Sample2   Sample3     GroupA         Member1      5          0         2     GroupA         Member2      3          4         10     GroupA         Member3      2          0         0     GroupB         Member4      10         3         7     GroupB         Member5      23         5         9      And so on. For each sample, the observations are counts of specific features of Member1,2... etc. If no specific features are present, the value is 0. I have two sets of observations between two different groups of samples (and groups have unequal size). My hypothesis is to determine if, at the `Major_group` level, there is a **preferential distribution** of non-0 counts for that major group in either of the two groups. I thought of using the Fisher's Exact Test for this, by building matrices like (for example, for groupA; numbers are only examples):                           Sample_group1     Sample_group2     Affected    10                4     Unaffected  9                 5      Where for each group `Affected` is the number of non-0 counting samples (and `Unaffected` the number of samples with observations = 0). My doubts of on this approach:   * is it "flattening" the contributions of individual members?   * is it actually correct for what I'm trying to do? As a disclaimer: my statistics background can be lacking, I'm a biologist by background.


title: Ordinal/continuous vs dummy variable for time series regression/data mining

text: Let's suppose I have a time series data that I would like to regress $y$ on $x$ and $Time$. See below for the dataset.               y   x   time     12  100 1     14  101 2     16  102 3     18  103 4     20  201 1     22  202 2     24  203 3     26  204 4      **Approach 1:** One approach is to do a multiple linear regression or neural networks or SVM directly on the dataset above by treating time as an ordinal (?) or continuous variable. I can do a time series regression as follows: $$y = a+\beta_1x+\beta_2 time+\varepsilon$$ where $\varepsilon$ is modeled as ARMA **Approach 2:** Alternatively I can rearrange data and create dummy variables for time as follows:               y   x   Time_1  Time_2  Time_3     12  100 1   0   0     14  101 0   1   0     16  102 0   0   1     18  103 0   0   0     20  201 1   0   0     22  202 0   1   0     24  203 0   0   1     26  204 0   0   0      and do a time series regression/neural network/SVM. For instance a time series would be $$y = a+\beta_1x+\beta_2 time_1+\beta_3 time_2+\beta_4 time_3+\varepsilon$$ where $\varepsilon$ is modeled as ARMA. Below are my questions:   1. What is the right approach - 1 or 2?    2. If we use the data mining approach such as a neural network or svm does it matter if we use either approach 1 or 2 ?   3. Is approach 1 more parsimonious since we have 1 variable representing time as opposed to approach 2 which has 3 variables?


title: Penalized cross validation in R

text: I have generated a L1 penalized Cox model using the `penalized` package in R. I used the `optL1()` function to generate the Breslow object (see below):               fit <- optL1(surv.obj, penalized = ..., etc)      In the reference manual, it says the `fit$predictions` are the cross-validated predictions for the left out samples. Does this mean these are predictions are based off the fold that was left out of the model building (what one usually thinks of as cross-validated) or is this just saying that these predictions are based off the optimal lambda which is established by cross-validation and are thus just refitting the data? I noticed that `predict(fit$fullfit)` and `fit$predictions` give two different values if that helps.


title: How to test if the mean of data collected over many days is significantly higher on a predicted day

text: I have some data examining blog posts on different days. Basically, about 2000 news articles pertaining a certain topic were sampled and each blog post was given a positivity percentage score depending on the amount of positive words in the text (from 0 to 1). This data was collected over a week. Therefore, there are 7 separate means, one for each day representing the average positivity of posts for a topic. I also have a sample of non-topic related news articles over a week that have been coded for positivity as well. I have a specific hypothesis that positivity in the articles about the topic should have increased on a specific day (e.g. on the 25th). However, this increase should be independent of any sort of non-topic related positivity increase that occurred on that day (i.e. it's not just because all articles are generally happier on that day, it's something specific to the topic). Is there a way to test that hypothesis with the data collected?


title: Drawing an ROC curve from one point

text: First, thanks for the wonderful site. Second, I am trying to use an ROC curves to check the reliability of some nondestructive testing of concrete slabs. As you may know, the Non-destructive tests(NDT) work same as any diagnostic system where the output either yes or no with some necessary transformation. I am planning to use some NDT technologies like (Impact echo, Ground penetration radar, infrared thermography, half cell potential, and hammer sounding). The problem is when I am trying to calculate the true positive rate (TPR) as the sensitivity and false positive rate (FPR) as 1- specificity, I will end up with a single point in the ROC curve. Can I depend on a single point to draw the ROC curve for this case.


title: Sample size for logistic regression?

text: I want to make a logistic model from my survey data. It is a small survey of four residential colonies in which only 154 respondent was interviewed. My dependent variable is "satisfactory transition to work". I found that of the 154 respondent 73 said that they have satisfactorily transited to work, while the rest did not. So the dependent variable is binary in nature and I decided to use logistic regression. I have seven independent variables (three continuous and four nominal). One guideline suggest that there should be 10 cases for each predictor/independent variable (Agresti, 2007). Based on this guideline I feel that it is ok to run logistic regression. Am I right? If not please let me know how to decide the number of independent variables?


title: Exercise 2.2 of The Elements of Statistical Learning

text: The text book first generates some 2-class data via: ![enter image description here](http://i.stack.imgur.com/rf6xF.png) ![enter image description here](http://i.stack.imgur.com/tXrP7.png) which gives: ![enter image description here](http://i.stack.imgur.com/OzIxZ.png) and then it asks: ![enter image description here](http://i.stack.imgur.com/EdClh.png) I try to solve this by first modelling this with this graphical model: ![enter image description here](http://i.stack.imgur.com/ZS2pv.png) where $c$ is the label, $h\,(1\le h \le 10)$ is the index of the selected mean $m_h^c$, and $x$ is the data point. This will give $$ \begin{align*} \Pr(x\mid m_h^c) =& \mathcal{N}(m_h^c,\mathbf{I}/5)\\ \Pr(m_h^c\mid h,c=\mathrm{blue}) =& \mathcal{N}((1,0)^T,\mathbf{I})\\ \Pr(m_h^c\mid h,c=\mathrm{orange}) =& \mathcal{N}((0,1)^T,\mathbf{I})\\ \Pr(h) =& \frac{1}{10}\\ \Pr(c) =& \frac{1}{2} \end{align*} $$ On the other hand, the boundary is $\{x:\Pr(c=\mathrm{blue}\mid x)=\Pr(c=\mathrm{orange}\mid x)\}$. With Bayesian rule, we have $$ \begin{align*} \Pr(c\mid x) =& \frac{\Pr(x\mid c)\Pr(c)}{\sum_c\Pr(x\mid c)\Pr(c)}\\ \Pr(x\mid c) =& \sum_h\int_{m_h^c}\Pr(h)\Pr(m_h^c\mid h,c)\Pr(x\mid m_h^c) \end{align*} $$ But later I found that the problem setting is symmetric so this may yield $x=y$ as the boundary. If the problem is asking the boundary when $m_h^c$ are conditioned, the equation will include $40$ parameters which I think it's unlikely to be the purpose of the exercise. So am I misunderstanding anything? Thank you.


title: Univariate non-significant results becoming significant on multivariate analysis

text: I'm having a look at x-ray waiting times in an emergency department and how these relate to the day of the week and the time of day that requests are ordered. (want to know if patients have to wait longer for their x-rays on weekends and at night compared with day shifts). The waiting times have a significant right skew when plotted on a histogram. I divided the continuous variable "waiting times" into 5 categories (0-30 mins being the reference category). I also dummy coded the 6 days of the week and the "time of day" (which was divided into the 3 shifts of the day) into 0-the day in question and 1-the rest of the days of the week. I performed univariate multinomial regression with waiting times as the dependent variable and each independent dummy variable in turn (omitting Saturday and Night shift from the analysis to act as "reference"). when I entered each variable independently, not one of them showed a significant relationship with x-ray waiting times. On multivariable multinomial regression, I found four variables were demonstrating p values less than 0.05, particularly with the waiting time category of 61-90 mins. I wondered if anyone knew why this was? if there's a plausible explanation or whether i've made a mistake somewhere? Please help if you have any suggestions as I'm a complete stats amateur....... Thanks very much.


title: Does the Holt-Winters algorithm for exponential smoothing in time series modelling require the normality assumption in residuals?

text: I'm working on a project to compare different approaches to time series modeling. In the model selection process, we perform residual analysis for the fitted models. For regression, we need to check the normality assumption in the residuals by performing Shapiro-Walk's test, or plotting QQ-plot, etc. Do we also need to check whether the residuals are normally distributed after fitting a model using HW algorithm?


title: Non-Bayesian alternatives to maximum likelihood estimators and method of moment estimators when there's only one observation

text: When trying to estimate the parameters of a known distribution, it might occur that the maximum likelihood estimator and the method of moment estimator don't work well when there's only one observation. What are some good non-Bayesian alternatives to maximum likelihood estimators and method of moment estimators if there's only one observation?


title: What does bandwidth mean?

text: I plotted density function in R and under the plot is number of bandwith. What does this number mean?


title: Multiple linear regression for classification

text: Recently I have been playing with linear regression for classification. For the binary case, it is straightforward. For the multiple case, I know there are issues with using this approach (like masking). However, I cannot fully understand how this masking arises. As understand in the multiple case one has to compute a discriminant function for each class and then classify according to the largest value obtained from each function. So why does MLR have issues with masking?


title: RSS examinations

text: Are the RSS examinations worth doing as a way to learn more statistics? I have been working a physicist for most of my life so I'm a little bit loathed to pay to learn maths as I can usually get my head around new statistical methods as and when I need them. The problem I have is that I only use what I know about and on more than one occasion after lurking on sites like this I have found myself thinking 'Oh, if only I'd have known that back then'. The RSS course seems like a way to force myself to learn and only have to pay for the examinations.


title: propagating errors through an ANOVA?

text: I'm sorry if this question has been answered before. I have a dataset:                 -1.5   -1.8    -1.625       -1.25  -1.525  -1.675       -1.15  -1.35   -1.65      etc Each column is a treatment, each number is a sample. I want to divide each treatment by another number of similar magnitude and units, but this number has a standard deviation that was calculated (from confidence intervals). There's no variation in that dataset because of the nature of the measurement. I'd like to run an ANOVA and maybe some pairwise comparisons on the resulting data (after subtraction), but I'm not sure what to do with the error term in the calculation parameter. Thanks! To be more specific:                  tmt 1    tmt 2   tmt 3        -1.5     -1.8    -1.625        -1.25    -1.525  -1.675        -1.15    -1.35   -1.65        -1.1     -1.45   -0.65        -0.9     -0.9    -0.9        -0.9     -1.15   -0.85        -1.35    -1.6    -1.6        -0.9     -1.8    -1.8        -1.35    -1.4    -1.5      This is the original dataset. Each column is a treatment, and each row within a column is an individual sample. I have an additional parameter for teach treatment (tmt). This additional parameter has a standard deviation associated with each treatment:                  tmt 1    tmt 2   tmt 3        -1.1     -1.5    -2.5             stdv 1   stdv 2  stdv 3        0.051    0.153   0.087      This is in the same units as the above dataset. I subtract the additional parameter from the original dataset to get a second dataset. The parameter (small table above) that was associated with each treatment (column) was subtracted from every sample in the original treatment:                  tmt 1   tmt 2   tmt 3        -0.4     -0.3    0.875        -0.15    -0.025  0.825        -0.05    0.15    0.85        0        0.05    1.85        0.2      0.6     1.6        0.2      0.35    1.65        -0.25    -0.1    0.9        0.2      -0.3    0.7        -0.25     0.1    1      This second dataset is what I want to run an ANOVA on, to look at differences in this calculated parameter between species. I'm not sure how (or if its even required) to take the error term from the additional parameter into account. Sorry about the nebulous question - I'm pretty new to all this


title: Extended COX model in R

text: Objective: To predict time to event of a customer using R I developed a Non parametric model but as per my understanding we cannot perform predictions from these models, so I have to build a model which can predict for any test data. I hope extended cox model will work for my case (we have time depended co variants) Questions:   1. How the data structure would be for extended cox model?   2. What we will predict from Cox model (In logistic model we will predict probability in the similar lines what is for survival model)?   3. Variable Importance   4. Can we predict time when a customer will leave from the study?   5. How to validate the model? It would be great if someone provide the R code for the above questions 2, 3, and 4 Thanks in advance.


title: Multi-class logarithmic loss function per class

text: In a multi-classification problem, we define the _logarithmic loss function_ $F$ in terms of the _logarithmic loss function per label_ $F_i$ as: $$ F = -\frac{1}{N}\sum_{i}^{N}\sum_{j}^{M}y_{ij} \cdot Ln(p_{ij}))=\sum_{j}^{M} \left (-\frac{1}{N}\sum_{i}^{N}y_{ij} \cdot Ln(p_{ij})) \right ) = \sum_{j}^{M}F_i $$ where $N$ is the number of instances, $M$ is the number of different labels, $y_{ij}$ is the binary variable with the expected labels and $p_{ij}$ is the classificiation probability output by the classifier for the $i$-instance and the $j$-label. The cost function $F$ measures the distance between two probability distributions, i.e. how similar is the distribution of actual labels and classifier probabilities. Hence, values close to zero are preferred. However, the cost function per label $F_i$ has any meaning? It seems that is measuring how good our classifier is doing per label, but it is affected by the number of instances $N$ that don't contain this label.


title: Probability of two aces given at least one ace (out of two cards) vs probability of two aces given one card is the ace of spades (out of two cards)

text: So, in the video lectures from Harvard's Statistics 110:Probability course that can be found on iTunes and Youtube, I encountered this problem. https://www.youtube.com/watch?v=JzDvVgNDxo8#t=566 I tried to summaries it here below: Suppose we are given a random two-card hand from a standard deck.   1. What is the probability that both cards are aces given that we have at least one ace? $$ P(both\ aces | have\ ace) = \frac{P(both\ aces, have\ ace)}{P(have\ ace)} $$ Since having at least one ace is implied if you have both aces, the intersection can be reduced to just $P(both aces)$ $$ P(both\ aces | have\ ace) = \frac{P(both\ aces)}{P(have\ ace)} $$ This is then just $$ P(both\ aces | have\ ace) = \frac{4C2\ /\ 52C2}{1-48C2\ /\ 52C2}=\frac{1}{33} $$   1. What is the probability that both cards are aces given that we have the ace of spades? $$ P(both\ aces | have\ ace\ of\ spades)=\frac{P(both\ aces, have\ ace\ of\ spades)}{P(have\ ace\ of\ spades)} $$ $$ P(both\ aces | have\ ace\ of\ spades)=\frac{(3C1*1C1)\ /\ 52C2}{2!*\frac{51}{52}*\frac{1}{51}}=\frac{1}{17} $$ Now, somewhere along these examples I got lost... The latter is obviously just the same as $\frac{3}{51}$, which makes a lot of sense (to me) that this would be the answer. If you are told that you have the ace of (say) spades, then you know that there are 3 more aces and 51 more cards. But in the former example, the math seem fine (and I believe the lecturer wouldn't give this example if it was incorrect...), but I can't wrap my head around this. How do I get some intuition for this problem?


title: Supervised Pattern Recognition with Probabilistic Labels

text: I am interested in supervised pattern recognition problems where the the label associated with each pattern gives the probability of membership for each of the $c$ classes, rather than assigning each pattern unequivocally to a single class. This can arise because the (possibly human) oracle providing the labels is unable to decide which class the pattern truly belongs to, or the labels are provided by another probabilistic classifier, and we don't want to throw away information by assigning each pattern to its most probable class. Note this is not a label noise or semi-supervised learning problem; the probabilistic labels accurately represent the oracles' state of knowledge. I've been experimenting with this type of model for sometime, but I can't find any papers where this type of model s explicitly discussed, so I expect I am not searching with the right terms. Can anybody suggest any suitable books/papers?


title: Appropriate PCA/ EFA rotation method

text: Question about appropriate PCA/ EFA rotation method. While I see that oblique rotation methods (e.g. promax) are suggested for correlated data (and I have correlated data) I see that the majority of researchers use Varimax rotation in PCA/ EFA. I intend to use factor scores in cluster analysis, where it is important that the considered variables are not correlated, therefore can orthogonal rotation (which I assume would impose greater distances) be justified on data that exhibit some correlation?


title: Questions about number of groups and group size in planned comparisons in ANCOVA

text: I have three questions.   1. I want to conduct a planned comparison in ANCOVA which compares three groups (n=2, n=7 and n=5) with two groups (n=11 and n=25). Is it ok to include the group with n=2 if I am only conducting a planned comparison as described above?   2. Is it ok to conduct an ANCOVA when I am only comparing two groups? E.g. if I wanted to compare left and right hemisphere groups, controlling for another variable.   3. If the DV in an ANCOVA has been transformed (because it was non-normal) does the covariate have to undergo the same transformation? E.g. If the DV is (log transformed) heart rate during an experimental condition, and the covariate is heart rate measured at baseline (which is normally distributed), does the baseline HR need to be transformed as well?


title: Percentile rank transforming z-scores

text: I have been looking through methods to convert z-scores to positive values, without taking absolute values and somehow taking the sign of the score into account (I avoided taking absolute value, because I want to differentiate between transformed value for z-scores with same magnitude but different sign). The z-score matrix I have is actually an z-score difference matrix calculated from two Pearson correlation coefficient matrices converted to z-scores via Fisher-r-to-z transformation. After searching around, I found out about percentile rank transformation of z-scores, but I do not have knowledge to determine whether this transformation can be applied in all the cases (especially the case I have). So my question is it this percentile rank transformation applicable in my case and are there any other z-score transformations that I could look into ? Thank you.


title: Interarrival time distribution of uniform arrival process

text: I am currently trying to model an uniform arrival process within my simulation model. However, I can only model it by means of an interarrival time (I can let the model wait for a certain amount seconds). I want to have: $$G(x) = {1 \over F(x)}$$ where $ F(x) \sim U(a,b) $. I started by using inverse transform sampling on the CDF of the inverse continious distribution (I should use a discrete distribution, but since the numbers are generally quite big, this is acceptable): $$ G(y) = { b - y^{-1} \over b - a }$$ $$ x = { 1 \over b - y(b-a) } $$ where $Y \sim U(0,1)$. However this gives me the wrong results. Lets take some number to clarify this. Let say I want to model an arrival process between 600 and 1800 units per hour. On average (1200/hr), the model has to wait 3 seconds. The boundaries of the waiting time are between 2 and 6 seconds. This means that 50% of the samples should lie in the interval $[2,3]$ and 50% in the interval $[3,6]$. This is exactly what $G(y)$ does, so I was expecting good results. However, if I calculate the average from this function, it does not match with what I want. $$g(y) = y^-2 { 1 \over b^* - a^*}$$ $$ E\left[G\right] = \int_{b^*}^{a^*} g(y) y dy = { b^* -a^* \over \ln \left| b^* \right| - \ln \left| a^* \right| } $$ where $a^* = {1 \over a}$ and $b^* = {1 \over b}$. When I plug in the numbers I get an average of 3.64, while I would like to have an average of 3. Do I make a mistake somewhere? Or, does somesome have a better approach in generating a random arrival process?


title: Explain the difference between type 1,2,3,4 sums of squares

text: Can someone explain the difference between type 1:4 sums of squares using examples with relatable variables?


title: If an ANOVA indicates no main effect and no interaction, should the lack of interaction be stated?

text: I measured response variable $Y$ at three levels of factor $A$ and four levels of factor $B$, $n=6$ reps/treatment. Results include   * $A$ has strong effects on $Y$.   * $B$ has no effect on $Y$   * There is no $A*B$ interaction I would like to report all three results (the second two are actually more interesting than the first since all were expected). So far I have: "The effect of A on Y was significant (ANOVA, $P<0.001$) and $Y$ was different at each level of $A$ (Tukey HSD, $P<0.001$). There was no effect of $B$ on $Y$, and there was no interaction between $A$ and $B.$ This implies that $B$ had no effect on $Y$ **and that the effect of $B$ on $Y$ was similar for all levels of $A$.** " (Text with bold emphasis is in question) This approach minimizes the double negative as in "the effect of $B$ on $Y$ was not different at different levels of $A$". Still, neither sounds 'right' to me, thus my post. In my field, as in many fields of science, non-significant responses are rarely mentioned in the text, and even more rarely are they interpreted, but in this case, previous studies find that $B$ almost always affects $Y$ and the lack of an interaction is notable.


title: Related sample t-test with missing values

text: I want to run a t-test on two conditions tested repeatedly on the same population. I'm doing this via a python implementation of a related sample t-test (`scipy.stats.mstats.ttest_rel`). The related sample t test requires my input arrays to have the same length. I do however have a number of missed measurements (which have a measurement value of "false", which python reads as zero and which distorts my t-test). What would be the less disruptive way to deal with missed measurements:   * replace each missed measurement with the **mean** of the rest of the data   * replace each missed measurement with the **median** of the rest of the data   * delete for each missed measurement the valid measurement from the other condition ?


title: Can this be modeled by a two-sample t-test?

text: Here is the question: A claim is made that men from two income families earn less than men who are the sole income earners in their household: To support this claim, data from a random sample of 22 married men was gathered. For this sample, the mean and standard deviation of salaries for 11 married men in two-income families was \$95,000 and \$15,000 respectively. The mean and standard deviation for the remaining 11 married men who were the sole source of family income was \$124,000 and \$18,000 respectively. (a) What are appropriate null and alternative hypotheses to test this claim? (b) State an appropriate test statistic to test this hypothesis and the null distribution of this test statistic, stating any assumptions that you use. (c) Calculate the test statistic chosen in (b), the corresponding p-value and form an appropriate conclusion. Here is my answer: (a) Essentially the null hypothesis is that there is no difference between married men from two income households versus one. So, I have: $H_{0} : U_{x}=U_{y}$, the alternative is that men who are sole earners do earn more, $H_{1} : U_{x}<U_{y}$ (b) Since the pop. variance is unknown, I use a pooled variance - a test statistic following a t-distribution is used : $\tau = \frac{\bar{X}-\bar{Y}}{S_{p}\sqrt{\frac{1}{n_{x}}+\frac{1}{n_{y}}}} \tilde{} t(n_{x}+n_{y}-2) = 20 $ My assumptions are that the samples have been obtained from a normal distribution population and that the populations have equal variance. (c) For the men from two income households we have: $\bar{x} = \$95,000 $ and $S_{x} = \$15,000 $, and the sole earners: $\bar{y} = \$124,000 $ and $S_{y} = \$18,000 $ The pooled variance is: $s_{p}^2 =\frac{(n_{1} - 1)s_{1}^2 + (n_{2} - 1)s_{2}^2}{n_{1}+n_{2} -2} $ $ = \frac{(11 - 1)15,000^2 + (11 - 1)18,000^2}{11+11 -2} = \$274,500,000$ so $ s_{p}= \$16,584.04$ Therefore the observed test statistic is: $\tau_{obs} = \frac{95,000-124,000}{16,568.04 \cdot \sqrt{\frac{2}{11}}} \approx -4.10$ With a 5% level of significance, the critical value is: $t_{0.05,20} = -1.725$ Since $-4.10 < -1.725$ we reject the null hypothesis. The alternative hypothesis is $H_{1}: U_{x}-U_{y} < 0$, so the p-value is $P(t_{20} \le \tau_{obs}) = P(t_{20} \le -4.1) \approx 0.00028$ Since the p value is less than the level of significance we can conclude that we reject the null hypothesis at the 5% level of significance. I.e. the average salaries of men who are sole earners are higher than those men in two income homes. Have I done this correctly? Also, what exactly is the 'null distribution' of this test statistic? I read somewhere that it is the distribution when the null hypothesis is true. So is that just the assumption that the samples are from a normal distribution? I don't quite understand this concept. Many thanks for the feedback in advance!


title: Phase space distribution of heads and tails (coin toss)

text: So I have the equation $$h(t) = 1 + vt - \frac{1}2 gt^2 \pm sin(\omega t) $$ to describe the motion of a flipped coin. It is just a kinematics equation with an angular component added to it, where $h$ is the height of the coin, $v$ is the upwards velocity of the coin, $\omega$ is the angular velocity, and $t$ is time. What I want to show is a 2D graph of angular velocity against velocity ($\omega$ vs. $v$) so that for all combinations of $v$, $\omega$, and $t$ such that $h = 0$ (when the coin lands), the graph is divided into regions that indicate whether the coin landed on heads or tails. What is the approach for doing this? I prefer Mathematica if possible.


title: Testing if treating a categorical variable as continuous is okay

text: Pardon me if this is a simple question, but I haven't found a great resource for this just yet. Yes, I know that when possible, we should try to treat ordinal variables as categorical rather than continuous, but I have a particular application where being able to interpret the ordinal variable (parental education) as interval would be really, really convenient. Pasta (2009) (bottom of page 5) uses the following to show that the deviations from nonlinearity of an ordinal variable are insignificant: You run models (in this case, OLS) of the outcome variable on the ordinal variable of interest. In this model, you include BOTH the version of the variable as continuous the version as categorical. The t-tests from this model on the categorical levels output from statistical software will give you a general idea if the nonlinearity is significant. I don't know how convincing that is seeing as the model drops two levels for collinearity reasons, so I was wondering if the following might work instead: You do the same, except you run separate models for each level of the categorical variables. So for example, you have the outcome variable _y_ and an ordinal independent variable _x_ , with 4 levels. As such, the categorical version of _x_ is _x.cat_ , and the 4 level dummies are _x.1_ , _x.2_ , _x.3_ , and _x.4_. You then run 4 models:   * _y ~ x.cat+x.1_   * _y ~ x.cat+x.2_   * _y ~ x.cat+x.3_   * _y ~ x.cat+x.4_ Could you then use these t-statistics on the dummies in each model to determine the extent to which nonlinearity is significant? Any feedback would be greatly appreciated.


title: Meaning of output of function "ar" in R

text: How should I read the output of the function `ar` in R. For example, take this VAR model:               library(tseries)     data(USeconomic)     US.ar <- ar(cbind(GNP, M1), method="ols",                 dmean=T, intercept=F)      (from the book `Introductory Time Series with R` by Cowpertwait) So `Us.ar` produces the following output:               > US.ar          Call:     ar(x = cbind(GNP, M1), method = "ols", dmean = T, intercept = F)          $ar     , , 1                   GNP    M1     GNP  1.27181 1.167     M1  -0.03383 1.588          , , 2                   GNP      M1     GNP -0.00423 -0.6942     M1   0.06354 -0.4839          , , 3                   GNP      M1     GNP -0.26715 -0.5103     M1  -0.02859 -0.1295               $var.pred            GNP    M1     GNP 618.69 16.38     M1   16.38 23.90      and `US.ar$ar` gives this other representation:               > US.ar$ar          , , GNP                     GNP          M1     1  1.271812104 -0.03383385     2 -0.004229937  0.06353801     3 -0.267154022 -0.02858942          , , M1                   GNP         M1     1  1.1674655  1.5876695     2 -0.6941813 -0.4838919     3 -0.5103451 -0.1294549      As I understand it, `,,1` refers to the order of the coefficient in the autoregression model. In the second representation, `, , M1` refers to the columns, whereas `GNP M1` describe rows and the numbers describe the order. I think this output is describing the following model: $$GNP_{t} = 1.27 GNP_{t-1} + 1.167 M1_{t-1} - 0.004 GNP_{t-2} - 0.6942 M1_{t-2} - 0.2671 GNP_{t-3} - 0.5104 M1_{t-3}$$ $$M1_{t} = - 0.03 GNP_{t-1} + 1.58 M1_{t-1} + 0.063 GNP_{t-2} - 0.48M1_{t-2} - 0.02 GNP_{t-3} - 0.12 M1_{t-3}$$ Therefore, in matrix notation, this should be equal to: $$ \begin{pmatrix} GNP_{t} \\ M1_{t} \\ \end{pmatrix} = \left[ \begin{pmatrix} 1.27 & 1.167 \\ -0.03 & 1.58 \\ \end{pmatrix}x + \begin{pmatrix} -0.004 & -0.6942 \\ 0.063 & -0.48 \\ \end{pmatrix}x^{2} + \begin{pmatrix} -0.267 & -0.5104 \\ -0.02 & -0.12 \\ \end{pmatrix}x^{3} \right] \begin{pmatrix} GNP_{t} \\ M1_{t} \\ \end{pmatrix}$$ However, the book says that the model is this: ![enter image description here](http://i.stack.imgur.com/lXHJA.png) As you can see, cross-terms are flipped. Is this a mistake?


title: Hypothesis testing books with problem sets and answers

text: I want to find hypothesis testing books with **problem sets** and answers covering the topics:   * Testing variance   * Type 1 and type 2 errors   * Test parameters for Poisson and binomial   * Chi-squared test I want to be able to read the question and determine on my own which test to be used.


title: Which variance inflation factor should I be using: $\text{GVIF}$ or $\text{GVIF}^{1/(2\cdot\text{df})}$?

text: I'm trying to interpret variance inflation factors using the `vif` function in the R package `car`. The function prints both a generalised $\text{VIF}$ and also $\text{GVIF}^{1/(2\cdot\text{df})}$. According to the help file, this latter value > To adjust for the dimension of the confidence ellipsoid, the function also > prints GVIF^[1/(2*df)] where df is the degrees of freedom associated with > the term. I don't understand the meaning of this explanation in the help file, so I'm not sure if I should be using $\text{GVIF}$ or $\text{GVIF}^{1/(2\cdot\text{df})}$. For my model these two values are very different (maximum $\text{GVIF}$ is ~$60$; maximum $\text{GVIF}^{1/(2\cdot\text{df})}$ is ~$3$). Could someone please explain to me which one I should be using, and what is meant by adjusting the dimension of the confidence ellipsoid?


title: unbiased estimates and MLE of central moments and of standardized moments?

text: I have heard of unbiased estimate and MLE of variance, and some about those of kurtosis. Are there general results about   * unbiased estimates of k-th order central moments?   * MLE of k-th order central moments?   * unbiased estimates of k-th order standardized moments?   * MLE of k-th order standardized moments? Thanks!


title: Can the standard deviation of non-negative data exceed the mean?

text: I have some triangulated 3D meshes. The statistics for the triangle areas are:   * Min 0.000    * Max 2341.141    * Mean 56.317    * Std dev 98.720 So, does it mean anything particularly useful about the standard deviation or suggest there are bugs in calculating it, when the figures work out like the above? The areas are certainly far from being normally distributed. And as someone mentioned in one of their responses below, the thing that really surprised me that it only took one SD from the mean for the numbers to go negative and thus out of the legal domain. Thanks


title: Statistical significance of betas in linear regression

text: This question is similar to another question I had recently posted but I have a follow on. In classical linear regression, we have $$\hat{\beta } \sim N(\beta,(X^{T}X)^{-1}\sigma^2).$$ Using this, one builds individual hypotheses of the significance of the coefficients, as done in the book by Tibshirani et al. My questions are two fold: 1) The book talks about a combined hypothesis built by proving that $$(\hat{\beta }-\beta)^T(X^{T}X)^{-1}(\hat{\beta }-\beta) \sim \sigma^2\chi_{p+1}^2.$$ I don't see how this formula can be derived from the equation I wrote above. I do see that $$\hat{\beta }-\beta \sim N(0,(X^{T}X)^{-1}\sigma^2).$$ How do we take the $X^TX$ matrix out and prove the above? I would be grateful if someone could outline the steps. 2) My second question is, how do we think about building the hypothesis? Do we think about building individual coefficient hypotheses or is it a good idea to view everything together? In other words, what is the difference/pros and cons of using the two different styles of hypotheses, the individual coefficient one or viewing everything together as per the above equation? Can we have an example of building a combined hypothesis? I am guessing that most statistical packages don't really take into account the correlation between different $\beta$ which is encoded in the matrix $X^TX$. Please clarify, any help will be much appreciated.


title: How to choose the tolerance parameter for ABC?

text: I have the following sorted data (sampling from parametric space [1,5]) with respect to their distances of parameter Theta. i.e., Let say N = 1000, Theta : 1.1, 1.7, 1.9, 2.4, 2.8, . . . , 4.9 Distance : 0.2, 0.3, 0.5, 0.9, 1.1, . . . , 1.9 From literature I know that more than 3% (i.e., 0.03*N) data are not useful, but I don't know where to cutoff? Could you suggest any re-sampling methods? How can I treat this problem? Classification or Regression? This data is basically output of rejection sampling (see http://en.wikipedia.org/wiki/Approximate_Bayesian_computation). Now I hope that would be clear.


title: How to choose number of lags for Zivot unit root test?

text: The fUnitRoots package considers `lag=2` as the default value in the `urzaTest()` function (Zivot & Andrews unit root test), while for the `ur.za()` function in the urca package the default is `lag=NULL`. What is the "criteria" to choose a good number of lags?


title: Reporting results of a logistic regression

text: I have the following logistic regression output:               Coefficients:                 Estimate Std. Error z value Pr(>|z|)         (Intercept)   0.5716     0.1734   3.297 0.000978 ***     R1           -0.4662     0.2183  -2.136 0.032697 *       R2           -0.5270     0.2590  -2.035 0.041898 *        Is it appropriate to report this in the following way: Beta coefficient, Odds ratio, Zvalue, P value. If yes, how can I obtain the Odds ratio?


title: What does a non positive definite covariance matrix tell me about my data?

text: I have a number of multivariate observations and would like to evaluate the probability density across all variables. It is assumed that the data is normally distributed. At low numbers of variables everything works as I would expect, but moving to greater numbers results in the covariance matrix becoming non positive definite. I have reduced the problem in Matlab to:               load raw_data.mat; % matrix number-of-values x number of variables     Sigma = cov(data);     [R,err] = cholcov(Sigma, 0); % Test for pos-def done in mvnpdf.      If err>0 then Sigma is not positive definite. Is there anything that I can do in order to evaluate my experimental data at higher dimensions? Does it tell me anything useful about my data? I'm somewhat of a beginner in this area so apologies if I've missed out something obvious.


title: Interpreting the change in random effect variance over time in a GLMM

text: I am measuring the change in random effect (random intercept) over time by running models such as `glmer( bull~p1+p2+p3+p4+p5+(1|school),subset=(yr=2011),data=bull, family=binomial(link=logit))` using different `subset=(yr=20??)` for each year. In the comments to my other question, it was stated that another way to do this was to use dummy variables for each year, and specify them in the model as random slopes, thus: `glmer(bull~p1+p2+p3+p4+p5+(i1+i2+i3+i4+i5+i6-1|school),data=bull, family=binomial(link=logit))` where `i1...i6` are the indicator variables of the different years. My question is: how is the output of this model to be interpreted in terms of the longitudinal change in random effects variance ?


title: How to calculate the correlation between smoking patients and heart attacks for each country?

text: I have a set of data on smoking patients and heart attack on different time points for 25 countries. How can I calculate the correlation between the smoking patients and heart attacks for each country? My data look like this table1:                        Smoking             Heart Attack     Country  1988  1984  2010    1988  1984  2010     -------  ----  ----  ----    ----  ----  ----     Congo    1200  1146   675     900   400   550     Nigeria  1100   786   765     950   568   590     Zimbabwe 1098   897   900    1098   769   865     ...      I am looking for a solution in R.


title: Frequentist /Bayesian hybrid for model selection: make sense?

text: I was told that this hybrid is wrong and I should never use it for any serious writing:               We estimated all of the ARMA models considered in this paper using the standard      procedure of maximum likelihood. We based model selection on the Schwarz Bayesian      Criterion, with the requirement that all parameters were significant (at the 5% level).      True or False? Why?


title: glmnet converges extremely slowly for multinomial with mass function concentrated at one point

text: If I have a response vector that is mostly one class label, convergence of Lasso in `glmnet` will be extremely slow, for example:                require(glmnet)      s <- matrix(rnorm(2000*1000),ncol=1000)      f <- rep(c(-10,0,10,20),2000)[sample(2000)]      f[200:2000] <- 0      fit <- glmnet(s,f,family='multinomial',nlambda = 2)      What settings should I tweak to improve this? I have hacked together a subsampling solution but would prefer a more robust method. So basically time complexity is infeasible at this number of features (or it's getting stuck on something).


title: Error in normal approximation to a uniform sum distribution

text: One naive method for approximating a normal distribution is to add together perhaps $100$ IID random variables uniformly distributed on $[0,1]$, then recenter and rescale, relying on the Central Limit Theorem. ( _Side note_ : There are more accurate methods such as the Box-Muller transform.) The sum of IID $U(0,1)$ random variables is known as the uniform sum distribution or Irwin-Hall distribution. **How large is the error in approximating a uniform sum distribution by a normal distribution?** Whenever this type of question comes up for approximating the sum of IID random variables, people (including me) bring up the Berry-Esseen Theorem, which is an effective version of the Central Limit Theorem given that the third moment exists: $$|F_n(x) - \Phi(x)| \le \frac{C \rho}{\sigma^3 \sqrt n} $$ where $F_n$ is the cumulative distribution function for the rescaled sum of $n$ IID random variables, $\rho$ is the absolute third central moment $E|(X-EX)^3|$, $\sigma$ is the standard deviation, and $C$ is an absolute constant which can be taken to be $1$ or even $1/2$. This is unsatisfactory. It seems to me that the Berry-Esseen estimate is closest to sharp on binomial distributions which are discrete, with the largest error at $0$ for a symmetric binomial distribution. The largest error comes at the largest jump. However, the uniform sum distribution has no jumps. **Numerical tests suggest that the error shrinks more rapidly than $c/\sqrt n$.** Using $C=1/2$, the Berry-Esseen estimate is $$|F_n(x) - \Phi(x)| \le \frac{\frac12 \frac{1}{32}}{\frac{1}{\sqrt{12}^3} \sqrt n} \approx \frac{0.650}{\sqrt n}$$ which for $n=10,20,40$ is about $0.205$, $0.145$, and $0.103$, respectively. The actual maximum differences for $n=10, 20, 40$ appear to be about $0.00281$, $0.00139$, and $0.000692$, respectively, which are much smaller and appear to fall as $c/n$ instead of $c/\sqrt n$.


title: How to determine the best number of weak classifiers to use in adaboost without overfitting the data

text: I was thinking by using validation but not quite sure how to go with it. Please list some papers or ideas on how. This is for multiclass problem (using one vs all approach). I think each class/category has its own speed in converging to a good strong classifier hence would need different number of weak classifiers to avoid overfitting...


title: Using BIC to estimate the number of k in KMEANS

text: I am currently trying to compute the BIC for my toy data set (ofc iris (: ). I want to reproduce the results as shown here (Fig. 5). That paper is also my source for the BIC formulas. I have 2 problems with this:   * **Notation:**     * $n_i$ = number of elements in cluster $i$     * $C_i$ = center coordinates of cluster $i$     * $x_j$ = data points assigned to cluster $i$     * $m$ = number of clusters 1) The variance as defined in Eq. (2): $$ \sum_i = \frac{1}{n_i-m}\sum_{j=1}^{n_i}\Vert x_j - C_i \Vert^2 $$ As far as I can see it is problematic and not covered that the variance can be negative when there are more clusters $m$ than elements in the cluster. Is this correct? 2) I just cannot make my code work to compute the correct BIC. Hopefully there is no error, but it would be highly appreciated if someone could check. The whole equation can be found in Eq. (5) in the paper. I am using scikit learn for everything right now (to justify the keyword :P ).               from sklearn import cluster     from scipy.spatial import distance     import sklearn.datasets     from sklearn.preprocessing import StandardScaler     import matplotlib.pyplot as plt     import numpy as np          def compute_bic(kmeans,X):         """         Computes the BIC metric for a given clusters              Parameters:         -----------------------------------------         kmeans:  List of clustering object from scikit learn              X     :  multidimension np array of data points              Returns:         -----------------------------------------         BIC value         """         # assign centers and labels         centers = [kmeans.cluster_centers_]         labels  = kmeans.labels_         #number of clusters         m = kmeans.n_clusters         # size of the clusters         n = np.bincount(labels)         #size of data set         N, d = X.shape              #compute variance for all clusters beforehand         cl_var = [(1.0 / (n[i] - m)) * sum(distance.cdist(X[np.where(labels == i)], [centers[0][i]], 'euclidean')**2)  for i in xrange(m)]              const_term = 0.5 * m * np.log10(N)              BIC = np.sum([n[i] * np.log10(n[i]) -                n[i] * np.log10(N) -              ((n[i] * d) / 2) * np.log10(2*np.pi) -               (n[i] / 2) * np.log10(cl_var[i]) -              ((n[i] - m) / 2) for i in xrange(m)]) - const_term              return(BIC)                    # IRIS DATA     iris = sklearn.datasets.load_iris()     X = iris.data[:, :4]  # extract only the features     #Xs = StandardScaler().fit_transform(X)     Y = iris.target          ks = range(1,10)          # run 9 times kmeans and save each result in the KMeans object     KMeans = [cluster.KMeans(n_clusters = i, init="k-means++").fit(X) for i in ks]          # now run for each cluster the BIC computation     BIC = [compute_bic(kmeansi,X) for kmeansi in KMeans]          plt.plot(ks,BIC,'r-o')     plt.title("iris data  (cluster vs BIC)")     plt.xlabel("# clusters")     plt.ylabel("# BIC")      My results for the BIC look like this: ![](http://s24.postimg.org/vwu97fr6d/figure_1.png) Which is not even close to what I have expected and also makes no sense... I looked at the equations now for some time and am not getting any further locating my mistake ):


title: Good references for time series?

text: I am wondering if anyone has book references for time series. I would like something comparable (in popularity) to the 'ESL' or to 'Machine learning' from Murphy in the machine learning field. Does anyone knows what are the most complete (in term of methods scope) books containing all about exponential smoothing (all of them), arima, sarima , arch, garch, neural network for time series, kalman filters, ... etc ?


title: How to interpret the p-value in this toy example (two coin tosses)

text: The p-value is said to be the probability of the respective event _assuming that $H_0$ is true_. The simplest possible toy example are two coin tosses. The 2-tailed $H_0$ would be that you consider the coin fair, $H_1$ in this case is that you consider it biased to one side or the other. For one 1-tailed test $H_0$ would be that it is either fair or biased towards e.g. heads. $H_1$ would then be that you consider it biased towards tails (the other side) in this case. You can experiment with this toy example in R, first the two-tailed case:               > binom.test(2,2,alternative="two.sided")              Exact binomial test          data:  2 and 2     number of successes = 2, number of trials = 2, p-value = 0.5     alternative hypothesis: true probability of success is not equal to 0.5     95 percent confidence interval:      0.1581139 1.0000000     sample estimates:     probability of success                           1      And now for the one-tailed case:               > binom.test(2,2,alternative="greater")              Exact binomial test          data:  2 and 2     number of successes = 2, number of trials = 2, p-value = 0.25     alternative hypothesis: true probability of success is greater than 0.5     95 percent confidence interval:      0.2236068 1.0000000     sample estimates:     probability of success                           1       _My question is: What is the correct interpretation of the p-value in the latter case?_ Here you assume under $H_0$ that the coin is either fair or biased towards one side but get a result which seems to be biased towards the other side. But the p-value says $0.25$ which is not the probability under the condition $H_0$ but again assuming a fair coin (one coin combination out of four possible results assuming a fair coin) - as in the two-tailed case! Is this in some way connected to the fact that you effectively halved your p-value by going from two-tailed to one-tailed _although_ you threw the exact same (!) combination in both cases?


title: Kalman filter conceptual question

text: I'm using the function `dlmGibbsDIG` (Gibbs sampler) in the `dlm`package from R to estimate the unknown variances. The output are the unknown variances together with the (saved) sampled states. Then the `dlmFilter` function performing the kalman filter only need as arguments the observations and the specified model (with the averaged variances gived by the output of `dlmGibbsDIG`) but I'm confused about the role of the thetas (states) simulated by the Gibbs sampler above.


title: Detecting structures in matrices of variable size (link spam)

text: I'm trying to create an algorithm that would detect link spam structures in directed graphs of nodes (converted to matrices). Each structure is represented by a sparse matrix of links between the central page and 2 levels of links. I was going to use neural networks, however the matrices that I have in my training set differ in size. In image analysis problem this is easy - you just resize the images to identical size. However in "link spam" problem I can't (can I?) resize the matrix. Any suggestions how to approach this?


title: Trying to use Cholesky decomposition of covariance matrix to sample error ellipsoid

text: I'm trying to construct an error ellipsoid from a covariance matrix (which exists for a 3D point) and then sample consistent xyz points in this region. In a previous question when I asked about this (on math.stackexchange, here) a user suggested that I compute the Cholesky decomposition of the matrix and then: > sample Gaussian with identity matrix as the covariance (easy to do by > sampling each coordinate as a Gaussian) and then apply the linear > transformation from the Cholesky decomposition matrix to get your sampled > point from the ellipsoid covariance Gaussian I've had a go at following this but I assume I've misunderstood to at least some extent as the result isn't consistent with what I'd expect. The steps I took were: 1) Calculate the Cholesky decomposition of the covariance matrix.   2) Sample each initial vertex point as a Gaussian with width 1 to generate (x', y', z')   3) Multiply (x',y',z') by the Cholesky decomposition matrix for the newly generated point. I'm certain this isn't correct, but don't have the experience to know exactly what is wrong. Thanks in advance! **Edit to add some specifics:** To add some more details, the covariance matrix of the initial point is: $ \begin{pmatrix} 10.0115 & -10.6835 & 5.18024 \\ -10.6835 & 11.4009 & -5.52798 \\ 5.18024 & -5.52798 & 2.77646 \\ \end{pmatrix} $ I then calculate the Cholesky decomposition, which is: $ \begin{pmatrix} 3.164095& 0& 0\\ -3.376478& 0.017131& 0\\ 1.637195& -0.001619& 0.309921\\ \end{pmatrix} $ And the initial xyz point is at (35.5361, -37.2661, 22.521). As a check before sampling the Gaussian that this multiplication would yield something sensible I multiply: $ \begin{pmatrix} 35.5361& 0& 0\\ 0& -37.2661& 0\\ 0& 0& 22.521\\ \end{pmatrix} $ by the Cholesky decomposition matrix (assuming that the Gaussian sampling would look something like this), and then I was expecting this to yield some new xyz point consistent with the original. The actual result is: $ \begin{pmatrix} 112.44& 0& 0\\ 0& -0.6384065& 0\\ 0& 0& 6.97973\\ \end{pmatrix} $


title: What is a good similarity measure to use when missing data is a significant issue?

text: I have a list of cities that I want to compare in terms of their similarity. Each city can described by a large but finite number of characteristics but most of them will have missing data for some random number of characteristics. If I consider each city to be a n-dimensional vector with with missing data in some dimensions, I'm thinking of using the cosine similarity to compare the cities so that only dimensions with data in both cities will be considered in the similarity. Is cosine similarity the most appropriate similarity measure to use when there's missing data in random dimensions? Since some cities may have a very high similarity result just by matching perhaps two dimensions while others aren't considered due to missing data - how can I show the relative "quality" of the results?


title: Different HC3 standard error estimates when normalising weights for weighted least squares fit using Python statsmodels

text: When I normalize the weights I use for fitting a line with weighted least squares, the parameters of the fitted line and the 'normal' standard errors stay exactly the same, as I would expect. The HC3 standard error estimates, however, change completely. I get the feeling that I am missing something quite important here, but... well, I am missing it... Here is some test code in Python:               import statsmodels.api as sm     import numpy as np          myData = np.array([1, 1, 2, 2, 3, 3, 4, 4], dtype=float)     myIndex = sm.add_constant(range(len(myData)))          myWeights = np.array([100,10,100,10,100,10,100,10], dtype=float)          fit1 = sm.WLS(myData, myIndex, weights=myWeights).fit()          print "Parameters: %s" % fit1.params     print "Normal standard errors: %s" % fit1.bse     print "HC3 estimates: %s" % fit1.HC3_se          # Normalise the weights     myWeights /= myWeights.sum()          fit2 = sm.WLS(myData, myIndex, weights=myWeights).fit()          print "Parameters: %s" % fit2.params     print "Normal standard errors: %s" % fit2.bse     print "HC3 estimates: %s" % fit2.HC3_se      Which produces:               Parameters: [ 0.9796748   0.49186992]     Normal standard errors: [ 0.09876738  0.02581648]     HC3 estimates: [ 0.00976334  0.00314918]     Parameters: [ 0.9796748   0.49186992]     Normal standard errors: [ 0.09876738  0.02581648]     HC3 estimates: [ 44.54587593   0.75750668]


title: When forcing intercept of 0 in linear regression is acceptable/advisable

text: I have a regression model to estimate the completion time of a process, based on various factors. I have 200 trials of these processes, where the 9 factors being measured vary widely. When I perform a linear regression of the 9 factors (and all 2 and 3 factor interactions), with no explicit intercept, I get an adjusted R${^2}$ of 0.915, if I force the intercept to 0 I get an adjusted R${^2}$ of 0.953. My intention of forcing the intercept to 0 was to ensure that trials that completed in very short amounts of time (< 1 second) did not result in predictions of < 0. Setting the intercept to 0 did not help with this. So my question is threefold. 1) When is it acceptable/advisable to force an intercept? 2) Does the improved R${^2}$ actually mean the model is a better fit (the plot of fitted vs measured did look better)? 3) Is there a way of ensuring the fitted values are all > 0?


title: Interpreting probit results in Stata

text: I have read the Stata manual but I am still having trouble understanding how to interpret these results. I am looking to determine how each parameter affects the outcome given a 1-unit change. The code I am using in Stata is as follows:               probit hhaffbyinvest ///     femintmig femdommig maleintmig maledommig ///     femeduc maleeduc lifecycl ///     fmeduc members ratio_workers ///     hmale ///     commun1 commun2 ///     if (fem1stmig >= 1991 | male1stmig >= 1991) & (femhomenow == 1 | malehomenow == 1) ///     & country == 8          est store guat      _hhaffbyinvest_ is whether a household has reported investing remittances. The line below that is the migration experiences of male and female heads of household, then education, lifecycle of household, size of household, etc... I am running this in three countries. The results I get are as follows: ![probit analysis](https://lh5.googleusercontent.com/-8OEt8tlwlQY/Trf1jfLsoLI/AAAAAAAAAAk/k2iNKrio- RE/s512/probit.jpg) So, for example, how do I interpret 1.4275 for _femintmig_ (a female head or spouse migrated internationally) with respect to the binary dependent variable (household invested remittances)? This is the first column in Guatemala.


title: R oneway.test with equal and non equal variance

text: can anybody tell me why the `oneway.test()` results in a p-value of NA when var.equal is set to FALSE but in a p-value of 0.7173 when it is set to TRUE Here is my example:               x <- c(3.921973, 5.703782, 3.921973, 3.921973, 3.921973, 6.075346, 3.921973, 3.921973, 3.921973, 3.921973)     y <- c(4.424847, 4.424847, 4.424847, 4.424847, 4.424847, 4.424847, 4.424847, 4.424847)          d <- c(x,y)     f <- c(rep("a", 10), rep("b", 8))          oneway.test(d ~ f, var.equal = FALSE)          oneway.test(d ~ f, var.equal = TRUE)      Thanks for you help!