title: cross validation vs shuffling

text: I have a classifier. when I shuffle input data (a few times), make training/test data and run classifier, the accuracy is high but when I use cross validation in sklearn it drops by 15% (both recall and precision). Wondering what is the problem? and how to fix it?


title: Interpreting distance between classes for binomial, multinomial, and ordinal regression

text: I would like to understand the difference between the binomial, multinomial and ordinal models when it comes to interpreting the "distance" between classes, in terms of response probabilities. To illustrate my question, I have prepared the following graph (just for the sake of example) ![enter image description here](http://i.stack.imgur.com/urPiq.png) Imagine I have use the logit link, and plotted the cumulative curves corresponding to 3 classes with a series of binomial regressions (first case), then I applied multinomial and ordinal regressions (cases 2 and 3, say for 4 classes, one of which is a reference class - not represented because the cumulative probability is 1). In every case, I will be able to estimate a mean response (probabilities) for each class, $\mu_i$. But the parameters will be different - understandably, as the formulation of the odds in the link function is different. Now what I dont understand, is how to interpret the distance between the classes (in green in my graph), in other words all the probabilities that are smaller than $\mu_2$ but greater than $\mu_1$, for example. What are the differences in interpretation according to the model I am choosing?


title: NN forecasting data source samples

text: I'm interested about the field of energy consumption forecasting, and because I'm in the learning stage I'm looking for some source codes and samples. More exactly I'm interested about forecasting with neural networks approach. Yesterday I've found on CV the link to the Neural Forecasting Competition but unfortunately I didn't find any results nor source codes. Also, I'm interested about some methods that were used to resolve the problems with the given data sets. Can you help me with some links or advices about where can I find some data sets&solutions in this energy forecasting field?


title: Count average on data for a long time period

text: I have a program system, where present a lot of tasks. Each task need some time for an execution and monitoring system writes time of execution for each task to a database. I need to calculate average execution time for each task for a long time period (much more than half year while it executes at least four times per month, in the worst case it executes once per day). Unfortunately I can't just calculate arithmetic average because task can be modified and task duration will increase or decrease. Also there may be some errors in task and it will fail and "spoil" time statistics. I thought that may be it is usefull to calculate average only on several last runs, but seems it won't give me information relevant enough. So my question is how to calculate average task running time with these clarifications? ps. not sure it should be here.


title: Standardizing single coefficient in multivariate analysis

text: I have panel data and for the I have a following equation $$ logY = \beta_1 + \beta_2 logX + \beta_3 m logW $$ Problem is with $\beta_3$ coefficient. Since m is outside the log and it is a very small value, the overall $mlogW$ is very small and the resulting coefficient is very large (in thousands). One alternative would be to use standardization but my results change drastically if I standardize the variables. A solution I thought of was to multiply the variable with a constant such as 10000. So the new equation will become: $$ logY = \beta_1 + \beta_2 logX + \beta_3 m logW*10000 $$ My question is, is this a right thing to do? It should only change the magnitude of the coefficients just for the sake of reporting the coefficients but actually it is changing how all the other coefficients behave when I multiply it say by 1000 or 10000. Or is there a way to standardize the coefficients after estimation? since I have panel data, the division and multiplication with standard deviations seem complex. Any idea how to make it work? Thank you


title: cubic relationship after linear relationship

text: Sorry for my title, but I really don't know how to describe this question. I am fitting a linear regression in R now, and I find that there is one parameter showing linear relationship before certain point and cubic relationship after that point. What I did is to separate the whole dataset with reference to that point and fit the model with parameter value and cubic value respectively. However, that point is determined manually. Does anyone here know how to determine that point with formal algorithm? Thanks in advance


title: How to interpret merits in Weka with ChiSquaredAttributeEval and SVMAttributeEval?

text: I want to interpret the goodness of attributes using feature selection with 10-fold cross validation. With ChiSquared I get something like this (deletet attributes with merrit was 0 in all folds):               average merit      average rank  attribute     13.632 +- 1.6        1   +- 0        8 E8      or               average merit      average rank  attribute      3.957 +- 0.306      2.1 +- 0.54     5 E5      8.015 +- 2.672      4.2 +- 8.94    23 E23      1.331 +- 3.994      7   +- 2.1     14 E14      1.095 +- 3.285      9   +- 2.68     8 E8      With SVM i get:               average merit      average rank  attribute     34     +- 0          1   +- 0        1 class (perfect attribute to compare merit)     31.5   +- 2.335      3.5 +- 2.33    13 E12     30.8   +- 1.887      4.2 +- 1.89     4 E3     30.6   +- 1.908      4.4 +- 1.91     9 E8     28.7   +- 3.407      6.3 +- 3.41    10 E9     27.6   +- 2.01       7.4 +- 2.01    20 E19      I couldn't find descriptions about what the numbers are telling me. Is it just a relative value to say which attribute is better than others or is there a meaning behind?


title: Pie charts vs. dot plots

text: I understand the critiques of pie charts as referenced here: Problems with pie charts However, the above response (and the R manual) always cite dot plots from Cleveland as an alternative. My question is why are dot plots considered an alternative? It seems to be that dot plots only apply when:   1. the data set is relatively small (so that you can perhaps even by eye count dots in the plot, each dot corresponding to a data point, and    2. they are not intended to display percentages / density estimates of data.  It seems to me like the main goal of the pie chart is (even if it does so badly) to highlight the fact that the categories must sum to 1 or 100%. with the dot plots you can show the cardinal value of each category, but it will not be obvious what fraction of all the categories a particular value is. Also, with a pie chart it doesn't matter if you have 10 data points or 1,000,000, where as a dot plot with a million points (even if the number of categories is small) seems odd and might just collapse to a bar graph. Could someone explain why dots plots are seen as an alternative and maybe provide a few examples of dot plots to other quantitative alternatives to pie charts (like bar plots?)


title: P-value correction in R: what is the parameter "n" used for and when should I use it?

text: The function `p.adjust()` in R accepts a parameter `n` that defaults to the length of the vector of p-values one wants to adjust: `p.adjust(p, method = p.adjust.methods, n = length(p))` In the documentation the explanation for `n` says: > n number of comparisons, must be at least length(p); only set this (to non- > default) when you know what you are doing! What does this mean? When could it be helpful to set it to a value different from the default?


title: Graphical representation of cross-validation errors for regression

text: What are some good ways of presenting/comparing cross-validated RMSE errors for regression using various models, graphically via plots? As of now, I have been presenting the quantitative results in tabular form.


title: How to find a sequence of random variables of infinite expectation that converges to zero?

text: I want to find an example such that $X_n\rightarrow0$ as $n\rightarrow \infty$ while $E(X_n)=\infty$. I am thinking of making use the divergence of $\sum\frac{1}{n}$, but fail to find suitable $X_n$ and $\Pr(X_n=k)$. Can you show me some examples or give me some hints?


title: How calculate average probabilities in MLP or SVM?

text: I have a system that find best model (best inputs and parameters of MLP/SVM) model in a financial problem for every inserted database and create a specific model for a specific data sample. I'm using 5-fold cross validation that repeats for 5 times (for increasing reliability of my system) for cost function of my optimization algorithm so In the end I will find best 25 models that now will use for out-sample data. When I insert out sample data to these 25 models in some cases I have slightly different outputs ( I have a classification MLP model that can generate probability). This is one of the outputs :               0.380032703, 0.356491904, 0.369914169, 0.377588996, 0.425321986 ,0.382467546, 0.341404809, 0.538460421, 0.353585952, 0.332340499, 0.36827261, 0.369171894, 0.329522399, 0.398054125, 0.360501952, 0.406897298, 0.480889919, 0.403619213, 0.406236551, 0.422529236, 0.393108809, 0.309808029, 0.357776713, 0.337373411, 0.397259506      As you know these outputs have difference parts of data as training sets and testing sets for creating model. Other things like inputs and parameters are the same.I should have only one probability for any out sample. My fist approach was using average. I want your suggestions to create one robust probability from these probabilities. When I check other out sample outputs, In a lot of cases difference is low but in some cases difference between probabilities is high. What is your idea about this problem? Thanks.


title: Can the coefficient of variation for two samples be added?

text: I have a quick question that should hopefully require a quick answer (my apologies if this is a rather simple question). I have two estimates of relative abundance from two sub-populations and want to sum them to obtain one overall number. Being that they're just estimates, I also have their corresponding measures of uncertainty (i.e. coefficients of variation, CVs). I know that if the output is a linear combination of the inputs with coefficients equal to 1, I can come up with an estimate of the variance for the output by adding together the variances of the inputs. For example: \begin{align} y &= x_{1} + x_{2} \\\ var(y) &= var(x_{1}) + var(x_{2}) \end{align} However, I'm obtaining the estimates from published reports and they express uncertainty in CVs and not variances. Can I add the CVs together like variances? Thanks in advance for your help!!


title: 2x3 ANOVA interaction no longer significant after including a covariate?

text: My experiment had a 2x3 design with one covariate. If I analyse the results using an ANOVA, I get a strongly significant interaction between the two main factors (p<.001). If I add the covariate into the analysis and do an ANCOVA, the interaction between the two main factors is not significant (p=.252). The covariate itself was not significant, nor did it have any significant interactions. There are similar questions on this website, but I was hoping to understand what does such a result pattern say about the covariate? How should I interpret the results? Adding the covariate was the novel feature of my study and therefore it is important to me to understand this interaction properly. Is the covariate influencing the results? EDIT: I've been thinking about this all day and I guess the possible interpretation would be: When controlling for the covariate, the interaction of the main effects is no longer significant. In which case I don't understand a couple of things: a) Does that mean the covariate is responsible for the interaction? b) Why is the interaction between FactorA*FactorB*Covariate not significant then? c) Why isn't the covariate significantly interacting with anything at all, when its presence seems to be influencing the results? EDIT 2: Maybe I am doing something wrong, because it is not making much sense to me. Here are my results: http://imageshack.us/g/580/ancova1.png/ The last two pictures show when I run an ANOVA without the covariate (repeated measures, Factor 1 (2 levels). Factor 2 (3 levels). The first three when I run an ANCOVA on the same factors, except I put in OptimismScore as a covariate.


title: Sampling using PCA and maximizing variance

