title: R - possible to use bootstrapping to test restricted/unrestricted model?

text: I'm trying to use the bootstrapping method to essentially do what an $F$-test does - determining if the restricted or unrestricted model is better. My two models are below:               lm3 <- lm(gamble ~ sex + status + income + verbal, data=teengamb)     lm4 <- lm(gamble ~ income, data=teengamb)      I'm having a bit of difficulty understanding the concept behind this bootstrapping. I can sample 1000 times and get 1000 linear models for each model (2000 total). What do I do with those? Do I calculate F tests for each pair? Thanks - a bit confused, and any guidance would be great.


title: Steps to write my own general equilibrium model in MATLAB

text: I'm facing a challenge: writing my own general equilibrium model in MATLAB. I would like to ask for:   1. the basic needed knowledge in mathematics   2. the basic needed knowledge in programming with references if possible. I'm a master's degree student in economics, but I'm willing to put in the necessary efforts to achieve my aim. Many thanks in advance.


title: GAM with negative binomial and temporal autocorrelation

text: I want to fit a GAM in R. My data are negative binomial. Plus they have a temporal autocorrelation.               m1 <- gam(Y~s(X,Day), s(fDate, bs="re"), family=negbin(c(1,10), link = log), data=data)       Resid <- resid(m1)        acf(Resid)      With the `acf()` function I have detected a temporal autocorrelation. Can I apply a correlation structure in `gam()`, or do I have to use `gamm()`, since in `gamm()` I can´t use a negative binomial function, or am I doing something wrong with that? Is it basically possible to apply a negative binomial family in `gamm()` like this?               m2 <- gamm(Y~s(X,Day), random=list(Date=~1), family=negbin(c(1,10), link = log),      correlation=corAR1(form=~Day), data=data)


title: Grouped and Ungrouped Data- What's the difference?

text: Would data in two columns such as Shoe Size (variate x) and Number of Pairs (frequency f) be considered grouped or ungrouped data? Does grouped data necessarily imply the existence or classification via class intervals?


title: Nonparametric mixture model and clusters

text: I have a question about clusters that I am contemplating to treat with a nonparametric mixture approach (I think). I am working on the explanation of human comportment.   Each row of my database contains:   1. the ID of someone    2. some parameters of the environment _X_ (example: the temperature, the wind, etc.)   3. a binary variable _Y_ representing the person's reaction to the parameters (example: get sick or does not get sick because of the weather). My idea (based on intuition and not on data) is that we can gather people in a finite number of groups so that in a group, people have the same reaction to temperature (some are easily sick, others are never sick...). In a given group, more formally, the law of _Y_ conditional to the parameters _X_ is the same. I have no idea of the law of _Y_ conditional to _X_. For the parameters _X_ , I can do some hypothesis if necessary. I would like to create some cluster of people "having more or less" the same reaction to parameter. Besides, I would like to predict the reaction of a given person to a given value of the parameters (even if this event has never happened in the database). It seems to me that we can treat the problem like a nonparametric mixture model. As I don't have hypothesis on the conditional law of _Y_ , I think I will have to create it with kernels method for instance. I have found this paper. Besides, it seems to me that, in this case, each row of observation $(X_i, Y_i)$ is not a simple realization of some random variable, but $X_i$ is a realization of a random variable, and $Y_i$ is a realization of a random variable conditional to $X_i$. I don't know if it makes a difference. I have around 100000 rows. The vector $X_i$ has some discrete components, and others are continuous. I am wondering:   * Is my approach correct?   * Would you advise another point of view for this problem? I would be very interested in any references about it.   Don't hesitate to ask me to reformulate the problem statement.


title: Outlier detection using clustering and dissimilarity matrix in R

text: I have some problems in finding the outliers using clustering. The data.frame is ~20000 observations and each row has mixed types of variables(numeric, nominal and binary). What I want to do is to detect the outliers by clustering. I have calculated the dissimilarity matrix using daisy() function in R:               diss = daisy(data,metric="gower")      And I know I can use pam() and hclust() functions to do the clustering. But how do I find the outliers after that? Here is my R code to find the outliers from pam():               kmedoid = pam(diss,k=10,diss=T)     centers = kmedoid$id.med     distMat = as.matrix(diss)        distances = rep(-99,20000)        for (k in 1:20000) {          distances[k] = min(distMat[centers,k])        }        outliers = order(distances, decreasing=T)[1:5]        outliers = data[outliers]        outliers      I don't know whether it is correct, because the result seems to be pretty different each time when I tried different value of k in pam(). So the main question is: **Once I have the "kmedoid" and "hc" calculated below, how do I find the outliers?**               kmedoid = pam(data,k=10,diss=T)     hc = hclust(data)      I did search Google, but there wasn't much info about this. I am not fluent in programming, so just using existing package and function in R would help me a lot:) And **is there any better method to find the outliers?** Thanks!


title: How does regression analysis help one understand how the typical value of the dependent variable change?

text: > regression analysis helps one understand how the typical value of the > dependent variable changes... \-- > http://en.wikipedia.org/wiki/Regression_analysis What does this mean? What "typical value" is it talking about? What kind of change is it looking for? Positive or negative or what? Change in what? Percentage? I'm confused. Lessons learnt:   * greened answer is good in that it had concrete clearness   * wikipedia documents a ton of info, most of which are unreadable and unusable   * most (80%) ppl, while im sure good meaning, just dont know how to give good answers   * and only human... have so many flaws, and thereby closes good questions


title: What is $s_x$ in Pearson correlation?

text: I read from my lecture notes that Pearson correlation can be computed as $\frac{s_{xy}}{s_xs_y}$. If I have been given a finite amount of data points $(x_i,y_i),i=1,\ldots,n$, how can I compute the Pearson correlation using formula $\tfrac{s_{xy}}{s_xs_y}$? I managed to find that $s_{xy}=\dfrac{\sum_{i=1}^n (x_i-\bar{x})(y_i-\bar{y})}{n-1}$ but I didn't found the similar formula for $s_x$ and $s_y$.


title: Explain regression to 7 years old

text: Please give a clean, simple, explanatory answer that any 7 yr old can understand. You can also link to a regression guide that is very good and simple. Should be fully explanatory all the way through. That's what makes a good, not bad, writing to learn from. Some sub-questions within my main question are:   * What is the purpose of regression?    * Why do we need it?    * What does it do?    * What specific category of math does it fall under, if any?    * Can you use a calculator on the Web to automatically calculate a regression? what's a good site that has this calculator or software to download and perform this function   * What does the result of a regression tell you? or suppose to tell you?  Questions asked that are related directly to this main question:   * How does regression analysis help one understand how the typical value of the dependent variable change?   * http://stats.stackexchange.com/questions/35831/what-does-it-mean-when-it-says-with-respect-to-a-conditional-probability-distri   * What does regression analysis estimate Exactly?


title: Probability distribution for transformation of a random variable

text: Let $g(x)=1$ if $x \leq c$ and $g(x)=(1-x)/(1-c)$, where $0 \leq x \leq 1$ and $0 <c <1$. So $g$ is an non-increasing function. Define $g^{-1}(y)=\inf\\{0 \leq x \leq 1 \mid g(x) \leq y\\}$. Now $X$ is a random variable that takes on value $0, 0.1, 0.2, \ldots, 1$. Is it true that $$ \mathbb P(g(X) \leq b) = \mathbb P(g^{-1} g(X) \geq g^{-1}(b)) \leq \mathbb P(X \geq g^{-1}(b)) \,? $$


title: How is the training set constructed for multi-class SVMs?

text: Support vector machines do binary classification. If there is more than two classes, it is possible to train several classifiers instead of one. Two common approaches are training _one vs. one_ (each class against each other class) and _one vs. all_ classifiers (each class against all other classes). Say that I have relatively small sample size, e.g. N=200. Do I split the data into training and test sets separately for each classifier?


title: Cumulative distribution of Binomial Random variable

text: How to prove this: **Proposition.** If $F_n$ is the distribution function of a $\textrm{Bin}(n,p)$ random variable, then, for every real fixed $t$, the sequence $\\{F_n(t)\\}_{n=1}^\infty$ is nonincreasing. P.S. There was some heavy editing here.


title: Had statisticians predicted 2008 financial crisis?

text: Are there any statistical or econometric studies before 2008 that predicted 2008 financial crisis? Note that there are some publications that attemp to predict contagion between markets using copula concept from probability theory.


title: Learning a multivariate polynomial with dependent coefficients

text: I have a polynomial of the form $K^2((a-i)^2 + (b-j)^2 + c^2) = (ct)^2$ where $a,b,c,t$ are unknowns. I have multiple observation points for the values of $i,j,\&\ K$.   * Can I use some technique to learn the values of $a,b,c,\&\ t$ from these?   * Could some machine learning technique help me with this?   * Would non-linear regression be of some help here?  I am new to machine learning and don't have an in-depth understanding of it. Thanks.


title: Is there a nice(r) Taylor expansion of the normal quantile function?

text: Letting $\Phi$ be the CDF of the standard normal, and $\Phi^{-1}$ be the quantile function of the normal, I am looking for the Taylor series expansion of $$ \Phi^{-1}\left(\Phi\left(x\right) + \epsilon\right). $$ Using the Taylor series expansion of the inverse error function, I believe I get something like $$ \Phi^{-1}\left(\Phi\left(x\right) + \epsilon\right) = x + \sqrt{\frac{\pi}{2}}\left[2\epsilon + \frac{\pi}{12}\left(3\left(2\Phi\left(x\right)-1\right)^2 2\epsilon \+ 3\left(2\Phi\left(x\right)-1\right) 4\epsilon^2 + 8\epsilon^3\right) + \ldots \right]. $$ This is a bit ugly for my tastes, and would prefer to collect the terms in powers of $\epsilon$. Is this a known expansion? (FWIW, I am getting the $\epsilon$ from the Edgeworth expansion.)


title: Dendrograms using clustergram vs traditional ways in Matlab

text: Good Morning, By comparing the dendrograms produced by the **clustergram** object and the "manual" approach i.e. **pdist** -> **linkage** -> **dendrogram** I found they are different, but cannot find an explanation for that difference.      Clustergram documentation says that the default distance used is **'Euclidean'** and the default linkage method is **'Average'** , same parameters I used for **pdist** and **linkage** functions.      I thought it might be related with the **standardization** performed by **clustergram** , so I used zscore for standardizing first my matrix and then each column separately, but the dendrograms are still different as compared with the one produced by **clustergram**.      Can somebody explain the reason behind that difference, and how can I use the "manual" approach for getting the same dendrogram?      Many thanks in advance. Edit: I compared the results obtained by the "manual" approach with the **SPSS** results and they are the same.   Also I tested all different **'data processing'** options with **clustergram** but the dendrogram keeps being different.


