title: Finding conditions on unspecified CDF that permit a solution to an equation

text: [A duplicate thread can also be found at http://mathoverflow.net/questions/131142/finding-conditions-on-unspecified- cdf-that-permit-a-solution-to-an-equation ] Let $F(\alpha) := \mathbb{P}(\tilde{\alpha} \le \alpha)$ be an arbitrary, strictly increasing and twice differentiable CDF that is defined on the interval $[0, \overline{\alpha}]$, where $\overline{\alpha}>1$ may be infinite. Moreover, let $\mathbb{E}(\tilde{\alpha}) = 1$. Let $N \ge 2$ be a natural number, and $\delta \in (0,1)$ real. I am looking for necessary and sufficient criteria on $F$ in order to find a solution to $(N-1)[1-\delta+\delta F(\alpha)] - \delta \alpha F'(\alpha) =0$, where $\alpha$ can be any real positive number not larger than $\overline{\alpha}$. Now, it is clear that the left term, $(N-1)[1-\delta+\delta F(\alpha)]$, is strictly increasing in $\alpha$ and bounded above by $N-1$. Hence, if the probability mass is sufficiently "concentrated" in some interval (implying that $F'(.)$ is large in that interval), a solution $\alpha^*$ must exist by continuity of $F'$ and $F$. However, it would be nice to have some sharper conditions on $F$ that are necessary or sufficient for a solution. Ideally, I'd wish to have a result that gives an upper bound on $F$'s variance, or similar. One thing that seems problematic is that the above equation never uses the fact that $\mathbb{E}(\tilde{\alpha})=1$. I've tried to apply Markov's inequality but it didn't help me much. Numerical simulations reveal that a solution can usually be found if the variance of the considered CDF is sufficiently low. Examples include $F$ log- normal or $F(\alpha) = \left(\dfrac{b\, \alpha}{b+1}\right)^{\displaystyle b}$ for $b>1$. The uniform distribution on $[0,2]$ is a borderline case that doesn't permit a solution. I would greatly appreciate any help or ideas to this (economics) research problem. Many thanks!


title: Do transformations prevent standardised coefficients from being comparable?

text: I assume that using standardized regression coefficients is a way to compare the relative effect of two or more predictors (measured on different scales/units) on the mean value of Y. **If you have several predictors that are comparable without standardization, and you transform one predictor to correct for skewness, would standardized coefficients of these two predictors be comparable?**


title: Finding the probability of given model

text: I'm trying to model a system given as a project, any help or advice is useful for now. Here our assumptions for a part of it:   * If I study for the finals and do assignments, I have $80\%$ chance to get a scholarship   * If I neither do assignments nor study, my chance reduces to $20\%$   * If I either do the assignments or study for the final, I have $60\%$ chance to get a scholarship.   * The chance that I will do assignments or study are both $50\%$. So, given that I got the scholarship what is the probability that I have did the assignments?


title: Models for calculating consumer behavior at coffee shop

text: I have the occasion to sit in a Starbuck's almost every day. I have noticed there are rush hours sometimes. It's like hundred of people decided to buy something at Starbucks at the very same time. They rush in, the line becomes longer and longer, and about 15 minutes later they disappear until the next rush hour starts about 45 minutes to 2 hours later. Also there are days when many people choose to stay inside of the Starbucks for like 30 minutes to 2 hours. One week later, the Starbucks might be totally empty under the very same conditions (depending on the day of the week, weather and temperature, school holidays, etc.). To me it seems that this process is totally random, at least I don't see any connections. I'd like to calculate how possible it is that my Starbucks will be full of people today. Are there any models for this or ways to calculate the consumers' behavior?


title: Does using box-cox transformation on individual data sets prevent these data from being comparable?

text: I have been using the box-cox transformation to normalise data for input to an Ecological Niche Factor Analysis software, as recommended by the software creators. However, it has occurred to me that the box-cox transformation method has (obviously!) been selecting different lambda values for each transformation. For example, I want to compare the influence of factors A,B and C on individual location over four sample dates. Each factor has been separately box-cox-ised (separately both from other factors and other dates). Does this mean that the results of each factor analysis will not be comparable (and that simple comparison of the transformed data, e.g. by ANOVA, will not be possible) because of the selection of different lambda values for each transformation? Many thanks in advance!


title: How are Kenward & Roger's degrees of freedom computed in PROC MIXED?

text: I run the following SAS model:               PROC MIXED DATA=dat order=formatted ;         CLASS LOT DOSE ;         MODEL TITRELOG10 = DOSE  / s DDFM=KR ;         RANDOM DOSE / subject=LOT type=UN ;     RUN; QUIT;      and it yields the following results: ![results_UN](http://i.stack.imgur.com/saTzJ.jpg) Now I run the same model except that I use a Toeplitz variance instead of an unstructured variance:               PROC MIXED DATA=dat order=formatted ;         CLASS LOT DOSE ;         MODEL TITRELOG10 = DOSE  / s DDFM=KR ;         RANDOM DOSE / subject=LOT type=TOEPLITZ ;     RUN; QUIT;      it yields: ![results_TOEP](http://i.stack.imgur.com/cq10C.jpg) The behavior of the Kenward & Roger degrees of freedom sounds strange: with the unstructured matrix they are all equal, whereas this is not the case with the Toeplitz matrix. The results are exactly the same with the Satterthwaite degrees of freedom (moreover, isn't it strange too that Satterthwaite and K&R degrees of freedom are the same?) I am not confident that I have mastered the degrees of freedom derivation and I am rather suspicious about the way they are derived in SAS, hence I would like to have your opinion about that. EDIT: My data are balanced


title: How to convert standardized coefficients to unstandardized coefficients?

text: My goal is to use the coefficients derived by previous research on the subject to predict actual outcomes given a set of independent variables. However, the research paper lists the Beta coefficients and t-value, only. I would like to know if it is possible to convert the standardized coefficients to unstandardized ones. Would it be useful to convert my unstandardized independent variables to standardized ones to calculate the predicted value? How would I return to an unstandardized predicted value (if that is even possible..) **Added** Sample row from paper: > Number of bus route (buslines) | 0.275(Beta) | 5.70*** (t-value) I am also given this regarding the independent variables: > Number of bus route (buslines) | 12.56(avg) | 9.02(Std) | 1(min) | 53(max)


title: Correlate bivariate Brownian bridges

text: Given two independently constructed Brownian bridges (from their marginal means and variances), is there a way to correlate the sample paths?


title: Lomax distributions - Kullback Leibler divergence

text: Does anyone know of a reference for an expression for the Kullback-Leibler divergence between two Lomax (Pareto II) distributions? Not really worried which way the Lomax is parameterized.


title: Can we trust the results of PROC MIXED when the G matrix is not positive definite?

text: PROC MIXED in SAS sometimes estimates a non-positive definite matrix of the random effects. However PROC MIXED nonetheless provides the estimates and the standard errors of the fixed effects in such a case. Can we consider that these estimates and standard errors are correct ?


title: How to plot multiple users' deviations from predictions of bandwidth consumption over time?

text: As it is said, a picture is worth 1000 words. I have some data, and I need to put them in the right chart representation. ### The working scenario is this: some users have to predict how much internet bandwidth they will consume during 13 weeks. The users can form groups in order to buy more cheaper internet bandwidth. If they consume more than they said, they will have to pay more. Idea is to be as near is possible by what was predicted. ### So, my data is like this:               user1 = {-0.075, -0.09, 0.32, -0.242, -0.368, -0.401, -0.73, -0.367, -0.294, -0.043, 1.296, 0.075, -0.373}       user2 = {-0.009, -0.013, -0.01, -0.008, -0.008, -0.01, -0.005, -0.02, 0.287, 0.345, -0.104, -0.324, 0.144}        user3 = {-0.197, -0.271, -0.153, -0.621, -0.549, -0.09, 1.745, 0.436, -0.271, 0.093, 0.085, 0.211, 0.331}       user4 = {-0.005, -0.005, -0.006, -0.006, -0.006, -0.005, -0.006, -0.086, -0.171, -0.15, -0.175, -0.067, 0.078}       user5 = {-0.223, -0.048, -0.129, 0.14, -0.535, -0.29, 0.51, 0.801, 0.521, 0.482, -0.105, 5.082, 5.516}       group = {-0.509, -0.427, 0.022, -0.737, -1.466, -0.796, 1.514, 0.764, 0.072, 0.727, 0.997, 4.977, 5.696}       Negative values represents that user consumed less that what he predicted. Positive values means that user consumed more that he said. Values from the group array is the sum of the difference between what was predicted and what was consumed by each user. ### Question: Do you have any idea, about how can I represent those data in a chart in a easy to read way? For example, as you can notice, difference between what was forecasted and what was consumed for user2 and user4 are very small, comparing to user2 and user5. A CVS file with the previous data is here.


title: identify different points

text: I have a large scatterplot, with about 100,000 (x,y) points. The x coordinate is the set of numbers from (1 to ~100,000) - in other words, no 2 points have the same x-coordinate. The y is mostly constant (around 50-70 in value), but there are key "regions" where the y value spikes to ~120 or drops to ~20. How would I statistically differentiate these regions? Clustering? Any other ideas? As a bonus, if you know R, if you could reference a certain method that would be helpful as well.


title: how to extract AIC(Akaike's Information Criterion) in LAR(Least Angle Regression) in R Studio?

text: I'm already done in conducting the whole LAR Algorithm using lars() function in R Studio. But my problem is how to extract or use AIC in R Studio for choosing enough the number of variable that will explain my model. If i will use the predict.lars() function, only the Cp Criterion will appear. Anybody? Thank you in advance.


title: I need the whole Lasso(Least Absolute Shrinkage and Selection Operator) Algorithm. Anyone?

text: Anyone could give me the whole Lasso Algorithm? The step by step procedure in determining Lasso solutions? I barely need this for my thesis. Thank You.


title: Controllong random deviation

text: I have a process in which X is a molecule going through birth ($x \rightarrow x + 1$) and death $(x \rightarrow x - 1)$. Rate of both reactions, namely synthesis rate for birth and elimination rate for death reaction, are proportional to number of molecules of X ($[X]$) at any given time. In this system, assume that synthesis rate itself increases with $[X]$; now I want to control the "random deviation" in $[X]$. Now I can't write down the expression which can tell me at what rate my elimination rate should increase to check the random deviation?


title: Is the Stoppa distribution a special case of any better-known, more general distribution?

text: The Stoppa distribution is a 3-parameter distribution that generalizes the Pareto distribution, adding a second shape parameter but no location term. The CDF is $$F(x) = \left[1-\left(\frac{x}{x_0}\right)^{-α}\right]^θ;\quad 0 < x_0 \leq x$$ Definition from Kleiber & Kotz (2003). This distribution is becoming increasingly important in the study of income distribution at the high end. Is there any better-known distribution for which the Stoppa is a special case, given some restriction on the parameters of the more general distribution? And if so, what is the distribution and the restriction? I suspect it is a special case of the five-parameter distribution that McDonald & Xu (1995) calls the Generalized Beta, Crooks (2012) calls the generalized beta prime, and some other authors call the Feller-Pareto (though all these names have also been used for the four-parameter distribution McDonald calls the generalized beta of the second kind). But I have not yet been able to find a set of restrictions that gets me to the Stoppa, so this may be wrong. And I have not been able to identify it with any member of the Pearson, Burr, Dagum, etc. distributional families, though I could easily miss it if it is just parameterized a little differently. Any information or helpful advice would be much appreciated. Crooks (2012). Survey of Simple, Continuous, Univariate Probability Distributions. Kleiber & Kotz (2003). Statistical Size Distributions in Economics and Actuarial Sciences. McDonald & Xu (1995). A Generalization of the Beta Distribution with Applications.


title: How to retrieve the prediction equation in R?

text: I have developed a prediction model prototype in R. The model uses Support Vector Regression to predict. But I need to develop the entire solution in Visual C++ for a real life implementation.   I have prior experience of coding normal regression equations for implementation.   * For Example:    In Linear regression after the modelling, if one gets a equation: Y = 0.5x + 0.8 . one needs to code out the equation for all the incoming 'x' values to get the 'Y' value. But here in I cannot retrieve the prediction formula to code it out. More ever in SVR has many support vectors created.   1. How to retrieve the prediction formula for a SVR in R(if possible as in normal Regression).   2. The model that i am using is: svm_model = svm(Y~., data= trainData, type= "eps-regression", kernel="radial")   3. How to construct the final prediction equation for SVR?


title: Calculating effect size (partial eta squared) for a planned comparison effect

text: I have a between-subjects ANCOVA (factors A and B, and covariate C), for which A*B is significant. I decomposed the interaction with a planned comparison described by the contrast A=(-1,1), B=(0,1), which was significant. I would like to find out the effect size (partial eta squared) for this effect, but unfortunately my stats package (Statistica) does not display this statistic for effects obtained via a contrast (planned comparison). So I set out to compute it manually. I used the definition np^2=SS_treatment/SS_total. Statistica only displays SS_effect and SS_error, so as SS_treatment I took the SS_effect reported in the planned comparison, and as SS_total I took the SS_effect reported for the A*B interaction. This, however, gives me a value of about 3.0, whereas I was expecting partial eta squared to be in the order of magnitude below 1, which is what I get for other similar effects. Is the way I computed the partial eta squared incorrect, or have I otherwise made other mistakes? Thanks for any help!


title: Exposing R functionality as Web services method

text: I have a R script which does some statistical modelling. I would like to expose this script as web service method which I would consume from a mobile application. How can I achieve this?


title: Visualise significant differences in a contingency table?

text: I'm working with pandas and scipy, and I would like to find a quick way to visualise unusual values in a contingency table. Say this is my contingency table:                                 wolves   lions   tigers   chickens    total     ate_zebra         14112    13891   12978    54          41035     ignored_zebra     6403     6377    6109     2637        21526     total             20515    20268   19087    2691        62561      It looks as though chickens are (mostly) vegetarian, unlike the other animals. If I run a chi-squared test, it tells me that the chi-squared value is significant and the expected values are:               wolves      lions       tigers      chickens    total     13456.19515 13294.18296 12519.54165 1765.080242 41035     7058.804846 6973.817043 6567.458353 925.9197583 21526     20515       20268       19087       2691        62561      But what I'd really like is some way to quickly, visually identify the most significant differences from the expected values, say something like this: ![contingency table with most significant deviations from expected values highlighted](http://i.stack.imgur.com/gEiwo.png) Excel would be good enough, but something Python-based would would be ideal :) Or do I need to write my own algorithm to do this? Any suggestions much appreciated.


title: Decision tree in R

text: I am new to machine learning in R. This is my data set.               channels<-sample(c("AFFILIATE","DIRECT","DISPLAY"),100,T)          booking<-sample(c("N","Y"),100,T)          placements<-sample(c("R","L","TR","TL"),100,T)          site<-sample(c("www.google.com","www.yahoo.com"),100,T)          campaigns<-sample(c("camp1","camp2","camp3"),100,T)          data<-data.frame(channels,campaigns,site,placements,booking)      I want to figure out which combination has higher booking probability. I have used the following R code to compute a classification tree.               tree<-rpart(booking~channels+campaigns+site+placements,data=data,method="class")      It generates one output, but not in the proper order (I want a tree where the order should be channels → campaigns → site → placements → booking). Also, it only gives two leaf nodes for each parent node, but I want to generate a tree where there will be more than two leaf nodes for each parent node. Are there any advanced methods which can fulfill the requirements?


title: Is linear discriminant analysis (LDA) more likely to overfit than support vector machine (SVM)?

text: I went to a short talk and the speaker quickly mentioned something like 'LDA (linear discriminant analysis) is more likely to be overfitted than SVM (support vector machine)'. Is this true? And why?


title: In ANOVA, why are $\text{MS}_\text{between}$ and $\text{MS}_\text{within}$ defined by divisions by k-1 and N-k?

text: I realize that the divisions are meant to take into account the number of factor levels (which will tend to push up $\text{SS}_\text{between}$), and the number of observations (which will tend to push up $\text{SS}_\text{within}$). However, I can't understand why we don't obtain these Mean Squares by dividing $\text{SS}_\text{between}$ by k and $\text{SS}_\text{within}$ by N. For example, in the following data there were three levels to the factor, and 18 observations in total. ![Example of an ANOVA table](http://i.stack.imgur.com/7MKDj.png)


title: How to know whether a set of Likert items form two scales using SPSS? How to assess the relationship between these two scales?

text: I have designed a questionnaire which has 16 Likert item questions; from 1 (strongly disagree) to 5 (strongly agree). The 16 Likert item questions are divided into 2 groups: 10 questions for customer loyalty, and 6 for customer satisfaction.   * **How do I know whether the Likert items should be grouped into the two scales?**   * **How can I calculate the relationship between the two scales?**


title: Cross-validation with dummy variables?

text: Does it make sense to use cross-validation with factor variables that have 3+ levels? When using bestglm, I get an error saying that it doesn't work with categorical variables. In the documentation the reasoning is "Cross-validation is not available when there are categorical variables since in this case it is likely that the training sample may not contain all levels and in this case we can’t predict the response in the validation sample." but I can't find anyone else discussing this issue.


title: Missing Value in Commodities Prices

text: I am trying to model the prices of four energetic commodities with ARIMA models in R. Unfortunately the price series is not regular, as for some days, like Christmas, no price is given. My series is made of 2277 days(no weekends) and I have 81 days where the prices of all the commodities don't exist. Moreover the missing prices happen altogether. I read that some fill in the gaps with averages of similar days. Would it be wrong to input for example the average of the previous and next week price (Monday with mondays and so on)? Or is it like imposing a greater seasonality than what the series really have?


title: How I can measure an error in a bar plot?

text: ![enter image description here](http://i.stack.imgur.com/Ruk10.png) This is my barplot. Model(red) is my estimated frequency at each X-axis vale (0-19) and actual is the frequency at each X-axis vale. How I can measure error of this estimating work in barplot?   If I use mean percentage error it does not consider the each size of bar. How should I solve that?? Of course, I have a data about expected frequencies and actual frequency at each value (X-axis).


title: Use lag operators with interaction term, is it correct?

text: I'm sorry I make reference to Stata here, but it's just to illustrate the point, the question is about whether it is correct or not using lagged variables in interaction terms. I have a linear model that I am estimating using FE and RE, in which my explanatory (continuous) variable is lagged, something like this: $$ y = \alpha + \beta x_{t-1} + \epsilon $$ Stata: `xtregar y L.variable, fe/re` For my analysis, I need to include a categorical interaction term. Stata, by using the graphical interface, does not allow me to use lags with my continuous variable, but it does give me results if I add 'L' manually, and these are very different to those without the L. What Stata's graphical interface allows me to do is:               xtregar y L.variable dummy#c.variable i.dummy, fe      And what I think I need is:               xtregar y L.variable dummy#cL.variable i.dummy, fe      I am confused since I imagine that if adding lags was OK, Stata would allow me to do so, wouldn't it?   1. Is it correct to add lags to an interaction term?   2. In case of that the answer to (1) is **no:** If I remove the variable (`L.variable`), is it ok to add lags to the interaction term? A related question: What to do if the variable is not significant but the interaction term is?


title: Problems in scale Bayesian network mode using R

text: The problem that we have is as follows. We have close to 60 discrete random variables each of which shall take on an average of 5 categorical values. We have developed a Bayesian network representation using our domain knowledge. We have the data for these 50 discrete random variables and how they interplay with each other using some logs. We are not able to compile/infer the conditional probability distribution and marginal probability distribution for this Bayesian network from the data/spreadsheet. The software libraries (gRain, bnlearn in CRAN) give up. As of now, we are trying to solve the problem by introducing some latent variables and by exploiting some local structure inherent in the problem. We are not successful so far. Any generic suggestions in terms of algorithms, tools to model, infer and solve problems of this scale shall be very useful.


title: Improving classification results

text: I have training data set with around 1500 positive set samples and 4500 negative set samples. All the features are numeric( floating or integer type values) and the data is specific to bio-informatics domain. I have tried using SVM and random forest, both give me an accuracy of around 88% on test set of size 700 and containing only positive samples. I have also tried LDA and adaboost in R, but they give worse results. I have tried feature selection which I observed improves the prediction accuracy by another 1.5%. Algorithm used for feature selection was random forest. I have two questions:   1. I need to improve my classification accuracy on the given test set.Please guide me how I can achieve the same.   2. Also, if there has to some step wise strategy that goes into selecting a classification algorithm. Thanks


title: Combining CATPCA with conventional PCA or factor analysis

text: I have some concerns regarding factor analysis and especially about combining the factor analysis for an ordinal scale (categorical data) - CATPCA with conventional PCA or with factor analysis. Basically, I need to enter my set of variables in CATPCA, take the transformed variables from CATPCA and introduce them in conventional PCA/factor analysis (because, for example, I need to rotate the factor to obtain a reasonable solution). Basically, I want to introduce first the variables in CATPCA because CATPCA does not require the data to be normally distributed and I also can replace the missing values with the mode (I have an ordinal scale so I think mode is better than mean). My concerns are: After introducing the transformed variables in conventional PCA/factor analysis, is it okay just to follow the steps for conventional PCA and not taking in consideration that the initial variables have been transformed (standardized, from what I have observed) in CATPCA? I mean, can I go through evaluating MSA, communalities, factor loadings and so on, as in a regular case? After obtaining the factors, if I want to assess the reliability, is it okay to take into consideration the initial variables (not the transformed ones in CATPCA which I also use in conventional PCA), assuming that I would like to generate the summated scales (for each factor) based on the initial variables and further use the summated scales for different constructs for correlation and regression? If I have a certain number of constructs and I would like to use the same methodology for each construct, is it okay to use CATPCA + conventional PCA/factor analysis even for the constructs for which I do not need to rotate the solution? (I want to use a combination of CATPCA+PCA for a number of reasons, including to rotate the solution, but some of the constructs I am analyzing are unidimensional or, even I they are multidimensional, the solution may be okay without rotation) Is it accepted to replace missing values with the mode in CATPCA while replacing the missing values with the mean (for the initial variables, not the transformed ones) for summated scales, correlation and regression? (This would also imply that I would use the initial variables and not the transformed ones for calculating summated scales, for correlation and for regression)


title: Non-Stationarity, first differences and Panel Data

text: I have build a sentiment index and am now validating its statistical properties and significance. I stumbled upon two problems. (1) `dfgls`indicates that my sentiment is non-stationary as the t-statistics are all above (less negativ) than the critical value. Accordingly, I use the first difference. My sentiment now granger-causes the (also non-stationary) market index, which is predicted by theory. But after taking the first difference it is not able to explain movements in the index using univariate regression, which it did before. Do I have to accept this as it is or do you distinguish between the applicability of first differences in the different tests? For example, would it be possible to use first differences for the granger-causality test and an estimator which takes into account the non-stationarity for the univariate analysis? (2) I want to apply my sentiment in panel data. Here, before taking the first difference, it applies very well. Afterwards, it does not. I also checked my dependend variable for unit root using `xtunitroot`and all p-values are 0.000. Therefore, I would judge that the data all panels are stationary. When I use `dfgls` for seperate panel IDs though, the tests imply non-stationarity for certain lags. Do I misunderstand the test? Thanks for any advice.


title: Stata: ml versus nl

text: I need to estimate an equation using some iterative method to optimally choose the coefficients. I won't go into the reasons that I can't use a simpler method like OLS, GLS or instrumental variables (it's a long, sad story). My question is about whether I can expect equivalent results from `ml` and `nl`. Say I want to estimate $y_t = \alpha + \beta y_{t-1} + \gamma x_t + u_t$ where I have autocorrelation in the error terms and they are correlated with $y_{t-1}$. Do I need to put in the effort to specify a likelihood function and ask Stata to maximize it, or will I get the same result if I ask Stata to minimize the sum of squares using `nl`?


title: Statistical analysis website postgresql

text: I'm currently working with postgresql on an amazon server ec2 and was wondering what is a good language to represent statistical analysis? I've been told that R is the best but heard from one experienced programmer that R runs slow when trying to display it on a website. I'm trying to perform some statistics with the help of my psql database and hopefully display some data on a webpage.


title: What is a "strictly positive distribution"?

text: I am reading Judea Pearl's "Causality" (second edition 2009) and in section 1.1.5 Conditional Independence and Graphoids, he states: > The following is a (partial) list of properties satisfied by the conditional > independence relation (X_||_Y | Z). > >   * **Symmetry:** (X_||_ Y | Z) ==> (Y_||_X | Z). >   * **Decomposition:** (X_||_ YW | Z) ==> (X_||_Y | Z). >   * **Weak union:** (X_||_ YW | Z) ==> (X_||_Y | ZW). >   * **Contraction:** (X_||_ Y | Z) & (X_||_ W | ZY) ==> (X_||_ YW | Z). >   * **Intersection:** (X_||_ W | ZY) & (X_||_ Y | ZW) (X_||_ YW | Z). > > > (Intersection is valid in **strictly positive probability distributions**.) > > (formula (1.28) given earlier in the publicatiob: [(X_||_ Y | Z) iff P(X | > Y,Z ) = P(X | Z) ) But what is an "strictly positive distribution" in general terms, and what distinguishes a "strictly positive distribution" form a distribution that is not strictly positive?


title: Evidential reasoning in Gaussian Bayesian Networks

text: I am working on Gaussian Bayesian Networks (GBN) i.e. the Bayesian Networks where all the random variables are continuous in nature. I am seriously trapped in the problem of evidential reasoning in case of GBNs. Can anybody help me by guiding me regarding   1. whether there are any established algorithm for evidential reasoning or inference in case of GBNs?   2. if yes then please provide me some references   3. if no then can anybody suggest how should I proceed in inference using GBNs.   4. can I directly apply the inference algorithms of discrete variable Bayesian Networks in case of GBNs?   5. lastly, my objective is to deduce some algorithm for multiply connected GBNs. in most of the cases the inference algorithms are for singly connected Bayesian Networks. For multiply connected Bayesian Networks what will be the changes required in inference algorithms? Thanks in advance.


title: Bootstrapping clusters in R

text: I am running a negative binomial regression of clinic counts in each county in the entire country (~3k counties). I'd like to at least partially account for the non-independence of neighboring counties by bootstrapping the confidence intervals in a "clustered" fashion--e.g. draw an entire state's (50 states total) worth of data at once. This has become standard practice, for better or for worse, in the econometric literature. I could write the code to do this myself, but the `boot` package seems like it should have the ability to do this somehow, and in general I prefer tested, general solutions to one-off hacks. Is there a way to coerce the `boot` package to do a clustered bootstrap? I tried the `strata` argument, but that randomizes _within_ strata rather than randomizing which cluster gets taken, as the following code confirms:               dat <- data.frame( cluster=rep(letters[1:5],each=10), x=runif(5*10), stringsAsFactors=TRUE )     boot.stat <- function(dat,idx) {         print(dat[idx,]$cluster)          print(table(dat[idx,]$cluster))         mean(dat[idx,]$x)         }         boot(           data=dat,           statistic=boot.stat,           strata=dat$cluster,          stype="i",          R=5      )


title: interrater reliability for multiple variables

text: i and another coder want to assess our agreement. our experiment involves coding the utterances of participants. there are 17 different codes that can be assigned, and we have independently counted their frequencies. participant 1 may have uttered code 1: 3 times, code 2: 2 times, code 3: 1 time, and so on until code 17. participant 2, ditto. on top of all this, i and the other coder want to compare the frequencies we found. so for participant 1, coder 1 and coder 2 might agree on the frequency of code 1, and not on code 3, and so on for participant 2 and 3 and 4. i don't know of a way to ascertain the overall interrater reliability between us- is there a way? and can it be done in SPSS? thanks for your time.


title: Should I use a paired or unpaired t-test for this problem?

text: I tested 5 objects with two different devices each, but on the last test for the second device, there was an error so the value was unusable (ie on objects 1-5, first device series was 1,1,2,2,3 and second series was 1,2,2,2,X). Should I use n=5 for the first device and n=4 for the second device and use an unpaired (heteroscedastic) test or should I throw out the good value for the last test of the first device and use a paired t-test with n=4 for both devices?


title: Survival analysis in SAS-is there a way to include a random effect with interval censored data?

text: I'm trying to use some form of survival analysis (e.g. accelerated failure time or proportional hazards) to study seed germination time. My data is interval censored (also called "grouped-time") and has a random effect (also called "frailty"). How can I do this, ideally in SAS? My seeds were on petri dishes that were censused on day 1, 3, 5, 8, and 14. These are interval censored, with a seed that germinated on d8 really germinated between day 6 and day 8. PROC LIFEREG can handle interval censoring, but as far as i know, it doesn't handle random effects. Seeds on the same petri dish are not independent, so I need to be able to include the dish as a random effect. I would like to make inferences about the mean or median time to germination for a seed grown under several treatment conditions (e.g. genotype, maternal temperature, seed temperature, & interactions) Is there an approach that allows interval censoring as well as random effects? If so, how should I modify the SAS code below to include random effects? Do I need to switch procedures to something like PROC NLMIXED? If this is something that cannot be done in SAS, what software is the best/easiest to do this? PROC LIFEREG DATA=survGeno2 plots=probplot; TITLE "Interval censored AFT"; CLASS geno; MODEL (start,fin)=geno; RUN; Thanks!


title: How to calculate zipf exponent of page degree distribution?

text: I want to calculate the exponent of the function which shows zipf's law in data. I'm working in the context of graph node degree analysis and want to comprehend results shown by Border et al. on a new dataset. I have pairs of values where the first value denotes the degree and the second one the number of nodes with that degree, like that:               1   188638380     2   20839145     3   9049892     4   5427599     [...]     641531  1     663136  1     820267  1     920489  1      I'm searching for a program that calculates the zipf exponent automatically and stumbled about zipfR of the R language. Unfortunately I'm unfamiliar with R and with statistic as well. Furthermore zipfR don't want to work with me. I needed to use `degrees<-data-frame(scan("path"), list(0,0),skip=1))`to get the data out of the file that's because I have about 13.000 data points. When I tried to use functions of zipfR I got errors that zipfR can't work with that. Afterwards I tried to adapt the algorithm shown when the question about frequency's shows up but the resulting exponent is far to small. There might be some reasons for it I'm not sure of(I tend to reason 2):   1. The code was written for single values and not for pairs. Did the algorithm assume iterative raising integers? Because in my dataset are gaps   2. Broder might truncated values. How to estimate a good point for truncation? To let you see what I've done:               dataset<-data.frame(scan("/data/dd_in.txt", list(0,0),skip=1))     p<-dataset/sum(dataset)     lzipf<-function(s,N) - s*log(1:N)-log(sum(1/(1:N)^s))     opt.f<-function(s) sum((log(p)-lzipf(s,920489))^2)     opt <- optimize(opt.f,c(0.5,4))      If I get it right, the values after c are the range for the exponent. I expect an value of around 2 for the exponent. The value in the opt.f function of 920.489 is the maximum value of the degrees. Is that correct? What about the number of pages with a certain degree? Where can I find them or add them? Is that implicitly done due to the p value or did I need to fill the gaps with 0?


title: Infinite p-value when checking normality of distribution

text: My data has following distribution: ![enter image description here](http://i.stack.imgur.com/Hgwys.png) I wish to determine if it is a Normal distribution, so:               > library(nortest)     Warning message:     package 'nortest' was built under R version 2.12.2      > sf.test(y)     Error in sf.test(y) : sample size must be between 5 and 5000     > ad.test(y)                 Anderson-Darling normality test     data:  y      A = 5487.108, p-value = Inf          > cvm.test(y)                Cramer-von Mises normality test     data:  y      W = 855.7627, p-value = Inf          > pearson.test(y)                 Pearson chi-square normality test         data:  y      P = 2456556, p-value < 2.2e-16          > qqnorm(y); abline(0, 1)      When I do `qqnorm`, I found ![enter image description here](http://i.stack.imgur.com/Sa9Ai.png) Can I conclude `y` does not have a normal distribution, since the p-value = INF and the qqnorm did not fit the abline? Should I expect all these p-values close to zero or one when my distribution is in normal distribution? How about the pearson.test? How should I interpret the numbers (A, W and P) before the p-values?


title: calculating the SE in a dice game

text: _I am trying to figure out why calculating the SE directly in 22-4 gives me a larger result than that given in the solution_ (examples change slightly each time the page loads): http://www.stat.berkeley.edu/~stark/SticiGui/Text/standardError.htm#SE_of_RV In 22-3, the method for directly calculating the SE seems to be $SE(X) = \sqrt{\sum \limits_{i=1}^{n}P(X=x_i)(x_i-E(X))^2} $, where you multiply each specific outcome by its probability, sum those values and then take the square root. The specific values I had for 22-4 were: You pay $\$2$ to roll a fair pair of dice. You win $\$3$ if you roll doubles. The SE of your net winnings after nine plays of the game is __? In 22-4, the solution calculates the SE of the binomial and then uses a linear transformation to get the standard error for playing with dollar amounts other than 0 and 1: $Y= -18+3X$ the linear relationship of net winnings for 9 plays, where X is the number of wins in 9 plays. $SE(X)=\sqrt{np(1-p)}=\sqrt{9\frac{1}{6}\frac{5}{6}}= 1.118$ The SE of a binomial distribution for 9 plays. $SE(Y)=|3|SE(X)=\color{blue}{3.354}$ The SE for 9 plays if you win $\$3$ by rolling doubles and lose $\$2$ otherwise. I tried to use the method I used for 22-3 (multiply the squared deviation of each outcome by its probability, sum those values and take the square root to get the SE), however I'm not sure why this doesn't work. $E(X) = x_0P(X=x_0)+x_1P(X=x_1)= -2(5/6)+3(1/6)=-7/6$ $SE(X) = \sqrt{P(X=x_0)(x_0 - E(X))^2+P(X=x_1)(x_1 - E(X))^2}$ $=\sqrt{\frac{5}{6}(-2+\frac{7}{6})^2+\frac{1}{6}(3+\frac{7}{6})^2}= 1.863$ for nine plays $\sqrt{9}SE(X)=\color{red}{5.5902}$


title: Problem with integrating exponential function in R

text: I'm trying to calculate the finite integral for the CDF of the normal distribution, when I plug the equation into Wolfram Alpha and integrate > exp((-t^2)/2) dt from -inf to 1000000 I get: 2.5 (for just the integral), when I run the same equation through R I get 0. On further digging I think the discrepency lies with the 'exp' function. I'm guessing that means I'm doing something wrong or I missed a step somewhere. Any help would be greatly appreciated.


title: Non-conjugate improper uniform priors in Bayesian data analysis: how to handle infinite sums?

text: I've been working on Griddy Gibbs sampler (paper: Ritter and Tanner) and I've implemented it in R. But I've faced a problem when I started thinking on its uses in other contexts. If I try to use an impoper uniform prior in a non conjugate setting, I'm supposed to perform numerical integration with the lower and upper bounds of the priors. When using flat priors for continuous values, the upper and lower values are simply minus and plus infinity. Thus, it is not possible to have a grid approximation to this interval. Are there any methods to handle this configuration? So far, I've been happy with the results of my sampler, but this feels like a situation that I can't handle with my current method.


title: Probability of pairwise difference of samples from distribution with finite support

text: I'd appreciate any help on the following problem: Let $X_1, X_2, \dots, X_N$ be i.i.d. continuous random variables with support $[0, 1]$. What is a reasonable bound on the probability that some pair of random variables is less than $\epsilon$ apart? Mathematically, I'd like to find some $\delta$ not too large such that $$\Pr \left( \exists_{i,j:\ i\ne j}\ |X_i - X_j| < \epsilon \right) < \delta.$$ If it helps, I'm not interested in the general case -- I'm only interested in the cases of a uniform distribution, a truncated Gaussian distribution, and a beta distribution. Again, I only need a reasonable probability bound, not the exact probability. A friend suggested that the Dirichlet distribution could possibly be helpful here, but I haven't seen a way to apply it here. Thank you in advance for help!


title: Regression analysis for more than one categorical variable in time series

text: I have a time series data for shipment with following variables:   1. Year: 2008, 2009, 2010, 2011, 2012, 2013   2. Month: jan, feb, ..., dec   3. Number of ordering days   4. Shipment Volume I want to know the effect of the number of ordering days on the shipment volume. Can this can be done using regression method? If so, should I consider month and year as categorical variables?


title: Proving the Accept-Reject Algorithm

text: Prove the following is equivalent to the Accept - Reject (AR) algorithm.   1. Generate $X \sim g$;   2. Generate $U|X = x \sim U_{[0; Mg(x)]}$;   3. Accept $Y = X$ if $U \le f(x)$   4. Return to 1. otherwise. Attempt: To prove the equivalency, I must prove that $f(x)$ is the same as $\frac{f(x)}{cg(x)}$ for the uniform distribution. If I can show that, I can then simply use the traditional AR proof (http://www.columbia.edu/~ks20/4703-Sigman/4703-07-Notes-ARM.pdf). I am stuck proving that $f(x) = \frac{f(x)}{cg(x)}$. Also, this is proof is equivalent to proving that $P(Y < y) = P(X < x | U < f(x))$. I think if I proved that, I could proove that the above algorithm is equivalent to the AR algorithm.


title: On univariate outlier tests (or: Dixon Q versus Grubbs)

text: In (most of) the analytical chemistry literature, the standard test for detecting outliers in univariate data (e.g. a sequence of measurements of some parameter) is Dixon's Q test. Invariably, all the procedures listed in the textbooks have you compute some quantity from the data to be compared with a tabular value. By hand, this is not much of a concern; however I am planning to write a computer program for Dixon Q, and just caching values strikes me as inelegant. Which brings me to my first question:   1. How are the tabular values for Dixon Q generated? Now, I have already looked into this article, but I'm of the feeling that this is a bit of cheating, in that the author merely constructs a spline that passes through the tabular values generated by Dixon. I have the feeling that a special function (e.g. error function or incomplete beta/gamma) will be needed somewhere, but at least I have algorithms for those. Now for my second question: ISO seems to be slowly recommending Grubbs's test over Dixon Q nowadays, but judging from the textbooks it has yet to catch on. This on the other hand was relatively easy to implement since it only involves computing the inverse of the CDF of Student t. Now for my second question:   1. Why would I want to use Grubbs's instead of Dixon's? On the obvious front in my case, the algorithm is "neater", but I suspect there are deeper reasons. Can anyone care to enlighten me?


title: Implementing a Metropolis Hastings Algorithm in R

text: Consider a univariate normal model with mean $µ$ and variance $τ$ . Suppose we use a Beta(2,2) prior for $µ$ (somehow we know µ is between zero and one) and a $log-normal(1,10)$ prior for $τ$ (recall that if a random variable $X$ is $log-normal(m,v)$ then $log X$ is $N(m,v))$. Assume a priori that $µ$ and $τ$ are independent. Implement a Metropolis-Hastings algorithm to evaluate the posterior distribution of $µ$ and $τ$ . Remember that you have to jointly accept or reject $µ$ and $τ$. Also compute the posterior probability that $µ$ is bigger than 0.5. Here are the data:               2.3656491 2.4952035 1.0837817 0.7586751 0.8780483 1.2765341     1.4598699 0.1801679 -1.0093589 1.4870201 -0.1193149 0.2578262      Attempt: (Heres to enter my data)               mu <- rbeta(1,2,2)     tau <- rlnorm(1,1,10)     x <-        c(2.3656491, 2.4952035, 1.0837817, 0.7586751, 0.8780483, 1.2765341,       1.4598699, 0.1801679, -1.0093589, 1.4870201, -0.1193149, 0.2578262)      This is a MH algorithm found online. I am not sure how to apply the above data to it.               # starting point     a0 <- -5      b0 <- -10           # length of chain     nit <- 1000     a <- rep(0,nit)     b <- rep(0,nit)          # initialize     a[1] <- a0     b[1] <- b0          # tuning parameter     s0 <- 2.0 # maximum step size in random walk proposal     # function. try different s0, e.g., 0.1, 1.0, 2.0          # start chain     counter = 0 # monitor number of acceptances.     for( i in 1:nit) {       s <- s0*runif(1)       theta<-2*3.1415926*runif(1)       anew <- a[i] + s*cos(theta) # random walk       bnew <- b[i] + 2*s*sin(theta) # random walk        r <- pdf(anew,bnew)/pdf(a[i],b[i])# acceptance ratio       test <- runif(1)       if(test < r ) # accept proposed moved.       {               a[i+1] <- anew          b[i+1] <- bnew         counter = counter + 1;       }       else # reject proposed move, stay put.       {         a[i+1] <- a[i]         b[i+1] <- b[i]       }     }     # acceptance rate =     counter /nit


title: Time series forecasting with CAAR model

text: Is there a package for Conditional autoregressive models in R (CARR). I need to use it for stock market forecasting and to compare with GARCH.


title: Interpreting a PCA biplot

text: I have a bilot of PCA and am not able to correctly interpret it. There are two type of whiskies and both have two unmatured treatments and two matured treatments, and compounds relating to each treatment are shown with their concentration. If you can see the snapshot of the data and plot below, I can see that some compounds like Furfural and Syringaldehyde seem to be related to each other. What I can't get is, for example, whisky1unmat2 has high furfural according to table. But in the biplot whisky1unmat2 and furfural are totally opposite to each other in different quadrants. Could you please tell me how to decipher this plot? ![pca biplot and data](http://i.stack.imgur.com/fjAOX.jpg)


title: Useful statistical functions for business -- for use by a newbie

text: I am a database developer working for a sales and manufacturing business. I am mostly ignorant about statistics. We need useful metrics. Our managers are tuned to accounting, and the vagaries of production tend to confound us. We do very little measuring of our production, and what we have is poorly formed. I should note we are a "job shop", not a "flow shop" -- we do a lot of engineer- to-order work, so the usual MRP standards are often hard to apply. Some tradition business metrics are known to us; for example "inventory turn- over rate." But we are unable to convert those to useful information. I believe our inability to qualify data statistically is a big reason why. Of course we perform averaging all the time. Rolling averages (smoothing data over time using a 3-week rolling average, for example) is a helpful extension. Recently I discovered how to apply standard deviation to labor costing with wonderful benefits. Now that I understand (A) averaging, (B) rolling averages, and (C) standard deviation, what are the next useful functions or functions I should seek to learn? I would love to have your insights on "business intelligence," I mean defining and using metrics. But it's the use of stats to get from raw data to usable information that I'm really after. No matter, give me whatever you've got.


title: Variance estimation from two gaussian distributions

text: Assume a stochastic process with observation $r$, and two hypotheses : $X \sim (0, \sigma^2)$ and $Y \sim (0, \sigma^2 + \tau^2)$. When we observe/receive $r$ we don't know which hypothesis $X$ or $Y$ it came from. Is there any way to have an estimate $\hat{\sigma}^2$ for the variance $\sigma^2$ under these conditions? We know that the variance in one distribution can be estimated as $\hat{\sigma}^2 = \frac{1}{N} \sum (x_i - \hat{\mu})^2$. But how to estimate $\hat{\sigma}^2$ when it is involved in two distributions?


title: Difference between Kernel classifier and linear classifier

text: I would just like to know what are the differences between kernel classifier and linear classifier? In what kind of problems the first is used and in what kind the second? What could be the advantages of the one over the other?


title: How to project test data onto the principal vectors produced by kpca?

text: So I have a dataset of dimension 6395x15 and when I apply kernel-PCA on this dataset, I get a rotated matrix of dimesion 6395x596. pcv() gives me the matrix of column-wise eigen-vectors and this matrix is also of dimension 6395x596. My question is how can I project a new data or my test data onto these eigen-vectors. My new data will be a vector of size 15 and I would want to project this onto the 596 dimensions produced by kpca so that I can test my classifier model build on a subset of 596 dimensions.


title: How are partial regression slopes calculated in multiple regression?

text: I'm trying to understand how multiple regression statistically controls for the effects of other predictor variables when calculating partial regression slopes. In a multiple regression of Y~X1+X2, would the partial regression slope of X1 be given by [Y]~[residuals of X1~X2], or by [residuals of Y~X2] ~ [residuals of X1~X2]? Different pages on the internet tell me different things. I've done some simulations to try and figure this out (see below), and it seems that both methods give the same estimates of slopes as multiple regression, but only the latter method has similar standard errors around those estimates. This makes me think the latter method is the one that multiple regression uses, but it would be nice to know for sure. Similarly, if I wanted to plot Y against X1 so that I could visualise how strongly the two were related while also controlling for any confounding with X2, would I plot [Y]~[residuals of X1~X2], or [residuals of Y~X2] ~ [residuals of X1~X2]? These two plots in the code below look very different in terms of the strength of the relationship. Thanks for your help, Jay               #1. simulate data, where x1 and x2 are correlated due to lurking variable,      #...and y is explained by both.     lurker <- rnorm(n=100)     x1 <- rnorm(n=100, mean=lurker*2, sd=1)     x2 <- rnorm(n=100, mean=lurker*5, sd=1)     y <- rnorm(n=100, mean=x1*2 + x2*5, sd=1)          #2. multiple regn model to estimate partial slopes:     summary(lm(y~x1+x2))      #partial slopes pretty close to simulated values          #3. calculate partial slopes manually, using either      #....(1) Y~[resids of X1~X2] OR (2) [resids of Y~X2]~[resids of X1~X2]            #3.a. based on (1) Y~[resids of X1~X2]       m.x1x2 <- lm(x1~x2)       resids.x1x2 <- m.x1x2$residuals       summary(lm(y ~resids.x1x2))    #slope pretty close to true value, but conf intervals MUCH larger than those for MR estimate            #3.b. based on (2) [resids of Y ~ X2]~[resids of X1~X2]       m.y1x2 <- lm(y~x2)       resids.y1x2 <- m.y1x2$residuals       summary(lm(resids.y1x2 ~resids.x1x2))  #also very close to true value, but conf intervals now similar scale to those from MR.            # plot the relationship between y and x1 after controlling for x2, based on the different methods:         op <- par(mfrow=c(2,1))         plot(y ~resids.x1x2)         plot(resids.y1x2 ~resids.x1x2)         par(op)


title: What, exactly, does the "coordinate gradient descent algorithm" do?

text: I knew about "Coordinate descent" and "Gradient descent" algorithms before; they are well-known and Wikipedia has articles for those. However I recently came across an algorithm called "Coordinate gradient descent". How does this algorithm exactly work, and how does it relate to "Coordinate descent" and "Gradient descent"? An example for usage of this term can be found on slide 7 at the link http://courses.cs.washington.edu/courses/cse546/12wi/slides/cse546wi12ClusteringEM.pdf. I think what K-means or EM does is "Coordinate gradient descent" (or "Coordinate gradient ascent", based on how you view the optimization), but I could not figure out the intuition behind the name.


title: Feature selection using mutual information in Matlab

text: I am trying to apply the idea of mutual information to feature selection, as described in these lecture notes (on page 5). My platform is Matlab. One problem I find when computing mutual information from empirical data is that the number is always biased upwards. I found about 3~4 different files to calculate MI on Matlab Central and they all give big numbers (like > 0.4) when I feed in independent random variables. I am not an expert, but the problem seems to be that if you simply use joint and marginal densities to compute MI, bias is introduced in the process because MI is by definition positive. Does anyone have practical advice on how to estimate mutual information accurately? A related question is, in practice, how do people actually use MI to select features? It is not obvious to me how to come up with a threshold value since MI is in theory unbounded. Or do people just rank the features by MI and take the top k features?


title: $ARIMA(p,d,q)+X_t$, Simulation over Forecasting period

text: I have time series data and I used an $ARIMA(p,d,q)+X_t$ as the model to fit the data. The $X_t$ is an indicator random variable that is either 0 (when I don’t see a rare event) or 1 (when I see the rare event). Based on previous observations that I have for $X_t$ , I can develop a model for $X_t$ using Variable Length Markov Chain methodology. This enables me to simulate the $X_t$ over the forecasting period and gives a sequence of zeros and ones. Since this is a rare event, I will not see $X_t=1$ often. I can forecast and obtain the prediction intervals based on the simulated values for $X_t$. Question: How can I develop an efficient simulation procedure to take into account the occurrence of 1’s in the simulated $X_t$ over the forecasting period? I need to obtain the mean and the forecasting intervals. The probability of observing 1 is too small for me to think that the regular Monte Carlo simulation will work well in this case. Maybe I can use “importance sampling”, but I am not sure exactly how. Thank you.


title: Trying to compute Gini index on StackOverflow reputation distribution?

text: I'm trying to compute the Gini index on the SO reputation distribution using SO Data Explorer. The equation I'm trying to implement is this: $$ G(S)=\frac{1}{n-1}\left(n+1-2\left(\frac{\sum^n_{i=1}(n+1-i)y_i}{\sum^n_{i=1}y_i}\right)\right) $$ Where: $n$ = number of users on the site; $i$ = user serial id (1 - 1,225,000); $y_i$ = reputation of user $i$. This is how I implemented it (copied from here):               DECLARE @numUsers int     SELECT @numUsers = COUNT(*) FROM Users     DECLARE @totalRep float     SELECT @totalRep = SUM(Users.Reputation) FROM Users     DECLARE @giniNominator float     SELECT @giniNominator = SUM( (@numUsers + 1 - CAST(Users.Id as Float)) *                                    CAST(Users.Reputation as Float)) FROM Users     DECLARE @giniCalc float     SELECT @giniCalc = (@numUsers + 1 - 2*(@giniNominator / @totalRep)) / @numUsers     SELECT @giniCalc      My result is (currently) -0.53, but it makes no sense: I'm not sure even how it could have become negative, and even in abs value, I would have expected the inequality to be much closer to 1, given how reputation grows the more you have it. Am I unknowingly ignoring some assumption about the distribution of the reputation/users? What do I do wrong?


title: Can co-integration be conducted if one of the variables is not of the same order?

text: All my independent variables are stationary at their levels except one variable which became stationary after first difference i.e. my dependent variable. Can I still conduct co-integration test ?


title: Accuracy of classifiers with Adaboost

text: Does Adaboost ensure that resultant accuracy is more than or at least equal to current accuracies? What happens if **Classifier A** performs badly and the weights are accordingly updated and the next **Classifier B** performs very well (better than **Classifier A** ) on the data. How does **Adaboost** handle this problem?


title: Covariance of a linear and quadratic form of a multivariate normal

text: Does anyone know of an explicit matrix expression for the covariance of a linear and quadratic form? That is, $\mathrm{Cov}[\mathbf{a' y},\mathbf{y' Hy}]$ where $\mathbf{y}\sim \mathcal N(\boldsymbol{\mu},\boldsymbol{\Sigma})$. I'm particularly interested in the case where $\boldsymbol{\mu}=\mathbf{0}$, and I think this simplifies (without the normal assumption) to $\mathbb E[(\mathbf{a'y})(\mathbf{y'Hy})]$. Since this involves cubic terms it probably isn't going to be simple.


title: How to determine variables that make certain cluster different from others?

text: Say clustering was performed. My concern is to understand what characterizes certain cluster. Which variables differ most for a certain cluster? Are there any methods for performing such kind of analysis? For example, say I have 50 variables and 10 clusters. The desired result would sound like: for a cluster 1 variables 5, 16, 23, 42 and 49 contribute to the biggest difference from other clusters.


title: How to calculate MAE and RMSE when using SVM light and libsvm?

text: Can SVM-LIGHT generate MAE (mean absolute error) and RMSE (root mean square error) for svm's model performance?If we using libsvm,can libsvm generate MAE and RMSE too by using windows command?


title: Alternative measures of prediction power?

text: What are alternative measures of the prediction power of a model, apart from the coefficient of determination $R^2$? What are their strengths and weaknesses, especially in comparison to the $R^2$?


title: Effect of treatment on correlation between two sets of variables

text: Imagine I wish to see the effect of some treatment on the correlation between two groups of variables, which I measure with Spearman's rank correlation coefficient. Treatment A gives a Spearman's rank of 0.4, 0.44, 0.43 for my 3 replicates. Treatment B gives a Spearman's rank of 0.48, 0.45, 0.46 for my 3 replicates Treatment C gives a Spearman's rank of 0.47, 0.50, 0.51 for my 3 replicates. Can I compare the groups using e.g. the Mann-Whitney U test? N.B. I'm not looking for the correlation between the sets of variables, I'm looking to see whether a specific treatment improves the correlation. Apologies if this question is similar to this Significance test on the difference of Spearman's correlation coefficient


title: Dixon test for outlier but which one is an outlier?

text: Following my previous question, I used Dixon test for outliers with the help of Michael Chernick answer. So now I have pvalues for say 10 numbers (basically 10 patients). But I have around 50k pvalues and I want to know that which patient has an outlier. Example-               chr     start.position    stop.pos     pvalue    pat1     pat2  pat3 ...      1              10         110            0.03      3     3      7      2              100        200            0.006   0.3     4      5               So I have around 50k rows in this format. So my first step is, I will take all the rows which have significant pvalues (lets say, less than 0.05). Ok, now I know which positions contain outliers but my next step is, which patient is an outlier? Like in aforementioned example, I can say that for first row (first position), patient3 is an outlier and in second row, I can see patient1 is an outlier.For less number of positions, I can easily see the outlier but what if I have so many rows? Is there any way to find this thing? I am also ready for applying any other test or any other procedure to detect this. Thanks in advance.


title: What are the different types of codings available for categorical variables (in R) and when would you use them?

text: If you fit a linear model or a mixed model there are different types of codings available to transform a categorical or nominal varibale into a number of variables for which paramaters are estimated, such as dummy conding (the R default) and effects coding. I heard that effects coding (sometimes called deviation or contrast coding) is preferred when you have interactions, but what are the possible contrasts and when would you use which type of contrast? The context is mixed modeling in R using `lme4`, but I think broader responses are fine. Sorry, if I missed a similar question. EDIT: Two helpful links are: effect coding & dummy coding explained.


title: How to calculate Mahalanobis distance in SPSS for an exploratory factor analysis?

text: I have a question regarding data screening for an exploratory factor analysis (EFA). I am conducting an EFA to identify the factor structure of 20 questions that I created on the topic of spirituality. I want to identify outliers in my sample using mahalanobis distances, and I am doing this on SPSS using a linear regression (`Analyze -> Regression -> Linear`).   * I entered the 20 questions in SPSS as the "Independents", but what would variable should be entered under the "Dependent" category?    * Or, if I can't use SPSS to find the mahalanobis distances, is there another (easy!!) way to find the m distances?


title: Showing $\mathbb{E}[T_n] = \theta \mathbb{E}_1[T_n]$ is scale equivariant?

text: This is question 5 is from Staudte and Sheather (1990), Robust estimation and testing. Let $X_1,\ldots , X_n$ be i.i.d with $$ F_\theta = F(\frac{x}{\theta}),\quad x>0;\theta>0.$$ Assume that $T_n = T_n(X_1,\ldots ,X_n)$ is scale equivariant. Show that $$\mathbb{E}[T_n] = \theta \mathbb{E}_1[T_n]$$. Using $\int_a^b f(u) du = \sum_{k=a}^b$ and $X\rightarrow \theta X$, my attempt at the question is as follows. \begin{aligned} \mathbb{E}[T_n] &= \int T_n dF(X_n) \\ &= \sum_x T_n P(X) &= \frac{1}{n}[\theta X_1 +,\ldots ,+ \theta X_n] \\ &= \theta\frac{1}{n}[ X_1 +,\ldots ,+ X_n] \\ &= \theta \sum_x T_n P(X) \\ &=\theta \int T_n dF(X_n) \\ &= \theta \mathbb{E}[T_n] \end{aligned} Is my working out correct? If not, where and why was I wrong? I should also point out that this is not an assignment question. I need to study this book to gain some background knowledge in robust statistics.


title: More interpretable measure of association than odds ratios for contingency tables with 0 counts

text: How strong/weak of a correlation is this?                                 Exposure                        No Yes  Total                 FALSE 139 467   606     Disease     TRUE   11 104   115                 Total 150 571   721          OR =  2.81      Exact 95% CI =  1.45, 5.97       Chi-squared = 10.49, 1 d.f., P value = 0.001     Fisher's exact test (2-sided) P value = 0.001       How about this?                                  Exposure                         No  Yes  Total     Disease     FALSE  301 1232  1533                 TRUE     0   17    17                 Total  301 1249  1550          OR =  Inf      Exact 95% CI =  1, Inf       Chi-squared = 4.14, 1 d.f., P value = 0.042     Fisher's exact test (2-sided) P value = 0.057        * Are there any better/other methods for calculating relative risk with a confidence interval for data like these?    * If there are other methods, what are their advantages/disadvantages? So, the first one is pretty easy, we could say "we estimate that those with the exposure are 2.81 [1.45, 5.97] times more likely to have the disease (p=.001)". That sounds like a pretty strong finding. Surprisingly strong for the hypothesis I'm working with. The second table doesn't work with the same formula, "we estimate that those with the exposure are infinitely more likely to have the disease [1, Inf], p < .05". What can be said about the second table? Shrug and say, "hrm, look at that, we are missing about 4 not exposed with the disease .. and there are none who have the disease who weren't exposed." I guess I already know the answers here. But, on the second table, I'm torn between thinking, "OMG, LOOK, UNEXPOSED HAVE 0, and p<.05!" and thinking, "that's interesting that that happened completely by chance". I guess that's when we say, "more research needed"?


title: Books on designing clinical trials for psychotherapy research

text: I'm designing a clinical trial that will evaluate a behavioral therapy against a control condition, and I'm looking for some good texts on this topic. There are a lot of books on clinical trials in medicine, but I've had a hard time finding one that also covers the topics specific to psychotherapy research. For instance, topics like choosing an appropriate control group and evaluating the therapist's adherence and competence are different in psychotherapy research than in drug trials. Any tips on a comprehensive book on this topic would be welcome.


title: Can GEE handle small unaccounted clusters in data?

text: A survey is being administered in a resource intensive setting. A 2-phase study design will be implemented using cases and controls identified by a certain behavior (say, smoking) and they will be screened for SNP information on a number of biomarkers for disease associated with the negative health outcomes associated with that behavior. The problem is that in this setting there are potentially large numbers of first degree relatives living in the same household, which would contribute to environmental factors influencing prevalence of the behavior or outcome risk independent of behavior (say, second hand smoking). Due to confidentiality reasons, it was decided not to collect residence. Assuming we could capture auxiliary exposure independent of behavior (such as "how many cigarettes are smoked in the same room as you in a week's time?), this could account for confounding by residence. Then relatedness would only be a unaccounted for cluster. How efficient is a GEE with independent covariance structure at accounting for small numbers of unaccounted for clusters in the data? If we use these data, can we come up with a recommended sample size difference accounting for and not accounting for relatedness to achieve a target power for determining association with disease status and behavior? I appreciate references to scholarly publications directly addressing this issue from related applied problems.


title: multicollinearity test

text: I use several control variables and   use independent variable $A$ for model 1,   use independent variable $B=1-A$ for model 2,   use independent variable $C=A-E-F$ for model 3,   use independent variable $D=A-C-F$ for model 4,   use independent variable $F=A-C-D$ for model 5. I have to report correlation coeffiecient matrix.   I wonder is it okay to report just one matrix with control variables and independent variable $A$? (Since other independent variables are just linear transformation of $A$)


title: Mixture model as a prior distribution

text: I've just started working with Bayesian models. My question is in the context of hierarchical Bayesian model. Suppose you have _n_ models to train. However, some of these models are similar to each other in some sense and I could assume that they their parameters are drawn from the same prior distribution. Another subgroup of that models are also similar to each other and their parameters share the same prior distribution. And it repeats for all subgroups (unknown in advance). So, my question is: Is it possible to use a mixture of distributions as a prior distribution for all these _n_ models? In such hierarchical model each component of the mixture model respond to a group of models. The underlying idea is to cluster these _n_ models by employing a mixture model as a prior distribution. Does it make sense? Thanks!


title: Why does one bootstrap sample help reduce bias in a biased variance estimator?

text: I am looking a bit on the effect of bootstrap on bias. I came by the following example in the estimation of variance (when dividing by n instead of n-1). Can someone explain to me why the bias when using the mean of a bootstrap sample comes out smaller than when using the simple mean? R code:               set.seed(3243)     n_sample <- 3     R <- 10000     var_values <- numeric(R)     var_with_mu_values <- numeric(R)     var_with_other_X_values <- numeric(R)     var_with_BS_values <- numeric(R)          var_with_mu <- function(x, mu=0) {         n <- length(x)         sum((x-mu)**2)/(n)     }          for(i in seq_len(R))     {         x <- rnorm(n_sample)              var_values[i] <- var_with_mu(x, mean(x))              var_with_mu_values[i] <- var_with_mu(x, 0)              other_x_mean <- mean(rnorm(n_sample))         var_with_other_X_values[i] <- var_with_mu(x, other_x_mean)              BS_mean <- mean(sample(x, replace = TRUE))         var_with_BS_values[i] <- var_with_mu(x, BS_mean)     }          # bias?     abs(mean(var_values)-1)     abs(mean(var_with_mu_values)-1)     abs(mean(var_with_other_X_values)-1)     abs(mean(var_with_BS_values)-1)      Results:               > abs(mean(var_values)-1) # known bias:     [1] 0.3422001     > abs(mean(var_with_mu_values)-1) # without bias:     [1] 0.007716617     > abs(mean(var_with_other_X_values)-1)     [1] 0.3199431     > abs(mean(var_with_BS_values)-1) ## Why is this bias smaller??     [1] 0.1202472        Thanks.


title: Cluster Big Data in R and Is Sampling Relevant?

text: I'm new to data science and have a problem finding clusters in a data set with 200,000 rows and 50 columns in R. Since the data have both numeric and nominal variables, methods like K-means which uses Euclidean distance measure doesn't seems like an appropriate choice. So I turn to PAM, agnes and hclust which accepts a distance matrix as input. The daisy method can work on mixed-type data but the distance matrix is just too big: 200,000 times 200,000 is much larger than 2^31-1 (the vector length limit before R 3.0.0.) The new R 3.0.0 released yesterday supports long vectors with length longer than 2^31-1. But a double matrix of 200,000 by 200,000 requires a continuous RAM larger than 16Gb which is not possible on my machine. I read about parallel computing and bigmemory package and am not sure if they are going to help: if I'm using daisy, it will generates a big matrix which cannot fit in memory anyway. I also read about the post about sampling: Is sampling relevant in the time of 'big data'? So in my case, is it relevant to use sampling on the data set, cluster on the sample and then infer the structure of the whole data set? Can you please give me some suggestion? Thank you! About my machine: R version 3.0.0 (2013-04-03) Platform: x86_64-w64-mingw32/x64 (64-bit) OS: Windows 7 64bit RAM: 16.0GB


title: Calculate the expected number of 95% confidence interval of binomial distribution

text: Can anyone show me how to calculate the expected number of 95% confidence interval of a binomial distribution using R, such as $\text{Bin}(100,0.5)$.


title: On Elo++ updating rule

text: I am not sure if this is the correct place to ask this kind of question, I hope it is. I am studying this paper on an improvement of Elo rating system called Elo++. On page 4 the author states that he want minimize a function called _total loss_ that is $$L=\sum_{i,j\in T}w_{i j}\left( \hat o_{i j}-o_{i j}\right)^2+\lambda\sum_{i\in D}\left(r_i-a_i\right)^2$$ where $r_i$ is the rank of the player that as to be estimated, $ \hat o_{i j}$ is the outcome of a game where player $r_i$ plays against player $r_j$, $a_i$ is the arithmetic average of the set $D$ of the player $r_j$ (the players $r_j$ which have played against $r_i$), $o_{i j}$ is the predicted outcome that is a function of $r_i, r_j$. The problem for me comes at page 5 where the author says that in order to minimize the total loss, he had used a stochastic gradient descent. The updating formula in the stochastic gradient descent is $$ r_i \leftarrow r_i-\eta\left(w_{i j}\left(\hat o_{i j}-o_{i j}\right)\hat o_{i j}\left(1-\hat o_{i j}\right)+\frac{\lambda}{N_i} \left(r_i- a_i\right)\right) $$ How is this updating formula is related to the problem of minimizing the first condition?


title: Attainable correlations for lognormal random variables

text: Consider the lognormal random variables $X_1$ and $X_2$ with $\log(X_1)\sim \mathcal{N}(0,1)$, and $\log(X_2)\sim \mathcal{N}(0,\sigma^2)$. I'm trying to calculate $\rho_{\max}$ and $\rho_{\min}$ for $\rho (X_1,X_2)$. One step in the given solution I have is: $\rho_{\max}=\rho (\exp(Z),\exp(\sigma Z))$ and $\rho_{\min}=\rho (\exp(Z),\exp(-\sigma Z))$, but they've made some references to comonotonicity and countercomonotonicity. I was hoping someone help me understand how they're relevant. (I know how to get this from the general expression but want to know specifically what the comonotonicity parts were saying.)


title: Does the probability of multiple independent events follow a normal distribution?

text: I'm heavily into role-playing systems, scripted a set of utilities in ruby for pen-and-paper games and I sort of understood statistics when I took it, but I could never for the life of me figure out the following: Given a varying-length series of independent 50/50 probabilities, i.e. 2, 3, 5 or 7 coin tosses, would the total number of successes (heads) have the same probability distribution as a fair roll of a die with number of sides equivalent to the number of coins + 1 (given that coins can produce all failures -- a zero)? This is a comparison, say, between the d20 system used by DnD and the multiple success-check dice used by Abberant. So in other words, would rolling a 20-sided die and flipping 20 coins and counting the heads result in the same probability distribution? I would be interested in a link that explains and visually represents the difference between distributions of independent events and a single event. EDIT: I was able to satisfy my basic curiosity with the above script, modified. No hard mathematical answer, but I'm not well-versed in this area anyway. Figured it was worth learning a bit though.


title: How to prove statistical significance between the 4 seasons

text: I have data regarding the number of deaths in a city over 35 years. I've collated the data, and now wish to prove the difference in total number of deaths between seasons is significant. Spring has 631 deaths, summer 540 deaths, autumn 502 deaths and winter has 605 deaths. Clearly there is a big difference between spring and autumn, but how do I prove that? Is it ANOVA? Or maybe t-test? I'm totally lost... Cheers


title: Patterns in binary sequences

text: I have data that looks like this:               011100111110100111     111111111111110010     111100001111000011      1D lanes of streams of data. Each row signifies the presence of a thing at that time with a 1 or a 0. Things tend to exist in consecutive chunks across the stream. The vertical order of the rows doesn't matter, patterns can occur between any set of rows. I am seeking patterns where dimensions coincide at two or more indices, both at the **start of '1' group's** , and also **where any '1' groups overlap** across all rows/dimensions. Pattern can be offset with nearby adjacent indices, if it's proportion is maintained. D = Dimension/Row N = Index in stream               (( D1(N), D25(N+4), D900(N-1) ), ( 3, 67, 90, 3000 ))          D1(3),    D25(7),    D900(2)     D1(67),   D25(71),   D900(66)     D1(90),   D25(94),   D900(89)     D1(3000), D25(3004), D900(2999)      An example of a pattern match that migh occur at multiple places. Dimension01 at N, Dimension25 at index N + 4, Dimension900 at index N -1 occur at indices 3, 67, 90 and 3000. ### The returned patterns:   * Are ordered by the number of matching indices or by the number of dimensions in the pattern.   * Match at least two dimensions at at least two points How can I go about this?


title: Hypothesis / Statsitical significance between continous variable and dependent binary variable

text: Sorry I am a stats newbie, I understood that I can use the search feature, I tried to search but I am afraid that I am not using the right terms and the results returned are not quite relevant to my question. I have a set of data readings (50 cases), A is a number between [0,1] and B is {yes,no}. My hypothesis / I am trying to show that A does not affect the outcome of B based on my data. May I have some advice what are the right technique / test I should use? Thanks


title: Measurement error in dependent variable

text: Let $Y_{t} \sim N(0,\sigma_{Y_t}^2) $ (independent var), $X_{t}\sim N(0,\sigma_{X_t}^2)$ (dependent var) be time series variables. However, suppose we don't observe $Y_t$ and instead measure it with some error. That is, we have access to $\hat{Y}_t = Y_t - \nu_t$. Assume that the measurement error is of the form $\nu_t \sim N(0,\sigma_\nu^2)$. Consider that I want to use $\{Y_t,X_t\}$ in a Linear Regression Model to be estimated with OLS for $t=1,2,...,T$. Then what I'll actually be estimating is: $(1a) \hat{Y}_t = a + b X_t + e_t$ $e_t \sim N(0,\sigma_{e_t}^2)$ _**Or, equivalently:_** $(1b) Y_t = a + b X_t + (e_t + \nu_t) = a + b X_t + \epsilon_t$ $\epsilon_t \sim N(0,\sigma_{\epsilon_t}^2)$ * * * What I want to know is whether the following is correct: Consider the residual variance co-variance matrix of $(1b)$, $\Sigma := E[\epsilon \epsilon']$, and the analogous population variance covariance matrix $\Omega := E[ee']$ (where "population" is taken to mean model $(1a)$ without any measurement error). **We have that $\textrm{diag}(\Sigma)_t = \textrm{diag}(\Omega)_t + \sigma_\nu^2$ for all $t$** because $\textrm{diag}(\Sigma)_t = E[\epsilon_t^2] = E[(e_t + \nu_t)^2] = E[e_t^2] + E[\nu_t^2] = E[e_t^2]+\sigma_\nu^2=\textrm{diag}(\Omega)_t+\sigma_\nu^2$.


title: Data imputation for statistical analysis

text: **Here is the situation:** I have an individual level data set $X$ where each row is a person $i$ and each column denote characteristics of $i$. The problem is that my data is missing an important variable, lets call this $z_{i}$. To resolve this situation I am considering using **data imputation**. Note, I have not done this before so I am: **(a)** a bit unsure of how it actually works in practice, **(b)** when it can be considered a valid thing to do (and when it is not appropriate). Information pertaining to each of these, or both would therefore be valuable. **My perception of the approach, and proposed way of proceeding is the following (please let me know if this is wrong):** I do not have data on $z$ in my dataset, but I have access to other data from which I could devise a empirical CDF $F(z)$. Next I would simply draw from this distribution and crudely impute these data into a new column in my dataset $X$. Is this the general approach taken? I understand that a more appropriate approach would be to devise the conditional empirical cdf $G(z|X’)$ and to sample form this. I can actually do that because I have a separate sample where a proxy for $z$ is available--call it $z'$, so I could probably use this to get $G(z'|X)$. Would I next simply draw from these conditional distributions in order to construct the new imputed variable in the old dataset $X$? I understand that these questions are rather naive, and I understand that any procedure of this kind must come with its own set of assumptions that will restrict when it is and is not appropriate, but the suggestions above are simply my take on the "intuition" behind the approach, and I firstly want to make sure that this is correct. Second, any information pertaining to the assumptions required on the data would be much appreciated too.


title: SVMs and solution

text: In SVMs, is the solution to the minimization problem $$\textbf{w} = \sum_{i=1}^{n} \alpha_i x_i y_i $$ and once we know $\textbf{w}$ we can get $\textbf{b}$? In plain English can somebody please describe the solution above? Basically the hyperplane that separates the data with the maximum margin is a linear combination of the training data? The quantity $\textbf{b}$ is the distance from the decision boundary to the closest training example? I have not been able to find a simple explanation of the solution above in the books and notes I have looked at.


title: Calculation of distribution of a gaussian process

text: currently finishing the last year of PhD in statistics, we wonder if you could help us with the following. Let $T = [0,1]$ and $X = \left( X_{t}, t \in T \right)$ be a gaussian process with mean function $m$ and covariance function $W$. Parentheses in functions are omitted in this notation. We write $X \sim GP(m,W)$. In a first time, suppose that the sample paths $t \mapsto X_{t}$ are continuous a.s. and let $Y$ be the stochastic process defined on $T$ by : $$Y_{t} = \int_{0}^{t} X(u) du,\mbox{ a.s.}.$$ The question is : what can you say about the law of $Y$ ? In particular, do you believe it to be $GP(m^{\star},W^{\star})$, given by : $$\begin{array}{ccc} m^{\star}(t) & = & \int_{0}^{t} m(u)du,\\ W^{\star}(s,t) & = & \int_{0}^{s} \int_{0}^{t} W(u,v)dudv. \end{array}$$ If the preceding statement is true, could you please give us some references? In a second time, suppose that the sample paths $t \mapsto X_{t}$ are differentiable a.s. and let $Z$ be the stochastic process defined on $T$ by : $$Z_{t} = \frac{dX_{t}}{dt},\mbox{ a.s.}.$$ The question is now : what can you say about the law of $Z$? Do you have any reference where we could find a result? We would be very grateful for any feedback you can give us. Thank you very much.


title: Calculating the Net Promoter Score for poll data on a 1-5 scale

text: From Wikipedia > The Net Promoter Score is obtained by asking customers a single question on > a 0 to 10 rating scale, where 10 is "extremely likely" and 0 is "not at all > likely": "How likely is it that you would recommend our company to a friend > or colleague?" Based on their responses, customers are categorized into one > of three groups: Promoters (9–10 rating), Passives (7–8 rating), and > Detractors (0–6 rating). The percentage of Detractors is then subtracted > from the percentage of Promoters to obtain a Net Promoter score (NPS). NPS > can be as low as -100 (everybody is a detractor) or as high as +100 > (everybody is a promoter). Now we want to scale it down to 1-5 for our purposes, because frankly our users will get confused on the 0-10 scale. `Now we are considering 1,2,3 to be detractors , 4 to be passive and 5 to be a promoter`. Now In the `Standard 0-10 scale, it is actually a 11 point scale – so promoters (9,10) is 2/11th of the distribution which is 18%. In the 1-5 scale we use promoters (5) is 1/5 which is 20%`. So more chances of getting a promoter score. `By same logic, 5 point scale will also show less detractors than 0-10 point scale (3/5=60% vs 7/11 = 63.6%)`. So a positive skew looks possible mathematically. But the fact is a user has only 1 option to choose from on the 1-5 scale to be a promoter, but 2 on the 1-10 scale. So someone who might be a 9 on a 10 scale might be a 4 on the 5 scale and become a passive and never show up and kill our NPS score to go to -100(if no one votes 5 , and a lot of people vote 4) . Three questions > > 1.) Which one of the two math arguments should be considered for scaling > down >> >> 2.) What should be the right passive pivot for the 1-5 scale? > > 3.) Are there other tests that can be done on this data to get an accurate > measure of customer feedback?


title: Significance of linear combination of two insignificant regression parameters

text: I'm running a model with an interaction term, of the general form: $Y=b_1X_1+b_2X_2+b_3X_1X_2$, whereby $X_2$ is a dummy variable and $X_1$ is continuous. Neither $X_1$ nor $X_2$ turn out as significant. However, a $t$-test of the hypothesis that $X_1$ is significant when $X_2=1$ (i.e. $H_0: X_1+X_2=0)$ turns out as significant. Can someone please make intuitive (or econometric) sense of this? Many thanks in advance, Barry


title: Most efficient vector construction for Dynamic Time Warping

text: I'm in the process of folding FastDTW into my SVM and the question now is how to best format my data (irrespective of normalization). Here's an example of what I'm attempting to do - given two 3d vectors: A:               10.3, 5.4, 2.3     7.2, 4.9, 1.4     5.9, 6.2, 3.1     9.2, 4.2, 2.5     2.5, 7.2, 4.2      B:               10.3, 5.4, 2.3     7.2, 4.9, 1.4     5.9, 6.2, 3.1     9.2, 4.2, 2.5     2.5, 7.2, 4.2      I'm not familiar with other DTW implementations but FastDTW allows you to define several columns per vector. If the columns are uneven it complains about it, stating that it cannot compare vectors of different lengths. This kind of confused me as I considered each 'row' to be one data point in a vector - perhaps it considers a vector of vectors? Anyway, as both of them are identical I am getting the expected result:               Warp Distance: 0.0     Warp Path:     [(0,0),(1,1),(2,2),(3,3),(4,4)]      So far so good - now I am modifying three data points horizontally               10.3, 5.4, 2.3     27.2, 14.9, 11.4 ** changed row     15.9, 6.2, 3.1     9.2, 4.2, 2.5     2.5, 7.2, 4.2      I leave vector B as is and this is the result:               Warp Distance: 17.320508075688775     Warp Path:     [(0,0),(1,1),(2,2),(3,3),(4,4)]      As expected the distance has changed (which is the value I'll be extracting further down the line). The only issue I see is that the warp path remains unchanged, but that's probably due to the short vector length. Okay, now we again change three values in vector A but this time vertically:               10.3, 5.4, 2.3     27.2, 4.9, 1.4 *     25.9, 6.2, 3.1 *     19.2, 4.2, 2.5 *     2.5, 7.2, 4.2      I have only changed the value in the first column of three sub-vectors. Similar amount of change but here's the result:               Warp Distance: 30.0     Warp Path:     [(0,0),(1,1),(2,2),(3,3),(4,4)]      The path if obviously the same but it seems that the warp distance has increased quite a bit. Similar amount of change and this suggests that sub- vector changes weigh less than changes across the vector's main axis (between rows). I hope this makes sense to you guys. I'm about to use this for financial time series projection and the plan is to use a window of price, volatility, and other type of market data. The idea it to simply take a window of e.g. 10 vectors and string them all together into one long multi-dimensional vector (just like LIBSVM likes it):               1x     1y     1z     –     2x     2y     2z     –     3x     3y     3z      But HOW to combine the data is the question. I could combine all the volatility data for instance and make it a series of 1-d volatility nodes followed by all the price data points, etc.               1x     2x     3x     1y     2y     3y     1z     2z     3z      Or I could put all the volatility into sub-nodes on one row, followed by all the price nodes, etc.               1x, 2x, 3x     1y, 2y, 3y     1z, 2z, 3z      I'll still get a different weight but the warp path will be much more similar, data normalization/scaling or not. If any of you have any pointers, tips, suggestions, insights - I'm all ears! Thanks in advance :-)


title: Signal dimension in regression model

text: Estimating Unknown Sparsity in Compressed Sensing is a paper about sparse signal. I am just learning the concepts. In the first paragraph, it says that when the number of observation data samples $n$ is less than the signal dimension $p$, then the desired signal $x$ is sparse. In this answer What is signal dimension, the answer is the number of components in a signal.   1. Can somebody please explain what is signal dimension in the context of regression model and what is meant by $p$? How is it different from model order?   2. An example of sparse AR or MA model will be helpful to visualize the difference if any between model order and dimension


title: Ridge and LASSO given a covariance structure?

text: After reading Chapter 3 in the Elements of Statistical Learning (Hastie, Tibshrani & Friedman), I wondered if it was possible to implement the famous shrinkage methods quoted on the title of this question given a covariance structure, i.e., minimize the (perhaps more general) quantity $$(\vec{y}-X\vec{\beta})^TV^{-1}(\vec{y}-X\vec{\beta})+\lambda f(\beta),\ \ \ (1)$$ instead of the usual $$(\vec{y}-X\vec{\beta})(\vec{y}-X\vec{\beta})+\lambda f(\beta).\ \ \ \ \ \ \ \ \ \ \ \ (2)$$ This was mainly motivated by the fact that in my particular application, we have different variances for the $\vec{y}$ (and sometimes even a covariance structure that can be estimated) and I would love to include them in the regression. I did it for ridge regression: at least with my implementation of it in Python/C, I see that there are important differences in the paths that the coefficients trace, which is also notable when comparing the cross-validation curves in both cases. I was now preparing to try to implement the LASSO via Least Angle Regression, but in order to do it I have to prove first that all its nice properties are still valid when minimizing $(1)$ instead of $(2)$. So far, I haven't seen any work that actually does all this, but some time ago I also read a quote that said something like " _those who don't know statistics are doomed to rediscover it_ " (by Brad Efron, perhaps?), so that's why I'm asking here first (given that I'm a relative newcomer to the statistics literature): is this already done somewhere for these models? Is it implemented in R in some way? (including the solution and implementation of the ridge by minimizing $(1)$ instead of $(2)$, which is what's implemented in the lm.ridge code in R)? Thanks in advance for your answers!


title: T-test on ranks vs Mann WHitney U-test equivalence

text: What does it mean when people say that a T-test performed on ranked data is equivalent to a Mann WHitney U -test? Does that mean they just test the same hypothesis/are useful in the same situations or are they are supposed to give the exact same p-values? The reason I ask is I tried both in R and compared two groups with very small sample sizes (3 and 4). I got completely different answers-one significant and one not. The two groups are: A=(1,2,3) B=(4,5,6,7) t-test p value=0.01, M_W test p value=0.06 Thanks!


title: Monte Carlo: generating autocorrelated data from empirical distribution

text: my problem is the following: having a distribution function of daily casfhlows resulting from electricity trading, I need to calculate a yearly 99% VaR, i.e. the 1% percentile of yearly casfhlows distribution. **Method 1/** Use CLT - fast and easy, albeit daily cashflows show some autocorrelation (0.3, 0.2 on the first two lags), hence assumptions are not satisfied. **Method 2/** Simulate as though daily cashflows are independent, i.e. 10 000 times generate 365 random values from the distribution function of daily cashflows and sum them. Like this, I get 10 000 yearly cashflows and get the respective percentile. **Method 3/** Simulate taking account autocorrelation of daily cashflows. This is what I don't know how to do. Having an empirical distribution function X(avg, sdev) how can I successively generate random realizations while ensuring that the generated data have (at least on the first lags) the same autocorrelation structure as the empirical daily cashflows that I have? Is there a way of doing this without estimating an ARMA model? The reason I want to avoid ARMA is that once I randomly generate one observation and model the following observations via the estimated ARMA process, the distribution of predicted cashflows will drift away from the empirical distribution of daily cashflows (due to the expected bad fit of the ARMA process). Any hint is highly appreciated, Daniel


title: Spatio-temporal version of Hammersley and Clifford theorem

text: Does anybody encountered or derived a spatiotemporal version of Hammersley and Clifford theorem? The original theorem (with the proof by Besag) fits just for the temporal models, but if I have a sample of data, that was collected at different times in the same locations, than I need also to incorporate a temporal effect, but still no luck to find any appropriate literature.


title: Reading Kmeans data and chart from R

text: I performed and plotted a kmeans analysis in R with the following commands:                km = kmeans(t(mat2), centers = 4)      plotcluster(t(mat2), km$cluster)      #from library(fpc)      Here is the result from the plot:![enter image description here](http://i.stack.imgur.com/sBWIS.png) What I want to know is how to make sense of the `km$centers` attribute and the plot. What I know is that `km$centers` is a 4 X 31 matrix. Each row represents the corresponding cluster. I think each column represents an iteration in the algorithm (correct me if I am wrong) so the final iteration and result of the algorithm for the centers would be given by:               km$centers[, 31]           0.008785652 -0.088641371 -0.012666252 -0.079348292       I must be wrong about a lot because this leads to the following questions:   1. The centers given by`km$centers` are not (x, y) coordinates. How do I get these (x, y) center coordinates?   2. The center for cluster 4 (according to the plot) must be something like (12, 2) but the above center numbers do not reflect any of these coordinates. In fact every number in the 4 X 31 matrix is less than 1. So, what is the relationship between `km$centers` and the plot? The ultimate goal here is to create a matching (not mentioned here) based upon the (x, y) coordinates. All help is greatly appreciated!


title: Understanding the p-value in Spearman's rank correlation

text: I have some pairs of datasets (`n=200` or thereabouts), of samples which are non-negative and not normally distributed. I think these pairs of variables are related, probably linearly. Calculating Spearman's rank correlation on these datasets gives some strange results. The correlation coefficients show that the pairs of variables are weakly, positively correlated (e.g. `rho` of around `0.4`), but the p-values are very low (e.g. `4.1e-10`). My vague understanding of this is that the variables are weakly, positively correlated but the probability of unrelated variables producing the same correlation is very low. Does this mean we can be reasonably certain that a positive correlation exists or have I misunderstood?


title: How to interpet residual plot?

text: As a part of a design of experiments course I'm taking, I ran an experiment at home. The experiment was checking how water boiling time changes under certain factors (5 overall factors) all which had 2 levels: with or without salt, with or without oil, type of pot, type of water and the diameter of the pot. The experiment was a fractional factorial of [2^(5-1)]*2. I ran the analysis an got this residual plot: ![residual by predicted](http://i.stack.imgur.com/cVGah.jpg) ![Residual by row](http://i.stack.imgur.com/JHwY2.jpg) I know it's supposed to be random, but it doesn't seem like it.


title: How do I convert unstandardized RMR estimates into standardized RMR estimates?

text: I am having difficulty turning unstandardized RMR estimates from AMOS software output in a SEM into standardized estimates to enable suitable interpretation. Can someone give me advice on where to start with this?


title: How many participants do I need to carry out an four-way Anova?

text: I am looking at attachment and cyber addiction. I'm carrying out an four-way Anova, attachment patterns being my ivs (secure, preoccupied, fearful and dismissing) with cyber addiction as my dependent variable. I am intending to use 100 participants in my study, I wanted to know if that is an appropraite amount to run a four-way anova?


title: How to convert molecular categorical variables to dummy variables for cluster analysis?

text: I would like to use a clustering method, e.g. 'mclust', in R to classify each individual in my dataset to k groups. I have 7 continuous and 3 categorical variables. These and other hierarchical clustering methods do not allow for use of categorical variables. Searching Google and this site it appears that converting the categorical variables to dummy variables is an option. Two of my categorical variables are SNP's (molecular markers). They are coded, 0 (fixed allele), 1 (heterozygous alleles), 2 (fixed allele). Do you have any suggestions on how to convert this to a dummy variable or am I incorrect in my thinking that the dummy approach is appropriate?