title: What is normality?

text: In many different statistical methods there is an "assumption of normality". What is "normality" and how do I know if there is normality?


title: Determine if a QTL interval region significantly overlaps with differntially expresed genes

text: I am trying to do an analysis of QTL regions and differentially expressed genes in Arabidopsis, and the trait I am interested in is hybrid vigor. I currently have 9 QTL regions (each containing many genes) associated with biomass heterosis, which were identified from a previous study. I also have a list of genes that are differentially expressed (DEGs) in hybrids. I then tested for overlap between the 9 separate gene lists (from the 9 genomic regions) and my list of DEGs. I see that there is significant overlap between the DEGs and 2 out of 9 lists of genes, using the hypergeometric test. I am wondering if there is a better way to determine if the overlap between the lists/regions is significant. I have seen in some QTL papers that a type of permutation or randomization test can be done, and I am wondering if that could be applied to my situation. The problem is that I am unfamiliar with these type of permutation tests, and am not sure what it entails. Any suggestions?


title: Sum of weights in portfolio theory is not equal to 1

text: I'm trying to understand basic portfolio theory using R. As far as I understood, the sum of the weights of assets must be equal to 1 . But in this link, that teaches how to compute the efficient frontier of a portfolio, the final sum of the weights is not equal to 1. What this means? Link of tutorial: http://economistatlarge.com/portfolio-theory/r-optimized- portfolio Weights final result: ![enter image description here](http://i.stack.imgur.com/rUc1k.png)


title: Repeated measures

text: I have administered an instrument (4 subscales, total items 32, no. of items in each subscale 12, 7, 9 & 4) three times over an year. i have worked out the factor scores. I am not sure I should compare means or go GENERAL LiNEAR MODEL, then REPEATED MEASURE. I need to know which way to go and why, also the steps in SPSS


title: Capitalization of n for sample size

text: Should the n for sample size be capitalized? Is there a difference between n and N?


title: How to incorporate prior knowledge in GPML?

text: I am using the MATLAB code for Rasmussen & Williams' book _Gaussian Processes for Machine Learning_. How can one incorporate prior knowledge in Gaussian process regression? Say, that the variance in one dimension of a two dimensional vector is greater. Is it only by considering the parameters of a normal distribution, or can it be more detailed?


title: Correcting standard errors when the independent variables are autocorrelated

text: I have a question about how to correct standard errors when the independent variable has correlation. In a simple time series setting we can use Newey- West covariance matrix with a bunch of lags and that will take care of the problem of correlation in the residuals. What does one do in a panel data setting? Imagine the situation where you observe firms over time: $$ Y_{i,t} = a + b\Delta{X_{i,t}} + \epsilon_{i,t} $$ where the $\Delta{X_{i,t}} = X_{i,t} - X_{i,t-n}$. It seems that clustering standard errors on $i$ and on $t$ should fix this problem. Am I correct?


title: Cross Validation folds in LibSVM

text: how the folds in cross validation of libSVM are done?, I mean, the splits are done aleatory? for example if I have 100 data and it is 4-fold then when it is an aleatory folding then can be splited like this?   * 20-45 : fold1   * 45-70 : fold2   * 70-95 : fold3    * 95-20 : fold4 or it make the folds always starting from the first one:   * 1-25 : fold1    * 25-50 : fold2    * 50-75 : fold3    * 75-100 : fold4 so how can I choose which kind on folding I want to use? Thank you for your answers


title: Z test of coefficients of output?

