title: Plot precision recall using a vector of ranks

text: I have developed a predictor which scores a datapoint, $x$ on how likely it is to belong to a particular class, $C$. The higher the score the more likely it is that the datapoint belongs to the class. The predictor takes as input a number of other datapoints (as training data). Typically I might have say 10000 non class members, and 50 class members. I am evaluating performance by taking my 50 class members, and hiding the class membership of each member in turn. I then obtain a score for my hidden class member and compare to the scores of all the other 99950 non-class members to produce a rank. So by repeating this I get 50 different ranks, one for each of the 50 class members respectively from when they were hidden from the predictor. I will then repeat this for a different class type, which may have e.g. 30 members, and get a list of ranks for these 30 genes. I repeat this for a number of different classes, of different sizes. So I end up with a vector of ranks for each class. Is there an easy way to plot this information as a precision-recall curve? So far I've done it by flattening out the rank vector to get a list of all ranks across all classes, and then for i in 1:10000 calculating precision as $\frac{(\text{number of ranks} < i)}{(i*10000)}$ and recall as $\frac{(\text{number of ranks} < i)}{ \text{total number of ranks}}$. However this makes for a slightly strange looking precision recall curver, which looks almost like two assymetric bell curves, which isn't how prec-recall curves should look! You may be able to see, unless the link gets blocked: http://picturepush.com/public/7618805 I hope that made some sense. I would appreciate any pointers anyone can give me on this. Let me know if I need to clarify anything, I fear I have waffled on a bit.


title: Conditional regression for 2 SD below / above the mean

text: I want to interpret an interaction. Usually a conditional regression analysis is done with 1 SD below/above the mean (see Aiken & West, 1993, p. 18). My result is a contrast with $p=.068$ if I take 1.5 SDs $p=.044$. with 2 SDs $p=.033$. It seems that contrasts get stronger by increasing values the number of SDs. Thus, is it appropriate to report contrasts for 2 SDs in a paper? Why does it have to be 1 SDs?


title: Correlation in space (influence of building function on number of people passing through some point)

text: I am in a need of solving such problem: I have two data sets: number of people passing nearby some point in space (two points were analysed) and the function of buildings around those points. How do I find the correlation (or test if influence exists) of building functions to number of people passing by. For example: The drawing below shows idealized situation. Color of the rectangle denotes building function (Bank, bakery etc.). In the middle of those two building groups (crossroads) number of passers by was counted. How to check if building functions influence people's choice of road crossing point? ![enter image description here](http://i.stack.imgur.com/sfbns.png)


title: Find and plot cumulative distribution function

text: We toss a fair coin. If a head appears then $X\sim N(0,1)$. If a tail appears then $X=1$ with probability $1$. Find the unconditional CDF of $X$ and plot it.


title: Cross Vaidation in Factor Analysis

text: I am trying to conduct an EFA with a sample size of 150 respondents. I would also like to use cross-validation but my professor says that the sample is not big enough for that. Is that true?


title: Measures of entropy/information: distinguish clustered configurations that would have the same information entropy

text: Let us consider the configuration of a 2D system and the standard definition of entropy $H=-\sum_{i=1}^{m}p_{i}\cdot \log(p_{i})$. Let us suppose that I can describe the state of my system by a 2D distribution over a square grid and suppose to consider two configurations (i.e. two distributions over these square grid ) which have the same entropy. I would like to know if there exist modified measures of entropy/information which take into account also clustering so that a clustered configuration of my system is no more degenerate with respect to a more sparse one.


title: Mean survival time for a Gamma survival function

text: I was not able to find the formula to allow me to calculate the mean and the median survival time for a Gamma survival function. Does anyone knows how to calculate these?


title: Energy function of RBM

text: The Hammersley-Clifford Theorem tells us that the distribution of a RBM must be Gibbs since it is Markov Random Field, but how to prove that its energy function must be of the form: $$E = -\sum_{i,j} w_{ij} \, v_i \, h_j -\sum_i \alpha_i \, v_i - \sum_i \beta_i \, h_i$$


title: Are latent variable models modelling causality?

text: Is the purpose of latent variable models to model causality, where the causes are not observable i.e. latent? Are latent variables modelling causes of the observable variables? Thanks and regards!


title: How to present the results of numerical comparisons

text: I ran a comparison where I estimated the parameters of a model using $S$ different methods several times. More specifically I simulated $M$ datasets $\boldsymbol{y}_1, \dots, \boldsymbol{y}_m$ from a model parametrized by $\boldsymbol{\theta}$, where $\boldsymbol{\theta}$ is $p$-dimensional. Then I have estimated the model parameters by using each dataset and each method. As a results I have $M \times p$ point estimates for each statistical method. So far I'm presenting the results by using a table where each element is an average squared error (over the $M$ runs):                               Model 1 Average Squared Errors               theta_1    theta_2    theta_3    theta_4           Method1     0.1        0.01      0.091      1.11      Method2     0.3       0.003      0.047      0.79      Method3     0.24       0.03      0.042      0.08     Method4     0.65       0.01      0.007       0.4      Best       Method1    Method2   Method4   Method3       In statistical methodology literature this is seems to be the standard way of presenting results, but I find it quite hard to interpret. In particular, it is not clear which method is the winner. Drawing conclusions is even more difficult when you try the methods on several models, i.e. you have a second table:                               Model 2 Average Squared Errors               theta_1    theta_2    theta_3    theta_4           Method1     0.5        0.14      0.01       0.3      Method2     0.1        0.03      0.007      0.7      Method3     0.14       0.023     0.1        0.8     Method4     0.05       0.1       0.02       0.4      Best       Method4    Method2   Method1   Method1       In my case I am using 5 methods and 6 model, so I can't ask the reader to go through page after page of tables. So if you really wanted to come up with a ranking of the methods, how would you pool the results together to come up with it? Finally this kind of tables are really boring to read, and there must be more entertaining ways of presenting such results!


title: Standarizing testing data to mean zero and unit variance on its own or from training data

text: In standarizing the data for supervised classification, one can do the following: 1) Extract mean and variance for each feature from training, use it to standarize both training and testing data. 2) Compute the mean and variance for entire data, training and testing, and standarizing both with the same statistics. 3) Separately standarize the training and testing data. These produce different classification results on a 10 fold cross validation. 1) and 2) are about the same results, with 2) being a little higher classification performance. 3) Produces the best results, by 5% improvement. Is it OK to do 3? Which one is better to pursue and why? Thank you


title: Gambler's ruin, fun problem

text: _This is mostly for fun and because I haven't found it on here. Hopefully nobody looks up the answer online for it ruins the challenge. Best answer will be chosen on who can solve this in their own way. If you look up the answer at least give it your own flavor._ ## Gambler's Ruin A gambler has a certain amount of money $M$ and is playing a game of chance with some win probability less than 1. Every time he wins, he raises his stake to a certain fraction, $1/N$, of his bankroll, where $N$ is a positive number. The gambler doesn't reduce his stake when he loses Every time he wins, he'll raise his stake to $M/N$, or his bankroll divided by $N$. When $M= 1000$ and $N=4$, for example, he'll gamble $250 M$ each time going forward. Should he win, he'll raise it again. Should he lose, he'll keep his stake at $250 M$. If he keeps at it, what are his expected winnings?


title: Can anyone help me with understanding ||w||

text: I am studying Machine Learning (SVMs), and I found this $||w||$ in front of me. I have not seen it again, and Google does not allow me to search for it. Does anyone knows what is $w$, what is $||w||$, how it is called, and why we need that? Many thanks


title: Sampling from a dynamic population

text: I need to create a sample of a given size from a population. However, the population is dynamic, that is, comes as a stream of items, and every item has a "time stamp" based on its location in this stream. I don't know beforehand how long is this stream, but I need that in the chosen sample, the time stamps will look as a random sequence. Practically, I can't "save" the whole population, and then choose a sample. I must somehow retain a sample throughout the process, and end up with the required sample of the gien size. Here is an idea I had: Suppose you need to sample m items, and the unknown population size total is n. \- choose a number k which is greater than m and smaller than n, and for which you can handle populations of size k \- for the first k elements in the stream, create a random sample of m elements from it and retain it as the "working sample" \- for the next k (that is the second chunk of k items), sample about m/2 elements from it, and replace randomly chosen m/2 elements from your working sample, wit the new sampled items. now you have an m sample from 2k items \- so on and so forth... for the i'th chunk of k items from the stream, choose about m/i elements randomly and substitute them randomly into your working sample... \- do it until the stream is over. your working sample is the result does this algorithm create a good sample? are there better ways to do it?


title: Computing Lorenz asymmetry cofficient (S)

text: Damgaard and Weiner in Describing Inequality in Plant Size or Fecundity discuss computing the Lorenz asymmetry coefficient, S as well as confidence intervals for the estimates of S for a sample. D A Seekell, P D'Odorico and M L Pace in Virtual water transfers unlikely to redress inequality in global water use also discuss the Lorenz asymmetry coefficient and reference the above article. I can not find a discussion anywhere on how to specifically compute the confidence intervals for the Lorenz asymmetry coefficient. Spicifically I am looking for: 1) The formulas to compute same. 2) An implementation in open source software that performs the calculation.


title: Beta Weights Equivalent?

text: I know that the beta weight tells you the relative importance of a predictor (IV) in predicting the criterion (DV). However, is it true that beta weights can not be compared across different samples? (i.e. if I run two different regressions, both with the same independent variables, but different dependent variables, I still cannot compare the beta weights for the same IV across the two samples)? If one cannot compare beta weights across samples, is there a method to do so? I'm really wanting to know if the relative importance of a predictor (IV) stays constant, grows, or diminishes over time. Thanks!


title: Most used name for Dirichlet Compound Multinomial dist (aka Dirichlet-Multinomial, Multivariate Polya)?

text: I'm doing a lot of work on the Wikipedia page for this distribution: http://en.wikipedia.org/wiki/Dirichlet_compound_multinomial_distribution The actual name of the page is "Multivariate Pólya distribution". However, in the text I've been writing, I actually use "Dirichlet compound multinomial distribution", abbreviated to DCM. I've also seen just "Dirichlet-multinomial distribution". Which one is the most used? The Google search test reveals:   * 20,000+ for "multivariate polya"   * 46,000+ for "dirichlet compound multinomial"   * 24,000+ for "dirichlet multinomial"   * 30,000+ for "dirichlet multinomial distribution" But this doesn't really tell us anything, esp. given the non-intuitive values of the last two searches. My sense is that "Multivariate Polya" (with or without the accent) is the least used, but this is just an intuition. Also, if we call it "Dirichlet-multinomial distribution", how do we abbreviate it, e.g. in a mathematical description? DM(), DirMult()? Or is it just left as-is? Note also that the article "Estimating a Dirichlet distribution" by Thomas P. Minka calls it simply the "Polya distribution". (Or rather, he introduces it as the "Dirichlet-multinomial distribution", with "Polya distribution" as a synonym, but then uses the latter term in the rest of the discussion.) However, the term "Polya distribution" seems to have lots of other possible definitions referring to various quite distinct distributions, so this name is not good as the canonical name.


