title: How to include prior knowledge that a model might be able to figure out itself

text: I have a problem where I want to predict the outcome of a sequence given another sequence online. Let $(x_1, x_2, ... x_T)$ be denoted by $x_{1:T}$, then I am estimating: $$ p(y_T|x_{1:T}) $$ where $y_t \in Y, x_t \in X~\forall t$. That is, given a history of observations I want to predict the outcome of another variable. In practice, I predict all the $y_{1:T}$, but since I do it online I predict $y_t$ independently of $y_{t'}$ where $t' > t$. I do that with a discriminatively trained nonlinear state space model, a recurrent neural net which is trained with gradient methods based and the negative log likelihood as a loss. (This can be thought of as a hidden Markov model with dirac distributions everywhere.) What I am wondering is, how I can incorporate any knowledge I have about $p(y_{t+1}|y_t)$. **Example:** I am modelling the position of a car $y_t$ given a history of motor commands $x_{1:t}$. I know the maximum velocity $v$ of a car, thus I could assume that $$p(y_{t+1}|y_t) \propto \begin{cases} 0~\text{if}~y_{t+1} - y_t > v \\ 1~\text{else} \end{cases} $$ My question now is: how to make use of that? One attempt would be to penalize those hypotheses which produce output probabilities which violate $p(x_{t+1}|x_t)$. I tried this with adding a regularizer to the log likelihood during parameter estimation, but it did not work out pretty well. Another guess is to just multiply that prior with my output probability: $$ p(y_{t+1}|x_{1:t}) = p(y_{t+1}|y_t) p(y_{t+1}|x_{1:t+1}) p(y_t|x_{1:t}). $$ But this is wrong since $p(y_{t+1}|x_{1:t+1})$ and $p(y_t|x_{1:t})$ are not independent.


title: Is unpaired t-test appropriate for 2-level random cluster sampling? Also, how does this affect the power of the test?

text: I'm working on power/sample size calculations for a non-profit and want to make sure my logic is correct, as I haven't worked with a situation quite like this before! This is a bit of a long post, as I want to make the situation clear, but the first two questions are the important ones. **The Situation** There is a project with a large number of groups (100-700) of 20 people in each group. The people are all undergoing the same treatment and there is strong reason to believe that there is cluster level homogeneity (the groups are very similar to each other). Because of geography, manpower, and financial considerations, we're trying to arrange for a more cost effective way to sample the population and calculate the effect of the treatment. It's a continuous variable on a 0-100 scale with data collected baseline (before treatment) and endline (after). **Sampling** Based on this, I've recommended a cluster sampling technique. The groups aren't all starting at the same time, so we've set up a 2-level cluster sampling design with systemic random sampling for clusters (sample Cluster 2, 5, 8, etc. with random start point) and simple random sampling within the clusters (randomly select a certain number of people, say 5, within the selected clusters). We WILL be sampling the same clusters baseline and endline, but we WON'T necessarily be sampling the same people within the clusters (that part will be random again). **Proposed analysis** I'd like to do a two-sample t-test on this data, comparing the sample means (endline vs. baseline) for the variable of interest. With our desired MDES, power, and significant level, R tells me I need a sample size of around 150-175 for both baseline and endline. We're going to sample around 200-250 for both baseline and endline. **The questions** 1) Cluster sampling tends to decrease the power in your analysis, since you aren't randomly sampling from the entire population. Hopefully this effect is small, as we do believe we have cluster-level homogeneity. I've found this resource: http://www.amstat.org/sections/srms/proceedings/y2005/Files/JSM2005-000817.pdf It implies that we can calculate an effective sample size based on anticipated ICC and # of samples per cluster. Is this appropriate here? Is there a way to calculate this directly in R or SAS (other program examples are fine, but those are my go-tos)? 2) If the effective sample size is appropriate, how do we then calculate the actual statistics? Do we just run the same unpaired t-test, but with the smaller n plugged in? 3) Since we _are_ sampling from the same groups for baseline/endline, we could pair up the group means and run a paired t-test on that. However, since we're randomly sampling the groups and not sampling the whole group this would decrease the power. Again, any good resources for figuring out how it would be affected and actually calculating the power? 4) How about the use of blocking here? We have fairly obvious blocks (the same groups that we're using as clusters for cluster sampling). Since this is just one treatment and a before/after measurement, blocking through ANOVA would be the same as the paired t-test, yes? But since we can't pair up the individuals (at least not naturally), we would need to pair up the blocks, and then we're not in a blocking situation anymore. Thanks for the help!


title: Calculate Variance for Exponential Distribution

text: I have a problem with calculating the variance of an exponential distribution. In my formulary there are these formulas for exponential distributions: $E(X)=\frac{1}{\lambda}$ $V(X)=\frac{1}{\lambda^2}$ Where $E(X)$ is the expected value and $V(X)$ the variance. I have a distribution function on the form $F_X(x)=C_1(1-e^{-\lambda_1x})+C_2(1-e^{-\lambda_2x})$ where the $C$s and $\lambda$s are constants. I calculated $E(X)$ in the following way: $E(X)=C_1*\frac{1}{\lambda_1}+C_2*\frac{1}{\lambda_2}$ But when I try to calculate the variance using the formula above for V(X) and the same method as I used when calculating $E(X)$ I don't get the right answer. I get the right answer if I use $V(X)=E(X^2)-(E(X))^2$, but both should generate the same answer since $E(X^2)=\frac{2}{\lambda^2}$ (from partial integration) which gives $V(X)=\frac{2}{\lambda^2}-(\frac{1}{\lambda})^2=\frac{1}{\lambda^2}$ Can someone explain this please?


title: I need help about EFA before CFA

