title: Testing periodogram "peaks": sine-like wave or AR/MA/ARMA noise?

text: I'm performing an harmonic fit to data I know (from physical constraints) come from a periodic source of the form $$\sum_j^M \sum_i^N a_{i,j}\sin(2\pi f_it)+b_{i,j}\cos(2\pi f_it)$$ using the Lomb-Scargle periodogram (I know this is not the optimal thing to do, but it's fast and reliable so far), where I'm also trying to estimate $N$ and $M$, the number of frequencies and the number of harmonics per frequency to fit. The number of harmonics is really not a very hard problem to solve once the number of frequencies is fixed, because I can perform methods like Ridge Regression, LASSO, etc. in order to find the best subset of harmonics, so my problem is really on estimating the optimal number of frequencies. The thing is that I know red noise is almost always present, and certain harmonic fits give AR-like spectra on the residuals (and, at least to me, it seems like I could even have ARMA-like spectra). The "classical" Lomb-Scargle periodogram, on the other hand, assumes white noise in order to test the significance of the peaks, so the significance test on the peaks doesn't seem like a good idea in order to test if there are any residual frequencies on the periodogram. I've been trying to derive significance tests assuming ARMA-like noise, but the problem of defining the order $p,q$ of the process arises: do you know of a way of differencing between correlated noise and sine wave-like spectra? This paper shows a way of doing this assuming AR(1) noise (in fact something more like a CAR-type noise?), but the problem of it is that it actually knows there is red noise. I think I can "see" how to do it: sine wave-like spectra has almost-symmetric decreasing peaks around the main peak because of the (uneven) sampling, which ARMA-like noise doesn't have (the amplitudes are random with an apparent upper limit dictated by the form of the spectra). I've also been thinking in performing fits with ARMA and sine-like waves, and perform model selection via AIC...any suggestions?


title: GWAS and Statistical theory - does the likelihood of a detectable main effect decrease with complexity?

text: I've been wondering recently about the difficulty of detecting statistically significant results from, say, genome-wide association studies. In these studies, many - ($10^{8}$) for example- statistical tests on tiny genetic changes are performed in order to detect an association with a trait - say heart disease. The difficulties are manifold. First, multiple testing. Simple Bonferroni correction for $10^{8}$ tests often leaves these entire experiments labelled "insignificant" if they don't find an association that withstands these mad Bonferroni corrections. (There's probably a case for the "Manhattan" plot being replaced with a Pareto plot, because of the hypothesis generating nature of these studies. It's a cultural issue and a little off my topic). However, using narrower sets of data and better defined phenotypes we can sometimes detect statistical interactions between particular genotypes which are associated with a particular (endo)phenotype. My questions is: is it in the nature of complex systems that main effects will "cancel out" but that statistical interactions will be more detectable? With many more "true" effects in your system - is it less likely you will detect those effects as main effects without controlling for others? This might be particularly true of the genome, but could be true of other complex systems. I don't really have a good structure in which to think about this problem and I wonder if it has been written about before. Thanks


title: What is the independent variable?

text: This question will more than likely be extremely simple for most of you, but the way the scenario is worded has me confused. Scenario: The communication quality of couples three month before and again three month after marriage were tested. It was assumed that the communication quality would decrease after marriage. One group studied was 19 couples who had received ordinary premarital counselling from the ministers who were going to marry them. The scores for the 19 husbands are listed as before and after. Question: What is the independent variable of the above? I understand the difference between IV and DV, however I'm unsure if the IV is marriage or the fact that 19 couples received premarital counselling due to the way the question is worded. All help will be appreciated!


title: When to use accuracy and precision to evaluate binary classifiers?

text: I came a cross two ways to evaluate the performance of binary classifiers: accuracy and precision. When to choose each? And what are the advantages and disadvantages of each one?


title: What is the probability of picking p primes from n random numbers?

text: The probability of that a number between 1 to x is prime is $\frac{1}{\ln{x}}$ as per prime number theorem and also the total number of prime numbers between $1$ to $x$ will be $\frac{x}{\ln{x}}$. But if we select $n$ (32 bit) random numbers, what is the probability that $p$ of them are primes? Or Simply put What is the probability of picking $p$ primes from $n$ random numbers (32 bit). TIA..


title: Should I re-shuffle my data?

text: We have a set biological samples that was pretty expensive to obtain. We put these samples through a series of tests to generate data that are used for building a predictive model. For this purpose we have divided the samples to training (70%) and testing (30%) sets. We have successfully created a model and applied it on the testing set to discover that the performance was "less than optimal". The experimentalists now want to improve the biological tests in order to create a better model. Provided that we cannot obtain new samples, would you suggest us to re-shuffle the samples to create new training and validation sets or to stick with the original division. (We don't have any indication that the division was a problematic one).


title: MCMC method -- good sources

text: Any suggestion for a good source to learn MCMC methods? thanks.


title: Exploratory factor analysis with fixed factors

text: Is it possible to conduct exploratory factor analysis using 'fixed number of factors' instead of going for 'factors having Eigenvalues more than 1'? Are there any technical problems in that?


title: Is it correct to combine PCA and NMDS axes in a multiple regression?

text: I am considering to do a multiple regression in which some of the predictive variables are PCA (principal components) axes whereas others are NMDS (nonmetric muliple dimension scaling) axes. I would like to know **if this is incorrect and why**. I'm considering this because I have two groups of many variables that I want to reduce before analyzing the relative effects of each of these groups over another variable. The first group is composed by anatomical variables from several species which are best represented by a Phylogenetic PCA. The second group is a bunch of environmental variables which are very non-normal and that I can't normalize through transformation (I tried several and the variables remained with exponential distribution). What would be a correct way to reduce all variables before doing the multiple regression via generalized least squares?


title: Logistic Regression with holdout sample

text: I am trying to run a logistic regression in R on my data where my independent variables are 13 continuous variables and my dependent variable is binary. I want to segment my data so that I train on the first 80% and test on the last 20%. I have a total of 3750 rows of data so I utilize the first 3000 for training. I have written the following:               mydata<-totaldata[1:3000,2:15]     mylogit<-glm(mydata$TARGET ~ mydata$VAR1+mydata$VAR2+mydata$VAR3+mydata$VAR4+ #$                                  mydata$VAR5+mydata$VAR6+mydata$VAR7+mydata$VAR8+                                  mydata$VAR9+mydata$VAR10+mydata$VAR11+mydata$VAR12+                                  mydata$VAR13, family="binomial")          predictdata=totaldata[3001:3751,3:15]     in_frame<-data.frame(predictdata)     predictions=predict(mylogit,in_frame,type="response")      However I get the following warning message: Warning message: 'newdata' had 751 rows but variable(s) found have 3000 rows Then when I look at predictions there are 3000 predictions not the 751 that I wanted. What can I do to fix this?


title: Solving a difficult equation for a variable?

text: I'm trying to obtain the maximum likelihood estimate of the parameters for a model I'm building. I have constants $\sigma$, $\mu$, and $q_0$; a boolean matrix $\alpha$; and vectors $A, \beta, r, d,$ and $Q$, all quantities non- negative. I set the derivative of the log likelihood function equal to zero with $$0=\sum_j \alpha_{kj} \Bigg[ \beta_j r_j \frac{\ln (A_j - r_j\sum_i \alpha_{ij} Q_i) + \sigma^2 - \mu}{\sigma^2 (A_j - r_j\sum_i \alpha_{ij} Q_i)} + (1-\beta_j)\frac{\ln(1-d_j)}{q_0} \Bigg]$$ How can I solve for $Q$?


Is there a Box-Cox like transformation for independent variables? That is, a transformation that optimizes the $x$ variable so that the `y~f(x)` will make a more reasonable fit for a linear model? If so, is there a function to perform this with `R`?


title: Comparing estimators in Cauchy distribution

text: Given $n$ observations with a Cauchy distribution (location=t,scale=1), I would like to compare estimators such as the mean and median by simulating sample means and mean square errors. This is an exercise given by my book. The exercise says I may assume t=0. Why is this? I thought I would try $n=10$ first so I typed in `x=rcauchy(10,location=0,scale=1)` and then median(x). Is this the right way to go about it? Also, what is the code for mean square error and coverage probabilities? Note: Advised to come here from Stack Overflow.


title: Getting started on BUGS, need help implementing a simple model

text: If I have a parameter with a prior $X\sim N(\mu,\sigma)$ and three data points, $x_1,x_2,x_3$, how can I calculate the posterior in R or BUGS? Here is an attempt at doing this in BUGS, that I have gotten stuck on:               data:     list(x=1.4, 2.1, 1.1, mu = 10, tau=1/10)          model {       for (i in 1:3) {            x[i] ~ dnorm(mu, tau)       }     }      Any help appreciated


title: Generation of distribution with given mean, standard deviation, lower and upper bounds

text: How to generate distribution in R, members of which are inside bounds and mean with SD are fixed? I just need a distribution which will fit in the supported creteria. This will be used for simulating results from other study for comparison with data I have. I need something like this:               rBootstrap<-function(n,mean,sd,lowerBound,upperBound){        ...     }      which result of `d <- rBootstrap(100,50,10,30,100))` had `min(d)` $\approx$ 30, `max(d)` $\approx$ 100, `mean(d)` $\approx$ 50, `sd(d)` $\approx$ 10. The best result I had achieved is:               rBootstrap<-function(n,mean,sd,lowerBound,upperBound){         data <- rnorm(n,mean,sd)         return(data)     }     summary(rBootstrap(100,50,10,30,100))        Min. 1st Qu.  Median    Mean 3rd Qu.    Max.        26.18   41.27   47.29   48.43   55.06   71.82      and it is very far from desired. I had tried several others but without success. **Is it possible to do this at all and in R in particular?**


title: Is this a two way ANOVA with replicates?

text: I have the following results from a biological experiment. I measure the size of the sample (2,4,8,16,32 or 64) and count its distribution in a population of 50 (sum of counts=50). I do those measurements in the wildtype (wt) and the experimental sample (KO) in three independent experiments (replicate 1,2,3). Am I in a setup for a two way ANOVA with replicates? How can I proceed rigorously, with the maximum power in R as I am not sure (1) my data is normally distributed and (2) how to code that in R anyway? Thanks!               strain  Replicate   size    count     wt  1   2   2     wt  1   4   6     wt  1   8   11     wt  1   16  12     wt  1   32  9     wt  1   64  10     wt  2   2   2     wt  2   4   7     wt  2   8   15     wt  2   16  11     wt  2   32  4     wt  2   64  11     wt  3   2   0     wt  3   4   5     wt  3   8   17     wt  3   16  9     wt  3   32  9     wt  3   64  10     KO  1   2   0     KO  1   4   0     KO  1   8   2     KO  1   16  15     KO  1   32  23     KO  1   64  10     KO  2   2   0     KO  2   4   2     KO  2   8   2     KO  2   16  20     KO  2   32  22     KO  2   64  4     KO  3   2   0     KO  3   4   0     KO  3   8   3     KO  3   16  20     KO  3   32  21     KO  3   64  6


title: How can one plot continuous by continuous interactions in ggplot2?

text: Let's say I have data:               x1 <- rnorm(100,2,10)     x2 <- rnorm(100,2,10)     y <- x1+x2+x1*x2+rnorm(100,1,2)     dat <- data.frame(y=y,x1=x1,x2=x2)     res <- lm(y~x1*x2,data=dat)     summary(res)      I want to plot the continuous by continuous interaction such that x1 is on the X axis and x2 is represented by 3 lines, one which represents x2 at a Z-score of 0, one at Z-score of +1, and another at a Z-score of -1, with each line a separate color and labelled. How can I do this using ggplot2? For example, it might look something like this (though of course with different colored lines rather than different line types): ![Example Image](http://i.stack.imgur.com/eUZJG.png)


title: Reading random text file in R

text: I have about 500 text documents stored in a directory. I want to randomly select one to read in the contents (this process will be repeated several times). The names of the files are not sequential. By this I mean they are not named something like 1.txt, 2.txt, 3.txt,... The names are basically random characters. How can I do this? Thanks


title: Test data has maximum number of unseen nGrams and prediction is failing

text: I have built a Machine learning model for twitter sentiment analysis. I have followed following steps to build my final model. Model Creation:   1. collected uni-bi-tri grams from entire training data after cleaning twitters removing URLs,stop words etc.   2. selected top nGrams based on it total frequency in the training data   3. Now used those selected nGrams to create features for each training tweets   4. Feature scaling   5. verified feature space doing PCA plot whether my features are good enough for machine to segregate positive and negative training tweets. ![PCA plot to describe feature space of training data](http://i.stack.imgur.com/atPZO.png)   6. Apply SVM (RBF) model to create final model for sentiment classification Testing/Prediction:   1. identify which nGrams which is part of nGram dictionary prepared during training are present in new tweet    2. created features using seen nGrams in the tweet   3. scale feature using calculated mean and standard deviation from training feature set.   4. finally called SVM prediction module for sentiment classification Result: training data prediction is very high. almost 95% test data prediction completely depend on seen nGrams. if number of seen nGrams is more then prediction is correct otherwise wrong. Now my question is that how can I increase seen nGram count in new tweet. increasing training data is only possible solution or is there any other approach? even if I increase training data then also no one can give guaranty that new tweets will have some number of seen nGrams. Is there any smart approach to solve this problem?


title: How is an ANOVA calculated for a repeated measures design?

text: The title says it all, and I'm confused. The following runs a repeated measures aov() in R, and runs what I thought was an equivalent lm() call, but they return different error residuals (although the sums of squares are the same). Clearly the residuals and fitted values from aov() are the ones used in the model, because their sums of squares add up to each of the model/residual sums of squares reported in summary(my.aov). So what are the actual linear models that are applied to a repeated measures design?               set.seed(1)     # make data frame,     # 5 participants, with 2 experimental factors, each with 2 levels     # factor1 is A, B     # factor2 is 1, 2     DF <- data.frame(participant=factor(1:5), A.1=rnorm(5, 50, 20), A.2=rnorm(5, 100, 20), B.1=rnorm(5, 20, 20), B.2=rnorm(5, 50, 20))          # get our experimental conditions     conditions <- names(DF)[ names(DF) != "participant" ]          # reshape it for aov     DFlong <- reshape(DF, direction="long", varying=conditions, v.names="value", idvar="participant", times=conditions, timevar="group")          # make the conditions separate variables called factor1 and factor2     DFlong$factor1 <- factor( rep(c("A", "B"), each=10) )     DFlong$factor2 <- factor( rep(c(1, 2), each=5) )          # call aov     my.aov <- aov(value ~ factor1*factor2 + Error(participant / (factor1*factor2)), DFlong)          # similar for an lm() call     fit <- lm(value ~ factor1*factor2 + participant, DFlong)          # what's aov telling us?     summary(my.aov)          # check SS residuals     sum(residuals(fit)^2)       # == 5945.668          # check they add up to the residuals from summary(my.aov)     2406.1 + 1744.1 + 1795.46   # == 5945.66          # all good so far, but how are the residuals in the aov calculated?     my.aov$"participant:factor1"$residuals          #clearly these are the ones used in the ANOVA:     sum(my.aov$"participant:factor1"$residuals ^ 2)          # this corresponds to the factor1 residuals here:     summary(my.aov)               # but they are different to the residuals reported from lm()     residuals(fit)     my.aov$"participant"$residuals     my.aov$"participant:factor1"$residuals     my.aov$"participant:factor1:factor2"$residuals


title: using outcome of clustering as an ordinal scale for regression - feasible?

text: Suppose I have access to five dental surgeries, who volunteered to collect data about patients and their regular check-ups. The dental surgeries are quite up to completely different from each other (size, sociodemographics of patients etc.). The collection of data for each patient runs as a small questionnaire (14 items on medical and sociodemographic preconditions). I want to identify, what has an influence on how long (= how many minutes) the patient needs for the check-up and potential ad-hoc therapy in need. There is _no_ reason to assume, that this duration in minutes is normally distributed, not even for the single dental surgery. Does the following idea hold - or is there a more direct way to work through? 1) I apply partitional clustering for each of the dental surgeries, but only on the duration of check-up plus therapy. I.e. I get one-dimensional values (in minutes), which I join to N groups (N a predefined number, e.g. 5) in such a way, that the mean of each group is as far away as possible from the means of the other groups. 2) Step 1) has established an ordinal scale for each dental surgery. I have heeps of data. Now I can start an ordinal regression... 2a) for each dental surgery separately 2b) for all dental surgeries together (?) My main questions are: is 1)+2) as establishing a scale of one's own feasible? Are there (easier, more correct) alternatives? One alternative I could imagine are quintiles instead of 1). Statos