title: How do you calculate standard error in a Dickey-Fuller test?

text: So in everything I've found, they tell you have to calculate $\rho$, or how to test for confidence interval for it. What I am trying to figure out is how to calculate the SE which would get us our critical value, i.e. if you get a $\rho$ of 0.8. To test it, you would do $(\rho-1)/(\text{SE}(\rho))$. To clarify, rho is our original coefficient in the AR model, it is not the transformed $(\rho-1)$ variable. Moving on, how does one get this $\text{SE}(\rho)$ manually?


title: $R^2$ correspondence for nonlinear time series

text: Is there a statistical measure for nonlinear time series data that is comparable to $R^2$ value in linear regression (giving an idea of how well the fit is)? The data is not monotonic, so I cannot use Spearman's $\rho$; and it's not linear, so Pearson's $\rho$ is not valid here either. Data looks like below (It's mostly monotonic decline. But since there is a spikes in the 11th year, it's not strictly monotonic): ![mostly monotonic decline. But since there are multiple spikes, it's not strict monotonic](http://i.stack.imgur.com/AJgys.png) Here is the fitted exponential curve (in blue) and the original data (in red) ![enter image description here](http://i.stack.imgur.com/DOToR.png) Here is the qqplot: (Doesn't seem to be a good fit although it should be a good one because of the image above?) ![enter image description here](http://i.stack.imgur.com/LsPkb.png)


title: Matrix completion approaches for healthcare big data

text: I am working on a prediction problem that leverage sparse clinical datasets. Missing data rate is in the range of 80%. 1- I am wondering if there is any example of application of matrix completion to clinical or other datasets with such a missing rates. 2- Currenlty exploring glmnet, pcaMethods and SoftImpute pacakages. I am also looking for R packages/SAS routines that can handle such sparse clinical data matrix and perform matrix completion. 3- I would like assess the reliability of my filled-in values, is there any metric or score to assess the quality of the matrix completion. Thanks in advance !


title: Are there concerns with doing a generalized additive model with a smoothing term of multiple variables?

text: In R, usually one performs a generalized additive model like so:               library(mgcv)     x <- gam(depvar ~ s(indvar1) + s(indvar2) + parametricvar)      Say I want to explore the non-parametric interaction of `indvar1` and `indvar2`. Is the following kosher? Are there issues of statistical validity with implementing such a model?               x <- gam(depvar ~ s(indvar,indvar2) + parametricvar)


title: Testing specific categories after Fisher's exact test in 2 x k tables

text: I have two variables, one with 2 categories and one with 6 categories. To test, whether there is any relation between these variables, I performed Fisher's exact test, which gave a significant result. Now I want to know for which of the categories of the second variable respective subjects have a differing distribution of the first variable (the one with 2 categories) with respect to subjects not in this specific category of the second variable. For this purpose I wanted to do the following: For each category c of the second variable (c in {1,...,6}): Group all subjects belonging to category c together and all subjects not belonging to category c, which results in a new binary variable variable2_c. Perform Fisher's exact test between variable2_c and the first variable (the one with 2 categories). Do you think this procedure is justified or do you see any problems with this? Or is there a better way to do this? Thank you in advance!


title: Is our queue processing speeding up?

text: We have been collecting statistics on a queueing and processing system that powers part of our business. We have multiple servers processing items in a queue. We recently added new servers, and we want to know whether we're processing faster. By "processing faster", I mean something like this: > Last week, all of our servers processed a combined average of X files/min > during our peak load. This week, all of our servers processed a combined > average of Y files/min during peak load. Are X and Y different enough to be > statistically significant? The data I have access to looks something like this graph, which shows the number of files processed per second for an hour-long time period. I have access to about 3 months of this type of data. ![Number of files processed per second](http://i.stack.imgur.com/Fyqxk.png) I have a minor in statistics, but I'm 10 years out of school. So please be gentle :-) Thanks so much! Alan


title: Do you need justification to transform non-normal data

text: First, a bit of background: I'm currently working on a project at work that produces variable width data. During a operation qualification I collected a bunch of data at low, nominal, and high levels. The data is all independent and I have no reason to expect anything but identical distributions part to part as the process is same with minor parameter tweaks (pressure and speed for example). The majority (15 out of 18 sample populations) passed a Anderson-Darling, but a few data sets did not. My quality engineer wants me to transform one of these data sets and use a Weibull distribution (which showed a P value of .079) just because there is a slight correlation. To clarify, the goal is to perform process capability analysis for each data set. My question(s): Is it okay to blindly transform data even though I have no rationale? I have always been taught that if you have independent and identically distributed data that the central limit theorem allows you to perform capability analysis on the data as if it were normal, is that incorrect?


title: What is lag length?

text: This might sound a bit basic and stupid, but I've been looking up all across the internet to see what IS lag length. I want to perform a EG ADF test, but for ADF, it always asks to specifiy a lag. It seems to me that whenever i look up lag length, most people just have trouble determining what it should be. I want to know what is lag length? What does it mean


title: Eigenvectors corresponding to eigenvalues

text: In R, the eigen() returns descending sorted eigenvalues. However, the eigenvectors do not correspond to these sorted eigenvalues. How do I identify the eigenvector corresponding to the ith sorted eigenvalue? One approach might be to sort by the norm of the eigenvectors (is there a procedure in R that can do this). I imagine there is a simple method I am missing.


title: How do I report Lasso results in an article?

text: I am using Lasso to reduce my number of variables and I have identified which variables had been retained after running my analysis (e.g. yrseduc, age). However, I don't know which numbers I should report in a scientific paper. I have been trying to find examples, but can't quite find the answer. I am wondering if I just need to report the parameter estimates of all the variables that were found significant? Should I just report them in a table or in the text? If in the text, what do these numbers really mean? ![Example](http://i.stack.imgur.com/8WQ3n.png) Thank you in advance, I am really at lost here.


title: Scale-invariant analysis of time series

text: When developing a general purpose time-series software, is it a good idea to make it scale invariant? How would one do that? I took a time series of around 40 points, and then multiplied by factors ranging from 10E-9 to 10E3 and then ran through the ARIMA functions of Forecast Pro and Minitab. In Forecast Pro, all resulted in the same answer (automatic modeling), whereas in Minitab, they were not. Not sure what Forecast Pro does, but they might just scale up or down all the numbers to a certain scale (let's say 100s) before running the model. Is this good idea in general?


title: How can I control the false positives rate?

text: I feel I'm pretty new to this, since some time passed since my last statistics assignment, so please bear with me. I am analyzing results of a biological experiment. Basically, I'm looking at a some graph over a genome, where each position in the genome has a value, and I'm looking for local minima (peaks). Now, I have to set some threshold since relatively high local minima also occur by chance. I can simulate the experiment computationally and get new data, but this is quite resource demanding (I can't run 1000 simulations, perhaps 100 or even only 20). What I'm currently doing is run some simulations; for each simulation: find all local minima, build a CDF for the local minima values. I then average all the simulation CDFs over all simulations so I have one 'average' CDF (CDF_simulations) which is supposed to show how would local minima distribute if everything in my genome is random. I do the same for the real data: find minima and build CDFs for their values, so I now have two CDFs - one for the real data and one for the average of the simulations. I now search for the max x such that CDF_simulations(x) / CDF_realdata(x) is < 10%. I report all minima in the real data with value < x as "true". I think this method should get me to a FP rate of 10%.   * Does this make sense?   * What is this method called and where can I find more about it?   * How should I scan the CDFs to find the right x? Sometimes for low x's CDF_simulations(x) > CDF_realdata(x).   * Where does the number of simulations come into play? Does it make sense to simply build an averaged CDF as I did? I think this is quite common, and the name FDR also pops to mind, but when reading about FDR I couldn't exactly figure how to apply it to my situation. Any comments and references (hopefully user-friendly ones) will be appreciated! Thanks Dave


title: Using the distance of p-value from alpha

text: I have two statistical tests which are inverse of each other, meaning that the null hypothesis are reversed. I want to use both the tests to take a decision. For this purpose, I am planning to do the following:   * If both the tests point to same result (by, say, rejecting the null hypothesis in test A and not rejecting the null hypothesis in test B, or vice versa), then I go ahead and take that decision   * On the other hand, if the results of the two tests are conflicting, I measure the difference of the p-values from alpha (say, 0.05) for each test, and go with the one having the largest deviation (I have some tie-breaking rules, but let's leave that here) It sounds reasonable to me, but is there some statistical ground in interpreting p-values like this? I certainly haven't seen similar application before (in my limited exposure). Edit: Let me clarify the question some more with the actual context and tests. I am testing for unit roots in time series, in order to determine if the series needs to be differenced in order to render it mean-stationary. The particular tests are KPSS test, with null being no unit roots, and ADF test with null being an unit root exists. Although the null hypotheses are related (and inverse of each other), but the test regression and the statistics are quite different in my opinion.


title: How many degrees of freedom are used in this model?

text: I have the table below, and I want to find the degrees of freedom for these data. Will my $k$ value from the formula $n-k-1$ include the dummy variables and `DIST^2`. That is, does $k=7$ or $3$, $4$, $5$, or $6$?               Dependent Variable: LFARE     Method: Least Squares     Date: 02/28/13 Time: 11:12     Sample: 1 4596     Included observations: 4596     Variable      Coefficient           Std. Error     C              4.770507             0.049427     LPASSEN        -0.080157            0.005603     BMKTSHR        0.301573             0.029724     DIST           0.000719             3.27E-05     DIST^2        -9.67E-08             1.21E-08     Y98            0.023120             0.013783     Y99            0.042357             0.013785     Y00            0.105950             0.013790          R-squared   0.428079     Adjusted R-squared 0.427206


title: The use of 'weights' in lmer to incorporate variance of individual data points?

text: This is a question regarding the use of the 'weights' function in lmer. Similar to the "weighted" approach discussed here for use in linear regressions: http://www.stat.cmu.edu/~cshalizi/uADA/13/lectures/ch07.pdf I wanted to use the "weights" function in lmer (using lme4 package) to account for variance of individual measurements in my dataset (ie.weights= 1/SEofY^2). That being said, I am not sure if this is correct and wanted to see if anyone had opinions on this. Thank you very much for your time.


title: Given more information, can a probability lessen?

text: Let $A$, $B$ and $C$ be events in the same probability space. Does $$\begin{align} \mathbb P(A\,|\,B\cup C) \ge \mathbb P(A\,|\,B) \end{align}$$ hold?


title: Orthogonal intersection in a Riemannian manifold

text: Let $S$ be the set of all probability distributions on $\mathbb{R}$ and $S_n=\\{p_\theta\\}$ be an $n$ dimensional submanifold of parameterized family of probability distributions on $\mathbb{R}$ where $\theta=(\theta_1,\dots, \theta_n)\in \Theta$ are called parameters of the probability density $p_\theta(x)$ and $\Theta$ is an open set homeomorphic to $\mathbb{R}^n$. It is also assumed that $p_\theta(x)>0$ for all $x$. For example, $$p_{\theta=(\mu,\sigma)}(x)=\frac{1}{\sqrt{2\pi}\sigma}e^{\frac{-(x-\mu)^2}{2\sigma^2}}$$ where $\mu\in \mathbb{R},\sigma>0$ forms a $2$-dimensional submanifold. On this manifold, a Riemannian metric is defined by the following, called Fisher information metric. \begin{eqnarray*} g_{i,j}(\partial_i,\partial_j) &=&\int \partial_i \log p_\theta(x) ~\partial_j \log p_\theta(x) ~p_\theta(x)~ dx\\\ &=& \int \partial_i p_\theta(x) ~\partial_j \log p_\theta(x)~ dx \end{eqnarray*} where $\partial_i=\frac{\partial}{\partial \theta_i}$. My question is regarding a point in page 384 from this Ann. Stat. paper by Amari. In the above mentioned paper, when he says a curve $\\{q_t(x)\\}$ in $S$, where $q_0(x)=p_\theta(x)$ intersects orthogonally $S_n$ at $q_0(x)=p_\theta(x)$, he means $$\int \partial_i \log p_\theta(x) ~\frac{d}{dt} \log q_t(x)\lvert_{t=0} ~q_0(x)~dx=0$$. i.e., $$\int \partial_i p_\theta(x) ~\frac{d}{dt} \log q_t(x)\lvert_{t=0} ~dx=0$$. Why? Could somebody make me understand by giving intuitive explanation of these concepts?


title: Good books/papers on credit scoring

text: I'm looking for recomendations of books on credit scoring. I'm interested in all aspects of this problem, but mostly in: 1) Good features. How to build them? Which have been proved to be good? 2) Neural networks. Their application to credit scoring problem. 3) I've chosen neural networks, but I'm interested in other methods as well.