text: GOOD DAY I am PHD student in management. I need some help about my analysis. My model is relatively complex and I need your advises. My model aims to examine the relationship between knowledge m management and managers' competencies. Knowledge management is an independent variable and consists of six constructs to measure it (4 observed variables for each construct). managers' competencies are divided into two groups (management competencies and personal competencies). Management competencies include nine competencies (observed variables are ranged between 4 and six). Personal competencies include three competencies (13, 10, and 5 observed variables). I collected 227 responses and the correlation was good as hypothesized. Now I am going to conduct CFA to measure the model. my questions are: 1- Do I have to conduct EFA before CFA? (actually I like to do this to reduce the items so I reduce the complexity of the model but I couldn't find support from literature) - (when I did so, I deleted the items with low commonalities (under .6) I found that the measurement model worked very good, is this right?) 2- if it is better to conduct EFA, do I have to conduct EFA for the full model? or I can conduct it for the three groups separately? 2- Do I have to examine the full structural model or I can examine only the path model (composite model as in one of your tutorials)? (if yes, please I need support from the literature). Thank you


title: Combining sensitivity and specificity to measure classification performance

text: I have 2-classes labelled data on which I'm performing classification using multiple classifiers. And the datasets are well balanced. When assessing the classifiers' performance, I need to take into consideration how accurate the classifier is in determining not only the true positives, but the true negatives also. Therefore, if I use accuracy, and if the classifier is biased toward positives and classifies everything as positive, I will get around 50% accuracy, even though it failed at classifying any true negatives. This property is extended to precision and recall as they focus on one class only, and in turn to F1-score. (This is what I understand even from this paper for example "Beyond Accuracy, F-score and ROC: a Family of Discriminant Measures for Performance Evaluation"). Therefore, I can use sensitivity and specificity (TPR and TNR) to see how the classifier performed for each class, where I aim to maximize these values. My question is that I am looking for a measure that combines both these values into one meaningful measure. I looked into the measures provided in that paper, but I found it to be non-trivial. And based on my understanding I was wondering why can't we apply something like the F-score, but instead of using precision and recall I would use sensitivity and specificity? So the formula would be $$ \text{my Performance Measure} = \frac{2 * \text{sensitivity} * \text{specificity}}{\text{sensitivity} + \text{specificity}} $$ and my aim would be to maximize this value. I find it to be very representative. Is there a similar formula already? And would this make sense or is it even mathematically fine? Any hint or direction is appreciated.


title: Extreme value distribution for multivariate normal

text: I have a series of data sets. Each data set represents a measurement in 3D space relative to a global origin. I want to model the extreme values of my data. If I were to calculate the extreme radius and use that would I still expect a gumbel distribution? Alternatively I could produce a separate model for each axis if necessary but this is not as elegant as a single model. Ideally I would like to be able to predict the probability that a future experiment will have an extreme value less then X. Is there anyway I can determine the parameters of the extreme value distribution from my normally distributed data? I have done some testing and it appears that I could reasonably assume that all my data sets are sample from the same population so I could assume this. Or, is it necessary to perform an MLE fit. I appreciate this functionality is available in most analysis software I was just curious how much one can tell about the properties of the EVD from the mean and standard deviation from which the extreme values are drawn.


title: dealing with categorical variables when using softImpute for missing data

text: I'm trying to impute missing values but I have problem dealing with categorical variables. The command softImpute calculate the missing values but they also turn categorical variables, which is inadequate for the analysis. For the missing values I did the following               >softImp = softImpute(as.matrix(train), rank.max = 60)     >data.comp = softImp$u %*% diag(softImp$d) %*% t(softImp$v)     >data.comp=data.frame(data.comp)      However when I look at categorical variables they are in decimal points               > head(data.comp$X91)     [1] 0.6037109 0.6263665 0.5373208 0.6092270 0.8796817 0.8643236      which is originally               > head(train[c(91)])       H0001600     1        0     2        0     3        1     4        1     5        0     6        1      Is there a certain way to impute missing values for categorical variables? I Any other suggestions dealing with missing values in categorical variables would also help a lot.


title: Information entropy of a mixture distribution

text: I'm working with a high-dimensional mixture distribution, and I'm interested in calculating its entropy. I think I could work it out if there were only two mixture components. Following @Daniel's suggestion in the comments, here's what I have so far: > If the two components had equal weight, I think the answer would be > $\frac{1}{2}(H_1 + H_2 - I_{1,2})$, where $H_1$ and $H_2$ are the entropies > of the component distributions and $I_{1,2}$ is the mutual information > between them. My logic is that subtracting the mutual information prevents > double-counting any entropy. With a bit more work, I could probably figure out how to generalize it to unequal weights. But my problem involves thousands of mixture components, and I don't know how I'd go about accounting for all of them.


title: Out-of-time testing (basic question)

text: I understand the importance of out-of-sample testing, but could you tell me **why I should (or shouldn't) do out-of-time testing** ? The only use that comes to mind is if the predictive model applies to economic activity and seeing whether it would work in both bull and bear markets. But more insight in the use and importance(?) of out-of-time testing would be very welcome.


title: How do you measure the accuracy of an inference hypothesis/procedure?

text: Take inference to mean reasoning/predicting the value of a hidden/laten variable $Z$ given some evidence/data $X$. For example, maybe you are trying to find out if your patient has Cancer (Z = 1 if he does, 0 if he doesn't) but you only have indirect ways of measuring it X. Say X could be the color of his skin, the ph of his blood, or whatever important features (or kernels of features) you might want to use. One thing that I was curious is, if you normally do not have the label $Z$ in your data set, then how are you suppose to know how well your predictions for $Z$ might be once you create a hypothesis for predicting? Another way to think about this question could be, if you have a data set $S_n$ (Set of size n) and you want to find hidden structure in it, how do you know that the hidden labels ($Z$) that you are learning are good just from the set $S_n$? The only way I can think to asses this data set is, to somehow get "force" to get a "test set" that happens to have the hidden labels that are normally hidden and then predict see if the labels match with whatever cost function you want to you (say maybe counting number of correct predictions). However, if "forcing" to get $Z$ out of your data set is not possible, how are you suppose to get a feel of how good your hypothesis might be? Without labels, k-fold validation is not possible, right? Note that by force I mean is to obtain the unknown label $Z$ for some set that you plan to use as test set (though, "forcing" is probably not possible in many cases or tractable, since if it were tractable/doable, we probably would not be doing learning, we would just "force" the answer out of the data set every time we need to!) * * * I am both interested in a **"practical" answer** , what professional machine learners do in practice and a **theoretical answer** , i.e. an answer that relies more on mathematical techniques/proof on how to measure accuracy in this case or even, a proof of why maybe unless with further assumptions, it might be impossible to do anything.


title: How to calculate the product of Gaussian variables in Matlab

text: I want to calculate the distribution of a product of two i.i.d. Gaussian distributed variables a and b. In principle, this should be possible by defining a new variable x with a dirac delta distribution ![new variable x with dirac delta distribution](http://i.stack.imgur.com/6TO55.png) To get the distribution over x (the product of a and b), a and b have to be marginalized out: ![calculation of p\(x\)](http://i.stack.imgur.com/bi6Fz.png) Since the dirac delta ties the variables a and b together via the variable x=ab, we can skip one integral and get ![one integral skipped](http://i.stack.imgur.com/7GwxI.png) So far so good. If I now want to calculate this in Matlab, I get a huge difference between the solution with `integral` and a simple sampling solution as approximation. Here is my MWE:               mu_a = 2;     sig2_a = 1;          mu_b = 0;     sig2_b = 1;          res = 250;     xmin = -10;     xmax = 10;     x_range = linspace(xmin,xmax,res);          normal = @(x, mu, sig2) 1/sqrt(2*pi*sig2)*...         exp((-1/2)*((x-mu)^2/sig2));          normal_a = @(a) arrayfun(normal,a,mu_a*ones(size(a)),sig2_a*ones(size(a)));     normal_b = @(b) arrayfun(normal,b,mu_b*ones(size(b)),sig2_b*ones(size(b)));          product_pdf1 = zeros(size(x_range));     product_pdf2 = zeros(size(x_range));     for i = 1:length(x_range)             x = x_range(i);              normal_ab = @(a) normal_a(a).*normal_b(x./a);         product_pdf1(i) = integral(normal_ab,-Inf,Inf);              normal_ab = @(b) normal_a(x./b).*normal_b(b);         product_pdf2(i) = integral(normal_ab,-Inf,Inf);     end          subplot(2,1,1,'replace'); hold all     plot(x_range, normal_a(x_range))     plot(x_range, normal_b(x_range))     plot(x_range, product_pdf1, '--')     plot(x_range, product_pdf2, '--')     axis tight     drawnow          % the sampling solution     N = 50000;     X_sequence = zeros(1,N);     for i=1:N         a = mu_a + sqrt(sig2_a)*randn();         b = mu_b + sqrt(sig2_b)*randn();         x = a*b;         X_sequence(i) = x;     end          subplot(2,1,2,'replace'); hold all     hist(X_sequence,250)     xlim([xmin xmax])      Note that I calculated the integral once for a and once for b, the output is completely different. ![green and red lines: input Gaussians, dashed lines: attempt to calculate the product](http://i.stack.imgur.com/AAhbb.jpg) Depending on the parameters, one, the other or both of the calculations are clearly wrong. In my example, the sampling solution is much peakier than the better looking integral attempt. The other one is completely wrong. I now wonder: Is my math wrong? Is my code wrong? Or is it a limitation of Matlab's `integral` function (e.g. having trouble integrating over discontinuities)?


title: Generalizing the data per respondent to group

text: I am currently conducting a survey for a class requirement but I do not know how to process the raw data. I am using the "Media Fitness Framework" by Higa and Gu which basically quantifies a channel's ability to perform via a Likert scale. The article discusses how to find the most fit media per respondent, but doesn't tell us how to analyze all the data to get _the_ most fit media for the whole group I sorry I couldn't disucss their data analysis in more detail; I am not really good at statistics. Here is the link to the article if you find it helpful https://docs.google.com/a/dlsu.edu.ph/file/d/0B9lX0U_MS2D1dzU0UmtKSVdMcGM/edit


title: 1D weighted Isotonic Regression (PAV): a simple description of the algorithm

text: I have to use IsotonicRegression class from scikit-learn with non-uniform point weights: in method `IsotonicRegression.fit` parameter `sample_weight!=None`. I roughly know how the weighted PAV algorithm works, but I don't want to insert the rough description into my article. I just want to give some reference. But I could not find any suitable description of weighted PAV simple enough to be understandable by non-mathematicians, best of all — in the form of some pseudocode. Can someone help me?


title: Which statistic to prove that one data point dominates the proportions in dataset?

text: I'm very new to statistics so apologies if this is a dumb question. I've got a grasp of the concepts but not real-life usage of them (but I am learning!). I have a skewed distribution where one datapoint (google, see below) dominates the dataset (the Visits metric). I can intuitively see that the %CR for the entire data set (calculated as Total Transactions/Total Visits) is driven by this data point. But what is the best statistical approach to prove that? I've tried to phrase my question in different ways in Google to point me in the right direction but I've had no luck. ![dataset](http://i.stack.imgur.com/wgP7z.png)


title: Simple question about variance

text: I always thought I understood the concept of variance, but this one confuses me. $ \begin{eqnarray} X & := & (1, 2, 3)\\ E(X) & = & (1 + 2 + 3) / 3 = 2\\ Var(X) & = & Cov(X, X) = E((X - E(X))^2)\\ & = & E((-1, 0, 1)^2)\\ & = & E(1,0,1)\\ & = & 2/3 \end{eqnarray} $ But `R` tells me that               > var(c(1,2,3))     [1] 1      **Which part of my calculation is wrong?**


title: BIC vs. Out of sample performance

text: I have two statistical models. Model 1 uses a GLM approach while model 2 uses a time series approach for fitting. I want to compare these two models. Model 1 (i.e. GLM) has a better out of sample performance. Model 2 has a better BIC criteria. So based on out of sample performance, I should pick up model 1 and based on BIC I should pick up model 2 as the preferred model. I should add that in this context and for the question I am trying to answer, Both the BIC and out of sample performance are important. The question is how to choose the best model in this case? Should I consider other criteria? Please let me know if you know any good reference with similar cases.


title: R cv.glm returns NaN for stepwise-generated regression model

text: I'm trying to run K-fold cross-validation on a multiple regression model that was generated via the `step` function in R. However, the call to `cv.glm` returns `NaN`. Here is a simplified script that replicates the problem:               # this works just fine     library(boot)     data(mammals, package="MASS")     mammals.glm <- glm(log(brain) ~ log(body), data = mammals)     cv.glm(mammals, mammals.glm, K=5)$delta          # but when I use the step function...     mammals.simple <- lm(("brain ~ 1"), data=mammals)     mammals.glm <- step(mammals.simple, direction="forward", test="F", scope="brain ~ body")     # it returns NaN     cv.glm(mammals, mammals.glm, K=5)$delta      Is this due to some fundamental misunderstanding, where it doesn't make sense to cross-validate a stepwise-generated model? I'm using R v3.0.1 (but 3.1.0 has the same problem).


title: Coefficient p-values for linear regression on large dataset

text: Please forgive me if this is a stupid question; I am still kind of new to statistics. I am running an experiment. There are two features: Gender and Area. Gender has two levels (Male/Female) and Area has four (North/South/East/West). The dependent variable is height. My sample size is really big (about a million people). My goal is to find out whether/how much of the difference in height between men and women is due to the area they live in. Because the data groups are of different sizes, and I can't guarantee normality of distribution of homogeneity of variance, I didn't think I could use ANOVA. Instead, I was trying to use linear regression with dummy variables and check whether the interaction terms were significant. I got this idea by reading https://en.wikipedia.org/wiki/Dummy_variable_%28statistics%29#Interactions_among_dummy_variables. Is this the right approach, and if so, how do I find out whether the individual coefficients are significant if I have a really large data set? I have read about some examples involving popular statistics packages that gave a std. error, t-statistic and p value, but given how much data I am using I wasn't sure they would be able to handle it. I know there are faster ways to get (or at least to accurately estimate i.e. stochastic gradient descent) just the coefficients, but I don't know about the standard error.


title: How to bootstrap respecting subject-level information?

text: This is my first post and I am not a skilled programmer, so please let me know if the question or code are unclear. I am trying to bootstrap an interaction (that is my test statistic) using the package "boot" in R. My problem is that for every resample, I would like the randomization to be done within subjects, so that observations from different subjects are not mixed. Here is the code to generate a dataframe similar to mine:               Subject = rep(c("S1","S2","S3","S4"),4)     Num     = rep(c("singular","plural"),8)      Gram    = rep(c("gram","gram","ungram","ungram"),4)     RT      = c(657,775,678,895,887,235,645,916,930,768,890,1016,590,978,450,920)     data    = data.frame(Subject,Num,Gram,RT)       This is the code I used to get the empirical interaction value:               summary(lm(RT ~ Num*Gram, data=data))      As you can see, the interaction between my two factors is -348. I want to get a bootstrap confidence interval for this statistic, which I can generate using the "boot" package:               #Function to create the statistic to be boostrapped     boot.huber <- function(data, indices) {      data <- data[indices, ] #select obs. in bootstrap sample      mod <- lm(RT ~ Num*Gram, data=data)      coefficients(mod)       #return coefficient vector     }          #Generate bootstrap estimate     data.boot <- boot(data, boot.huber, 1999)          #Get confidence interval     boot.ci(data.boot, index=4, type=c("norm", "perc", "bca"),conf=0.95) #4 gets the CI for the interaction      My problem is that I think the resamples should be generated without mixing the individual subjects observations: that is, to generate the new resamples, the observations from subject 1 (S1) should be shuffled within subject 1, not mixing them with the observations from subjects 2, etc... I don't know how "boot" is doing the resampling (I read the documentation but don't understand how the function is doing it) Does anyone know how I could make sure that the resampling procedure used by "boot" respects the subject level information? Thanks a lot for your help/advice!


title: How to determine a statistical model to fit into data based on distribution

text: I am trying to learn to fit data into models.The data is multi dimensional with `n` observations for each of the `k`dimensions.I went through the link Link but in my case the data is multi dimensional having several columns of data. The data is obtained from each temporal frame of a video representing features such as pixel coordinate (x,y) ,speed,acceleration along x and y component for 20 videos each of 200 frame duration.Can somebody please explain how to do this and whether it is possible to do linear or nonlinear regression models.   1. Is it necessary to find the distribution of the data before fitting into a model? Is there a criteria for selecting a model based on distribution?   2. Should the number of variables in the model equal to the number of columns of data?If so then it will not be an AR model, I think. I should then go for higher dimensional polynomials where each variable will represent a feature vector.Or are there other linear or nonlinear regression models?I am particulary interested in NARMAX modelling but do not understand the complexities.   3. After fitting into the model,how do I test the "goodness of fit" ? An illustration with an example of a model will be really helpful.


title: Classifying data in R

text: I'm using R and I want to use some type of machine learning technique (SVM, perhaps) to classify data. I have data that looks like this, for example                    [,1] [,2] [,3] [,4] [,5]     [1,]    3    4    5    6    0     [2,]    2    6    7    3    1     [3,]    7    1    4    6    1     [4,]    3    4    5    9    0      Where columns 1, 2, 3, and 4 are variables with some provided data, and column 5 is the classification. So, data entry 1 has values 3, 4, 5, and 6 for the four variables, and is classified as a 0. Classifications can be either 0 or 1. I want to be able to build an SVM so that given a new set of variables, it will classify it as "0" or "1" What is the best way to do this?


title: Errors in Variables and Deming's multivariate regression: Assumptions

text: There has been extensive literature that puts forth a standard set of assumptions for the Ordinary Least Squares (OLS) estimator. I am very interested in working around the two classical problems of a)Errors in Variables regression, b)Deming's Regression in a multivariate setting where the error's in measurement is only in the independent variables and not in the response. My problem is that I have not been able to find any literature that gives a standard set of assumptions for these two problems in a multivariate setting within the literature. Can someone put forth the assumed models and the distributional/probabilistic/independence assumptions clearly and succinctly for one or both of these problems in a multivariate setting with a matrix notation? A majority of my problem's/doubts are usually resolved once the assumptions are put forth clearly.


title: Probability model question

text: Counts of the number of broken bones in college athletes during the season would be best represented by which of the following probability models? Question options:   1. Binomial --Thinking this is the best answer since they either yes broke a bone or no did not break a bone   2. Poisson    3. Uniform    4. Normal (This is not the correct answer)


title: Power of chisquared test for large sample size

text: I would like to calculate the power of a Chi-Square test for goodness-of-fit as a function of sample size for a specified alpha-value (say 0.01). Specifically, I am referring to power as the probability that the test will correctly reject the null hypothesis. In most places I look, I can find only vague references that this can be done and/or am referred to software. It seems to me that this should not be difficult. The specific question I would like to be able to answer (without using some sort of opaque software) follows: Given an α value, for what sample size, n, can I expect a particular power, β? Even better, how might I be able to calculate an ROC?


title: Binary classifier testing with minority positive data

text: I have a question about a testing methodology for a binary classifier, that I'm not entirely sure how to describe. So apologies if this has been answered many times before but I haven't had much luck searching for it. I'm attempting to test a binary classifier where the positive class is a tiny minority of the data, approximately 1% of the actual data. My question is: what is the best way of testing a classifier so that a failure to classify something as positive actually counts towards the accuracy score. For instance, if the classifier always said Negative, it would still be 99% accurate. I've been trying 5-fold cross-validation but obviously the above case occurs and it always give an amazing score. Any suggestions would be greatly appreciated.


title: Are Bayesian approches used for classification (supervised) or for clustering (unsupervised)?

text: Are Bayesian approaches (static and dynamic) used for classification (which is supervised) or for clustering (which is unsupervised)? or can they be used for both ? I even see that for instance to compute the likelihood they need the class labels of data, so I was thinking that it is only convenient for supervised cases where we have the class labels


title: Equation for a logit link function for a series of events

text: I have modeled some data using generalized linear modeling with a binomial distribution and logit link function. However, my data is not dichotomous, it is actually a series of events. I have fixed the number of trials at 20 (because there were 20 trials) so it is not a binary outcome (i.e., logistic regression). I believe the logit link function looks like this: exp(X * beta)/[1 + exp(X * beta)] or In(natural log) (probability of the event happening/ 1 - the probability of the event not happening). I am trying to do something on excel so I need the actual equation for the link function I used but because this is not a traditional logit I am a bit confused. Can someone tell me based on this what the equation would be in a form I could use? In other words, what does the logit link function equation look like when instead of a binary DV it is viewed as some number of events, for a fixed amount of trials (20 in my case)? Thank you!


title: What happens to an (AR)MA model when doing out-of-sample forecasting?

text: What happens to the error terms in an (AR)MA model when doing out-of-sample forecasting? As I understand it, when doing an in-sample fit, the estimate is simply the residual of the ground truth data. However, in an out-of-sample setting I don't know the residual. What happens? Do I take zero as the expected value as an estimate (which would make a MA model produce constant forecasts if I understand correctly)? Could somebody please confirm or correct this? Thanks!


title: Stationary processes I(0) cointegration, does it make sense?

text: I'm looking at two time series and I would like to determine how they move together. They both are however stationary. Would it make sense to test for cointegration? Linear relationship between both series has a very low R^2 and correlation but I think that there might be a lag or something and that they move together.


title: Conditional expectations by conditioning on functions of random variables

text: I have conjectured the following: Let $f:\mathbb{R}\supseteq A \rightarrow B \subseteq \mathbb{R}$ be an injective function. Let $X$ be a random variable with support $A$ and $Y$ be some random variable that is not independent from $X$. Then, $$E[Y | X=x]=E[Y|f(X)=f(x)].$$ Is that correct? If it is correct, is there any "weaker" assumption (weaker than $f$ being injective) that would make this true? Thanks.


title: Asymptotic distribution of maximum order statistic of IID random normals

text: Is there a nice limiting distribution of $\max( X_1,X_2,...,X_n) $ as $n$ goes to $\infty$, assuming that they are _iid_ normal distributions with variance $\sigma^2$. This is almost certainly a well known problem with a clever proof and nice solution, but I've been digging around and haven't found anything.


title: Optimization of models (ANN, radial basis, etc.) in R to target predictor levels to produce a desired response

text: I've learned how to use some packages like `neuralnet` and `caret` to create models based on designed experiments and am getting to the point where I think my models are relatively decent. There's a program called ISight that some others at my company use for radial basis function analysis. I prefer R for a bunch of reasons, but one of the features it has is an "optimizer." This lets you input a desired value for your response(s) and it chugs away and produces it's best suggestion of what input variable values will achieve this output. It also allows you to weight each output response importance (if you have multiple responses) and will drive toward a weighted/balanced response set based on your ranking. I'm wondering if there's something like this for R? I could just generate a huge data frame of various combinations, run it through the `predict` functions of my models and then look for the responses around my target value, which will tell me at what approximate input values I can achieve that value? That will probably be fine for one response, but when I start trying to model multiple responses it will be much more difficult to figure out which input levels optimize across the whole array of desired output values. Feel free to provide any terminology used to describe this. I'm new to the area and only know that ISight calls it the "optimizer." My google searches may have been poor due to not knowing what to call what I'm looking for.


title: Specifying the linear predictor, survival function and PDF of a log-normal survival distribution

text: Continuing on my exploration of the log-normal distribution, I'm working on reimplementing some code originally written for a Weibull/Exponential model for a log-normal model. Among the things it does is use SAS's PROX NLMIXED to hard-calculate some likelihoods. While this is technically SAS code, I think it's pretty readable in a code-agnostic fashion. For example, the Weibull Version:               *Linear Predictor of Event A     lam=exp(-(f0*alpha+f1*alpha*X));          *Density of Event A;     ff1=alpha*lam*t1**(alpha-1)*exp(-lam*t1**alpha);     ff2=alpha*lam*t2**(alpha-1)*exp(-lam*t2**alpha);          *Survival function of Event A;     sf1=exp(-lam*t1**alpha);     sf2=exp(-lam*t2**alpha);      The two survival functions are just at two different points in time for interval censoring. I've implemented, I believe, the appropriate survival and density functions for a Log-Normal as follows:               *Density of EventA;     ff1 = exp(-0.5*((log(t1)-mu)/sigma)**2)/((t1*(2*CONSTANT('PI'))**0.5)*sigma);     ff2 = exp(-0.5*((log(t2)-mu)/sigma)**2)/((t2*(2*CONSTANT('PI'))**0.5)*sigma);          *Survival function of EventA;     sf1= 1 - CDF('Normal',((log(t1)-mu)/sigma));     sf2= 1 - CDF('Normal',((log(t2)-mu)/sigma));      I guess the first part of the question is does that appear correct? My second question is where I'm getting caught up, actually defining the linear predictor for the log normal. I think it should be something along the form of: `mu = exp(f0+f1*X)` ala the Weibull above and with sigma estimated outside the linear predictor, but I'm not positive. Can someone point me in the right direction? Edit: The goal PDF once the syntax errors are fixed is: ![enter image description here](http://i.stack.imgur.com/us46d.gif)


title: How to generate subsets with the highest T-value

text: I following this  data mining project. It explains how to do it but I don't understand the T-value part. **"For each class, generate subsets with top 2,4,6,8,10,12,15,20,25, and 30 top genes with the highest T-value"** So I calculate the T-value like this: $$(Avg1 - Avg2) / \sqrt(Stdev1^2/N1+ Stdev2^2/N2)$$   * $Avg1$ is the average for one class across the gene sample and $Avg2$ is the average for the other 4 classes.   * $Stdev1$ is the standard deviation for one class and $Stdev2$ is the standard deviation for the other classes.   * $N1$ is the number of samples that have the class whose $T-value$ we are interested in, and $N2$ is the number of samples that does not have the $T-value$ that we are interested in. * * * How do I use that T-value to create those subsets? ![The Data](http://i.stack.imgur.com/vhwEX.jpg) The bold values in the data represent Average and Stdev per column.


title: Philosophical justification for possible but "implausible" (probability = 0) events in Bayesian Inference?

text: It violates some intuition in me that no matter what particular outcome has occurred of a variable described by a continuous distribution, it must have had probability 0. In Frequentist inference, I feel ok (but not great) by recognizing that the asymptotic value of the number of occurrences of any outcome divided by n (number of trials) is 0 as n approaches infinity. In Bayesian, I don't feel satisfied (yet) with just saying that you can (correctly) have had 0 expectation of certainty for an outcome that has indeed happened. I'm just hoping to find a perspective to understand this better. Thanks! Edit: Due to some confusion in the comments, I just want to emphasize that this is not a question about why point events must always have probability 0 in a (pure) continuous distribution. That is a strictly mathematical consequence because, for a given probability space, there cannot exist an uncountable number of disjoint events with probability greater than 0. I'm just having trouble reconciling how we can interpret that with a Bayesian approach.


title: Out-of-core data analysis options

text: I have been using SAS professionally for close to 5 years now. I have it installed on my laptop and frequently have to analyze datasets with 1,000-2,000 variables and hundreds of thousands of observations. I have been looking for alternatives to SAS that allow me to conduct analyses on similar sized data sets. I am curious what other people use for situations such as this. This certainly isn't "Big Data" in the manner that is used today. Nor are my datasets small enough to hold in memory. I need a solution that can apply algorithms to data stored on a hard drive. These are the things I have investigated to no avail:   1. R - BigMemory can create matrices stored out of memory, but the elements have to be the same mode. I work with data that is almost a 50/50 split between character and numeric. The FF package gets closer to what I need, but I don't quite understand which procedures are compatible with it. I think support is somewhat limited.   2. Pandas - I was very excited about a Pythonic alternative to R. However, it too has to hold all of the data in memory.   3. Revolution R - This one shows quite a bit of promise. I have a copy on my home computer (free if you sign up for Kaggle) and have yet to test it as viable alternative to SAS. Comments on Revolution R as a SAS alternative are much appreciated. Thanks **UPDATE 1** Editing to add that I am looking for real-life, practical solutions that people have used successfully. For the most part, SAS lets me chug through big files without worrying one bit about memory constraints. However SAS is implemented, they figured out how to make memory management transparent to the user. But, it is with a heavy heart that I used SAS for my job (I have to) and would LOVE a FOSS alternative that allows me to work on "large" data without having to think too hard about where the data is located at a specific time (in memory or on disk). The closest things I have come across are R's FF package and something on the horizon for Python called Blaze. And yet, these problems have existed for many years so what have analysts been doing in the mean time? How are they handling these same issues with Memory limits? The majority of solutions on offer seem to be:   * Get more RAM -- This isn't a good solution, imo. It's easy to find a dataset that can exceed RAM yet still fit on a hard-drive. Furthermore, the work flow has to accommodate all of the structures that are created during exploratory data analysis.   * Subset the data -- This is fine for exploration but not for finalizing results and reporting. Eventually, whatever processes are developed on a subset will have to be applied to the entire dataset (in my case, anyway).   * Chunk through the data -- This is what I would like to know more about from people who actually implement this work-flow. How is it done? With what tools? Can it be done in a way that's transparent to the user? (i.e., create some on-disk data structure and the frame-work takes care of the chunking under the hood).


title: Intuitive feel in mixed effect models

text: This is something I have been thinking about for sometime. Consider the random effects model $y= Zu + e$ where $u \sim N(0, \sigma^{2}I)$ and $ e \sim N(0, \epsilon^{2}I)$ I now want to compare 2 models. In the first model I use an overfit matrix $Z_{1}$ in which one singular value is large and the rest are close to zero. In the second model I use a matrix $Z_{2}$ which is a rank-1 approximation of $Z_{1}$ (i.e. de-noised). Intutively speaking, how would we expect the estimate of the means and standard errors of Variance($Z_{1}u_{1} + e_{1}$) and Variance($Z_{2}u_{2} + e_{2}$) to compare. Is there some paper which makes this hand wavy argument concrete? Although this is a little far fetched, I have been thinking if the estimates of $u$ show something equivalent to the bias-variance tradeoff seen in linear regression


title: Confusion regarding Likelihood Ratio Tests (LRT)

text: I am trying to construct a LRT to test the hypothesis $H_0: p \ge p_0$ and $H_1: p < p_0$ where $\alpha = .1$ and $p_0 = .6$ and give a critical region. Attempt: $\lambda(x) = \frac{"restricted" MLE}{"unrestricted" MLE} = \frac{L(\hat{\theta_0}|X)}{L(\hat{\theta}|X)}$ I am looking for $P(X \in \Re) \le \alpha$ and $P(\lambda(x) \le c) \le \alpha$ where $C \in (0,1)$ $Reject: \{ x: \lambda(x) \le c \}$ My confusion, is with this little information, how can I calculate $\frac{L(\hat{\theta_0}|X)}{L(\hat{\theta}|X)}$.


title: Literature for prediction models where each training example has a different amount of data?

text: This could be a machine learning question as much as a statistics question, but I think this is the best place to put the question. Here are three different examples of problems where each observation can have a different amount of data:   * Users surfing a website will click around from page to page. Based on the pages the user has visited, how likely is she to leave the website?   * Customers shopping in a supermarket move from aisle to aisle filling their cart. Based on the products they have stopped to look at (and their order), how much longer do we expect them to continue shopping for?   * Baseball players often spend time in the minor leagues before making their major league debut. Some players will be in the minors for 2-3 years, while others can have 5-10 years of minor league experience. Based on all of a player's minor league stats, how can we predict his rookie statistics if he were to jump to the majors? In each case, we'll have some examples that consist of 1 datapoint (a user only went to one site and left, a customer picked up bread and went straight to the register, a baseball player went straight to the majors after one minor league season), while other examples can have 2 or 3 or 5 or 10 individual datapoints. What are some literature examples, or keywords that I can search for, that explains how to model this type of problem?


title: Clear steps to calculate coherence between two time series

text: _I originally posted this on stackoverflow.com and then deleted it and moved it here_ My question is similar to Similarity of two discrete fourier tranforms (specifically the selected answer). I've also gleaned some helpful information from this R help thread. I'm stuck on the actual implementation of this process, however. I have two Fourier series; my ultimate goal is to calculate a sum of weighted coherence. (In actuality I have millions of series corresponding to pixels, but, let's just call it two for now).                   TS1 <- c(5051.29, 5355.31, 5602.18, 5784.4, 5896.45, 5934.9, 5898.6,      5788.64, 5608.37, 5363.27, 5060.78, 4710.09, 4321.86, 3907.89,      3480.75, 3053.42, 2638.89, 2249.75, 1897.82, 1593.81, 1346.93,      1164.71, 1052.67, 1014.21, 1050.52, 1160.47, 1340.74, 1585.84,      1888.33, 2239.02, 2627.25, 3041.22, 3468.36, 3895.69, 4310.22,      4699.36, 5051.29)         TS2 <- c(4192.83, 4532.62, 4836.41, 5094.96, 5300.41, 5446.53, 5528.88,      5544.95, 5494.25, 5378.33, 5200.71, 4966.78, 4683.65, 4359.93,      4005.45, 3630.98, 3247.9, 2867.85, 2502.38, 2162.59, 1858.8,      1600.25, 1394.8, 1248.67, 1166.33, 1150.26, 1200.95, 1316.87,      1494.5, 1728.43, 2011.55, 2335.28, 2689.76, 3064.23, 3447.31,      3827.36, 4192.83)      From the above links I see that I need to calculate coherence of the cross- spectrum at each frequency, and weight the frequencies that have higher power spectral density. I know I can extract spectral density from               spectrum(data.frame(TS1, TS2))$spec      I know I can extract coherence using:               spectrum(data.frame(TS1, TS2))$coh      I can see that the value of coherence extracted is always equal to 1. It seems the "solution" to this is to lag the time series somehow. I'm not sure how to 1) calculate optimal lag, especially because I'd like it to be comparable across all the coherence calculations, and 2) how to specify lag in the spectrum(). spec.pgram takes "spans", but those seem to have more to do with smoothing than offsetting. I tried a maximum ccf function defined here, but in many cases (not in the example set I've provided, you'll have to trust me) the "optimal lag" returned is 0, which makes sense because the ccf is 1, but isn't actually what I'm looking for. Should I be using `acf()` or `ccf()` to manually calculate coherence? What steps (R functions would be really helpful) do I need to perform to get from where I am to having a single weighted coherence value?


title: Bootstrap residuals: Wild vs Block Bootstrap?

text: (This is related to my (so far unanswered) last question) I want to use residual bootstrap to examine uncertainty and robustness of a model that fits a series of environmental measurements. I know close to nothing about the structure of my data's noise. It appears to be heteroskedastic and may (!) be autocorrelated in some cases (I have a few different data sets). Would the most simple form of the wild bootstrap (switching signs of the residuals with 50 % probability) be valid in this context, or would I have to rely on harder to implement methods like (moving) block bootstrapping because of possible autocorrelation? Please stop me if I'm completely off the track; while the idea of bootstrapping is more or less easy to understand, the many available variants are extremely confusing to a newcomer. /edit: Just to be clear: I want to resample the residuals (as an estimate for the error), not the time-series itself.


title: Statistical formula for intra-rater reliability

text: I designed an assessment rubric for a class assignment and wish to calculate my own intra-rater reliability. My goal is to use the rubric to evaluate the same ordinal data at two different points in time. Would a Cohen's kappa be fine for this or is there another formula in SPSS you would recommend? Thank you.


title: How to calculate confidence intervals using subsampling after a nonparametric estimator about the empirical distribution function?

text: I have a problem where I think subsampling is more appropriate than the bootstrap. (Reason in another post.) However, I found no quick reference on subsampling CIs, and my naive inversion of the basic theory led to something that looks like the analogue of the "basic bootstrap" CI which seems to be used less often itself. Also, another question got a reply that subsampling applications usually "rely on asymptotic normality of the estimate", and "a bootstrap of some sort is more appropriate eg for asymmetric confidence intervals." My (limited) understanding of subsampling does not seem to imply that. I think I can get CIs as $$(\theta -\frac{\sqrt{b}}{\sqrt{n}}(\theta^{*}_{(1-\alpha)}-\theta);\theta -\frac{\sqrt{b}}{\sqrt{n}}(\theta^{*}_{(\alpha)}-\theta))$$ where $\theta^{*}_{(1-\alpha)}$ denotes the $1-\alpha$ percentile of the subsampled coefficients $\theta^{*}$ and $\theta$ is original estimate (so yes, everything could use $\hat{\theta}$ here). So, if I have a nonparametric estimator for the "excess mass" at some parts of the distribution (implicitly using kernel density estimates usually unfit for the bootstrap), and I "simply" repeat the process with subsamples without replacement, how do I get a CI from the distribution of estimates? Take the raw percentiles (á la "percentile bootstrap") or invert the statistic from the general theory ("basic bootstrap"). Monte Carlos would be prohibitively expensive to compare these. By the way, I know subsampling is not really a smart move without deriving something about the rate of convergence (simple square root above) and the right method to pick the size of the subsample (b). But is a "percentile bootstrap" any smarter in this case? Finally, maybe a question worthy of its own posting, but if "b=n" is not required for the common bootstrap either (nor is no-replacement, thus subsampling could be a subset of bootstrap methods), why isn't the $\frac{\sqrt{b}}{\sqrt{n}}$ correction part of the bootstrap CI formulas (basic or BCa)?


title: Outlier detection/imputation - discussion

text: _**Introduction_** : I'm working with heart rate data. IBI (InterBeat Interval) is defined as the time period between any two consecutive heart beat and is usually measured in millisecond. I have followed a subject for 6 days and using a device, I have all his IBI measures. My dataset has about 450,000 IBI measure. The device I'm using to measure IBI's is not 100% accurate. Therefore, I may get IBI measures that our out of range and are not consistent with previous and future measures. Consider the data of: `1500,580, 590, 570, 580, 1450, 560, 590, ...` for sure, 1450 is an outlier. _**Question_** : I'm trying to find a reasonable way to:   1. Detect outliers (like 1500, 1450 in the example above)   2. Impute a reasonable value for detected outliers _**What I have done so far_** : My idea is to fit a cubic natural spline to the data. Then, any point outside of the 95% confidence interval around the N.S. fit, is considered as outlier. I can then impute for the outliers by using the fitted value from the N.S. fit. _**Issue_** : Consider the data above, the outliers `1500` and `1450` have substantial effect on the N.S fit. They might even cause the true values right next to them (580, 590, ...) being marked as outlier. This means that even those values are valid, they need to be imputed with fitted value which in case are lot larger than those valid values. _**How can I solve the issue above?_** Is it reasonable to first fit an initial N.S. fit and mark all outliers. Then, fit a secondary N.S. (natural spline) fit with all outliers excluded from the dataset and use this secondary fit to impute for outliers? Would you have any better idea? Thanks a lot for your help. I believe this message can initiate a very good discussion on outlier detection/imputation methods.


title: Partitional vs. hierarchical clustering with distance matrix

text: I am exploring the flexibility of partitional clustering algorithms. In particular, I would like to introduce more general distances than the ones which are used by default. Let us consider, for simplicity, `dbscan` contained in the R-package `fpc`. It allows the user to specify a "data matrix, data.frame, dissimilarity matrix or dist-object". My idea would be to compute the distance matrix of the given data w.r.t. my chosen distance, and run `dbscan`. Here comes the point where I am stuck. Is it true that specifying a distance matrix should lead inevitably to a hierarchical clustering? My intuition says that a hierarchical clustering in presence of a distance matrix makes more sense that a partitional one on the elements of the matrix itself. As I am no expert in clustering, I cannot judge the above statement properly. Would you use a partitional algorithm on the distance matrix of a given dataset? Is this correct? Thank you all!


title: Confusion in interpreting precision matrix

text: I have this confusion related to interpreting the precision matrix. Suppose I have four variables forming a multivariate gaussian distribution with the following precision matrix               Q =  1 -1  0   0         -1  2 -1   0          0 -1  2  -1          0  0 -1   1      If I invert this precision matrix and get the covariance matrix from there if I calculate the correlation matrix, I get                  R = 1    0.998006970135374   0.996015933289707   0.994030843789017     0.998006970135374   1   0.998004987033912   0.996015933289707     0.996015933289707   0.998004987033912   1   0.998006970135374     0.994030843789017   0.996015933289707   0.998006970135374   1      So even though the variables are conditionally independent, they are highly correlated. Can anyone explain the intuition behind this? I mean they are conditionally independent since there are 0 entries. But still they are highly correlated.


title: How to compute the prediction variance from Relevance Vector Machine regression

text: I am using a relevance vector machine in R, rvm(), to solve a regression problem. I need to know the variance of the fitted values for each identified RVs. Does anyone know how to program to compute the prediction variances with R's Kernlab package?


title: How is the sample size determined for a test?

text: In general, how is the required sample size derived for a specific power and a specific significance in a test rule?   1. Since the alternative hypothesis may not be simple and the power is a function of the sample's distribution, what does the specific power require? Is it the minimum power over the sample's distributions in the alternative hypothesis?   2. If the power function is continuous at the border between null and alternative hypotheses, no matter how many the sample size is, the minimum power over the sample's distributions in the alternative hypothesis is always the significance. So isn't it meaningless to find a sample size to achieve a specific power? For example, for a t test for $\mu=0$ versus $\mu \neq 0$ with significance $\alpha = 0.05$, the power function is a continuous function of $\mu$ at $\mu=0$. Thanks and regards!


title: Hausman test: Include or not year effects and/or interaction variables

text: I have a problem when performing a Hausman test. I have a panel dataset that has five panels. I am estimating the same model twice, once using quarterly and another using half-year data. My dependent variable and some of my explanatory variables contain time-series data, which changes across individuals and time. However, I also have one time dummy, year, and two interaction terms with a dummy. The year dummies do not change across panels, since I consider the same years for each panel, and do not systematically change along time, since in the case of quarterly data I have the same year 4 times per panel, and for half-years I have year twice per panel. Something like this               Year | Quarters          1998   1998q1     1998   1998q2     1998   1998q3     1998   1998q4     (...)      In the case of the interaction terms, the values do change across individuals and time, but only in 4 of the 5 panels, since for the first panel the variable is multiplied by zero. Therefore, all the values for the first panel are zero. According to Wooldridge (2010, p.329) "Econometric analysis of cross sectional and panel data", in the section about comparing FE and RE, he says: > "Because the FE approach only identifies coefficients on time-varying > explanatory variables, we clearly cannot compare FE and RE coefficients on > time-constant variables. But there is a more subtle issue: we cannot include > in our comparison coefficients on aggregate time-effects--that is, variables > that change only across t. (...) the problem with comparing coefficients on > aggregate time effects is not one of identification; we know RE and FE both > allow inclusion of a full set of time period dummies. The problem is one of > singularity in the asymptotic variance matrix of the difference between FE > beta estimate and RE beta estimate." After experimenting I have the following problems: 1) If I regress only using the 'pure' variables (no interaction), with/without year effects I get the error I asked about here. 2) If I include the interaction terms, everything seems OK. But, **is it OK including these interaction terms when at least in one panel its values do not change along t?** 3) The result of the test from including/not including year effects are different, in the sense that in one case it's significant and in the other it's not. Independent of these results, **should I include year effects (year dummies) in the model from which I get the estimates I use for the Hausman test?**


title: Markov models with conditional transition probabilities

text: First, let me acknowledge up front that I'm not as well versed in statistics and mathematics as I'd like to be. Some might say have just enough knowledge to be dangerous. :D I apologize if I'm not using terminology correctly. I am trying to model the probabilities of a system transitioning from one state to another. A simple Markov model is a good start. (Set of states, set of initial state probabilities, set of transition probabilities between states.) However, the system I'm modeling is more complex than that. The transition probabilities leading to a state at time T are most certainly dependent on variables other than the state at T-1. For example, S1 -> S2 might have a transition probability of 40% when the sun is shining, but S1 -> S2 probability goes to 80% when it is raining. **Additional info from commenters' questions:**   1. The states are observable.   2. There will only be 5-10 states.   3. There are currently about 30 covariates that we want to investigate, though the final model will certainly have fewer than this.   4. Some covariates are continuous, others are discrete. **Three questions:**   1. How can I incorporate conditional transition probabilities into my Markov model?   2. Or, is there another perspective entirely from which I should approach this issue?   3. Also, what keywords/concepts should I search for online to learn more about this? I've already been around the web searching for things like "markov models with conditional transition probabilities," but so far nothing has slapped me in the face and said, "This is your answer, dummy!" Thank you for your help and patience.


title: Statistical significance of cluster validity

text: Hi I m working on a unsupervised problem to partition my dataset. I have access to the class labels for this dataset. Now I am trying to use Jaccard coefficient to compute correlation between cluster and class matrices of my data to perform cluster validation. I have two questions: a) Is Jaccard coefficient (JC) better than Random Statistic based cluster validation since the class and cluster matrices are symmetric.   b) Statistical significance of JC measure obtained. I studied this paper "On Clustering Validation Techniques" Journal of Intelligent Information Systems, 17:2/3, 107–145, 2001, where they talk about clustering validation techniques and significance testing approach using an example on page 21. I am unable to understand what they are trying to say when they say the following: rule _*_ "We accept or reject the null hypothesis whether (1−ρ)*r = (1−0.05)*100 = 95 values of Ri , Ji, FMi , i are greater or smaller than the corresponding values of R, J , FM." Here ρ=0.05 and r=100 different random datasets generated. In a nutshell,first they are generating what is the similarity based metric for this cluster validation using Jaccard coefficient, Random stat and others. Next they are generating 100 random datasets under uniform distribution to test the significance of the Jaccard measure, Random stat and others. Then,they are applying same clustering parameters on all the random datasets individually and computing the validity based measures. Finally, they are comparing the random dataset measures against the initial cluster validity using the above rule _*_. I am having trouble understanding that rule. For example, lets say my Jaccard coefficient was 0.68, are they saying if 95 random Jaccard coefficient values are >=0.68 then they reject the null hypothesis. But to me that's a huge error acceptance. Please comment. Any help is appreciated.


