title: How can one construct a cumulative probability distribution function from 2 others?

text: I dip into project time estimation and can't find intuition. What is the cumulative probability distribution of an event when two independent tasks both complete successfully (when performed in parallel)? We have two cumulative probability distribution functions on 2 independent random variables that represent duration from project start to possible finish date (more time – more likely that a team completes a task): $$F_1(x) = P(\xi_1 \le x)$$ $$F_2(x) = P(\xi_2 \le x)$$ What is: $$F(x) = P(\xi_1 \le x\ \&\ \xi_2 \le x)$$ How can it be derived or represented from/by $F_1$ and $F_2$? **UPDATE** I agree that $$P(\xi_1 \le x\ \&\ \xi_2 \le y) = P(\xi_1 \le x) \cdot P(\xi_2 \le y) $$ for event $(\xi_1, \xi_2) \in R^2$ But misconception come from expected value of joint event: what is average time to complete both task? I think that it is: $$E(X) = \int_0^\infty xf_1(x)f_2(x)dx$$ Previously I think that it will be $E(\xi_1)\cdot E(\xi_2)$ - but this is wrong (you get **time^2** unit!!) and so doubt about multiplication low as think that something wrong in that part.


title: Proof of conditional probability

text: $X$ and $Y$ are real-valued random variables such that the distribution of $(X,Y)$ is absolutely continuous with density function $p$ and let $p_x$ denote the marginal density function of $X$. Suppose that there exists a point $x_0$ such that $p_x(x_0) > 0$, $p_x$ is continuous at $x_0$, and for almost all $y$, $p(\cdot,y)$ is continuous at $x_0$. Let $A$ denote a subset of $\mathbb{R}$. For each $\epsilon \gt 0$, let $$d(\epsilon)=\Pr(Y \in A | x_0 \leq X \leq x_0 + \epsilon).$$ Show that $\Pr[Y \in A|X = x_o]$ = $\lim_{\epsilon \to 0} d(\epsilon).$ **Attempt:** $$\eqalign{ \lim_{\epsilon \to 0} d(\epsilon) &= \lim_{\epsilon \to 0}\Pr(Y \in A|x_0≤ X≤ x_0+\epsilon)\\\ &= \Pr(Y∈A|x_0≤X≤x_0+0)\\\ &= \Pr(Y \in A|x_0≤ X ≤x_0)\\\ &= \Pr(Y \in A|X = x_0). }$$ I realize this attempt is wrong, but I'm not sure how to approach this problem.


title: Is it OK to use correlated variables for cluster analysis?

text: I know there is a series of regression diagnostics procedures (correlation, beta, residual, etc.) before, during, and after regression analysis. But, is there any common procedure to follow for cluster analysis (like, Ward)? What are the R commands? Thanks!


title: Efficient method/technique to update covariance matrix