title: Post-PCA analysis phase

text: Using PCA analysis, I was able to reduce the initial 23 variables into 10 principal components. But I do not understand what to do with this insight. I mean, how do I validate this information on, say, another set of data? Is there any way to identify which of the original set of 23 variables are more important, so that in the future I do not have to measure those variables and hence reduce the data overload?


title: How to identify a bimodal distribution?

text: I understand that once we plot the values as a chart, we can identify a bimodal distribution by observing the twin-peaks, but how does one find it programmatically? (I am looking for an algorithm.)


title: How do I parameterize a Weibull distribution in JAGS / BUGS?

text: Based on the answer to a previous question, For which distributions are the parameterizations in BUGS and R different? I have been transforming R parameterizations to JAGS parameterizations, but I have been getting errors so I am asking a separate question to clarify that the transformation is correct. R's `?dweibull` states the pdf (with $a = \text{shape}$ and $b = \text{scale}$): $$(^a/_b)(^x/_b)^{a-1} \text{exp}(- (^x/_b)^a)$$ And the JAGS manual states the pdf as: $$\nu \lambda x^{\nu-1}\text{exp}(-\lambda x^\nu)$$ I can see that the parameterization used by R has $f(b)^a$ but the JAGS parameterization does not have an equivalent $f(\lambda)^\nu$. But I can't bring myself to pull out a pencil when I have all of this computational power at my fingertips. So I did the following empirical demonstration that the JAGS parameterization is not to simply transform shape to rate as $\text{rate} =1/b$: The following plot represents two samples (code below),   * R (red) from $\small{Y\sim\text{Weibull}(a = 2, b = 50)}$    * JAGS (blue) from $\small{Y\sim\text{Weibull}(\lambda = 2, \nu = 1/50)}$  ![enter image description here](http://i.stack.imgur.com/CiPMy.png)               library(rjags)     shape <- 2     rate  <- 50     set.seed(0)     model.string <-      writeLines(paste("\nmodel{\nbeta ~ dweib (2,",1/rate,")\nY <- beta\n}"), con = 'weibulltest.bug')     j.model  <- jags.model(file = "weibulltest.bug", data = list(x=NA)) #hack     mcmc.object <- coda.samples(model = j.model, variable.names = c('Y'), n.iter = 10000)     Y.jags <- as.matrix(mcmc.object)     Y.r <- rweibull(10000, 2, 20)     plot(density(Y.r), col = 'red', ylim = c(0,0.15))     lines(density(Y.jags), col = 'blue')


title: Should a predictor, significant on its own but not with other predictors, be included in an overall multinomial logistic regression?

text: I constructed a model via multinominal logistic regression analysis. The final model contains three predictors. All predictors are significant when they are the only predictors. However, the coefficient of one of the predictors is not significant when included all three predictors are included in the model. Should I include this predictor in the final multinomial logistic regression equation?


title: detector confidence

text: the situation is as follows: Say, I am going to answer a multiple choice question based on suggestion of my two friends. One gives correct answer 60% percent of the times, the other gives correct answer 40% of the times. Is there any way I can answer a question correctly more than 60% of the time? (in a long run average leaving aside natural variability)


title: What does the "different populations" result of a significant group t-test mean exactly?

text: Why does finding a significant difference on a given measure, say mathematical ability, between two groups (e.g., via an independent samples t-test) have to be equivalent to the two groups coming from different underlying populations - can "population" really be defined as simplistically as "the population of people good at math" and "the population of people rubbish at math"? And assuming the two groups tested represent, say, the two halves of a classroom, what does it mean that there's a significant difference between the groups? That the populations from which each of the two halves were sampled are different in terms of how maths scores are distributed within that population, or different in other ways?


title: Interpretation of dummies (several variables and categories!)

text: I am trying to figure out how to interpret the coefficients when using several qualitative independent variables of which some have more than two categories. I found this short article that touches upon that issue (the one page preview is enough to understand what I mean but I will explain here too). The article says that the interpretation is messed up in this case. Let me use the example variables of the article to make it clear what I mean. Wage is my dependent variable. The independent variables are all dummies that can be grouped in: _educational attainment_ E1 = postgraduate (value 1 if postgraduate, 0 otherwise) E2 = bachelor (value 1 if bachelor's degree, 0 otherwise) E3 = high school (value 1 if high school, 0 otherwise) _marital status_ R1 = married/in relationship (value 1 if married, 0 otherwise) R2 = divorced (value 1 if divorced, 0 otherwise) R3 = single (value 1 if single, 0 otherwise) _sex_ M = male (value 1 if male, 0 otherwise) F = female (value 1 if female, 0 otherwise) Now, to avoid the dummy trap I delete E3, R3 and F. My regression equation would then look like this: wage = b0 + b1*E1 + b2*E2 + b3*R1 + b4*R2 + b5*M The article states that the baseline to compare against would thus be a single female that went to high school and that would mess up interpreting the individual dummy coefficients. My questions:   1. Do I understand it correctly that it is indeed using several qualitative variables that have more than 2 categories (in this case educational attainment and marital status have each 3) that messes up the interpretation? I get that if I for example only used the marital status dummies (again excluding R3 to not face the dummy trap), I could easily state that married people (R1) earn more or less compared to single people (R3). But is this no longer possible as soon as I include the educational attainment dummies in my regression as well?   2. Are dummy variables that only have 2 categories and therefore are represented in the regression equation with 1 dummy only (e.g. my sex variable -> M in the equation) affected by this at all or could I still state that men earn more or less compared to females, despite having educational attainment and marital status in my equation? Thank you for your time, I greatly appreciate any feedback.


title: Outliers and Linearity for EFA

text: I am having 56 likert scaled items for IV and 28 likert scaled items for dv. As to fulfill the assumption for EFA, outliers and linearity need to be checked. Can anyone tells me what method/analysis is the best for identifying outliers and linearity? For outliers, can I use Z-scores?


title: Viterbi training vs Baum-Welch algorithm

text: I'm trying to find the most probable path (i.e. sequence of states) on an HMM using the Viterbi algorithm. However, I don't know the transition and emission matrices, which I need to estimate from the observations (data). To estimate these matrices, which algorithm should I use: the Baum-Welch algorithm or the Viterbi Training algorithm? Why? In case I should use the Viterbi training algorithm, can anyone provide me a good pseudocode (it's not easy to find)?


title: PCA for large matrix that don't fit in memory

text: I'm looking for some matrix library that can do PCA and can work with large matrix that can't whole fit in memory at once. is there exist such library ?


title: How was this deviance formula for Cox Proportional Hazard derived?

text: The gbm package in R documents the following deviance formula for CoxPH models http://gradientboostedmodels.googlecode.com/git/gbm/inst/doc/gbm.pdf $-2 \sum w_i \delta_i (f(x_i) - \log(R_i / w_i)$, where $w_i$ is a weight for the observation, $\delta_i$ is the right censor indicator and $R_i = \sum_{j=1}^N 1(t_j >= t_i) e^{f(x_i)}$ is the total hazard of the risk set However this formula does not seem to agree with the author's implementation. https://github.com/harrysouthworth/gbm/issues/22 I'm wondering how this was derived so I can figure out if there is a bug in the implementation or if there is a mistake with the way the formula was written


title: Random Balls in Random Buckets: What are the characteristics of the distribution?

text: I have N buckets, numbered 1 to N. I draw k random integers, uniformly distributed in the range 1 to N, with replacement, and for each integer I drop a ball into the corresponding bucket. k can be any size; specifically it can be any value from 2 to N, or larger than N.   * What are the statistical characteristics (probability distribution, etc) of the numbers of balls in each bucket at the end?    * What are the statistical characteristics of the numbers of buckets with 0,1,...k balls? This question arises from a need to measure the 'goodness' (in some sense) of a hashing algorithm. Given a sample of distribution of keys among buckets, I need a measure of where it lies on a scale from 'Very good' to 'Very bad', and to be able to work out things like 'what are the chances of more than x balls in a bucket given k and N?' Obviously, from the way I've written this, and from the random names I've assigned to my variables, I have absolutely no statistical sophistication. Please be gentle; I want to learn. Please feel free, for example, to change the variable names to something more conventional, or anything else.


title: Which data mining packages support anomaly detection?

text: I aim to have some anomaly detection process on my data but Weka, Rapidminer or Knime do not support anomaly detection algorithms. How would I take care of the process?


title: Online resources to aid in understanding SAS results

text: Every time I look at SAS's result page I get lost among the numbers. Is there any online material that can help with understanding SAS's result page?


title: How should I use prop.test function?

text: I did an experiment with 12 participants. They were asked to rate 3 stimuli, with the following 3 possible answers: "increased", "decreased", "no difference". I have as a result this table:                         Increase Decreased No_difference           StimulusA   10      1           1     StimulusB   12      0           0     StimulusC    8      2           2      Now,if I want to understand in one shot if the "Increase" answer is chosen significantly more than the other two answers in all the Stimuli, is it correct if I simply use the binomial distribution summing the successes along all the 3 stimuli? I mean, is it correct doing this?               > prop.test(30,36)              1-sample proportions test with continuity correction          data:  30 out of 36, null probability 0.5      X-squared = 14.6944, df = 1, p-value = 0.0001264     alternative hypothesis: true p is not equal to 0.5      95 percent confidence interval:      0.6652978 0.9303666      sample estimates:         p      0.8333333       Or it is better to use the function in this way instead?               > prop.test(c(10,12,8),c(12,12,12),c(0.5,0.5,0.5))              3-sample test for given proportions without continuity correction          data:  c(10, 12, 8) out of c(12, 12, 12), null probabilities c(0.5, 0.5, 0.5)      X-squared = 18.6667, df = 3, p-value = 0.0003204     alternative hypothesis: two.sided      null values:     prop 1 prop 2 prop 3         0.5    0.5    0.5      sample estimates:        prop 1    prop 2    prop 3      0.8333333 1.0000000 0.6666667


title: Local polynomial (linear) regression of binary data --- logit transformation?

text: I got a bit confused about how to fit a local polynomial to binary outcomes if I would rather approximate the underlying index (within a link function) instead. (Basically for the same reason why someone would estimate a logit or a probit instead of a linear probability model.) Of course, I would like to plug the variable in the inverse of the link function --- but of course, 0 and 1 will give me $-\infty$ and $+\infty$ in common link functions. To clarify, my problem is to flexibly predict $y(x)$, where $x \in [-1,1]$ but $y \in \\{0,1\\}$. I was not thinking about splines, as I would be more interested in a smooth prediction of the index $f(x)$ if the prediction is $\hat{y}=g(f(x))$ with a link function g. I do not have a prior on the order of the splines and the location of the breakpoints. I would be interested in general advice on the appropriate approach. As a shorthand, I substituted 0 with $\frac{1}{N}$ and 1 with $\frac{N-1}{N}$, and generated the index with these. (Using a logit link, thus I was smoothing/predicting $f(x)=\Lambda(\frac{1}{N})$ if y(x)=0, e.g.) But I am uneasy about this. I coded up leave-out cross-validation to pick the right bandwidth for smoothing f with the short-hand f, and even used the delta method to get back to confidence intervals of the original binary variable from my smoothing of the index. But that might be fake sophistication if I approach the prediction wrong in the first place. Thanks for any thoughts on this!


title: Problems with extremum of two uniform random variables

text: Here is the problem from the book: > Let $X = \min(U,V)$ and $Y = \max(U,V)$ for independent > $\text{uniform}(0,1)$ variables $U$ and $V$. Find the distributions of a) > $X$; b) $1-Y$; c) $Y-X$. I know that each has density $2(1-x)$ for $0<x<1$. What I'm hoping for is an intuitive explanation of how to understand this problem and arrive at that solution without knowing the beta distribution.