text: In 6 forests, I've 3 variables that I know the values for all trees: $height$, $diameter$ and $density$. I can only cut some trees to know their $mass$. The question is how to choose the trees to cut? As I supposed a linear model such as: $Mass = a_1 Diameter + a_2 Height + a_3 Density + a_4$ (or as mention in an answer $log(Mass) = a_1 log(Diameter) + a_2 log(Height) + a_3 log(Density) + a_4$), I thought to maximize the variance of the sampling to have the best fitting afterward. Is it correct to use a Principal Component Analysis over the 3 variables and then choose in the 2 first axes of the PCA the most spread out sampling ? Is there another method to choose sampling by maximizing the variance over variables that I already know ? Thank you.


title: Estimating parameters of inifinite scale mixture from data

text: Suppose that I have an infinite scale mixture of zero-mean normal distributions, whose mixing distribution is gamma with parameters $\alpha$ and $\beta$. The data is thus distributed according to a generalized Student's $t$ distribution. I am handed a pair of samples, $(z,y)$, generated through the following procedure, and asked to estimate $\alpha$ and $\beta$. First, $x$ is sampled directly from the Gamma distribution. Then it is corrupted by gaussian noise, so that $z = x+N(0,\delta)$. Finally, $y \sim N(0,1/x).$ Are there known estimators for $\alpha$ and $\beta$?


title: How to validate "equal slopes" (proportional odds) in ordinal regression

text: I would like to fit an ordinal regression model using proportional odds. I learned to test for "equal slopes" in order to say something about the model's validity. Therefore, I fit a model _with_ equal slopes and I fit a model _without_ equal slopes (e.g. with R using the VGAM package). After that, I do a Chi-square test to compare both results. If this test _is not_ significant, we have equal slopes. If this test _is_ significant, we do not. A colleague states that those tests lack power. With such low power, it could happen that the test is _not_ significant (e.g. $p=.2$) while there _are_ different slopes. My questions are:   1. Is there a more powerful way to test for equal slopes?   2. Are there other ways to test the validity of the model? Thx in advance!


title: How to estimate Relative Risks in Multivariate Binary Logistic Regression Models, instead of Odds Ratios?

text: All software programs I have tried report only odds ratios (ORs) for binary logistic regression predictors (as exponential of the betas). I am interested to know how can I compute the relative risk (RR) from a binary logistic regression model? My reason is that RRs (besides ORs) would make my variables more understandable. Also I see many journals (even many of the top ones) confuse these two, and I want to hopefully have a good report with scientific merit, and not repeating previous mistakes. I found some methods over the net (which I have not tried yet). Here and Here for example. However in a forum, I saw some experts had opposed reporting RRs for logistic regressions. For example: "It is incorrect to use the relative risk as a measure of association in a logistic regression. The measure of association in a logistic regression is the odds ratio. The odds ratio is an approximation of the relative risk. The approximation becomes progressively better as the disease becomes progressively rarer. Regardless of whether the disease is rare or not, inferences drawn from a logistic regression are valid. Please do not report a logistic regression using relative risk. It is not correct to do so." --John Sorkin Or "I am curious why one would want risk ratios. Unlike odds ratios, they are not interpretable without reference to the base risk. For example a risk ratio of 2 cannot possibly apply to anyone with a starting risk exceeding 1/2." --Frank Harrell So my questions are:   1. Is it a good way to estimate RR for binary logistic regression? Or do you agree with the above quotes, disagreeing with RR for logistic regressions?   2. Could you please let me know why, if you agree or disagree?   3. Do you know of any implemented algorithms (e.g., a macro or an R function) which can do it without needing to manually computing it? Many thanks.


title: Real world challenge: how to predict without training data

text: For one campaign, I have X as independent variable, Y as dependent variable, i.e., X is income, Y is credit score. I can use X to predict Y with proper model. Now I want to know how Y would perform if I added one more independent variable Z into the model, i.e., age. The challenge is I don't have any data on Z, that is to say, I don't know the relationship between Y and Z for this campaign. But I do have data that illustrates the relationship between Y and Z for other campaigns. Is there any model that can be used to answer this question? Any thought is helpful.


title: Martingale difference vs. sub-independence

text: We know that martingale difference is an assumption weaker than dependence, and stronger than uncorrelated. The same is true for sub-independence. My question is how sub-independence and martingale difference are related.


title: Can two random variables have the same distribution, yet be almost surely different?

text: Is it possible that two random variables have the same distribution and yet they are almost surely different?


title: How to understand standardized residual in regression analysis?

text: According to _Regression Analysis by Example_ , the residual is the difference between response and predicted value, then it is said that every residual has different variance, so we need to consider standardized residuals. But the variance is for a group of values, how could a single value have variance?


title: Estimate the population variance from a set of means

text: I have a set of measurements which is partitioned into M partitions. However, I only have the partition sizes $N_i$ and the means $\bar{x}_i$ from each partition. Because all measurements are assumed to be from the same distribution, I believe I can estimate the mean of the population, $\bar{y}$, and standard deviation of the mean, $\sigma_{mean}$: $$ N=\sum_{i=1}^M N_i $$ $$ \bar{y} = \frac{1}{N}\sum_{i=1}^MN_i\bar{x}_i $$ $$ \sigma_{mean}=\sqrt{\frac{1}{N}\sum_i N_i(\bar{x}_i-\bar{y})^2} $$ My questions:   1. Am I right in my assumptions, that the mean $\bar{y}$ can be computed as above?   2. How can I find the standard deviation for the population, given only the means? I read that the standard deviation of the population and standard deviation of the mean is related with $$ \sigma_{mean}=\frac{\sigma}{\sqrt{n}} \mbox{[1]} $$ where $n$ is the number of samples used in the computation of $\bar{x}_i$. So is it actually as simple as just multiplying $\sigma_{mean}$ with $\sqrt{n}$ if $n$ for all means are the same?    3. If it's that simple, what do I do if each $\bar{x}_i$ is computed using a different number of samples?  [1] Wikipedia:Standard Deviation


title: Notation of estimators (tilde vs. hat)

text: 1\. Is there any naming convention regarding the hat and the tilde symbol in stats? I found $\hat{\beta}$ is describing an estimator for $\beta$ ( Wikipedia ) But I also found $\tilde{\beta}$ is describing an estimator for $\beta$ (Wolfram ). Is there any difference in the meaning?   On the web I found kind of a difference but I am not sure about the meaning Reference for Stats Symbols. There it is distinguished between "estimates of parameters" and "estimates of variables". Could someone be so kind to explain in which case to use tilde and hat? 2\. Regarding the expectation operator, is there any difference in $E(X)$ and $E[X]$ and $E\\{X\\}$ regarding the brackets? I got the advice to use the curly brackets. But I am not sure about the meaning. I used to use the brackets only for reading/visualization rather than pointing on some meaning. Any advice on that?


title: Command-line tool to calculate basic statistics for stream of values

text: Is there any command-line tool that accepts the flow of numbers (in ascii format) from standard input and gives the basic descriptive statistics for this flow, such as min, max, average, median, RMS, quantiles etc? The output is welcome to be parseable by the next command in command-line chain. Working environment is Linux, but other options are welcome.


title: How to normalize Poisson distributed data before PCA

text: I am a newbie to principal component analysis (PCA). I will have to do PCA to data sets consisting of count statistics: all data are positive integers. Before PCA the data needs to be normalized. It is more or less standard to do that by subtracting the mean and dividing by the standard deviation in the variable over the sample set. I wonder whether this is appropriate for data sets that are very skewed.


title: Combining experimental error in a mean

text: I understand the rules for combining experimental errors in sums, differences and ratios (as explained here), but what happens to an experimental error when you average it? Say, ruler measurements of the length of beetles that are all values like 12.3 cm $\pm$ 2mm.   I can find plenty of explanation on the web about how to add and multiply values like that together, but what happens to the error when you take an average of (say) 10 measurements?


title: Do you know if I can find database with volleyball statistic?