text: I have fitted a certain ARMA model and got the following output: ![sp](http://i.stack.imgur.com/TedV8.png) As you can see, some coeff were set to zero and the focus is on the ma5 coeff and the intercept. To be exactly, the ma5 coeff is `-0.05077979` and the intercept is `-0.00044509`. Also the exact s.e. are `0.02078887` (ma5) and `0.00032599` (intercept). I want to calculate the p-values: From the answer with 7 thumbs up on this thread (from mpiktas) I calculated it with:               (1-pnorm(abs(-0.05077979)/sqrt(0.02078887^2)))/2      which should give me the p-value for the ma5 coeff, and the output is               0.003645037      For the intercept I do               (1-pnorm(abs(-0.00044509)/sqrt(0.00032599^2)))/2      and get                0.04303588      Now I control this via coeff test of the lmtest package:               library(lmtest)     coeftest(myarmamodel)      and I get ![ot](http://i.stack.imgur.com/kZFax.png) As you can see, the estimates and the s.e. are the same, but the p-value is different. My question is: Why? What is wrong? Is the p-value of the lmtest wrong or the solution in the other post?


title: Predicting next character based on few previous -- how to combine predictions?

text: I asked this question at math.stackexchange.com first, but nobody answered. Perhaps statisticians can help me. The question is like this: I have a signal generator, which each second generates one of the letters 'a', 'b' or 'c'. I don't know anything else about this signal generator, but I suspect that there are some patterns to it. The signal generator is started at time 0. Given output of the signal from time 0 to time n, I need to forecast its output at time n+1. I suspect that there are some patterns to the signal, so I create a frequency table: For each sub-string (up to length 7) of the first n symbols of the signal generated data, I calculate three values: the number of occurrences of 'a','b', and 'c' after that sub-string. So, for example for sub-string "abc"(if it exists in the data), I store:   * the number of cases when symbol 'a' comes after string "abc"   * the number of cases when symbol 'b' comes after string "abc"   * the number of cases when symbol 'c' comes after string "abc" Now that I have all these data, I have 7 predictions as to what the next symbol could be. If for example the last 7 characters of the signal are "acbccba", then: If I look just at the one-character frequency table, then I will have a certain prediction, which would look like: "Since the last character of the string is 'a', and since coming directly after character 'a' there were 40 cases of letter 'a', 25 cases of letter 'b', and 130 cases of letter 'c', I predict that the next character will be letter 'c'" Similarly for the last 2-letter ("ba"), 3-letter ("cba"), ... , 7-letter ("acbccba"). So in the end I have 7 predictions. The question is, how do I find which next character is actually the most probable for this signal generator? Different predictions are based on different sample sizes, so how do I combine them effectively?


title: probability of interest

text: The model that I created in R is: > fit <\- lm(hired ~ educ + exper + sex, data=data) what I am unsure of is how to fit to model to predict probability of interest where p = pr(hiring = 1). Any help would be appreciated thanks, Clay **Edit:** This is the computer output for what I have computed so far. I am unsure if this is even a step in the right direction to find the answer to this question. What I am trying to do is, Fit a logistic regression model to predict the probability of being hired using years of education, years of experience and sex of job applicants.                > test<-glm(hired ~ educ + exper + sex, data=data, family=binomial(link="logit"))      > summary(test)           Call:      glm(formula = hired ~ educ + exper + sex, family = binomial(link = "logit"),           data = data)           Deviance Residuals:           Min       1Q   Median       3Q      Max        -1.4380  -0.4573  -0.1009   0.1294   2.1804             Coefficients:                  Estimate Std. Error z value Pr(>|z|)        (Intercept) -14.2483     6.0805  -2.343   0.0191 *      educ          1.1549     0.6023   1.917   0.0552 .      exper         0.9098     0.4293   2.119   0.0341 *      sex           5.6037     2.6028   2.153   0.0313 *      ---      Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1            (Dispersion parameter for binomial family taken to be 1)               Null deviance: 35.165  on 27  degrees of freedom      Residual deviance: 14.735  on 24  degrees of freedom      AIC: 22.735           Number of Fisher Scoring iterations: 7


title: Good Internet resource for preparation for studying mathematical statistics

text: I'm going back to university this September and I picked up a course which has some prerequisites > Mathematics & Statistics (Sci): Sample space, events, conditional > probability, independence of events, Bayes' Theorem. Basic combinatorial > probability, random variables, discrete and continuous univariate and > multivariate distributions. Independence of random variables. Inequalities, > weak law of large numbers, central limit theorem. I've studied statistics many years ago but as I don't use it in my work I completely forgot all of it. Since there're still a month before September, I think I should make some preparation for my study if I don't want to drop the course. So is there any good Internet studying resources that I can use to refresh my statistic knowledge? (I thought about buying textbooks, but they are really expensive)


title: Detrended line still has a trend

text: I have a line and I have detrended it but it seems to still have a trend. Visually it has one (see jpg) and it has a unit root (Dickey-Fuller) and fails the Durbin-Watson test, so it is not stationary. The line that I am calling 'detrended' is a plot of the residuals (errors) of a regression of the data on time. My question is why would a supposedly detrended line still have a trend? ![data](http://i.stack.imgur.com/TQ3wF.jpg)![data detrend](http://i.stack.imgur.com/UsTtn.jpg)


title: Text Classification in R

text: New to R, and am trying to do text classification. I am using R package tm to convert raw txt data into matrix. Here's the relevant code snippet.                col <- Corpus(DirSource(path),                     readerControl = list(reader = readPlain,                     language = "en",                     load = TRUE)))           tdm <-  TermDocumentMatrix(col, control = list(tokenize = NGramTokenizer))      I have the following questions: 1) Feature selection I need to do chi-squared or information gain based feature selection on my data. Which R packages can I look at? I came across at caret and boruta but they do not seem to be appropriate for what I am wanting to do. 2) Handling new (unseen) instances Let's say I have trained my model using my training set. When the test set comes in, I would need to pass it through same filters (stemming, stopword removal, tf-idf weighting, feature selection etc.). I have no idea how to do this ! Any hint/help will be much appreciated. Thanks in advance.


title: Are t tests of coefficients in multiple regression post hoc tests?

text: In multiple regression, if a global F test is significant, then are t tests (or Wald tests) for the coefficients considered to be multiple comparisons and post hoc tests and should they be adjusted?


title: Choosing the best featureset for prediction

text: I have this N sets of features F each with $F_i$ number of features. All the feature sets have 20000 examples and we have 20,000 labels for them. Lets say feature set $F_1$ has 10 features and feature set $F_2$ has 20 features. Each of these feature sets have 20000 examples. Now all these examples have labels. So I want to decide which of these feature sets is useful for the classification. I don't want to create a classifier and decide because my labels are highly unbalanced. I was thinking of doing this with canonical correlation. I will take one feature set $F_i$ and the label vector and do cannonical correlation, get the score. Do that for all the feature sets and then sort them based upon the scores obtained. How would that be?


title: Picking the smallest number of i. i. d. exponential variables to satisfy a condition

text: _The lifetime of a brand of lamp is exponentially distributed with parameter $\lambda=0.01\;h^{-1}$. How many lamps are needed so that there's a 95% probability that a room is lit for 1,000 hours, if_ _a) when a lamp burns, the next one is used?_ _b) all lamps are turned on simultaneously?_ * * * a) The lifetimes $X_1,…,X_n$ are independent, and I know that the sum of i. i. d. exponential variables has a gamma distribution. $$X_1+...+X_n\triangleq Y\overset\perp\sim \mathrm{Gamma}(n,\lambda)\\\ \mathbb P(Y\ge 1000) \ge 0.95\\\ \mathbb P(2\lambda Y\gt 2000\lambda) \ge 0.95\\\ 2\lambda Y\triangleq T\sim \mathrm{Gamma}(n,0.5)\equiv \mathrm{\chi^2}(2n)\\\ \mathbb P(T\gt 20) \ge 0.95\quad,$$ and a chi-squared table (and the restriction that $2n$ is even) gives $n=16$. Is there a way to answer this without a table (or just a normal table)? * * * b) Since the previous item uses one lamp at a time, I assume one is enough to light the room. Therefore, I need the largest lifetime to be 1,000 h or more. $$\max\\{X_1,...,X_n\\}\triangleq Z\\\ \begin{align}\mathbb P(Z\ge1000) &= 1-\mathbb P(Z\le1000) =\\\ &=1 - \mathbb P(X_1\le1000,...,X_n\le1000)\overset{\mathrm{i.i.d.}}=\\\ &=1 - [\mathrm{cdf_{Exp(0.01)}}(1000)]^n = 1 - (1-e^{-0.01\cdot1000})^n \end{align}\\\ 1 - (1-e^{-10})^n \ge 0.95\\\ (1-e^{-10})^n \le 0.05\\\ n\log(1-e^{-10}) \le \log 0.05\\\ n\log[(1-e^{-10})^{-1}] \ge \log 20\\\ 4.5401\cdot10^{-5}\cdot n\ge 2.9957\\\ n\ge 65983.9\quad,$$ and so I need 65,984 lamps, a suspiciouly large number, even though I want one of them to last ten times its average lifetime. Is something wrong?


title: How to set the radius value in self organizing map?

text: I'm training the self organizig map, I need to set the value for the radius of it. is there any method to find the optimum radius size ?


title: "Exposure" and "Outcome" can be switched without changing the odds ratio?

text: I'm enrolled in an introductory course in statistics. In one of my lecture notes, the instructor explains why Odds Ratio is popular and he mentions: "Exposure" and "Outcome" can be switched without changing the odds ratio. Could anyone explain what that means? I appreciate your help.


title: assumption GARCH (1,1)

text: We have a problem with parameters’ estimation of GARCH(1,1). Our return equation is $$ r_t = \mu + h_t^{1/2} \cdot z_t, $$ where $z_t$ has normal distribution and $h_t$ is the conditional variance and residuals $\epsilon_t = h_t^{1/2} z_t$. For GARCH(1,1) we have $$ h_{t+1} = \omega + (\alpha + \beta)h_t + \beta(r^2_t - h_t), $$ and we couldn’t explain why $r_t^2 $correspond to $\epsilon^2_t$ if $v_t=\epsilon^2_t - h_t$ for the ARMA assumption. We usually find in general for GARCH (1,1) $$ h_t = \omega + (\alpha + \beta) h_{t−1} + \alpha(\epsilon^2_{t−1} – h_{t−1}). $$


title: R plugin for SPSS

text: Can anyone tell me how can I integrate the R plugin for SPSS version 17?? I need to do some statistics using the R plugin for SPSS, unfortunately I can't find the plugin on the internet, so can anyone help me


title: Negative binomial distribution mixture model with R

text: I have two data vectors of observed count data: $A$ and $B$, where count $A_n$ and $B_n$ refer to the same observation point. $A$ is assumed to follow a negative binomial distribution. $B$ is assumed to be the result from two underlying processes. For one again an independent negative binomial distribution and additionally a kind of shadowing effect from $A$, where a certain fraction $f$ of counts from $A$ are also observed in $B$. Accordingly, one can simulate such data by:               set.seed(13579)     f <- runif(1)     A <- rnbinom(100, mu = 100, size = 1)     B <- floor(f*A) + rnbinom(100, mu = 20, size = 1)      **Now to my question: Since $B$ is a mixture model of two negative binomial distributed variables, how can I estimate the value of factor $f$ explaining the underlying data best?** The final result should be a vector $B^\prime$, correcting $B$ for its shadow portion: $B^\prime = B - (A/f)$. Since I know $A$ and $B$, to gain this I need to estimate $f$ from the data and the underlying assumption of a negative binomial distribution. I'm neither a statistics nor an R expert, therefore any help is highly appreciated.


title: What would be a good way to work with a large data set in Excel?

text: I have a large data set which is in .dbf format right now and what I would like to do is be able to manipulate it easily in Excel and do something like subtotal and calculate stdev and ratios. Details of the data set; This data set contains shopper information. It has 1.2 million rows and 20 columns where the rows are each a unique shopper and the columns hold their shopping data (what they bought). I am using Office 2007 programs, I know Excel the best but was wondering what alternatives I could use to accomplish my goals (subtotal, calculate stdev, and ratio's).


title: How can I use credibility intervals in Bayesian logistic regression?

text: I'm using Gibbs sampling to learn the distributions of coefficients for a multinomial logistic regression model. At the end, I end up using the mean values of distributions of coefficients, and the resulting logistic regression is used as a classifier. I'm trying to find out advantages of having probability distributions for coefficients and the response variable, but I can't really see the way to leverage credibility intervals. What can I do with these distributions that I end up with, other than using their mean values? ok, I think I've failed to express my question clearly. Here is the update Let's assume that I have y = b1*x1 + b2*x2 + b3*x3 in my hands. for all b1,b2,b3 I have normal distributions, with mean and quantiles. These distributions come from the mcmc results. I can use the means of b1,b2,b3 and that'd give me a classifier. I can also use %2.5 quantile values and 97.5 quantile values for b1,b2,b3 and so on. What would be the probabilistic interpretation of these other equations? Can I produce a smoother classifier this way, rather than using only means? The use of credibility intervals is quite clear to me when they are used for a single variable, but in this case, I'm talking about N variables (5 actually) each with their own credibility intervals. I'm having trouble getting the semantics of this setup. I have not seen any papers etc that discusses this, and any pointers would be appreciated.


title: Addition of Probability Intersects

text: Good day everyone, I am currently looking at a self-help excercise which goes about exhibiting the following behavior which really puzzled me. \begin{equation} > P(D) = P(A \cap D) + P(B \cap D) + P(C \cap D) \end{equation} After drawing out the Venn diagram, it just doesnt seem to make sense. ![enter image description here](http://i.stack.imgur.com/at76u.jpg) The original question that prompted my question as follows: ![enter image description here](http://i.stack.imgur.com/aKqzJ.jpg)


title: Equivalence of AIC and p-values in model selection

text: In a comment to the answer of this question, it was stated that using AIC in model selection was equivalent to using a p-value of 0.154. I tried it in R, where I used a "backward" subset selection algorithm to throw out variables from a full specification. First, by sequentially throwing out the variable with the highest p-value and stopping when all p-values are below 0.154 and, secondly, by dropping the variable which results in lowest AIC when removed until no improvement can be made. It turned out that they give roughly the same results when I use a p-value of 0.154 as threshold. Is this actually true? If so, does anyone know why or can refer to a source which explains it? P.S. I couldn't ask the person commenting or write a comment, because just signed up. I am aware that this is not the most suitable approach to model selection and inference etc.


title: Is there an ANOVA table generalization for two nested CART models?

text: The last line is an example of what I'm looking for:               data(airquality)     attach(airquality)          lm1 <- lm(Ozone ~ Solar.R+Wind)     lm2 <- lm(Ozone ~ Solar.R+Wind+Temp)     anova(lm1 , lm2)          require(rpart)     rp1 <- rpart(Ozone ~ Solar.R+Wind)     rp2 <- rpart(Ozone ~ Solar.R+Wind+Temp)     anova(rp1 , rp2) # this doesn't exist - is there something like it? some sort of anova.rpart function?      Thanks!


title: Explaining regression performance differences

text: **Context** I have a regression framework and two sets of data. Using leave-one-out cross- validation, the first set gives very good performance and the second set gives rather poor performance. I need to explain the reason for this difference in performance. Having looked at the data, it is clear that the first set is a much easier test case. The data in the first set does not vary as much as in the second set. In addition, the test and training data of the first set are strikingly similar, which is not the case in the second dataset. To explain the reason for the performance difference, I can show these differences in data qualitatively using pictures of the physical phenomena being modeled. However, I would also like to quantify these dataset differences and relate it to the performance discrepancy. **Questions** What metrics can I use to show that the first dataset is a trivial test case and that the second set is more challenging? As a statistician, what else would you want to know in order to be confident that the performance discrepancy was caused by these differences in data and not something else?


title: Optimal blind poll construction

text: This is a question strongly related to Cauchy "characters". I'm constructing a 4 question canvassing questionnaire that will tell the likely voter being contacted which of the presidential candidates most closely matches them. The advantage of this approach for a dark horse presidential candidate is obvious, presuming, of course, that most of the likely voters match him. I do have the verbal interest in this from the state-level executive director for such a presidential candidate. One might be entitled to think that this work has been done umpteen times by the thousands of political science PhDs and/or major polling organizations -- at least using the General Social Survey data if nothing else -- and one would be wrong. Moreover, we don't have much time to deploy. Ideally the questionnaire construction would result in a kind of decision tree where the door-to-door canvassing volunteer could have a mobile device app providing the next question to ask based on the answers to prior questions. Also ideally, the construction process, itself, would minimize the likely- voter contact as this drives the expense. Using GSS data zeros out that cost and would be optimal if we could get access to the raw GSS data, but we can't. We have to do a survey to gather the data for construction of the 4-deep tree of questions. On the result side, as a practical compromise, I've proposed falling back to finding just 4 questions rather than finding a 4-deep tree of questions. On the construction side, as a practical compromise, I've proposed a prize- fund backing a tournament where:   1. The contestants each submit 4 questions.    2. The submitted questionnaires are paired up for the contests.   3. Campaign volunteers each get a pair of contested questionnaires, and get ten likely voters to completely answer all 8 questions.   4. The winning questionnaire of a contest is the one whose author can, from the answers to his own questionnaire, best-guess the answers to the opponent's questionnaire.   5. Award prizes after the log2(N) contests have selected a winner of that tournament.   6. Publish the rankings of the questionnaires and, by permission, their authors. As resources permit, this tournament is iterated for multiple rounds. We really have to place weight on the value of up-front volunteer time, so minimizing the construction labor is crucial. I know this is fairly far from a pure mathematics question but I've brought it as close as I can to some kind of weighted figure of merit involving high- value labor in the construction phase and the expected accuracy of the resulting 4 questions answered by likely voters contacted during blind poll canvassing by lower value labor. The question: About how far from optimal would be the proposed practical compromise of the 4-question questionnaire (constructed as described) from the ideal 4-deep decision tree of questions constructed from an infinite number of samples during the construction phase? A secondary question: Is there a better way to make use of the same up-front volunteer time?


title: Need algorithm to compute relative likelihood that data are sample from normal vs lognormal distribution

text: Let's say you have a set of values, and you want to know if it is more likely that they were sampled from a Gaussian (normal) distribution or sampled from a lognormal distribution? Of course, ideally you'd know something about the population or about the sources of experimental error, so would have additional information useful to answering the question. But here, assume we only have a set of numbers and no other information. Which is more likely: sampling from a Gaussian or sampling from a lognormal distribution? How much more likely? What I am hoping for is an algorithm to select between the two models, and hopefully quantify the relative likelihood of each.


title: Imputation in normalized signals

text: I'm currently analyzing a variety of signals. The problem I have is that I have several "missing" values. These "missing" values represent the absence of signal, they are not errors in sampling or something else. So it seems important to keep them when performing the data analysis. However, I'm unsure about the best way to include them in the analysis. I have devised 2 main options, both using standardized data:   1. To transform continuous data into categorical data using the standard deviations of the normalized data, and adding a category for "no-signal".   2. Replacing all missing values by a z-score of -4 or -6. I would like to ask about your opinion (pros and cons) of using one of these methods. PS: Maybe is worth to say that my data are right skewed ranging from ~ (-2) to (+5) SD.


title: Different hierarchical clustering results

text: I'm running a hierarchical clustering on a sample of data using the steps below:               library(RODBC)          setwd('D:/r/cluster2')     channel <- odbcConnectExcel('cluster.xls')     data <- sqlFetch(channel, 'clust9')          y9 <- data.frame(inf=data$infest, faible=data$faible, moyen=data$moyen, fort=data$fort, lon=data$Lon, lat=data$Lat)          y9 <- na.omit(y9) # listwise deletion of missing     y9.use <- y9     y9 <- scale(y9) # standardize variables          wss <- (nrow(y9)-1)*sum(apply(y9,2,var))     for (i in 2:15) wss[i] <- sum(kmeans(y9, centers=i)$withinss)          plot(1:15, wss, type="b", xlab="Number of Clusters", ylab="Within groups sum of squares")          # K-Means Cluster Analysis     fit <- kmeans(y9, 5) # 5 cluster solution          aggregate(y9,by=list(fit$cluster),FUN=mean)          y9 <- data.frame(y9, fit$cluster)          # Ward Hierarchical Clustering          d <- dist(y9, method = "euclidean") # distance matrix          fit <- hclust(d, method="ward")           plot(fit) # display dendogram          rect.hclust(fit, k=5, border="red")      and i got this results: ![enter image description here](http://i.stack.imgur.com/65YdS.png) But when I did the same steps the next day I got different results: ![enter image description here](http://i.stack.imgur.com/clS5h.png) They are not different in everything, but there are individual that they now belong to another cluster! so I don't know why this behavior? i'm interested in interpreting and explaning the results, so when i get different results each time, that will make my previous interpretation wrong, what can I do for now ?


title: Checking factor / covariate independence in ANCOVA

text: One of the ANCOVA assumptions, as I read it here, says that: _**Independent variables orthogonal to covariates**. In traditional ANCOVA, the independents area assumed to be orthogonal to the factors. If the covariate is influenced by the categorical independents, then the control adjustment ANCOVA makes on the dependent variable prior to assessing the effects of the categorical independents will be biased since some indirect effects of the independents will be removed from the dependent._ See also the third assumption here. Could you tell me how should I verify that this assumption is met, in a given model containing three factors and one covariate? You can refer to SPSS or R ways of doing this. Thank you very much


title: Regression with standard deviation dependent on predictor?

text: I got this question of how to deal with a regression model involving a set av predictors. The model wanted was one where the response should be divided by the (non-constant) standard deviation of the error, but this standard deviation is a simple linear function of one of the predictors (σ(x) = a+bx) so the response would be Y/σ(x). The person wo asked me had been suggested to set up a maximum likelihood estimation procedure, an iterative one. But even if that would be possible, it seemed as a problem principally simple enough to warrant the search for a solution which has already been found. Any suggestions? Robert


title: The variance-covariance matrix of the least squares parameter estimatimation

text: I'm learning Linear Regression for Regression from "The Elements of Statistical Learning". Why The variance-covariance matrix of the least squares parameter estimates is easily derived( from(3.6) and is given by   $$ Var(\hat{\beta}) = (X^TX)^{-1}\sigma2. $$   Typically one estimates the variance $\sigma^2$ by   $$ \hat{\sigma^2} = \frac{1}{N-p-1}\sum_{i=1}^N(y_i-\hat{y}_i)^2 $$ Could someone explain the above formulas in detail, like the following:   $$ E(\hat{\beta})=E((X^TX)^{-1}X^Ty) \\\ =(X^TX)^{-1}X^TE(y) \\\ =(X^TX)^{-1}X^TX\beta \\\ =\beta $$


title: How to find a fitted statistical model to a series of data?

text: There are several outputs of a numeric computations, I'm seeking a method to model them with a _general statistical model_ (e.g., distribution function). For example:               ...     0.044423     0.127713     0.010692     0.019980     0.013242     0.071450     0.054766     0.013618     0.067680     ...      Each series has different length. Intuitively, Statistica has an option such as "Distribution Fitting" which works well but only on a series of data at a time. It provides a graph and model as follows for each series:               Exponential, Chi-Square test = 72.94, df = 12 (adjusted) , p = 0.00      ![enter image description here](http://i.stack.imgur.com/ljy0e.png) There are at least 30 such a series. How to fit a general (best) model to each of them but finally select the best overall fitted distribution model? How to avoid doing the job manually? A short comment about the above output is also appreciated. How can the information be interpreted? * * * Here is what I recently did. I just put together all histograms of 30 series (black dots) and then used a kernel-based smoothing algorithm to find a curve (red curve). Now it should be easy to fit a known distribution to that curve (because it is only one). **Any comments about this solution?** ![enter image description here](http://i.stack.imgur.com/Bo5FD.jpg) ** please note the updated comments put on answers.


title: Can hazard ratios be undefined?

text: Can a hazard ratio be undefined because you are comparing a "lower value" to a "higher value?" For example, suppose you are interested in some variable $x \in [0, 100]$. Say you are interested in comparing 5 values of $x$. If you choose the highest value of $x$ to be the reference level... Would all the hazard ratios of the other values be undefined? **Edit.** Suppose you have the following groups: $A,B,C$ and $D$ with $D$ as the reference level. You are looking at the outcome death for each of these groups. If there are no outcome events in group $A$ for example, would the hazards ratio for group $A$ be $0$ or undefined?


title: Leave-one-out cross validation and boosted regression trees

text: Colleagues of mine recently presented a work where they calibrate boosted regression trees (BRT) models on small data sets ($n= 30$). They validated the models using leave-one-out cross validation (LOOCV) using R2, RMSPE and RPD indices. They also provided these indices computed by training and validating the model on the full dataset. The R2, RMSPE and RPD values obtained through LOOCV were almost strictly equal to the R2, RMSPE and RPD values obtained when validating on the training data set. My questions are :   * Is such a results expected for LOOCV on BRT?   * Is this because BRT is relatively insensitive to outliers (and to single individuals?) that excluding one individual during LOOCV does not make a difference, providing nearly similar calibrated models with same performance metrics on the excluded individuals?   * In that case does LOOCV for BRT makes any sense, compared to repeated k-fold CV with $k < n$? Thank you in advance


title: Bootstrap test for correlation coefficient

text: I want to test the hypothesis that correlation coefficient between X and Y is 0 with a bootstrap, however I don't know which is a correct way to construct bootstrap samples. I have several ideas, which may be wrong: 1) randomly permute Y among X's as in permutation test, but sample Y with replacement 2) permute Y among X's without replacement and sample pairs (X,Y) with replacement Could you please tell me the right way to do bootstrapping in this case and justify it? Thank a lot!


title: How to draw random samples from a negative binomial distribution in R?

text: Let $X$ have the Negative Binomial distribution with parameters $r$ and $p$. The Negative Binomial distribution is a mixture distribution or compound distribution. That is $X$ is $\text{Poisson}(\lambda)$ where $\lambda$ is randomly chosen from a $\text{Gamma}(r, p/(1-p))$. Use this relation to write an R function for randomly drawing $n$ Negative Binomial random samples. _Hint_ : Use the R functions `rpois` and `rgamma`.


title: Transforming the distance value from a center, to a probability value

text: Let $c_i$ be the center of a micro-cluster (i.e. we have many centers representing some fragments of clusters). Let $c_1$ be the center which is the closest to a new data-point $x$, such that $d_{x,c_1} = distance(x, c_1)$. (1) Is it possible to express d by a probability $d_{x,c_1}$ by a probability $p_{x,c_1}$ which tends to 1 as $d_{x,c_1}$ is higher, and to 0 as $d_{x,c_1}$ is smaller ? I think it is more intuitive to manipulate a probability instead of a distance. (2) Same question as (1), however this time I suppose that $p_{x,c_1} = 0$ if $d_{x,c_1} < radius$.


title: Efficacy of an intervention

text: I am new to the site so apologies if this (probably quite straightforward) question has been answered before- I have carried out research on an intervention with 2 groups (control and intervention), and collected pre and post data for a number of measures (stress, mood, etc.). To measure the efficacy of the intervention (or lack thereof), I think I need to use a mixed between-within ANOVA to compare the 2 groups on each of these dependent variables. Would I be correct in saying that? I'm fairly sure, but thought it best to check with those in the know! I plan to carry out further analyses but this is my main focus, so any help would be much appreciated.


title: ANOVA using day numbers

text: I am looking at differences in dates of maximum counts of geese at 3 sites in the UK. My data look like these created in R:               df <- data.frame(day.number=c(sample(70:75, 30, replace=T),                                    sample(77:82, 30, replace=T),                                    sample(80:84, 30, replace=T)),                      site=paste0('site', sort(rep(1:3, 30))))      The variable day.number is the day number on which the maximum count of geese was recorded in each year. So day number 70 would be 70 days from 1st January in that year. Each row represents a year and there are 30 years worth of data from each site. I am considering carrying out an ANOVA to test to see if the date of the maximum count is different among sites. Example R code:               summary(aov(day.number ~ site, data=df))      Is an ANOVA appropriate here? Is there a better way of looking at differences in day number of maximum count among sites?


title: Is a "split plot" ANOVA with two factors the same as two-way ANOVA with repeated measures in one factor?

text: Is a "split plot" ANOVA with two factors identical to two-way ANOVA with repeated measures in one factor? if not, what is the distinction?


title: Categorizing transportation modes with k-nearest neighbor?

text: I'm working on a classifier to sort out transportation modes based on certain attributes of an activity. I use a nearest neighbor algorithm like this : 2 example sets of training data (let's assume they are the only ones I have) This set represents a car ride:               Attribute  |  AvgSpeed   TopSpeed   GyroValue      ---------------------------------------------     Value      |  40         80         20      This set represents a walk:               Attribute  |  AvgSpeed   TopSpeed   GyroValue      ---------------------------------------------     Value      |  5          8          90      AvgSpeed, TopSpeed and GyroValue have different weights (e.g. 0.3, 0.3, 0.4) If I feed the program something like `A=25, B=60, C=40`, it should be able to predict a car ride. I'm currently using a k-neighbors classification system for this (more specifically the KNeighborsClassifier algorithm with scikit-learn in python). Would there be a stronger case for another type of ML algorithm here? I'm not convinced this is the best fit for my data, and this is not exactly my field of expertise. Thank you.


title: Normalizing Data in Multiple Categories

text: I couldn't decide if this question belonged here or on the Math stack exchange, but I'm posting it here since it seems more directly relevant. Basically, I have a bunch of categories, labeled, say, A-Z. Each of these categories has a bunch of "tasks" inside it. Tasks are specific to a category - they only go into one. For the purposes of abstracting the situation, we'll say there are infinite tasks. Each category has a different level of difficulty - a task in category A may take 30-40 seconds, whereas a task in category X might take 10-15 minutes. I have, for the sake of nice numbers, 100 users working in these categories. They choose a pseudo-random category (so you can't rely on their selection of category to normalize things) and may switch multiple times in a day. Some users will only work categories A-C, some D-F, some might work in all of them - it's completely inconsistent. What I am trying to do is compare these users to eachother, to find who is the best/fastest. However, I can't just look at who had the fastest tasks-per- hour, because they might have been doing those tasks in an easy category. I need to come up with a coefficient for each category to calculate a sort of "weighted total tasks"-per-hour for each user, in order to fairly rank them. I was previously doing something like: `(Total (all users) Completed Tasks in A/Total (all users) time spent in A)/(Total (all users) Time Spent in ALL categories/Total (all users) tasks completed in ALL categories)` I hope that makes sense. This gave a coefficient (with a mean of all coefficients at '1') for each category that I could use to normalize things. However, in examining this more closely, I don't think it's fair. If user 23 is really fast, but they only ever work in Category A, that will make everybody who works in category A look bad except for user 23 - user 23 drives the coefficient down. I think I need a way to normalize a user first to themselves, then normalize the categories coefficients using that user-normalized data. If that makes any sense at all - I'm hurting my OWN brain, here. Can anybody help me out, or point me in the right direction? PS: The categories difficulty is changing pseudo-randomly on a live basis, so benchmarking users in each category is a no-go.


title: Validate cluster analysis in R

text: I am trying to validate hierarchical cluster analysis result following a paper by Guy Brock, et al. _clValid: An R Package for Cluster Validation_ (pdf). Do I have to use all these methods? What are the most common validation methods? Also I am stuck at the following code:               emp <- as.matrix(emp)     library(clValid)     intern <- clValid(emp, 4:10, clMethods = c("hierarchical", "kmeans", "diana", "fanny",                                                 "som", "pam", "sota", "clara", "model"),                        validation = "internal")     summary(intern)          Error in if (xx != 0) xx/10 else z/10 : argument is of length zero     In addition: There were 50 or more warnings (use warnings() to see the first 50)      How to deal with the zeros in the matrix?


title: Multivariate Weighted Linear Regression

text: Very simple. I am looking for a package that does Multivariate Linear Regression with weights on the observations. Does anyone know of a package that does this? I am shocked that I have not been able to find any. **NOTE:** R does _NOT_ do multivariate regression. The `lm()` help page specifically states: "If response is a matrix a linear model is fitted separately by least-squares to each column of the matrix. " This means independent regression models for each response variable. Thus `lm()` does **NOT** do multivariate linear regression. It merely does several univariate linear regressions for convenience.


title: Pretest Posttest with 2 groups

text: Help! I have 2 groups (New Professionals and Experienced Professionals) would have received the same training. Both groups were given a survey before and after the training. I have not been given the survey, but I do have 18 variables provided in the data set as clues to what was is the survey. I am suppose to report on the usefulness of the survey and draw conclusions for future trainings. There are 59 total participants. There are 35 new professionals and 24 experienced professionals. Should I use a principal component factor analysis, then chi square to test the differences between groups or run MANOVA with the pretest as a covariate? Will the different number of professionals in the groups hurt my analysis?


title: Example of strong correlation coefficient with a high p value

text: I was wondering, is it possible to have a very strong correlation coefficient (say .9 or higher), with a high p value (say .25 or higher)? Here's an example of a low correlation coefficient, with a high p value:               set.seed(10)     y <- rnorm(100)     x <- rnorm(100)+.1*y     cor.test(x,y)      cor = 0.03908927, p=0.6994 High correlation coefficient, low p value:               y <- rnorm(100)     x <- rnorm(100)+2*y     cor.test(x,y)      cor = 0.8807809, p=2.2e-16 Low correlation coefficient, low p value:               y <- rnorm(100000)     x <- rnorm(100000)+.1*y     cor.test(x,y)      cor = 0.1035018, p=2.2e-16 High correlation coefficient, high p value: ???


title: Standard deviation of weighted sums of random distributions with weights random but fixed

text: Let $a,b,c,d$ be independent normally distributed random variables. I'm aware that the following distributions: $$c a + d b$$ $$(c+d)a$$ both have the same standard deviation ($=\sqrt{2}$ if $a,b,c,d$ have unit variance). I'm interested in the statistics generated by the following process: Let $a,b$ represent vectors containing a large number of samples, while $c,d$ each represent a single sample. Then what is the expected standard deviation of the two distributions above? Using the MATLAB code below, I get mean standard deviations of approximately $1.25$ and $1.125$, respectively.               N=100;     for m=1:100000         a = randn(N,1);         b = randn(N,1);         c = randn(1);         d = randn(1);              x(m)=std(c*a+d*b);         y(m)=std((c+d)*a);     end          mean(x)     mean(y)      If possible, I'd like to characterize the distribution of standard deviations fully, but a way of deriving the mean would be helpful too. My apologies if this question is too easy (I hope so!).


title: How to choose best value of 2 parameters simultaneously?

text: I have devised an algorithm that takes as input a dataset $s$ and two parameters, $k$ and $d$. $k$ is discrete while $d$ is categorical. My algorithm outputs a vector $v$ of length 20 which contains real numbers between 0 an 1. I run my algorithm on three different datasets $s_1$, $s_2$, $s_3$. For each dataset, I vary $k$ in the range {5, 10, 20, 30, 40, 50, 60}, and I vary $d$ from {1, 2, 3}. I run each combination and save the resulting $v$. My goal is to find out which values of $k$ and $d$ show the best average performance (in terms of $v$, larger values are better) across all three datasets. I'm ideally looking for a test or method that will tell me exactly which value of $k$ and $d$ are the best overall values for these three datasets. How can I do this?


title: How to combine p values?

text: I have three datasets A, B and C. A and C can construct a contingency table, and I can get a p value (pValue1) with Fisher Exact Test. B and C can construct another contingency table, and I can get another p value (pValue2) with Fisher Exact Test. I want to know how to combine these two p values and get a single one. One possible method is using Fisher method, but it tends to generate a new value close to the minimum one. I want to use some other method to tend to generate one close to the maximum one, namely both reject the null hypothesis. Can anyone tell me which method I can use? Thanks.


title: How is mean date calculated?

text: I was surprised to discover yesterday that various summary statistics can be calculated on dates. Googling has found lots on syntax required to calculate summary statistics on dates, but very little on how calculations are done. How are summary statistics (mean, sd etc) calculated on dates? Example R code:               myDates <- as.Date(paste0(sample(1:28, 100, replace=T), '/',                                sample(1:12, 100, replace=T), '/',                                2013),                         '%d/%m/%Y')          mean(myDates)     median(myDates)     sd(myDates)     var(myDates)


title: Can data processing time be estimated?

text: I'm curious if there are any analytical or empirical methods to estimate the time it would take to process a data set (grouping or various prediction algorithms). Known variables would obviously include the size of the dataset, but also the types of variables (continuous/discrete), the number of factors for discrete variables, the amount of memory available, and processor power. This sounds like a long shot to me, but I thought I'd throw it out there anyways.


title: How to do this test?

text: Suppose I have a sample $\\{(x_{1i}, x_{2i}, \dots, x_{mi}), i=1,\dots, n\\}$of $m$ unknown random variables $X_1, X_2,\dots, X_m$. How can I test if $X_1 = X_2 =\dots = X_m$? Furthermore, if there is a nonrandom explanatory variable $Y$ such that $Y=y_j$ for $X_j, j=1, \dots, m$, how can I test if $Y$ has an influence on the difference between $X_1, X_2,\dots, X_m$?


title: Is Heckman two-step overly sensitive to dummy variables?

text: I am doing a Heckman probit in Stata and whenever I try to include region dummy variables (which in normal probit are all significant), the model either doesn't converge, or does converge with rho ~ 1, which in my model makes no sense. So is it something about dummy variables that Heckman struggles to estimate? Thanks a lot!


title: KPCA in R proportion of variability explained

text: I'm using `kpca` function from `kernlab` and try to get the proportion of variance explained by each component as in standard pca. I don't select the number of features a priori since I would like to check their contribution. However, I get 124 components which is much more than my original dataset which has 10 covariates. On the other hand, there are no large eigenvalues so I can cut at some level. Is there any alternative to select a number of features? here are my eigenvalues:                    Comp.1      Comp.2      Comp.3      Comp.4      Comp.5      Comp.6      Comp.7      Comp.8      Comp.9     Comp.10     Comp.11      0.040155876 0.031142483 0.029499281 0.024417995 0.023562053 0.021718055 0.020310953 0.018984183 0.017789973 0.017311123 0.015484136          Comp.12     Comp.13     Comp.14     Comp.15     Comp.16     Comp.17     Comp.18     Comp.19     Comp.20     Comp.21     Comp.22      0.015346860 0.015005007 0.013791102 0.013291260 0.012670090 0.012180261 0.011882593 0.011523336 0.011107896 0.011045174 0.010477924          Comp.23     Comp.24     Comp.25     Comp.26     Comp.27     Comp.28     Comp.29     Comp.30     Comp.31     Comp.32     Comp.33      0.010251907 0.009882142 0.009606943 0.009529857 0.009362611 0.009340580 0.009062668 0.008987593 0.008699146 0.008670243 0.008549814          Comp.34     Comp.35     Comp.36     Comp.37     Comp.38     Comp.39     Comp.40     Comp.41     Comp.42     Comp.43     Comp.44      0.008398879 0.008214842 0.008091366 0.008029260 0.007924718 0.007857977 0.007771030 0.007691762 0.007657295 0.007582320 0.007510590          Comp.45     Comp.46     Comp.47     Comp.48     Comp.49     Comp.50     Comp.51     Comp.52     Comp.53     Comp.54     Comp.55      0.007470620 0.007404477 0.007285695 0.007246443 0.007134445 0.007087406 0.006956178 0.006935525 0.006898103 0.006864934 0.006653101          Comp.56     Comp.57     Comp.58     Comp.59     Comp.60     Comp.61     Comp.62     Comp.63     Comp.64     Comp.65     Comp.66      0.006605607 0.006585557 0.006513107 0.006395417 0.006207376 0.006171564 0.006043939 0.006023738 0.005955121 0.005894706 0.005788706          Comp.67     Comp.68     Comp.69     Comp.70     Comp.71     Comp.72     Comp.73     Comp.74     Comp.75     Comp.76     Comp.77      0.005714630 0.005700367 0.005601950 0.005550044 0.005441031 0.005410300 0.005367971 0.005246899 0.005161450 0.005093251 0.005026677          Comp.78     Comp.79     Comp.80     Comp.81     Comp.82     Comp.83     Comp.84     Comp.85     Comp.86     Comp.87     Comp.88      0.004984414 0.004866770 0.004674961 0.004655324 0.004644769 0.004529852 0.004447371 0.004411176 0.004338879 0.004258299 0.004135511          Comp.89     Comp.90     Comp.91     Comp.92     Comp.93     Comp.94     Comp.95     Comp.96     Comp.97     Comp.98     Comp.99      0.004078752 0.003985330 0.003902527 0.003838939 0.003734150 0.003582305 0.003547204 0.003485095 0.003440328 0.003397815 0.003363152         Comp.100    Comp.101    Comp.102    Comp.103    Comp.104    Comp.105    Comp.106    Comp.107    Comp.108    Comp.109    Comp.110      0.003246147 0.003223031 0.003184239 0.003091351 0.002938476 0.002868938 0.002765338 0.002645138 0.002572225 0.002544704 0.002466896         Comp.111    Comp.112    Comp.113    Comp.114    Comp.115    Comp.116    Comp.117    Comp.118    Comp.119    Comp.120    Comp.121      0.002419687 0.002298704 0.002187789 0.002089151 0.002019031 0.001957721 0.001908535 0.001887064 0.001760442 0.001705021 0.001587056         Comp.122    Comp.123    Comp.124      0.001536336 0.001228544 0.001079629


title: Non-parametric estimators for time-varying binomial proportion

text: I have a bunch of count data associated with time intervals (potentially overlapping and of variable lengths), say $(s_i, t_i, n_i, N_i)$ where $N_i$ is a count of the total number of events observed by a particular observer over a time interval $[s_i, t_i]$, and $n_i$ is the number of those events of some class of interest. I'm interested in modelling how the _proportion_ of interesting events varies over time, assuming that this proportion is independent of the observer. The total number of events observed in each interval may vary somewhat arbitrarily based on the observer, although if necessary we can safely assume that in aggregate over all observers it gives a reasonable picture of the overall event density. The simplest model I could think of for this data would be to assume $n_i \sim Binomial\left[N_i, \frac{1}{t_i-s_i}\int_{s_i}^{t_i}\theta(w) dw\right]$ For some unknown smooth function $\theta(t)$ which lies in $[0,1]$ for all $t$, and look for non-parametric methods to estimate $\theta$. A slightly better model could perhaps be to assume an inhomogeneous Poisson process for the overall events with rate $\lambda(t)$, and that the events of interest are a thinning of this process in a time-varying proportion $\theta(t)$ yielding another inhomog. Poisson process with rate $\lambda(t)\theta(t)$. We can further assume that both these processes are thinned further by each observer, independently to the "event of interest" thinning, in some unknown fixed proportion representing the proportion of events of both types which the observer sees. This proportion cancels out leaving: $n_i \sim Binomial\left[N_i, \frac{\int_{s_i}^{t_i}\theta(w)\lambda(w) dw}{\int_{s_i}^{t_i} \lambda(w) dw}\right]$. Here $\lambda(t)$ only matters up to a constant factor, so perhaps it would be better to model it as a density rather than a Poisson process rate? Anyway the simpler model falls out as a special case of this if you assume $\lambda(t)$ is constant/uniform. I was wondering if anyone has any pointers on how to estimate the proportion function $\theta(t)$ in this kind of model? A naive ad-hoc guess would be to compute a separate estimate for the overall rate and the rate of interest using kernel smoothing methods, then compute a ratio of these two estimates (perhaps adding on some small constants akin to what you'd get using a Beta prior to estimate a standalone binomial proportion e.g. $\frac{n+\alpha}{N+\alpha+\beta}$). To bring out the heavy guns, I suppose one could put a logistic Gaussian Process prior on $\theta$ and estimate look for a posterior estimate. I imagine this would require some fairly hairy approximate Bayesian methods. I was wondering there were any simpler frequentist estimators for this (based for example on kernel density estimators), and if anything is known about hyper-parameter selection, convergence rates, confidence intervals and so on. If it helps to consider a simpler problem first, we could perhaps assume that we have point observations for the events (of interest and not). Thanks!