title: What is the best way to treat popularity for ranking?

text: I have items for purchase and I have statistics how much each were purchased. I want to use this amount for ranking them. The simple guess is to `Log(2 + num_purchased)`. I thinking about better formula for that, and most important I want to have some probabilistic background and meaning behind it.


title: Is $f(Y | X)$ in the same family as $f(Y)$?

text: Is it the case that for random variables $Y$ and $X$, $f(Y | X)$ in the same family as $f(Y)$? If so how can I prove it? If not are there any situations (some families of distributions) where they can be the same? Or is it related to the link function(s)? Thanks


title: Cross validation for variable selection and coefficient shrinkage?

text: Is cross validation an appropriate technique for variable selection and regression coefficient shrinkage? A former colleague of mine used 10-fold CV to compare the regression coefficients from the 10 training models. Variables were kicked out if the coefficients changed sign. The average coefficient values for the remaining variables was then used to estimate the coefficients of the final model. Does this homebrewed methodology sound like a valid technique?


title: Establishing a dynamic threshold for partitioning data

text: Given a sorted set of data points, where each point represents the length of a common substring between two files, I want to systematically decide a cut off point where I can start considering them to be "short" and no longer care about them. Basically my goal is to exclude the data points that represent commonalities between the two files that exist because we’re looking at such small substrings the probability that they will repeat is large, rather than because they are actual common substrings. My attempt to solve this was to first make a simple histogram using the Freedman-Diaconis rule for choosing a bin width. An example here shows the bucket number on the left and the frequency on the right.               0  : 68     1  : 14     2  : 5     3  : 6     4  : 1     5  : 0     6  : 1     7  : 1     8  : 1     9  : 1     10 : 2     11 : 2      You can see from above that the first couple buckets are littered with these small commonalities I am wishing to avoid. I don't want to have a "hard coded" cut off point but would like to decide this for each data set. Any ideas or pointers based on what I mentioned as far to solve this? Also I would like some input into my idea of using a histogram to solve this. If any other methods are better/more practical or any improvements to what I wrote I would love to hear.


title: Interpreting logit regression with both continuous and categorical variables

text: I have a logistic regression like this:               Y = a1 + b1*(number positive scores) + b2*(number negative scores) + b3*Z             + b4*Z*(number positive scores) + b5*Z*(number negative scores)             + additional non-interaction terms      `Y` is the probability of an outcome that takes on binary values (0,1). `Z` is a continuous variable; I am trying to determine if it is significant Number of positive and negative scores take on integer values between (0,5) When I run the regression using `glm` in `matlab`, I find that all coefficients are significant except for b3; b1 and b4 are positive whereas b2 and b5 are negative. I would like to draw conclusions about whether Z is in fact a significant factor in the outcome variable via interactions b4 and b5, but I understand that in a logistic regression all coefficients and particularly interactions need to be evaluated in the context of specific values of the independent variable x. So this is what I have done so far: Let's say bhat is the estimated vector of coefficients `b=(a1,b1,b2,b3,b4,b5...)` and `xhat` is the sample mean of `x=(num pos scores, num neg scores,Z,Z*num pos scores, Z*num neg scores)`. Also say that bi_sigma is the estimated standard error on the ith element of bhat and `xi_sigma` is the sample standard deviation of the ith element of x. Let's say `L(bx)` is the logistic cdf evaluated at bx. Am I right to evaluate the odds ratios `exp(bhat*xhat-bi_sigma*xi_hat)` and `exp(bhat*xhat+bi_sigma*xi_hat)`, then determine whether each element bi_hat is significant based on whether the range of these odds ratios is strictly greater than or less than one? In other words, my thinking is that if the odds ratios don't include one then they significantly improve or decrease the odds. For instance an odds ratio range of (1.3,2) reflects a bi that improves the odds for mean levels of x. Yes? Secondly, am I correct to evaluate the first differences, `L(bhat*xhat)-L(bhat*xhat-bi*xi_sigma)` and `L(bhat*xhat+bi*xi_sigma)-L(b*xhat)` as two measures of the size of the impact of each `xi`. Can I do this for the interaction also? Thanks very much for any advice or help. If there are references that you would suggest for this, I would appreciate that too. All the examples I've found online involve binary or categorical variables only, and no continuous variables.


title: Estimating non-stationary Markov chain transition probabilities from data

text: I have some call centre outcome data: each person is called several times, and the result of the call is recorded as one of a few discrete outcomes ("no answer", "wrong number", "third party", "correct party", etc.). The goal of the analysis is to try and work out how many times to call before giving up. My current plan is to consider the outcomes as a Markov chain. If I assume that the data represents a stationary state, then it is easy to get the transition probabilities. The problem is, I don't believe that they are stationary: having "no answer" 20 times is a different situation to be in than having "no answer" once. How do you calculate the transition probabilities if your data is non- stationary? For bonus points, is there a function in R to do this?


title: I have 3 dimension reduction methods and 3 classifiers, how to select the best combination?

text: I am trying to classify two types of objects, which unfortunately have high- dimensional features with few samples. (230 features 12 samples from each group). As a first step: To reduce the dimension, I have tested three different approaches of PCA (each with slightly different parameters) and use the scores of the first 3 PCs as the features instead of the original 230 features. As a second step: in order to classify those objects, I have tested three different classifiers (SVM, 3-NN, Naive Bayes) where I use the leave one out method, i.e.: training the classifier on 23 objects and tested it on the one leave out, for testing. Sum up: I have use 3 different dimensional reduction methods and three use different classifiers (nine combinations overall). In one combination (out of 9), I got excellent classification performance. My question: Is there some statistic approach I can use here in order to prove that this option was indeed robust and the correct one, and didn't happen by chance?


title: What is the equivalent confidence interval for a directional hypotheses?

text: Testing the null hypothesis $H_0: \theta = 0$ against the alternative $H_1: \theta \neq 0$ is equivalent to constructing a confidence interval and checking whether or not zero belongs to it. But what about the alternative $H_1: \theta > 0$. What is the "equivalent confidence interval?".


title: NAs in summary of linear model generated from significant genes for continuous response variable identified using SAM