title: Visualizing many p-values

text: I have a dataset, on which I am doing some data mining. My mining method has certain parameters which I can fiddle with, presumably leading to better or worse quality results. I measure the quality of these results with some statistical tests, which gives me a low P-value for "good" results. This generates a lot of P-values over a large range. I want to visualize how the P-value changes in response to my parameter. Currently I do this by plotting the P-value over the parameter, however the plot doesn't come out smooth at all. For my personal use, I tried "smoothing" by doing a moving average of the P-values. I know this is a nonsense method but it helps me see the general trend. I have two questions:   1. What is the "right" way to smooth such a curve, ie. a curve generated by a function `f(x)=p` which returns a P-value corresponding to each parameter x?   2. Does the answer to 1 change if I was plotting the log of P-values, rather than the values themselves?


title: Cyclostationary time series

text: http://en.wikipedia.org/wiki/Cyclostationary_process What are the methods in modelling and forecasting such time series? It is mentioned in the link above that there is a deterministic approach to treat such a time series. Can anyone explain what does it mean?


title: Missing rates and multiple imputation

text: Is there a limit which is the least acceptable when using multiple imputation (MI)? For example can I use MI if the missing values in a variable are the 20% of the cases while and other variables have missing values but not to such a high level?


title: Non-independent correlation, jacknife on subjects?