title: How to compare distributions of two groups where each group has multiple observations from a small group?

text: I have data on continuous measurements (length of time of a behaviour) from two groups of individuals which differ in their phenotype (phenotype A or B, the variable used to group them). Each group is made up of a small number of individuals (4-5), and I have multiple (>100) measurements from every individual within each group. The number of individuals per group and the number of measurements per individual are not equal. How would I go about statistically comparing the mean and variances of these measurements between groups (phenotypes), accounting for the fact that each measurement is not independent (i.e. there are multiple observations from each individual within the groups, rather than each measurement being from a different individual? I basically want to know if the length of time spent in a particular behaviour differs between phenotype A and phenotype B - do the means differ and is one phenotype more consistent (have a tighter distribution) than the other. Many thanks in advance!


title: Neural network and dynamical system

text: Dynamical systemsare those whose evolution can be described by a rule, evolves with time and is deterministic. In this context can I say that Neural networks have a rule of evolution which is the activation function $f(\sum of product of weights and features)$ ? Are neural networks (A) dynamical system (B) Linear or Nonlinear dynamical systems? Can somebody please shed some light on this?


title: How to derive $P$ in continuous Kalman filter?

text: I want to derive the equation $\dot P = K'P+PK'^T + Cov(Kv-Gw)$ from the system $\dot x=Fx+Gw \ \ \ \ \ \ w \sim N(0,Q)$ $z=Hx+v \ \ \ \ \ \ v \sim N(0,R)$ given $\dot {\hat x}=K'\hat x+Kz$ $K'=F-KH$ $\dot {\hat x}=F\hat x+K(z-H\hat x)$ $d\tilde x/dt = K'\tilde x + Kv -Gw$ $\tilde x = \hat x - x$ $P(t)=E[\tilde x\tilde x^T]$ where $K'$ and $K$ are to be chosen to optimize the estimate $\hat x$. Here is how I tried (and failed) to show this: $P(t)=E[\tilde x\tilde x^T]$ $\dot P(t)=E[(K'\tilde x+Kv-Gw)\tilde x^T+\tilde x(K'\tilde x+Kv-Gw)^T]$ $=E[(K'\tilde x)\tilde x^T+\tilde x(K'\tilde x)^T]+E[(Kv-Gw)\tilde x^T+\tilde x(Kv-Gw)^T]$ $=K'P+PK'^T+E[(Kv-Gw)\tilde x^T+\tilde x(Kv-Gw)^T]$ I cannot figure out how to proceed further. Any help will be appreciated.


title: How to convert an unrestricted co-occurrence network to a restricted one?

text: I have built an unrestricted co-occurrence network of words from a songs corpus. To convert it to a restricted network, Ramon Ferrer Cancho and Ricard V. Sole describe the following approach in their paper The small world of human language: > The technique can be improved by choosing only pairs of consecutive words, > the mutual co-occurrence of which is larger than expected by chance. This > can be measured with the condition $p_{ij} > p_i*p_j$, which defines the > presence of correlations beyond that expected from a random ordering of > words. If a pair of words co-occurs less than expected when independence > between such words is assumed, the pair is considered to be spurious. Graphs > in which this condition is used will be called restricted (unrestricted > otherwise). Word pairs are considered as co-occurring if they occur within a maximum distance of 2 within a song. Here, $i$ and $j$ are the two words; $p_{ij}$ is the probability that the two words co-occur. If $i$ and $j$ are statistically independent, then the probability that they co-occur is given by the product $p_{i}\cdot p_{j}$. If they are not independent, and they have a tendency to co-occur, then $p_{ij}$ will be greater than $p_{i}\cdot p_{j}$. I am confused on how to calculate $p_{ij}$, $p_{i}$, and $p_{j}$. For example, to calculate $p_{i}$, can I take the count of all occurrences of $i$ divided by total number of songs in my corpus? Alternatively, $p_{i}$ can also be calculated by dividing count of all occurrences of $i$ by total number of words in the corpus. Which would be right approach out of the two? Similarly, to calculate $p_{ij}$, can I take the count of co-occurences of $i$ and $j$ divided by total number of words or total number of possible co- occurring word pairs? **EDIT:** Moreover, If $i$ and $j$ are statistically independent, then shouldn't the probability that they co-occur be given by $4∗p_{i}⋅p_{j}$ since word pairs are considered as co-occurring if they occur within a maximum distance of 2 within a song.


title: Is it possible to differentiate in expectation?

text: If $$E[f(x)]=0$$ can we derive that $$E[f'(x)]=0?$$ For example $f(x)$ means some noise with zero mean, gaussian distribution.


title: How to compute multidimensional omega with R

text: I need to perform a computation of reliability of a 5-point Likert scale having 6 items. From a factor analysis I found that my scale is a multidimensional scale (3 factors), so I cannot use Cronbach's alpha to compute reliability. I have seen in several papers that it is possible to use the multidimensional extension of the McDonald's omega. Does the `omega` function in the **psych** package allow to do this? Also, is there an R function to compute the Stratified alpha?


title: My data passed Johansen test, but they couldn't pass Phillips-Perron test. Why?

text: I want to choose two-pairs for pairtrading. For pairtradings, two pairs need to pass two tests, the Johansen test (for cointegration) and the P.P. test (for stationary). As I knew, if they related in cointegration, they would be stationary. But, My data passed Johansen test, but they couldn't pass P.P. test. Why does this result occur?


title: Forensics in wireless networks, anomaly detection and beyond?

text: first i'de like to apologize if this is not the right place. Next year i'm gonna be working on my final project in computer security, i have to build a wireless forensics tool that can analyse a data capture and scan it for attack possibilities or weird behavior, ( WEP attack, WPA, evil twin, attack on radius, DOS, NMAP Scan, and i'm open to any ideas/suggestions you have). My goal for in this project being :   1. Take on as many attacks as possible.   2. Get an optimal execution time. I was asked to do my research and come up with a plan for the project, so i thought it would be a good idea to include some Machine learning features and GPU computing techniques, i wonder if anyone here has ever worked on something similar o this, if so, some guidance would be nice, where to start ? any books or resources i should check ? how useful can this be ? what are the possibilities ? is building a neural network for attack detection a good idea or should i stick to the manual packet analysis method ? As to the programming, i'm gonna be using python with Wireshark and TCPDump for packet analysis and filtering. thank you


title: Expectation of multivariate function of random variables

text: How would I define a random variable that is:   * Some random value (uniform distribution) if the magnitude of the difference of the parameters is greater than some threshold $t$    * Otherwise the random variable takes on the value of the first parameter  I'm thinking of something like this: $$ g(x,y) = \left\\{ \begin{array}{l l} x & \quad |x-y| < t\\\ \text{Random real between $0$ and $N$} & \quad |x-y|>t \end{array} \right. $$ But I have no idea how to actually use this function, for example I want to calculate the expectation of $g(X,Y)$ but I can't figure out a way to do so using my current definition. Really my question is if it is possible to evaluate this integral given my definition of $g(X,Y)$: $$ E[g(X,Y)] = \int_{-\infty}^\infty \int_{-\infty}^\infty g(x,y)f_{X,Y}(x,y)\mathrm{d}x\mathrm{d}y $$ where $f_{X,Y}(x,y)$ is the joint PDF of two continuous random variables $X$ and $Y$. Any suggestions, or recommendations would be greatly appreciated!


title: EM algorithm for mixture of Negative Binomial distributions

text: I am trying to derive the EM-algorithm of mixtures of negative binomial distribution $Neg\;Bin(r,p)$. I have the updating equations for updating the E-step as well as $p$ and the mixing coefficients $\phi_j$ s. However, when I try to derive the M-step for $r$, I got stuck with an equation that involves the digamma function. $$\frac{dQ}{dr} = \sum_{i=1}^m w^{(i)} \left( \psi(x^{(i)} + r) - \psi(r) + log(p) \right) = 0$$ Note that $Q$ is conditional expectation of log-likelihood and $w^{(i)}$ comes form the E-step, and hence can be treated as a constant. And $\psi$ is the digamma function. The superscript of $w$ and $x$ are just the index of data. I read a paper here , which primarily deals with mixtures of gamma distribution. And it uses gradient ascent algorithm to estimate $\alpha$. If I follow the paper's approach to estimate $r$, then I will have the following gradient ascent algorithm: $$r^{(t+1)} = r^{(t)} + \frac{a_t}{n} \frac{dQ}{dr}$$ I am not sure if this is the right way to go about this problem. Does anyone what I can do solve for $r$ ?


title: Predicting Y using X for the following data

text: I have observed the following pattern of behavior in my data-set and I'm wondering how I can explore whether I can predict Y using a value of X. This plot is a 2D histogram of values. ![enter image description here](http://i.stack.imgur.com/PP95y.png) I have previously explored interquartile regression to identify the upper bound, but I'm wondering if there are any other approaches which I can apply which would enable me to identify the range of Y, depending on the X value. I assumed with heteroskedastic data like this that I would need to apply a transformation, but I'm not sure if this is the correct approach. Could I just fit an interquartile regression to the upper and lower bounds and then provide an estimated value (with an error) for any value of Y that is predicted? **EDIT:** Following some suggestions from @whuber and @jbowman I have included some more graphics of the dataset. These are plots of the mean and standard deviation binned by groups of `8` (e.g. 1 - 8, 9 - 16), and a histogram of slices through the data at 4 intervals. ![enter image description here](http://i.stack.imgur.com/CsByA.jpg) Histograms of slices through the data use a binwidth of 400. For example, for `1000`, I included all values of `X` that are 200 above or below so the slice ranges from `800 - 1200`: ![enter image description here](http://i.stack.imgur.com/O8sPf.jpg)


title: Evaluate deviation from negative binomial model

text: I'm trying to figure out how to determine to what extent a sample deviates from a negative binomial model fitted to a larger population. As an example, I generated counts of doctor visits for a population of males and females aged between 20 and 60:               library("MASS")     library("ggplot2")          set.seed(8)          n <- 1000     dispersion <- 1     baseline <- 1     ageeffect <- 0.02     gendereffect <- 0.5          gender <- sample(c("male", "female"), size=n, replace=TRUE)     age <- round(runif(n, min=20, max=60))     mu <- baseline + age * ageeffect + (gender=="female") * gendereffect     visits <- rnbinom(n, size=dispersion, mu=mu)     barplot(table(visits), space=0)          data <- data.frame(visits, gender, age)     summary(m <- glm.nb(visits ~ age + gender, data = data))          newdata <- data.frame(       age = rep(seq(20, 60), 2),       gender = c(rep("male", 41), rep("female", 41))     )          p <- predict(m, newdata, type="link", se.fit=TRUE)          newdata$visits <- exp(p$fit)     newdata$ll <- exp(p$fit - 1.96 * p$se.fit)         newdata$ul <- exp(p$fit + 1.96 * p$se.fit)          ggplot(newdata, aes(age, visits)) +       geom_ribbon(aes(ymin = ll, ymax = ul, fill = gender), alpha = .25) +       geom_line(aes(colour = gender), size = 2) +       labs(x = "Age", y = "Prediction")      Now I would like to quantify how much a specific group of people deviates from the model. Groups can be of different sizes and composition, and I'm only interested in groups which exhibit a above-normal number of doctor visits. The idea is to calculate some kind of scoring for each group (let's say the inhabitants of a specific neighbourhood) and follow this up over time. What would be the preferred way of doing this? I did calculate p-values for each observation in the sample, but I'm not sure how to proceed.               samplesize <- 10     s <- data[sample(1:n, samplesize),]     pp <- predict(m, s, type="link", se.fit=TRUE)     s$mu <- as.vector(exp(pp$fit))     s$theta <- rep(m$theta, samplesize)     s$p <-  1 - pnbinom(s$visits, m$theta, mu=s$mu)     s              visits gender age       mu     theta          p     487      0 female  48 2.385442 0.9899275 0.70307422     873      2 female  38 2.115261 0.9899275 0.31264813     19       3   male  29 1.440554 0.9899275 0.12164221     78       2 female  33 1.991872 0.9899275 0.29476639     362      1   male  43 1.704576 0.9899275 0.39648278     120      1 female  48 2.385442 0.9899275 0.49534441     658      1   male  53 1.922300 0.9899275 0.43182233     14       5 female  22 1.745163 0.9899275 0.06636807     663      2 female  57 2.657989 0.9899275 0.38294175     77       1   male  24 1.356523 0.9899275 0.33088457      **Edit 1:** After reading http://mikelove.wordpress.com/2012/03/12/combining- p-values-fishers-method-sum-of-p-values-binomial/, would Fisher's method be a good approach? **Edit 2:** I'm exploring Pearson residuals now, but something's not clear to me. In the code below I generate a dataset and construct a model. I then calculate Pearson statistics for samples of 20 and 1000 counts. From what I read I expected these statistics to follow a Chi-squared distribution with n-2 degrees of freedom, but that doesn't really seem to be the case. Am I doing something wrong?               library("MASS")          n <- 1000     dispersion <- 1     baseline <- 1     ageeffect <- 0.02     gendereffect <- 0.5          dataset <- function(n) {       gender <- sample(c("male", "female"), size=n, replace=TRUE)       age <- round(runif(n, min=20, max=60))       mu <- baseline + age * ageeffect + (gender=="female") * gendereffect       visits <- rnbinom(n, size=dispersion, mu=mu)       data <- data.frame(visits, gender, age)       return(data)       }          model <- function(data) {       m <- glm.nb(visits ~ age + gender, data = data)     }          pred <- function(data, m) {       p <- predict(m, data, type="link", se.fit=TRUE)       data$mu <- as.vector(exp(p$fit))       data$theta <- rep(m$theta, length(data$visits))           data$pr <- (data$visits - data$mu) / sqrt(data$mu + data$mu^2 / data$theta)       return(data)     }          pearson.stat <- function(data) {       return(sum(data$pr^2))     }          test.chisq <- function(data) {       return(pchisq(pearson.stat(data), df = length(data$pr) - 2, lower=FALSE))     }          # population of 100000 and model          dp <- dataset(100000)     m <- model(dp)          # pearson statistic (samples of 1000)          ps <- NULL     for (i in seq(1, 1000)) {       d <- dataset(1000)       p <- pred(d, m)       ps <- c(ps, pearson.stat(p))     }          # pearson statistic (samples of 20)          ps2 <- NULL     for (i in seq(1, 1000)) {       d <- dataset(20)       p <- pred(d, m)       ps2 <- c(ps2, pearson.stat(p))     }          # plot          dplot <- function(res, df) {       x <- seq(min(res), max(res), length.out=100)       y <- dchisq(x, df=df)       h <- hist(res, prob=T, breaks=30, main=paste(df, "degrees of freedom"), ylim=c(0, max(y)), xlab="Pearson statistic", ylab="density")       lines(x, y, col="cornflowerblue", lwd=2)     }          par(mfrow=c(2, 1))     dplot(ps, df=1000-2)     dplot(ps2, df=20-2)      ![enter image description here](http://i.stack.imgur.com/hZyOM.png)


title: Hypothesis testing with a response scale of "1, 2, 3, 4, or not applicable"

text: Let's imagine I'm testing a new drug that purports to reduce the effect of a hangover when the drug is administered shortly before an episode of heavy drinking. On Friday night, I stake out the entrance to a frat party and, using random assignment, give each new guest either a placebo or the active drug. The next morning, I ask each participant to complete this two-item questionnaire: >   1. Did you get drunk last night? > >     * no >     * yes >   2. If "yes," how hungover do you feel this morning? > >     * 1: not at all >     * 2 >     * 3 >     * 4: extremely > Normally I'd perform a simple t-test for the difference in mean response to the Likert scale. But in this case, **the drug could have some effect in discouraging people from getting drunk in the first place**. Anyone who answers "no" to Question 1 necessarily has missing data for Question 2. How can I do the usual hypothesis test -- that the drug does/does not differ from placebo -- given these idiosyncrasies?


title: Help translating the following PCA results into layman's terms?

text: As a non-maths doctoral candidate, I would be absolutely thrilled to bits if anyone here might be willing to take a stab at translating the following PCA results into layman's terms. There are only two sets of raw data: TOEIC test scores and DMIS stages. Additionally, I also need to be able to state the significance, if any, of this PCA analysis (e.g. I have to answer the perennial dissertation question "so what?") Lastly, I must assert that it was definitely not me who performed this PCA, but rather a friend who is conversant with SPSS. Thank you in advance, even if you're not actually in a position to assist me:               Communalities         Initial Extraction     TOEIC   1.000   .546     DMIS    1.000   .546     Extraction Method: Principal Component Analysis.               Total Variance Explained     Component       Initial Eigenvalues                    Extraction Sums of Squared Loadings                Total   % of Variance   Cumulative %   Total   % of Variance   Cumulative %     1          1.092      54.579         54.579       1.092      54.579         54.579     2           .908      45.421        100.000              Extraction Method: Principal Component Analysis.


title: Convergence results for block-gibbs sampling?

text: Suppose you have some complex model you want to sample from by Markov chain Monte Carlo. There are many types of situations where you can divide your variables into, say, two groups, and efficiently and exactly sample the variables in one group conditional on the others. That is, we can sample from $ p(x | y)$ and $ p(y | x)$ and so we would run our Gibbs sampler by repeatedly drawing $x$ from $p(x|y)$ and $y$ from $p(y | x)$. Examples include the restricted Boltzmann machine, scale mixtures of Gaussians in certain image models, and some Bayesian models (one block for the state, and one for parameters). Anyway, _in practice_ , being able to implement block Gibbs sampling seems to lead to much faster mixing of the Markov chain. But, are there any theoretical results on this? Presumably some conditions will need to be satisfied by the two distributions, which will lead to a bound on conductance, or a coupling of some sort. Thanks in advance!


title: Predicting the time until an expected event occurs

text: I have a situation where an event is supposed to occur every x minutes for a number of different sites (each site could be configured for a different time interval x). From time to time the event may fail. I'd like to be able to predict an interval of time that I should worry if an event hasn't occurred. (e.g. for site A 1 missed event indicates a problem that needs attention, but site B misses a lot, so I wouldn't worry about 5 missed events). To make this a little more concrete. Let's say that I have 3 stores that are supposed to upload their sales data every 30 minutes. If store 1 misses an upload, it may be a temporary problem and they'll be able to upload at the next 30 minute interval. This is especially true if this happens frequently for store 1. On the other hand store 2 may almost always upload at the 30 minute interval mark, so that one failure indicates a problem that someone needs to look at. Assuming I have a lot of historical data on how frequently each store fails uploading, how would I go about creating a model for each store based on their history that I can use to determine an estimated prediction of when each store should next upload their data?


title: What graphical model to capture conditional relationships between tokens?

text: I have a list of tokens/words which co-occur with each other. Some of these tokens co-occur more than the others. Moreover, there's a conditional relationship between many. For example, "bat" occurs conditionally with "sport" & "baseball". Such conditional relationships can be on multiple variables. I want to represent these "conditional relationships" (depending on presence/absence of values) using a DAG. Is a graphical model the correct way to do so? Mathematically, I have a large symmetric matrix r x c. The values of this matrix are binary signifying the collocation of r,c viz. if r,c = 1 it means that r,c are co-occurring and may have a dependence between them.


title: Why would I use ANOVA instead of a Rank-Sum test?

text: A colleague of mine with little statistics experience is trying to perform an experimental evaluation of a computer program. He created a between subjects design and solicited test subjects. 11 people were given his "new and improved" computer program to use. 10 others got the "old and boring" computer program to use, for the same task. He asked me, and several other people around the lab how to analyze his data. I told him he should examine the data for normality. If it was normally distributed, he should use a t-test. If it was not, he should use a Wilcoxon rank sum test. One of my colleagues told him he should use ANOVA, even though he only has 2 groups. Apparently using ANOVA on non-normal data in R produces a new degree of freedom measure which can be put into a t-test. I've never heard of such a thing. Is this true? Is it statistically valid? Why would anyone use it instead of just doing a rank-sum test?


title: Why is gender typically coded 0/1 rather than 1/2, for example?

text: I understand the logic of coding for data analysis. My question below is on the use of a specific code.   * Is there a reason why gender is often coded as 0 for female and 1 for male?    * Why is this coding considered 'standard'?    * Compare this with Female = 1 and Male = 2. Is there a problem with this coding?


title: Does it make sense to regress the residuals with rep(0, N) vector?

text: I need to use **Breusch-Pagan test** , to check whether two series are homoscedastic. First I calculate the residuals of two lists of prices after [...] and then I test residuals for homoscedasticity. My doubt is, does it make sense to regress the residuals with a rep(0, N) line? I thought of it because the residuals should be 0 (to have a perfect fit) so i think it makes sense to regress the residuals with a 0,0,0,0,0 etc etc line, what do you think? The only problem that if I use `bptest(residuals ~ rep(0,500))` I get a `NA`


title: Estimator for a binomial distribution

text: How do we define an estimator for data coming from a binomial distribution? For bernoulli I can think of an estimator estimating a parameter p, but for binomial I can't see what parameters to estimate when we have n characterizing the distribution? Update: By an estimator I mean a function of the observed data. An estimator is used to estimate the parameters of the distribution generating the data.


title: Book review of "A Practitioner's Guide to Resampling for Data Analysis, Data Mining and Modeling"

text: I am looking at the text, _A Practitioner's Guide to Resampling for Data Analysis, Data Mining, and Modeling_ , by Phillip Good (Amazon).   * I would like to see if anyone has read this book (or seen a review) and can offer objective advice on it (rating and level of material). (I have a suspicion this question - being subjective - should be placed in a community location. Problem is, not sure how to do that. So, I will beg forgiveness in advance and look to a moderator to help direct me.)


title: How can I check if two numbers are equal (with some allowed error)?

text: I have a lot (about 10000 or more) measurement results. I measured the performance of different algorithms (it doesnt really matter which algorithms for now on). I want to check if my measurement results are stationary or not. I can simply check it with SIMPLE MOVING AVERAGE, and I did it. The thing is, I have to compare my moving averages I calculated, check if they're equal or not. If they're equal, my data set is stationary. But having so much data I need to consider a small error, with I can assume my means are equal. Example:   1. I have a data set with the measurement results: `X = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0};`, so `N = 10`   2. I chose the window, `w = 3`   3. I calculated means with the simple moving average method, stored means in `SAM` table: (1+2+3)/3 = 2, (2+3+4)/3 = 3, (3+4+5)/3 = 4, ..., (8+9+10)/3 = 9 so `SAM = {2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0}`   1. For my data to be stationary, `SAM` table should look like this: `SAM = {9.0, 9.0, 9.0, 9.0, 9.0, 9.0, 9.0, 9.0}` (moving means are equal)   2. Having `X` set with some randomness, my SAM will not look quite that neat, it could be `SAM = {9.2, 9.7, 9.3, 9.6, 9.4, 10.0, 9.1, 9.7}` and heres my question: how can I calculate the acceptable error with I can assume my means are equal? When `SAM` looks like this: `SAM = {9.2, 9.7, 9.3, 9.6, 9.4, 10.0, 9.1, 9.7}` and the error is, lets say +/-0.7 then I can consider my means as equal and my data stationary. The main problem is, _how can I calculate this acceptable error?_


title: How to prove the significance of features in classification?

text: I have a binary classification problem. I have extracted 500 features from a set of 5000 samples using my domain knowledge. In other words, I have got **hand crafted features**. > I wish to prove that these features actually are enough for performing > classification and they make the 2 classes of samples separable. i.e. When > the samples are represented with these features, there is exists a > (reasonable) decision boundary. Please advice on how I can prove this. Is there any statistically appropriate way of measuring the significance of the set of features as a whole (NOT the significance of individual features)? EDIT: Suppose if I observe that different families of classifiers such as kNN, DT, NB, ANN, SVM etc. (with 10 fold cross validation) offer reasonable accuracy, can I conclude that the features are meaningful?


title: What is the best analysis to perform on brute force tests?

text: Basically, I have run a network simulation for video streaming a number of times. I have assigned weights to particular performance metrics of the network measuring their performance each time. Depending on the weight configuration, the outputted quality of the video (measured in a mean opinion score or MOS) is different each time. What I want to know is which metric was most important for a good MOS and which metric was the least important for a good MOS score. I am not very good at statistics so if I wanted to show this which test should I perform in Minitab? The test was run 16 times, with the weights selected at random. EDIT: To go into more detail, my results showed that high values for at least two of three particular performance metrics resulted in one of the highest results. How do I represent this?


title: 2 intercept cooficients in glmnet output

text: I am not sure I understand why glmnet returns 2 intercepts in its result. If model matrix already has intercept column, why do we need another intercept for the model?               > mm = model.matrix(...)     > colnames(mm)      [1] "(Intercept)"  "v1" "v2"....          > model.fit = glmnet(mm, resp, standardize=FALSE, family='binomial', nlambda = 10)     > coef(model.fit, s= 1e-05)       48 x 1 sparse Matrix of class "dgCMatrix"                        1       (Intercept)        5.9487109234       (Intercept)        .                  v1        0.2580378293       v2       -0.4693849121       .....


title: Can I use simple regression with variables that are not significant in multiple regression?

text: As I am testing a number of models, I notice that none of my demographic variables are ever signficant. For example, I am testing a model to predict the dependent variable 'perceived substitutability'. The independent variables are contentment about the content on traditional television and contentment about the content on online television. I run this model together with gender, age, professional status, income (yes/no), family composition and marital status. I used a hierarchical method where I used the demographic variables in the first block and the other variables in the second block. As I run the regression, all demographic variables were not significant (p>0,05). This is of course due to correlations between the independent variables. If I use one demographic variable at a time to explain the dependent variable (so six times a simple regression where no other demographic variables are used), professional status and income become significant in the simple regression (not together because of course they correlate). Can I use this in my result-section: that in the simple regressions (1 dependent, 1 independent) they are significant while in the multiple they were not? I am using SPSS


title: Redundancy in Performance Measures

text: I have a performance rating that correlates positively with a variety of other parameters used individually to evaluate the performance rating of a given operation. However, I suspect that many of these other parameters are incorporated in the performance rating and that these parameters 'compose' the performance rating, whose composition is not known. To rephrase, I know nothing about the functional dependence of the performance rating, but do know how these other parameters correlate with both the performance rating and the performance itself (the rating is noisy so these are not necessarily the same). To define these terms more clearly let me give an example: Performance Rating: 4.5 Performance: Estimated RPM = 1000 In other words, the rating creates a unit-less way to compare how something is expected to perform. It turns out that the performance rating predicts performance in a way that is similar to the way specific parameters predict performances, suggesting that these parameters might be components of the rating. Is there a best way to go about testing whether the reason the performance rating gives similar results is because it is composed of some of these parameters? That is, is there a way to test whether the predictive capacity of the performance rating is entirely due to these parameters? Would this just follow a basic regression?


title: confidence interval for co-variance

text: is there a way to compute the confidence interval of a co-variance matrix in Matlab.. Suppose I have 2 vectors of stock returns. If I compute a co-variance matrix (2x2) the sample co-variance will not differ from zero. However, I need to know if the estimate co-variance is different from zero.. In order to do so, I will need the confidence interval of the co-variance. I searched a lot on google and other forums. First result was a Wald test, but I not sure if this one is correct. The second one is to normalize the co-variance with the std of both assets to obtain the correlation coefficient.. The second one is not appropriate for my research.. Has someone a suggestion? Thanks in advance for your reply.. Could this be the answer? cov(x,y)= s_hat the confidence interval of s_hat is ((s_hat/(Chi^2_n(1-\alpha/2)) , (s_hat/Chi^2_n(alpha/2))) where Chi^2_n is the chi-square distribution with n degrees of freedom


title: Equations in the news: Translating a multi-level model to a general audience

text: The New York Times has a long comment on the 'value-added' teacher evaluation system being used to give feedback to New York City educators. The lede is the equation used to calculate the scores - presented without context. The rhetorical strategy appears to be intimidation via math: ![alt text](http://graphics8.nytimes.com/images/2011/03/07/education/07winerip_graphic/07winerip_graphic- articleLarge-v2.jpg) The full text of the article is available at: http://www.nytimes.com/2011/03/07/education/07winerip.html The author, Michael Winerip, argues that the meaning of the equation is beyond the capacity of anyone other than, um, Matt Damon to understand, much less an average teacher: > "The calculation for Ms. Isaacson’s 3.69 predicted score is even more > daunting. It is based on 32 variables — including whether a student was > “retained in grade before pretest year” and whether a student is “new to > city in pretest or post-test year.” > > Those 32 variables are plugged into a statistical model that looks like one > of those equations that in “Good Will Hunting” only Matt Damon was capable > of solving. > > The process appears transparent, but it is clear as mud, even for smart lay > people like teachers, principals and — I hesitate to say this — journalists. > > Ms. Isaacson may have two Ivy League degrees, but she is lost. “I find this > impossible to understand,” she said. > > In plain English, Ms. Isaacson’s best guess about what the department is > trying to tell her is: Even though 65 of her 66 students scored proficient > on the state test, more of her 3s should have been 4s. > > But that is only a guess." How would you explain the model to a layperson? FYI, the full technical report is at: http://schools.nyc.gov/NR/rdonlyres/A62750A4-B5F5-43C7-B9A3-F2B55CDF8949/87046/TDINYCTechnicalReportFinal072010.pdf Update: Andrew Gelman offers his thoughts here: http://www.stat.columbia.edu/~cook/movabletype/archives/2011/03/its_no_fun_bein.html


title: Estimate how coefficient of variation changes as a function of $x$ without replication

text: Suppose I have data points of the form (x,y) and I know that CV should be constant across all x values. (i.e. heteroskedastic data). That being said, how can I compute the estimated CV if I don't have multiple y values for each x? Typically, CV is computed on discrete data (ex. repeatedly measuring y values for some given x). But in this case, not all x values have multiple y values. Note that constant variance can be assumed throughout the domain of x. If anyone could direct me to literature that handles scenarios as this, that would be great.


title: Regression: finding p-value and F

text: I'm working on stats homework, and I've come across a problem I cannot figure out. It's a multi-part problem, and I've correctly answered a number of the questions, but there are a few that are not coming up right. Problem: > Healthy bones are continually being renewed by two processes. Through bone > formation, new bone is built; through bone resorption, old bone is removed. > If one or both of these processes is disturbed, by disease, aging, or space > travel, for example, bone loss can be the result. Tartrate resistant acid > phosphatase (TRAP) is a biochemical marker for bone resorption that is also > measured in blood. The data below gives values for this biomarker and a > measure of bone resorption `VO--`. Analyze these data using the questions in > the previous exercise as a guide. (Give answers accurate to 2 decimal > places.) Data > > Analyze these data using the logs of both TRAP and `VO--`. Run the > regression using TRAP to predict `VO--`. Summarize the results. (Round your > answers to two decimal places.)   >  $\log($VO--$) = \mathbf{5.62} + .42 \cdot \log($TRAP$)$   >  $F = \mathbf{1590.97}$   >  $P = \mathbf0$   >  $R^2= .1$ I wasn't quite sure how to figure out the answers, but I entered the data into list 3 and list 4 on my TI-83 calculator. Then I entered `ln(L3)` into list 1 and `ln(L4)` into list 2. From there, I did `STAT > TESTS > LinReg(a+bx)` to get the regression equation: $y = 5.62 + .42x$ and $R^2= .1$. I input these numbers, and the $R^2$ and slope were correct, but the $y$-intercept was not. I then used `STAT > TESTS > ANAOVA(L1, L2)` to get $F = 1590.97$ and $P = 0$, but neither of those answers were correct. I've listed my answers above with incorrect answers bold. Can anyone help let me know what I'm doing wrong? I'd be very appreciative!


title: SVM confusion matrix whose dimensions are more than two

text: While I am using SVM, I train it with a train data and then I try to predict a sample if its label is -1 or +1. However, I see some confusion matrice for SVM like below. Mine are 2x2 matrice but their dimesions are larger, e.g. 15x15. Do they have more than one SVM? How people do such things? It is easy to write one diagonal but how they decide other values? In the figure below, for 2, it is easy to write 38 but how they write 2? BTW, I am new to these topics. ![15x15 confusion matrix](http://i.stack.imgur.com/zsZn9.png)


title: How to NOT use statistics

text: This is sort of an open ended question but I wanna be clear. Given a sufficient population you might be able to learn something (this is the open part) but whatever you learn about your population, when is it ever applicable to a member of the population? From what I understand of statistics it's never applicable to a single member of a population, however, all to often I find myself in a discussion where the other person goes "I read that 10% of the world population has this disease" and continue to conclude that every tenth person in the room has this disease. I understand that ten people in this room is not a big enough sample for the statistic to be relevant but apparently a lot don't. Then there's this thing about _large enough_ samples. You only need to probe a large enough population to get reliable statistics. This though, isn't it proportional to the complexity of the statistic? If I'm measuring something that's very rare, doesn't that mean I need a much bigger sample to be able to determine the relevance for such a statistic? The thing is, I truly question the validity of any newspaper or article when statistics is involved, they way it's used to build confidence. That's a bit of background. Back to the question, **in what ways can you NOT or may you NOT use statistics to form an argument**. I negated the question because I'd like to find out more about common misconceptions regarding statistics.


title: How to quickly sample X if exp(X) ~ Gamma?

text: I have a simple sampling problem, where my inner loop looks like:               v = sample_gamma(k, a)      where `sample_gamma` samples from the Gamma distribution to form a Dirichlet sample. It works well, but for some values of k/a, some of the downstream computation underflows. I adapted it to use log space variables:               v = log(sample_gamma(k, a))      After adapting all of the rest of the programme, it works correctly (at least it gives me the same exact results on test cases). However, it is slower than before. Is there a way to directly sample $X, \exp(X) \sim \text{Gamma}$ without using slow functions like $\log()$? I tried googling for this, but I do not even know if this distribution has a common name (log-gamma?).


title: Issues with sequential feature selection

text: I am trying to do some feature selection in gene expression data with 22215 features. I followed the tutorial here. I initially applied filter method(ttest) to select the features having the best p values. I selected 100 features from them initially. Then I tried to apply sequential feature selection method on them with SVM classifier. However, when I do               [fs1, history] = sequentialfs(@SVM_class_fun, reducedL, yS1, 'cv', c);      it always returned me the 1st feature only. I mean in fs1 every other feature except the first one is 0. If I try to force it to give me 10 features with               [fs1, history] = sequentialfs(@SVM_class_fun, reducedL, yS1, 'cv', c, 'nfeatures', 10);      Here is my SVM_class_fun               function err = SVM_class_fun(xTrain, yTrain, xTest, yTest)       model = svmtrain(xTrain, yTrain, 'Kernel_Function', 'rbf', 'boxconstraint', 10);       err = sum(svmclassify(model, xTest) ~= yTest);     end      it will give me the first 10 selected by the filter method having lowest p values. So this mean using sequentialfs is not helpful in this case. To let you know I have just 12 examples. So my data matrix is of dimension 12x22215. Might this be the issue? Can anyone provide some insights?


I know correlation does not imply causation but instead the strength and direction of the relationship. Does simple linear regression imply causation? Or is an inferential (t-test, etc.) statistical test required for that?


title: AIC or p-value: which one to choose for model selection?

text: I'm brand new to this R thing but am unsure which model to select.   1. I did a _stepwise forward regression_ selecting each variable based on the lowest AIC. I came up with 3 models that I'm unsure which is the "best".              Model 1: Var1 (p=0.03) AIC=14.978     Model 2: Var1 (p=0.09) + Var2 (p=0.199) AIC = 12.543     Model 3: Var1 (p=0.04) + Var2 (p=0.04) + Var3 (p=0.06) AIC= -17.09      I'm inclined to go with Model #3 because it has the lowest AIC (I heard negative is ok) and the p-values are still rather low. I've ran 8 variables as predictors of Hatchling Mass and found that these three variables are the best predictors.   2. My next forward stepwise I choose Model 2 because even though the AIC was slightly larger the p values were all smaller. Do you agree this is the best?              Model 1: Var1 (p=0.321) + Var2 (p=0.162) + Var3 (p=0.163) + Var4 (p=0.222)  AIC = 25.63     Model 2: Var1 (p=0.131) + Var2 (p=0.009) + Var3 (p=0.0056)                  AIC = 26.518     Model 3: Var1 (p=0.258) + Var2 (p=0.0254)                                   AIC = 36.905      thanks!


title: proper statistical test to check differences

text: I have 3 replicas for a value in different individuals. Each of these values are ratios $ab$, and $a$ and $b$ are means from $n=20$ sample pool each. Thus there are 3 times each ratio $ab$ for each individual. When comparing for differences is student t-test correct?                                1st sampling -   2nd         -   3rd         -   4th     individual 1  -  0.7164165213 -  0.6057539083 -  0.5242174359 -  0.7670756899     individual 2  -  0.6540839702 -  0.8140762612 -  0.6057645321 -  X     individual 3  -  0.611493629  -  0.7270260938 -  0.5255522645 -  0.9964242368      as said each sampling value comes from A/B, whereas A and B are means for treated and non treated for a given variable with n ranging from 15-25.


title: How to interpret 95%CI's in conjunction with significance tests for differences between group means?

text: I have data for two groups, and I am investigating if the two groups perform differently on a task. I have 95% CIs calculated by bootstrapping that give lower and upper bounds to mean point estimates for each group. I have done a t-test comparing the means and it suggests that the difference between the means is significant, suggesting an effect. I am having trouble interpreting what the CIs mean, in conjunction with the sig. test. For example, the CIs do overlap to a small degree. Therefore, I interpret that the groups may have similar means that the population level and a significant different may not be observed. Is this correct? I am also wondering the best way to report significance tests and CIs as complimentary in a report. Is it alright to discuss a significant difference as important, even though the CIs indicate a degree of overlap?


title: When can I suppress the intercept using treatreg?

text: Can I suppress the intercept if I know the treatment will be zero if the independent variables are zero. Also, can I suppress the intercept if I know the right hand side of the primary regression equation CANNOT equal zero. Example:               treatreg y1 x1 x2 x3, treat(y*=x4 x5, noconstant) noconstant      Can I do this if I know for sure `x1`, `x2` and `x3` cannot combine to equal zero? Also if `x4` and `x5` combine to zero then `y*` MUST be zero. Is this an appropriate use of noconstant?