title: How to analyse the relationship between many categorical variables

text: I have a number of categorical variables, each with varying numbers of levels: gender of client (m,f) age group (0-8, etc.) diagnosis (mild, moderate, severe, profound) location (community, residential) restrictive interventions used (physical restraint, chemical restraint, seclusion, etc.) (and more) I want to see if certain group characteristics, e.g. a client's age or gender (or both combined) are more predictive than others of receiving a particular restrictive intervention. E.g. Does gender predict type of restrictive intervention received more strongly than range of diagnosis? What is the strongest predictor of seclusion? I'm pretty much at a loss in terms of what test to use to answer these questions. I'm at a basic level of statistics, so any help would be appreciated. Cheers!


title: Tukey-Kramer test

text: Can anybody point me out a good reference for the Tukey-Kramer test? Secondly, suppose that we want to compare the mean of the group called "male" and the group called "female" (each group is composed of 10 sample). I executed this test with JMP and it gives to me a comparison table like the following:                      female   male        male   2.18  -3.15     female  -3.11   2.18      What do these numbers represent? From the JMP manual, I know if they are positive there is a significant difference, if they are negative there is no significant difference. If I read a negative value in the cell male-male or female-female, does it means that among the ten sample of male (or among the 10 samples of female) there is no significant difference in term of their mean?


title: How to detect changepoint with a bit of tolerance?

text: I'm using the **changepoint** package in R because I need to check if the variance is stable or if there are some breaks inside. The function **cpt.var** does the job well but there is not a parameter to set the "lags". With lags I mean a parameter to specify how many errors can be accepted before should be consider it a breaks. Take a look at the image below. The function correctly finds the break(below -2 on the plot) but as you can see the variance of the follow residuals continues constant. ![Residuals](http://i.stack.imgur.com/9l8sk.jpg) Is there a similar function with the possibility to set the number of "errors" that should be found before considering it a break? Thank you!


title: Forecasting binary time series

text: I have a binary time series with 1 when the car is not moving, and 0 when the car is moving. I want to make a forecast for a time horizon up to 36 hours ahead and for each hour. My first approach was to use a Naive Bayes using the following inputs: t-24 (daily seasonal), t-48 (weekly seasonal), hour of the day. However, the results are not very good. Which articles or software do you recommend for this problem?


title: High precision with low recall SVM

text: I'm classifying a data set using SVM and those are the precision and recall values for two classes.                    precision    recall  f1-score   support          H       0.91      0.99      0.95      1504     R       0.81      0.23      0.36       192          avg/total 0.90      0.91      0.88      1696      Following is the confusion matrix.               Confusion Matrix:     [[1494   148]      [ 10    44]]      ![enter image description here](http://i.stack.imgur.com/2hTZS.png) Can I say this is a good classifier for my dataset based on average? I'm not sure because I'm getting a low recall value for class 'R'?


title: Compare patterns of text in R

text: I would like to compare a number of patterns within a huge dataset. I have one variable that is in textform, i.e.               1 der tus wif     2 fiu seg suf erf utz tus wif      What I would want to find is how many of these snippets differ from one line to the next one. In this case, "wif" and "tus" overlap to the next line, but the other text snippets are new or different.? How can I create such variables? I have already assigned all snippets into one variable each. Does it help?


title: Testing for paired subjects AND two timepoints

text: How should I best analyze data that have "2x pairing"? E.g., _N_ set of twins, each twin randomly assigned to active drug or placebo, then measure indicator _X_ at timepoints _A (before administration)_ and _B (after administration)_. My understanding here is that there are 2 pairings: between time _A_ and _B_ , and between drug and placebo. Does it make sense to manually calculate the difference in _X_ between the timepoints and just perform a simple paired t-test / Wilcoxon on the differences? If that's OK, then is there a difference if I do it the other way around, i.e. manually calculate the difference in X between drug and placebo, and perform a simple paired test across the timepoints? (It feels less logical, but I can't place why exactly)


title: How do I report my ranked ordinal data?

text: I passed out a survey to the community and I am unsure how to properly report the data. The question is: "Rate from 1-5: one being the most important and five being the least important. I will attend the City Heights Seniors Activity Group because I want to:   * (__) Be more social   * (__) Exercise More   * (__) Be creative   * (__) share new cultures   * (__) Learn new skills" I am wondering how to record the ranked results on Microsoft Excel, and also how to show it in a table. Any advice would be great!