text: I have an array of 12 subjects each of whom have 2 variables (i.e.right now I have 24 x 1 array..each subject repeated twice..). When I correlate this, I get a significant correlation. But because each subject is counted twice, I was told to jackknife on the subjects. While I understand the theory, I'm not sure how to implement this in matlab.. So I did the following: p=1:24; [rnhat pnhat]=corr(a,b); for i=1:12 [rn(i) pn(i)]=corr(a(setxor([i i+12],p)),b(setxor([i i+12],p))); end bias=11*(mean(rn)-rnhat); =-0.0073 bias=23*(mean(rn)-rnhat); =-0.0153 Q: Is the above implementation correct? Q: How should bias be calculated? i.e. 11 or 23 or something else? Q: What does this mean? i.e. not sure what the bias should be to consider the correlation significant. Thank you!


title: Differential Entropy

text: Differential entropy of Gaussian R.V. is $\log_2(\sigma \sqrt{2\pi e})$. This is dependent on $\sigma$, which is the standard deviation. If we normalize the random variable so that it has unit variance its differential entropy drops. To me this is counter-intuitive because Kolmogorov complexity of normalizing constant should be very small compared to the reduction in entropy. One can simply devise an encoder decoder which divides/multiples with the normalizing constant to recover any dataset generated by this random variable. Probably my understanding is off. Could you please point out my flaw? Thank you