title: Errors when fitting a SARIMAX model in R

text: Currently working on big data sets, I would like to fit a SARIMAX model to it in order to forecast future values. Here are my questions : **Q1.** As a matter of information, does anybody know how does R compute the fitting of an ARIMA/SARIMA model ?               themodel = arima(auch[0:700,1], order = c(0,1,1), seasonal = list(order = c(0,1,1), period = 144), xreg=tmp[0:700])      For instance, this formula works well with an average calculation time - around 30 sec/1 minute - while : **Q2.** If I try to fit another set of parameters to the same data :               themodel = arima(auch[0:700,1], order = c(2,1,5), seasonal = list(order = c(1,2,8), period = 144), xreg=tmp[0:700])      I get the following error which I don't really understand :               Error in makeARIMA(trarma[[1L]], trarma[[2L]], Delta, kappa, SSinit) : maximum supported lag is 350      Lag is supposed to be `144`. Unless `D = 2` makes it so the second differentiation exceeds the `maximum supported lag` ? If so, is there anyway to work around this limitation ? **Edit :** From the SARIMA general equation which is $\phi(B)\Phi(B^s)W_t = \theta(B)\Theta(B^s)Z_t$ we can guess that the maximum lag for a given model `SARIMA(p,d,q)(P,D,Q)144` is `max((p+Px144), (q+Q*144))` (correct me if I'm wrong because I am not a specialist). My problem is I need to try orders of `(p,d,q,P,D,Q)` greater than 1 to find a better model to fit to my data set. Any idea how can I overcome this limit of 350 ? **Q3.** While doing the same command line as above, i.e. :               themodel = arima(auch[0:700,1], order = c(2,1,5), seasonal = list(order = c(1,2,2), period = 144), xreg=tmp[0:700])      I also sometimes get the following error :               Error in makeARIMA(trarma[[1L]], trarma[[2L]], Delta, kappa, SSinit) : cannot allocate memory block of size 7.0 Gb      But I tried `memory.size()` and `memory.limit()` which gave me :               memory.size()     [1] 109.27     memory.limit()     [1] 1535      I think I understood that the memory space may not be contiguous, which would cause the memory allocation problem, but it seems strange to me since I worked with `arima` and bigger data sets without problems before.


title: Independent t-tests. One- or two-tailed?

text: For an independent t-test, if my hypothesis states there will be no significant difference between the two groups, is that a one- or two-tailed test?


title: Compare fits of model to transformed and untransformed response

text: I want to compare data that proportions among three different groups e.g.:                ID Group Prop.Nitrogen      1    A     0.89      2    A     0.85      3    B     0.92      4    B     0.97      Following Wharton and Hui (doi:10.1890/10-0340.11) I though I'd see if these data would be better dealt with using a logit transformed. When I look at diagnostic plots for linear models on the transformed and un- transformed data they look very similar with no obvious problems, and there are only small differences in estimated parameters. However, I'd still like to be able to say something about how well the model fits the transformed and untransformed versions of the data - I know I can't compare AIC values directly. Is there a correction and I can make to examine this? Or should I be taking a different approach?


title: Principal Components and Noise

text: When doing principal components it is intuitively clear that noise accumulates more towards the last components. What would be the formal explanation of how and why this happens?


title: Simulating Convergence in Probability to a constant

text: Asymptotic results cannot be _proven_ by computer simulation, because they are statements involving the concept of infinity. But we should be able to obtain a sense that things do indeed march the way theory tells us. Consider the theoretical result $$\lim_{n\rightarrow\infty}P(|X_n|>\epsilon) = 0, \qquad \epsilon >0$$ where $X_n$ is a function of $n$ random variables, say identically and independently distributed. This says that $X_n$ converges in probability to zero. The archetypal example here I guess is the case where $X_n$ is the sample mean minus the common expected value of the i.i.d. r.v.'s of the sample, $$X_n = \frac 1n\sum_{i=1}^nY_i - E[Y_1]$$ **QUESTION:** **How could we convincingly show to somebody that the above relation "materializes in the real world", by using computer simulation results from necessarily finite samples?** Please do note that I specifically chose convergence to a _constant_. I provide below my approach as an answer, and I hope for better ones. **UPDATE:** Something in the back of my head bothered me -and I found out what. I dug up an older question where a most interesting discussion went on in the comments to one of the answers. In there, @Cardinal provided an example of an estimator that it is consistent but its variance remains non-zero and finite asymptotically. So **a tougher variant of my question becomes: how do we show by simulation that a statistic converges in probability to a constant, when this statistic maintains non-zero and finite variance asymptotically?**


title: Practical problems with difficult posteriors

text: I'm looking for difficult Bayesian inference problems to test out different Monte Carlo sampling methods. I've mostly been looking at Hamiltonian Monte Carlo based algorithms and in particular, I've been using the trick of using the hessian at the Maximum A Posteriori point for the mass matrix. What are some practical inference problems with difficult posteriors to sample from?


title: How to test differences in budget allocations

text: Suppose a grant program allows city administrators to allocate the grant across three mutually exclusive categories `1, 2, 3`. City administrators are surveyed to elicit their preferred allocations. A response by administrator `i` is a triplet `R_i = (p1, p2, p3)`, where by definition: `0 <= pi <= 100` for all `i` in `{1,2,3}` and `p3 = 100 - p1 - p2`. In words, `p1, p2, p3` represent the proportion of grant money he or she would allocate to each of three categories `1, 2, 3`. Suppose 80 administrators are surveyed. Of these 34 selected at random had previously received a prompt meant to shift the distribution of allocations in favor of `p1`, at the expense of the other two items. The motivation is that the grantor of money has preference ordering `1 > 2 > 3`. I want to test the null that elicited budget allocations across treatment and control are the same, against the one-sided alternative of a location shift towards `p1`. Due to small sample size, I am mostly interested in non- parametric approaches. An obvious possibility is to test `p1_T = p1_C`, for treatment and control, against `p1_T > p1_C`. But I am not sure this makes the most effective use of the information, as any such difference must also imply a decline in `p2_T` or `p3_T` or both relative to control. In addition, I am not sure marginal test of homogeneity (like Chi square test) offer a directed alternative. For these reasons I am considering a weighted response like `Y = p1*1 + p2*2 + p3*3` for each administrator, as the categories are ordered. The alternative would then be `Y_T < Y_C`. Is this correct? Any suggestions and rationale?


title: Quartiles in Excel