text: I have a gene expression data-set with log2-transformed expression values (no NAs) for 495 genes for 59 samples for which values of a continuous response variable (r) are also known (no NAs). I want to use leave-one-out cross validation to test if r of the test sample can be predicted from the sample's gene expression. For this, I intend to use the samr R package for Significance Analysis of Microarrays to identify significant genes associated with r in the training set of samples. Then, I want to generate a linear model using the significant genes as variables, which will then be used to predict r of the test sample. I have tried the following code to begin with, but when I generate the model and examine it, I see many NAs in the model summary, which makes me suspect that I am doing something wrong. Can someone tell me what I might be doing wrong? Secondly, I will appreciate any comment on the use of nperms (in SAM) with a value of 100. Is it too low for an expression data-set for 495 genes.               # rVals with the r values is read as a vector from a row of a table for phenotypic data read from a tab-delimited file with sample-names as column names and phenotype features as row-names     # geneVals is the log2-transformed gene expression data-set read as a matrix from a tab-delimited file with sample-names as column names and gene-names as row-names          # Perform SAM with FDR of 5% and obtain list of significant genes          sam <- SAM(x=geneVals, y=rVals, resp.type=c("Quantitative"),     testStatistic=c("standard"), regression.method=c("standard"), logged2=TRUE,      fdr.output=0.05, eigengene.number=1, knn.neighbors=10, nperms=100,      genenames=as.vector(rownames(geneVals)))          sigGenes <- rbind(sam$siggenes.table$genes.up, sam$siggenes.table$genes.lo)          # Generate linear model     toModel <- data.frame(t(rbind(rVals, geneVals)), check.names=FALSE)     myModel <- lm(toModel[c('r', sigGenes[,c("Gene ID")])])          # Examine model     summary(myModel)          ...output...          Call:     lm(formula = toModel[c("rVals", sigGenes[, c("Gene ID")])])          Residuals:     ALL 59 residuals are 0: no residual degrees of freedom!          Coefficients: (58 not defined because of singularities)                Estimate Std. Error t value Pr(>|t|)     (Intercept)   -18.29363         NA      NA       NA     `let-7e`       -1.70545         NA      NA       NA     `miR-125a-5p`   2.43177         NA      NA       NA     `miR-151-5p`    2.67439         NA      NA       NA     ...


title: Probability of total events occurring given that one or more events occur in specified number of months

text: As a non-statistician, I have a real world statistical/probability problem that I'm having trouble framing. The software I rely on in inventory management interprets the 'movements' (number of times inventory is used) in a strange way. It offers the number of months, out of 24 months, that the item has been used. For example, a moving code of 20 means that the item was used at least once for 20 out of 24 months.. What I need to be able to do is translate that to find **the most probable number of movements over that 24 month period.** If movements randomly fall into 20 out of 24 months with no limitations, obviously the number of movements that really occur is likely to be much greater than 20. How much greater? Sorry, this question is extra challenging because I have no ideas on how to begin tackling this. Any help is much appreciated.


title: Is there a "pure R" implementation for loess? (with no C code?)

text: I would like to manipulate the work done by the loess function in R. However, the main workhorse of this function is written in C. Is there some pure R implementation of the code? Thanks.


title: Is the maximum entropy distribution consistent with given marginal distributions the product distribution of the marginals?

text: There are generally many joint distributions $P(X_1 = x_1, X_2 = x_2, ..., X_n = x_n)$ consistent with a known set marginal distributions $f_i(x_i) = P(X_i = x_i)$. Of these joint distributions, is the product formed by taking the product of the marginals $\prod_i f_i(x_i)$ the one with the highest entropy? I certainly believe this is true, but would really like to see a proof. I'm most interested in the case where all variables are discrete, but would also be interested in commentary about entropy relative to product measures in the continuous case.


title: Calculate the probability of being in spatial zone

text: I know that the probability of a point on a 2D map follows a normal distribution dnorm(x), where x is the distance to another point. Lets call it A. so x is calculated by a distance function               A = c(2,1)     dist = function(mypoint){        sqrt((mypoint[1]-A[1])^2+(mypoint[2]-A[2])^2)     }      I want to calculate the probability to be in a certain zone of the map. For example $X=(x,y)$ and $S\in[S_{x1},S_{x2}]*[S_{y1},S_{y2}] $ $P(X\in S) = P(S_{x1}<X_x<S_{x2} \cap S_{y1}<X_y<S_{y2})$ $P(X\in S) = \int_{min(S_{x})}^{max(S_{x})}\int_{min(S_{y})}^{max(S_{y})} dnorm(dist([x,y]))dx.dy$ Where $dnorm$ is the normal density. I have two questions:   1. Does this make any sens? Is there another (better) way to acheve this?   2. How can I calculate such a probability with R?


title: households survey data - determine 'valid' product categories

text: I'm working with household survey data where households report their expenditure for various products. The goal is to look into income elasticities for various products thus in essence I'm interested in households' expenditure given their incomes. The problem is to find a consistent approach in determining which product categories are invalid for this type of analysis. By invalid I mean the following: suppose we have two products - milk and lobsters. It would be quite reasonable to assume that as households consume much more milk than lobsters, the given sample will be more likely to depict the households' behaviour towards their purchases of milk rather than lobsters throughout the income spectrum. I could try to look at for example variation in expenditure among households for a given product by analysing coefficients of variation or somehow look at variance of the dependent variable in linear regression like `expenditure ~ income` or even more simply I could look at the sample sizes for each product category (that is, how many households reported buying a particular good). However, I wonder if there is a consistent statistical approach to this kind of problem in determining the 'validity' of product categories? Any suggestions or references would be very helpful. Thanks!


title: Multiple Correspondence Analysis percentages, and points near to one axis

text: I have done a multiple correspondence analysis based on 10 factors: every factor is a `Yes`or `No`question. The first factor give 22.3% of the total variance and the second 13.6%. I am used to Principal Component Analysis (PCA) and to having a high percentage (>70%). Question1: Is a total percentage of 35% (22+13) enough for an MCA? Question2: base on the image attached, is there any implication when the variables are close to the second axis? Any other comment about the MCA I have done is appreciated. ![enter image description here](http://i.stack.imgur.com/SAlzm.png)


title: Measuring Cronbach Alpha for very large number of items

text: How can I calculate Cronbach Alpha for two tests that I have administered to 165 people? One test contains 110 items and the other contains 67 items. The scoring is done on the hard copies and the scores of each subscale rather than individual items are entered for further calculations. There are 11 sub-scales in the one with 110 items, and four sub-scales in the one with 67 items. Presumably, the subscales are uni-dimensional in construct, so should I be calculating alpha for them. I mean do the results become invalid if I do not find alpha for my data. Is there any other more practical way of ensuring the internal consistency, or do the original alpha values given by the developers of the scales suffice? The task of putting scores of 165 people for 177 items is not practical. My tables show scores of each subscale (not items) for all the participants. Please help me understand what is the best way to go about it.


title: Getting unexpected lme results

text: I am evaluating an experiment in which subjects had to rate quality of audio files. The design of the experiment is repeated-measures. So I'm trying to figure out which factors have an impact on quality. I am using the `lme` function in order to control within subject variance. `amplitude` and `number` are my two factors to predict quality.               lme(quality ~ amplitude*number, data = some_data, random = ~ 1 | subject)      This is the result:                Correlation:                       (Intr) ampltd number     amplitude        -0.246                   number           -0.234  0.331            amplitude:number  0.165 -0.704 -0.666      That's not what I want. If I change the forumla to `quality ~ amplitude*number` I get                Correlation:              (Intr)     quality -0.845      which is what I wanted. I am confused as it seems that the values to predict are on the right side of the formula but I thought that I have to specify the formula like `y ~ factor1*factor2` but obviously `factor1*factor2 ~ y` is correct. What am I getting wrong here?


title: What type of test and how many participants?

text: I'm looking at collecting information on whether participants playing an educational game on a particular topic will achieve lower/comparable/superior gains between pretest/post-test scores compared to participants who learn that subject matter without a game (lecture only). Thus, I intend to randomly assign participants to one of two groups: game or no game. Each participant will take a pretest and post-test. Assume the results are as follows: all participants from both groups score an average of 50% on the pretests, the participants who play the game have an average post-test score of 90%, and the participants who don't play the game (lecture only) have an average post-test score of 70%. With the correct statistical technique and "enough" participants, it might suggest that the game made a difference. Additionally, I also plan on examining whether gender affects pretest/post- test scores (thus, should I have four groups? Male/game, Male/no game, Female/game, Female/no game) I'll likely also look at other demographics like age (categories to be determined, but assume: age 20-40/game, 20-40/no game, 41-60/game, etc.) Unfortunately, it's been over 20 years since I had a stats course. What technique do I use for these types of calculations (ANOVA?), and how many participants do I need in each category for the results to be valid? Additionally, what confidence interval is normally considered "valid"? 95%? 99%?


title: Compute a confidence interval for Bernoulli distribution

text: Given a Bernoulli distribution with a success probability of p = 0.02. Let's say we have N=3000 samples. How can I compute a confidence intervals for the expected number of successes (e.g., with 5% significance level)?


title: Quantile regression: Which standard errors?

text: The `summary.rq` function from the quantreg vignette provides a multitude of choices for standard error estimates of quantile regression coefficients. What are the special scenarios where each of these becomes optimal/desirable?   * "rank" which produces confidence intervals for the estimated parameters by inverting a rank test as described in Koenker (1994). The default option assumes that the errors are iid, while the option iid = FALSE implements the proposal of Koenker Machado (1999). See the documentation for rq.fit.br for additional arguments.   * "iid" which presumes that the errors are iid and computes an estimate of the asymptotic covariance matrix as in KB(1978).   * "nid" which presumes local (in tau) linearity (in x) of the the conditional quantile functions and computes a Huber sandwich estimate using a local estimate of the sparsity.   * "ker" which uses a kernel estimate of the sandwich as proposed by Powell(1990).   * "boot" which implements one of several possible bootstrapping alternatives for estimating standard errors. I have read at least 20 empirical papers where this is applied either in the time-series or the cross-sectional dimension and haven't seen a mention of standard error choice.


title: Finding Max, Min, and Joint Density of max and min

text: X, Y are two independent R.V.s uniformly distributed over the interval [0,1] 1)Z = Max[X,Y]. Find the C.D.F. and density 2)W = Min[X,Y]. Find the C.D.F and density 3) Verify that the joint density of (W,Z) is 2 for 0 I thought I knew how to do it, but since the joint density equals 2, I don't think I did it right. For 1, I got x^2 as the C.D.F. and 2x for the P.D.F.


title: About Covariance

text: I'm coding in python for a project. I have all my data and run an "optimize.leastsq". I want the uncertainty in cov[1,1] and cov[0,0]. Is it just the square root of these matrix elements or is there more involved? Picture attached below ![enter image description here](http://i.stack.imgur.com/VHSxj.png) ![enter image description here](http://i.stack.imgur.com/pChXM.png)


title: Probability of game type lottery outcomes