title: Methods for temporal patterns extraction

text: For example a video or series of images, or usage patterns data on a website, or a univariate time series, is there some flexible methods for extracting patterns of any length, such as head movements/gestures for videos, trends and similar repetitive patterns for time series? I've learned some simple clustering methods like K-means and PCA for reducing dimensions but I don't understand how those could work with the basic dataset of temporal data, since we need to group only consecutive samples to form patterns. The only solution I know is windowing $\begin{bmatrix}a_1 b_1\\\ a_2 b_2\\\ a_3 b_3\\\ a_4 b_4\\\ a_5 b_5\\\ a_6 b_6\end{bmatrix}$ becomes with w=3 $\begin{bmatrix} a_1 a_2 a_3 b_1 b_2 b_3\\\ a_2 a_3 a_4 b_2 b_3 b_4\\\ a_3 a_4 a_5 b_3 b_4 b_5\\\ a_4 a_5 b_6 b_4 b_5 b_6 \end{bmatrix}$ but it doesn't look flexible Finding an optimal w is hard, and wouldn't fit for different sizes of patterns.   A too large window won't be able to capture the small patterns, too small windows could do a better work, but this is not so well, multiple overlapping patterns would be found for large patterns Is there a more flexible way to address this problem of time-series patterns?


title: From one-dimension-SVR to multi-dimension-SVR (?) using MATLAB/LibSVM

text: this is my first post here. I already checked the Thread-Tags but I did not find any answer... My problem/question is how I can upgrade my SVR from one-dimension- to multi- dimension-regression. Its more a basic question. So let me introduce you into this. First I built a one-dim-SVR using MATLAB/LibSVM recording to an already existing Thread: LibSVM and matlab for time series prediction               x = [1:10]'; //full data set     y = [11:11:110]'; //full label set     TrainingSet = x(1:end-3);  //training data for SVR (70%)     TrainingSetLabels = y(1:end-3);  //training labels for SVR (70%)     TestSet = x(8:end);  //test data (30%)     TestSetLabels = [0;0;0]; //test labels      options = ' -s 3 -t 0 -c 100 -p 0.001 -h 0';     model   = svmtrain(TrainingSetLabels, TrainingSet, options);          [predicted_label, accuracy, decision_values] = ...         svmpredict(TestSetLabels, TestSet, model);      This is a simple code for one-dim-regression. With a data amount of 10 values. To train the machine it uses the first 7 values with its labels. After this the last 3 values are used for regression. With this code I get pretty good regression and accuracy. This works. Fine! * * * Now, I want to simply upgrade this into a multi-dim-regression. To understand this I want a 2-class-regression. Lets say each class has also 10 values like the first code above. And the values are the same (from 1 to 10). The labels are the products of coresponding values like:   1. labels X1 X2   2. 1 1 1   3. 4 2 2   4. 9 3 3   5. 16 4 4   6. 25 5 5   7. 36 6 6   8. 49 7 7   9. ? 8 8   10. ? 9 9   11. ? 10 10 The questionmarks representing the "unknown labels". (Of course we know them; not the machine) So I created this code... just by reading the LibSVM-Readme.               x = [1:10; 1:10]'; //full data set with two classes     y = [1; 4;9;16;25;36;49;64;81;100]; // full label set     TrainingSet = x(1:end-3, :); //training data set     TrainingSetLabels = y(1:end-3); //training label set     TestSet = x(8:end, :); //test data set     TestSetLabels = [0;0;0]; //test label set     options = ' -s 3 -t 0 -b 1 -c 100 -p 0.001 -h 0';     model   = svmtrain(TrainingSetLabels, TrainingSet, options);     [predicted_label, accuracy, decision_values] = ...         svmpredict(TestSetLabels, TestSet, model);      Here I get no single good regression and accuracy. Of course I used several Kernels. But there was no exceptable change. So here is my question/problem: Is my code wrong? Did I miss something? Why do I not get better regression? Please help me to get into this. Thank you!


title: Maximum likelihood for number of events in Poisson process

text: I have a Poisson process with parameter $\lambda$ known. How do I compute the maximum-likelihood estimator for $N$, ie. the number of events over a specific time spell $T$. To repeat, I know $\lambda$ and $T$, and I need to find the ML estimator for $N$ I am asking this because, given that the average time between Poisson events can be modelled as an exponential distribution with parameter $\lambda$, I am not even sure which of the two distributions I should use to write down the likelihood function and even if I solve this problem, I would not be able to take derivatives with respect to $N$.


title: The correct use of 'statistical significance'

text: Can chi-squared tests be used for making statements like this or is this an inaccurate analysis of my data? > The majority of the groups (n=9) affected by disease ranged in prevalence > from 5.56% to 12.73% with the exception of three groups (groups A, B and C) > which produced markedly higher prevalence rates. Chi-squared tests were > performed on these however, and it was determined that neither the 39.13% of > group A (χ2= 1.08, df=1, P=0.297) nor the 21.05% from group B (χ2= 6.36, df= > 1 P=0.12) were statistically significant. The 23.19% from group C (χ2= > 59.52, df=1, P=<0.05) was significant thereby marking it out as unique in > this assemblage. Is this acceptable? Or would a more accurate way be to do a chi-square on a 2x9 contingency table? And if it is, how then do you know which group is significantly different from the others?


title: The "mode" of sum of dependent random variables

text: Is there some study on finding the mode of sum of dependent random variales? Or give some bounded zone for the mode? Million thanks for your helps!


title: Setting significance levels for multiple testing

text: I have data from a number of data sets I am looking to find a correlation between measured volatiles. I have performed a t-test and have a p-value for each set of data. I now need to determine if I have a significant number of p<0.05. For example I have 40 sets with p<0.05 and 5 sets p>0.05, I now need to be able to say if this is significant, can you help?


title: Conditional Expectation of Poisson Random Variables

text: Suppose $X_1,X_2,\ldots,X_n$ is a random sample from a Poisson Distribution with mean $\theta$. How can I find the conditional expectation $E \left( X_1+X_2+3X_3 |\sum_{i=1}^n X_i \right)$? I know that $\sum X_i $ has a $poisson (n\theta$) distribution. Similarly the random variable $X_1+X_2+X_3$ has a $poisson (3\theta)$ distribution. I get confused with the required summations afterwards though. Thank you.


title: What does it mean to "fit a regression function" and then use it to update other functions?

text: Referring to the algorithm on page 11 in this paper on boosting algorithms, I really don't understand step 2, (ii) and (iii). What does this mean: > (ii) Fit the regression function $g_j^h (x)$ by weighted least-squares of > the working response $z_{ij}$ to $x_i$ with weights $w_{ij}$ on the training > data. > > (iii) Set $g_j(x) \leftarrow g_j(x) + g_j^h (x)$ For convenience, I'm pasting the algorithm here: ![gentleboost.c](http://i.stack.imgur.com/d5C6P.png) I do understand the theory and general concept of least squares regression. I understand we have data whose graph doesn't have a clear line or curve formula, and so we try to come up with a function that'll help us predict a response given an input. Using "least-squares" means coming up with a formula whose parameters are such that the square of the errors of the prediction is minimized. But once this regression function is "fitted" to data, how is it "added" to other functions like itself? I'm having a hard time understanding the application of the above quoted statement. Not just in this algorithm, but in others, wherever they mention "fitting" a regression function and then using that function with others. At best I understand using regression functions to actually predict values as output. But how do we use them to predict/generate more functions? Could someone please explain this simply, maybe using arrays of data or something?


title: Reporting contrasts between binary probability parameters in Bayesian data analysis - odds ratios or difference in probability?

text: In a bayesian data analysis, if one is modeling differences in binomial/bernoulli probability parameter differences between populations, is it still standard to report the difference in the binary probability as an odds ratio / relative risk as opposed to a probability difference? Take for example, estimating a set of binomial proportions and computing posterior contrasts between them via a beta-binomial model. The use of the odds ratio in the context of logistic regression makes sense since there's a direct correspondence between the coefficient and the change in the log odds ratio. However, in Bayesian data analyses, I commonly see the posterior distribution shown on a 0-1 scale and the posterior contrasts on the same scale as a posterior distribution. What's the preferred scale on which to present binary probability contrasts in a Bayesian data analysis (and why is it preferred)?


title: How to interpret coefficients in a regression with ARIMA errors?

text: I've got some time-series business data that I can fit relatively well with a `ARIMA(2,1,0)(1,1,0)[12]` model (using R's excellent `forecast::Arima` \-- thanks Prof. Hyndman!). The series is dominated by seasonal effects, but has trends as well, thus the differencing. I'm not an expert in forecasting. I'm exploring a future experiment (power analysis sorts of things) by simulating the effect of some sort of intervention that may increase (or decrease) the values in the series, probably multiplicatively. To do this, I'm scaling the numbers for the last N months by X%, and using the `xreg` parameter to change the model from a differenced AR(1) to a regression with time-series errors. The vector I'm using as the regressor looks like `[0, 0, ..., 0, 1, 1, 1]`, where 1s represent the months with the intervention in effect. The coefficient I get from the model appears to be an additive effect, which makes sense, but is _much_ smaller than the actual effect (4000 vs 100,000). However, when I use `forecast`, with and without 1s in the regressor, the difference is of the expected magnitude -- if anything, too high. So, my questions:   1. How do I interpret that coefficient. Is it additive?    2. Is it correct to use 1s in the regressor vector for all time periods that the treatment is in effect, or should I be thinking of this as an impulse that offsets the trend in the ARIMA model, and using a pattern like `0, 1, 0, 0, -1, 0`?   3. Any other advice? This is related to these questions, which either don't answer my question or I don't fully understand:   * Variables importance on regression with ARIMA errors model   * Intervention analysis in time-series regression with seasonal ARIMA errors   * How to detect a significant change in time series data due to a "policy" change? Thanks!


title: Computing autocorrelation from weighted samples

text: Suppose I've used Metropolis-Hastings to get samples from a candidate distribution, then associated with each an importance weight. (Perhaps the target distribution was too rowdy to use MH directly, but it seems to work fine on something similar.) I now have weighted samples from a stationary process. How would I compute time-lag autocorrelation? (I don't actually have weighted samples like that. I'm writing a statistical library in which everything can be done using weighted samples, and I ran into a snag.) I can think of a lot of hacks involving averaging the weights, but I'd like something better motivated than "it was the only thing I could think of". We'll assume that the weights are relevant, so I can't just chuck them.