text: I'm student and I have project from Data Warehouses's course. I have to create cube from database (or I should create database from raw date). It is required that the data are real (I don't mean real number). Maybe someone with your experience knows source of date on statistics in volleyball


title: Implication for OLS of not using time component with panel data?

text: I have a subset of data from a panel covering 100 cities and 20 years but the subset for some reason does not include the time variable or a variable to distinguish one city from another. What are the implications for OLS? Would this violate some assumptions? How can I get around it (if at all)? Can I treat the data as cross-sectional data? Even though there are 20 observations per city across the variables. **Additional Info:** To be more specific my dependent variable is year-on-year change in population and the main regressors are year-on-year changes in other variables. Does the fact I am using changes get around these problems at all? (I am thinking it would as long as I can assume the city specific factors are constant).


How can reliability and validity of content analysis be quantified when there is only one person coding the data?


title: Confidence bands for difference of time series

text: Assume that I have two time series $Y_{1t}$ and $Y_{2t}$ that are sampled at the same frequency. Is there a way to quantify the uncertainty in their difference $Y_{1t} - Y_{2t}$? That is, can we get confidence bands on $Y_{1t} - Y_{2t}$? My thought is that some sort of a dependent bootstrap should apply. This question also seems related to this question.


title: How to calculate the probability $P(A' \cup B')$ and related properties?

text: Let $A$ and $B$ denote two events, and suppose that $\text{P}(A) = 0.2$, $\text{P}(B) = 0.3$, and $\text{P}(A\cap B) = 0.1$. Are the following computations correct?   1. $\rm{P}(A \cup B) = \rm{P}(A) + \rm{P}(B) - \rm{P}(A\cap B) = 0.4$   2. $\rm{P}(B') = 1 - \rm{P}(B) = 0.7$, where $B'$ denotes the complement of $B$   3. $\rm{P}(B\cap A') = \rm{P}(B)\times \rm{P}(A') = 0.24$   4. $\rm{P}(A' \cup B') = \rm{P}(A') + \rm{P}(B') - \rm{P}(A'\cap B') = 0.94$ [NB: $\cup$ = "or" and $\cap$ = "and"]


title: Expected number of times the empirical mean will exceed a value

text: Given a sequence of i.i.d. random variables, say, $X_i \in [0,1]$ for $i = 1,2,...,n$, I'm trying to bound the expected number of times the empirical mean $\frac{1}{n}\sum_{i=1}^n X_i$ will exceed a value, $c \geq 0$, as we continue to draw samples, that is: $$ \mathcal{T} \overset{def}{=} \sum_{j=1}^n \mathbb{P} \left(\left\\{ \frac{1}{j}\sum_{i=1}^j X_i \geq c\right\\}\right) $$ If we assume that $c = a + \mathbb{E}[X]$ for some $a > 0$, we can use Hoeffding's inequality to arrive at \begin{align} \mathcal{T} & \leq \sum_{j=1}^n e^{-2ja^2} \\\ & = \frac{1 - e^{-2 a^2 n}}{e^{2 a^2}-1} \end{align} Which looks nice (maybe) but is actually quite a loose bound, are there any better ways of bounding this value? I expect there may be a way since the different events (for each $j$) are clearly not independent, I'm not aware of any way to exploit this dependence. Also, it would be nice to remove the restriction that $c$ is greater than the mean. **edit** : The restriction on $c$ being greater than the mean can be removed if we use Markov's Inequality as follows: \begin{align} \mathcal{T} & \leq \sum_{j=1}^n \frac{\frac{1}{j}\mathbb{E}[X]}{c} \\\ & = \frac{\mathbb{E}[X]H_n}{c} \end{align} Which is more general, but much worse than the above bound, although it's clear that $\mathcal{T}$ must diverge whenever $c \leq \mathbb{E}[X]$.


title: How to simulate a Cox proportional hazards model with change point and code it in R

text: I have a model that has the following characteristics:   * The covariate $X$ follows a $Be(1/3)$.   * If $X=0$, survival time $Y$ follows an $E=Exponential (1)$.   * If $X=1$, survival time $Y$ is generated as $E$ if $E\le\Psi$, and as $\Psi+E_\lambda$ if $E>\Psi$, where $\Psi$ is the 'change point' (let $\Psi=1$) and $E_\lambda=Exponential(\lambda)$ independent from $E$. Does anyone know how to simulate this model and translate it in R code?


title: How to project a vector onto matrix of rotation in PCA

text: I want to project a vector onto space transformed by PCA. I've calculated PCA in R language using `prcomp`. Now I should be able to multiply my vector by matrix of rotation. Should principal components in this matrix be arranged in rows or columns? Vector that I want to project:               attr1 attr2 ...       Rotation matrix:               PC1 PC2 PC3 ...      ... ... ...


title: Arithmetic relation between expectation values of transformed random variables

text: Let $X_1$ and $X_2$ be discrete non-negative random variables with $E(X_1) \le E(X_2)$. Assume a function $f$ with the following properties:   * $f$ is positive, i.e. $f(x) \ge 0$   * $f$ is (strictly) monotone increasing   * $f$ is concave  I conjecture that $E(f(X_1)) \le E(f(X_2))$ but couldn't find an appropriate reference. Could you provide a reference for a proof of this (or a more general property implying this case) and/or sketch the basic proof idea? Remark: The concrete case for which I require this property is for random variables $X_1,X_2$ that evaluate to $[0,1]$ and $f(x) = 1-(1-x)^n$ for $n\in\mathbb{N}, n \ge 1$.


title: Probability to find the treasure five times in row

text: Let's pretend there are 10 doors in front of me. Behind one of them, there is a treasure and nothing behind all the others. Intuitively, I can easily determine I have 1 chance out of 10 the pick open the right door (the one with the treasure). I repeat this 5 times (i.e. all door closed, open one door). What are my chance to find the treasure the five times? In brief, I am wondering how to determine the probably that I find the threasure 5 times in row. Thansk


title: What is the difference between test of independence and test of randomness in linear regression?

text: In section 3.3 of Kutner's Applied Statistical Linear Models: > **Nonindependence of Error Terms** > > Whenever data are obtained in a time sequence or some other type of > sequence, such as for adjacent geographic areas, it is a good idea to > prepare a sequence plot of the residuals. > > ... When the error terms are independent, we expect the residuals in a > sequence plot to fluctuate in a more or less random pattern around the base > line 0. ... In section 3.4 > **Tests for Randomness** > > A runs test is frequently used to test for lack of randomness in the > residuals arranged in time order. Another test, specifically designed for > lack of randomness in least squares residuals, is the Durbin-Watson test. > This test is discussed in Chapter 12. What is the difference of testing independence of error terms and testing of randomness in the residuals? I guess they are the same thing? Thanks!


title: Investigating evidence of autocorrelation in a regression model

text: Is there anyone who can explain as simple as possible how to go about doing this? I have been watching different youtube videos and using google as well, however, I have not understood how to go about approaching this. I know of the Durbin-Watson Test, of dl, du, 4, 0, 2 but it just doesn't make sense to me. I have problems analyzing and applying this procedure for investigating whether there is evidence of autocorrelation. I would appreciate your help.


title: Over represented gender in linear modelling data

text: I am fitting a model to some sales data and looking to accurately represent male/female behaviour (e.g basket analysis). For example, I know the proportion of females/males in the population buying these products are 50:50, and my products are identical to competitor products, however, my data contains 30% females and 70% males. I believe this may bias my model, and from what I can learn, I need to apply a weighting function. I am using R and lm() to model. If using weighting is the right way to go, how is it represented in lm(x,y,weighting=)? E.g. do males receive a weighting of .3 to counter their over representation? or what would be the appropriate weighting in this case?


title: How to mathematically prove that we are sampling from same distributions?

text: The content of this question is about rigorously proving something which is otherwise considered easily correct intuitively. Let's assume we have a multivariate distribution $g(x_1,x_2,...,x_n)$ over the variables $x_{1:n}$. Let's assume that we know how to sample from that distribution, too. We draw samples $x^{1}_{1:n},x^{2}_{1:n}, ... ,x^{N}_{1:n}$ from this distribution. Then we assume that we have the distributions $f_{1}(x_1), f_{2}(x_2|x_1),f_{3}(x_3|x_2,x_1),...,f_n(x_n|x_{1:n-1})$ which is equal to $g(x_1,x_2,...,x_n) = f_{1}(x_1)f_{2}(x_2|x_1)f_{3}(x_3|x_2,x_1),...,f_n(x_n|x_{1:n-1})$ by the chain rule of probabilities. Now, for each sample $x^{i}_{1:n}$ we first sample $x^{i}_{1}$ from $f_{1}(x_1)$, then $x^{i}_{2}$ from $f_{2}(x_2|x_1)$ up to $x^{i}_{N}$ from $f_n(x_n|x_{1:n-1})$. We again obtain $N$ samples. Intuitively we know that the first $N$ samples which come from $g(x_1,x_2,...,x_n)$ and second $N$ samples each of which sequentially come from $f_{1}(x_1),f_{2}(x_2|x_1),f_{3}(x_3|x_2,x_1),...,f_n(x_n|x_{1:n-1})$ are identically distributed. But how can we show this fact in a mathematically rigorous way? I could not think of any procedure and became stuck. Thanks in advance


title: What test of significance should I use to analyse visits to a clinic by various groups?

text: Apologies for the noob question. After a day on Google and Wikipedia I still can't quite work out what to do. So here I am. A small private healthcare clinic in the UK has asked me to look at their patient appointments data. Patients hear about the clinic from various sources (Google, word-of-mouth, leaflets, etc etc). Plotting [patient numbers] against [number of appointments per patient] appears to reveal intriguing patterns in patient "loyalty": ![plot of patient numbers and average appointments per patient](http://i.stack.imgur.com/bSxBA.png) The X-axis is each source:   * saw the clinic while visiting another business in the same building   * recommended by another patient   * friend of the clinician   * etc   * Google   * saw our street sign   * etc   * etc etc etc The pink bars show total number of patients per source. Not too surprising that Google is up high, as is the street sign. But intriguingly the patients from Google/street-sign only come for 3 or 4 appointments - see the blue bars - whereas the patients to the left from sources including word-of-mouth appear to be more "loyal", visiting for a higher mean number of appointments. I'm curious if this is a real effect, or a glitch in the smaller sample sizes of the left-most data. The left-most 4 sources have sample sizes of 7, 12, 6 and 9. Exactly which test should I use to determine significance here? I have month- by-month totals for each source for the past year, like this: ![month by month total patient numbers per source](http://i.stack.imgur.com/qoPAL.png) As you can see, I've started calculating standard deviations and standard errors. However I'm a little stuck now. Should I be using ANOVA here? If so, how exactly? I should add I'm a novice at this so feel free to use small words and type slowly. I've looked through "similar questions" but I'm still no wiser. Many thanks for your help.


title: Estimating Failure Rate from Observed Data

text: I recently read the excellent book Probabilistic Programming & Bayesian Methods for Hackers and I'm trying to solve some problems on my own: > I perform an experiment to estimate the reliability of a device (e.g. CPU, > car, fan, etc.) by purchasing 1000 devices and running them over a 1 yr > period to see how many fail. My experiment shows that 0 of them have failed > over the period. What is my 95% confidence interval of the failure rate? If the number of devices that failed were non-zero, I could estimate the mean and confidence interval with standard techniques like calculating the standard error or bootstrapping. The fact that 0 devices failed makes it much trickier. I try and use PyMC to solve it:               import numpy as np     import pymc as pm          data = np.zeros(1000)  # observed data: zero failures out of 1000 devices     p = pm.Uniform('p', 0, 1) # model the failure rate as a uniform distribution from 0 to 1     obs = pm.Bernoulli('obs', p, value=data, observed=True) # each device can fail or not fail. i.e. the observations follow a Bernoulli distribution     model = pm.Model([obs, p])          mcmc = pm.MCMC(model)     mcmc.sample(40000, 10000, 1)          print np.percentile(mcmc.trace('p')[:], [2.5, 97.5])          > [3.3206054853225512e-05, 0.0037895137242935613]        1. Is my model correct?   2. Is my usage of PyMC correct?   3. Can someone confirm my answer with another method? Maybe analytically by using a Poisson distribution? P.S. I practically know nothing about the theory behind MCMC but the application-first and mathematics-second approach of the book is excellent.


title: Optimal feature selection for MAPE criteria with RandomForest cross-validation

text: I am trying to optimize my set of features against random forest cross- validation using MAPE criteria. I tried forward selection with Univariate linear regression test (f_regression in sklearn), I calculate MAPE for each set of variables selected by SelectKBest:               for i in range(1,len(X.columns)):         selektor = SelectKBest(f_regression, k = i)         clf = RandomForestRegressor(n_estimators=10, max_depth=None)         pred = clf.fit(X, y).predict(T)         MAPE = mean(abs(pred-y_real)/y_real)      I also tries backward selection with the RandomForest feature_importance attribute. I start with full set and in each iteration I remove the least important feature, and I calculate MAPE. I remove then, the feature with the least important MAPE:               while X:          clf = RandomForestRegressor(n_estimators=n, max_depth=None)         pred = clf.fit(X, y).predict(T)         imp = dict(zip(list(X.columns.values), clf.feature_importances_))         todrop = min(imp.iteritems(), key=itemgetter(1))[0]      First technique returns decent results but it is not optimized because I select K best features with liner regression score criteria. Second technique returns noisy results (MAPE doesn't converge). I am looking for a technique for feature selection to use in sklearn that will measure MAPE directly. The size of my full set is 150 features and 80,000 observations. Thank you in advance for any suggestion...


title: Autocorrelation in the presence of non-stationarity?

text: Does the autocorrelation function have any meaning with a non-stationary time series? The time series is generally assumed to be stationary before autocorrelation is used for Box and Jenkins modeling purposes.


title: Best way to handle unbalanced multiclass dataset with SVM

text: I'm trying to build a prediction model with SVMs on fairly unbalanced data. My labels/output have three classes, positive, neutral and negative. I would say the positive example makes about 10 - 20% of my data, neutral about 50 - 60%, and negative about 30 - 40%. I'm trying to balance out the classes as the cost associated with incorrect predictions among the classes are not the same. One method was resampling the training data and producing an equally balanced dataset, which was larger than the original. Interestingly, when I do that, I tend to get better predictions for the other class (e.g. when I balanced the data, i increased the number of examples for the positive class, but in out of sample predictions, the negative class did better). Anyone can explain generally why this occurs? If I increase the number of example for the negative class, would I get something similar for the positive class in out of sample predictions (e.g., better predictions)? Also very much open to other thoughts on how I can address the unbalanced data either through imposing different costs on misclassification or using the class weights in LibSVM (not sure how to select/tune those properly though).


title: Multiple dependent variables in factorial design

text: I have a 2x2x2 factorial design with two dependent variables (lets say height and weight). I can examine the effect of the three factors for each dependent variable separately. But I also want to check if height has anything to do with weight? I.e., is the tallest subject the heaviest?   * Any idea on how to analyse this in spss??   * Also what would the assumptions be?


title: Alternatives to bag-of-words based classifiers for text classification?

text: Most of the text classifiers are based on the bag-of-words approach where you loose the context that a particular word appears. As a solution (or simple solution?) we can use n-grams as features. But are there any classifiers which "gist" the idea and model it in someway before training?


title: R time series analysis in Hadoop

text: need to do time series analysis(TSA) and forecasting on over 100 customers in Hadoop environment(RHadoop), for each customer, it has its own time series which I can analyze individually, but when I do it in Hadoop, is it possible that each customer's time series might be broken into pieces(i.e it doesn't have continuous time series from begin date to end date) - if possible, is there a way to prevent this from happening? any example code of time series using R in Hadoop is highly appreciated


title: How to detect outliers with longitudinal data?

text: I am running a pooled OLS and Random Effects (RE) model and I would like to test for whether there are any outliers. I know how to do this for OLS, but I just dont know how to do it for Random Effects model? If there is outliers when running pooled OLS does it mean that outliers are also present when running a Random Effects model? And if there is outliers, what to do in RE?


title: Natural language processing application in algorithmic trading and sports betting

text: I'm want to do a PhD in NLP and I'm defining the topic at the moment. I've heard that NLP can be used to track online sentiment which in turn can be used for algorithmic trading on stock markets. See this hedge fund I was wondering if NLP techniques, maybe with online sentiment tracking, could be used to predict sport events on online betting sites? Would that be feasible? If yes, what research has been done on the subject?


title: Can a coin be biased?

text: If I have a coin that is not necessarily expected to turn up heads half of the time, then is it correct to call the coin biased? Or does _bias_ in statistics only mean the bias of an estimator in which case it would be inappropriate to describe the coin itself as biased? This question does not mean to be about issues like this, so maybe a different example would have been better.


title: Chi-square contrast tests with R

text: I have Group A, B, C. My hypotheses are" 1) proportions of C are different from proportions of A (=independance)   2) proportions of C are different from proportions of B (=independance)   3) proportions of A are NOT different from proportions of B (=NO independance)   Is there something like post-hoc contrasts for chi-squared, so I can test all my hypotheses in ONE modell and not have alpha inflation? Intuitively I would test C against a composite of A and B but that still doesnt test if A is different/independent There is a post on Stackoverflow regarding this topic, but there wasn't an answer.


title: Document classification with Bayes

text: I want to build a document classifier in R, using the Naive Bayes approach. Here are steps, that I've done so far:   * I have corpus with about 30 documents from 2 authors (Classes are: "target author" and "other author").   * "Vocabulary" (training set) has been pre-processed (removed numbers, removed punctuation, words to lower case, removed stop words, stem documents, strip whitespace), and I am considering only frequent words (top 700).   * Now I have matrix which looks like: ![enter image description here](http://i.stack.imgur.com/1bVnX.png) Then I trained my classifier using Bayes using some existing R library, e1071. Here are my questions: I want to test my classifier on other documents that were not part of the training set.   * How to prepare my data matrix? What if those other documents don't contain all the words (attributes) from my training set? Should I put dummy columns there (e.g., with `value=0`)?   * Does the position of the words (columns order) matter?  Here is an example: Training attributes:               "wild"  "wind"  "woman"      Testing attributes:               "woman" "wind" "wild"        Is this ok, or should columns be in the same order as in training matrix?


title: How to define what a "sample" is?

text: If I give you three numbers that are independently and identically drawn from a standard normal distribution, then have I given you three samples or one sample?   If the answer is one sample, then is there a short name for what I have given you three of?


title: Getting started with cluster analysis in R

text: I have a huge dataset which contains 20 columns and many rows. I have done clustering in SAS, Knime and SPSS, but I am new to R. I have to do clustering on my dataset. I have imported my data into R.   * What are some suggestions for getting started with cluster analysis in R?


title: Bootstrap confidence intervals on parameters or on distribution?

text: Excuse what may be an obvious question about bootstrapping. I got sucked in the Bayesian world early and never really explored bootstrapping as much as I should have. I ran across an analysis in which the authors were interested in a survival analysis related to some time to failure data. They had about 100 points and used regression to fit a Weibull distribution to the data. A result of this they obtained estimates of the scale and shape parameters. A very traditional approach. However, they next used bootstrapping to sample from the original data set and, for each new sample, performed a regression and came up with a new Weibull distribution. The results of the bootstrapping was then used to construct confidence intervals on the survival distribution. My intuition is a bit conflicted. I'm familiar with bootstrapping confidence intervals on parameters, but not seen it used for constructing distribution confidence intervals. Can anyone point me toward a reference/source that might provide some insight? Thanks in advance.


title: GLMM - test of significance

text: _Adapted from myprevious, unanswered question:_ I'm analyzing the results of a hormone manipulation experiment. I measured a number of variables at three times in three groups. The groups are different sizes and not all individuals were measured every time, so I'd like to use GLMM rather than a repeated-measures ANOVA. I created the model then tested the significance of the terms (time, treatment, and time x treatment) with ANOVA. I'm quite new to GLMM, but after doing the tests, further reading suggests that my approach may be inappropriate, particularly with small data sets (I have ~seven animals per group). One reason given against using this method is that it seems that there is disagreement about what the degrees of freedom should be. Is this method of testing the significance of factors in the model appropriate? I'm aware that the place for me to go is probably the Pinheiro and Bates book, but I don't have access to it at this time. Thanks in advance for any advice.               library(nlme)          datums<-data.frame(id=rep(1:20,each=3),var1=runif(60,4,6),var2=runif(60,25,30),var3=runif(60,0,1),var4=runif(60,10,15),var.time=rep(1:3,times=20),var.treatment=rep(c('a','b','c'),each=20))          datums$var.time<-as.factor(datums$var.time)     datums$id<-as.factor(datums$id)               #and now the GLMMs on each variable - I'll show just one here          var1.glmm<-lme(var1~var.time + var.treatment + var.time*var.treatment, data=datums, random = ~1| id)          summary(var1.glmm)          anova(var1.glmm)


title: After the fact hypothesis testing

text: I'm not really a statistician, but rather in need of statistical guidance, so I hope this is not too much of an off-topic question. I'm writing a master's thesis (computational linguistics/NLP), and I've got several result sets I'm comparing. Now, I didn't really formulate a null and an alternative hypothesis before I ran the experiments, which I understand means that ideally I shouldn't really go about using the T-test on my datasets. But some of the different results are so close that I've given in to the temptation and T-tested them. How unclean is this? Is it bad enough that I should leave it out of my thesis entirely, or is it less severe? In case it matters, the data are the error rates of different language models, by 10-fold cross-validation.


title: Binary interaction terms using lmer

text: I am trying to create a model using the `lmer` function. The model will contain the continuous response term "Average.profit" and explanatory terms "Type", "OtherType," and "Game" all 3 of which are binary. I also would like to include the two way interaction terms, these should have 4 possible levels `{00,01,01,11}` all of which I am interested in, As I would like to know if any of the 4 combinations increase or decrease "Average.profit". The model is below.               > m1<-lmer(Average.payoff~Game+Type+Others.Type+Type:Others.Type+Game:Others.Type+                           Game:Type+(1|Subjects),REML=FALSE, data=Subjectsm1)          > m1         Linear mixed model fit by maximum likelihood          Formula: Average.payoff ~ Game + Type + Others.Type + Type:Others.Type +           Game:Others.Type + Game:Type + (1 | Subjects)             Data: Subjectsm1             AIC   BIC logLik deviance REMLdev          54.47 69.67 -18.24    36.47   49.83         Random effects:          Groups   Name        Variance Std.Dev.          Subjects (Intercept) 0.025177 0.15867           Residual             0.122703 0.35029          Number of obs: 40, groups: Subjects, 20              Fixed effects:                            Estimate Std. Error t value         (Intercept)          0.6545     0.1746   3.749         Game1                0.2492     0.2088   1.193         Type1               -0.1675     0.2088  -0.802         Others.Type1        -0.5242     0.2088  -2.511         **Type1:Others.Type1  -0.2983**     0.2261  -1.319         Game1:Others.Type1   0.4500     0.2685   1.676         Game1:Type1          0.1733     0.2261   0.767              Correlation of Fixed Effects:                     (Intr) Game1  Type1  Oth.T1 T1:O.T G1:O.T         Game1       -0.718                                            Type1       -0.718  0.435                                     Others.Typ1 -0.718  0.531  0.435                              Typ1:Oth.T1  0.389 -0.108 -0.542 -0.542                       Gm1:Othr.T1  0.461 -0.643 -0.129 -0.643  0.000                Game1:Type1  0.389 -0.542 -0.542 -0.108  0.000  0.000        From my understanding the boldsection shows that when Type=1 AND OtherType=1 the response variable decreases by -0.2983 compared to any of the other combinations of type and others type ie. 00,01 or 10. How then do I know about the effect of any of the other combinations of Type and OtherType ie. Does 01 result in a higher response variable than 10? If this can't be shown from this model what information can I infer about the interaction between Type and OthersType?


I understand that in the dual form of the model for support vector machines, the feature vectors are expressed only as a dot product. Mapping the feature vectors to a higher dimensional space can accommodate classes that are not linearly separable in the original feature space, but calculating this mapping and working with the higher-dimensional feature vectors is computationally prohibitive. Instead, kernels can be used to efficiently calculate the same value as the dot product of the mapped vectors. How do support vector machines avoid overfitting? Is maximizing the margin of the decision boundary the only trick that they use, or am I missing something?


title: Can a trend stationary series be modeled with ARIMA?

text: I have a question / confusion about stationary series required for modeling with ARIMA(X). I am thinking of this more in terms of inference (effect of an intervention), but would like to know if forecasting versus inference makes any difference in the response. **Question:** All the introductory resources I have read state that the series needs to be stationary, which makes sense to me and that is where the "I" in arima comes in (differencing). What confuses me is the use of trends and drifts in ARIMA(X) and implications (if any) for stationary requirements. Does the use of either a constant/drift term and/or trend variable as an exogenous variable (i.e. adding 't' as a regressor) negate the requirement of the series being stationary? Is the answer different depending on if the series has a unit root (e.g. adf test) or has a deterministic trend but no unit root? OR Does a series always have to be stationary, made so via differencing and/or detrending before using ARIMA(X)?


title: Comparing Hazard Rate Ratios for different outcomes on the same sample

text: Suppose I have survival data coming from a cohort study. The outcome of interest is death and a subject can die either because of cancer or because of other causes. The two outcomes under study are therefore mutually exclusive. Now, suppose that I'm interested in assessing whether the magnitude of the associations between the exposure (for simplicity let's imagine it's dichotomous) and the outcomes is different for the two outcomes considered. The measure of association could be, for example, the Hazard Rate Ratio. I have found this paper by Lunn and McNeil that provides a possible solution to this problem. For those of you who don't have access to jstor or don't want to read it, here's the idea behind it: I duplicate my data (since in this case I have 2 possible outcomes), and create an indicator variable $\delta_i$ which represents the failure type for the subject $i$ (let's say $\delta_i=0$ for cancer death and $\delta_i=1$ for other cause of death). This means that each subject is included twice in my dataset, for example:               subject time status failure_type exposure     i       t_i  1      0            e_i     i       t_i  0      1            e_i      as the 2 outcomes are mutually exclusive, for each subject, the `failure_type` vector can be either $(0,0)'$ (censored) or $(0,1)'$ (died because of other causes) or $(1,0)'$ (died because of cancer). The time to event/censoring (`time`) and the exposure (`exposure`) are the same for the two records of each subject. At this point I fit a Cox regression model (using a robust variance/covariance estimator): $$\lambda_i(t) = \lambda_0(t)\exp(\beta_0*failuretype_i+\beta_1*exposure_i+\beta_2*failuretype_i*exposure_i)$$ and test $\hat{\beta_2}$ against the null $H_0: \beta_2 = 0$. (Of course I can also relax the assumption of proportionality of the hazards by stratifying over `failure_type`, but the bottom line doesn't change). Do you know of other possible ways of performing this analysis? I've read a paper where the authors fit two different models (one for each cause of death) and then use an Hotelling test to compare the two HRRs, but they give no references or additional information and this method seems unsound to me.


title: Calculating confidence intervals for prevalence differences

text: I would like to calculate the confidence intervals for prevalence differences. If I have the CIs for the prevalence, can I simply subtract the upper limit of one Ci from the upper limit of the other, and do the same for the lower limits as well? Or is calculating the CIs more complicated than this?


title: Multiple FDR corrected experiments using the same data

text: I am testing many (500,000) genetic variants, and the tests are FDR corrected and give me a q-value. Normally I would just call everything with q < .05 significant. But in this case I am testing those same genetic variants in two other related experiments (not using exactly the same individuals, but the samples may overlap). What to do? Would changing the significance threshold for q to .05/3=.0167 be an option? Many thanks!


title: Trying to fit a model after detrending

text: I have data for Hydrogen Sulfide Series, see here http://www.wikiupload.com/Y4WAZJ4Z0IMTK7V I applied a Box-Cox Transformation with $\lambda =1/3$ to try to stabilize the data. I plotted a few sample PACF/ACF to show that the series is not stationary and does not demonstrate constant second order properties with time. I have eliminated a possible trend and seasonal component by assuming a model of the form $y_t = m_t + s_t + x_t$ where $x_t$ is the stochastic process I am trying to model. I removed a possible seasonal component and a fluctuating mean and I get the following ACF/PACF for my series. Does that look like something that is known? How can I fit a stationary model in R to this? Is it even stationary? Maybe the above decomposition is not really applicable in this case, i.e. the time dependence is more complicated. Here are the resulting ACF-PACF after detrending ![ACF & PACF](http://i.stack.imgur.com/tjQnt.png)   Larger image


title: SVM without offset

text: I would like to know if the linear-SVM-without-offset solver: $$\min \frac{1}{2}\|w\|^2+C\sum_{i=1}^m \xi_i, \quad \mbox{s.t.}\quad y_iw^\top x_i \geq 1-\xi_i, \quad \xi_i\geq 0 \quad \forall i=1,\ldots,m.$$ can be applied to classify linearly separable data, where the hyperplane does not pass through the origin. Maybe a change of coordinate system would help? Thank you in advance.


title: Reliability test, Cronbach's alpha or else

text: I have a set of questions for my research on Importance Menu Design towards Customer Purchasing Decision. After running Cronbach's Alpha for reliability test, it turned out less than 0.5 These are the set of questions in the survey. > Psychology Effects > >   1. I read the menu thoroughly before I make an order >   2. I always order the same / familiar item without looking at the menu >   3. I always end up order more items than my initial plan after I look at > the menu >   4. The menu is attractive and persuade me to order >   5. The menu is user friendly and interested > > > Layout > >   1. The menu has clear dish classification >   2. I am attracted to big pictures on the top or left upper corner >   3. The pictures are appealing and encouraging me to order >   4. I always look at pictures before I read the menu descriptions >   5. The grid format used in the menu help me to read better > > > Dynamic of Menu Design > >   1. The menu book is sturdy and has clean appearance >   2. The menu is too heavy to hold and make me uncomfortable to make > decision >   3. The menu size is too big which make it difficult to flip over and take > too much space on the table >   4. The menu color is encourage me to order >   5. The book format is suitable to cover menu range offer by restaurant > > > Menu Information > >   1. The dish description is too complicated >   2. I prefer rounded nominal of price without currency symbol (eg. 10, 20) > rather than decimal nominal of price with currency symbol (eg. RM 9.99) >   3. I am aware with the price given >   4. I ask waiter for more details before I order >   5. The font type and size can be read easily > All the question is answered on scale 1=totally disagree, 5= totally agree Can someone help me to pin point my mistake? Did I use the wrong reliability test or something wrong with the question I constructed? Thank you. PS: I am a newbie in research


title: Creating a CI for the mean from an approximation of x-bar

text: The "classical" $100(1-\alpha)\%$-confidence interval starts from the Student's t statistic $$t=\frac{\bar{x}-\mu}{s/\sqrt{n}}.$$ Then, one obtains the desired result, e.g. $\bar{x}\pm t_{1-\alpha/2;n-1}\frac{s}{\sqrt{n}}$ by simple algebraic manipulation. Here's my problem. I am using an approximation for $\bar{x}$, which leads me to the following modified $t$ statistic: $$t'=\frac{a(\bar{x}-\mu)^2+(\bar{x}-\mu)+b}{s/\sqrt{n}}$$ where $a$ and $b$ are parameters depending on sample size, and some central moments. I would like to derive a CI for $\mu$ using $t'$, but the squared term gives me a little trouble. Is it as simple as treating the enumerator as a quadratic equation in $\bar{x}-\mu$, solving it and moving terms around? That is, let $$a(\bar{x}-\mu)^2+(\bar{x}-\mu)+b=t' \frac{s}{\sqrt{n}},$$ find the roots $r_1, r_2$ and then write, e.g., $\bar{x}-\mu=r_1$, yielding $\mu = \bar{x}-r_1$ and thus the CI is $(\bar{x}-r_1) \pm t' (s n^{-1/2})$ ? If both roots are real, how do I decide on which to pick?


title: How to display a matrix of correlations with missing entries?

text: I'd like to obtain a graphic representation of the correlations in articles I have gathered so far to easily explore the relationships between variables. I used to draw a (messy) graph but I have too much data now. Basically, I have a table with:   * [0]: name of variable 1   * [1]: name of variable 2    * [2]: correlation value The "overall" matrix is incomplete (e.g., I have the correlation of V1*V2, V2*V3, but not V1*V3). Is there a way to graphically represent this ?


title: Which terms are more standard: precision and recall or correctness and completeness?

text: I am familiar with terms like precision rate and recall rate from the machine learning literature, but recently I was reviewing a thesis which used the terms correctness and completeness. A quick Google search showed that these are used in some fields, such as image processing of remote sensing data. I'm interested to know which set of terms are more widely used, or is this just something that varies by application area?


title: Is mixed logit simulation appropriate for categorical data which include likert ratings?

text: I have data from a study where subjects were shown videos and asked to identify the contents as one of two options, as well as to rate their confidence of correct classification on a likert scale of 1 to 7. Each subject was asked to perform this classification task multiple times, and given different video material each time. I am using mixed logit because I think this data is categorical. I am encoding it as each choice situation having 14 alternatives. That is, the subject identifies the content as type 1 or type 2 then provides a rating of their confidence from 1 to 7, so I view that as 2 * 7 = 14 alternatives of which only 1 may be given as the answer. Does this seem a reasonable approach? Or is this an inappropriate way to use mixed logit estimation? I'm new to discrete choice methods, and my general statistical knowledge is rusty, so I welcome reading recommendations. Thanks very much for your time!


title: How does cross-validation overcome the overfitting problem?

text: Why does a cross-validation procedure overcome the problem of overfitting a model?


title: One sample $t$-test with range values

text: I need to compare an average value for an experimental group ($n=5$) with a reference value given by legislation. But this reference value in given as [min - max] range. Do I have to compare my average (with SD) with the two extremity of the range or could I use the mean value of the range ?


title: Doubts in linear regression

text: If a linear regression model has a constant term say 1 or 0.2, for example if the original model is $y(t) = 0.2 + ay(t-1) $, then what does this constant term imply? Will it hamper the estimates if the constant term is ignored? The question is in terms of estimation of linear models using Maximum Likelihood estimation or any other estimation technique, where in most of the examples I have seen that the parameters are estimated and not the constant terms.


title: Shape of confidence interval for predicted values in linear regression

text: I have noticed that the confidence interval for predicted values in an linear regression tends to be narrow around the mean of the predictor and fat around the minimum and maximum values of the predictor. This can be seen in plots of these 4 linear regressions: ![enter image description here](http://i.stack.imgur.com/NhMqj.jpg) I initially thought this was because most values of the predictors were concentrated around the mean of the predictor. However, I then noticed that the narrow middle of the confidence interval would occur even if many values of were concentrated around the extremes of the predictor, as in the bottom left linear regression, which lots of values of the predictor are concentrated around the minimum of the predictor. is anyone able to explain why confidence intervals for the predicted values in an linear regression tend to be narrow in the middle and fat at the extremes?


title: Need help on Gibbs sampling with truncated normal and gamma

text: I am trying to use Gibbs Sampling to simulate a random sample from a joint distribution $f(\beta ,{{Z}_{1}},...,{{Z}_{75}},{{\lambda }_{1}},...,{{\lambda }_{75}})$, where the fully conditioned distribution function are: $$\beta |{{Z}_{1}},...,{{Z}_{75}},{{\lambda }_{1}},...,{{\lambda }_{75}}\sim N\left( \frac{\sum\limits_{i=1}^{75}{{{\lambda }_{i}}{{Z}_{i}}}}{\sum\limits_{i=1}^{75}{{{\lambda }_{i}}}},\frac{1}{\sum\limits_{i=1}^{75}{{{\lambda }_{i}}}} \right)$$ For i=1,...,24, ${{Z}_{i}}|\beta ,{{\lambda }_{1}},...,{{\lambda }_{75}}$ ~ left truncated normal at 0 $${{f}_{L}}\left( t;\beta ,0,\frac{1}{\lambda } \right)=\left\\{ \begin{array}{*{35}{l}} 0 & \text{if t}\le 0 \\\ \frac{{{e}^{-\frac{{{\lambda }_{i}}}{2}{{(t-\beta )}^{2}}}}}{\sqrt{2\pi /{{\lambda }_{i}}}[1-\Phi (-\beta \sqrt{{{\lambda }_{i}}})]} & \text{if t}>0 \\\ \end{array} \right.$$ $$$$ For i=25,...,75, ${{Z}_{i}}|\beta ,{{\lambda }_{1}},...,{{\lambda }_{75}} \sim$ right truncated normal at 0 $${{f}_{R}}\left( t;\beta ,0,\frac{1}{\lambda } \right)=\left\\{ \begin{array}{*{35}{l}} \frac{{{e}^{-\frac{{{\lambda }_{i}}}{2}{{(t-\beta )}^{2}}}}}{\sqrt{2\pi /{{\lambda }_{i}}}[\Phi (-\beta \sqrt{{{\lambda }_{i}}})]} & \text{if t}\le 0 \\\ 0 & \text{if t}>0 \\\ \end{array} \right.$$ $$$$ For i = 1,...,75, $${{\lambda }_{i}}|\beta ,{{Z}_{1}},...,{{Z}_{75}}\sim {\rm Gamma}\left( \frac{5}{2},\frac{2}{4+{{({{Z}_{i}}+\beta )}^{2}}} \right).$$ I am having trouble implementing this. My algorithm: choose some fixed initial $Z_{i}$'s and $\beta$, and generate $\lambda_{i}$ from Gamma distribution.   Let $Z_{i}=0$ and $\beta=0$, I get $\lambda_{i} \sim {\rm Gamma}(5/2, 2/4)$. I then use these $\lambda_{i}$'s to generate $Z_{i}$. Next step: how do I generate truncated normal, how do I know if my t is greater than zero or not? Any experts would like to do what I am doing right or wrong? Maybe suggestions? All welcome! Thank you very much in advance, I have been very confused about this! Note: I am using Java.


title: Should I add noise to my truth data before before training my classifier?

text: **My task is to develop a system that will take in a series of measurements and return the probability that an object is a type 1, type 2,... type n.** I will refer to the system I have to create as a the classifier. I know all the possible types, so I know it must be one and only one of the types. I have to develop my system based on some data I have been given from a simulation. The physics-based simulation generates both the truth data, what is actually occurring in the simulation, and the simulated measurement data, what the simulated sensors are measuring. Unfortunately, I don't have as much data as I would like, and obtaining more data by running the physics-based simulation is not an option at this time, so I must make do with what I have. My superior believes that using a Gaussian Mixed Model (GMM) would produce the best classifier, so my solution must use GMM. His proposed solution is to use the Expectation-Maximization Algorithm to develop the classifier, only using the truth data. To measure how likely a given measurement is to be a given type, we measured the Mahalanobis Distance between the mean and covariance given by the GMM and the measurement. The type with the smallest Mahalanobis Distance, we thought would be the most likely type of the object. Once we went through our initial testing, we could measured the effectiveness of our classifier by plugging in the simulated measurement data and seeing how we did. (We simply looked at which type the classifier decided was most probable). Our initial performance was dismal, we found that nearly everything was being classified as type 1. The only time any other type was being correctly classified was when the truth and the measurement were nearly identical. After examining the data more closely, the truth data is tightly clustered for all types except type 1, causing the covariance associated with the GMM to be very small. For type 1, the covariance is large, because the particular data set I have shows a lot of variability. I realize that the "best" solution would be to get more simulation data, but that's not an option at this time, so I'm trying to figure out how else I can improve my classifier performance. My question is, should I artificially add noise into the training data to give the classifier a "preview" of the measurement data or should I try to find an alternative method of classification? I'm not a mathematician, I'm a engineer, so if I sound like I don't know what I'm talking about, its because I am in unfamiliar territory. So please let me know if I'm using the wrong terms or if I'm being unclear and I'll do my best to clarify. * * * The data has a total of _three dimensions_ , one of which is time. I can't seem to find a way to use time dimension effectively because some simulation runs will last 500 seconds and others will last only 200 seconds. I realize that I could increase the number of dimensions by taking the derivative of one of the dimensions with respect to another, or by using a nonlinear equation to create a new dimension based on one or two of the variables, but I'm not sure how helpful that would be.


title: Communicating Regression Model Results

text: I am concerned about how unequipped most people are (both within and without academia) to properly employ standard model building methods such as linear regression and to interpret the results of these models. Both from my own observations and the literature, it is clear that most people are being poorly served by the standard statistical tools they have available. To improve this situation, I would like to propose a different set of model evaluation output tables (examples shown below in Table 2 and Table 3) to replace the standard output table that are almost universally used today (an example shown in Table 1). This new output format is robust to human error (both in model specification and analysis) to a much greater extent that is our current format. I would like ask people for their critique of this proposed method and whether there was interest in using it. (Sorry for the length of this, it got much longer than expected) ## Approach Overview First, calculate an "honest $R^2$" using leave one out cross validation (LOOCV) to get an overall measure of the model's fit. Simply put, remove each data point in turn, estimate the model without that data point, then use the reduced model to calculate the value of the removed data point and use these errors to estimate your $R^2$. Then for each coefficient in the model, remove that parameter from the model and calculate an "honest $R^2$" for that submodel. Calculate the effect of the adding the coefficient to the model as the overall model's $R^2$ minus the $R^2$ of the submodel when the coefficient is removed. Thus a positive value for this coefficient $R^2$ indicates the fit improves when it is added to the model; a negative value indicates the fit decreases. Results would be reported as a simple table of honest $R^2$ values (see examples in Table 2 and Table 3 below). ## Discussion This approach is designed to be straightforward and to do three things:   1. Provide end users with the information they need to interpret the (practical) significance of models and coefficients   2. Provide end users with this information in a form they are already familiar with   3. Minimize the effect of user error both in developing the models and interpreting results I think this approach does all three. The use of the universally understood $R^2$ metric should mean it is readily understandable and can plug right into users' existing mental frameworks. Second, the use of this approach is robust: both to misspecified models and to poor interpretation. Most models depend on certain assumptions that users often fail to check are met. Of course, p-values for coefficients in linear regressions depend on certain parametric assumptions that are generally violated to a greater or a lesser extent. This proposed approach is resilient to such issues (the main issue it is susceptible to is, as with regular output statistics, correlation between observations, which could lead to an overestimation of the honest $R^2$). In regards to the interpretation of results, p-values or often misinterpreted (as has been extensively noted). I believe this honest $R^2$ approach is much more resilient to these issues as it is focused on effect size (which is often what people incorrectly take the p-value as a proxy for). It emphasizes practical significance in place of statistical significance which is a concept people have a lot of trouble with. Furthermore, this method allows the direct comparison of results to those generated by other methods (even those that don't generate likelihoods) such as machine learning techniques such as random forests. One issue with this approach is computational burden. Where LOOCV is computationally infeasible I would suggest 10- or 5-fold CV. Since these methods should result in higher error (as the model is training on smaller data sets), the honest $R^2$ values reported by them would be _conservative_. Therefore they can be used and reported in placed of the LOOCV with any risk of mis-comparison being of a conservative nature. Another issue is the LOOCV estimate is known to have high bias. I'm not really sure how this could be dealt with or if it is a serious problem. One last point is the LOOCV is asymptotically equivalent to AIC, so this could fit into that paradigm a bit. ## Applied Example Taking a set of housing data ( http://archive.ics.uci.edu/ml/datasets/Housing ) we try to predict average house value in a suburb based on average house age, number of rooms in the house, pollutant levels, pupil to teacher ratio, and proximity to a highway. We start with a linear regression. Below is similar to what most software will currently output (this specific table was generated by R). ### Table 1. Standard regression table.               Coefficients:                    Estimate Std. Error t value Pr(|t|)         (Intercept)     7.76739    4.98881   1.557 0.120112         AGE            -0.01509    0.01378  -1.096 0.273773         ROOMS           7.00565    0.41172  17.015    2e-16 ***     NOX           -13.31418    3.90262  -3.412 0.000698 ***     PUPIL.TEACHER  -1.11645    0.14799  -7.544 2.17e-13 ***     HIGHWAY        -0.02487    0.04257  -0.584 0.559341         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 5.819 on 500 degrees of freedom     Multiple R-squared: 0.6037, Adjusted R-squared: 0.5997       The more honest statistics results are shown below. Note that these results aren't surprising given what we saw in table 1. The parameters that were not significant in Table 1, result in worse models as measured by the change in the honest $R^2$ (again the coefficient Honest $R^2$ value indicate how the honest $R^2$ changed when that coefficient was added to the model). ### Table 2. Proposed "honest" statistics table.                         Item Coefficient Honest.R2         -Full Model-              0.593080          (Intercept)      7.7674 -0.000998                  AGE     -0.0151 -0.000409              ROOMS      7.0056  0.228901                NOX    -13.3142  0.008123      PUPIL.TEACHER     -1.1165  0.045506            HIGHWAY     -0.0249 -0.002018      Finally, this approach allows us to directly compare more exotic algorithms such as Random Forests. Unlike Bayes factor, AIC, BIC or some other methods, likelihoods are not required for comparison. Anything that creates predictions can be compared. Classification algorithms can also fit into this scheme if you use one of the common pseudo $R^2$ approaches for calculating $R^2$. ### Table 3. Comparison between algorithms. (In this case Support Vector Machines and Random Forests)                       SVM               |         Random Forest               Item Honest.R2  |             Item Honest.R2       -Full Model-  0.71052   |     -Full Model-   0.7450                AGE  0.00834   |              AGE  -0.0128              ROOMS  0.34722   |            ROOMS   0.2068                NOX  0.02338   |              NOX   0.0654      PUPIL.TEACHER  0.01813   |    PUPIL.TEACHER  -0.0127            HIGHWAY -0.00381   |          HIGHWAY  -0.0255      Side note: It is interesting to interpret these results as we can see in the linear model pupil per teacher has a statistically significant effect on housing prices. While in the random forest model, it does not. Since the random forest model is the most predictive of the models, I would have to conclude that this was evidence that this coefficient did not significantly (practical significance) affected housing prices. This is a small illustration of the weakness of blindly applying linear models to everything and using that to carry out hypothesis testing. ## Implementation I have uploaded (rough-hewn, unoptimized) code in R to implement this technique at: http://dl.dropbox.com/u/94002/HonestStats.R If there is interest, I can optimize this to be much faster for linear models and implement things like k-Fold CV.


title: Model averaging in prediction -- "Wisdom of the Crowd"

text: Suppose I'm trying to predict $Y$ (a real number) and I have $n$ experts with guesses $Y_1,...Y_n$. Each prediction is a reasonable guess as to the value of Y in itself (hence the name "expert"), but we should be able to come up with a better prediction for $Y$ by combining the experts' guesses. We could take our guess $Y_o$ to be $$Y_o=\sum_{i=1}^n\frac{Y_i}{n}$$ but this isn't optimal if some experts are better than others. I have a few questions related to this topic. Question 1: Is this equal weights method used much? Are there any other simple methods? Can you refer me to (preferably relatively basic) information on the topic? Question 2: I've seen model averaging discussed in this as a possible alternative to other methods of model selection for linear regression which may fail, for example, if there is multicollinearity. That's led me to wonder whether the following solution might be justifiable: We want (perhaps out of naïveté) to regress $Y$ against predictors $X_1$, $X_2$ and $W$. Suppose, in fact, that $X_1$ and $X_2$ themselves are predictors of $Y$. Using multilinear regression for $Y$ with all three variables is problematic because, being predictors of $Y$, $X_1$ and $X_2$ exhibit multicollinearity, but what about creating two multilinear regression models $$Y=\alpha_o + \alpha_1X_1 + \alpha_2W$$ and $$Y = \beta_o + \beta_1X_2 + \beta_2W$$ and then using some average-like combination of these two regression-based predictors to predict $Y$? What weights might we use? Again, is there any literature on this topic? Thank you. EDIT: Since I posted this, I've read about stacked regression, as described here. Does the situation I've described sound suitable for applying stacked regression?


title: Is the Binomial Effect Size Display (BESD) usually an unsound, misleading technique?

text: It's hard for me to accept that Donald Rubin would ever come up with a true lemon of a technique. Yet that's my perception of the BESD [1, 2], 3]3. He and Robert Rosenthal (1982*) claimed that there was value in showing "how to recast any product-moment correlation into such a [2x2] display, whether the original data are continuous or categorical." The table below is from p. 451 of the 2nd link above: ![enter image description here](http://i.stack.imgur.com/2ac1b.jpg) This technique seems to overstate the magnitude of almost any effect size. Here, RSQ from the original data = .01, but when "translated" into a 2x2 crosstab, we seem to be faced with a much stronger effect. I don't deny that, when the data are recast into categorical format in this way, phi indeed = .1, but I feel something has been very distorted in the translation. Am I missing something truly valuable here? Also, I have the impression that in the past 10 or so years the statistical community has by and large rejected this as a legitimate method; am I wrong on that? *Rosenthal, R., & Rubin, D.B. (1982). A simple general purpose display of magnitude of experimental effect. Journal of Educational Psychology, 74, 166–169.


title: Similarity measures for point processes

text: I have multiple measurements of a point process: vectors of 0's and 1's. I'm trying to gauge the similarity of the measurements, but have no idea how to proceed. Any suggestions? Thanks!


title: Relationship between correlation and sample variance

text: A correlation coefficient is lower if there's a low variance in the characteristic of the sample. For example, the correlation between IQ and school achievement follows this pattern. The correlation is lower if you only include students with similar school achievement. It gets higher if you include students from very different school types. Is there a (mathematical) explanation?


title: Heading in plots by mvpart()?

text: How can I put in a heading in the tree output when using `mvpart()`? I am using `xv="1se"`. e.g. using the spider data in mvpart fit<-mvpart(zora.spin~water+sand+moss+reft+twigs+herbs,data=spider) from this I get a tree in R. And I am wondering how to "put in" a heading in this figure.


title: Computing document similarity in latent semantic analysis

text: I have a question regarding Latent Semantic Analysis - after performing SVD decomposition of term-document matrix and choosing some number of dimensions, I get the set of new document vectors. Now, how can I calculate similarity between two documents? New document vectors contain negative values, and results produced by cosine similarity make no sense.


title: more variables than cases in PCA/factor analysis

text: I'm trying to run a PCA on 15 variables for which I have a sample size of 10 (i.e. 15 columns x 10 rows, all continuous data). I have two problems: 1) more variables than cases and 2) some highly correlated variables (r > 0.9). I don't want to make arbitrary decisions about excluding variables prior to the PCA, but if I go ahead with them all included, are my results invalid? I've done lots of reading about singular matrices and matrices that are not positive definite, etc, but still can't get a straight answer on whether the results are completely invalid or not. Also, the most highly correlated variables are the ones I am most reluctant to remove from a scientific/functional point of view.


title: Reparameterizing the binomial link for psychometric data

text: Consider an experiment where a square and diamond appear on a screen, one before the other, and participants are required to judge which came first. Manipulating the time interval between the two shapes and averaging responses to obtain the probability that participants respond "Square came first" typically yields s-shaped curves like: ![enter image description here](http://i.stack.imgur.com/1a0sf.png) Now, I know that averaging responses to proportions within each of multiple participants is a bad idea, so presumably I'd want to analyze data from such experiments using a generalized mixed effects model with a link function that recognizes the binomial nature of the raw response data. However, I find the logit and probit links dissatisfying because they have a slope+intercept parameterization that confounds properties of the curve that have independent psychometric interest. Specifically, while the slope is indeed a parameter of interest (as it indexes perceptual sensitivity), the intercept is not. More interesting than the intercept is the shift of the curve, usually indexed by the point on the x-axis at which the curve hits p("square first")=.5. To elaborate, consider an extension of the above related experiment where there are three groups of participants:   1. A "control" group with no distractions and equal proportions of "square-first" and "diamond-first" trials   2. A "distracted" group listening to distracting music   3. A "biased" group who experience more "square-first" trials than "diamond-first" trials Presumably, the control group and distraction group should differ in terms of their slope, evaluated via:               #fit a model that contains no group-by-slope interaction     m1 = lmer(         response ~ (1|participant) + interval + group         , family = binomial         , REML = FALSE     )     #fit a model that contains the groups-by-slope interaction     m2 = lmer(         response ~ (1|participant) + interval*group         , family = binomial         , REML = FALSE     )     (AIC(m1)-AIC(m2)) #AIC-corrected log-base-e likelihood ratio      Furthermore, presumably the biased group will differ from the other groups on the shift of the function. _If_ I am willing to assume that the groups _don't_ differ in terms of slope, then the shift can be assessed in a manner similar to above by evaluating a group effect on the intercept:               #fit a model with no group effect on the intercept     m3 = lmer(         response ~ (1|participant) + 1         , family = binomial         , REML = FALSE     )     #fit a model with a group effect on the intercept     m4 = lmer(         response ~ (1|participant) + group         , family = binomial         , REML = FALSE     )     (AIC(m1)-AIC(m2)) #AIC-corrected log-base-e likelihood ratio      However, this approach to evaluating effects on shift breaks down if you cannot assume that slopes are not also affected by your predictor variable. Sure you can evaluate the assumption of no-effect-on-slopes as described above, but this not only leaves you in the lurch when you _do_ find an effect on slopes, it is susceptible to the same criticisms of other approaches to assumption tests. Finally, even with the assumption of equal slopes, presumably this is an assumption that the slopes are equal _on average_ ; variability in the slopes will cause analysis-of-intercepts-as-a-proxy-for- shift to have lower power relative to a direct analysis of shifts (because the variance of the intercepts will represent the combined variance of the slopes and shifts). So, my questions are:   1. Is there a link function that handles binomial data while also implementing the slope+shift parameterization?   2. Failing a "yes" answer to #1, any suggestions on how might I go about coding such a link function? I know how to code an MLE search given a single participant's response data across intervals in a single experimental condition, but I'm not sure how to scale this up to a link function that could be used to analyze data across groups of individuals and conditions. P.S. I am aware that, given a set of maximum likelihood estimates for the slope and intercept parameter values, it is simply a matter of arithmetic to compute the ML estimate for the shift value. However, this does not pertain as it cannot be applied to the inferential procedures described above.


title: Can we run a chi squared test on a normalized function?

text: Hi I am fit a maxwell distribution and attempt to find the chi squared value in two cases: >   1. When the data is normalized. >   2. When the data is un-normalized. > My problem is that the two different cases give completely different values of chi-squared! So my question is: is it mathematically correct to use the chi squared test on normalized data? Notes on my data: Number of counts: 41,000 Expected fit: Maxwell distribution. Method to calculate Chi-squared: Find the best maxwell fit to data, then use scipy.stats.chisquare function in python programming language to calculate the chi-squared value using the experimental and expected counts. Chi-squared documentation For reference: http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.chisquare.html#scipy.stats.chisquare I have attached the two graphs here for you to see the difference in chi squared values:   1. The un-normalized graph: please ignore the y-axis which says "Normalized" ![Please ignore the "Normalized" on the y-axis here](http://i.stack.imgur.com/RkdFG.png)   2. The un-normalized graph: ![enter image description here](http://i.stack.imgur.com/PSXeT.png)


title: PyMC: Which parameter is passed to the model by sampler? How to find the value passed by custom sampler?

text: I am relatively new for PyMC and as a whole I have some less knowledge in probability. As far as I understand MCMC algorithm, its main task is to propose samples of parameters to the model (in iterations) until the probability of model is maximized (Please do correct me if my understanding is wrong). For example let me define a model as \-----------------------------------testmodel-----------------------------------------------               import numpy as np     import scipy.stats as ss     from pymc import *     N = 50     data = np.array([ss.gamma.rvs(2, scale=50) for i in range(N)])     labels = np.random.randint(0, 2, 50)     @ stochastic     def logprob(samples=[2, 50, 2, 50],  data=data, labels=labels, value= -10000):         a, b, c, d = samples         data1, data2 = data[(labels == 0)], data[(labels == 1)]         pdf1 = ss.gamma.logpdf(data1, a, scale=b)         pdf2 = ss.gamma.logpdf(data2, c, scale=d)         pdf1[(pdf1 == -np.inf)] = -500         pdf2 [(pdf2 == -np.inf)] = -500         return np.sum(pdf1)+np.sum(pdf2)      ## data is the observation we have and it is known that it follows a gamma distribution. Lets suppose all four the parameters follows the normal distributions as a, c~N(2,0.2) b, d~N(50, 5) I want to create a custom step method that proposes the random variable drawn from these distribution. I tried to create some but ran into error while application.               class StandardNormal(Gibbs):         def __init__(self, stochastic, verbose=None):                  Gibbs.__init__(self, stochastic, verbose=verbose)              def step(self):             self.stochastic.value = [ss.norm.rvs(2, scale=0.2),                                      ss.norm.rvs(50, scale=5),                                      ss.norm.rvs(2, scale=0.2),                                      ss.norm.rvs(50, scale=5)]     import testmod     M=MCMC(testmod)     M.use_step_method(StandardNormal, samples)      it says 'samples' are not defined. If I define the 'samples' in the model as given normal distributions, are the traces of the 'samples' the same as proposed by sampler ? I tried without creating custom step method it worked, But I would like to get the values of the samples proposed to the model by default sampler. As next part of my program will need those values, that maximizes the pdf of the data model to classify the two data.


title: Where to find csv files of census and IRS data

text: Google searches for "census data by zipcode" have come up mostly empty of anything useful and I'm going in circles on the census's website, so I'm wondering where people commonly find this data for use in their analyses.


title: Estimating hidden transfers of market share

text: Suppose we have yearly data representing the market share of three companies, say A, B and C. In other words, we have observations: $$ A_t, \; B_t \;\; \text{and} \;\; C_t \;\; \text{where} \; \; A_t+B_t+C_t = 1 $$ for $t = 1, \dots,T$. Suppose that in year $t$ the market share of company A has changed by $\Delta A_t = A_t - A_{t-1}$. Is there any way of estimating how that change can be sub-divided into market share lost to or acquired from companies B and C? My actual problem includes 5 companies, but I guess that the solution shouldn't change too much.   Thanks!


title: Estimating cumulative effects of lags in ARIMA transfer Functions

text: What is the difference between and AR(1) term in an arima results with only and advertising spend regressor with 1 lag and a Koyck specification with only an advertising spend variable with 1 lag. If I specified a model in autobox or auto.arima with only adspend and it's lag as predictors, would I calculate the cumulative effects of adspend the same way as a Koyck or Autoregressive DL model?


title: Non-Orthogonality in PCA?

text: i) What is the main role of "only" trying to find orthogonal components in PCA? I can understand, that we would not want a zero-solution as well as find directions that are orthogonal in order to explain most of the variance. When we look at the problem in terms of finding a projection matrix that preserves the gram matrix of mean centered data, ii) What would an optimization under a non orthonormal constraint produce, as long as we make sure the solution is non-zero? Why would it be or not be useful? I am aware of the concept of non-orthogonal eigen-functions as well, just in case you want to generalize the answer to a Kernel PCA, that would be fine as well.


title: Bayesian Information Criterion (BIC) for large samples

text: The Bayesian information criterion is defined as $BIC = -2 \text{ln}(L) + k\text{ln}(n)$, where $L$ is the maximized likelihood of the data, and where $n$ is the sample size. In case of a huge sample size, BIC tend to $\infty$. Is there any transformation that needs to be done in order to compute the BIC for large samples? Thank you,


title: What percentage of values (compared to response) must be present in a variable to qualify as a predictor

text: Consider a sparse dataframe with lots of NA and a large number of explanatory variables and a response variable. The response variable doesn't contain any NA. I want to clean this dataframe by discarding some of the explanatory variables having most of their entries NAs. Is there any study that suggests the minimum percentage of non-NA values that must be present in a variable to qualify it as a predictor? For example, if the total percentage of non-NA values in column X is less than 5% of total records, then drop column X.


title: Weka batch filtering on random projection

text: I am new to Weka. I have used the following statement in Weka CLI:               java weka.filters.MultiFilter -b -F     "weka.filters.unsupervised.attribute.Remove -R 4,62,83,85,86" -F     "weka.filters.unsupervised.instance.RemoveWithValues -S 0.0 -C 9 -L 7 -V" -F     "weka.filters.unsupervised.instance.NonSparseToSparse" -F     "weka.filters.unsupervised.attribute.NumericToNominal -R first-83" -F     "weka.filters.unsupervised.attribute.RandomProjection -P 80.0 -R 42 -D Sparse2" -     "weka.filters.unsupervised.attribute.Wavelet" -F     "weka.filters.supervised.attribute.PLSFilter" -i train_and_test.arff -o train_and_test_o.arff -r test.arff -s test_o.arff      the `train_and_test.arff` is my train set, but I appended my test data of `test.arff` at the end of it to see the result. The statement runs without error message. I expect to see the "appended" test data at the end of train set file output same as that of test set output file. However, it turns out not the same. Do I make some error? Or do I have a miscoception about "batch filtering" in Weka?


title: Prey preference - can I use a paired t- or wilcoxon test?

text: I have done an experiment in which I looked at the prey preference of an aquatic invertebrate. I used 10 containers all with one predator but with different densities of prey. I used two different prey species of which the densities were similar for every container. For example container 1 contained 5 individuals of prey A and 5 of prey B. Container 2 contained 10 individuals of prey A and 10 of prey B and so on. Now I want to see which prey is prefered by the predator. My data looks as followed: > prey.a <\- c(5,5,10,10,9,12,13,8,15,17) ## number of captured prey (A) > > prey.b <\- c(5,5,4,2,8,8,4,9,4,6) ## number of captured prey (B) > > prey.density <\- c(10,10,20,20,20,30,30,30,40,40) ## total number of prey > available (A and B) My initial thought was to use a paired t-test. The data violate the assumption of independent samples because they are kind of paired that is why I thought a paired t-test. However, my data violate the assumption of normality so perhaps a Wilcoxon test for paired measurements would be suitable. The thing is that both these test are designed for a before and after treatment which I don't have. Thus, I'm not really sure if I can actually use these tests. The problem is tha it's a paired rank test of the differences and my low density treatments will always be ranked lower than my high density treatments. My question is; What would be a suitable test to analyse these data? I also had this wild thought about a GLMM in which I use the container as a random effect and add prey species as a fixed effect, nevertheless, normality. btw I used different prey densities because I also looked at the effect of prey density on the total number of captured prey, which has nothing to do with my question.


title: Features selection by filter methods for multivariate time series

text: I have a data set in which the samples are multivariate (about 30 variable/features) time series. These samples refer to two classes. I would like to select the variables more relevant to discriminate between the classes and also to remove redundant variables. I am thinking about a filter method for feature selection (I am not interested on a specific predictor), but I have never worked with time series. I know Fisher score and AUC method, but I don't know how I could apply to these data. Anyone could provide me with some help? thank you


title: What would be a reliable way to check if a square matrix is computationally singular?

text: To diagnose computational singularity, one could, for example, check rank deficiency or compute singular-value decomposition. I have the following questions: (1). Do the two methods listed above generally produce results that coincide with each other? (2). If the answer to (1) is negative, which method is considered more reliable? _**ADDITIONAL INFO:_** I don't think I should've used the word "reliable" here. I am writing a C++ function using the Armadillo library. The function will be called in `R`. Part of the C++ function uses the Armadillo function `solve()`. Under singularity, the function would throw an exception. But I would like to catch the problem before `solve()` is executed and output a more meaningful error message to R. So I was thinking perhaps I could check rank deficiency or svd to accomplish that.


title: When actual analysis need to be done differently from designed

text: For example, you calculate the sample size with given power and significant level for a, e.g. two-way balance ANOVA design study. But at the end, its assumptions are not met and cannot be remedial. One switch to a more suitable method, e.g. Poisson regression and there are debatable in calculate power for the new method (you do not know the true power now). How can one explain about the power of the study now? with what evidences?


title: What is the expected number of steps before the spider exits? Markov Chain

text: A little spider lives in a rectangular box of which the sides are 3 and 4 cm long. It can only sit in one of the four corners marked with the numbers 1,2,3,4 (clockwise). Assume in Corner 2 there is a bigger spider ready to eat the little spider and in Corner 3 there is a hole leading to the outside through which the spider can escape. If the spider in Corner 1, what is the expected number of steps before the spider exists? Here is the Markov chain transition matrix $${\begin{bmatrix} 0 & 15/47 & 12/47 & 20/47\\\ 15/47 & 0 & 20/47 & 12/47\\\ 12/47 & 20/47 & 0 & 15/47\\\ 20/47 & 12/47 & 15/47 & 0 \end{bmatrix}}_{4\times 4}$$ In my notes my lecturer says we have to solve               e1 = 15/47+12/47+20/47(1+e4)     e2 = 20/47+15/47(1+e1)+12/47(1+e4)     e3 = 20/47+12/47(1+e1)+15/47(1+e4)     e4 = 12/47+15/47+20/47(1+e1)      I don't really understand. Anyone can explain why?