text: This is about a lottery type football. It consists of 14 games. 13 of them have 3 possible outcomes: home team wins, tie, or away team wins. The 14th game is about 1 of the 13 games in which one chooses which team scores first (or no score), and in what quarter. So it can be home team scores in 1q, 2q, 3q or 4q, and the same for the away team, or no score at all. So in this game there are 9 possible outcomes. I already calculated the probability of winning 1 game, 2 games, 3 games, etc. using the binomial distribution. I want to calculate the same but with an alternative way of play. You can play the lottery by choosing doubles and triples (this option is only valid for the first 13 games). Doubles allow to choose 2 of the 3 outcomes, and triples allow you to choose all three, securing that result. Questions. If I use, let's say, 3 doubles and 1 triple: How many combinations would I be using? and, what is the probability of choosing correctly 0 games, 1 game, 2 games, etc.? (the 2 questions include choosing 1 outcome of the 14th game). Another part of the problem. I have the data of the past 500 lotteries (only for the first 13 games) and did a frecuency table of each outcome, calculated means and standard deviations. I want to know the probability of the outcomes based on their frecuencies. I did a normal distribution but I don't know if it's the correct way, because it doesn't add up exactly. For example, it shows a very small probability of 15 home wins, and this should be zero because there are only 13 games. And also the probabilities of each outcome don't add up to 1, it's very close, but it should be exactly 1. Based on the frecuencies, the most common results are 5 home wins, 5 ties, and 3 away wins. Final question. How many combinations can be made with that array? Also including 1 outcome of the 14th game.


title: Estimate center and radius of a sphere from points on the surface

text: If we assume that our data points were sampled from the surface of a sphere (with some perturbation), how can we recover the center of that sphere? In my searching, I found papers on something labeled "spherical regression", but it didn't quite seem like that was doing the same thing. Maybe I just didn't understand it. Is there a straightforward formula, similar to linear regression, that finds a sphere center point and radius that minimize the sum-squared distance of a set of data points from the surface of the sphere? * * * Edit 1: We can assume that the noise will be 2 or 3 orders of magnitude smaller than the radius of the sphere and uniformly spherically Gaussian. However, the samples themselves will definitely not be drawn uniformly from the sphere's surface, but will likely be clustered in a few patches on the surface, likely all within one hemisphere. A solution that works for data in $\mathbb R^3$ is fine, but a general solution for arbitrary dimensionality is great too. * * * Edit 2: What are the chances that I might get a sensible answer if I were to use linear regression, $y = X\beta + \epsilon$, in the 7 dimensional space pretending that the squared components are independent from the other parameters: $\begin{align} X &= \begin{array}{ccccccc}[-2x& -2y&-2z&1&1&1&-1]\end{array}\\\ \beta &= \begin{array}{ccccccc}[x_0 & y_0 & z_0 & x_0^2 & y_0^2 & z_0^2 & r^2]'\end{array}\\\ y &= x^2+y^2+z^2\end{align}$ At best, I suppose that my error metric will be a bit wacky. At worst the solution won't be even close to consistent.   ...or that's silly because with four identical columns, we get a singular matrix when we try to do regression. * * * Edit 3: So, it looks like these are my options:   1. Non-linear numerical optimization using some cost function: $f(x_0,y_0,z_0,r|X) = \frac{1}{2}\sum_{i=1}^n \left(r - \sqrt{(x_i-x_0)^2+(y_i-y_0)^2+(z_i-z_0)^2}\right)^2$   2. Hough-transform: discretize the plausible space or possible centers and radii around the data points. Each point casts a vote for the potential centers that it could be part of at each specific radius discretization. Most votes wins. This might be okay if there were potentially an unknown number of spheres, but with just one it's a messy solution.   3. Randomly (or systematically) select groups of 4 points and analytically compute the center. Reject the sampling if ill-conditioned (points are nearly co-planar). Reject outliers and find the mean center. From that we can find the mean radius. Does anyone have a better method?


title: Differences between MANOVA and Repeated Measures ANOVA?

text: * What is the difference between a repeated measures ANOVA over some factor (say experimental condition) and a MANOVA?   * In particular one website I stumbled across suggested that MANOVA does not make the same assumption of sphericity that repeated measures ANOVA does, is that true?       * If so, why would one not just always use MANOVA?    * I am trying to conduct a repeated measures ANOVA with multiple DVs, what is the appropriate approach? Thank you.


title: MANOVA or repeated-measures ANOVA?

text: thank you very much for your time! I am currently working on a neuroimaging study comparing lateralization indexes between the two hemispheres in three different groups. I have 4 lateralisation indexes (two different independent parameters for two different fiber tracts) and am now wondering whether it's better to use MANOVA or repeated-measures ANOVA. In the literature I found some studies looking for differences in diverse lateralization indexes via repeated measures ANOVA but I am not quite sure if this is the correct approach since the two indexes are independent and are not describing the same fiber tract. Thank you very much for your help in this matter!


title: curve fitting with extremely "large" numbers

text: here is an unconventional but apparently workable idea. wondering if anyone has tried something like it, esp looking for references, examples, or nearby related work. am working on a mathematical induction problem with a sequence $x_i$ of strings which are "generally increasing" in size as $i$ goes to infinity. from inspection the function is loosely related to the Fibonacci sequence. these strings could be converted to very large multidigit base-N numbers (ie with something in the range of dozens of digits or more) and one could do a curve fit using standard algorithms eg Marquardt-Levenburg nonlinear least squares algorithm, powells algorithm, etc. of course its preferable to get an _exact_ curve fit such that it matches $x_i$ for some finite $i \leq k$ and then generalizes to unseen data $i>k$. ie _induction_. > so what are some of the largest magnitude numbers that curve fitting has > been applied to, and esp cases where exact formulas have been obtained? are > there some researchers, refs, or some software that specializes in this? etc my question should not be taken to be asking about datasets with many datapoints. this question is a much different regime of interest where there may be "not so many datapoints", but each datapoint has a very large numerical magnitude with lots of digits. so its "big data" in a sense but with a key twist. note however equivalently this question _is_ the same as fitting floating point values to very high degrees of precision and the tight fit is not arbitrary but instead meaningful to the problem in some way. (this question is partly inspired by the AMS ref on _experimental mathematics_.[1]) [1] Experimental Mathematics: Examples, Methods and Implications by Bailey/Borwein, notes of the AMS, 2005


title: How to find the closest distribution of a given data?