title: Principle of factor analysis

text: Please is anybody can explain me the principle of factor analysis (exploratory and confirmatory) that I will use in the validation of a questionnaire translated. Do I have to make a survey using the original and experimental Version or only the experimental version? And what is the sample size for this survey?


title: How to use standard deviation to place items in categories?

text: I would like to somehow use the standard deviation of the following set of numbers to categorize values into five different groups. How would I do this? Do I take the max value of the set and divide by four, to get five equal parts? **Set of numbers** 1   1   1   15   22   9   99   100   240   3   7 **Standard Deviation**   75 **Divided by 4 (creates 5 groups)**   * 60 - numbers less than 60 equal category 1    * 120 - numbers less than 120 and greater and equal to 60 equals category 2    * 180 - numbers less than 180 and greater and equal to 120 equals category 3   * 240 - numbers less than 240 and greater and equal to 180 equals category 4   * numbers greater than 240 equals category 5


title: What is a good optimization algorithm for finding the MAP?

text: I'm working on a Bayesian inference package and I have a function that finds the maximum a posteriori point. Right now I'm using scipy's implementation of BFGS to find that minimum of the log posterior, but I've had some trouble with handling large negatives and -infinities. I'd like to see if there's a better algorithm, particularly something that works well on problems that   1. Will be stiff in some dimensions (like scale parameters)   2. Have boundaries (like uniform distributions)    3. Have a largish number of dimensions (say at least 100) I'm looking for suggestions on what optimization algorithms to look at. I'm okay with implementing something myself.