title: Multiple Regression Model determination suggestions? (Time-series** or Cross-section)

text: What I have to model is the trend that occurs after a transient event (water pump failure, water pump repair etc.) where pump flow goes to 0 in order to forecast the pump flow after a planned event in the future. The most important things to capture are the area under the curve (Total water filtered) and then the peak water filtration. One way for me to describe the trend after these transient events is that the rate starts from 0, increases sharply to about the pre-transient event rate value, then continues to ramp up to a peak (since more water has built up around the pump since it was not pumping for some time) then trends back down and stabilizes at the pre-transient event rate. To describe this in numbers : Pre-transient event - Rate = 60 1:00 - Rate = 0 2:00 - Rate = 20 3:00 - Rate = 40 4:00 - Rate = 60 5:00 - Rate = 75 6:00 - Rate = 85 7:00 - Rate = 92.5 8:00 - Rate = 97.5 9:00 - Rate =~ 92.5 10:00- Rate =~ 87 11:00 - Rate =~ 85 ... 20:00 - Rate =~ 60 (Pre-transient rate) 21:00 - Rate =~ 60 ... **Data I have:** ~20 Filters ~3-4 individual filter events (repair, replacement etc) basically this mentioned trend 3-4 times. ~hourly data for the life of the water filter 1 Dependant variable ( Water Pump Rate ) 7-10 Independant variables ( pump rate before transient event, number of hours it takes to replace filter, temperatures, pressures, etc.) 1 common event that leads to a similar trend \- Can be a filter replacement \- Filter repair \- Filter All data is captured on a hourly basis of which the trend usually takes 15-20 hours to ramp up and return to it's pre-trasient position. I want advice on my approach listed below please. I am trying to best model the trend in order to predict the total water that is pumped during this post- transient event period, and also the length of time that this pump-flow period lasts (can be deduced if I can get the trend properly) for forecasting future events. My approach: 1) Since the pump rate data is noisy, smooth pre or post regression model creation with Sovitzky-Golay (retains peak and total pump volume). 2) Put all data into Matlab 3) Determine multiple linear regression equation with all variable combinations. 4) Use 5 comparisons to determine best model (SSE, adjusted R^2, Akaikes Information Criterion, Correlated AIC, Schwarz Bayesian info criterion) 5) Attempt to optimize by adding polynomials of these variables (^2,^3 maximum), or logs of variables. 6) Since independent variables (such as the temperatures/pressures) follow linear or otherwise predictable trends after a transient period, determine equations to fit these variables based on their values before the event and the length of time to repair. 7) Use these values in the model for forecasting the water pump rates going forward 100 hours into the future (4-5 days) I have just began to learn more about time-series and where this is a recurring event after each transient and often based on the pump rate before the transient and the length of time of the transient, I may be able to model it that way as well. I think that would be less reliable but could provide type-curves for this event. SORRY IN ADVANCE FOR LENGTH OF THIS POST. If anyone would better understand it with help of a spreadsheet of sample data, please don't hesitate to ask. @irishstat


title: PCA, ICA and Laplacian eigenmaps