text: I have inter-arrival times of vehicles recorded by a vehicle detection algorithm. I want to find the closest distribution (e.g., Poisson or other) of this data. How can I do that? Here is a graph of the inter-arrival times from a crosspost on SO. ![enter image description here](http://i.stack.imgur.com/gxQwj.jpg)


What are the meanings of "stochastically dependent" and "functionally dependent"? What is the difference? (I saw the usage of the above terms in the paper "Maximum entropy sampling and optimal Bayesian experimental design") Thanks a lot.


title: Is AdaBoost less or more prone to overfitting?

text: I have read various (seemingly) contradicting statements whether or not AdaBoost (or other boosting techniques) are less or more prone to overfitting compared to other learning methods. Are there good reasons to believe one or the other? If it depends, what does it depend on? What are the reasons that AdaBoost is less/more prone to overfitting?


title: DLM package, State equation

text: I would like to know how can I write in dlm form ( dlm PAckage by Petris) a state equation like this: $$a(t+1) = k(t) + F*a(t) + T(t)w(t)$$ where $w$ is $N(0,1)$. To make it easier I can assume $k$ and $T$ constant. My question is how can I rearrange it in the way described by "An R Package for Dynamic Linear Models"? Thanks


title: Significant improvement in education - can I show it?

text: I have a data set of observations that show an observed effect following education. I have a baseline position (taken prior to the education) which shows the number of students showing particular knowledge and then post- education an increased number of students showing particular knowledge. How do I show that the gain is significant? Unfortunately sample size is relatively small (60 students). Also I have data on the genders of the students. Looking at the data I have I would say the impact of gender is not relevant but how do I show that is true? Thank you!


title: Second moment of draws from a multivariate normal covariance matrix

text: Suppose I have an $N\times N$ covariance matrix that describes a multivariate normal joint distribution. Now take 100,000 draws of the covariance matrix. I measure the variance of values for each draw -- call this the "cross sectional dispersion". If I average the result of cross-sectional dispersion across draws what is the analytical expression for the mean and variance of this dispersion? Here is some code if you'd like to toy with this or test your idea:               # create 3x3 covariance matrix     m1 <- c( .1 , .15 , 0 )     S1 <- matrix( c( .1 , .05 , .02 , .05 , .1 , .03 , .02 , .03 , .1 ), nrow = 3 )          # draw 10,000 times from covariance matrix     set.seed(1000)        library(MASS)     draws <- MASS::mvrnorm(10000,mu=m1,Sigma=S1)          # measure the "cross-sectional" dispersion at each row     rowDispersion <- apply( draws , 1 , var )          # measure the mean cross-sectional dispersion     mean( rowDispersion )          # measure the variance of cross-sectional dispersion     var( rowDispersion )


title: Two negative main effects yet positive interaction effect?

text: I have two main effects, V1 and V2. The effects of V1 and V2 on the response variables are negative. However, for some reason I am getting positive coefficient for the interaction term V1*V2. How can I interpret this? is such situation possible? thank you


title: How to calculate the conditional variance for AR(2) model?

text: I need to find the conditional variance of $Y_t$ given information up to time t. $Y_t$ = $\mu $ + $\phi_1$$Y_{t-1}$ + $\phi_2$$Y_{t-2}$ +$\epsilon_{t+1}$ Need to find the conditional variance of this equation. I need to find three step ahead forecast. Was able to do the first, but got stuck at second step. Var($Y_{t+1}$|$\Omega_t)$ = $\sigma^2$ But need to find Var($Y_{t+2}$|$\Omega_t)$? Need some guidance.


title: Limiting Distribution of $W_n=\frac{Z_n}{n^2}$ , $Z_n \sim \chi ^2 (n)$

text: My try ended in an awkward result. I thought it best to use the moment generating function (MGF) technique. We can derive the MGF of $W_n$ as follows: $$ E \left[ e^{tZ /n^2} \right]= \left(1-\frac{2t}{n^2} \right)^{-n/2}$$ from the chi-squared MGF. But the problem is that the limit of that as $n \to \infty$ leaves $1$ and I am left puzzled whether I did everything right. Have I missed something? Thanks.


title: Sort lists using statistics

text: I want to make a ranking of music genres. I have a list of artists and top genres with a value (1 - 100) which describes how much their music is from that genre. Now I need to organise the genres so the most popular is on top. Bellow I have genres from 4 artists (3 genres per each). We can see that `Rock` and `Metal` are both in 3 groups, but when looking at the value it is clear that `Rock` is more popular than `Metal`. I was thinking about summing the values from the genres, and order the by the summ, but I think in statistics it must be an algorithm that will give better result when doing the rank. The question is how to order the genres (btw, I'm using PHP). * * *   * metal (71)   * industrial (65)   * rock (36) * * *   * rock (100)   * alternative (67)   * metal (36) * * *   * metal (100)   * Melodic Death Metal (72)   * modern metal (68) * * *   * alternative rock (100)   * rock (78)   * alternative (63) * * *


title: Interpretation of log likelihood and covariate significance in Cox regression

text: Another newbie question here (probably piece of cake for you guys). When I run a Cox Regression and one of my covariates come out as significant:                                                                   95,0% CI for Exp(B)             B       SE      Wald    df  Sig.    Exp(B)  Lower   Upper     risk    2,224   1,107   4,036   1   ,045    9,244   1,056   80,950     <5 covariates removed for clarity>      But my Omnibus Tests of Model Coefficients(a) isn't significant. Change from previous step: Chi-square=10,290. df=6. Sig.=0,113. Does this mean my model is simply too bad to use for anything? Can I say anything credible about my data after running this analysis? Thanks for any clues that'll take me closer to the beauty of applied statistics... ;) Regards Alex


title: How to handle age at measure-start in Cox regression?

text: I'm a statistics newbie (medical student) trying my luck with a Cox regression for a survival analysis on the outcome of a specific type of operation. And I'm trying to determine which variables to control for. And how to handle the age of the patient at operation-time (measure-start)... Currently I've set the following as covariates in my analysis (PASW/SPSS):               Risk (Low/High)     Gender (Female/Male)     Shunt used in operation (False/True)     Artery stenosis (False/True)     Age at operation (Integer value)      I've set the 4 first variables as categories (dummy coded as 0 and 1), and specified the first value (0) as the reference value. But what about the patients age at the operation? I assume this must be controlled for? But the age at measure-start (operation) isn't time-dependent (is it?), but are always the same value, right? Does this mean I have to control for the patients real age aswell? And does this mean I have to use the `COX with Time-dependent Covariables` function instead of the regular `COX Regression`? Thanks for any help! Regards Alex


title: Akinator.com and Naive Bayes classifier

text: Context: I am a programmer with some (half-forgotten) experience in statistics from uni courses. Recently I stumbled upon http://akinator.com and spent some time trying to make it fail. And who wasn't? :) I've decided to find out how it could work. After googling and reading related blog posts and adding some of my (limited) knowledge into resulting mix I come up with the following model (I'm sure that I'll use the wrong notation, please don't kill me for that): There are Subjects (S) and Questions (Q). Goal of the predictor is to select the subject S which has the greatest aposterior probability of being the subject that user is thinking about, given questions and answers collected so far. Let game G be a set of questions asked and answers given: $\\{q_1, a_1\\}, \\{q_2, a_2\\} ... \\{q_n, a_n\\}$. Then predictor is looking for $P(S|G) = \frac{P(G|S) * P(S)}{P(G)}$. Prior for subjects ($P(S)$) could be just the number of times subject has been guessed divided by total number of games. Making an assumption that all answers are independent, we could compute the likelihood of subject S given the game G like so: $P(G|S) = \prod_{i=1..n} P(\\{q_i, a_i\\} | S)$ We could calculate the $P(\\{q_i, a_i\\} | S)$ if we keep track of which questions and answers were given when the used have though of given subject: $P({q, a} | S) = \frac{answer\ a\ was\ given\ to\ question\ q\ in\ the\ game\ when\ S\ was\ the\ subject}{number\ of\ times\ q\ was\ asked\ in\ the\ games\ involving\ S}$ Now, $P(S|G)$ defines a probability distribution over subjects and when we need to select the next question we have to select the one for which the expected change in the entropy of this distribution is maximal: $argmax_j (H[P(S|G)] - \sum_{a=yes,no,maybe...} H[P(S|G \vee \\{q_j, a\\})]$ I've tried to implement this and it works. But, obviously, as the number of subjects goes up, performance degrades due to the need to recalculate the $P(S|G)$ after each move and calculate updated distribution $P(S|G \vee \\{q_j, a\\})$ for question selection. I suspect that I simply have chosen the wrong model, being constrained by the limits of my knowledge. Or, maybe, there is an error in the math. Please enlighten me: what should I make myself familiar with, or how to change the predictor so that it could cope with millions of subjects and thousands of questions?


title: Classical measurement error model in R

text: I'm not sure if this is precisely a measurement error model or not. I'm working on performing meta analysis, and the model I'm starting with is fairly basic. \begin{aligned} X_i &= \mu_i + e_i \\\ Y_i &= \beta \mu_i + g_i + \delta_i \end{aligned} The random components are $e_i$, $g_i$, and $\delta_i$,and the variance is known for $e_i$ and $\delta_i$. This falls under a measurement error model with measurement in both the predictor and response. How would I fit this model in R?


title: Categorical Missing Value Replacement

text: I want to ask how I could replace categorical missing values with reference please?


title: R-squared result in linear regression and "unexplained variance"

text: I did a linear regression in R and got the following result:                                       Estimate Std. Error t value Pr(>|t|)         (Intercept)        192116.40    6437.27  29.844  < 2e-16 ***     cdd                   272.74      26.94  10.123 1.56e-09 ***     pmax(hdd - 450, 0)     61.73      22.54   2.738   0.0123 *       ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 16500 on 21 degrees of freedom     Multiple R-squared: 0.8454, Adjusted R-squared: 0.8307      F-statistic: 57.41 on 2 and 21 DF,  p-value: 3.072e-09       My question regards the R-squared value, 0.83 and what it means if I want to specify approximate percentage contributions of each (monthly) variable. EDIT: See data, below. Say I take the first 12 hdd and cdd data points, and calculate the sum of the 12 predictions (i.e. the first year's total prediction), using the coefficients, above. The baseline (intercept) contribution to the year would be approximately 12 * 192116.40 = 2305397, right? Similarly, the cdd contribution to the year would be approximately 1608 * 272.74 = 438565.9, and hdd would be (after my hand-made hinge function) approximately 1329 * 61.73 = 82039.17. Summing the three values yields 2826002, which is within 1.3% of actual total usage (2862840, the sum of the first 12 elec's). Can I then say that cdd contributes 438565.9/2826002= 0.1551895, or approximately 16% of the yearly total? Or do I need to take that and compensate for the adjusted R-squared: 0.1551895*0.8307= 0.1289159 (i.e. multiply by the adjusted R-squares), for approximately 13% of the total? Or is none of this correct reasoning? My data is:                    elec  hdd cdd     1  235940  880   3     2  205380  772   4     3  211780  551   9     4  192220  281  68     5  221440  165 119     6  304840   15 364     7  283160    4 434     8  300440   11 339     9  272900   42 214     10 204220  322  44     11 201060  592   8     12 229460  784   2     13 214520 1064   0     14 197900  719   2     15 186660  618  15     16 195340  332  88     17 241200  109 159     18 260700   18 282     19 299940   29 367     20 293240    2 426     21 268740   51 159     22 208380  319  36     23 183820  452   7     24 231360  903   0      (The monthly billing cycle for elec can be anywhere from 29 to 32 days, so that injects a lot of variance right there. I do not yet have all of the billing cycle lengths to to a trading day kind of adjustment.)


title: Expected values for chi-squared test on binned paired counts

text: I want to do a chi-squared test on data that looks like this:               A B     0 0      1 0     0 1     1 1     8 0     3 4     ...      You can think of each pair as one trial with two participants. In each trial, there are a different number of observations for each participant. I have binned each data pair like so: I counted how many pairs have 0 for both pairs (e.g. 0-0), how many have exactly one 0 (e.g., 0-1, 1-0, 8-0, etc.), and how many have greater than 0 for both (e.g., 1-1, 3-4, etc.). This gives me the following counts: Two zeroes: 227   One zero: 277   No zeroes: 146 The problem is that I am not sure how to calculate expected values here. Each pair represents the number of times something appears within an arbitrary number of observations; each pair represents a different number of observations, and each element of the pair does too. So, for instance, for a given 0-1 pair the first participant might have given 200 observations (with no hits), and the second might have yielded 150 observations (with 1 hit). Another 0-1 pair might have yielded 100 observations and 50 observations, respectively. So, in this case, just totaling the overall number of hits and dividing that by the total number of observations won't get me the right expected values... Am I missing something obvious here? For instance if we know that the relative frequency of hits per observation is 0.01217, and if there are a different number of observations for each trial, is there a simple way to get expected values in a problem like this? This is what I want to test: whether the hits between participants in each trial are independent. I expect that as the hits of one participant in one pair increases, so will the other. I realize that there might be better ways to test this, but my committee have asked me to use a chi-sq test, if possible. So what I expect to find is that the number of trials in the 0-0 and in the no-zero bins will be higher than would be expected if hits are distributed evenly.


title: How to graphically compare predicted and actual values from multivariate regression in R?

text: I'm trying to write a function to graphically display predicted vs. actual relationships in a linear regression. What I have so far works well for linear models, but I'd like to extend it in a few ways.   1. Handle glm models   2. Deal with NAs in the predicted values Does what I have so far seem like a good solution, or is there an existing package somewhere that's already implemented this?               DF <- as.data.frame(na.exclude(airquality))     DF$Month <- as.factor(DF$Month)     DF$Day <- as.factor(DF$Day)          my_model <- lm(Ozone~Solar.R+Wind+Temp+Month+Day,DF)          PvA<- function(model,varlist=NULL,smooth=.5) { #Plot predicted vs actual for a model              indvars <- attr(terms(model),"term.labels")              if (is.null(varlist)) {             varlist <- indvars         }              Y <- as.character(as.list(attr(terms(model),"variables"))[2])         P.Y <- paste('P',Y,sep='.')              DF <- as.data.frame(get(as.character(model$call$data)))         DF[,P.Y] <- predict.lm(model)              par(ask=TRUE)         for (X in varlist) {             print(X)             A <- na.omit(DF[,c(X,Y)])             P <- na.omit(DF[,c(X,P.Y)])             plot(A)             points(P,col=2)             lines(lowess(A,f=smooth),col=1)             lines(lowess(P,f=smooth),col=2)         }          }     PvA(my_model)


title: How to correct uneven sampling distribution when calculating the mean?

text: Suppose I have a function `f`, and I want to sample it at 100 points in the interval `[0, 100]`. For some reason (that seemed smart to me at the time), I decided to not sample at equidistant intervals, but rather use the following function to determine the sample points:               log2(x)*(100/log2(100))      This gives me a sequence of sample points that becomes denser as it approaches 100. The problem is, that now I need to calculate the mean over the values I have sampled, but due to the bad sampling that would be heavily biased. I cannot resample the data, this would take way too long (several days), and I am on a very tight schedule. So, the solution that comes to mind is to calculate a weighted average to correct the error. My question is, how do I determine the weights?


title: Should the average prediction for a given attribute value equal the rate for that value?

text: Let's say I'm predicting the likelihood that someone will buy a widget, using their age, eye color, and gender as input attributes. I split my data into a training set and a test set, and train up my predictive model using the (large) training set. Let's say that 15% of people in the test set with blue eyes buy the widget, should I expect the average predicted buy probability of people with blue eyes would be around 0.15 also (assuming a good predictive model)? Intuitively it seems like it should, but I don't want to assume that it's true.


title: Cycling in k-means algorithm

text: According to wiki the most widely used convergence criterion is "assigment hasn't changed". I was wondering whether cycling can occur if we use such convergence criterion? I'd be pleased if anyone pointed a reference to an article that gives an example of cycling or proves that this is impossible.


title: More efficient plot functions in R when millions of points are present?

text: I find R can take a long time to generate plots when millions of points are present - unsurprising given that points are plotted individually. Many of the points overlap and form a black mass and a lot of time is spent plotting more points into that mass. Are there any plotting functions that trade off processing time to work out in advance what plotting is redundant? I am not looking for a density plot (although those are often useful), I would want the same output as a simple plot call but much faster than millions of overplots if possible.


title: New mean calculation in multivariate gaussian

text: Given a Gaussian distribution $N(\mu_1,\sigma_1^2)$, i would like to choose another mean $\mu_2$ which is $2\sigma_1$ away from $\mu_1$. In this case our new mean $\mu_2=\mu_1\pm 2\sigma_1$. How do we calculate the new mean($\mu_2$) in multivariate case? I mean to say, when your multivariate Gaussian distribution is $N(\mu_1,\Sigma_1)$ and my $\Sigma_1$ is symmetric positive definite matrix. i.e $\left[ \begin{array}{cc} \sigma_x^2 & \sigma_{xy} \\\ \sigma_{yx} & \sigma_y^2 \end{array} \right]$.


title: Should I use Chi Squared/Fisher's Exact Test or McNemar's Test?

text: Say I did a survey of 100 people asking if they eat vegetables daily. Out of my sample population, 60 say they do, 40 say they don't. Now, I also have the breakdown of males and females. Say in the 60 who eat veggies daily, 45 are females and 15 are males. In the 40 who don't eat veggies daily, 15 are females and 25 are males. Now, my initial study didn't control for gender but I now want to see if gender has any effect on whether they eat veggies daily or not. Can I use Chi-Squared/Fischer's exact test with [45,15;15,25] as my contigency table? Or would I need to use a McNemar's Test and do paired matching? Also if i did chi-squared would this be considered a "Chi-Square Goodness of Fit Test" or a "Chi-Square Test of Independence"?


title: Calculating Stocks/Shares Odds

text: I am interested in the calculation of odds for stocks. For example, sports odds are like so:               Turkey vs Ukraine     hwin, draw, awin     2.20  3.40  3.20      2.20 for a Turkey Win, 3.40 for a Draw and 3.20 for a Ukraine Win. I am interested in the calculation of odds in such a format but for stocks in this case. Can you point me or show me a way to achieve this such that:               GOOGLE     1.20 4.50 3.51      That is: 1.20 for a rise in value, 4.50 for the value to remain the same, and 3.51 for the share value to increase. What will be needed in order to calculate such odds? Share price history etc?


title: Why does Random Forest Regression perform worse than autoregression

text: I have a dataset of NFL games. Each game has one row for each team in the game. Each team's row contains the team's statistics in that game (such as points scored, passing yards, red zone attempts, etc.), the team's statistics over the team's previous four games, and the team's opponent's statistics over their previous four games. So, for example, the one row might contain a team's passing yards in game `gt`, that team's total passing yards over games `gt - 4` through `gt - 1`, and the team's opponent's passing yards allowed over games `go - 4` through `go - 1`. What I have found is that, for some statistics, a team's previous four games' total (and possibly their opponent's previous four games' total) gives a better prediction of what they will do in a given game than a Random Forest Regressor trained on all the variables. Using the example from above:               >>> # keys beginning with 'tm_' represent the team's statistics in this game     >>> # keys beginning with '_tm_' represent the team's statistics over the previous 4 games     >>> # keys beginning with '_op_' represent the opponent's statistics (allowed) over the previous 4     >>>     >>> from scipy.stats.stats import pearsonr     >>> from sklearn.ensemble import RandomForestRegressor     >>>     >>> # Calculate r-squared between passing yards over the previous 4 and passing yards     >>> r, p = pearsonr(df['_tm_py'], df['tm_py'])     >>> r**2     0.10695998799573359     >>>     >>> # Use a random forest regressor to predict passing yards using all the features in the data set     >>> rgr = RandomForestRegressor(n_jobs=-1, max_features=None, n_estimators=500, oob_score=True)     >>> rgr.fit(df[features], df['tm_py'])     RandomForestRegressor(bootstrap=True, compute_importances=None,                criterion='mse', max_depth=None, max_features=None,                max_leaf_nodes=None, min_density=None, min_samples_leaf=1,                min_samples_split=2, n_estimators=500, n_jobs=-1,                oob_score=True, random_state=None, verbose=0)     >>> rgr.oob_score_     0.099712596456913757     >>>     >>>     >>> # print the top-10 most important features     >>> for f, i in sorted(zip(features, rgr.feature_importances_), key=lambda x: x[1], reverse=True)[:10]: print f, '  \t', round(i, 4)     ...     _tm_py      0.0929     _tm_rzpy    0.0098     _op_py      0.008     _tm_pts     0.0075     _tm_pc      0.0071     _tm_sfpy    0.007     _tm_pfd     0.0068     _tm_td      0.0067     _op_p2y     0.0065     _tm_p2y     0.0065      As you can see, a team's passing yards over their previous four games is the most important feature for predicting their passing yards (followed by the team's red zone passing yards over their previous four games and their opponent's passing yards allowed over their opponent's previous four games). However, the `oob_score_` for the Random Forest Regressor turns out to be lower than the r-squared between `tm_py` and `_tm_py`. Why is this? Is it because all of the other features are somehow "polluting" the direct relationship between these two variables? Is it possible that using a linear correlation between a single feature and the target can be better than using a random forest regressor between many features and the target?


title: Moving average ARIMA error term

text: Lets say that I have a time series data $Y_{t}$. I'm trying to forecast using am moving average MA(1) model using Box Jenkins methodology. The following is the equation for an MA(1) obtained from a textbook. $$ Y_{t} = \theta_{0} + e_{t} + \theta_{1} e_{t-1} $$ I assume we would use some type of nonlinear optimization to estimate coefficients $\theta_{0}$ and $\theta_{1}$ via maximum likelihood or conditional least squares. My question is about the very first error term $e$ when $t = 1$.   1. Does the very first error term $e_{1}$ is also automatically determined by the optimization ?   2. What would be the value for $e_{0}$ ? Thanks so much


title: Got a good Random Forest predictor. What's next?

text: Okay, I have created a Random Forest, I get a very good prediction of response. What's next? How can I improve it further? I ask this because I often see competitions in Kaggle, where RF is used as a benchmark, and many people manage to obtain scores much higher than that of the RF based solution. In some cases, the data is also unlabelled, preventing you to apply manual judgment. How do they beat RF by such decent margins?


title: Modeling user satisfaction with increasing sampling rate

text: Standard Deviation is inversely proportional to $\sqrt(Sampling\\_Rate)$ (MLE Consistency). Can user satisfaction with the reduction of error by increasing the sampling rate be modeled somehow other than the obvious way of performing a user study? Is there any literature that discusses similar ideas?


title: A regression decomposition for dummy variables

text: Suppose that I have a regression model \begin{align} Y = X\beta_X + Z_1 \beta_{Z_1} + Z_2 \beta_{Z_2} + \varepsilon \end{align} where $X$ are some controls, and $Z_1$ and $Z_2$ are two sets of dummy ($\\{0,1\\}$-valued) variables. I assume that $Z_2$ does not depend upon $Z_1$, but $Z_1$ might depend upon $Z_2$ in the sense that the probability that an element of $Z_1$ might be affected by the value of $Z_2$. I would like to decompose the information (variation) in $Z_2$ into a component that affects $Y$ only through $Z_1$, and a residual that affects $Y$ directly. If $Z_1$ were continuous, I could do something like regress (each coordinate of) $Z_1$ onto $Z_2$, then regress $Y$ on the fitted values. The natural analogue of this, say fitting a _logistic_ (or similar) regression model for $Z_1$ gives up linearity, losing interpretability, and takes much longer to compute. I would appreciate any insight on literature addressing this or similar problems.


title: Cross-section definition in a pooled regression

text: I'm trying to understand the results of a pooled regression model coming from this paper http://repository.cmu.edu/cgi/viewcontent.cgi?article=1010&context=sds It's a behavioral economics experiment made to measure the "offers" made in different groups of people across 20 time periods. Four "treatments" are being compared (called UG, PC2, RC2 and RC5, it doesn't matter how they differ, but the authors are basically interested in seeing how the evolution of the offers in the 20 time periods will differ between treatments). There is a VERY small sample size in each treatment (always less than 6 groups in each treatment). Here is the table summarising their results: ![enter image description here](http://i.stack.imgur.com/ug9bZ.png) There are a few things I don't understand about those results :   * the authors say they "pool the four treatments UG, RC2, RC5 and PC2". But the four treatments do not really correspond to different "cross-sections" or "units", as they are experimentally very different. What should be pooled instead are the time-series for each group within each treatment, am I right?   * authors say they pool the four treatments "with UG being the omitted category". What does that mean exactly? I've never seen this in a pooled regression before.   * why do you think they used a Wald test to test the significance of the results? If my questions are too general, I'd appreciate any online resource that you could point at me to understand more what was done in this analysis.


title: Automatic labeling of training set

text: I have once meet the following question, given a training set, is that possible to do the automatic labelling? In addition, if this training set consists of plain text files, is that possible to know what is the best description of positive class, and what is the best description of negative class? My feeling is that, automatic labelling can be treated as an unsupervised learning; finding the best description from a set of text documents can be treated as an information extraction problem. But I would like to know, is there any formalized methodology or paper to discuss this problem? What’s the state-of-art for this issue?


The question is simply what is stated in the title: **When does the law of large numbers fail?** What I mean is, in what cases will the frequency of an event not tend to the theoretical probability?