title: Significant Statistical difference One-sample t-test

text: I am trying to test if one algorithm is statistically significantly faster than another using the One-sample t-test These are the results I have and I am trying to prove Algo2 is significantly faster than Algo1                       Mean (ms)  Variance     Algo1     5171       782      Algo2     3753       1920      with a sample size of 78,869 I've tried using the equation: ![one-sample t-test](http://i.stack.imgur.com/wtIEI.png) but I don't seem to be getting a value that makes sense. How do I work this out?


title: AIC and deviance don't agree

text: I'm trying to test the fit of two logistic regression models. I understand that deviance and AIC are good methods for this, but what if one model has better AIC but worse deviance? This is what I think has happened in my case. Here's the summary/anova table from the first logistic regression:               Call:     glm(formula = winperc ~ prevwinperc, family = binomial(link = "logit"),          data = New_Data, weights = gp_col)          Deviance Residuals:         Min       1Q   Median       3Q      Max       -4.9750  -1.0090   0.0503   1.0859   3.7078            Coefficients:                     Estimate Std. Error z value Pr(>|z|)          (Intercept) -1.02740    0.07812  -13.15   <2e-16        prevwinperc  2.04401    0.15453   13.23   <2e-16           (Dispersion parameter for binomial family taken to be 1)              Null deviance: 710.52  on 235  degrees of freedom      Residual deviance: 534.69  on 234  degrees of freedom      AIC: 1842.8          Number of Fisher Scoring iterations: 3          Analysis of Deviance Table         Model: binomial, link: logit         Response: winperc         Terms added sequentially (first to last)                    Df Deviance Resid. Df Resid. Dev     NULL                          235     710.52     prevwinperc  1   175.83       234     534.69      and here's the same info for the second:               Call:     glm(formula = winperc ~ prevPyth, family = binomial(link = "logit"),          data = New_Data, weights = gp_col)          Deviance Residuals:          Min       1Q   Median       3Q      Max       -4.9997  -1.1134   0.0745   1.1310   3.7486            Coefficients:                 Estimate Std. Error z value Pr(>|z|)         (Intercept) -1.13383    0.08564  -13.24   <2e-16 ***     prevPyth     2.25682    0.16969   13.30   <2e-16 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           (Dispersion parameter for binomial family taken to be 1)              Null deviance: 710.52  on 235  degrees of freedom     Residual deviance: 532.75  on 234  degrees of freedom     AIC: 1840.8          Number of Fisher Scoring iterations: 3          Analysis of Deviance Table         Model: binomial, link: logit         Response: winperc         Terms added sequentially (first to last)                     Df Deviance Resid. Df Resid. Dev     NULL                       235     710.52     prevPyth  1   177.77       234     532.75      So the first model has higher AIC but lower Deviance.


title: Online material to learn time series analysis

text: My question is if there are any good online materials for learning this. Something that introduces things well, especially ARMA models and the related math. **Edit:** I'm looking for something of the high-end undergraduate level. Something like in Brockwell and Davis' _Introduction to Time Series and Forecasting_


title: What exactly is the "proportion of variability explained"?

text: One often hears to say _"more than 70% variability is explained by ..."_ What exactly is meant by this? Th proportion of the sum of squares (SSE), or mean sum of squares (MSE)? For example in the following anova table:                                                   Df Sum Sq Mean Sq F value Pr(>F)         as.factor(site)                    444   8357   18.82   163.1 <2e-16 ***     as.factor(year)                     12    569   47.43   410.9 <2e-16 ***     as.factor(month)                     5    863  172.53  1494.8 <2e-16 ***     as.factor(year):as.factor(month)    60    769   12.82   111.1 <2e-16 ***     Residuals                        34188   3946    0.12                        ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1     7176 observations deleted due to missingness      could we say that most of the variability was explained by `site`? We see that site covers most of the SSE but as there is a lot of sites, the MSE for site is almost the lowest in the table. And how would I interpret this in practice? I want to know where is the variability, whether it varies mostly accross time or space. Is the `site` actually the biggest source of variability, or is it a `month` and `year`? Shall I read SSE or MSE column for this purpose? PS: please note I am not a professional statistician, so if you are about to respond with a lot of math then please make also some simple summary for dummies :-)


title: Interpreting F, p, & partial eta squared from an ANOVA