text: A covariance matrix of multivariate random variable can be constructed given a time-series random variables. Eg. If you observe a student's performance in different objects (Math, English, Physics, etc) for a period of time; then you can construct the covariance matrix for those objects for that specific student. However, the random variable can be statistical unstable. Hence, the covariance matrix must be updated each time a new value of randome variable is observed. I'm looking for an efficient method/technique to update that covariance matrix. Although just the name for the method is enough, if you have a link to a tutorial, lecture note or a paper on the topic, it will be much more helpful too. Thanks, PS: Please correct me if I use some terminologies incorrectly. (I'm not a mathematician or statistician).


title: How to calculate a distribution out of a series of min-max values?

text: For each minute of an hour, I have a minimum and maximum observed value. for example:               time   min    max     00:00  12.13  15.10     00:01  14.23  17.02     00:02  11.12  12.10     ...     00:59  09.11  09.89      We can assume that, during that minute, the observed values were uniformly distributed between the minimum and the maximum values.   1. How can I calculate the hourly distribution of the observed values out of the 60 min and max values. My goal is to calculate the x-percentile values.   2. Do you know if this type of calculation has a name I can search for or if it is implemented in a statistical software, such as R?


title: Create n-variate normal distributed random numbers

text: I estimated the sample covariance "C" of a population and get a symmetric matrix. With "C" I like to create n-variate normal distributed r.n. but therefor I need the Cholesky-Decomposition of "C". What should I do if "C" is not positive semidefinite?


title: 95% CI or standard error of predicted risk in Cox regression using R

text: I am doing a multivariate Cox regression. I got the predicted risk for a new scenario based on the existing Cox model in R using function `predictSurvProb` in the `pec` package. Now I want to figure out the 95% CI of this predicted risk. How can I do it? Here are my codes:               # fit a Cox model     f.coxph <- coxph(Surv(recurtime, recurrence) ~ age + pgs.f + svi.f + ece.f                 +  margin.f, data=prostate)      # predicted 5-year recurrence risk for a 55-year old, pathological Gleason score 6,      # positive surgical margin, negative SVI, and negative ECE     case1 <- data.frame(age=55, pgs.f="GS 5-6", margin.f="Pos",               svi.f="Neg", ece.f="Neg")     predictSurvProb(f.coxph, newdata=case1, times=60)      Or how to calculate the standard error of this predicted recurrence risk? So I can use `fit+/-1.96*se.fit` to get the 95% CI. Thanks!


title: Is it possible to construct a discrete-time multilevel hazard model in R?

text: I'm trying to run a discrete-time multilevel hazard analysis comparable to the model proposed by Barber et al. I am attempting to model the hazard of migrating internationally using predictors at the individual, household, community, and regional levels. (Most of the variables of interest are at the individual, community, and regional levels--I just need to account for clustering by household) Is there a package that will allow me to do this in R? I've used lme4 for regular multilevel modeling, but can it also be used for multilevel survival models? How would I go about coding such a model? (And if this can't be done in R, can it be done in Stata, or do I have to bite the bullet and buy and learn HLM or MLwiN?) Please let me know if you need any additional information. Thanks! ETA: Barber et al. refers to: Barber, Jennifer S., Susan A. Murphy, William Axinn, and Jerry Maples. 2000. "Discrete-Time Multilevel Hazard Analysis." Sociological Methodology, 30: 201-235.


title: Generate random numbers following a distribution within an interval in R

text: I need to generate random numbers following Normal distribution within the interval $(a,b)$. I know the function `rnorm(n,mean,sd)` will generate random numbers following normal distribution,but how to set the interval limits within that? Is there any particular R functions available for that?


title: How to predict the amount of data needed for modeling?

text: Is there a way to estimate the amount of data (or the number of records) required to build a statistical model? I read few blogs and I feel that most of the responses concur that there is no way or it is extremely hard to predict sample size for an application. One one blog asks for 10 times the total number of features, less if I am using regularized version of the ML model. A sample size of 50 seem to be the minimum. One method suggests to build the model and check the generalization error. If the error is unacceptable, then gather more data and iterate. Considering the data collection for my application to be time consuming and extremely costly, what options do I have in time-constrained business organization?


title: Horseshoe distribution

text: Studying nonlinear dimension reduction topics I came across a "horseshoe distribution" or "horsheshoe phenomena". Figure (b) [source, p. 12] depicting horseshoe distribution is included below. ![Horseshoe distribution.](http://i.stack.imgur.com/DfT0V.png) I do not have a good idea how to simulate this type of distribution. Thanks in advance for any pointers.


title: Is there any free software for Oden's Ipop?

text: I'm trying to use software that has an Oden's Ipop Statistics function for testing spatial autocorrelation adjusted for population. I prefer knowing if something is available for FREE. I have searched through the internet, but unfortunately I did not find nothing. I have one option for this, which is for me to use `R` and make a function of Oden's Ipop Statistics using the formula for it. But just in case someone knows of any free software for the Oden's Ipop, that would be a helpful to me with my research study.


title: PCA on independent but non-identically distributed multinormal data

text: In the empirical phase diagram approach the data are given by a set ${\boldsymbol X}$ of input (controlled) variables (such as pH and temperature) and a set ${\boldsymbol Y}$ of output variables. The empirical phase diagram aims to visualize (at least, roughly) the relation between the input and the output. The strategy runs as follows in case when there are only two input variables:   * firstly we ignore the input ${\boldsymbol X}$ and we run a PCA on ${\boldsymbol Y}$   * keep the first three principal components    * transform the three coordinates on the principal components into a color using the red-green-blue coding   * plot the color in function of the two input variables, and visualize I am rather new in PCA analysis (so please do not hesitate to tell me my questions are stupid if they are). I know PCA has a better interpretation when it is applied with a sample of i.i.d. random multivariate normal variables. But I don't know what are the possible pitfalls when this assumption does not hold. Assume for instance a multivariate regression model for which the distribution of a single multivariate response ${\boldsymbol Y}$ is assumed to be: $${\boldsymbol Y}_i = f({\boldsymbol X}_i) + \epsilon_i \quad \textrm{with } \epsilon_i \sim {\cal N}({\boldsymbol 0}, \Sigma).$$ I think that we ideally should run the PCA on the centered responses ${\boldsymbol Y}_i - \hat{f}({\boldsymbol X}_i)$ in such a situation. So what are the possible pitfalls if we run the above strategy in such a situation ?


title: Can unsupervised evaluation measures for clustering replace a supervised evaluation measure?

text: Is it possible to have the same evaluation performances when comparing some clustering algorithms using many unsupervised evaluation measures instead of a supervised one ?


title: Software for advanced visualization

text: This question might not be directly connected with the statistics but nevertheless as many people here deal with data visualization, I would like to ask if anyone of you uses some special tools for visualizing data or could recommend any resources for creating special effects with data. Here is an example video with some possible cases Thank you!


title: How to solve/compute for normal distribution and log-normal CDF inverse?

text: > **Possible Duplicate:**   >  Calculating percentile of normal distribution Edit: I think the hint is Box-Muller transformation. Any ideas? I am trying to find the inverse of normal CDF. This is normal CDF I got exactly from Wiki, $\Phi (x)=\frac{1}{\sqrt{2\pi}}\int\limits_{-\infty }^{x}{{{e}^{-{{t}^{2}}/2}}}dt=\frac{1}{2}\left[ 1+erf\left( \frac{x}{\sqrt{2}} \right) \right]$. This is the log-normal CDF, ${{F}_{x}}(x;\mu ,\sigma )=\frac{1}{2}erfc\left[ -\frac{\ln x-\mu }{\sigma \sqrt{2}} \right]=\Phi \left( \frac{\ln x-\mu }{\sigma } \right)$. I was told there is an easy way using radian, but don't know how. Any experts would know? I have seen this: http://home.online.no/~pjacklam/notes/invnorm/ This is part of the inverse transformation I am doing. R or Java will be used, an explicit formula would be helpful enough


title: What is the difference between the Mann-Whitney and Wilcoxon rank-sumtest?

text: First of all, what is the difference between the Mann-Whitney and Wilcoxon rank-sum test? How do I choose between the two? I know the latter is implemented in R, is there a way to do the Mann-Whitney test? Is there a way I can include my (very large) datasets in this post so you can get an idea of what I'm trying to do?


title: How to convert the wilcox.test W value output for two independent groups to U in R?

text: I need the U values of Mann-Whitney test, but the wilcox.test in R outputs only W values. I do not think they are identical, so here is my question - how to convert W to U?


title: conditional expectation of squared standard normal

text: Let $A,B$ independent standard normals. What is $E(A^2|A+B)$? Is the following ok? $A,B$ iid and hence $(A^2,A+B),(B^2,A+B)$ iid. Therefore we have $\int_M A^2 dP = \int_M B^2 dP$ for every $A+B$-measurable set $M$ and hence $E(A^2|A+B) = E(B^2|A+B)$. We obtain $2 \cdot E(A^2|A+B) = E(A^2|A+B) + E(B^2|A+B) = E(A^2+B^2|A+B) = A^2+B^2$ where the last equation holds since $A^2+B^2$ is $A+B$-measurable. Finally we have $E(A^2|A+B) = \frac{A^2+B^2}{2}$.


title: Doubt on Covariance Matrix in Weighted Least Square Estimation

text: ![this is from a book of gilbert strang..](http://i.stack.imgur.com/jlF6b.png) This is a page from the book linear algebra, geodesy and gps by Gilbert Strang. The page explains about the justification of the inverse of the of the variance-covariance matrix of measurement vector $b$ in the overdetermined system $Ax = b$ as the best weight matrix for best estimate of x.   1. what does the line referring to "errors... contained in the matrix" (underlined in blue) mean?    2. how does that substitution take place marked with a blue curve...   3. is $E[rr^T] = E[bb^T]$ ?? any explanation is welcome...


title: Question about standard deviation and central limit theorem

text: I have a quick question about the central limit theorem. Lets say I measure some value that comes from an arbitrary distribution N times and I repeat this M times. I understand that if I calculcate the mean from the N values I will have a set of M values that follows a normal distribution. But what if I measure the sample standard deviation from N, will my resulting distribution also be normal? Following the derivation of the CLT I do not see this to be the case, but intuitively I think this is true, at least for some distributions. Any light on the issue would be greatly appreciated. First I will quote the CLT from wiki: > the central limit theorem (CLT) states that, given certain conditions, the > arithmetic mean of a sufficiently large number of iterates of independent > random variables, each with a well-defined expected value and well-defined > variance, be approximately normally distributed. My question then is a variant on the quote from the wiki page: will the the central limit theorem (CLT) state that, given certain conditions, the STANDARD DEVIATION of a sufficiently large number of iterates of independent random variables, each with a well-defined expected value and well-defined variance, be approximately normally distributed? Here is a series of plots I made from random numbers that follow a beta distribution. I have generated 1000 sets of 5000 points each. The first plot is a histogram of the first set. The second is a histogram of the 1000 calculated means, and the 3rd is a histogram of the 1000 calculated std. ![Simulation results](http://i.stack.imgur.com/Ozb5q.png) ![Results 2](http://i.stack.imgur.com/Wk8rj.png) ![Results 3](http://i.stack.imgur.com/1ab1U.png)


title: Estimating mean and st dev of a truncated gaussian curve without spike

text: Suppose I have a black box that generates data following a normal distribution with mean m and standard deviation s. Suppose, however, that whenever it outputs a value < 0 it does not record anything (can't even tell that it's outputted such a value). We have a truncated gaussian distribution without a spike. How can I estimate these parameters?


title: Calculating linguistic diversity of a text with respect to size

text: I'm doing linguistics work and I'm trying to make a measure of linguistic diversity in a text. A simple calculation seems to be linguistic uniformity ( or diversity) = number of words / unique words But this must "privilege" shorter texts like poetry for diversity, since the author is less likely to repeat themselves in so short a space. For example, Blake's poetry gets a score of 5, whereas the KJV bible scores 79 and this has a lot to do with the size of the bible compared to Blake's poetry. If we analysed more of the Blakes's corpus, his score would definitely go higher. How can I normalize the calculation, taking into account the inevitable repetition that comes with longer texts, and the fact that shorter texts are so small that we just haven't seen an appreciable part of an author's lexicon? I know that there's a name for this, I just can't remember what it is. I've looked at the wiki article on diversity index but it has not helped me.


title: The distribution of the linear combination of Gamma random variables

text: If $X_i\sim\Gamma(\alpha_i,\beta_i)$ for $1\leq i\leq n$, let $Y = \sum_{i=1}^n c_iX_i$ where $c_i$ are positive real numbers. Assume all the parameters $\alpha_i$'s and $\beta_i$'s are all known, what is $Y$'s distribution ?


title: Recommend some books/articles/guides to enter predictive analytics?

text: What learning material would you suggest for a CS person / novice statistician / novice mathematician to get into predictive analytics?


title: Books on statistical ecology?

text: I know this question was asked before: Reference book for ecological studies but it is not what I am looking for. What I am looking for is if anyone could recommend a good book (or a canonical reference) on statistical ecology? I have a very good understanding of statistics so the book could really be at any level. I would be using the book to teach myself more about the application of statistics in ecology than anything else so even an introductory book with good/interesting examples would be much appreciated. Also, my research tends to be geared towards Bayesian statistics so a book incorporating Bayesian statistics is even better!


title: Confidence interval based on time series

text: I have a timeseries of data which was gathered by driving a car around randomly. One data point was gathered every minute and each data point is either "Yes" or "No". Yes = temperature is above a threshold, No = it is below the threshold. Assume the car never stops. I'm interested in the proportion of Yeses and related confidence intervals. I'm inclined to use binomial confidence intervals except that it seems to me that the data points are not independent.   Would the fact that the data was gathered by driving around randomly buy me independence? How could I derive meaningful confidence intervals?


title: How would I calculate the expected change?

text: Based on the info below, how would I calculate the expected change in the test scores for males and females? The company decides to test for differences in performance across genders. To do so, they run a regression and obtain the following results:               Predicted Change in Score = 3 + .25*Female      R2 = .7      Female is a variable whose value is $1$ for females and $0$ for males. The t-statistic on the constant is $10$ and the t-statistic on the coefficient for female is $.9$. The company elects to use a $5\%$ level of significance. The class is very large, in excess of $30$ students.


title: One vs All appraoch to multiclassification in LIBSVM

text: I'm working on multiclass classification problem (precisely 4 classes). I want to use the most simple approach for this problem: one-vs-all! I have 4 different test set (only labels are different). Using one-vs-all approach, during test, for each input pattern, I have to compute 4 different objective function values from 4 different SMVs. So, the pattern will belong to the class with the greatest objective function value. So, I tried this: ./svm-train -s 0 -t 5 -c 16 -g 0.05 -b 1 'traindata' ./svm-predict -b 1 'testdata' 'traindata.model' 'outfile' I'm not sure about the means of **b parameter**. Opening the file 'outfile' generated from svm-predict, I get probabilistic values. I suppose that comparing these probabilistic values between all SVMs, I'm able to classify each pattern in the class with the greatest probabilistic value. My assumption is correct? In other words, I want to be sure that these probabilistic values are directly proportional to the relative objective function values.


title: How to predict Ozone concentration in few years time?

text: my friend is a chemist and his problem is to predict the level of ozone concentration in a single site. We have the data for the last 12 years. We want to predict the concentration for the coming years (as much as possible). I know this data set is small, so my question is is this possible? And how much of data is required if not? What is the tool/method to use? Here is a plot of my data: ![alt text](http://i.stack.imgur.com/KP3IQ.png) any suggestions/ help?


title: Kappa for Predictive Model

text: The "standard" way to compute Kappa for a predictive classification model (Witten and Frank page 163) is to construct the random confusion matrix in such a way that the number of predictions for each class is the same as the model predicted. For a visual, see (right side is the random): ![alt text](http://i.stack.imgur.com/oetej.jpg) Does anyone know why this is the case, instead of truly creating a random confusion matrix where the prior probabilities drive the number of predictions for each class. That seems the more accurate comparison against "a null model". For example, in this case, the number of actual and predicted classes would coincide (in the image uploaded, this would mean that the columns of the random confusion matrix would be 100, 60 and 40 respectively). Thanks!   BMiner


title: Determining sample size for a comparative study involving two unequal groups

text: I will like to know the formula for calculating the minimum sample size required for a comparative study involving 2 groups with unequal size. Group 1 is 10 times greater than group 2.


title: Hyperprior density for hierarchical Gamma-Poisson model

text: In a hierarchical model of data $y$ where $$y \sim \textrm{Poisson}(\lambda)$$ $$\lambda \sim \textrm{Gamma}(\alpha, \beta)$$ it appears to be typical in practice to chose values ($\alpha, \beta)$ such that the mean and variance of the gamma distribution roughly match the mean and variance of the data $y$ (e.g., Clayton and Kaldor, 1987 "Empirical Bayes Estimates of Age-Standardized Relative Risks for Disease Mapping," _Biometrics_ ). Clearly this is just an _ad hoc_ solution, though, since it would overstate the researcher's confidence in the parameters $(\alpha, \beta)$ and small fluctuations in the realized data could have large consequences for the gamma density, _even if_ the underlying data generation process remains the same. Furthermore, in _Bayesian Data Analysis_ (2nd Ed), Gelman writes that this method is "sloppy;" in the book and this paper (starting p. 3232), he instead suggests that some hyperprior density $p(\alpha, \beta)$ should be chosen, in a fashion similar to the rat tumors example (starting p. 130). Although it's clear that any $p(\alpha, \beta)$ is admissible so long as it produces a finite posterior density, I have not found any examples of hyperprior densities that researchers have used for this problem in the past. I would greatly appreciate it if someone could point me to books or articles which have employed a hyperprior density to estimate a Poisson-Gamma model. Ideally, I am interested in $p(\alpha, \beta)$ that is relatively flat and would be dominated by the data as in the rat tumor example, or a discussion comparing several alternative specifications and the trade-offs associated with each.


title: Test similar to Kruskal-Wallis test but for multidimensional vector

text: I need to do the similar kind of test as Kruskal-Wallis test, but for a multidimensional vector. Would you please kindly suggest such kind of test or provide some references for N-dimensional case


title: log-probit model: Calculation of confidence intervals for ED50 data

text: I calculate log-probit models according to the following method (many thanks to COOLSerdash, Aniko, whuber):               # data: textbook from Erna Weber     dosis <- c(2.90, 3.36, 3.90, 4.52, 5.52, 6.09,6.90)     nges <- c(10,10,10,10,10,10,10)     nok <- c(1, 2, 3, 4, 7, 6, 8)     edx.data <- data.frame(dosis, nges, nok)          # glm     glm.logit <- glm(cbind(nok,nges-nok) ~ log(dosis), family=binomial('probit'), data=edx.data)          # calculation of EDx data and confidence intervals     r <- dose.p(glm.logit,p=seq(0.1,0.9,0.2))     se <- attr(r,'SE')     xr <- data.frame(as.matrix(exp(r)))     names(xr) <- c('Dose')     xr$Lower <- exp(r-se)      xr$Upper <- exp(r+se)      The output is as follows:               > print(xr)                  Dose    Lower     Upper     p = 0.1: 2.827593 2.444115  3.271238     p = 0.3: 3.920316 3.606605  4.261315     p = 0.5: 4.915905 4.576043  5.281009     p = 0.7: 6.164330 5.580168  6.809647     p = 0.9: 8.546535 7.219686 10.117236      EDx data ("Dose") are according to the textbook by E. Weber. However, confidence intervals should be as follows:                 p   ED(p)   lower    upper     0.1  2.8272  1.6574   3.4642     0.3  3.9203  3.0380   4.5064     0.5  4.9159  4.2478   5.8820     0.7  6.1643  5.2973   8.6082     0.9  8.5476  6.7830  16.0300      The question is: where is my mistake? How to calculate these confidence intervals? Thank you very much in advance.


title: finding an optimal subgroup of binary indicators

text: My dependent variable is continuous variable that measures the (potential) success of a person in some activity. I have hundreds of binary indicators, each indicates about the existence of a specific beneficial behavior or property (e.g. the person typically wakes up before 8am). each of these indicators is positively linked to the outcome both theoretically and empirically (although some of the empirical differences are insignificant). Of course that the indicators are correlated. My task is to find a subset of 10-20 indicators that best predict the outcome. Ideally, what I want to have is a subset of indicators from which I'll construct a single variable - the number of positive indicators from that subset for a person, and wish that this variable has the highest possible correlation with the outcome. What I did first is a linear regression with all the indicators in. It gave me lots of negative coefficients (which makes sense mathematically, but is not intuitive as you'd expect that each variable will contribute something), and it didn't help much in choosing the subset. I have also tried a sort of a genetic algorithm optimization to help me choose such subset. It helped, but I feel I'm missing some statistical way of doing that (PCA perhaps?), and would be happy to get any leads/references/suggestions about approaching this problem, which looks like a model selection problem. Any ideas? Thanks.


title: Model Formulae. correct model?

text: I'm trying to build the following model in R, however I'm quite confused about the model formulae to use to include an interaction (x1 and x2) $Y_{}=a+b*x_1+cx_2+d(x_1*x_2)$ this intuitive formula seems to be wrong model:               y~x1+x2+(x1*x2)       is this the appropriate formula?:               y~x1*x2 ?


title: ACF and PACF plot analysis

text: I am new to ARIMA, and I am trying to understand these lag plots. Are the following ACF and PACF suggesting that the lag of my time series is 4? If I am wrong, please help me understand these plots. > ![enter image description here](http://i.stack.imgur.com/N5x9U.png)


title: Lower classification rate than expected by chance

text: I'm using scikit-learn for a small sample (36) classification problem with three features and three outputs (one output is binary and the other two are ternary). I'm using separate classifiers for each output, and I'm achieving reasonable classification rates for two outputs, but the remaining one, which is a 3-class problem, is abysmal. I've tried many different classifiers (SVM, random forest, logistic regression, naive bayes), and using leave-one-out cross-validation the classification rate fluctuates between 0% and 33% (this only happens when the classifier always outputs the same class). What I don't understand is how it's possible for the classifier to achieve a lower rate then expected by chance. A thought experiment says that when it achieves 0% classification rate, I could ignore the predicted class and randomly choose between the other two, achieving 50% classification rate, which would be better than chance. Why is this not happening automatically?


title: How to assess the importance of the features which come from intersection of features of the two models?

text: I have two models from two different data sets. Model 1 contain 50 features and model 2 contain 40 features. the intersection of features of model 1 and 2 is 10. so how can I assess the relative importance of the features which comes from the intersection of two models ?


title: How to generate random categorical data?

text: Let's say that I have a categorical variable which can take the values A, B, C and D. How can I generate 10000 random data points and control for the frequency of each? For example: A = 10% B = 20% C = 65% D = 5% Any ideas how I can do this?


title: How to model categorical (discrete-valued) time series?

text: Just want to make a little survey, What are, according to you, the best approach to model categorical time series? I'm building a model able to generate time series reproduicing the characteristics of a set of data. The available data is a discrete time, discrete-valued time serie. I have generated synthetic time serie using Markov chain, semi-markov chain and moving block bootstrap, I would like to know if you think of better way to do that job.


title: Obtain the asymptotic distribution of X/Y if X and Y are iid and independent of each other

text: ![example](http://i.imgur.com/gfshR1C.png). We just covered large sample theory and I'm going through the examples in our textbook. This one kinda confused me and I was hoping someone could help me understand why and how the book got, $W_i=\mu_X*Y_i - \mu_Y*X_i$ and how the book derived 6.4.3. I tried asking around but the response I keep getting is that each steps are explained but that doesn't help because they seem vague to me.


title: Definition of "degree of interaction" in the MARS model

text: What exactly is the meaning of the "degree of interaction" or "interaction degree" in the MARS model? e.g. R: `earth(..., degree=2)`


title: Structural changes in the observed standard deviation

text: Supposing to have a stochastic process with known normal distribution N(m,s) we can observe samples from. For quality test purpose, we define a threshold for good/bad samples equal to 3 standard deviation from mean value. Hence we accept a bad sample ratio equal to 0.28 %. We collect samples originated from such underlying stochastic process and we like to test if they match the expected distribution in terms of bad sample ratio or there is a structural change in the underlying process (i.e. the standard deviation of their underlying probabilistic distribution is different from the expected one). Please, what is the correct approach ? I know about "cusum", however it is applied to structural changes to the mean value, not the standard deviation.


title: Getting heads in an infinite series of coin flips- certain or almost certain?

text: A while ago I was introduced to the idea of something with probability happening "almost certainly", (or "almost surely"). Conceptually as I understand it, the probability of an event may be 1, but that does not mean it is certain, only _almost_ certain. Wikipedia gives this example, based on throwing a dart at a square: > Consider the event that "the dart hits the diagonal of the unit square > exactly". The probability that the dart lands on any subregion of the square > is proportional to the area of that subregion. But, since the area of the > diagonal of the square is zero, the probability that the dart lands exactly > on the diagonal is zero. So, the dart will almost never land on the diagonal > (i.e. it will almost surely not land on the diagonal). Nonetheless the set > of points on the diagonal is not empty and a point on the diagonal is no > less possible than any other point, therefore theoretically it is possible > that the dart actually hits the diagonal. So my basic question is: if I flip a fair coin an infinite number of times, am I **certain** , or merely **almost certain** , to get heads at least once? **Thought experiment** My initial reaction to the question was that the answer must be almost certain. However, I came up with a thought experiment which made me doubtful: Two immortal scientists are tasked with monitoring a fair coin which will flip once a minute. They both share the same coin, but are given different instructions. The first scientist's instructions are: If the coin has not come up heads after infinite flips, report "false". If the coin comes up heads on any single flip, report "true". The second scientist is given this flow chart to describe what he should do: ![coin flipping](http://i.stack.imgur.com/kcj8M.png) Given this, I would say the following:   * There is clearly no flow through the second scientist's diagram where there is any possible result apart from "true".    * There is no sequence of coin flips where, at any point either:      * The two scientists report different results OR     * One scientist has reported a result while the other is still yet to report a result. The first of those points seems to mean that it is certain, rather than almost certain, that the second scientist will report true. The second of those points seems to mean that the second scientist being certain to report true implies that the first scientist is certain to report true. Putting those together, the first scientist is certain to report true, giving a counterintuitive answer to the basic question. So to wrap up, is it certain or almost certain that an infinite series of flips of a fair coin will result in heads at least once? If **certain** , then what's the explanation for that counterintuitive result given that there's nothing forcing any single coin to come up heads? If **almost certain** , then what is the flaw in the reasoning of the thought experiment?


title: How can I estimate a weighted linear regression with fixed effects and many groups in R?

text: I have a data set derived from administrative registers covering the population of a small European country, containing a large number of defined groups (15 000+). For each of these groups I have multiple observations (different years), and two explanatory variables that characterize the exposure of each group to a policy intervention. In STATA I would estimate this using the regression command, with aweight=group_size and using absorb(group_identifier). I've done this before, and it is extremely fast. How do I do this in R? In R, my problem is that I cannot find a package that allows me to both weight by group size AND include fixed effects. Using the standard lm function I could include the group-size weight, but would have to estimate the group-effects explicitly (simply adding the groups identifier as a factor variable). This is extremely time consuming (I cancelled it after it had run for several minutes). Alternatively, I could turn to fixed-effects packages - such as the lfe package. These allow for fixed effects, but I have not found any such package that also accepts weighting by group size. Any suggestions?


title: Taking log of variables for my GMM study.

text: I am assisting a professor in implementation of Arellano Bond estimator for a study on Factors, which are crucial determinants for Bank profitability. Someone suggested me to take log of the variables to remove the Heteroskedasticity when implementing Arellano Bond estimator. Here are my Doubts:   1. Can't i just use the method available in Stata to give me the robust estimators, thus taking care of Heteroskedasticity. I am talking about xtabond method.   2. Even if i have to take log of variables, some of them are negative, or zero. Is there any way i can circumvent this problem? Kindly Help


title: Combining several a/b tests

text: Imagine I have two ads (banners), A and B. They can be placed at several different positions on the page (only one banner is shown on the page). Each position can be viewed as a separate a/b test. I do not control showing of banners, I only have access to clicks / shows data. I need to determine which ad is more clickable. If there is only one possible position, I only need to compare binomial proportion confidence intervals, which is easy enough. But consider the following example: Position 1:   A: 3 clicks out of 100 shows   B: 1 click out of 5 shows Position 2:   A: 0 clicks out of 20 shows   B: 25 clicks out of 200 shows As you can see, both tests favor ad B, but significance is low. How do I combine two tests to improve statistical significance?


title: Polynomial model with unpaired data

text: I'm trying to model data as a 2nd degree polynomial, but the data is unpaired and each data point of average values has a standard error for each axis. My data:   1. A time series in minutes (time tissue is exposed to cryoprotectant)   2. A % value of tissue survival after liquid nitrogen storage (3 replicates for each time point)   3. A % value for the uptake of cryoprotectant in tissue directly after exposure to cryoprotectant (3 replicates for each time point) The tissue used for the uptake is of the same type but is from a different source, as the uptake measuring protocol is destructive. So there are 2 sets of 3 values for each time point. However, a scatterplot of this data does not make sense because the survival and the regrowth values are not paired - I could mix and match my values within each time point. At low levels of cryoprotectant, survival is low because of ice formation. At high levels of cryoprotectant survival is low because of cryoprotectant toxicity. Therefore:   * I would like to model survival vs. uptake to find the maximum survival value's corresponding cryoprotectant uptake. I can:   * Perform a polynomial fit on survival vs. uptake using the average values for each time point   * Use that fit to find a predicted optimum uptake % I cannot:   * Get a p-value for the fit that is valid (i.e. one that assesses the error at each time point for both variables) Therefore is it possible to model such data with a valid p-value, or does the fact that it is unpaired mean it doesn't make sense? I'm currently carrying this out in R.


title: Skew of log-normal distribution using sciPy

text: I am doing bioinformatics and I am trying to fit some values to a log-normal distribution with `python's sciPy` version 0.11. According to the skew of the resulting distribution, I would like to make a 0-1 decision, i.e. if skew is positive, give a "0" value, as most of the data is on the left, while, if skew is negative, give a "1" value, as most of the data is on the right. While trying to implement that, I am confused as to which formula to use for the skew. From Wikipedia, the skew is equal to: $(e^{\sigma^2}+2)\sqrt{e^{\sigma^2}-1}$, but this can never be a negative number. Using `sciPy's` `moment` method like this: `scipy.stats.moment(data,3)/(std**3)`, seems to yield negative numbers, but then I found out this is the non-central moment, so I am afraid it is wrong. Secondly, I have looked into this very useful question here, on how to actually fit data to the log-normal distribution. I have tried to make a toy example, where I generate 50 numbers on the interval $0.1-9$ and 1000 numbers on the interval $9-12$. So, I would expect that this distribution should have a negative skew. The code is as follows:               import numpy as np     import scipy.stats as s     from matplotlib import pyplot as plt     test1 = np.linspace(0.1,9,50)     test2 = np.linspace(9,12,1000)     test = list(test1) + list(test2)     shape,loc,scale = s.lognorm.fit(test,floc=0)     distro = s.lognorm(shape,loc,scale)     x = np.linspace(0.1,20,200)     plt.plot(x,distro.pdf(x))     plt.plot(x,distro.cdf(x))     plt.axvline(x=scale,color='r',ls='--')     plt.legend(('PDF','CDF',r'$\mu$'), loc='upper left')     plt.hist(test,bins=100,normed=True,alpha=0.4)      And the corresponding image: ![A toy example histogram and fitted log-normal distribution](http://i.stack.imgur.com/xzFIo.png) My question is, does it really take that much to shift a log-normal distribution to the right? What I have done here seems to spread it wide, but shift it only a little. Should I use a different measure for my purposes?


title: Random sampling the data again to keep equal size across groups

text: I have an experiment with a between-group variable. It has 3 groups in total. When conducting the experiment, although we did random assignment for participants, the size of the 3 groups end up differently: 15, 21 and 27, respectively. With the data I want to conduct one-way ANOVA. I also have a repeated-measure variable, with which I also want to conduct mixed ANOVA. Currently my data violates the homogeneity of variance assumption. It seems with balanced group size ANOVA is pretty robust against the violation but not if sizes are unbalanced (reference). My questions is, could I randomly sample the 2 bigger groups again (without replacement) into size of 15, so that to keep the size equal across groups? If so, is there any theoretical support for that? Or I should just use robust non-parametric methods instead?


title: How can i perform Cosinor analysis in Matlab?

text: I have monthly and daily data from births in Greece and i want to to do a cosinor analysis to find the rhythm, (circadian,circannual,half yearly etc etc ) How can i use Matlab to perform the cosinor analysis? I found this link http://www.mathworks.com/matlabcentral/fileexchange/20329-cosinor- analysis/content/html/cosinor.html Not sure if it can help. Is there anything else that i could do with these data apart from cosinor analysis? I did a fourier spectral analysis and i found a strong 12 month cycle and a strong 2 month cycle.


title: Can't replicate the results of Golub's Generlized Cross Validation procedure (GCV)

text: I'm trying to understand this paper: Generalized Cross-Validation as a Method for Choosing a Good Ridge Parameter Author(s): Gene H. Golub, Michael Heath and Grace Wahba Source: Technometrics, Vol. 21, No. 2 (May, 1979), pp. 215-223 It's quite a well known paper and I can't think it would be incorrect (but obviously my understanding of it is!) The trouble is that when I try to implement their main formula (formula 1.4) on the raw data I don't get the same answer as when I use what is supposed to be a different form of the same formula shown here: Use of the Bootstrap and Cross-Validation in Ridge Regression Author(s): Nancy Jo Delaney and Sangit Chatterjee Source: Journal of Business & Economic Statistics, Vol. 4, No. 2 (Apr., 1986), pp. 255-262 Their formula shown in section 2.2.5 (is more intuitively appealing and I understand it) but it doesn't give the same answer as the one in the Golub et al paper. In the Golub paper GCV is defined to be PRESS (predicted residual sum of squares) carried out on transformed data (as shown in the paper). The transformed version is 2.3 on p217. They also state that the GCV is the value that minimizes either 1.4 or 2.3 (i.e. the models based on the original and transformed data respectively). So if I understand this correctly even though the minimum value will not be the same between the original and transformed models, the parameter value which minimizes them be the same (is this correct?) However when I implemented both methods they were not minimized by the same value. The first was minimized by approx 2 and the 2nd by approx 1.5. So my questions are: 1.) Are equation 1.4 in Golub and the one shown in section 2.2.5 of Delaney and Chatterjee equivalent? If so can you show me how one can be derived from the other. 2.) Is my coding of each model correct? I can post these papers up if someone tells me how.               function [GCV]=GCV(y,x,k)          [n, p]=size(x);          % the Delaney and Chatterjee method          h=(x'*x+k*eye(p));     h=h\x';     h=x*h;          den=(1-diag(h));     den=sum(den);     den=den/n;     den=den^2;               sum_e_i_2=0;          for i=1:n          y_e=y;     y_e(i,:)=[];          x_e=x;     x_e(i,:)=[];               pseudo = sqrt(k) * eye(p);     x_ridge  = [x_e;pseudo];     y_ridge  = [y_e;zeros(p,1)];               b2 = x_ridge\y_ridge;               e_i=y(i)-x(i,:)*b2;          sum_e_i_2=sum_e_i_2+e_i^2;          end          num=sum_e_i_2/n;          GCV=num/den;  %Delaney and Chatterjee answer          % the Golub method          A=h;          I_mat=eye(n);          num2=(I_mat-A);     num2=num2*y;     num2=norm(num2);     num2=num2^2;     num2=num2/n;               den2=(I_mat-A);     den2=trace(den2);     den2=den2/n;     den2=den2^2;          GCV2=num2/den2;               end


title: How could I confirm if the time series is stationary?

text: I'm using KPSS method to check if the time series is stationary, but I would also like to use another test to **confirm** the hypothesis that the time series is stationary. I'm already using Unit root tests as ADF or PP, but I need a more specific test to check the stationarity. What method should I use?


title: Device age vs. downtime relationship

text: **The Problem:** Retailer Y uses (for client operation purposes) a network of about 3000 devices. The whole network has 10 different models (from distinct manufacturers), aged from 0 years old (purchased in the current year) to 12 years old. Each device may suffer from malfunctions for which immediate repair is requested. The total time of unavailability for each device is accounted at the end of month and we call it downtime. The same device may have no malfunction during that month period or may have one or more unavailability periods during that month period. Common sense says that older equipment should have higher downtimes. **The Question:** (a) How can I check whether the Age variable has effective impact over the Downtime variable? Which statistical model should I use? (b) How can I account for the fact that each type of equipment should behave differently to aging regarding downtime? Thanks, Eduardo


title: Fitting a multivariate ARIMA model with R

text: I have 4 correlated time series, and I want to predict one of them, from the other 3. There is a clear seasonal effect in the 4 time series, so my first thought was to fit a multivariate ARIMA model, but I can not seem to find an R-function for this.


title: Decision boundaries and Gaussian density functions

text: This is for my hw, and if anyone can solve the first part of the question it will be great. Here is the question: Assume a two-class problem with equal a priori class probabilities and Gaussian class-conditional densities as follows: $$p(x\mid w_1) = {\cal N}\left(\begin{bmatrix} 0 \\\ 0 \end{bmatrix},\begin{bmatrix} a & c \\\ c & b \end{bmatrix}\right)\quad\text{and}\quad p(x\mid w_2) = {\cal N}\left(\begin{bmatrix} d \\\ e \end{bmatrix},\begin{bmatrix} 1 & 0 \\\ 0 & 1 \end{bmatrix}\right)$$ where $ab-c^2=1$.   1. Find the equation of the decision boundary between these two classes in terms of the given parameters, after choosing a logarithmic discriminant function.   2. Determine the constraints on the values of a, b, c, d and e, such that the resulting discriminant function results with a linear decision boundary.


title: Estimators - OLS & other

text: Question: ![Problem](http://i.stack.imgur.com/V2FVy.png) I am assigned this problem, but I don't quite understand what it's referring to. Is $\alpha$ the same as the intercept of the regression, and the coefficient on the variable u hat is 1? Would the OLS estimator for $\alpha$ come from the normal equations, then? Why do we need to square u hat in part b - what could this be related to? I'm not sure how I can compute anything just given this information. Thank you.


title: Best practice when analysing pre-post treatment-control designs

text: Imagine the following common design:   * 100 participants are randomly allocated to either a treatment or a control group   * the dependent variable is numeric and measured pre- and post- treatment Three obvious options for analysing such data are:   * Test the group by time interaction effect in mixed ANOVA   * Do an ANCOVA with condition as the IV and the pre- measure as the covariate and post measure as the DV   * Do a t-test with condition as the IV and pre-post change scores as the DV **Question:**   * What is the best way to analyse such data?    * Are there reasons to prefer one approach over another?


title: Which test should I use for this?

text: I have two groups (A and B), observed on day 0 and day 7. On each group there are 10 observations. So 2 x 2 x 10 = 40 results on my dataset. I want to find out if the control factor between group A and B had any effects over the week. Which statistical test should I use ?


title: Pre-post test with control group

text: > **Possible Duplicate:**   >  Best practice when analysing pre-post treatment-control designs I'm doing an action research study with 2 of my classes: 1 is the 'treatment' class and the other is the 'control'. Both classes undertook a pre-test prior to the 'treatment'. Then, after the 'treatment' stage (the control class didn't receive any 'treatment'), both classes did a post-test. What statistics should I use? Paired t-tests? ANOVA with repeated measure? Help!


title: Methods for analyzing pre/post test data

text: What methods of analysis are utilized for pre/post test data when outcome are knowledge scores? In addition what methods of analysis would be used for Likert behavioral score data?


title: still trying in R with CI and predictions

text: I looked over the tutorial in R for prediction and decided to reword my question. I need the prediction on average for male with a 95%CI then repeat with maximal values using the variables status, income and verbal. Now If I did the predictions correctly shown below with the mean then the max values. But doesn't that give 2 prediction intervals and 2 confident intervals for both the average and max. values                g2<-data.frame(status=75, income=15, verbal=10, sex=0)      g1<-data.frame(status=43, income=4.64, verbal=6.66, sex=0)      ie..Average               > g1<-data.frame(status=43, income=4.64, verbal=6.66, sex=0)     > predict(g,g1, interval='confidence', level=.90)     fit lwr upr     1 28.11506 19.7603 36.46983     > predict(g,g1, interval='prediction', level=.90)     fit lwr upr     1 28.11506 -10.95281 67.18293      Max Values used:               > predict(g,g2, interval='confidence', level=.90)            fit      lwr      upr     1 71.30794 47.07516 95.54072     > predict(g,g2, interval='prediction', level=.90)            fit      lwr      upr     1 71.30794 26.10037 116.5155      **Do** I need both the prediction and confidence interval? **Question:** how do I plot the CI to determine which is wider?


title: Confidence intervals for proportions (prevalence)

text: I have some data on a sample of $n=1776$ hospitals. For each hospital there is a total number of patients (`patients`), and a number of patients diagnosed with a particular condition (`diagnosed`). Do I take the mean of this proportion,`diagnosed`/`patients`, for all hospitals in the sample, $\hat{\mu}$, and calculate a 95% confidence interval as $\hat{\mu} \pm 1.96\sigma / \sqrt{n}$ or as $\hat{\mu} \pm 1.96 \sqrt{\hat{\mu}(1-\hat{\mu})/n}$ ? Or.... ? **Update** [Following comments from whuber]. Additionally, the data are broken down into 2 age groups (young and old) and 3 risk scores. That is, all 1776 hospitals have total numbers of patients as follows:                              younger patients       older patients                       Low risk            A                      D          Medium risk         B                      E          High risk           C                      F      ...and similarly for the numbers of patients with the condition. So, for each combination of age group and risk score, I would like to estimate the mean prevalence and a confidence interval for it. Here is some summary of the data               Risk   age    mean   sd      n     1      u50    0.37   0.19    1776     2      u50    0.49   0.25    1776     3      u50    0.54   0.26    1776     1      o50    0.45   0.36    1776     2      o50    0.52   0.42    1776     3      o50    0.67   0.41    1776


title: Has "Understanding Monte Carlo Markov Chain" by Roberts and Tweedie ever been published?

text: In Monte Carlo Statistical Methods by Robert and Casella, I saw there is a reference to Understanding Monte Carlo Markov Chain (Springer Series in Statistics) by Gareth O. Roberts and Richard L. Tweedie published in 2004. I searched this book on the internet and didn't find its availability. On Amazon, one link says it was published in 2008 and another link says it will be in 2015. Thanks!


title: Robust regularized regression

text: I've been using elastic net implemented in R (via glmnet) for some modeling, but I was wondering, due to the number of outliers in my data, if there was some sort of modeling approach for regularized robust regression? e.g. something like elastic net applied to robust regression. if there's something already in R, even better. Just curious to know what's out there. Thanks!


title: How to interpret p-values for factors

text: I have fitted a Poisson model to my data in R which includes two factors as independent variables. Each factor has 5 levels. I have used the `contrasts = list(FactorA='contr.sum', FactorB='contr.sum')`command to change the constraints of my model. So, the parameter estimates for each factor sum to 1. My problem is that I find it very difficult to understand what exactly the p-values of the estimates mean. For instance, if the p-value for one level of `FactorA` is significant, does this mean that it's generally significant for explaining the response variable or only in comparison to the other levels of the same factor? Are the constraints that I used involved in this interpretation that I have to make for the significance of each level? Also, I do not get any output for the last level of each factor. Is there any way to interpret those two levels?


title: CART (rpart) balanced vs. unbalanced dataset

text: I am fitting a tree (CART) to the olives-dataset. The training data has 436 observations (test data: 136). I have 3 responses (the 'Region' variable) which splits the training data into 116 / 74 / 246 observations. If I plot the variables eicosenoic and linoleic, I can see an almost perfect classification. I used a balanced dataset with 74 observations for each response (btw, is that correct or should I use a smaller size than 74 observations?) and got almost the same prediction results of the testdata as for the unbalanced dataset. That is why I am wondering if a balanced dataset is required in this case? I assume that balancing is not requried but I am not sure and would like to know other opinions.


title: Resources about forecasting stock returns with SVR

text: I'm currently working on SVR to predict stock market returns and **I would like to know if anyone of you could give me some advise regarding interesting papers on it.** These are the articles I'm working with:   1. Cao, L., Tay, F.E.H. (2001). Financial forecasting using support vector machines. Neural Computing and Applications, 10:184-192.   2. Wang,L., Zhu, J. (2010). Financial market forecasting using a two-step kernel learning method for the support vector regression. Annals of Operations Research, 174:103-120. [DOI] The second paper is really interesting but quite hard to put into practise as the author is using a two-step kernel learning method to find the optimal kernel function for the SVR. All the papers I found does not provide many details about tuning parameters (C and sigma for the Gaussian kernel). I'm using R for my computations (`Kernlab` and `e 1071` packages).


title: Statistical test to compare multiple treatments on the same 'individual' over time

text: I'm looking for a statistical test or tests to compare the effects of different 'treatments' on a single 'individual' over time. The single 'individual' could be one person or a group acting as a single entity (e.g. a team, a company, etc.). A simple example is comparing the effects of two different exercise regimes on one person's weight loss. The main concern I have with using a 'standard' statistical test (in this case, something like a t-test, Mann-Whitney U test or Kolmogorov-Smirnov Test seems appropriate) is the obvious non-independence of the data. I could check for independence by running some form of ABAB study (i.e. alternating the regimes several times to check that they are not interacting). This seems to be the approach that the author of this 'Self Experimentation' paper (Roberts, S.) seems to take. He then uses repeated chi-squared tests (first A against first B, second A against first B, etc.). Is it a valid statistical approach to perform multiple tests on the same data or does this increase the chance of type I errors as indicated in this paper? > The more often one analyses the accumulating data, the greater the chance of > eventually and wrongly detecting a difference, thereby drawing incorrect > conclusions from the trial. Should I worry about the non-independence? Can I simply perform multiple tests against the same data? Is there a specialised test for this kind of experiment?


title: Strategies for teaching the sampling distribution

text: **The tl;dr version** What successful strategies do you employ to teach the sampling distribution (of a sample mean, for example) at an introductory undergraduate level? _**The background_** In September I'll be teaching an introductory statistics course for second year social science (mainly political science and sociology) students using _The Basic Practice of Statistics_ by David Moore. It'll be the fifth time I've taught this course and one issue I've consistently had is that **the students have really struggled with the notion of the sampling distribution**. It's covered as the background for inference and follows a basic introduction to probability with which they don't seem to have trouble after some initial hiccups (and by basic, I mean **basic** \-- after all, many of these students have been self-selected into a specific course stream because they were trying to avoid anything with even a vague hint of "math"). I would guess that probably 60% leave the course with no to minimal understanding, about 25% understand the principle but not the connections to other concepts, and the remaining 15% fully understand. _**The main issue_** The trouble students seem to have is with the application. It's difficult to explain what the precise issue is other than to say they just don't get it. From a poll I conducted last semester and from exam responses, I think that part of the difficulty is confusion between two related and similar sounding phrases (sampling distribution and sample distribution), so I've don't use the phrase "sample distribution" anymore, but surely this is something that, while confusing at first, is easily grasped with a little effort and anyway it can't explain the general confusion of the concept of a sampling distribution. (I realize that it might be _me_ and my teaching that's at issue here! However I think ignoring that uncomfortable possibility is reasonable to do since _some_ students do seem to get it and overall everybody seems to do quite well...) _**What I've tried_** I had to argue with the undergraduate administrator in our department to introduce mandatory sessions in the computer lab thinking that repeated demonstrations might be helpful (before I started teaching this course there was no computing involved). While I think this helps overall understanding of the course material in general, I don't think it's helped with this specific topic. One idea I've had is to simply not teach it at all or to not give it much weight, a position advocated by some (e.g. Andrew Gelman). I don't find this particularly satisfying since it has the whiff of teaching to the lowest common denominator and more importantly denies strong and motivated students who want to learn more about statistical application from really understanding how important concepts work (not only the sampling distribution!). On the other hand, the median student does seem to grasp p-values for example, so maybe they don't need to understand the sampling distribution anyway. _**The question_** What strategies do you employ to teach the sampling distribution? I know there are materials and discussions available (e.g. here and here and this paper which opens a PDF file) but I'm just wondering if I can get some concrete examples of what works for people (or I guess even what doesn't work so I'll know not to try it!). My plan now, as I plan my course for September, is to follow Gelman's advice and "deemphasize" the sampling distribution. I'll teach it, but I'll assure the students that this is a sort of FYI-only topic and will not appear on an exam (except perhaps as a bonus question?!). However, I'm really interested in hearing other approaches people have used.


title: Determine range that contains given ratio of data points

text: I'd like to determine the "representative" range for a series of data points. For example, I'd like to know the range where 98% (or any other given ratio) of points are located. I could have used the inter-quantile range (1%-99% in this case), but in case of seriously skewed distributions this gives "too much attention" to long tails of the distribution. I'd better consider the following idea: imagine a histogram built for that series of points. Maximum of the histogram is considered to be 100%. Now draw the horisontal line at 1% mark and find the leftmost and rightmost intersections of the histogram curve with this line. Those intersections would give the left and right boundaries of the scale that contain most of the points. Does this approach have a commonly defined name? Is this approach appropriate for determining the "proper" range of the data?


title: What kinds of variables should we use the normality test for?

text: When should we apply normality tests? For which types of the variables should we apply the normality test? For example dependent variables, independent variables, or control variables, etc?


title: Summing up two matrices of distances among centroids

text: I am ecologist, using mainly the 'vegan' R package. I want to combine ('cbind') 2 matrices (sample x species) but samples (row) are different. I sampled within 24 sites, by using 2 distinct methods, fish assemblages (6 replicates) and macro-invertebrates assemblages (3 replicates). (See data below) I want to compute the pairwise distances among sites considering both assemblages(using braycurtis dissimilarity), let's say I want to get the overall distances among sites. I see 2 options: option 1, averaging at the site scale fishes and macro-invertebrates abundances , 'cbind' the two mean abundances matrix and calculating bray- curtis option 2, for each assemblage, computing bray-curtis and distances among sites centroids, and then summing up the 2 distance matrix. Please, could you tell me if the option 2 is correct and more appropriate than option 1. thank you in advance. Pierre                library(plyr)      library(vegan)          #assemblage 1: 15 fish species, 6 replicates per site             a1.env=data.frame(Habitat=paste("H",gl(2,12*6),sep=""),Site=paste("S",gl(24,6),sep=""),Replicate=rep(paste("R",1:6,sep=""),24))     summary(a1.env)          a1.bio=as.data.frame(replicate(15,rpois(144,sample(1:10,1))))     names(a1.bio)=paste("F",1:15,sep="")     a1.bio[1:72,]=2*a1.bio[1:72,]          #assemblage 2: 10 taxa of macro-invertebrates, 3 replicates per site     a2.env=a1.env[a1.env$Replicate%in%c("R1","R2","R3"),]     summary(a2.env)          a2.bio=as.data.frame(replicate(10,rpois(72,sample(10:100,1))))     names(a2.bio)=paste("I",1:10,sep="")     a2.bio[1:36,]=0.5*a2.bio[1:36,]               # "environmental" data at the sit scale     env=unique(a1.env[,c("Habitat","Site")])     env=env[order(env$Site),]          ######## OPTION 1, averaging abundances and cbind          a1.bio.mean=ddply(cbind(a1.bio,a1.env),.(Habitat,Site),numcolwise(mean))     a1.bio.mean=a1.bio.mean[order(a1.bio.mean$Site),]          a2.bio.mean=ddply(cbind(a2.bio,a2.env),.(Habitat,Site),numcolwise(mean))     a2.bio.mean=a2.bio.mean[order(a2.bio.mean$Site),]          bio.mean=cbind(a1.bio.mean[,-c(1:2)],a2.bio.mean[,-c(1:2)])          dist.mean=vegdist(sqrt(bio.mean),"bray")               #######OPTION 2, computing for each assemblage distance among centroids and summing the 2 distances matrix     a1.dist=vegdist(sqrt(a1.bio),"bray")     a1.coord.centroid=betadisper(a1.dist,a1.env$Site)$centroids     a1.dist.centroid=vegdist(a1.coord.centroid,"eucl")          a2.dist=vegdist(sqrt(a2.bio),"bray")     a2.coord.centroid=betadisper(a2.dist,a2.env$Site)$centroids     a2.dist.centroid=vegdist(a2.coord.centroid,"eucl")          dist.centroid=a1.dist.centroid+a2.dist.centroid     coord.centroid=cmdscale(dist.centroid,k=23) # for further distance-based analysis (is it correct ?)               ### COMPARING OPTION 1 AND 2     pco.mean=cmdscale(vegdist(sqrt(bio.mean),"bray"))     pco.centroid=cmdscale(dist.centroid)          comparison=procrustes(pco.centroid,pco.mean)     protest(pco.centroid,pco.mean)


title: Calculate probability for LibLinear classification results

text: I am using **LibLinear** for a document classification task, in which I would like to calculate the probability of correctness for each prediction. In fact, in the **LibLinear** , it does provide probability output for **logistic regression** , but not for default **support vector classification** task. Furthermore, based on 10-fold cross validation, the **logistic regression** is nearly 10% worse than **support vector classification**. So can anyone tell me, if I continue to use the **support vector classification** for the solution, is there a method to calculate the probability separately from the program? I greatly appreciate your help !


title: Expectation Maximization Clarification

text: I found very helpful tutorial regarding EM algorithm. The example and the picture from the tutorial is simply brilliant. ![enter image description here](http://i.stack.imgur.com/mj0nb.gif) Related question about calculating probabilities how does expectation maximization work? I have another question regarding how to connect the theory described in tutorial to the example. > During the E-step, EM chooses a function $g_t$ that lower bounds $\log > P(x;\Theta)$ everywhere, and for which $g_t( \hat{\Theta}^{(t)}) = \log P(x; > \hat{\Theta}^{(t)})$. So what the $g_t$ in our example, and it looks like it should be different for every iteration. In addition, in example $\hat{\Theta}_A^{(0)} = 0.6$ and $\hat{\Theta}_B^{(0)} = 0.5$ then applying them to the data we get that $\hat{\Theta}_A^{(1)} = 0.71$ and $\hat{\Theta}_B^{(1)} = 0.58$. Which is for me looks counter intuitive. We had some prior assumptions, applied it to the data and get new assumptions, so the data somehow changed the assumptions. I don't understand why $\hat{\Theta}^{(0)}$ doesn't equal to $\hat{\Theta}^{(1)}$. In addition, more questions emerge when you see Supplementary Note 1 to this tutorial. For example what is $Q(z)$ in our case. It's not clear to me why the inequality is tight when $Q(z)=P(z|x;\Theta)$ Thank you.


title: Laplace smoothing and Dirichlet prior

text: On the wikipedia article of Laplace smoothing (or additive smoothing), it is said that from a Bayesian point of view, > this corresponds to the expected value of the posterior distribution, using > a symmetric Dirichlet distribution with parameter $\alpha$ as a prior. I'm puzzled about how that is actually true. Could someone help me understand how those two things are equivalent? Thanks!


title: Excessive splitting in Random Forest

text: I've noted that quite a few of the node splits in my Random Trees are redundant since both child nodes have identical classifications. I suspect this is caused by an unfortunate choice of parameters, but which? I'm using the OpenCV implementation to generate the trees; for prediction I use my own code. That's how I discovered these redundant nodes. The OpenCV trainings parameters are:   * max_depth = 20    * min_sample_count = 7   * nactive_vars = 2   * max_num_trees = 250   * forest_accuracy = 0   * term_crit = CV_TERMCRIT_ITER (i.e. 250 trees) I have +/- 50.000 samples, 25 variables per sample and 7 classes. No missing variables, no categorical variables. The problem with redundant nodes appears with nodes much shallow than `max_depth` so that's probably not relevant.


title: How to compute for the verification rate at 0.001 FAR?

text: I want to test my pair matching method based on this evaluation metric but unsure how to do this? I only have positive (matching) test samples.


title: Multiple correlation: PCA and behavior

text: This is follow-up from a previous question (How to separate groups using PCA?). I have 25 normals and 12 patients.For each of them I have a vector representing a spectrogram (2000). So i have a matrix Z = [25x2000;12x2000]. I calculate [coeffZ, score, latent, tsquared, explained, mu]=pca(Z); bar(explained) shows me that the 1st 3 PCs explain most of the variance. I also have a behavioral score (from testing) for each of the subjects [37x1]. It was suggested that I see if the 1st 3 PCs can predict the behavioral score using multiple correlation coefficient. Specifically this - http://en.wikipedia.org/wiki/Multiple_correlation Does this make sense? Does anybody have an idea of how I can implement this in Matlab?


title: Distribution of dot products between two random unit vectors in $\mathbb{R}^D$

text: If $\mathbf{x}$ and $\mathbf{y}$ are two [ **Update:** independent] random unit vectors in $\mathbb{R}^D$, what is the distribution of dot products between them in the limit of large $D$? I guess it quickly (?) becomes almost normal (does it?) with zero mean and variance decreasing in higher dimensions $$\lim_{D\to\infty}\sigma^2(D) \to 0,$$ but is there an explicit approximate formula for $\sigma^2(D)$? **Update** Well, I ran some quick simulations. First, generating 10000 pairs of random unit vectors for $D=1000$ it is easy to see that the distribution of their dot products is perfectly Gaussian (in fact it is quite Gaussian already for $D=100$), see the subplot on the left. Second, for each $D$ ranging from 1 to 10000 (with increasing steps) I generated 1000 pairs and computed the variance. Log-log plot is shown on the right, and it is clear that the formula is very well approximated by $1/D$. Note that for $D=1$ and $D=2$ this formula even gives _exact_ results (but I am not sure what happens later). ![dot products between random unit vectors](http://i.stack.imgur.com/hC3FM.png) Why?


title: Why are all regression predictors in balanced multi-way ANOVA orthogonal?

text: As far as I understand, in fully balanced ANOVA (with any number of factors) understood as general linear model all categorical predictors (all main effects and all interactions of all orders) are mutually orthogonal, which leads to many nice properties like unique decomposition of variance. Orthogonality is completely obvious for 1-way ANOVA with several levels (from the way dummy coding works), and _sort of_ clear for higher-way ANOVAs as well, but I would like to see either a formal proof, or some intuition, or, ideally, both. I consulted Rutherford, 2001, "Introducing Anova and Ancova, GLM approach" (pdf), but weirdly could not find this discussed.


title: Intuition about a joint entropy

text: I am having trouble building some intuition about joint entropy. $H(X,Y)$ = uncertainty in the joint distribution $p(x,y)$; $H(X)$ = uncertainty in $p_x(x)$; $H(Y)$ = uncertainty in $p_y(y)$. If H(X) is high then the distribution is more uncertain and if you know the outcome of such a distribution then you have more information! So H(X) also quantifies information. Now we can show $H(X,Y) \leq H(X) + H(Y)$ But if you know $p(x,y)$ you can get $p_x(x)$ and $p_y(y)$ so in some sense $p(x,y)$ has more information than both $p_x(x)$ and $p_y(y)$, so shouldn't the uncertainty related to p(x,y) be more that the sum of the individual uncertainties?


title: Does the optimal number of trees in a random forest depend on the number of predictors?

text: Can someone explain why we need a large number of trees in random forest when the number of predictors is large? How can we determine the optimal number of trees?


title: Fitting a time-varying coefficient DLM

text: I want to fit a DLM with time-varying coefficients, i.e. an extension to the usual linear regression, $y_t = \theta_1 + \theta_2x_2$. I have a predictor ($x_2$) and a response variable ($y_t$), marine & inland annual fish catches respectively from 1950 - 2011. I want the DLM regression model to follow, $y_t = \theta_{t,1} + \theta_{t,2}x_t$ where the system evolution equation is $\theta_t = G_t \theta_{t-1}$ from page 43 of Dynamic Linear Models With R by Petris et al. Some coding here,               fishdata <- read.csv("http://dl.dropbox.com/s/4w0utkqdhqribl4/fishdata.csv", header=T)     x <- fishdata$marinefao         y <- fishdata$inlandfao          lmodel <- lm(y ~ x)     summary(lmodel)     plot(x, y)     abline(lmodel)      Clearly time-varying coefficients of the regression model are more appropriate here. I follow his example from pages 121 - 125 and want to apply this to my own data. This is the coding from the example               ############ PAGE 123     require(dlm)          capm <- read.table("http://shazam.econ.ubc.ca/intro/P.txt", header=T)     capm.ts <- ts(capm, start = c(1978, 1), frequency = 12)     colnames(capm)     plot(capm.ts)     IBM <- capm.ts[, "IBM"]  - capm.ts[, "RKFREE"]     x <- capm.ts[, "MARKET"] - capm.ts[, "RKFREE"]     x     plot(x)     outLM <- lm(IBM ~ x)     outLM$coef         acf(outLM$res)     qqnorm(outLM$res)         sig <- var(outLM$res)     sig          mod <- dlmModReg(x,dV = sig, m0 = c(0, 1.5), C0 = diag(c(1e+07, 1)))     outF <- dlmFilter(IBM, mod)     outF$m         plot(outF$m)     outF$m[ 1 + length(IBM), ]          ########## PAGES 124-125     buildCapm <- function(u){       dlmModReg(x, dV = exp(u[1]), dW = exp(u[2:3]))     }          outMLE <- dlmMLE(IBM, parm = rep(0,3), buildCapm)     exp(outMLE$par)         outMLE         outMLE$value     mod <- buildCapm(outMLE$par)         outS <- dlmSmooth(IBM, mod)         plot(dropFirst(outS$s))     outS$s      I want to be able to plot the smoothing estimates `plot(dropFirst(outS$s))` for my own data, which I'm having trouble executing. **UPDATE** I can now produce these plots but I don't think they are correct.               fishdata <- read.csv("http://dl.dropbox.com/s/4w0utkqdhqribl4/fishdata.csv", header=T)     x <- as.numeric(fishdata$marinefao)         y <- as.numeric(fishdata$inlandfao)     xts <- ts(x, start=c(1950,1), frequency=1)     xts     yts <- ts(y, start=c(1950,1), frequency=1)     yts          lmodel <- lm(yts ~ xts)     #################################################     require(dlm)         buildCapm <- function(u){       dlmModReg(xts, dV = exp(u[1]), dW = exp(u[2:3]))     }          outMLE <- dlmMLE(yts, parm = rep(0,3), buildCapm)     exp(outMLE$par)             outMLE$value     mod <- buildCapm(outMLE$par)             outS <- dlmSmooth(yts, mod)             plot(dropFirst(outS$s))          > summary(outS$s); lmodel$coef            V1              V2             Min.   :87.67   Min.   :1.445        1st Qu.:87.67   1st Qu.:1.924        Median :87.67   Median :3.803        Mean   :87.67   Mean   :4.084        3rd Qu.:87.67   3rd Qu.:6.244        Max.   :87.67   Max.   :7.853        (Intercept)          xts      273858.30308      1.22505       The intercept smoothing estimate (V1) is far from the lm regression coefficient. I assume they should be nearer to each other.


title: Evaluating the clustering of a Kohonen UMatrix

text: Given a converged Kohonen feature map, how would one evaluate the clustering in terms of intra- and inter-cluster distances? Assuming that both the trained codebook vectors and Unified Distance Matrix are available. Currently I have been evaluating them empirically using heatmaps of the UMatrix. However, I was wondering if there was a more reliable, statistical method.


title: How to deal with missing categorical data in logistic regression models?

text: I am participating in a Project with data from a complex survey. We are going to analyze data from a national fertility survey. Some of the questions in the questionnaire were only asked by a subsample. For example, a question about the job was only asked to employed people. This variable and others like that are important in our study and we would like to include it in our logistic model. What approach do you recommend? I have seen in some papers of this area to set up another category “else”, that is, they combine other people (to which the question was not raised) creating an additional category. If this approach is acceptable, should i attributed random values to the income variable in order to reduce bias? I will be very grateful if you could help me.


title: How to perform hypothesis testing for comparing different classifiers

text: I am trying to classify a small dataset (around 500 records) into two classes. I used various methods like SVM, Naive Bayes and k-nn classifier. Now, I would like to set the results from one of the classifiers are my baseline and perform a statistical hypothesis testing. I am not much familiar with this field of statistical testing, and I wonder how to proceed on this. I have been thinking of setting the SVM classifier as my baseline, but I am not sure how to perform a t-test (or similar) on the data. The input dataset has 10 attributes. Should I use the classification results from two classifiers and do a paired t-test on them? For example, I could take the result from Naive Bayes and perform the paired t-test with the SVM results (which is the baseline). Is this the right approach? Also, I am confused with the explanation for null and alternate hypothesis. Could someone be willing to give an idea about how to fix the null and alternate hypothesis.


title: Specification of logical node (with distribution?) in WINBUGS

text: For a piece of homework I have an assignment using WINBUGS which I must admit confuses me to say the least. Tangential to my question but I have a few stochastic nodes that are to be gamma distributed. Taking one example, the parameters of the gamma distributions (call `a` and `b`) also feed in to lambda which we know the prior distribution for.               d ~ gamma(a, b)     z = a/b ~ gamma(0.01, 0.01)     w = b/a ~ gamma(0.01, 0.01)      At first look I thought `z` and `w` should be specified as a logical node but then I can't add the distribution information (gamma). My code runs along these lines (although I've tried numerous versions):               d ~  dgamma(a, b)     z ~  dgamma(0.01, 0.01)     w ~  dgamma(0.01, 0.01)     a <- z*b     b <- a/z     a <- b/w     b <- a*w      But then I get errors along the lines of "multiple definitions of node...", if I remove one of each of the assignments for `a` and `b` to resolve the multiple definitions, WINBUGS stops responding. Would really appreciate any words of wisdom of where I'm going wrong.


title: Under what conditions should Likert scales be used as nominal data?

text: Similar to this question, Under what conditions should Likert scales be used as ordinal or interval data?, I am curious under what conditions should/could Likert scales be used as nominal data?


title: Cholesky decomposition and confidence ellipsoid

text: I'm trying to construct an error ellipsoid from a covariance matrix (which exists for a 3D point) and then sample consistent xyz points in this region. (This question succeeds this one.) What I'm currently doing is: 1) Calculate the Cholesky decomposition of the covariance matrix.   2) Sample each initial vertex point as a Gaussian with width 1 to generate (x', y', z')   3) Multiply (x',y',z') by the Cholesky decomposition matrix for the newly generated point.   4) Add this result to a matrix of the mean values. To add some more concrete details, the covariance matrix of the initial point is: $ \begin{pmatrix} 10.0115 & -10.6835 & 5.18024 \\\ -10.6835 & 11.4009 & -5.52798 \\\ 5.18024 & -5.52798 & 2.77646 \\\ \end{pmatrix} $ I then calculate the Cholesky decomposition, which is: $ \begin{pmatrix} 3.164095& 0& 0\\\ -3.376478& 0.017131& 0\\\ 1.637195& -0.001619& 0.309921\\\ \end{pmatrix} $ And the initial xyz point is at (35.5361, -37.2661, 22.521). Generating three random numbers from a Gaussian of mean 0 and width 1 yielded: -0.377495, -0.933623, 0.241011, and at this point I'm not completely clear on the correct procedure. My assumption is that it is correct to multiply the Cholesky decomposition matrix by a matrix containing the randomly generated Gaussian points, and then add this to a matrix containing the initial values. This would look like this: $ \begin{pmatrix} 35.5361& 0& 0\\\ 0& -37.2661& 0\\\ 0& 0& 22.521\\\ \end{pmatrix}+ \begin{pmatrix} 3.164095& 0& 0\\\ -3.376478& 0.017131& 0\\\ 1.637195& -0.001619& 0.309921\\\ \end{pmatrix} \begin{pmatrix} -0.377495& 0& 0\\\ 0 & -0.933623 & 0 \\\ 0 & 0 & 0.241011\\\ \end{pmatrix} $ For which the result is: $ \begin{pmatrix} 34.3417�& 0& 0\\\ 0 & -37.281 & 0 \\\ 0 & 0 & 22.5957\\\ \end{pmatrix} $ And I could definitely believe that a new xyz point at (34.3417,-37.281,22.5957) could be consistent with what I want to generate, but I'm still not completely confident I have followed a valid procedure to generate this. Any comments would be much appreciated!


title: How to organise the variable names in R without messing up?

text: When I work in R, every time I create a new variable I try to name it simple so that it iss easy to type like `x1, x2, ...` or `a1, a2, ...` Sometimes I use temporary variables with names temp, and again I continue making variables `temp1, temp2, ...`. As a result finally I mess up everything to a point impossible to understand the data if looked at after some days. Is there any easy way or rule to follow in creating new variables and organizing them?


title: Issues with regression and prediction

text: I am having a lot of fun with regression analysis at the moment, and by fun I mean bashing myself repeatedly over the head. I have a set of 200 data points, by filtering on a property of interest, I end up with 153 points of use. I initially used these 153 points to generate a linear regression, with an excellent R${^2}$ and a plot of fitted vs actual variables of almost a perfect diagonal. Great! However, it was suggested that this might only be an internally predictive model (which as I understand it means the model fits the data, rather than the opposite). So, I then tried this: I randomly selected a sample of 100 of the 153 results, and built the same model, it still gave a relatively good fit. I then used the predict function in R to try to predict the outcome of the other 53 records. It did not go well. What I got was one of 2 things.   1. the predictions made no sense at all, not even on the same scale as the actual values.   2. most of the predictions made sense (although weren't very accurate) and one or two, were on an entirely different scale (orders of magnitude larger, or smaller). Since the model I am fitting has time as the response variable, it was suggested I use a Gamma fit regression instead of a plain old linear regression. I tried this and ended up essentially with the result. So, am I using R correctly, was Gamma a good choice for this? I'm pretty sure my data is good (non biased) so if I am unable to predict, despite the good model - does this mean my model is useless? I've been working on this for some weeks now, and it would be great if I could salvage something. The R commands I have used:               modelSet<-sample(1:nrow(myData),100)     modelData<-myData[modelSet,]     predictData<-myData[-modelSet,]          fit<-lm("time~(x1+x2+x3+x4+x5+x6)^3", data=modelData)     pred<-predict(fit, predictData)     plot(predictData$time, pred) <- gives a really not useful plot               fit2<-glm("time~(x1+x2+x3+x4+x5+x6)^3", data=modelData, family=Gamma) # tried with link=log too     pred2<-predict(fit2, predictData)     plot(predictData$time, pred2) <- gives an even less useful plot


title: Quantify possible model improvement given samples X in Y dimensions and Z% accuracy in a random forest model

text: I've been doing some Machine Learning as a hobby and am playing around with Random Forest models and basically trying to apply that towards predicting things in my own life. So I have a series of product data with, say, 60 dimensions (which PCAs down to 5 factors) over 400 data points. I generate a few models by basically sampling 200 random data points at a time, and end up being 60%ish accurate in guessing which price bin of 10 clusters ($0-$10, $10-$30, $300-$500, etc.) an item is on a market place. So what I want to know: how do I express the possible improvement of accuracy if given more data points? Or more useful: how many data points do I need to maintain a high, say 95% or arbitrary X, confidence level? (Edit: I reread some old college texts and literature about confidence levels but am not sure where specifically to go for my use case. I may have some bad assumptions about how to do this since I haven't mathed in years, so please let me know how I should rephrase/rethink this!)


title: augmenting test set

text: In training a neural network for an image classification task, is it ok to augment the test set by affine transforms when the dataset is small (~1000 images), rather than doing k-fold cross validation. Augmenting the training set is standard practice to reduce overfitting, but very costly with k-fold CV. But how about augmenting the test set too, then I can do standard CV. (appreciate papers references applying this practice)


title: What is the interval that relates to the mean as the equal tailed interval relates to the median and the highest density interval relates to the mode?

text: When summarizing a one dimensional continuous distribution (e.g. a posterior distribution) it is common to use either an equal tailed interval (aka quantile-based) or a highest density interval. The 95% equal tailed interval conceptually corresponds to the median in that when the coverage of the interval $\rightarrow 0\%$ the interval converges to the median. In the same way the highest density interval corresponds to the mode as the mode is the point of the highest density. But another popular point summary of a continuous distribution is the mean and my question is: **What is the interval that corresponds to the mean?** That is:   * What is that interval called?   * How is it defined/calculated? If someone has a comment on the conundrum why it is the case that the mean is a really popular way of summarizing a distribution while the corresponding interval is not as popular (as is my impression), I would be gratefull!


title: Kalman Filter forecasting

text: I am trying to finish my Masters Dissertation in the UK and have to submit it in less than 2 weeks. Part of it is to forecast a bivariate regression using Kalman Filtering. It would be based on the lines of Time-Varying Parameters(TVP) so that Kalman Filter allow the parameters to vary over the data sample. My equation is quite simple: yt = a + bx1t + cx2t where yt is the dependent variable. a = constant. x1 = independent variable. b = coefficient for x1. x2 = independent variable. c = coefficient for x2. I am trying to use Kalman filter to forecast yt in a time-varying parameter (TVP) format. I have gone through some of the packages and the articles but I am not sure how to fill in some of the arguments (FF, V, GG, W, m0, C0, dV etc.) for the **dlm** package commands (dlmFilter, dlmForecast). I understand that in the package **sspir** , **SS** and **kfilter** commands could be used for this. Would someone tell me step-by-step how to use these commands esp. the arguments, for my data please. This is for my masters dissertation and I am running short on time.


title: Probability theory books for self-study

text: Are there any good books that explain important concepts of probability theory like probability distribution functions and cumulative distribution functions? Please avoid referring books like "Mathematical Statistics and Data Analysis" by John Rice which start with simple permutation concepts and then suddenly (in 2nd chapter) take a leap assuming knowledge in real-analysis, multiple and surface integrals and begin describing CDFs and PDFs and illustrating them in 3-dimensional figures. One is left scratching head as to how everything is connected. I am looking for self-study books and any book in the same category as "Calculus for the Practical Man" would be of great help?


title: Nonextensive statistics

text: I know that the Tsallis($S_q$) entropy is called nonextensive in the sense that if $P$ and $Q$ are two probability distributions then $S_q(P\times Q)=S_q(P)+S_q(Q)+(1-q)S_q(P)S_q(Q)$. Can someone explain what is meant by "nonextensive statistics" or "nonextensive statistical mechanics"? I searched a bit. I couldn't find anything which explains clearly. I hope someone in stat exchange can explain it in mathematical point of view.


title: Interpretation of Shapiro test

text: I'm pretty new to statistics and I need your help. I have a small sample looking as follows: > H4U 0.269 0.357 0.2 0.221 0.275 0.277 0.253 0.127 0.246 I run the Shapiro test using R: `shapiro.test(precisionH4U$H4U)` and I got the following result: > W = 0.9502, p-value = 0.6921 Now, if I assume the significance level at 0.05 than the p-value is larger then alpha (0.6921 > 0.05) and I cannot rejest the H0 hypothesis about the normal distribution, but does it allow me to say that the sample have normal distribution? Thanks!


title: Mixed joint distribution

text: Let $Y = \max\\{0,Y^*\\}$, where $Y^*$ follows a continuous distribution. I found in the paper which I am reading that $$ f_{Y,Y^*}(Y, Y^*) = I(Y = Y^*)f_{Y^*}(Y^*), $$ where $I(\cdot)$ is an indicator function. I do not know how to derive this expression. If I use the Bayes rule, I get $$ f_{Y,Y^*}(Y, Y^*) = f_{Y|Y^*}(Y|Y^*)f_{Y^*}(Y^*), $$ but I do not know how to deal with this object.


title: How to fix the threshold for statistical validity of p-values produced by ANOVAs?

text: I have run experiments on a group of users under two conditions, measuring the time it took users to finish their experiments. I used a cross-over design where half of the users started in the first conditions to end with the second, and the other half of the users did the other way around. I analyze the data provided in a few different ANOVAs and find different p-values for my hypotheses. Some are below 0.05, some are below 0.01, some are over 0.05. Do I need to fix an alpha level of statistical significance to be used in all my analysis, or can I report something like 'Hypothesis A is proven true at alpha level 0.05, while Hypothesis B is true at alpha level 0.01 (thus, possibly a stronger proof)'? I don't know if I am being clear enough here. Let me know and I'll add details if needed. Thanks.


title: Inference with one-sample Likert-type data

text: In case, one wants to do some inference with one-sample Likert-type data, what tests could one use? Signed Rank Test?


title: How many lags to use in the Ljung-Box test of a time series?

text: After an ARMA model is fit to a time series, it is common to check the residuals via the Ljung-Box portmanteau test (among other tests). The Ljung- Box test returns a p value. It has a parameter, _h_ , which is the number of lags to be tested. Some texts recommend using _h_ =20; others recommend using _h_ =ln(n); most do not say what _h_ to use. Rather than using a single value for _h_ , suppose that I do the Ljung-Box test for all _h_ <50, and then pick the _h_ which gives the minimum p value. Is that approach reasonable? What are the advantages and disadvantages? (One obvious disadvantage is increased computation time, but that is not a problem here.) Is there literature on this? To elaborate slightly.... If the test gives p>0.05 for all _h_ , then obviously the time series (residuals) pass the test. My question concerns how to interpret the test if p<0.05 for some values of _h_ and not for other values.


Is it possible to do meta analysis of only two studies. What will be limitation of such analysis.


title: Calculating % unsampled in sampling with replacement

text: You sample `N` of `N` items **_with_** replacement. How do you calculate the expected percent not sampled from original population `N`? * * * **Extra Credit:** Generalize to sampling **`k`** of `N` items **_with_** replacement.


title: Sample Size Problem with Tukey HSD Test

text: I was planning to use a Tukey HSD test for the following data, as it was used in a study with a similar approach. I regressed stock data for three time ranges (1 year daily data, 2 years weekly data, 8 years monthly data from today) and got the (beta and) alpha values. Dependent variable: stock returns of S&P 500; independent variable: always the returns of one of 13 hotel real estate investment trusts. Beta is calculated out of the REIT stock return and the S&P500, which is used to calculate Jensen's alpha (as per CAPM). After that I conduct a one-way ANOVA of the alphas (three alphas in total, one for each time period) and THEN eventually the Tukey HSD test. It is supposed to show whether REITs performed significantly better in one of the periods. I am unsure what my sample size is for the test specifically (or if it is equal), because the number of data points I have in the three time ranges to calculate the alphas differs, although the number of companies obviously stays the same. Does this makes sense?


title: Handling large data sets in R -- tutorials, best practices, etc

text: I am an R noob who is required to do various kinds of analysis on large data sets in R. So while looking around this site and elsewhere, it appeared to me that there are a lot of esoteric and less well known issues involved here - like which package to use when, what transformations to (not) apply on the data etc. I am just wondering if there is a book/tutorial/guide that demystifies all this and presents the information in a systematic way? I prefer doing this instead of looking around, and collating information from different sources online. Thanks in advance.


title: Logistic models as variables in a logistic model

text: I have a binomial logistic model I'm working on, with roughly 6 significant variables. I'd like one of those variables to be how likely something else is to happen. Lets say that I'm primarily interested in the probability that cars get sold from a dealers lots, with sub-interest in whether the dealers are golfers. The fact that they are/aren't golfers could affect whether they're there when the right buyer comes along, and thus it could significantly affect the probability of the car being sold. In its own right, the sub-model might be interesting, but more important is that if I can get it to predict being a golfer well enough, and add it in the main model, it can significantly improve the prediction power of the main model. Can something like this be done, or are there problems with using probabilities as a variables? Should I use the logit from the sub-model instead? One other question is what would happen if I used variables in the sub-model, and then again in the main model? Perhaps what I'm trying to do has a specific name and I'm just not searching for it correctly.


title: Asymmetric confidence intervals on bootstrap estimates

text: I've performed bootstrapping on my leastsq parameters and now I have a load of data from which I can get the mean and standard deviation for each parameter. Lovely. But when I look at a histogram of this data, I find (I suppose unsurprisingly) that the data is skewed ![histogram](http://i.stack.imgur.com/wN5me.png) (The line is the mean) How do I report the clearly asymmetric variance and get the proper centre (in this ~0.75)?


title: Using lme with a fixed beta (slope), and estimating the intercept only

text: lme.1.combo <- lme(ComboRate ~ p_w, random = ~1 | Rat,data=x)      The line above will return a fitted intercept term, and a fitted beta (slope) term given these variables (allowing for random intercepts for each rat). What I would like to do, is run the same regression, but with a FIXED slope. That is, if I keep the slope at a given value, what is the best fitting intercept going to be? So my output would be the best intercept value, given the slope value I assigned. Is there a way to do this using `lme` or any other function in R?


title: Given loads of data, can we always model it with polynomials?

text: Given Taylor series and enough data so as to not risk over-fitting, do you actually need to think about if your phenomenon is following an exponential, quadratic, logarithmic, ..., behaviour? I'm sure you can come up with counter- examples of when this isn't a good idea, but if we're very general and pragmatic, is it usually 'fine' to just fit the data to some degree nine polynomial hoping whatever bizzare pattern is hiding in the data will be very well approximated by its power series?


title: Steps required to choose a suitable model for Poisson regression

text: I have data on claim frequency and several predictors like age, gender of insured person and characteristic of vehicle. I have checked for the assumption of mean equal to variance of a Poisson model. Now I want to test for multicollinearity and homocedasticity. How am I supposed to proceed so that I choose a suitable model and fit it using the Poisson regression model?


title: t-test correction for serial correlation

text: I have a time series of serially correlated data (differences between spatial means of meteorological variables). I calculated the lag-1 autocorrelation (rho), and corrected the number of degrees of freedom for my t-test by (1-abs(rho))/(1+abs(rho)). When I calculate the standard error, do I divide by N, or by the adjusted number of degrees of freedom?