text: I am interested in the definition of quartile that is usually used when you're in basic statistics. I have a Stat 101 type book and it just gives an intuitive definition. "About one quarter of the data falls on or below the first quartile..." But, it gives an example where it calculates Q1, Q2, and Q3 for the set of data               5, 7, 9, 10, 11, 13, 14, 15, 16, 17, 18, 18, 20, 21, 37      Since there are 15 pieces of data, it chooses 15 as the median, Q2. It then splits up the remaining data into two halves, 5 through 14, and 16 through 37. These each contain 7 pieces of data and they find the median of each of these sets, 10 and 18, as Q1 and Q3, respectively. This is how I would calculate it myself. I looked at Wikipedia's article and it gives 2 methods. One agrees with the above, and one says you could also include the median 15 in both sets (but you wouldn't include the median if it was the average of the two middle numbers in the case of an even number of data points). This all makes sense to me. But, then I checked Excel to see how Excel calculates it. I am using Excel 2010, which has 3 different functions. Quartile was available in 2007 and previous versions. It seems they want you to stop using this in 2010 but it's still available. Quartile.Inc is new but agrees exactly with Quartile as far as I can tell. And, there is Quartile.Exc as well. Both of the last 2 are new in 2010 I believe. This time, I just tried using the integers 1, 2, 3, ..., 10. I'm expecting Excel to give median of 5.5, Q1 of 3, and Q3 of 8. The method from the statistics book, as well as both methods on Wikipedia would give these answers, since the median is the average of the middle two numbers. Excel gives               quartile number, Quartile.Inc, Quartile.Exc     1,               3.25,         2.75      2,               5.5,          5.5     3,               7.75,         8.25      Neither of these agree with what I have previously talked about. The descriptions in the help file for Excel are: Quartile.Inc - Returns the quartile of a data set, based on percentile values from 0..1, inclusive. Quartile.Exc - Returns the quartile of the data set, based on percentile values from 0..1, exclusive. Can any one help me understand this definition Excel is using? Thanks


title: If I have good reason to expect a certain value as the sample mean and I obtain something very different, should I obtain another sample?

text: In my specific case, I am measuring the mean of a sample for different values of a factor (let's say the possibilities are integers from 1 to 100). I have very good reason to expect the mean of each sample to be roughly proportional to the value of the factor; in any case, related values for the factor should give related values for the outcome. However, a few sample means deviate quite a lot from the expected value (and the neighbouring means), and I suspect that it is simply due to chance. How I do I deal with the deviating samples? I cannot just obtain new samples until I get a value which fits better; that would interfere with the underlying distribution (besides simply being dishonest). But maybe I could obtain just one other sample and replace the old one regardless of the outcome? Or is it possible to obtain a new and _bigger_ sample? If so, how much bigger must it be? Can I pool the old and the new sample? Or maybe I should just present the data that I got the first time around and note that based on predictions, that particular sample is likely a fluke?


title: When to use proportional odds model?

text: Suppose we fit a proportional odds model to some data with $Y$ the response variable being ordinal. If we get 0 estimated for all the coefficients (except the intercepts), what does that mean? How would we interpret the model?


title: Unfold degrees of freedom of interaction in factorial design analysis in R

text: I've been trying to analyse data from a factorial design with missing cells without much success until now. I'm using Anova from car package, that, according with I've read, can handle unbalanced data with type III sum of squares. But when I detect an interaction effect I don't know how to unfold the degree of freedom of the interaction and discover what is the best combination of factors. I think I have to do something like in this article > da Rocha, F et al. An approach to the decomposition of interaction in a > factorial experiment with five factors, Acta Sci., Agron. 34(1), 2012. Some one know's how to unfold the degrees of freedom? Here follows a reproducible example               fdata <- read.table("http://labpib.ffclrp.usp.br/~rsilva/fdata.txt")     library(car)     options(contrasts=c("contr.sum","contr.poly"))     fdata.hi <- aov(log(resp+1) ~ PIE*NS*La*Li*Lf, data=fdata)     # trying to analyse as a full factorial design     Anova(fdata.hi, type="III", singular.ok=TRUE)     # trying to analyse as an one-way anova     library(agricolae)     model1 <- aov(log(resp+1) ~ ALL, fdata)     HSD.test(model1, "ALL")          x11(); par(mfrow=c(3,2))      apply(combn(2:6,2)[,1:6], 2,            function(x) interaction.plot(fdata[,x[1]], fdata[,x[2]], fdata$resp))         x11(); par(mfrow=c(2,2))         apply(combn(2:6,2)[,7:10], 2,                function(x) interaction.plot(fdata[,x[1]], fdata[,x[2]], fdata$resp))      par(mfrow=c(1,1))       That is what I'm doing until now, try to verify graphically what is the best factor combination, knowing that I have a significant interaction among four factors.


title: What is the minimum sample size for panel data analysis?

text: I am working with a panel data and want to perform a regression analysis (OLS, fixed effects model or random effects model). The panel data is unbalanced with $T$ varying from 7 to 13, and $N=11$. I am confused about having a small number of observations (around 100). What is the minimum time observations and cross-sectional observations to perform a statistically significant panel analysis. And what are the alternatives if my dataset is too small?


title: EM maximum likelihood estimation for Weibull distribution

text: **Note:** _I am posting a question from a former student of mine unable to post on his own for technical reasons._ Given an iid sample $x_1,\ldots,x_n$ from a Weibull distribution with pdf $$ f(x) = k x^{k-1} e^{-x^k} \quad x>0 $$ is there a missing variable representation and hence an EM (expectation-maximisation) algorithm to find the MLE of $k$?


title: Finding uniform minimum variance unbiased estimators

text: $X_1, X_2, X_3, \ldots{}, X_n$, i.i.d., follow $\mathcal{N}(0, \theta^2)$, $\theta > 0$. What are the UMVUEs (Uniform Minimum Variance Unbiased Estimators) of $\theta$ as well as $\Phi(\tfrac{1}{\theta})$, where $\Phi$ is the cumulative distribution function of $\mathcal{N}(0,1)$?


title: Computing overhead of statistical models for training?

text: Could someone provide overhead of the following model for training (With respect to input size or if there are any relevant parameters). Overhead I mean somewhat like asymptotic time complexity form.   1. Multiple linear regression models (Least Sqaures)   2. Decision Tree (C 4.5 OR ID3)   3. K nearest neighbor  Thank you.


title: Presenting the error term in a quantile regression specification

text: Let $Y_i$ be the response and $X_i$ be the independent variables. Whenever I've seen a quantile regression specification they'll go: $Q_{\tau}(Y_i | X_i) = a(\tau) + b(\tau) X_i$ Or, alternatively: $Q_{Y_i}(\tau|X_i) = a(\tau) + b(\tau) X_i$ where $X_i$ is a vector of independent variable observations for cross-section $i$ and $b(\tau)$ is a vector of coefficient estimates for some $\tau \in [0,1]$. These are the two alternatives I see throughout the different examples in "Economic Applications of Quantile Regressions". However, I want to specify an error term in my quantile regression equation. How would I go about doing this? Further, I'd like to ask why those worked examples in that book (all from different papers) don't specify an error term?


title: Distribution hypothesis testing - what is the point of doing it if you can't "accept" your null hypothesis?

text: Various hypothesis tests, such as the $\chi^{2}$ GOF test, Kolmogorov-Smirnov, Anderson-Darling, etc., follow this basic format: $H_0$: The data follow the given distribution. $H_1$: The data do not follow the given distribution. Typically, one assesses the claim that some given data follows some given distribution, and if one rejects $H_0$, the data is not a good fit for the given distribution at some $\alpha$ level. But what if we don't reject $H_0$? I've always been taught that one cannot "accept" $H_0$, so basically, we do not evidence to reject $H_0$. That is, there is no evidence that we reject that the data follow the given distribution. Thus, my question is, what is the point of performing such testing if we can't conclude whether or not the data follow a given distribution?


title: Information gain and information gain ratio: Do I have to pick just one?

text: When splitting attributes while constructing a decision tree, i can use information gain or information gain ratio to try and determine the best value to split the tree on.   * I'd use information gain when I have smaller numbers of distinct values (mostly nominal and ordinal values that are not continuous)   * I'd use information gain ratio for attributes which have large numbers of distinct values (mainly continuous values) That said, do I necessarily have to choose _one_ method for use with all attributes across the entire tree? My gut feeling is no, because information gain or information gain ratio for comparison to other attributes is not a relevant comparison; it's only relevant for comparison to a different value of the same attribute for the purposes of producing a split. I'm thinking that one might be able to go so far as to swap them for different splits in the tree, depending on the training data. If the subset of values for that attribute is small compared to previous splits where the distinct values were large, one could use the value that's more biased to producing a more useful result.


title: How to find probability with multiple pieces of evidence?

text: I'm trying to find the probability of P5 having C1, C3, or C4. Here are the constraints of the problem. C1 - C5: Are cards that can be held by 1 person. P1 - P5: Are players that each have 1 card. X means the person does not have the card. O means the person has the card ? means it's unknown who has the card.                P1 P2 P3 P4 P5          C1 ?   X   X   ?   ?          C2 X  O   X   X   X          C3 ?   X   X   ?   ?          C4 ?.  X.  X   X.  ?          C5 X.  X. O.  X.  X      I've calculated by looking at each row: C1 has a 1/3 chance of being in P5 C3 has a 1/3 chance of being in P5 C4 has a 1/2 chance of being in P5 I've calculated by looking at P5 column: There is a 1/3 chance that they have C1, C3, or C4 I'm thinking there is a way to combine these probabilities to give me an answer as to the chance of P5 having C1, C3, or C4. Will you please show me the math and let me know what this is called in statistics?


title: Bayes' classification in R

text: For a machine learning class I am taking, on our first homework assignment we are given the following problem that has me stuck:               Consider the following simulated data set:          set.seed(123)     n<-100     X<-runif(n)     Y<-rbinom(n,1,exp(0.5+X)/(1+exp(0.5+X)))          a) Find the Bayes' classifier     b) Construct an empirical version of the Bayes' classifier using MLE (you can use      the glm function)      I don't understand how to find the Bayes' classifier using R. I can find it algebraically, but how do you implement Bayes' classifiers in R? When I search around, the only sources I can find are on "Naive Bayes' classifiers", which don't appear to be the same thing. This: http://en.wikipedia.org/wiki/Bayes_classifier Is the Bayes' classifier I want to find, but I can't find any sources on it for R. Further, even if I did know how to find the Bayes' classifier, I don't understand what the difference would be between finding it and constructing an empirical version using MLE. The question doesn't even make sense to me. How do I use the glm function to use MLE to construct a classifier? I imagine it has something to do with fitting a logistic model, but I don't understand how to use the glm function in the way I am being asked to? I suspect I might just be getting caught up in the terminology/notation and confusing myself unnecessarily. Anyone have any pointers for how to get started on this? I'm not asking anyone to code it for me, but it would be nice if someone could point me in the right direction.


title: When has a Bayesian approach been critical to addressing a theory, hypothesis or problem?

text: A question was recently posted to an listserve that I subscribe to asking (perhaps cynically?) when a Bayesian approach has been crucial to "getting the job done" in addressing questions in the field of ecology. I'm wondering, in general, about when a Bayesian approach has been essential to progress in a particular field. In ecology, Bayesian methods seem to be most frequently used in applied situations with big, complicated data sets, so I'd be especially interested in circumstances that relate to important or classic theories or hypotheses in a field. For example, in ecology, Bayesian methods appear to be the only way to fit complex hierarchical models and obtain accurate estimates of things like the size of a population of animals or the survival rate of an individual in a population of critters. I'm not familiar with instances when progress was made on burning theoretical questions because a Bayesian approach was used, though this could be because ecological theory is often addressed with reductionistic experiments in an ANOVA-like framework where p-values are the historically valued currency.


title: What is the effect of scaling down covariance in the formula for rho?

text: The rho formula divides the covariance by SD of X and SD of Y. Does it result in a valid estimate of the true correlation (i.e., in the population)?


title: What is the best way to detect trends in item popularity?

text: Assume I have kept track of the number of views per page and I want to show trending pages. I am thinking about applying some statistical tests like Log-likelihood/chi- squared to detect breakthroughs from the popularity expectation. Can somebody point me to some resources related to this matter?


title: Use matrix feature for machine learning or cluster analysis

text: I have a bunch of features that I would like to use for classification/machine learning and cluster analysis. Normally I use single point values or transformations of values for features and everything is fine Now I would like to use a matrix as a feature. The matrix is probably going to be a fairly big (say 50x50) but will only be filled with 1's and 0's. It is pretty much an 'image' matrix. It is the shape/pattern of the matrix entries which is important. Is there anyway I can **easily** use the matrix as a feature for machine learning? I know I could use each matrix entry, say Row1Column1 as a feature and then give it a value, but then I would have 2500 features from my 50x50 matrix, which is what I am trying to get away from. Any ideas would be greatly appreciated.


title: Evaluation of dynamic model

text: This is kind of a beginner level question. I want to know how to assess the performance of a model I have created. Suppose that we have an identification system aboard a spaceship. The job of the system is to identify the objects detected by the sensors of the ship. For example, the sensors may detect asteroids, UFOs, satellites, and so on. The identification system must be able to provide a response at any time, and so it is refreshed with new sensor inputs every couple seconds. More information about a detected object may become available as time passes. The 'lifespan' of a detected object may be for a couple of hours, or only for a few minutes. It is important to correctly identify the objects, because how the spaceship interacts with them may vary from one classification to another. The goal is to assess how well the identification system performs. What techniques exist to do this? Normally, for a classification problem, we can look at the number of correct and incorrect classifications and just look at the accuracy or confusion matrix. But the temporal aspect in this case is important - Maybe the identification system can produce a correct result after an hours' worth of data inputs. But it would be much better if it found the correct answer after a much shorter time span. How can I do this? What metrics are established for this kind of problem?


title: Combining instance-based learning with regression analysis to improve predictions

text: I have a table that contains items. Each day its possible that an incident can occur to some items but equally no incident may happen. I wish to see if its possible to create a prediction for the next day based on the information below               **Date      Item 1    Incident       Location**     20130701    40          12           ES     20130702    50           6           ES     20120701    60          10           ES     20120702    20           8           ES       This could be a regression problem (at the moment I'm only interested in seeing the next days prediction) but I have an idea that taking the same time period for the previous year and combining it with this years data based on week number (e.g. week 37 this year with week 37 last year) would enrich the prediction because the weather could be similar or the conditions for the incident could be similar. This might make it possible to get a better prediction. I'm using rapidminer for this but would know if anyone would know how to implement this correctly or even point me to some research papers that may have covered this. I can get up to five years of data so in essence week 37 five times but for the different years.


title: Comparision of slopes

text: How do i compare (p-value) regression slopes for 2 groups fitted with the same model but with different sample size? Model: y=cov_1*X1+...+cov_5*X5+slope_group_1*X6 and similarly for group 2. I followed What is the correct way to test for significant differences between coefficients? example and computed the p-value as: deltab = stats_group_1.beta(7)-stats_group_2.beta(7); t=deltab/se; p = 2*tcdf(-abs(t), se); Can anyone please help me verify this test?


title: When is brownian covariance less appropriate than linear covariance?

text: I've just been introduced (vaguely) to brownian/distance covariance. It seems particularly useful in many non-linear situations, when testing for dependence. But it doesn't seem to be used very often, even though cov/corr is often used for non-linear/chaotic data. That has me thinking that there must be some draw-backs, so what are they, and why doesn't everyone just always use distance covariance?


title: S shaped curve for logistic regression

text: I was trying to fit logistic regression with a binary response variable and a continuous predictor variable. But when I plotted the fitted probabilities vs the predictor, I got an almost straight line instead of the expected sigmoid curve. What can be the reason for this?


title: How to calculate the posterior distribution given Inverse Gamma conjugate prior?

text: I have a state-space model (actually belonging to a Kalman Filter) as in the given graphical model: ![Kalman](http://i.stack.imgur.com/bJgXb.png) This is a typical 1 dimensional Linear Dynamic System model. The variances $Q$ and $R$ have Inverse Gamma priors and the coefficient $A$ has a Gaussian prior. We know the exact values of $x_1$ and $B$ a priori. Now, what I want is to find a closed-form formula for the posterior probability of $Q$ which is $$P(Q|x_1,x_2,x_3,y_1,y_2,y_3,A,R,B)$$ From the full joint distribution of the model, I have managed to reduce this posterior distribution to the formula below: $$P(Q|x_1,x_2,x_3,y_1,y_2,y_3,A,R,B) = \frac{IG(Q|\alpha,\beta)\prod_{i=1}^3 N(y_i|Bx_i,Q)}{\int_{Q}IG(Q|\alpha,\beta)\prod_{i=1}^3 N(y_i|Bx_i,Q)dQ}$$ $IG$ is the inverse gamma and $N$ is the Gaussian distribution. Since we know the means of the Gaussians exactly, Inverse Gamma should be a conjugate prior for the variance $Q$ and the posterior should be distributed as Inverse Gamma as well, but I have failed to show that. How can I manage to construct the posterior for $Q$ in this case and is my conjugate prior choice is true? Thanks in advance.


title: Given 3 discriminant functions I was able to classify 2D patterns, but how do I plot decision boundaries via matplotlib?

text: I have implemented the discriminant function and was able to classify the 2D patterns (via Python), but I have troubles thinking about an approach to plot the decision boundaries. Hope anyone has an idea how I can achieve this. More specifically, I am wondering if I need to derive the separate decision boundaries via g1(x) = g2(x), g1(x) = g3(x), g2(x) = g3(x) from an analytical solution or if there is a way around that in matplotlib that would allow me to use the equation of the discriminant function (which I use to classify the patterns below) directly to plot the decision boundaries and or decision regions. The discriminant function would be defined as: $ g_i(\pmb{x}) = \pmb{x}^{\,t} \bigg( - \frac{1}{2} \Sigma_i^{-1} \bigg) \pmb{x} + \bigg( \Sigma_i^{-1} \pmb{\mu}_{i}\bigg)^t \pmb x + \bigg( -\frac{1}{2} \pmb{\mu}_{i}^{\,t} \Sigma_{i}^{-1} \pmb{\mu}_{i} -\frac{1}{2} ln(|\Sigma_i|)\bigg)$ So that I ended up with the 3 discriminant functions $\quad g_1(\pmb{x}) = \pmb{x}^{\,t} \bigg( - \frac{1}{2} \Sigma_1^{-1} \bigg) \pmb{x} + \bigg( \Sigma_1^{-1} \pmb{\mu}_{\,1}\bigg)^t \pmb x + \bigg( -\frac{1}{2} \pmb{\mu}_{\,1}^{\,t} \Sigma_{1}^{-1} \pmb{\mu}_{\,1} -\frac{1}{2} ln(|\Sigma_1|)\bigg) \\\ \quad g_2(\pmb{x}) = \pmb{x}^{\,t} \bigg( - \frac{1}{2} \Sigma_2^{-1} \bigg) \pmb{x} + \bigg( \Sigma_2^{-1} \pmb{\mu}_{\,2}\bigg)^t \pmb x + \bigg( -\frac{1}{2} \pmb{\mu}_{\,2}^{\,t} \Sigma_{2}^{-1} \pmb{\mu}_{\,2} -\frac{1}{2} ln(|\Sigma_2|)\bigg) \\\ \quad g_3(\pmb{x}) = \pmb{x}^{\,t} \bigg( - \frac{1}{2} \Sigma_3^{-1} \bigg) \pmb{x} + \bigg( \Sigma_3^{-1} \pmb{\mu}_{\,3}\bigg)^t \pmb x + \bigg( -\frac{1}{2} \pmb{\mu}_{\,3}^{\,t} \Sigma_{3}^{-1} \pmb{\mu}_{\,3} -\frac{1}{2} ln(|\Sigma_3|)\bigg)$ And the following parameters are known: $p([x_1, x_2]^t |\omega_1) ∼ N([0,0]^t,4I), \\\ p([x_1, x_2]^t |\omega_2) ∼ N([10,0]^t,4I), \\\ p([x_1, x_2]^t |\omega_3) ∼ N([5,5]^t,5I),$ I tried to solve the discriminant functions analytically, but I may have made some mistakes, because g1(x) - g2(x) = 0 would lead to weird solutions. So here is what I got so far: * * * **Let:** $\pmb{W}_{i} = - \frac{1}{2} \Sigma_i^{-1}\\\ \pmb{w}_i = \Sigma_i^{-1} \pmb{\mu}_{\,i}\\\ \omega_{i0} = \bigg( -\frac{1}{2} \pmb{\mu}_{\,i}^{\,t}\; \Sigma_{i}^{-1} \pmb{\mu}_{\,i} -\frac{1}{2} ln(|\Sigma_i|)\bigg)$ * * * $ \pmb{W}_{1} = \bigg[ \begin{array}{cc} -(1/8) & 0\\\ 0 & -(1/8) \\\ \end{array} \bigg] $ $ \pmb{w}_{1} = \bigg[ \begin{array}{cc} (1/4) & 0\\\ 0 & (1/4) \\\ \end{array} \bigg] \cdot \bigg[ \begin{array}{c} 0 \\\ 0 \\\ \end{array} \bigg] = \bigg[ \begin{array}{c} 0 \\\ 0 \\\ \end{array} \bigg]$ $ \omega_{10} = -\frac{1}{2} [0 \quad 0 ] \bigg[ \begin{array}{cc} (1/4) & 0\\\ 0 & (1/4) \\\ \end{array} \bigg] \cdot \bigg[ \begin{array}{c} 0 \\\ 0 \\\ \end{array} \bigg] \- ln(4) = -ln(4)$ with $ \quad g_1(\pmb{x}) = \pmb{x}^{\,t} \bigg( - \frac{1}{2} \Sigma_1^{-1} \bigg) \pmb{x} + \bigg( \Sigma_1^{-1} \pmb{\mu}_{\,1}\bigg)^t \pmb x + \bigg( -\frac{1}{2} \pmb{\mu}_{\,1}^{\,t} \Sigma_{1}^{-1} \pmb{\mu}_{\,1} -\frac{1}{2} ln(|\Sigma_1|)\bigg) $ $ \Rightarrow g_1(\pmb{x}) = \pmb{x}^{t} \bigg[ \begin{array}{cc} -(1/8) & 0\\\ 0 & -(1/8) \\\ \end{array} \bigg] \pmb{x} + [0 \quad 0 ] \; \pmb x - ln(4) = \pmb{x}^{t} (-(1/8) \; \pmb{x}) - ln(4) $ $ \pmb{W}_{2} = \bigg[ \begin{array}{cc} -(1/8) & 0\\\ 0 & -(1/8) \\\ \end{array} \bigg] $ $ \pmb{w}_{2} = \bigg[ \begin{array}{cc} (1/4) & 0\\\ 0 & (1/4) \\\ \end{array} \bigg] \cdot \bigg[ \begin{array}{c} 10 \\\ 0 \\\ \end{array} \bigg] = \bigg[ \begin{array}{c} 2.5 \\\ 0 \\\ \end{array} \bigg]$ $ \omega_{20} = -\frac{1}{2} [10 \quad 0 ] \bigg[ \begin{array}{cc} (1/4) & 0\\\ 0 & (1/4) \\\ \end{array} \bigg] \cdot \bigg[ \begin{array}{c} 10 \\\ 0 \\\ \end{array} \bigg] \- ln(4) = -12.5-ln(4)$ $ \Rightarrow g_2(\pmb{x}) = \pmb{x}^{t} \bigg[ \begin{array}{cc} -(1/8) & 0\\\ 0 & -(1/8) \\\ \end{array} \bigg] \pmb{x} + [2.5 \quad 0 ]\;\pmb x - ln(4) = \pmb{x}^{t} \cdot( -(1/8) \; \pmb{x}) + 2.5 x_1 - ln(4) $ $ \pmb{W}_{3} = \bigg[ \begin{array}{cc} -(1/10) & 0\\\ 0 & -(1/10) \\\ \end{array} \bigg] $ $ \pmb{w}_{3} = \bigg[ \begin{array}{cc} (1/5) & 0\\\ 0 & (1/5) \\\ \end{array} \bigg] \cdot \bigg[ \begin{array}{c} 5 \\\ 5 \\\ \end{array} \bigg] = \bigg[ \begin{array}{c} 1 \\\ 1 \\\ \end{array} \bigg]$ $ \omega_{30} = -\frac{1}{2} [5 \quad 5 ] \bigg[ \begin{array}{cc} (1/5) & 0\\\ 0 & (1/5) \\\ \end{array} \bigg] \cdot \bigg[ \begin{array}{c} 5 \\\ 5 \\\ \end{array} \bigg] \- ln(5) = 5 -ln(5)$ $ \Rightarrow g_3(\pmb{x}) = \pmb{x}^{t} \bigg[ \begin{array}{cc} -(1/10) & 0\\\ 0 & -(1/10) \\\ \end{array} \bigg] \pmb{x} + [1 \quad 1 ]\; \pmb x - ln(4) = \pmb{x}^{t} \cdot ( -(1/10)\; \pmb x )+ \; {x_1} + {x_2} + 5 - ln(5) $ As I mentioned above, the results may not be correct :(   I would really appreciate it if someone has so time to check whether this is correct! * * * So what I have done so far: I was given a dataset of 2D patterns that can stem from 3 classes, 1, 2, or 3. The model is Gaussian and the parameters are known: When I plot the training data set (with class labels known) it would look like this:               # Training Dataset     f, ax = plt.subplots(figsize=(7, 7))     ax.scatter(train_set[train_set[:,2] == 1][:,0], train_set[train_set[:,2] == 1][:,1], \                marker='o', color='green', s=40, alpha=0.5, label='$\omega_1$')     ax.scatter(train_set[train_set[:,2] == 2][:,0], train_set[train_set[:,2] == 2][:,1], \                marker='^', color='red', s=40, alpha=0.5, label='$\omega_2$')     ax.scatter(train_set[train_set[:,2] == 3][:,0], train_set[train_set[:,2] == 3][:,1], \                marker='s', color='blue', s=40, alpha=0.5, label='$\omega_3$')     plt.legend(loc='upper right')      plt.title('Training Dataset', size=20)     plt.ylabel('$x_2$', size=20)     plt.xlabel('$x_1$', size=20)     plt.show()      ![enter image description here](http://i.stack.imgur.com/MVtK7.png) So I came up with the discriminant functions (based on Bayes' rule for statistical pattern classification): then I implemented the discriminant function to classify the patterns (pretending I don't know the labels). The equation is $ g_i(\pmb{x}) = \pmb{x}^{\,t} \bigg( - \frac{1}{2} \Sigma_i^{-1} \bigg) \pmb{x} + \bigg( \Sigma_i^{-1} \pmb{\mu}_{i}\bigg)^t \pmb x + \bigg( -\frac{1}{2} \pmb{\mu}_{i}^{\,t} \Sigma_{i}^{-1} \pmb{\mu}_{i} -\frac{1}{2} ln(|\Sigma_i|)\bigg)$               def discriminant_function(x_vec, cov_mat, mu_vec):         """         Calculates the value of the discriminant function for a dx1 dimensional         sample given the covariance matrix and mean vector.              Keyword arguments:             x_vec: A dx1 dimensional numpy array representing the sample.             cov_mat: numpy array of the covariance matrix.             mu_vec: dx1 dimensional numpy array of the sample mean.              Returns a float value as result of the discriminant function.              """         W_i = (-1/2) * np.linalg.inv(cov_mat)         assert(W_i.shape[0] > 1 and W_i.shape[1] > 1), 'W_i must be a matrix'              w_i = np.linalg.inv(cov_mat).dot(mu_vec)         assert(w_i.shape[0] > 1 and w_i.shape[1] == 1), 'w_i must be a column vector'              omega_i_p1 = (((-1/2) * (mu_vec).T).dot(np.linalg.inv(cov_mat))).dot(mu_vec)         omega_i_p2 = (-1/2) * np.log(np.linalg.det(cov_mat))         omega_i = omega_i_p1 - omega_i_p2         assert(omega_i.shape == (1, 1)), 'omega_i must be a scalar'              g = ((x_vec.T).dot(W_i)).dot(x_vec) + (w_i.T).dot(x_vec) + omega_i         return float(g)          import operator          def classify_data(x_vec, g, mu_vecs, cov_mats):         """         Classifies an input sample into 1 out of 3 classes determined by         maximizing the discriminant function g_i().              Keyword arguments:             x_vec: A dx1 dimensional numpy array representing the sample.             g: The discriminant function.             mu_vecs: A list of mean vectors as input for g.             cov_mats: A list of covariance matrices as input for g.              Returns a tuple (g_i()_value, class label).              """         assert(len(mu_vecs) == len(cov_mats)), 'Number of mu_vecs and cov_mats must be equal.'              g_vals = []         for m,c in zip(mu_vecs, cov_mats):              g_vals.append(g(x_vec, mu_vec=m, cov_mat=c))              max_index, max_value = max(enumerate(g_vals), key=operator.itemgetter(1))         return (max_value, max_index + 1)      And then I classified the data to calculate the error:               # Empirical Error of the training set          import prettytable          class1_as_1 = 0     class1_as_2 = 0     class1_as_3 = 0     for row in train_set[train_set[:,2] == 1]:         g = classify_data(             row[0:2],              discriminant_function,             [mu_vec_1, mu_vec_2, mu_vec_3],             [cov_mat_1, cov_mat_2, cov_mat_3]         )         if g == 2:             class1_as_2 += 1         elif g == 3:             class1_as_3 += 1         else:             class1_as_1 += 1          class2_as_1 = 0     class2_as_2 = 0     class2_as_3 = 0     for row in train_set[train_set[:,2] == 2]:         g = classify_data(             row[0:2],              discriminant_function,             [mu_vec_1, mu_vec_2, mu_vec_3],             [cov_mat_1, cov_mat_2, cov_mat_3]         )         if g == 2:             class2_as_2 += 1         elif g == 3:             class2_as_3 += 1         else:             class2_as_1 += 1          class3_as_1 = 0     class3_as_2 = 0     class3_as_3 = 0     for row in train_set[train_set[:,2] == 3]:         g = classify_data(             row[0:2],              discriminant_function,             [mu_vec_1, mu_vec_2, mu_vec_3],             [cov_mat_1, cov_mat_2, cov_mat_3]         )         if g == 2:             class3_as_2 += 1         elif g == 3:             class3_as_3 += 1         else:             class3_as_1 += 1      Which yielded:               +--------------+----------------+----------------+----------------+     | test dataset | w1 (predicted) | w2 (predicted) | w3 (predicted) |     +--------------+----------------+----------------+----------------+     | w1 (actual)  |      327       |       3        |       20       |     | w2 (actual)  |       0        |      323       |       27       |     | w3 (actual)  |       7        |       10       |      333       |     +--------------+----------------+----------------+----------------+     Empirical Error: 0.06 (6.38%)


title: How does the expected value relate to mean, median, etc. in a non-normal distribution?

text: How does the expected value of a continuous random variable relate to arithmetic mean, median, etc. in a non-normal distribution (eg. skew-normal)? I'm interested in any common/interesting distributions (eg. log-normal, simple bi/multimodal distributions, anything else weird and wonderful). I'm looking mostly for qualitative answers, but any but quatitative or fomulaic answers are also welcome. I'd particularly like to see any visual representations that make it clearer.


title: How to do a Hypothesis test for 10 samples?

text: I have 10 samples with sample size $N=30$ in each case. Null Hypothesis = (population mean) $\mu\le14$ Alternative Hypothesis= (population mean) $\mu>14$ confidence level = 0.95 there is an equation for hypothesis testing $Z$ $$ Z=\frac{(\bar{x}-\mu)}{\frac{SD}{\sqrt{N}}} $$ can some one tell me please how to apply this equation for the 10 samples that I have selected. can the hypothesis be tested using the 10 samples at once rather than having seperate hypothesis testing for each sample?


title: Fitting panel data with both variable and constant coefficients using R

text: I'm working with panel data and want to fit a model of the form: $y_{ij} = \alpha_{i} + \beta_{1i} x_{1ij} + \beta_{2i} x_{2ij} + \beta_{3} treatment_{ij} + \epsilon_{ij}$ The data consist of (unbalanced) time series for a number of individuals, each of which received the treatment at a different time. For each individual in the data, the coefficents on $x_1$ and $x_2$ should be different, as well as the intercepts. However, I want to estimate the coefficient on treatment across all individuals. I'm having trouble figuring out which functions to use to get this done. One possibility would be "least squares dummy variables", using lm:               lm(y ~ x1*individual + x2*individual + treatment)      where "individual" is a factor variable. My data has a large number of individuals and so this takes a lot of time and memory to run. Other than that, is there anything wrong with this approach? I've also looked into the plm package, which would have the advantage of time- demeaning the data instead of fitting for the intercepts for each individual. But I can't figure out how to make it fit a model of the form that I want. The closest I can get is:               pvcm(y ~ x1 + x2, index=c("individual","date"), model="within")      That will fit for different slopes for each individual, but I can't find a way to include the treatment variable such that the coefficient on that is the same for all individuals. I'm also reading up on lme4, but am finding the documentation difficult. Is there a way to do this using "lmer"? Any advice on which packages are applicable and how they might differ would be much appreciated!


title: Is the binomial sign test really the exact McNemar test?

text: I am calculating two different $p$ values for the McNemar test and what I believe is the corresponding _exact_ test, even for combined counts larger than 25. Take, for example the contingency table                                  Posttest     Pretest      Outcome A   Outcome B     Outcome A       xx          50     Outcome B       65          xx      The chi-square statistic is $$\chi^2=\frac{(50-65)^2}{50+65}=1.96,$$ with $$p=0.1619. \tag{R}$$ According to Sheskin (2011, Test 20, VI.3, pg 844), the exact test for these situations is essentially a binomial sign test (for a single sample) with parameter $\pi=0.5$ and the two counts equal to the two the cells of interest in the contingency table. My understanding of this is that I'm calculating the probability that, out of $50+65=115$ trials with probability of success 0.5, at least 65 of them are successes. My attempt to calculate this for the particular situation above is: $$ p = \sum_{i=65}^{115} \begin{pmatrix} 115 \\\ i \end{pmatrix} 0.5^{i}\ 0.5^{115-i} = 0.0957 \tag{Mathematica} $$ Why do these two $p$-values differ so much? Am I making a simple mistake?


title: Time and space complexity of Deep Belief Nets (DBN)

text: What is the time and space(memory) complexity of DBNs? given d:number of dimensions(attributes), n:number of records, and l:number of hidden layers.


title: fitted() function in R vs adding the residuals to the original data

text: I've found a discrepancy between the output of the fitted() function and adding the residuals to the original data set. Is the fitted() function not doing what I think it should be doing?


title: Is AR(1) a Markov process?

text: Is AR(1) process such as $y_t=\rho y_{t-1}+\varepsilon_t$ a Markov process? If it is, then VAR(1) is the vector version of Markov process?


title: What is the distribution of the product of a Bernoulli & a normal random variable?

text: I am struggling to find a solution on finding the distribution of following random variable: $$Y = Z \cdot |X|$$ here, $Z$ is a random variable takes 1 or -1 with equal probability, and $X$ is a standard normal variate, and $|\cdot|$ denotes absolute value. Can somebody help me with some pointer from where I should start? Thanks for your help.


title: Confusion related to convex optimization of a function in a svm paper

text: I was reading this paper http://www.ist.temple.edu/~vucetic/documents/wang11kdd.pdf related to adaptive multi-hyperplane machine for non linear classification In that paper, they have mentioned about multiclass SVM, with multiple weights for each class. The loss for any classification is $l(x_n,y_n) = max_{i\epsilon y\\\\\y_n}(0,1 + max g(i,x_n) - g(y_n,x_n))$ where $y_n$ is the label for the nth example and $x_n$ is the features. I have this confusion when they do the training of this algorithm. They call this SVM MM(Multiple Hyperplane). They say the convex-approximated problem is defined as $min_{W}P(W|z) = \frac{\lambda}{2}||W||^2 + \frac{1}{N}\sum_{n=1}^{N}l_{cvx}(W;(x_n,y_n);z_n)$ where they have the concave term $-g(y_n,z_n)$ replaced with the convex term $-w^T_{y_n,z_n}x_n$. I am not sure if I have described it clearly. But I am going to attach the screenshot of the paper as well. The thing is I didn't get what's the difference between $-g(y_n,z_n)$ and $-w^T_{y_n,z_n}x_n$. They seem the same term to me. I might be asking a lot. But can anyone provide some info? ![enter image description here](http://i.stack.imgur.com/33iB6.png) ![enter image description here](http://i.stack.imgur.com/TNtk0.png) I have marked by the red rectangle the part that I didn't understand. I might be asking a lot. But I didn't get that part. Why is it so?


title: QQ plot for Generalized linear model in R

text: I want to draw a QQ plot for Negative Binomial distribution to see if the Generalized linear model fit well or not. But since there is no normality assumption exist in my model, normal QQ plot will not be suitable in this case. What kind of plot should I do? And how can I do it in R? I can not simply use qqplot(negative binomial, data). I am writing a function including this plot, and the data for observations and variables are the things I will put in after I write down the function, and they are not fixed values. So I do not have a fixed form for the negative binomial part. However, what I want to know is the fit of the model, so I need to create something else to check that.


title: contrapositive of probability

text: If P(A|B) = 95%, then is P(B'|A') also 95%? The subject is hypothesis testing. If the null hypothesis is true and there is a 95% probability that the data should pass the test, then does failing the test imply the null hypothesis is wrong with 95% chance?


title: Detecting a trend in 4 samples

text: I have 4 levels for a single factor, which is called "Compliance". Compliance has an ordinal scale (0, 1, 2, 3) and it is computed as the number of rules a team complies with. My hypothesis is the following: "the more you comply, the more effective you are". Effectiveness is a real positive number. I have 4 independent samples, one for each compliance. My questions are: 1) Can I apply Jonckheere-Terpstra? 2) If I apply J-T and it succeeds, I can only conclude the existence of a weak monotony. I would like to check whether there are significant differences among pairs. In particular, I would like to check the following: C=0 vs C=1, C=1 vs C=2 and C=2 vs C=3. What test should I run? Should I apply any correction? Thanks.


title: Random Forest or other machine learning techniques with paired samples

text: I am trying to select features from paired samples and was wondering if there are methods out there that adapt random Forest or other machine learning algorithms to paired samples i.e. if randomForest was a machine learning analog to a standard FDR-corrected ANOVA approach, what would be the machine learning analog be of a paired t-test between two categories?


title: How is finding the centroid different from finding the mean?

text: When performing hierarchical clustering, one can use many metrics to measure the distance between clusters. Two such metrics imply calculation of the centroids and means of data points in the clusters. What is the difference between the mean and the centroid? Aren't these the same point in cluster?


title: Area to the right or left when computing test statistic

text: Can someone generally explain when the area to the right or area to the left is sufficient in solving a statistical problem (using z, t, or chi)? For example, in a chi-squared problem, if we're looking at a table of p-values we are generally interested in the area to the right. Why is this? What is the area to the right defined as?


title: Best model for change in scores over three time points

text: I have a response variable measured at three time points per individual (week 0, 18, and 36). I am interested in differences in the change of the response over the 36 weeks within some categorical variable X. I see two ways of modeling this.   1. One way ANOVA with response = `week_36_score - week_0_score` (this seems like the simplest option)   2. Repeated Measures ANCOVA with response = `week_18_score` and `week_36_score`. Covariate = `Week_0_score`. ### Question   * Which would you prefer (if any, maybe there is a better choice)?  I know there has been a lot of articles on this and they seem to say each has its own strengths and weaknesses. I believe the second model would have more power. I am not worried about a ceiling effect here.


title: Given a small pilot study, how to compute number of samples needed for a certain power?

text: I ran a small pilot study and computed both a p-value and Cohen's d. Now how do I compute the number of samples needed for the full study if I want a given power (e.g. 95%). I hoped I could use this table, but I can't figure out how it's calculated. **Edit:** More details on the study: We're manipulating a single variable (condition A and B). Each participant will have N readings from each condition. The question is, how do we compute how many participants we'll need given the results of a pilot study with one participant.


title: For bootstrapping, why does a higher subsample size lead to lower variance?

text: I've been working on a bootstrapping problem that's left me a little confused and wondering whether I'm doing things correctly. We have around 200 samples from a population of about 3,400, we want to bootstrap a value to estimate the total value for the 3,400 cases. Myself and a colleague both took slightly different approaches. I took the 200 samples and randomly selected 3,400 observations to create a new sample the same size as the population, then got the sum of the 3,400 values. I repeated this 10,000 and took the average and the standard deviation of all 10,000 totals. This gave me an estimate for the total value with a 95% confidence interval. My colleague did almost the exact same thing, but instead of taking a sub- sample of 3,400, his 10,000 subsamples each only had 200 in them. He got the average and standard deviation of the whole lot and multiplied it by 3,400 to get the estimate for the total. When we compared results, we found we got the exact same answer for the estimate - which is good. However, the standard deviation from his method was much bigger. From doing some research it seems like his method of resampling to the same number as the original sample is correct, but can anyone explain exactly why the difference in the standard deviation? The difference got me wondering if this is how we calculate the standard deviation at all. Should we be calculating the standard deviation of each subsample 10,000 times and estimating that the same way as the sum? Also, can anyone point to any resources/tutorials to clear things up? Thanks!


title: Kernel PCA (in R)

text: I am attempting to use the kernel PCA features in `kernlab` but am having trouble understanding the output. In particular, it's unclear what scale the results are in when I try to project the original data unto the resulting eigenvectors. Here is an example:               require(kernlab)     x = data.frame(a=1:10, b=10:1, c=sample(1:10,10), d=sample(1:10,10))     x = as.data.frame(lapply(x, scale, scale=F))     kpc = kpca(~., data=x, kernel="polydot", kpar=list(scale=1, offset=0, degree=2), features=1)      Printing `x` and the data from `x` projected unto the eigenvectors shows:               print(x)           a    b    c    d     # 1  -4.5  4.5  1.5 -2.5     # 2  -3.5  3.5 -4.5 -3.5     # 3  -2.5  2.5 -1.5 -4.5     # 4  -1.5  1.5  0.5 -1.5     # 5  -0.5  0.5 -0.5 -0.5     # 6   0.5 -0.5 -2.5  4.5     # 7   1.5 -1.5  3.5  1.5     # 8   2.5 -2.5  2.5  3.5     # 9   3.5 -3.5  4.5  2.5     #10   4.5 -4.5 -3.5  0.5          print(rotated(kpc))              [,1]      #1    68.18682      #2  -109.41619      #3   -24.34769      #4    30.26417      #5    24.81883      #6    41.32570      #7   -20.17489      #8   -37.57738      #9   -93.93368     #10   120.85431      I'm having a tough time understanding the scale of this result. Does `rotated` not scale the projection? In particular, the eigenvectors that PCA returns should be length 1. We don't know the eigenvectors in feature space, however, so well we can return is the projection of our data in that feature space using the kernel trick. The results below don't seem to be on the correct scale, however.


title: What is the test statistic in Fisher's exact test?

text: For a 2 by 2 contingency table, some said Fisher's exact test uses the count $X_{1,1}$ in the (1,1) cell in the table as the test statistic, and under null hypothesis, $X_{1,1}$ will have a hypergeometric distribution. Some said its test statistic is $$ |X_{1,1} - \mu| $$ where $\mu$ is the mean of the hypergeometric distribution (mentioned above) under null. It also said that p-values is determined based on the tabel of hypergometric distribution. I was wondering if there is some reason to subtract mean and then to take absolute value? $|X_{1,1} - \mu|$ doesn't have a hypergeometric distribution under null, does it?


title: Multi-group SEM with a categorical outcome (using lavaan)

text: I'd like to ask the SO community for some help in regard to the interpretation of a structural equation model with three groups, featuring a categorical outcome. I have found a lot of sources treating multi-group CFA with categorical factor indicators, but nothing accessible for models with continuous predictors but a categorical outcome. In a simplified form, my model has one factor, F1, with four continuous predictors, X1-X4. I am interested in the effect of F1 on a categorical variable U, which has three categories. Additionally, I am interested whether the effect of F1 differs for three different educational groups, people with high, secondary and low education levels. Variable U has two thresholds, a variance and a scale. I've established that for the three groups the condition of invariance holds for:   1. factor loadings   2. intercepts   3. means    4. thresholds   5. regression coefficients   6. residuals I am also able to impose the condition of equality on scale of y* for U for the first two groups, while keeping the scale of y* for the third group free. Variances for the latent factor F1 are non-invariant. Based on my understanding of SEM I would draw the following tentative conclusions. (A) I'd imagine since factor loadings and intercepts are equal, that F1 is measured in the same fashion across the three groups. (B) The equality of means would imply in my understanding that the average levels of F1 are equal across groups. (C) Means of F1 are equal, but the variance of F1 differs between groups. I am not sure whether I can say that the effect of F1 on U is equal between groups. The corresponding regression coefficient is equal, as are thresholds for U. Yet the scale of y* for you is only equal between the first two groups, and the variance of U differs between the three groups as well. Finally, I am unable to impose equality on the variance of U by stating the following in the lavaan model:               u ~~ c(v1,v1,v1)*u       The label `v1` appears in the model, but the coefficients still differ. My question is (1) whether I am wrong in my initial conclusions, (2) how to interpret the equality of the effect of F1 on U and (3) whether someone is familiar with lavaan to specify the syntax for fixing the variance of U. Thank you so much in advance!


title: What is the purpose of and how to use the xreg argument when fitting ARIMA model in R?

text: I am really new with R and time series. But, I have understood most concept. Part where I am (very) confused is the `xreg=` argument in `arima()` from tseries package and forecast package. As I have read most of the related threads on this site, I understand that `xreg` is used for exogenous data. See for example, How to fit an ARIMAX-model with R?. Thus, I conclude that I can fit SARIMAX or ARMAX method using `arima()`. However, I am currently reading Shumwhay's book and his web tutorial. If you check the website (close to the bottom of the page), you can found that he said > `xreg` in `arima()` does not fix ARMAX model. It is also explained in the R issues on his website. In addition, he proposed fitting ARMAX model in state space model. So, the answers found on this site and the explanation in the website are completely different. Can anyone explain me why there is such discrepancy? Which one is correct? If Shumwhay is correct, is there any function/package that can fit ARMAX or SARIMAX model in R? Sorry if it turns out I only missed some points, but I hope you can point me to the right direction.


title: What sort of problems is backpropagation best suited to solving, and what are the best alternatives to backprop for solving those problems?

text: I am developing a neuroscience-inspired training algorithm for feed-forward neural networks. The natural benchmark for comparison is backpropagation. So I need to know to know what sort of classification problems backpropagation (and thus feed-forward NNs) is most suited for. Also, since I ultimately want to move beyond backprop, what non NN algorithms are the best alternatives to backprop, when it comes to the sort of problems that backprop is good at. Or, in other words, what problems should I benchmark a feed-forward network on, and what are the algorithms to beat for those problems?


title: Negative Expectation of a constant is positive or negative

text: Say, I have a normal distribution $ X \sim N(\mu, \sigma)$, where $\mu$ unknown and $\sigma^2$ known Then, $-E[2/\sigma]=2/\sigma$ or $-2/\sigma$ This is in reference to the Fisher Information where you have $-E(\cdot)$ of something. Edit:   I was taking expectation of Fisher Information. $$-E[\frac{1}{2\sigma^4}-\frac{(X-\mu)^2}{\sigma^6}|\mu, \sigma^2]=\frac{1}{2\sigma^4}$$   where did the negative sign go?


title: Pre-processing time series data for data mining / predictive modeling input

text: What are some ways to prepare/pre-process time series data to use the series data as a predictor(s) in a predictive model (classification or regression)? Specifically, what are the methods to be considered in order to   * Get the most predicatively useful signal from the data?   * Reduce the dimensionality of the series As a concrete example, I have 90 days of ending balance data (amount on deposit in a checking account). I want to use that data to predict if the owner of the account will close it in the next 2 weeks (I have an indicator for this occurring or not 2 weeks after the end of the series). ADDITION: After reviewing the responses, I think I was looking for a list of commonly used techniques.   * Certainly there are the feature creation that Matt Krause wrote about (each customers balance series are treated separately in all these methods): Things like differences and % changes in the series values each day or week. There are sliding window aggrgations like weekly averages, min, max and standard variations. Also number of increases and decreases, indicators for changes in the balances (absolute or relative of a certain size).    * I have considered fitting a linear or polynomial regression to each serie and using the coefficients in the model as predictors.    * Other things I've wondered about is for each balance series compute the auto correlations for a maximum number of lags and use these values as predictors.   * Cluster the time series into a relatively small number of values and use the indicators for cluster membership as predictors (use dynamic time warping distance and hierarchical clustering for example). Are there others?   * For example do Fourier Transforms work here? I may post a separate question about them.


title: Confidence interval for a relative risk when one proportion is 0

text: I am analysing data from a small study. There are two treatment groups and for each, one notes success and failures. For one of the groups, the number of successes is 0. I did an exact significance test, but reviewers of the paper want a confidence interval for the relative risk - in other words, a confidence interval for $p_1/p_2$, where $p_1=0$. What should I do? The usual formula assumes that neither proportion is 0 or 1.


title: Cluster analysis of boolean vectors in R

text: I have 114 vectors with 6 boolean attributes. I saw that might be several distinct clusters in a simple visualization. K-means clustering on the transformed vectors (true = 1, false = 0) results in roughly the clusters that I had seen in the visualization. However, I am not sure what the most appropriate clustering method for this kind of data is, and how to determine the confidence in those factors (the k-means results change every time due to randomization). Should I treat the data as nominal or as numerical data? What would be the best way to do a cluster analysis on this kind of data in R?


title: Are logistic regression coefficient estimates biased when the predictor has large variance?

text: I'm simulating data from a logistic regression model:               log(p/1-p)= 0 + X      where $X \sim N(0,\sigma^2)$. After I simulate the data, I fit a logistic regression model to the data and compare the fitted regression coefficients to the actual regression coefficients. I've noticed that as I increase $\sigma$ (i.e. the variance of the original $X$ data) the fitted regression coefficient for $X$ (i.e. $\beta_1$) is consistently greater than 1 (however, the sd of the estimate also increased so 1 is still contained in the confidence interval for beta1) I was wondering why when you increase the variance, the fitted $\beta_1$'s tend to be greater than the actual $\beta_1$ (i.e. $\beta_1 = 1$), not less than? Is there a statistical explanation for this? Thanks!               beta0 = 0     beta1 = 1     sigma = 1     number_samples = 10000     genLogit = function(pos_prop,sd){         generated_data = c()          xtest = rnorm(10000,0,sd)     linpred = beta0 + (xtest * beta1)     prob = exp(linpred)/ (1+exp(linpred))          runis = runif(10000,0,1)     ytest = ifelse(runis<prob,1,0)          pos = sample(xtest[ytest ==1],floor(pos_prop*1000))     neg = sample(xtest[ytest == 0], floor((1-pos_prop)*1000))          generated_data = rbind(cbind(pos,rep(1,floor(pos_prop*1000))),cbind(neg,rep(0,floor((1-pos_prop)*1000))))     colnames(generated_data) = c('X','Y')     generated_data = data.frame(generated_data)          fit = glm(Y~X,data =generated_data, family=binomial(link="logit"))          return(fit)     }      If you run `genLogit(.5,1000)` this is generating balanced (50/50) data with X distributed normal(0,1000). Running it multiple times, I get a beta0 estimate much greater than 0.


title: Oddity in simulating Weibull survival times

text: Based on this guide to simulating survival times in SAS and R, I used the following code to generate two survival functions:               data simulation;         call streaminit(321);         beta_death = 0.9163; *Simulated ln(HR) of ~ln(2.5);         lambda_death = 0.002; *Baseline hazard of death;         beta_discharge = 0.4055; *Simulated ln(HR) of ~ln(1.5);         lambda_discharge = 0.001; *Baseline hazard of discharge;         do i = 1 to 100;             x = rand('BERN', 0.10); *10% exposure rate;             linpred_death = exp(-beta_death*x);             linpred_dis = exp(-beta_discharge*x);             t_death = rand('Weibull', 1, lambda_death* linpred_death); *Time to death;             t_dis = rand('Weibull', 1, lambda_discharge* linpred_dis); *Time to discharge;             output;         end;     run;      The weird part is when I look at the generated survival times - they're quite small. And when I _increase_ what's being called the baseline hazard, the survival times get longer. That doesn't make sense to be - decreasing the baseline hazard should, I feel, increase the survival times, should it not? What am I missing?


title: How to use SPSS output for calculating additional variable by syntax?

text: I need to use the output from certain tests like Factor Analysis to compute additional variable like socioeconomic status.


title: Inequality involving the regression function and its approximation

text: In the proof of Theorem 2.2 from the book by Devroye et al., I don't understand why $g(x) \ne g^*(x)$ implies $|\eta(x) - \tilde\eta(x)| \geq |\eta(x) - \dfrac{1}{2}|$. Note that $\eta(x) = \mathbb{E}\\{Y|X=x\\}$ is the regression function, $\tilde\eta(x)$ is an approximation of $\eta(x)$, $g^*(x)$ is the Bayes classifier $$ g^*(x) = \begin{cases} 0 & \text{if } \eta(x)\leq\dfrac{1}{2} \\\ 1 & \text{otherwise} \end{cases}$$ and finally, $g(x)$ is defined like $g^*(x)$ with $\tilde\eta(x)$ replacing $\eta(x)$.


title: Measuring correlation of point processes

text: There is a huge literature on time series analysis. My data does not seem to fit into the standard model in that it consists of event times, that is the times at which an event occurs. What is a good way to measure the correlation between two sequences of event times? Is there any literature on analyzing event time sequences? The data might be the times at which a particular neuron fires, for example. So in this case I would like to say that the firing times of two different neurons are highly correlated ( or not). I found this question looks related Analysis of cross correlation between point-processes.


title: Bootstrap vs numerical integration

text: My understanding of the bootstrap approach is based on Wasserman's framework (almost verbatim): > Let $T_n = g(X_1, ..., X_n)$ be a statistic ($X_i$ is the iid sample drawn > from distribution $F$). Suppose we want to estimate $V_F(T_n)$ - the > variance of $T_n$ given $F$. > > The bootstrap approach follows these two steps: > >   1. Estimate $V_F(T_n)$ with $V_{\hat{F}}(T_n)$, where $\hat{F}$ is the > empirical distribution function. > >   2. Approximate $V_{\hat{F}}(T_n)$ using simulation. > > Do I understand correctly that the simulation in step 2 could be replaced with a precise calculation, except that it is infeasible for practically useful values of $n$? Here's my thinking: $V_{\hat{F}}$ precisely equals an integral of $T_n(X_1, ..., X_n)d\hat{F}(X_1)d\hat{F}(X_2)...d\hat{F}(X_n)$. $\hat{F}$ is a step-function, with a finite number $n$ steps; so we can ignore all points except the $n$ points where $d\hat{F}(x)$ has non-zero mass. So the integral is precisely equal to a sum of $n^n$ terms. Once $n$ exceeds 14, a simple direct calculation is impossible. But all we're trying to do is calculate an integral. Why not replace the brute-force bootstrap simulation with any of the traditional numerical algorithms for taking integrals? Wouldn't it result in much higher precision for the same computational time? Even something as simple as splitting the sample space in sections (perhaps with smaller volumes where sample statistic varies faster), and estimating the value of the statistic in each section by using the middle point, seems to be better than blind bootstrap. What am I missing? Perhaps bootstrap works so well and so fast that there's no need to do anything more complicated? (For example, if the loss of precision in step 1 is so much greater than in step 2, then improvements to step 2 are rather useless.)


title: Estimating Markov transition probabilities from sequence data

text: I have a full set of sequences (432 observations to be precise) of 4 states $A-D$: eg $$Y=\left(\begin{array}{c c c c c c c} A& C& D&D & B & A &C\\\ B& A& A&C & A&\- &-\\\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\ B& C& A&D & A & B & A\\\ \end{array}\right)$$ **EDIT** : The observation sequences are of unequal lengths! Does this change anything? Is there a way of calculating the transition matrix $$P_{ij}(Y_{t}=j|Y_{t-1}=i)$$ in Matlab or R or similar? I think the HMM package might help. Any thoughts? eg: Estimating Markov chain probabilities


title: Comparing structural equation models with multivariate non-normality

text: I am aware that multivariate non-normality in SEM can inflate the chi-square statistic and deflate standard errors of parameter estimates. I can deal with the latter (in AMOS) using bootstrapping procedures. But here's my specific question: **while I can't correct for chi-square inflation in AMOS (as far as I am aware), does this inflation matter if I am mostly interested in comparing several models using the chi-square change, each of which will presumably be inflated in the same manner, rather than simply testing whether a single model is well-fitting?** Moreover, my best fitting model fits well by all sensible criteria, so the chi-square inflation per se is not a major concern. As this query implies, I use AMOS and haven't currently got access to M-Plus or another similar package to apply the correction to chi-square inflation they provide.


title: Which statistical method should be used in the study of stability of an enzyme?

text: I did a study monitoring the stability of an enzyme exposed to different conditions over time. Each day for 30 days, an aliquot was taken from each storage sample and analyzed in duplicate. I wanted to use statistical analysis to show mean changes during the 30 day experiment between the different conditions for analysis. Here are the statistical methods I was thinking to use:   * For condition, mean value and standard deviation (SD) values at each day will be calculated and plotted versus day    * Day 1 values of mean value will be compared across the 8 conditions using a one factor analysis of variance (ANOVA) to determine if the eight methods had comparable starting value.    * If the starting value was significantly different, then the change from Day 1 will be calculated and used for comparison of the conditions. A plot of means and SDs of change in value from Day 1 by method and day will be made    * Comparison of values from the 8 conditions across the 30 days will be made using a two factor ANOVA, with conditions and days as the two factors.    * If a significant interaction of condition by day was observed, then a one-factor ANOVA will be used for each day to compare the conditions. A p value ≤ 0.05 will be used as statistically significant    * To determine deterioration in the sample based on the condition, a one-factor ANOVA will be used to make comparisons among days for each condition. This will be followed by a one-tailed Dunnett’s test to compare Days 2-30 to Day 1. There will be up to 29 comparisons for each method (30 days-1). Some of the conditions may have less than 30 days due to the condition the sample is placed in making it unable to be tested. A reasonable experiment-wise error rate of p value < 0.01 will be used.


title: Are there popular dimensionality reduction tools for classification type supervised learning?

text: I am familiar with PCA as a linear transformation in order to align the axis of the IV space with the directions of maximal variance in order to possibly be able to reduce the dimensionality of the problem. Are there similar techniques (or can one use PCA itself) when facing a binary classification problem? Like a logistic regression? I.E. Suppose we have a dependent variable we want to predict which only takes values 0 and 1 and we have a lot of features from which to predict, but maybe there are "principal components" in the high dimensional feature space which carry the most explanatory power.


title: What do the bias units represent in a restricted boltzmann machine?

text: I'm reading up on RBMs and this is not obvious to me. I'm imagining RBMs being used for something like the Netflix Prize (since that was one of the papers I read on it). So you have a bunch of movies in the visible nodes, and the genre of those movies in the latent variables. But in the energy function, there are bias units for all of the nodes. What do they represent?


title: References about univariable vs multivariable variable selection

text: Suppose I have variables $X_j$, $j=1,\ldots,p$, some of which are correlated, and some continuous output $y$. I want to rank the variables by importance. One way is to do an association test of each variable j with the output (linear regression for each variable). This ignores the correlations between the variables, and potentially leads to selection of redundant variables in terms of prediction. Another way is to build a multivariable model of all variables (potentially with some penalisation, l1 or l2, not material here). The multivariable model implicitly accounts for the correlations between the variables, so if two variables are correlated with the output and also high correlated with each other, usually the one more correlated with the output will get a higher regression coefficient, and the weaker one will get a low coefficient (or zero in case of lasso). This is contrast with the univariable approach where both will get a high coefficient. Another way to say this is that the multivariable model includes all the potential confounders of the association of each variable j with the output but the univariable models is confounded most of the time unless it happens to look at the causal input. Therefore, I would expect the estimates from the univariable methods to be inflated relative to the true effects, with a high false positive rate of detecting the causal variables. Is anyone aware of literature discussing these intuitions more rigorously?


title: LDA, Significance of orthonormality- Trace Ratio Maximization

text: The objective of fisher linear discriminant analysis can be formulated as maximizing $\frac{Tr[X^TAX]}{Tr[X^TBX]}$ over $X$ where $A$ and $B$ are positive semi-definite with orthonormality constraints over X. What is the significance of the orthonormality constraint in here? I am also curious if reformulating this objective as $Tr[X^TAX]-\alpha Tr[X^TBX]$ is justified? $\alpha$ belongs to positive reals and the objective is a function of $X$ and $\alpha$ in this case. Was reading section II in : http://www.eecs.berkeley.edu/~jiayq/papers/tnn_traceratio.pdf and was wondering about the orthonormality constraint. I believe that orthonormality constraints in some trace minimization problems prevent a trivial solution, where the entries in X can all go to zero. The objective for Fisher's LDA seems to be having the orthonormality constraint for quite a different reason, and am not sure about it. Has it got do to something with a unique solution over rotation matrices? I don't believe that to be true either..


title: A question about a formula for lmer

text: I have been asking a number of questions about fitting mixed effects models. Thanks to a number of wonderful statisticians on here I have come far, and I thank you for you continued patience. I was going to ask this question live in chat but when I officially registered my username I lost my reputation and it has reverted from 23 to 1, so I apologise for yet another question. I now understand how to formulate a simple mixed-effects model, enter into the R package using lme4 and how to look at AICc etc. to compare models. However, one response variable may need a more complex design which I would appreciate help with if at all possible. I wish to see if habitat choice of broods at one site is different from the other. Each observation may have multiple habitats used so I have a proportion of use for each habitat type for each observation. I have six broods, three at each site. Maybe an alternate test is better, but as I only have six broods that I wish to compare (three at each beach) so using means makes the sample very small, and I cannot use chi-sq due to it being proportional data. So, I am not sure if a formula is available in a model that would then allow me to compare it to a null model to see if site has an effect on the response variable 'habitat'. However I have 7 habitats so have a variable (proportion of use) for each habitat within each observation. More importantly I would like to see how the use of each habitat affects a response variable such as mean foraging rate. Another study did this but their model contain one variable 'habitat' which made me think that they had only one habitat per observation. However they ended up calculating the Estimate for each habitat by averaging the two top models (AIC delta<2), so it seems that it may have somehow contained a formula showing each habitat use within?) The easy way would be just to add a variable for each habitat, but feel that they should be 'related' in some way. However I would still like to get an Estimate for each habitat. I hope that I am making sense. Thank you all for your ongoing patience. Rachel


title: Questions on bag of features methods

text: I've been studying the bag of features (BoF) methods. I understand the basics and how it's done but I think I still am not good at the details. So here's my understanding of such methods: Suppose you have a data base of images and you need to classify them. OK, then you must select local features in the first place. I read that there's no general rule on how they're selected, it highly depends on the application. Then, considering each of these local features is called "visual vocabulary", the same steps as bag of words is applied here, too. In bag of words, we count the important words of a document using a histogram. Here are my questions:   * When your data base of image is relatively large and you have no information of the images individually, how are those local features selected? How do you know which parts of the image contains information and are thus important to us?   * I read that these methods are simple because spatial features are ignored. Don't we lose information this way, or even, do we have enough information? If you could suggest some good reads on the topic, I'd really appreciate it but please don't publish Google results, I've read most of them (if not all of them!). In case someone is interested to know why is the word "bag" used; that's because they're unordered collections.


title: Algorithm implementation help

text: I was reading a paper by Cui et al., "Discover Breaking Events with Popular Hashtags in Twitter". I was trying to implement their hashtag instability algorithm, ![enter image description here](http://i.stack.imgur.com/2dQbs.png) I do not understand the actual algorithm, can some one please help me to convert it to a pseudo code, it will be really helpful, thanks!


title: t-test on highly skewed data

text: I have a data set with tens of thousands of observations of medical cost data. This data is highly skewed to the right and has a lot of zeros. It looks like this for two sets of people (in this case two age bands with > 3000 obs each):                Min.  1st Qu.   Median     Mean  3rd Qu.     Max.       0.0      0.0      0.0   4536.0    302.6 395300.0      Min.  1st Qu.   Median     Mean  3rd Qu.     Max.       0.0      0.0      0.0   4964.0    423.8 721700.0       If I perform Welch's t-test on this data I get a result back:               Welch Two Sample t-test          data:  x and y     t = -0.4777, df = 3366.488, p-value = 0.6329     alternative hypothesis: true difference in means is not equal to 0     95 percent confidence interval:      -2185.896  1329.358     sample estimates:     mean of x mean of y       4536.186  4964.455       I know its not correct to use a t-test on this data since its so badly non- normal. However, if I use a permutation test for the difference of the means, I get nearly the same p-value all the time (and it gets closer with more iterations). Using perm package in R and permTS with exact Monte Carlo                   Exact Permutation Test Estimated by Monte Carlo          data:  x and y     p-value = 0.6188     alternative hypothesis: true mean x - mean y is not equal to 0     sample estimates:     mean x - mean y            -428.2691           p-value estimated from 500 Monte Carlo replications     99 percent confidence interval on p-value:      0.5117552 0.7277040       Why is the permutation test statistic coming out so close to the t.test value? If I take logs of the data then I get a t.test p-value of 0.28 and the same from the permutation test. I thought the t-test values wold be more garbage than what I am getting here. This is true of many other data sets I have like this and am wondering why the t-test appears to be working when it shouldn't. My concern here is that the individual costs are not i.i.d. There are many sub-groups of people with very different cost distributions (women vs men, chronic conditions etc) that seem to voilate the iid requirement for central limit theorem, or should I not worry about that?


title: How to validate a random sample

text: I have been given two sets of data, one with 10,000 observations which is a subset of the other with 100,000 observations. It is **claimed** that the small dataset is a random sample of the large dataset. How do I verify whether this is true? I could compare the distributions, e.g. mean/median/interquartile ranges for different variables in both sets of data. However, is there a general way to test the hypothesis: "dataset x is not randomly sampled from dataset y."