text: I am trying to interpret data from previous research and I don't understand what the individual statistics (F, p, and np2--partial eta squared) are telling me about the data, or how they relate to each other to convey information about the data. For example, here are some of the information from Spruyt et al., 2009: > The 2 groups are people with synesthesia and nonsynesthetes They are each > given 2 Stroop Type color priming tasks, one synesthetic task and one > standard task The levels for each task are congruent and incongruent What I don't understand is how they make inferences based on the following data:   1. A significant color priming effect emerged in the standard task, F(1,22)= 43.87, p<.001, np2= .53   2. Color priming effect was reliable in the nonsynesthetic F(1,22)=13.10, p<.005, np2= .37 as well as the synesthetic group F(1,22)= 24.39, p<.001, np2= .53 I have no idea what the partial eta squared means, but I have come to the understanding that if F is bigger than 1 it is significant.   * But then, what does the p value have to do with it?    * Would a big F and p<.001 be more significant than a big F and p<.005?    * If the p values are the same and both Fs are large but one is 3 times larger than the other, is the larger F more significant?    * How do I know what effect or interaction is more significant or prominent than another?


title: identifying events (patterns) that occur before an event of interest using sequence of events data

text: I'm interested in **identifying events** (patterns) that **occur before an event of interest**. For example, a customer calls in to complain or a customer checks their balance online then the customer closes the account (event of interest). Data is in a form of a sequence of events with a time stamp. I would like to know what methodologies (and software) people are using for such a task. I'm open to **exact patterns** (sequence of events), **frequency of events** (calls in multiple times to complain), and **timing of events** (80% of event A leads directly to event of interest within X days). Basically, open to any approach people are using to identify patterns leading up to an event. Thus far I found the `CSPADE` algorithm as available in the `arulesSequences` package in R. It seems to be able to identify sequence of patterns and which items co-occur. However, I don't think one can set a target event for it to find a pattern. I'm open to algorithms available in R, Python, or SAS. Thanks so much!


title: What to do if a covariate is collinear in one model but not in another

text: I have 10 different species with presence/absences data, as well as 6 different covariates relating to the design of marinas, including 3 continuous (lengths of walls, pontoons and groynes) and 3 factor variables (distance from freshwater, type of entrance and marina location). I have made binomial GLMs to assess which design elements of a marina which most influence the presence of each species. However, when I run variance inflation factor (VIF) values to assess if there is collinearity among the covariates, the covariate ‘location’ only comes up as collinear in some of my models, even tho they all contain the same covariates... why is this? Should I just exclude 'location' from all my models (I can see why it would be collinear)? Or should I include/exclude it from different models based on the models VIF? The following are two of my 10 models, the remainder follow the same format but with different response variables (all binary presence/absence data)               bin.fit1<- glm(S_clava ~ as.factor(fresh) + wall + groyne + pontoon + mooring + as.factor(entrance) + as.factor(location), family= binomial, data= data1)          bin.fit2<- glm(C_mutica ~ as.factor(fresh) + wall + groyne + pontoon + mooring + as.factor(entrance) + as.factor(location), family= binomial, data= data1)


title: What would be the logical inverse of a hypothesis test? (If any)

text: I just thought I'd ask this question to see if there is a method or reasoning of proving a hypothesis test inversely, whilst logic remains the same. I will use the following example, which I recently got confused over: > Question: Test on a 5% level of significance if the model is significant. > State the null hypothesis, alternative hypothesis, the p-value and > conclusion.   1. The correct answer:      * H0(null): Model is not significant     * H1(alt) Model is significant     * p-value = 0.4 > alpha (0.05)     * Do not reject H0     * model is not significant I understand the logic of this answer, but I alternatively wrote an answer like this:   1. My answer:      * H0(null): Model is significant     * H1(alt) Model is not significant     * p-value = 0.4 > alpha (0.05)     * Reject H0     * model is not significant After reviewing my error, I understand that we cannot reject H0 if the p-value is greater than the 5% level of significance, and this is where my question comes in. How would it be possible to reject MY H0? Would I use an alternate value? My own suggested answer to this (based on the 5% level of significance) is that if I have inverted everything, then I would invert my 5% level of significance and say that: p-value = 0.4 < 0.95 Therefore, reject H0 on a 5% level of significance. My reasoning is deduced from the fact that (based on the 5% level of significance) it must hold that the p-value be greater than 0.95 to accept the null hypothesis. Can anyone confirm or reject this?


title: How to estimate the parameter of this time series?

text: The time series is governed by the equation $S(T)=S(0)e^{(\mu-\frac{\delta^2}{2})T+\delta(w(T)-w(0))}$, in which $w(t)$ is a standard Brownian motion. Now given the data $\{S(t)\}_{t=0}^{t=T}$， how to estimate $\delta$ and $\mu$?


title: Setting up data for MANOVA in SPSS

text: I do not understand how to arrange my data in SPSS for running MANOVA.   I have 6 groups of participants that have rated 100 pictures on 3 different dimensions (dim1, dim2, dim3). What I don't understand is this: Do I have to use the mean value for the 3 dimensions (dim1, dim2, dim3) (say, mean of 120 participants) and then make every picture a variable (p1, p2......p100)?   Or do I have to again get the overall mean of all the pictures?


title: Cluster standard error _versus_ fixed effects

text: I am working with experimental data -correspondence test in the rental housing market-. On which ground can I choose between:   1. F.E. for the region where the apt is and cluster s.e. on the day the application to the apt was sent    2. cluster s.e. on region where the apt is and F.E. for the day the application was sent (day is too much, perhaps, better w-e or month)? More in general. On what ground can I choose whether to control for a variable or to cluster the s.e. on that variable? The first option changes the estimate of the coefficients while the second changes the estimate of the s.e.'s....both of them should be used when I suspect there is something which is substantially non observable. So, I know the practical difference between option `1` and `2`, but on what ground do I choose?


title: t-test to compare two means

text: I have asked question concerning comparison between two groups. Got really nice answers and I want to perform downstream analysis using really good gui11aume's suggestion.   But not everything is clear to me, hence I need help to clarify significance test I should apply.   **What I want:** To compare two means of percentage and find out if both compared groups contain same percentage of analyte.   **What I have:**               Group   Mean   Size       A     0.1     5       B     2.3    124      **Question:** How can I do it? gui11aume suggested to compare two means using t-test just like in this example. However I found that for t-test one can also use full string (not only pre-calculated mean). So, if it's possible to calculate t-test statistics from `what I have` - how can I do it (in R)?   And what is the difference between such t-test and cbeleites answer using `prop.test`? They seem similar to me.


title: Doing CFA on a known theoretical model, but having problems with convergent and discriminant validity