text: I am very interested in the Laplacian eigenmaps method. Currently, I am using it to do dimension reduction for my medical data sets. However, I have run into a problem using the method. For example, I have some data (spectra signals), I can use PCA (or ICA) to get some PCs and ICs. The problem is how to get similar dimension reduced components of the original data? According to the Laplacian eigenmaps method, we need to solve the generalised eigenvalue problem, which is $L y = \lambda D y$ Here y is the eigenvectors. If I plot the eigenvectors, e.g., top 3 of the y vectors (set the solution according to 3 eigenvalues), the results are not interpretable. However, I can always plot the top 3 PCs and top 3 ICs, which somehow represent the original data x. I assume the reason is because the matrix L is defined by the weighting matrix (adjcent matrix W), and the data x has been fitted with the heat kernel to create W, which is using an exponential function. My question is how to retrieve the reduced components of x (not the eigenvector y of matrix L)? Thanks very much and I am looking forward to your reply. * * * Thanks very much for your guys replied. My data set is restricted and not easy to demonstrate the problem. Here I created a toy problem to show what I meant and what I want to ask. Please see the picture, Firstly, I create some sine waves A, B, C showing in red curves (first column of the figure). A, B, and C have 1000 samples, in other words, saved in 1x1000 vectors. Secondly, I mixed the sources A, B, C using randomly created linear combinations, e.g., $M = r1*A + r2*B + r3*C$, in which r1, r2, r3 are random values. The mixed signal M is in very high dimensional space, e.g., $M \in R^{1517\times1000}$, 1517 is randomly chosen high dimensional space. I show only first three rows of signal M in green curves (second column of the figure). Next, I run PCA, ICA and Laplacian eigenmaps to get the dimension reduction results. I chose to use 3 PCs, 3 ICs, and 3 LEs to do a fair comparison (blue curves showed as 3rd, 4th, and last column of the figure respectively). From the results of PCA and ICA (3rd, 4th column of the figure), we can see that we can interpret the results as some dimension reduction, i.e., for ICA results, we can recover the mixed signal by $M = b1*IC1 + b2*IC2 + b3*IC3$ (I am not sure if we can also get $M = a1*PC1 + a2*PC2 + a3*PC3$ with PCA results but the result seems quite right for me). However, please look at the results of LE, I can barely interpret the results (last column of the figure). It seems something 'wrong' with the reduced components. Also, I want to mention that eventually the plot of the last column is the eigenvector $y$ in formula $L y = \lambda D y$ Have you people got more ideas? Figure 1 using 12 nearest neighbours and sigma in the heating kernel is 0.5: ![Columns from left to right: original signal, mixed signal, PCs, ICs, LEs](http://i.stack.imgur.com/qahl5.png) Figure 2 using 1000 nearest neighbours and sigma in the heating kernel is 0.5: ![Columns from left to right: original signal, mixed signal, PCs, ICs, LEs](http://i.stack.imgur.com/z0XOi.png) Matlab codes with required package are uploaded to http://www.mediafire.com/?0cqr10fe63jn1d3 Thanks so much.


title: Image Clustering with K-means - Postprocessing

text: I did some clustering on an image (each pixel is an observation that has 5 variables associated with it), I get pretty detailed results but they are a little bit noisey... I think. I used K-means. Does anyone have a nice idea on how to reduce the noise a bit? anyone know some postprocessing for K-means. I would usually just apply a median filter to the image or something of the sort but I want to know if there is something a little nicer out there. Thank you in advanced. Not sure if posting this question here was the correct decision. Let me know. Julian. p.s. this was all done in python by the way, only brightly colored pixels were clustered. ![enter image description here](http://i.stack.imgur.com/29e3L.jpg)


title: Computing ridge regression with prior different from 0

text: I compute ridge regression results with Matlab, not using their implementation but simply computing `(trans(X)X)+kI)^-1+trans(X)y` as seen here. The given formula equals the ridge formula given in the link, just above "examples". "trans" stands for "transposed", "^" stands for "taken to the power of". "k" is equivalent to what is usually referred to as lambda. Ridge regression penalizes for moving away from 0, but I want to penalize for moving away from a certain prior. The prior is different for each coefficient that has to be computed. Is there a simple way modify the formula above to take the different priors into account? Thanks


title: Multiple linear regression on a data set with Python?

text: I'll preface my question with the fact that I'm just learning about linear regression so I may be thinking about this wrong. I have a set of data. In this set I have one dependent variable and about 10 independent variables and the data set is growing regularly. It's rows of data in database with 10 columns of independent variables and one column of the dependent variable. You can see my previous question for an example of what I'm trying to do: Variables importance: who can do the most pushups? The output of a linear regression is a formula right? Now I want to write a python script (I could use R also but I'd greatly prefer python) to take this data as input and output the linear regression formula. Is there a python method to do this? Do I need to run a regression comparing each independent variable with the dependent variable one at a time? Or is there a python method to feed in the data with all 10 independent variables and come out with a formula?


title: Logistic regression - the model is significant in predicting the DV, yet the percent correct decreases

text: How can this be? My thoughts are that while the percent correct decreased from 80.8% to 80% with the model, perhaps the model is regarded to be a significant predictor due to the specificity having increased. Any thoughts?


title: How to select best model out of linear, quadratic, and a model involving an exponential? Are the models nested?

text: I'm examining the growth rate of species en function of ion carbonate. For each specie, I have different measured of ion carbonate, it can be in different temperature or pH. For each measure, I'm doing graphics:   * x - ion carbonate;   * y - growth rate.  I know that for each measure I can have one of this models:   1. y=ax+b   2. y=ax^2+bx +c   3. y=a*(1-exp(-x/b)) + c  so I have to choose which one to select. I'm doing this with AIC. But AIC is only for nested models, is my 3 model is nested? 1 and 2 are nested I think. If not, how can I select what is the "best" model for my data?


title: Different F-ratios for within subjects effects when using SPSS and R's aov

text: I've just compared the ANOVA tables generated by SPSS and Statistica with the aov table provided by summary(aov.model). They yield identical between Subject effects (e.g., NativeLanguage(English vs Other), but different F ratios for within subject effects (e.g., Class(animate vs inanimate), aov F-ratios being consistently smaller and more conservative. The interaction of the within- between factors again yields identical terms. I am stumped. HOw can this be? Any suggestion? Here are more details: I was analyzing the RTs lexdec data base from the languageR package. I did some minor data adjustments. E.g., To get a sense for RTs I reversed the log transform exp(lexdec$RT) then removed Error responses and RT outliers. Using ddply I obtained condition means for NativeLanguage and Class for each subject (the data frame is shown at the bottom ob my post). Analyzing these data, I obtained different aov and STATISTICA summaries. Specifically for the within- subjects factor Class was p~.18 with aov and p~.06 with STATISTICA and SPSS, with larger Class SS values shown by the two commercial packages (569) than by aov (301). I've tried to make the two anova outputs (shown below) look transparent but it appears that the posting format does not match the format shown in the question window.               > C1.anova <- aov(RT ~ (Class * NativeLanguage)     +  + Error(Subject/Class) + (NativeLanguage), data=C1 )     > summary(C1.anova)          Error: Subject                    Df Sum Sq Mean Sq F value  Pr(>F)       NativeLanguage  1  81413   81413  6.2973 0.02131 *     Residuals      19 245637   12928                       ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Error: Subject:Class                          Df  Sum Sq Mean Sq F value   Pr(>F)        Class                 1  301.51  301.51  1.9661 0.176994        Class:NativeLanguage  1 2175.86 2175.86 14.1880 0.001305 **     Residuals            19 2913.83  153.36                         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1       STATISTICA                                      SS  DF        MS         F          p     Intercept        15387240   1  15387240  1190.202   0.000000     NativeLanguage      81413   1     81413     6.297   0.021311     Error              245637  19     12928          Class                 569   1       569     3.709   0.069217     Class*NativeLanguage 2176   1      2176    14.188   0.001305     Error                2914  19       153           Data:                  Subject  Class NativeLanguage       RT     1       A1 animal        English 557.6410     2       A1  plant        English 548.4687     3       A2 animal        English 533.4737     4       A2  plant        English 511.7941     5       A3 animal          Other 598.9545     6       A3  plant          Other 602.4118     7        C animal        English 562.8864     8        C  plant        English 560.0588     9        D animal          Other 630.1464     10       D  plant          Other 604.0286     11       I animal          Other 542.1219     12       I  plant          Other 533.1666     13       J animal          Other 565.4324     14       J  plant          Other 513.2333     15       K animal        English 492.4500     16       K  plant        English 517.7333     17      M1 animal        English 481.8372     18      M1  plant        English 497.9687     19      M2 animal          Other 671.6666     20      M2  plant          Other 655.8750     21       P animal          Other 640.7209     22       P  plant          Other 610.0286     23      R1 animal        English 552.9744     24      R1  plant        English 545.4242     25      R2 animal        English 636.8864     26      R2  plant        English 675.1714     27      R3 animal        English 607.8572     28      R3  plant        English 614.9428     29       S animal        English 599.9285     30       S  plant        English 586.6286     31      T1 animal        English 580.0500     32      T1  plant        English 583.2857     33      T2 animal          Other 892.5526     34      T2  plant          Other 862.1000     35       V animal          Other 736.2619     36       V  plant          Other 718.3529     37      W1 animal        English 517.0465     38      W1  plant        English 539.2727     39      W2 animal        English 639.1363     40      W2  plant        English 666.7143     41       Z animal          Other 725.3750     42       Z  plant          Other 706.2069


title: Whether and how to perform VIF and Condition Index on formative measurement model?

text: I found this very interesting paper which I think provides some good guidance on testing the reliability of my **formative measurement model** \- see earlier post. The relevant information is on page 6 of the paper and I am quoting it below (because I know external links are unstable and should not be used for ongoing reference): **Extract from:** Phatcharee Toghaw Thongrattana, "Assessing reliability and validity of a measurement instrument for studying uncertain factors in Thai rice supply chain" (October 1, 2010). SBS HDR Student Conference. Paper 4. http://ro.uow.edu.au/sbshdr/2010/papers/4 > _**_*Reliability test for Formative constructs*** > > As formative constructs composited of different aspects of a construct that > their indicators are not necessary to correlate with each other > (Diamantopoulos and Winklhofer, 2001). Straub, Boudreau et al. (2004, p.400) > state that “it is not clear that reliability is a concept that applies well > to formative constructs”. This statement is also supported by Diamantopoulos > and Siguaw (2006, p.270) and Rossiter (2002, p.315) that no dimensionality > and reliability test are performed on formative indicators because factorial > unity in factor analysis and internal consistency are not relevant. > Although, low item-to-total correlation should be dropped from measurement > scales to increase internal consistency reliability for reflective > measurement model because the scales are from the same content construct, > the removal of measurement scales in formative measurement model can lead to > change the empirical and conceptual meaning (MacKenzie et al., 2005). > Andreev, Heart et al. (2009) conclude that construct reliability of > formative should be performed by multicollinearity, test of indicator > validity (path coefficients significance), and optionally, if appropriate, > testretest (Petter et al., 2007). > > *On the other hand, reflective constructs that multicollinearity among items > in the same construct is desirable such as high Cronbach’s alpha, but > reliability of formative construct in term of multicollinearity is not > present because if multicollinearity is present, it means that indicators > are tapping into the same aspect of the construct (Petter et al., 2007, > p.641). Likewise, formative measurement model is based on a multi-regression > that multicollinearity should not exist (Diamantopoulos and Winklhofer, > 2001). Thus, reliability evaluation for formative constructs is to assess > the assumption of no multicollinearity (Diamantopoulos and Siguaw, 2006). > Variance Inflation Factor (VIF) is evaluated. There are some guidelines that > can be applied: > > • VIF is less than 3.3 that shows a excellent value (Diamantopoulos and > Siguaw, 2006). • VIF is less than 10 that no collinearity is commonly > accepted (Hair et al., 1995). > > As collinearity also can be harmful effects to formative constructs, > condition index is the standard diagnostics that measure the relative amount > of variance associated with an eigenvalue. Its threshold value should be > less than 30 to find no support for the existence of collinearity (Hair et > al., 1995). If multicollinearity exists, Petter, Straub et al. (2007, p.642) > recommended that at first, the model construct may have both formative and > reflective measures. Secondly, the correlated measurement items can be > removed, if content validity is not affected. Thirdly, the correlated > measurement items can be collapsed into a composite index. Lastly, it can be > converted into multidimensional construct. ### My questions   * How can I perform VIF and Condition Index using SPSS?   * Is VIF and Condition Index suitable for assessing reliability of formative constructs?


title: Determine separation between two modes from distribution

text: I’ve got a sample of pairwise distances between points in a 2D picture. Some of these points lie within the same object. Their distance to each other is thus smaller than some well-defined threshold (the object’s diameter). Points that lie in different objects (predominantly) have a pairwise distance greater than said threshold. Points that lie within the same object are however rare (<10%). I would like to determine this distance threshold empirically from my sample. For “appropriate” parameters (well, herein lies the rub, doesn’t it?) the threshold is visible in the density plot: ![enter image description here](http://i.stack.imgur.com/qqd75.png) The threshold is marked by the arrow. This is the _objectively right_ cut-off for my application: it is the dip after the first tall plateau which corresponds to the distribution of the few points lying within the same object, and it corresponds to the object diameter that can be individually verified in the original picture, but not easily automatically deduced from my data. Unfortunately, I have no idea how to determine it in an automated fashion. Even the `adjust` argument / bandwidth for the density function has been found by trial and error, and a different input data set I’ve tried requires a different bandwidth. Is there _any_ hope? Or should I just give up?


title: Fitting distributions, goodness of fit, p-value. Is it possible to do this with Scipy (Python)?

text: INTRODUCTION: I'm a bioinformatician. In my analysis which I perform on all human genes (about 20 000) I search for a particular short sequence motif to check how many times this motif occurs in each gene. Genes are 'written' in a linear sequence in four letters (A,T,G,C). For example: CGTAGGGGGTTTAC... This is the four-letter alphabet of genetic code which is like the secret language of each cell, it’s how DNA actually stores information. I suspect that frequent repetations of a particular short motif sequence (AGTGGAC) in some genes are crucial in a specific biochemical process in the cell. Since the motif itself is very short it is difficult with computational tools to distinguish between true functional examples in genes and those that look similar by chance. To avoid this problem I get sequences of all genes and concatenated into a single string and shuffled. The length of each of the original genes was stored. Then for each of the original sequence lengths, a random sequence was constructed by repeatedly picking A or T or G or C at random from the concatenated sequence and transferring it to the random sequence. In this way, the resulting set of randomized sequences has the same length distribution, as well as the overall A,T,G,C composition. Then I search for the motif in these randomized sequences. I permormed this procedure 1000 times and averaged the results.   * 15000 genes that do not contain a given motif    * 5000 genes that contain 1 motif    * 3000 genes that contain 2 motifs    * 1000 genes that contain 3 motifs    * ...    * 1 gene that contain 6 motifs So even after 1000 times randomization of true genetic code, there aren't any genes which have more than 6 motifs. But in the true genetic code, there are a few genes which contain more then 20 occurrences of the motif, which suggest that these repetition might be functional and it's unlikely to find them in such an abundance by pure chance. PROBLEM: I would like to know the probability of finding a gene with let's say 20 occurences of the motif in my distribution. So I want to know the probability to find such a gene by chance. I would like to implement this in Python, but I don't know how. Can I do such an analysis in Python? Any help would be appreciated.


title: How to test the predictive power of a model?

text: I want to build a model to predict the outcomes of experiments. My predictive model gives out scores with an range 1 to 100 values. I want to test if my predictive scores can be used to classify experimental outcomes as "good" or "bad" groups. Experimentally, we did the 1000 experiments. Using my predictive model, I have 1000 scores. To test if my predictive model statistically acceptable, what should I do? I have done ROC and sensitivity test for these 1000 X 2 data. ROC were plotted for all 1000 experimental data and predictive scores. By looking at the AUC values for the plot (sensitity vs 1-specificity), AUC=0.64. Let's said if my predictive score has a cut off value of 5, i.e. it is likelihood that the experimental outcome will be "good", score > 5 are likely to have "bad" experimental outcome. I calculate the enrichment of my predictive model, i.e. no. of real "good" results / no. of predictive score < 5. Did I do anything wrong here? What else should I do to check the predictive power of a model?


title: How to create one synthetic variable from 5 measured variables

text: I have 5 more less correlated variables measuring conceptually same thing - size. For example, height, weight, shoe size, jacket size and age. I just want to summarize the information from all variables into one and use that single (synthetic) variable in further analysis. What are the options for doing that? One of the requirements is to be able to Measure and "tweak" weights of each of the components of the synthetic variable based on exPert judgement... lats say the requirement may be to increase weight of shoe size from 11.34% to 20%. Or reduce weight of another variable because of known measuremen error or some other valid reason. Any ideas?


title: Does the estimated overdispersion parameter of Negative Binomial depend on mean

text: Negative Binomial distribution can be parameterized using mean, $\mu$, and overdispersion $\psi$, so that the variance of NB is $\mu + \frac{\mu^2}{\psi}$. We know there is no analytical solution for estimating $\psi$. I understand the variance of NB depends on mean $\mu$. But does the estimated $\psi$ also depend on mean. Say using MLE or any other commonly used method.


title: How can an improper prior lead to a proper posterior distribution?

text: We know that in the case of a proper prior distribution, $P(\theta \mid X) = \dfrac{P(X \mid \theta)P(\theta)}{P(X)}$ $ \propto P(X \mid \theta)P(\theta)$. The usual justification for this step is that the marginal distribution of $X$, $P(X)$, is constant with respect to $\theta$ and can thus be ignored when deriving the posterior distribution. However, in the case of an improper prior, how do you know that the posterior distribution actually exists? There seems to be something missing in this seemingly circular argument. In other words, if I assume the posterior exists, I understand the mechanics of how to derive the posterior, but I seem to be missing the theoretical justification for why it even exists. P.S. I also recognize that there are cases in which an improper prior leads to an improper posterior.


title: Aggregating time series data to one variable

text: I have got various time series data for different observations. I want to represent one time series observation by **one aggregated variable**. The easiest one is to take the `mean` of the time series, but `mean` does not capture the variations in the observation over time. I thought of computing the distance (`dist`) between an arbitrary reference and the time series data. I am curious to know if there is other technique that I could apply to fairly represent the time series observation by **one aggregation variable**? Is there any R package that I can use? Thanks


title: Clustering and Feature Selection For Audio Data

text: I'm very new to stats and I'm at a loss as to where to start with this problem. I'm not sure what tools or methods to use to extract meaningful answers from my data. I'll try and describe the problem generally Basically, I've extracted about 22 features or so from individual audio tracks of varying instrumentation(drums, guitars etc) on a frame by frame basis of each track, so they're effectively a time series of features in MATLAB. I've selected these features due to various literature I've read. I know through listening that some of these tracks are similar, so they must therefore have common features. I want go through all my audio and identify which tracks have similar features to each other and identify what features makes them so similar. For example I want to know what are the common features to all the drum tracks I have. I effectively want to identify which are the best features to group the instruments by. What kind of stat analysis do I need to do try and figure this out? Linear regression, time series analysis, Dynamic Time Warping?? I have no idea where to start. My knowledge of statistical methods is poor.


title: Applying ARMAX model from r output

text: I'm trying to apply R output to generate a scenario using external data, I'm not sure how exactly to use the coefficients in each from the R output. I have an ARMAX(1, 1) model   * Coefficient of AR1: $A$   * Coefficient of MA1: $M$   * Intercept = $k$   * Coefficient of external regressors: $B_1$, $B_2$   * External regressors: $X_1$, $X_2$   * Y actual = $Y_a$   * Y predicted = $Y_p$ The formula I'm using is $Y_p = k + A*(Y_a(t-1)-k) - M*(Y_a(t-1) - Y_p(t-1)) + B_1*(X_1(t) - X_1(t-1)) + B_2*(X_2(t) - X_2(t-1))$ Is there anything wrong with the formula I'm using? I'm asking because of two things, one is that while AIC is much better for this model, the sum of squared errors is actually higher than if I just used AR(1) model; another thing is that when I adjust the AR term, the sum of squared residuals gets reduced, but I thought the model is supposed to minimize the sum of squared residuals for each of the terms? I'm using the same data set for generating the model and to fit the model, I need to test to see how much improvement this model offers compared to simpler models. I'm doing the fitting in Excel. Any help would be greatly appreciated, and if additional information is needed to answer this please ask. This is my R output               ARMA11R2R30               Call:     arima(x = YieldReOLD[, 5], order = c(1, 0, 1), xreg = cbind(YieldReOLD[, 2],          YieldReOLD[, 4]))          Coefficients:              ar1      ma1  intercept  cbind(YieldReOLD[, 2], YieldReOLD[, 4])1           0.9872  -0.1970    -6.2862                                   -0.1867     s.e.  0.0072   0.0489     0.3743                                    0.0566                cbind(YieldReOLD[, 2], YieldReOLD[, 4])2                                            -0.3999     s.e.                                    0.1140


title: Finding the level of significance of a sample

text: I have an exercise with the following information; $ 20\% $ of the set are youth, my sample of $ 139 $ persons have $ 39 $ youths. The exercise is to test if this sample is a representative sample of the set and to find the level of significance of the result. I have two hypothesis $ H_0 = 0.20 $ and $ H_1 \ne 0.20 $, sample statistics $ \hat{p} = 0.2955 $ and the test statistic $ z = 2.74 $ the p-value: $ 2 * P(z\geqslant 2.74) = 0.0062 $ How do I go from this to finding the level of significance? Help is much appreciated.


title: Guess a probability distribution function that take latest or more recent values from data

text: Is there any good source or link that has several probability distribution functions (pdf) we can use based on the simple observations or patterns we observe from data. The question I have in my mind is to adopt some distribution function of prospective customers of car sales based on the quotations one enter from car- dealer website (a car sale is dependent on several features eg: Car Type -sedan or SUV, Engine kind, fully loaded etc). If I hypothesize there is more likely chance that customer with latest quotation (or recent quotation) to buy the car with the last entered set of features. Is there any pdf that can describe the population behavior?


title: Why do I always fail to reject the null hypothesis on a binomial test for `p=.5`?

text: When performing a two-tailed test on a binomial process with p=.5, why are my p-values all over the place, and I virtually never reject the null hypothesis? I answered a question on stackoverflow about testing a random process (i.e., that a function designed to yield a binomial random variable was working correctly). I got output like this:               can reject the hypothesis that the following tests are NOT the     results of a binomial process (with their given respective     probabilities) with probability < .01, 1000000 trials each     p = 0.01 {False: 10084, True: 989916} 4.94065645841e-324 reject null     p = 0.1 {False: 100524, True: 899476} 1.48219693752e-323 reject null     p = 0.33 {False: 100633, True: 899367} 2.96439387505e-323 reject null     p = 0.5 {False: 500369, True: 499631} 0.461122365668 fail to reject     p = 0.66 {False: 900144, True: 99856} 2.96439387505e-323 reject null     p = 0.9 {False: 899988, True: 100012} 1.48219693752e-323 reject null     p = 0.99 {False: 989950, True: 10050} 4.94065645841e-324 reject null      My p-values are always all over the place for p = 0.5, but fairly consistent (almost too consistent) on the others. I'm sure if I went back into the math, I could prove to myself why this is, but it's been too long since I've studied this, and I just don't have the time right now. So who can tell me why this is, or if perhaps I've made a mistake here. ## Edit: Conclusion I will consider this fully answered: it appears the p=.5 _is_ the null hypothesis for the test I was performing. So, semantically my problem now is how do I test a null p != .5 (i.e. two-tailed)? To those who need a analogy, say I want to test a fair coin. I have found a specific answer on the site addressing this new problem of mine. Thank you so much for participating in this exercise. Here's the link: Testing if a coin is fair


title: How to normalize ranked data

text: I am doing some machine learning and need help with the stats aspect of my problem. I have a number of addresses of webpages and some features for these webpages. I am running TF-IDF on the webpage text to create more features. One of the features I have extracted is, for each web address, I retrieve the Google Page ranking. This value can be any value in the world,but obviously, the lower the "better quality" Google has deemed it to be. How can I normalise this figure, given that I have 7,000 web addresses and the ranks can vary enormously (www.google.com, for instance, may be ranked #1, while www.bbc.co.uk may be #1,117). We can assume my Logistic Regression is merely trying to predict whether a webpage is "good" or not. I am doing this just as a learning exercise. I'd also appreciate a quick explanation of how TF-IDF creates more features, I cannot seem to find this anywhere (though I understand the general concept of putting more weight on less common terms and less weight on more common terms). Thanks very much!


title: Strength of association test with binary variables

text: I have a dataset with different purchases for two different items from the same users. So the users purchased the two items at different points in time. I also have 3 different variables: `High`, `Mid`, and `Low` twice (Low, Mid,High for each product) within this dataset. So the product are broken up into these categories. For example, suppose I bought item X and its classified as a High price item. I then want to see if the same user bought an item Y at the same classification level (High). My end goal is to see if there is an association of a customer buying items at the same classification level across different items. My first thought was to use a Chi square test but I am not sure this is the best way to do it. I am doing this in R. All of the data are binary variables. Here is an example of what the dataset looks like.               User      Type             High  Mid  Low             1         Product 1        0     1    0            1         Product 2        1     0    0            2         Product 1        0     1    0     2         Product 2        0     1    0     3         Product 1        0     0    1     3         Product 2        0     1    0


title: Methods for 3D histograms comparison

text: When I wish to compare 2D image histograms I can use methods like Chi-Square and Intersection but what are my options if I wish to compare 3D histograms (e.g, based on R,G,B values)? Thanks.


title: Lasso cross validation

text: I want to perform cross validation to find the regularization parameter for Lasso. I am using scikit-learn library in python. I first generate the dataset and then perform k-fold cross-validation. Here is my code (most of it from an example at scikit-learn website):               # generate some sparse data to play with     import numpy as np     n_samples, n_features = 5000, 200     X = np.random.randn(n_samples, n_features)     coef = 3 * np.random.randn(n_features)     coef[10:] = 0  # sparsify coef     y = np.dot(X, coef)          # add noise     y += 0.01 * np.random.normal((n_samples,))          # Split data in train set and test set     n_samples = X.shape[0]     X_train, y_train = X[:n_samples / 2], y[:n_samples / 2]     X_test, y_test = X[n_samples / 2:], y[n_samples / 2:]          ###############################################################################     # Lasso     from sklearn.linear_model import Lasso     from sklearn.cross_validation import KFold     from matplotlib import pyplot as plt          kf = KFold(X_train.shape[0], n_folds = 10,)               alphas = np.logspace(-16, 3, num = 50, base = 2)          e_alphas = list()     e_alphas_r = list()  #holds average r2 error     for alpha in alphas:         lasso = Lasso(alpha=alpha)         err = list()         err_2 = list()         for tr_idx, tt_idx in kf:             X_tr , X_tt = X_train[tr_idx], X_test[tt_idx]             y_tr, y_tt = y_train[tr_idx], y_test[tt_idx]             lasso.fit(X_tr, y_tr)             y_hat = lasso.predict(X_tt)             err_2.append(lasso.score(X_tt,y_tt))             err.append(np.average((y_hat - y_tt)**2))         e_alphas.append(np.average(err))         e_alphas_r.append(np.average(err_2))          plt.figsize = (15,10)     fig = plt.figure()          ax = fig.add_subplot(111)     ax.plot(alphas, e_alphas, 'b-')     ax.plot(alphas, e_alphas_r, 'g--')     ax.set_xlabel("alpha")     plt.show()      The graph of error is show in the figure at below: ![Lasso CV](https://dl.dropboxusercontent.com/u/943457/lasso_cv.png) I know that there are other ways in scikit-learn to do a lassoCV but I just want to know how do you select the parameter given that kind of graph I am getting. Thanks for your reply.


title: Maximum likelihood factor analysis fails

text: When trying to use maximum likelihood Factor analysis instead of Principle component analysis in SPSS, what are the reasons that I am unable to complete a factor extraction using maximum likelihood. When I attempt to perform it, it says `"Attempted to extract factors, no local minimum was found. Extraction terminated"`.


title: Compute percentage variance explained with NIPALS

text: I am using NIPALS (Non-linear iterative partial least squares) to compute first $K$ components of a $N\times M$ data matrix $X$ with $N$ samples and $M$ features. Along with the principal components NIPALS returns the eigenvalues of the first $K$ components. I would like to compute the percentage variance explained given by $$ \frac{\lambda_k}{\sum_{j=1}^M \lambda_j} $$ However I do not know the eigenvalues for $j>K$ and consequently I can't compute the denominator. Is there a solution to this? Here is my rough idea how to compute the denominator. $$ \sum_j \lambda_j = \mathrm{Trace}(X X^T) = \sum_{i=0}^M x_i x_i^T $$ where $x_i$ is the $i$th column of $X$. Can someone please tell me whether this idea is correct or not? I'm concerned because I saw that the eigenvalues computed by the SVD-PCA are different in magnitude than those computed by NIPALS. However SVD-PCA works with $X^TX$ covariance matrix while Nipals works with $X X^T$, right? So this is to be expected and as long as I compute the trace of $X X^T$ the above idea should work fine.


title: Translating R code about treatment into Effect size. Expected mean and variance

text: Using R, I created groups of individuals with trait values. Then I simulated a treatment that modified their trait value (see below). Finally I run a one-way Anova on them using the individuals traits value as dependent variable and the groups as independent variable. I want to know how does my **treatment magnitude can be translated into Effect Size** in order to compare my results with what I would expect using the pwr package. I don't want to measure the effect size afterward, I want to know what effect size did I simulated. I'm not sure of which words I should use, I might have better to call it the expected effect size given the treatment I simulated. Here is the definition of effect size for Anova as defined by Cohen.               Effect.size = sqrt(Epsilon.over.all.i(p.i *(mu.i - mu)^2)/sigma.square)         mu.i = the mean of the group i         mu = the mean of the overall         p.i = n.i / N         n.i = the number of individuals in group i         N = the total number of individuals         sigma square = the error variance within groups      You'll find it more explicitly written directly on my source (one third of the page in the ANOVA section) Here is my R script:               # Create individuals and groups.         N=100 # sample size     nb.groups = 4     sd=1     My.data = data.frame(     group=rep(1:nb.groups,n/nb.groups),     Trait=rnorm(n=N, mean=0,sd=sd)     )          # Simulate a treatment          Treatment = 0.55      Way 1 of simulating a treatment               # Adding a fixed value if Treatment = 4, nb.groups=4:       # group 1: nothing is added       # group 2: 1/3 * 4 is added       # group 3: 2/3 * 4 is added       # group 4: 3/3 * 4 is added          for (i in 2:nb.groups){          My.data[which(My.mat[,1]==i),2] =           My.data[which(My.mat[,1]==i),2] + Treatment*((i-1)/(nb.groups-1))     }      Way 2 of simulating a treatment               # Adding a value drawn from a random distribution which mean equals the fixed value we added in Way 1 and sd always equals my parameter "sd"          for (i in 1:nb.groups){          My.data[which(My.mat[,1]==i),2] =          My.data[which(My.mat[,1]==i),2] + rnorm(N/nb.groups,mean=Treatment*((i-1)/(nb.groups-1)),sd=sd)     }      I repeat my question: Both for way 1 and way 2, I'm trying to find a function that "translate" the parameter "Treatment" into the common "Effect size". In order to use the above formula for effect-size, I must be able to infer the expected mu.i, mu and sigma.square but I don't know how to do that ! Here is the ANOVA I ran:               aov(my.data$Trait~my.data$group)      P.s. I already posted this question but did not get any answer or comment. So I try again. Please let me know if my question is unanswerable.


title: Count data as an independent variable in OLS- using a dummy variable+ the variable linearly to account for skewness

text: I am using OLS to model the relationship between amount of foreign aid (dependent variable, logged) and media coverage (number of newspaper articles, count variable). I assume a linear relationship between the two and use the "media coverage" variable as a continuous predictor. There is a spike in 0, however, in the count variable which makes it highly skewed to the right. Although this problem was addressed in different threads, I would like to better understand why using a dummy variable alongside the original variable would improve the model. What does the dummy variable do for the skewness? how can it be interpreted in parallel with the continuous part of the model? Moreover, do the values of 0 in the count variable have to be kept or replaced as "missing values" ? I am also using a second IV, measuring the proportion of negative articles in the total amount of articles, which naturally spikes at the value of 0 as well. Does the dummy variable control for the zeros in this variable as well? Thank you in advance!


title: Negative predictions for binomial predictions from gbm in R

text: I've just fit a binomial model (training y = 0 or 1) using R's gbm package. When I calculated predicted values using my validation data, some of the predicted values were less than 0. Is this normal behavior? If so, how do I handle the negative values, floor them at 0?


title: Modelling access to data in a storage and replication system

text: I need to model data access for a storage system where I have:   * $S = $ Number of Server Nodes (that hand out data to clients and store data)   * $D = $ Number of Data Nodes (that do not hand out data to clients and only store data)   * $F = $ Number of Frontend Nodes (that hand out data to clients and do not store data)   * $R = $ Replication factor (copies of the data, no more than one copy per node) For my purposes, either $D$ or $F$ will be $0$ (not sure if that affects the formula but it makes thinking about it easier). What I am trying to model is the probability that on any given I/O operation (read or write), the $S$ or $F$ node being accessed by a client will have NO local copy of the data.


title: Ljung-Box Statistics for ARIMA residuals in R: confusing test results

text: I have a time series I am trying to forecast, for which I have used the seasonal ARIMA(0,0,0)(0,1,0)[12] model (=fit2). It is different from what R suggested with auto.arima (R calculated ARIMA(0,1,1)(0,1,0)[12] would be a better fit, I named it fit1). However, in the last 12 months of my time series my model (fit2) seems to be a better fit when adjusted (it was chronically biased, I have added the residual mean and the new fit seems to sit more snugly around the original time series. Here is the example of the last 12 months and MAPE for 12 most recent months for both fits: ![fit1, fit2 and original data](http://i.stack.imgur.com/kkUOb.png) The time series looks like this: ![original time series](http://i.stack.imgur.com/twNkT.png) So far so good. I have performed residual analysis for both models, and here is the confusion. The acf(resid(fit1)) looks great, very white-noisey: ![acf of fit1](http://i.stack.imgur.com/gyIv3.png) However, Ljung-Box test doesn't look good for , for instance, 20 lags:                   Box.test(resid(fit1),type="Ljung",lag=20,fitdf=1)      I get the following results:                   X-squared = 26.8511, df = 19, p-value = 0.1082      To my understanding, this is the confirmation that the residuals are not independent ( p-value is too big to stay with the Independence Hypothesis). However, for lag 1 everything is great:                   Box.test(resid(fit1),type="Ljung",lag=1,fitdf=1)      gives me the result:                   X-squared = 0.3512, df = 0, p-value < 2.2e-16      Either I am not understanding the test, or it is slightly contradicting to what I see on the acf plot. The autocorrelation is laughably low. Then I checked fit2. The autocorrelation function looks like this: ![acf fit2](http://i.stack.imgur.com/JZ7Sc.png) Despite such obvious autocorrelation at several first lags, the Ljung-Box test gave me much better results at 20 lags, than fit1:                   Box.test(resid(fit2),type="Ljung",lag=20,fitdf=0)      results in :                   X-squared = 147.4062, df = 20, p-value < 2.2e-16      whereas just checking autocorrelation at lag1, also gives me the confirmation of the null-hypothesis!                   Box.test(resid(arima2.fit),type="Ljung",lag=1,fitdf=0)         X-squared = 30.8958, df = 1, p-value = 2.723e-08       Am I understanding the test correctly? The p-value should be preferrably smaller than 0.05 in order to confirm the null hypothesis of residuals independence. Which fit is better to use for forecasting, fit1 or fit2? Additional info: residuals of fit1 display normal distribution, those of fit2 do not.


title: What is the literal interpretation for having a very low power (~30%)?

text: I am going to be conducting an experiment where I have two discrete samples. I have a limited number of samples, and when I calculate the Power of the test, (given a set number of samples, and difference in the proportions that would be significant), it is low (~30%). I am not a statistics guru, but it seems to me this experiment isn't even worth running. What are the practical implications of having very low Power? I understand that Power is the probability that the test will reject the null hypothesis, given that the null hypothesis ought to be rejected, but I am still having a hard time wrapping my head around what this means. Is it equivalent to saying that any difference we see between the two samples is with the "noise" of the experiment?


title: Testing a model

text: So I hear the term 'testing a model' used a lot. There seems to be a variety of ways to find whether a model has explanatory power (R^2) or if there is an overall significant relationship between the variables (g-test). However, let's say I have data which I acquire each day. Let's say I ran a logistic regression on the past three days of day. Now, I want to know if that model has any explanatory power in explaining data from the newest day. So does this model still hold up if I look at just the newest day of data? I'm just wondering if there's a better way to test a model this way than by simply running that same model on the new dataset?


title: issue in using GP_fit() of GPfit package in R

text: I am trying to work with Gaussian process through `GPfit` package of R but every time I use `GP_fit()` of this package, either it gets hanged or throws error which is like this:               Error in chol.default(Sig) :        the leading minor of order 2 is not positive definite.       My input matrix looks like this:                    visitLength       hour        [1,] 0.0241057543 0.73913043        [2,] 0.0413780103 0.91304348        [3,] 0.0432866771 0.30434783        [4,] 0.0000000000 0.17391304        [5,] 0.0639756822 0.60869565        [6,] 0.1239219567 0.86956522        [7,] 0.0000000000 0.34782609        [8,] 0.0389038126 0.17391304        [9,] 0.0886469674 0.04347826          .      .             .          .      .             .      and so on. There are about 200 rows in this input matrix. My output matrix looks like this:                    [,1]        [1,]    0        [2,]    0        [3,]    1        [4,]    0        [5,]    1        [6,]    0        [7,]    0        [8,]    1        [9,]    0      and this too has same number of rows as input matrix. Now when I use:               model = GP_fit(input_mat,outpt_mat)      either it gets hanged or the above mentioned error is thrown. I have no idea why it gets hanged when there are only 2000 rows and also I have scaled values in range [0,1] or it throws as error which does not clarify the meaning or reason. Can anyone provide some information about what and why is it happening? **UPDATE:** The signature of GP_fit() is:               GP_fit(X,Y)       where               X is the (n x d) design matrix     Y is the (n x 1) vector of simulator outputs.      One of the **example code** given **in the package documentation** is like this:                library(GPfit)       n = 5; d = 1;         computer_simulator <- function(x){         x = 2*x+0.5;         y = sin(10*pi*x)/(2*x) + (x-1)^4;         return(y)         }         set.seed(3);         library(lhs);         x = maximinLHS(n,d);         y = computer_simulator(x);         GPmodel = GP_fit(x,y);         print(GPmodel)


title: Nonparametric credible intervals

text: I have a small ( _n_ ≈ 15) IID sample from a continuous distribution that I have no reason to believe (a priori or a posteriori) belongs to a familiar parametric family. I do know a priori that the distribution is bounded, and I know its bounds, but that's about all I know. I'd like to construct 95% credible intervals for the population mean and standard deviation. I'm at a loss as to how to do this because it appears that credible intervals require a posterior PDF, which requires a likelihood function, which requires a parametric model. The term "nonparametric Bayes" seems to refer to techniques that use a parametric model whose number of parameters is allowed to grow with the sample size. I don't think that's what I want. I can imagine using something like the set of all piecewise linear functions with integral 1 and _n_ pieces as a parametric model, but a bit of preliminary work in this direction makes me think it would be very messy. If what I want here isn't really possible, then what would a conscientious Bayesian do in my situation? Is there perhaps a way to frame the problem so that I can say that I'm 95% sure the population mean is in some 95% _confidence_ interval?


title: Sum of Products of Rademacher random variables

text: Let $x_1 \ldots x_a,y_1 \ldots y_b$ be independent random variables taking values $+1$ or $-1$ with probability 0.5 each. Consider the sum $S = \sum_{i,j} x_i\times y_j$. I wish to upper bound the probability $P(|S| > t)$. The best bound I have right now is $2e^{-\frac{ct}{\max(a,b)}}$ where $c$ is a universal constant. This is achieved by lower bounding the probability $Pr(|x_1 + \dots + x_n|<\sqrt{t})$ and $Pr(|y_1 + \dots + y_n|<\sqrt{t})$ by application of simple Chernoff bounds. Can I hope to get something that is significantly better than this bound ? For starters can I at least get $e^{-c\frac{t}{\sqrt{ab}}}$. If I can get sub-gaussian tails that would probably be the best but can we expect that (I dont think so but can't think of an argument)? Thanks in advance


title: When should you consider using time series in analysis?

text: I have a number of projects that I am working on where data is collected over time periods of months to years. My first question is when should time be considered as a factor in conducting analyses and tests on the data? E.g. I have a list of 100 individual products and track view and purchase events which occur infrequently over a period of 3 months. I wish to compare the ratio of views to purchases to see if one product is purchased more frequently than another based on its popularity (views) then should I consider time as a factor in this analysis? My second question is what is the best test to use to check for difference between ratio of views and purchases of these products? Many thanks in advance for any assistance.


title: Nonparametric/nonlinear regression

text: I am looking for a survey/book on some state-of-the-art non-parametric (or nonlinear regression) methods, preferably with an inclination towards sequential data. Till date I have used `gaussian process` regression and some others such as `knn` regression, `Random Forest` regression, etc. I am looking for other possible regressors for example `particle filtering` based regression, etc. As mentioned earlier a recent survey would be great


title: What is the minimum sample size for using exploratory factor analysis to reduce a pool of questionnaire items?

text: **Context:** In my real experiment I am planning to include a questionnaire. I aim to measure 4 different constructs with multiple questions per construct. The questionnaire now consists of 24 items. I created more questions than I need with the idea of doing a pre-test and then reducing the number of questions to 12 based on the results of factor analysis. I was thinking of having a sample size of 30 people to participate in my pre- test, but is just what I came up with. **Question:** What is the minimum sample size required for exploratory factor analysis in order to reliably refine a multi-factor questionnaire?


title: Constructing a specific statistical hypothesis?

text: I'm going to pretty much rephrase my problem in a straightforward terms. Let's say one has $N$ friends and every one of them has either an odd or even number of cats with probability $0.5 \pm \epsilon$. If $kN$ ($k \in [0;1]$) different friends were randomly chosen and given a one additional cat then how likely is the fact that the probability of having odd\even number of cats is within $[0.5-\delta;0.5+\delta]$ range? How to evaluate this? I am really confused because this is essentially a probability of some other probability... Originally I wanted to evaluate how distorting can this effect of giving away cats be, but simply calculating the possible probabilities is not enough because e.g. if I had $10$ friends ($5$ of each category) and gave half of them an additional cat then the outcome can be everything from all friends with odd number of cats to all friends with even number of cats.


title: Replacing kernel for a Poisson distribution in an equation

text: In Example 4.4.2 on page 163 of Statistical Inference, the following is implied: $\sum_{t=0}^{\infty}\frac{((1 - p)\lambda)^t}{t!} = e^{(1-p)\lambda}$ with a note that says "[the] sum is a kernel for a Poisson distribution". I understand the Poisson distribution, but I don't fully grasp the kernel concept yet. I tried searching for this identity, but couldn't find it. Can someone help me understand this, or provide references that I can use for background information? Update: I initially had the incorrect equation above. I corrected it after Anand's comments.


title: Random variables mapping to a state space that isn't a set of real numbers?

text: A random variable usually has a definition like "a function that maps from the original sample space to a new sample space, typically a set of real numbers". Is it ever useful to map into a new space that isn't a set of real numbers?


title: Bandwidth value change each time the script in run in hdrcde function

text: The following code give different values of bandwidth everytime I run it. I use a Mackey-Glass time series 1500 data, double, negative and positive data.               library(hdrcde)     myfile = ("C:\\Users\\ernest\\Desktop\\A Rcode\\testfile.txt")     y <- read.table(myfile, header=FALSE,col.names=c("A"), dec=".")          HDRlevelVal <- 0.99      x <- y$A     hHDR <- hdrbw(x,HDRlevelVal, gridsize = 100)     HDRhat <- hdr.den(x,prob=100*(1-HDRlevelVal),h=hHDR)     print (hHDR)     print(HDRhat)      Any idea of what happen? Thank you.


title: Comparing classification performance over similar datasets

text: My original dataset (training set $n=291$) is multiclass and I evaluated the performance of 5 classifiers on a test set ($n=64$). I then converted the multiclass problem to a binary one and ran the same classifiers. The accuracy of SVM improved from 56.25% to 76.56%. I'm not completely sure which test to use for comparing the performance of SVM between the two datasets (the data are exactly the same except the classes are reduced from five to two). WEKA allows the use of a paired t-test; do you suggest anything different?


title: ARMA - GARCH fitting

text: I have an ARMA model which has these terms `Constant + AR23 + MA1 + MA3 + MA24`. The residual of this model has heteroskedasticity (from `ArchTest`) and I am trying to model these using a GARCH. While modelling the residuals through GARCH, I do want to retain the `AR23 + MA1 + MA3 + MA24` terms, although they can have different coefficients after optimization. Can you please let me know how to do this.   1. Which function (& package) shall I use to model GARCH if I have to retain all the above mentioned parameters?    2. I have tried following in `ugarchspec` function but it is not working.  2.1 When I use `fixed.pars` argument, I get this error: > Error in pars[zf, 1] = unlist(fixed.pars) : NAs are not allowed in > subscripted assignments function call as follows:              garch_spec <- ugarchspec(variance.model=list(model="sGARCH", garchOrder=c(1,1)),                               mean.model=list(armaOrder = c(1, 4), include.mean = TRUE),                               fixed.pars=list(ar23=-0.026922,ma1=-0.497002,ma3=-0.022282,                                               ma7=0.015574, ma24=-0.445093))      2.2 If don't use `fixed.pars`, then I get a different set of AR and MA terms which I don't want and get all nulls in robust standard errors function call:              garch_spec <- ugarchspec(variance.model=list(model="sGARCH", garchOrder=c(1,1)),                               mean.model=list(armaOrder = c(1, 4), include.mean=TRUE),                               start.pars=list(ar23=-0.026922,ma1=-0.497002,ma3=-0.022282,                                               ma7=0.015574,ma24=-0.445093))      2.3 In the mean.model, one can mention the order of terms but is there a way to mention the positions of terms, e.g., in my case `AR23 + MA1 + MA3 + MA24`? 2.4 Is there any other argument of `ugarchspec` that I need to use?


title: Posterior density

text: Absolutely stuck on section A of this question I found. I've been trying over 3 hours and still failing, can anyone here guide me through? Suppose a random sample $x = (x_1, x_2,\ldots, x_n)$ is taken from a normal N($\Theta$, 1) distribution. It is desired to estimate the mean $\Theta$. A normal distribution with zero mean and variance $\frac{1}{t^2}$ is used as a prior distribution for $\Theta$. (a) Show that the posterior density $\pi(\Theta|x)$ satisfies $$ \pi (\Theta |x)\propto\exp\left \\{ -\frac{1}{2}\left ( \Theta ^{2}\left ( n +t^{2} \right ) - 2n\mu \Theta \right ) \right \\} $$ where $$ \mu = \frac{1}{n}\sum_{i=1}^{n}x_{i} $$ By writing this posterior density in the form $$ \pi (\Theta |x)\propto exp\left ( -\frac{\left ( \Theta -m \right )^{2}}{2\nu ^{2}} \right ) $$ deduce the posterior distribution of $\Theta$. (b) Using your posterior distribution for $\Theta$, what is your estimate for $\Theta$? By recalling that $X_i - N(\Theta, 1)$, obtain the mean and variance of your estimate. (c) Discuss what happens to your estimate (i) if $n$ is large, (ii) if $t$ is large, (iii) if $t$ is small. (d) Discuss why someone might choose the case (i) large $t$, (ii) small $t$. $T_{post} = T_{prior} + T_{data}$ but I'm not sure how to utilize this. For Ci. is it correct that the estimate will become closer to the sample mean when $N$ is large? Cii. make $t$ converge towards the mean? Ciii. the larger the ratio it is towards the estimate?


title: Nonparametric test

text: To get the power of the t-test I use this R code:               pwr.t.test(n=25,d=0.35,sig.level=0.05,type="one.sample",alternative="greater")       and for the test of the Pearson coefficient of correlation I use this R code :               pwr.r.test(r=0.3,n=36,sig.level=0.05,alternative="greater")      I would like to know if it is possible to get the R code/script to get the power of the 3 nonparametric test (Wilcoxon, Wilcoxon/Mann-Whitney and the Kruskal-Wallis)?