text: I did a survey based on the UTAUT model / theory with standard questions used. The N was 150. I am trying to validate the model for my survey with structural equation modeling (via the AMOS software). But I have heard that before testing the model fit, "it is absolutely necessary to establish convergent and discriminant validity, as well as reliability, when doing a CFA" (confirmatory factor analysis). Doing CFA, I found out that factor loadings are OK (using AMOS, loadings > 0,7). I did CR for reliability testing, everything OK - over 0,7. But (not like other published papers) I have problems with convergent and discriminant validity (note results below). ![enter image description here](http://i.stack.imgur.com/8FUnX.png) Is there some way to deal with this or do I have to restate the questions and repeat the survey? (I think doing exploratory factor analysis is not logical if so many others have proven the theory works, so does that mean something must be wrong with my data?)


title: How can I set a maximum and minimum level for a dependent variable?

text: I have to make sure that a dependent variable I explain using linear regression ranges between a minimum of 0% and a maximum of 30% (it is an investment weight in a portfolio). How should I proceed ?


title: What is the criterion of model validation?

text: Suppose I have a model that boasts given 1 year's daily data to calibrate the parameters, it could predict the behavior of future 1 month. What should be the criterion for back-testing or, model validation? Currently I calculate the error residuals and do Hypothesis testing:   1. Calculate $\epsilon = x - \hat x$ , where $x$ is the realized value, and $\hat x$ is the estimated value based on the model $X(\theta)$ , while $\theta$ are those parameters calibrated based on 1 year's data.   2. Since the model declares $\epsilon$ follows normal distribution: $$E\{\epsilon\} = 0$$ and $$Var\{\epsilon\} = \sigma^2 $$ I do the Hypothesis testing on the next month's data (28-31 days). Now, suppose I did such Hypothesis testing on the previous 5 years, during the 60 times hypothesis testing (significance level $\alpha=5\%$), 2 times failure in testing $E\{\epsilon\} = 0$, and 3 times failure in testing $Var\{\epsilon\} = \sigma^2 $, where there is a "bad month" both mean test and variance test failed. Should I accept this model or reject it and request improvement? Is there some guidelines or common practice on such issue?


title: Can I select a subset of predictor variables without using step()?

text: I have a set of 20 predictor variables, and I want to formulate a regression model by applying my own variable selection technique basically with backward approach (just for an experiment purpose.) I'm rather new to `R`, and as far as I have studied so far, if I were to do it automatically I would use `step()` like this:               y2y <- read.csv(my CSV location)     y2y.lm <- lm(myTarget~., data=y2y)     y2y.step <- step(y2y.lm)      However, I'd like to apply my own technique to select variables using the calculation results based on p-values, and so far I have figured out that I can obtain the largest p-value and the corresponding variable name by               maxV <- which.max(summary.lm(y2y.lm)$coefficients[,"Pr(>|t|)"])     maxV     attr(maxV,"names")      If I do the calculation repeatedly (and stepwise) in a loop until a condition is met, using `for()` or `while()`, then I can probably achieve my goal. However, I just can't figure out the way to formulate the model dynamically to feed to `lm()` in the loop, e.g. `myTarget ~ x1 + x2 + x3 + x4.....` using the set of variables I selected. For instance, in the first loop, if `x2` was selected as a variable to drop, I want the equation to look like `myTarget ~ x1 + x3 + x4.....` Is there a way to dynamically formulate the equation to feed to `lm()`? I'd appreciate your suggestions and insights.


title: Variable not significant in log specification

text: I am estimating a model which gives me significant parameter estimates for my variables of interest. However, when I use the log version of the model log(dependent)=log(independent), I don't get significant estimates at the 95% confidence level. I am replicating two other studies, and the authors do get significant results under both specifications. Can someone please explain why does this happen and what does it mean? Any way of solving the problem?


title: How to express the difference between two percentages?

text: If I have some share of a product in some sector as:               2012 2013     50%  60%      and I find share change for 2012-13 then should it be 10% or 0.1?


title: Difficulty in understanding a vector quantization algorithm

text: Neural gas for vector quantizationpaper explains a technique for symbolizing or quantizing data. Algorithm presents the algorithm in Section 4. An application in EEG data symbolization is presented Application. In the Application it is shown that an n- dimensional data like EEG recordings is vectorized to 1D. But Neural gas is basically an unsupervised clustering algorithm where data can be assigned to more than one cluster. The definition of Vector Quantization (VQ) is the process of quantizing n-dimensional input vectors to a limited set of n-dimensional output vectors referred to as code-vectors. The set of possible code-vectors is called the codebook. The codebook is usually generated by clustering a given set of training vectors (called training set), the codebook is then used to quantize input vectors. What I do not understand is that how a multi dimensional data is converted to single dimension as in the Application paper and the module where the symbols are assigned. In the Algorithm, I did not understand what is "reference vector $w_{c_i}$, the set $A$ and Step 3 where it says to order all elements of $A$ according to distance. What is the distance measure here? Can somebody please explain the algorithm in simple terms how the symbolization is done? Thank you


title: Can I use data for factor analysis with very low KMO (less than 0.4), determinant negative, and data not significant according to Barnett's test?

text: I am analyzing rainfall water quality with cluster analysis first and then principal component analysis (PCA). There are 7 variables. Results from cluster analysis that extracting with PCA, almost of them show negative determinant, KMO very low (less than 0.3) and according to Barnett's test on sphericity, those data truly not significant (only 10% of group samples have p value <0.05). Then, I checked on the correlation matrix and found that some of them are very highly correlated and some very low. To raise good KMO and more significance, I dropped the least variable, but the results are not very different. Although KMO became little higher, the significance test showed p value over 0.400.   1. Can I use these data to get factors?   2. What I must do with results like that?


title: WinBUGS: Multiple definitions of a node

text: So this question is about the BUGS modeling language. So you either know it or have no clue. I'm a newbie to this so it's been driving me mad. I want to define a simple two-state hidden Markov model (HMM) where the emission of each state follows a Normal distribution. I have an array data of Nxl dimension where each row is a subject and column is a continuous variable Time associated with the subject. So in WinBUGS, I would define the model as:                for(i in 1:N) { # for each subject         # Sample the initial observation         Time[i, 1] ~ dnorm(mu[State[1]], tau[State[1]])              for(j in 2:l[i]) { # for each observation             # Sample the observed variables             Time[i, j] ~ dnorm(mu[State[j]], tau[State[j]])             # Sample the hidden states             State[j] ~ dcat(P[State[j-1], ])         }     }          # P is the transition matrix of the hidden Markov chain     P[1, 1:2] ~ ddirch(Pinit)     P[2, 1:2] ~ ddirch(Pinit)          # Set the initial states     State[1] ~ dcat(Pinit)          # Sample the initial params     mu[1] ~ dnorm(200, 1.0E-6)     mu[2] ~ dnorm(400, 1.0E-6)          # The precision params     tau[1] ~ dgamma(0.001, 0.001)     tau[2] ~ dgamma(0.001, 0.001)      Initial values have also been provided via R2WinBUGS. The model is syntactically correct. But when I run, I got this error: " **multiple definitions of node State[2]** " Can you please tell me why and how to solve this? I've searched around on the error but each has their own specific case and there's not a generic solution or explanation as to why this arises.


title: How to model the probability of an event occurring on a given day when it must occur once in a set of 10 days?

text: I am a muslim and there is ramadan going on. The last 10 days of ramadan are really special, and in one of these 10 days there is one special night that no one knows about.   The special day is more likely to happen on the odd nights, but it can be on any day. **How can the probability of a given day being the special day be modelled statistically?** **My logic:** Since we know it needs to fall in one of the 10 days and these all equally likely to land on any day. So the probability of that happening is 10%. I would think the binomial theorem would play a role because either it happens that day or it doesn't. But when I calculated it all the way to n =9 (0 included) my probability got really low as I approached the last day.


title: How to run a Granger Causality Test with Stata

text: I'm doing a study on the determinants of FDI (Foreign Direct Investment) in the ASEAN countries. Before doing a panel data analysis, I'd like to run a Granger Causality Test between the potential FDI determinants time series (GDP, exchange rate, ecc.) and the FDI's one to support the choice of these variables for the above panel analysis. My questions are:   * Does it make sense?   * If the time series are non-stationary (I've already discovered it with a ADF test) could I (in any case) run the Granger CT or should I have to make time series stationary with some cointegration process before? (and in this case, which one?)   * ( _the most important_ ) How can I do it with Stata? (I think I have to run a VAR or SVAR before doing the Granger CT, is it correct? does it make sense in non-stationary (unit root) time series?) how should I set up the test? (dependent variable, lag periods, ecc.)   * ( _and finally_ ) How can I interpret the (Stata Granger CT) results?  thanks in advance to everybody


title: Analysis with complex data, anything different?

text: Say for example you are doing a linear model, but the data $y$ is complex. $ y = x \beta + \epsilon $ My data set is complex, as in all the numbers in $y$ are of the form $(a + bi)$. Is there anything procedurally different when working with such data? I ask because, you will end up getting complex covariance matrices, and test statistics which are complex valued.. Do you need to use a conjugate transposes instead of transposes when doing least squares? is a complex valued covariance meaningful?


title: Is it correct that no statistical simulation is done using $2^{256}$ (or more) outputs from a MWC prng?

text: In a conversation about Marsaglia's Multiply-With-Carry pseudo-random number generator and its potential use in generating random data for statistical purposes, someone said: > There is no statistical simulation ever done that will use $2^{256}$ outputs > from such a generator. Is that correct?   * If _"yes"_ — Is there an explanation which might help me understand why that is correct? And is there any other statistical purpose where more than $2^{256}$ of a MWC prng output would be used?   * If _"no"_ — Can you give me a name of a statistical simulation I could mention as an example?  The reason why I'm asking is because (being a stats-noob) I trusted in the fact that Marsaglia (being a stats-guru) had a good reason to develop the MWC prng with its _"near-record period"_. But now, I have a hard time understanding for what statistical purposes the MWC prng might come handy if it's indeed true that no statistical simulation will ever use $2^{256}$ (or more) outputs from such a generator.


title: How can root mean square error be used to predict logistic regression model accuracy?

text: The problem which I have at hand consists of a logistic regression model for risk evaluation that has been made on some credit card data of Quarter-1'12 (Jan'12 - Mar'12). Now I use the same model to evaluate risk for the data of Quarter-2'12(Apr'12-Jun'12). I want to devise an accuracy score that gives an insight into how "accurately" my model fits the new data (Quarter -2 data). I have used Hosmer Lemeshow Statistic and Balanced Accuracy Method till now but none served the purpose. What further can be done in this regard?


title: Introduction to statistics for mathematicians

text: What is a good introduction to statistics for a mathematician who is already well-versed in probability? I have two distinct motivations for asking, which may well lead to different suggestions:   1. I'd like to better understand the statistics motivation behind many problems considered by probabilists.   2. I'd like to know how to better interpret the results of Monte Carlo simulations which I sometimes do to form mathematical conjectures. I'm open to the possibility that the best way to go is not to look for something like "Statistics for Probabilists" and just go to a more introductory source.


title: Using AdaBoost on multi-class in R on unbalanced data

text: I have a data set which is highly imbalanced and I have used the SMOTE algorithm (using the R package DMwR) to balance the binary class in the data set. I have been using the R Ada package to then train an Ada Boost model on this data set to predict the binary class, with very good results. In the same data set, I have another class variable which has multiple values (6 in total). In this case I realise that I can't use the AdaBoost algorithm as implemented in the ada package as it only deals with the binary case. I therefore have 2 problems:   1. I'd like to use the SMOTE algorithm on second class variable but this also only works with binary classes. Is there an algorithm or package I can use in R to "rebalance" a data set based on a class with multiple values in a similar way to SMOTE?   2. I'd like to use a classifier to predict the multiple class variable. I have tried using the one-vs-all approach with AdaBoost but I cannot get this to work well (my approach is below). Boosting seems to work well with this data set. Are there any other boosting algorithms or other approaches I could use in R that handle classes with multiple values. I have tried using Random Forest but one of my nominal inputs has too many discrete values to use it. **Approach for AdaBoost one-vs-all**   * Build a vector with a binary variable for each discrete class value   * Train one AdaBoost model against each binary class vector   * Generate probability prediction for each AdaBoost model   * Select the class with the highest probability Many thanks


title: How to include a linear and quadratic term when also including interaction with those variables?

text: When adding a numeric predictor with categorical predictors and their interactions, it is usually considered necessary to center the variables at 0 beforehand. The reasoning is that the main effects are otherwise hard to interpret as they are evaluated with the numeric predictor at 0. My question now is how to center if one not only includes the original numeric variable (as a linear term) but also the quadratic term of this variable? Here, two different approaches are necessary:   1. **Centering both variables at their individual mean.** This has the unfortunate downside that the 0 now is at a different position for both variables considering the original variable.   2. **Centering both variables at the mean of the original variable** (i.e., subtracting the mean from the original variable for the linear term and subtracting the square of the mean of the original variable from the quadratic term). With this approach the 0 would represent the same value of the original variable, but the quadratic variable would not be centered at 0 (i.e., the mean of the variable wouldn't be 0). I think that approach 2 seems reasonable given the reason for centering after all. However, I cannot find anything about it (also not in the related questions: a and b). Or is it generally a bad idea to include linear and quadratic terms and their interactions with other variables in a model?


title: Does applying ARMA-GARCH require stationarity?

text: I am going to use the ARMA-GARCH model for financial time series and was wondering whether the series should be stationary before applying the said model. I know to apply ARMA model the series should be stationary, however I'm not sure for ARMA-GARCH since I'm including GARCH errors which imply volatility clustering and non-constant variance and hence non-stationary series no matter what transformation I do. Are financial time series usually stationary or non-stationary? I tried applying ADF test to a few volatile series and got p-value<0.01 which seems to indicate stationarity but the principle of volatile series itself tells us that the series isn't stationary. Can somebody clear that up for me?I'm getting really confused