title: Identifiability in nonlinear mixed-effects models

text: I am interested in the identifiability of linear mixed effects models. Let's assume $p$ subject are observed at different instants in time. Let $\mathbf{y}_{i}$ $(1 \leq i \leq p)$ the vector containing the $n_{i}$ observations of the $i$-th subject. We assume the following model : $$ \mathbf{y}_{i} = \mathbf{f}_{i}(\beta,\mathbf{b}_{i}) + \varepsilon_{i} $$ where $\mathbf{f}_{i}$ is a known function, $\beta$ is a vector of unknown fixed parameters, $\mathbf{b}_{i}$ (resp. $\varepsilon_{i})$ is a random vector whose distribution is assumed to be multivariate normal with mean $0$ and covariance matrix $\mathbf{D}$ (resp. with mean $0$ and covariance matrix $\sigma^{2}\mathrm{I}_{n_{i}}$, independent from $b_{i}$). I searched on Internet and I didn't find much about identifiability for such models. All I found is this article : When is a nonlinear mixed-effects model identifiable ? but I'm a bit confused because the article doesn't seem to make a difference between fixed and random effects. I also found a document in which they write : > Let $\mathcal{L}(\mathbf{y}_{i} \vert \mathbf{b}_{i}) = > \mathcal{N}(\mathbf{f}_{i}(\beta,\mathbf{b_{i}}),\sigma^{2}\mathrm{I}_{n_{i}})$. > By the deconvolution theorem of Kagan, Linnik and Rao (1973), as long as > $\mathbf{f}_{i}(\beta,\mathbf{b}_{i})$, as a random variable for any fixed > $\beta$, does not contain any Gaussian components, and > $\mathbf{f}_{i}(\beta,b)=\mathbf{f}_{i}(\beta,b')$ implies $b=b'$, then > $\beta$ is identifiable. Are there other documents dealing with this problem ? Do you know of other ways to test for the identifiability of such a model ?


title: Testing the random slope with correlated random effects

text: I have a mixed/random effects model $$\mathbf{y}_i=\mathbf{X}_i\boldsymbol\beta+\mathbf{Z}_i\mathbf b_i+\boldsymbol\epsilon_i,$$ where random effects $\mathbf b_i$ has variance- covariance matrix $$D=\text{cov}(\mathbf b_i)=\left(\begin{array}{cc} \sigma_{1}^{2} & \sigma_{12}\\\ \sigma_{12} & \sigma_{2}^{2} \end{array}\right).$$ I would like to test $H_0:\sigma_2^2=0$. I know that under the null $\sigma_{12}=0$. But I am not willing to restrict that $\sigma_{12}=0$ under the alternative, i.e., $\sigma_{12}$ and $\sigma^2_1$ are unspecified. I know there are some tests developed assuming independent random effects. For example, the `RLRsim` package can do the simulation-based exact (restricted) likelihood ratio test for testing the presence of a single variance component. Besides, according to @StasK, "the bootstrap on the boundary of the parameter space has its own set of issues and problems leading to inconsistency." So I wonder if there is any potential solutions or references for the individual variance component test without assuming independent random effects. It would be perfect if it is readily available or can be easily implemented. Thanks!


title: Quantile regression standard error and the OLS standard error

text: I was asked by a non-statistician what is the driving factor for the standard error estimates of the parameters. Here are my thoughts:   1. From asymptotic theory, both OLS standard error and Quantile regression standard error depends on the error/residual distribution. The OLS standard error mainly depends on the variance of the error distribution, while the quantile regression depends on the density at the $\tau$th quantile.   2. More often the iid assumption does not hold for the residuals, so bootstrap methods are used. There is no explicit formular to explain what should be driving the standard error estimation in this case. But it should perform similar to the iid case when there are lots of observations. Anyone has a more right, clear and intuitive explanation?


title: Can you see any S-curve in the scatter plot

text: Can you see any S-curve in the below scatter plot? If there is the S-curve is it correct to use a scatterplot plot? If the linear correlation is 0.85, does it mean the S-curve is unlikely? I am trying to understand the relationship between the two variables. ![enter image description here](http://i.stack.imgur.com/4L0Dr.png)


title: How can I test the same categorical variable across two populations?

text: I have data that looks a little something like this:               ID         Status     01         A     02         G     03         E     ...        ...     100        G      You get the idea, I think. I have this data from two separate populations (cohorts) and I want to compare the distribution of the status variable in one population to the distribution in another. The question I'm answering is something like this: If you didn't know any better, could these be from the same population? I think this means I should perform a Person's chi squared, though I am not sure. I am also not sure on the methodology for transforming the variable in a way that allows you to run the test. (I'd especially like to know how to do this in R.)


title: How to analyze ratings/purchasing data to find purchasing drivers/probabilities?

text: Given survey data where respondents rated 10 different questions on a scale from 1-10, and whether or not they purchased (binary variable): We are trying to determine which of the survey questions are most strongly correlated to the person's purchasing behavior by calculating a probability or an odds ratio value per question. The ratings (1-10) seem to follow a bell curve distribution, but what type of analysis can/should we use to determine the correlations/probabilities? Conceptually, if we were to plot this out on a graph, what would it look like? Would we have 1 and 0 (purchased and not purchased) on the Y axis and the ratings 1-10 on the X axis?


title: Mean absolute deviation vs. standard deviation

text: In the text book _"New Comprehensive Mathematics for O Level"_ by Greer (1983), I see averaged deviation calculated like this: > Sum up absolute differences between single values and the mean. Then get its > average. Througout the chapter the term **mean deviation** is used. But I've recently seen several references that use the term **standard deviation** and this is what they do: > Calculate squares of differences between single values and the mean. Then > get their average and finally the root of the answer. I tried both methods on a common set of data and their answers differ. I'm not a statistician. I got confused while trying to teach deviation to my kids. So in short, are the terms _standard deviation_ and _mean deviation_ the same or is my old text book wrong?


title: How to estimate training time prior to training? (Machine Learning & Neural Networks)

text: I'd like to know ahead of time if my training will take 8 hours, 8 days or 8 weeks. (The 8 was an arbitrary number I chose obviously). Is there a reliable way to estimate the time it will take? Can I extrapolate the time it takes to train 200,000 as roughly double the time it takes to train 100,000? It would be helpful to be able to estimate whether it will take a couple hours or a couple days or even weeks because then I can tweak the parameters ahead of time.


title: Scaling sensor output before training a machine learning model

text: I'm working with the Emotiv EEG headset. My goal is to collect the signal from its 14 sensors while a subject who is wearing it performs a certain task, and use a machine learning model (a neural network in the present case) to identify certain patterns of interest in the signal from the device. I would like to scale the sensor data before providing it as input to the model. I'm collecting this data in separate sessions for each individual subject. I observed that initially, the sensor data varies between 7,500 and 8,500. I'm considering scaling the data across all the subjects to zero-mean and unit-variance, followed by PCA or an equivalent decorrelation step. Then I'll train the model. My questions: (1) Would my steps for scaling the data be ok, or are there any additional precautions I need to take while dealing with the sensor data? (2) When I'm using the trained model in real setting with a new subject, and I want it to make predictions while I collect the sensor information about the subject, how should I scale the data then? Should I scale it by the same factor with which I scaled the data used to train/test the model? (3) And what if I would like to update my model as I'm collecting new data from other subjects in an online manner? How should I scale the sensor data then? I will be trying different possibilities next and I was hoping to get some advice here. Thanks!


title: A little problem to achive this by coding in convex constraint

text: IN convex optimization : min f(x), I have two constraits : one is sum_i|x_i| <= t (t can be any numbers) , another constraints is |x_j|=|x_j-1 k_j|(set k_1 can not equal to zero and k_j can be thought as instrumental variable). As we know , if we do not have the second constraint , the models can be solved by LASSO using coordinate descent of LARS ; however if have the second constraint , how can we solve this problem ? using coordinate descent ?


title: How do I compare two groups by several dependent variables and get 1 p-value as a result?

text: Suppose I have two groups and several dependant variables. I would like to show statistically whether the two groups are different or not, based on the combination of these dependant variables. The dependant variables will be correlated to some extent. When doing t-test or ANOVA each of my dependant variables are almost significant based on my chosen p-value cut-off. I would like to harness the power of all of my dependant variables together and get a single p-value out from the comparison. I could do this by a repeated measures ANOVA but perhaps this wouldn't be appropriate since it would be more accurate to say that my measurements are of related but seperate dependant variables rather than an individual dependant variable with repeated measurements. It would probably be more appropriate to model the situation with a MANOVA, but with that test I don't think I can get a single value of significance to summarise over all dependant variables?


title: Why is a Bayesian not allowed to look at the residuals?

text: In the article "Discussion: Should Ecologists Become Bayesians?" Brian Dennis gives a surprisingly balanced and positive view of Bayesian statistics when his aim seems to be to warn people about it. However, in one paragraph, without any citations or justifications, he says: > Bayesians, you see, are not allowed to look at their residuals. It violates > the likelihood principle to judge an outcome by how extreme it is under a > model. To a Bayesian, there are no bad models, just bad beliefs. Why would a Bayesian not be allowed to look at the residuals? What would be the appropriate citation for this (i.e. who is he quoting)? Dennis, B.   Discussion: Should Ecologists Become Bayesians?   _Ecological Applications, Ecological Society of America_ , **1996** , 6, 1095-1103


title: Closed form recurrence formula for getting N consecutive heads on a coin

text: I want to find the expected number of coin tosses to get $N$ heads in a row, where $p$ is the probability of getting a head in a single toss. Let $F(N)$ be the expected number of tosses to get $N$ heads consecutively, so $$F(N) = 1 + p F(N-1) + (1-p) F(N)$$ which gives $$F(N) = 1/p + F(N-1)$$ With base condition : $F(1) = 1 + 1/p$ My logic is as follows: > if we get a head, in current toss we need to get N-1 more heads > consecutively, but if we get a tail, we have to start over This is what I thought, but is not correct. Can you please help me?


title: Hypothesis Testing with Big Data

text: How do you perform hypothesis tests with big data? I wrote the following MATLAB script to emphasize my confusion. All it does is generate two random series, and run a simple linear regression of one variable on the other. It performs this regression several times using different random values and reports averages. What tends to happen is as I increase the sample size, p-values on average get very small. I know that because the power of a test increases with sample size, given a large enough sample, p-values will become small enough, even with random data, to reject any hypothesis test. I asked around and some people said that with 'Big Data' its more important to look at the effect size, ie. whether the test is significant AND has an effect large enough for us to care about. This is because in large sample sizes the p-values will pick up of very small differences, like it is explained here. However, the effect size can determined by scaling of the data. Below I scale the explanatory variable to a small enough magnitude that given a large enough sample size, it has large significant effect on the dependent variable. So I'm wondering, how do we gain any insight from Big Data if these problems exist?               %make average     %decide from how many values to make average     obs_inside_average = 100;          %make average counter     average_count = 1;          for average_i = 1:obs_inside_average,                                   %do regression loop     %number of observations     n = 1000;          %first independent variable (constant term)     x(1:10,1) = 1;           %create dependent variable and the one regressor     for i = 1:10,              y(i,1) = 100 + 100*rand();              x(i,2) = 0.1*rand();          end                              %calculate coefficients     beta = (x'*x)\x'*y;          %calculate residuals     u = y - x*beta;          %calcuatate sum of squares residuals     s_2 = (n-2)\u'*u;          %calculate t-statistics     design = s_2*inv(x'*x);          %calculate standard errors     stn_err = [sqrt(design(1,1));sqrt(design(2,2))];          %calculate t-statistics     t_stat(1,1) = sqrt(design(1,1))\(beta(1,1) - 0);     t_stat(2,1) = sqrt(design(2,2))\(beta(2,1) - 0);          %calculate p-statistics     p_val(1,1) = 2*(1 - tcdf(abs(t_stat(1,1)), n-2));     p_val(2,1) = 2*(1 - tcdf(abs(t_stat(2,1)), n-2));                                   %save first beta to data column 1     data(average_i,1) = beta(1,1);          %save second beta to data column 2     data(average_i,2) = beta(2,1);          %save first s.e. to data column 3     data(average_i,3) = stn_err(1,1);          %save second s.e. to data column 4     data(average_i,4) = stn_err(2,1);          %save first t-stat to data column 5     data(average_i,5) = t_stat(1,1);          %save second t-stat to data column 6     data(average_i,6) = t_stat(2,1);          %save first p-val to data column 7     data(average_i,7) = p_val(1,1);          %save second p-val to data column 8     data(average_i,8) = p_val(2,1);          end          %calculate first and second beta average     b1_average = mean(data(:,1));     b2_average = mean(data(:,2));          beta = [b1_average;b2_average];          %calculate first and second s.e. average     se1_average = mean(data(:,3));     se2_average = mean(data(:,4));          stn_err = [se1_average;se2_average];          %calculate first and second t-stat average     t1_average = mean(data(:,5));     t2_average = mean(data(:,6));          t_stat = [t1_average;t2_average];          %calculate first and second p-val average     p1_average = mean(data(:,7));     p2_average = mean(data(:,8));          p_val = [p1_average;p2_average];          beta     stn_err     t_stat     p_val


title: Bootstrapping approach: how is the resampled sample chosen?

text: Bootstrapping approach is based on the assumption that the sample is the population. Then, a resampling of the sample is done, replacing the individuals of the population by chance. The number of the resampled individuals is the same as the 'real' sample. So, how is the new sample different from the original one? I don't have a huge background in statistic so please be plain.


title: Problems with bounded response variable

text: My response variable is a questionnaire score ranging **from 1 to 5** (which can take intermediate values like 1.3, 4.7, etc.). However, when the plot is made, I can see the density line **reaching values under 1 and over 5** , something like: ![enter image description here](http://i.stack.imgur.com/CtGr7.jpg) So my multiple question is: Is there a name to explain this phenomenon?? Why is this happenning and, the most important, is there any solution? I would appreciate any help.


title: What method should I use to analyze data?

text: I failed to pass my regression analysis course as the lecturer didn't say what kind of analysis he wants on the data. The problem was: > The index of biotic integrity (IBI) is a measure of the water quality in > streams. IBI and land-use measures for a collection of streams in the Ozark > Highland ecoregression of Arkansas were collected as part of a study. Table > below gives the data for IBI and the area of watershed is square kilometers > for streams in the original sample with area less than or equal to 70 > $\text{km}^2$. > > Use numerical and graphical methods to describe the variable IBI. Do the > same for area. Summarize your results.               Area IBI       21 47       34 76        6 33       47 78       10 62       49 78       23 33       32 64       12 83       16 67       29 61       49 85       28 46        8 53       57 55        9 71       31 59       10 41       21 82       26 56       31 39       52 89       21 32        8 43       18 29        5 55       18 81       26 82       27 82       26 85       32 59        2 74       59 80       58 88       19 29       14 58       16 71        9 60       23 86       28 91       34 72       70 89       69 80       54 84       39 54        9 71       21 75       54 84       26 79      I have no idea what method should I use and what kind of analysis would be useful. Is that kind of questions normal in basic studies? Should I try to find some book on regression analysis which contains more rigorous exercises?


title: $\bar{X}$ versus $\mathbb{E}(\bar{X})$?

text: I was not able to find this question here, so I am going to ask this:   * What is the difference between $\mathbb{E}(\bar{X})$ (expected value of $X$ bar) and the actual $\bar{X}$? I am very confused about these two concepts. How come $\bar{X}$ is one of the estimators of normal distribution (the other one being $S^2$) and then what is the point of $\mathbb{E}(\bar{X})$?   * Another question that I have is, what is the relationship between $\bar{X}$ and $E(\widehat{\mu})$ ? I understand that mu is a true mean, but then what is $E(\widehat{\mu})$? I would appreciate any explanation on these concepts and since I'm only a beginner in math stats, I am struggling with notation, so any simplified explanation will be much appreciated! And then, I just discovered that there is also an expected value of $\sigma^2$ which, why lie, completely blew my mind! So I suspect I'm struggling with the whole concept of expected value and how it relates to the the population, sample distribution, and sampling distribution, so I would definitely appreciate an explanation of expected value of $\bar{X}$ in light of this.


title: Calculation AIC without loglikelihood-function

text: I want to calculate the AIC without calculating the loglikelihood-function (which seems complicated). If the residuals are normally distributed, this can be done, according to wikipedia, as follows: $AIC_{\sigma} = n log(\sigma^2_Z) + 2k $, where $\sigma^2_Z$ is the variance of the residuals, $n$ the number of samples and $k$ the number of parameters. In order to test this I wrote a simple R-script:               library(MASS)     data(cats)     l<-lm(Hwt ~ Bwt,data=cats)          l.resi <- resid(l)     qqnorm(l.resi)     qqline(l.resi)     shapiro.test(l.resi)          n<-nrow(cats)     k<-2 # one cefficient + intercept     aic_s<- n*log(mean(l.resi^2) ) + 2*k           print(paste("aic: ",AIC(l)))      print(paste("aic_s: ",aic_self))      The script returns:               Shapiro-Wilk normality test          data:  l.resi     W = 0.9845, p-value = 0.1046          [1] "aic:  520.121593704369"     [1] "aic_s:  109.467296141423"      So, the residuals seem to be normally distributed. However, the AICs differ. Why is that?


title: Existence of the moment generating function and variance

text: Can a distribution with finite mean and infinite variance have a moment generating function? What about a distribution with finite mean and finite variance but infinite higher moments?


title: How to prove three properties of the moment generating function?

text: The moment generating function of a random variable $X$ is defined to be the function $$M_{X}(t)=E(e^{tX})=\sum_{n=0}^{\infty}\frac{E(X^n)}{n!}t^n.$$ Let $I=\\{t\in\mathbb R:M_{X}(t)<\infty\\}.$ I wish to show that   1. $I$ is possibly a degenerate interval and $0\in I$. (Degenerate means the interval includes only one real number.)    2. $M_{X}(t)$ is a convex function on $I$.    3. If $0$ is an interior point of $I$, then $E(X^k)\lt\infty$ for all $k\in \mathbb N$; i.e., $X$ has finite moments of all orders.


title: I got a dataset, there are about 400 variables,most of them are quite similar

text: I got a dataset, there are about 400 variables, most of the variables are quite similar. Actually, I would like to say their number are pretty close, but I think they won't too correlated, since the graph for the variables is like a heart beat monitor. But just at some point it varies too much. ![pic for different variables](http://i.stack.imgur.com/HgtbA.png) At first, I was trying to perform a PCA to it. Here is the out put.But it doesn't turn out very well. ![PCA](http://i.stack.imgur.com/gwKQi.png) Maybe it is because the similar variables cover too much of the dataset. I wonder if there are any other way to analyze it? to reduce dimension?


title: Using test data sets in simulation study

text: I'd like to know the correct way to simulate test data for a simulation study. For simplicity, suppose that I want to test a linear regression model. Chapter 7 of ESL explains that the average test error is equal to the expected prediction error, where the average is over the training data. $\text{Err}_{\mathcal{T}}$ is the test error for a specific training set $\mathcal{T}$. They have an example stating that $\text{Err}_{\mathcal{T}}$ was calculated for 100 simulated training sets. $\sum_{\mathcal{T}=1}^{100}\text{Err}_{\mathcal{T}}/100$ is then an estimate of the expected prediction error $\text{Err} = E(\text{Err}_{\mathcal{T}})$. My question is: must 100 test sets be simulated for the calculation or is each $\text{Err}_{\mathcal{T}}$ calculated using the same test set? Here is some R code to demonstrate:               library(MASS)     n = 50     m = 200     b0 = 0.5     b = c(1,2,0,0)     rho = 0.5     sigma = 1     iter = 100     p = length(b)     r = matrix(rho, p, p); diag(r) = 1          # test set option A     # x.new = mvrnorm(m, rep(0,p), r)     # y.new = b0 + x.new %*% b + rnorm(m, 0, sigma)          err.t = rep(0, iter)     for (i in 1:iter) {       x = mvrnorm(n, rep(0,p), r)       y = b0 + x %*% b + rnorm(n, 0, sigma)            # test set option B       x.new = mvrnorm(m, rep(0,p), r)       y.new = b0 + x.new %*% b + rnorm(m, 0, sigma)            mod = lm(y ~ x)       pred = predict(mod, x.new, type="response")       err.t[i] = sum((y.new - pred)^2)/m     }     err = mean(err.t)      Should the test set be generated outside the loop as in option A, or inside the loop as in option B? The same question goes for a validation set. If I want to use a separate data set for model selection, should the validation set be outside or inside the loop?


title: Tokenizer for text categorization

text: In text categorization, it is common to use tokens as features. But, there are several different ways to convert a sentence to tokens. For example, see the different types of Tokenizers here: https://github.com/NaturalNode/natural . What is the most commonly used tokenizer for categorizing short texts?


title: Unbiased estimator of variance for samples *without* replacement

text: This is a follow-up question on that one:   Could Bessel's correction make sample variance estimation even more biased? I understand that you need Bessel's correction to get an unbiased estimate of variance in case you sample _with_ replacement. Now in practice most sampling is done _without_ replacement (e.g. asking voters about their preferences or studying the effects of new drugs). My questions are:   \- What is the unbiased estimator for variance in cases of samples _without_ replacement?   \- Choosing only between the classical estimators for _population variance_ and _sample variance_ wouldn't it be better in most cases to take the _population variance_ for samples because it is less biased than _sample variance_ in cases of samples _without_ replacement?


title: Classification between most probable classes - what is its proper name?

text: I have to classify data between N classes. At the first step, I make a classification between N classes with low threshold (for example, 10%). When first step is over, I have $M$ labels assigned to my datum $(M \leqslant N)$. At the second step, I make normal classification with normal threshold between the $M$ classes to get the final decision. What is the proper name of such policy?


title: Non-random distribution of SNPs (single-nucleotide polymorphisms)

text: I have a bunch of groups of DNA sequences in which I have observed changes in different regions (block D and H) vs. the rest of the molecule. The data is available as a table like the one below. The first row means the total number of changes in cluster 1 were 4, 0 in block D and 0 in block H, meaning they were found in the rest of the molecule.               cluster TotalSNPS blockD blockH     c1          4        0      0     c2          4        2      0     c3          3        1      2     c4          4        2      0      I would like to test if the distribution of changes within the molecule is random. My expectation is that SNP changes should occur just about anywhere. What statistical test should I apply or what is a suitable approach to test this?


title: When is binomial distribution function above/below its limiting Poisson distribution function?

text: Let $B(n,p,r)$ denote the binomial distribution function (DF) with parameters $n \in \mathbb N$ and $p \in (0,1)$ evaluated at $r \in \\{0,1,\ldots,n\\}$: \begin{equation} B(n,p,r) = \sum_{i=0}^r \binom{n}{i} p^i (1-p)^{n-i}, \end{equation} and let $F(\nu,r)$ denote the Poisson DF with parameter $a \in \mathbb R^+$ evaluated at $r \in \\{0,1,2,\ldots\\}$: \begin{equation} F(a,r) = e^{-a} \sum_{i=0}^r \frac{a^i}{i!}. \end{equation} Consider $p \rightarrow 0$, and let $n$ be defined as $\lceil a/p-d \rceil$, where $d$ is a constant of the order of $1$. Since $np \rightarrow a$, the function $B(n,p,r)$ converges to $F(a,r)$ for all $r$, as is well known. With the above definition for $n$, I'm interested in determining the values of $a$ for which \begin{equation} B(n,p,r) > F(a,r) \quad \forall p \in (0,1), \end{equation} and similarly those for which \begin{equation} B(n,p,r) < F(a,r) \quad \forall p \in (0,1). \end{equation} I have been able to prove that the first inequality holds for $a$ sufficiently smaller than $r$; more specifically, for $a$ lower than a certain bound $g(r)$, with $g(r)<r$. Similarly, the second inequality holds for $a$ sufficiently larger than $r$, i.e. for $a$ greater than a certain bound $h(r)$, with $h(r)>r$. (The expressions of the bounds $g(r)$ and $h(r)$ are irrelevant here. I will provide the details to anyone interested.) However, numerical results suggest that those inequalities hold for less stringent bounds, that is, for $a$ closer to $r$ than I can prove. So, I'd like to know if there is some theorem or result which establishes under which conditions each inequality holds (for all $p$); that is, when the binomial DF is guaranteed to be above/below its limiting Poisson DF. If such theorem doesn't exist, any idea or pointer in the right direction would be appreciated. Please note that a similar question, phrased in terms of incomplete beta and gamma functions, was posted in math.stackexchange.com but got no answer.


title: How to predict or extend regression lines in ggplot2?

text: I have a data frame that contains two time series: the dates and version numbers of Emacs and Firefox releases. Using one ggplot2 command it's easy to make a chart that uses loess (in a way that looks a bit amusing, which I don't mind) to turn the points into lines. How can I extend the lines into the future? I want to determine where and when Emacs and Firefox version numbers will cross, and if there's a way to show an error range, all the better. Given that ggplot2 is plotting the lines, it must have a model, but I don't see how to tell it to extend the lines, or to get the model out and do something with it.               > library(ggplot2)     > programs <- read.csv("http://www.miskatonic.org/files/se-program-versions.csv")     > programs$Date <- as.Date(programs$Date, format="%B %d, %Y")     > head(programs)       Program Version       Date     1   Emacs    24.1 2012-06-10     2   Emacs    23.4 2012-01-29     3   Emacs    23.3 2011-03-10     4   Emacs    23.2 2010-05-08     5   Emacs    23.1 2009-07-29     6   Emacs    22.3 2008-09-05     > head(subset(programs, Program == "Firefox"))        Program Version       Date     18 Firefox      16 2012-10-09     19 Firefox      15 2012-08-28     20 Firefox      14 2012-06-26     21 Firefox      13 2012-06-15     22 Firefox      12 2012-04-24     23 Firefox      11 2012-03-13     > ggplot(programs, aes(y = Version, x = Date, colour = Program)) + geom_point() + geom_smooth(span = 0.5, fill = NA)      ![Emacs and Firefox versions by date](http://i.stack.imgur.com/Mswg9.png) (Note: I had to fudge the early Firefox versions and turn 0.1 onto 0.01, etc., because "dot one" and "dot ten" are equal arithmetically. I know Firefox is releasing every six weeks now, but they don't exist yet, and I'm interested in a general answer to this prediction question.)


title: Converting a list of partial rankings into a global ranking

text: I'm working on something like the following problem. I have a bunch of users and N books. Each user creates an ordered ranking of all the books he's read (which is likely a subset of the N books), e.g., Book 1 > Book 40 > Book 25. Now I want to turn these individual user rankings into a single ordered ranking of all the books. Are there any good or standard approaches to try? So far, I'm thinking of Bradley-Terry models applied to pairwise comparisons, but I'm wondering if there's anything else.


title: Calculate goodness-of-fit (with deviance) to compare averaged models?

text: I need to compare the goodness of fit of several averaged logistic regression models by calculating the deviance explained. I'm using the `MuMIn` package in R to average many logistic regression models into a single averaged model. I ultimately want to compare the explanatory power of several averaged models, in part by using the deviance explained. My questions are:   1. Does deviance explained apply to averaged models as a strong measure of the goodness of fit?   2. How does one calculate the deviance explained (calculated as the null deviance less the residual deviance as a proportion of the null deviance) from the averaged model output from the `model.avg()` command in `MuMIn`?  Examining the structure of the averaged model object indicates that the null and residual deviances are calculated on each individual model that contributes to the averaged model, but I'm not sure how to extract them and then calculate the overall deviance explained by the averaged model.


title: Linear regression with a R squared of 0.5

text: I am trying to fit a linear model to my data, which consist of a continuous dependent variable and several predictors, some continuous, some discrete and one dummy. I used the R package glmnet, and the best R squared that I get is of 0.5. Is there something that I could try to improve the model? Perhaps, I could scale the variables, but since I have a dummy variable, it seems that I should not.


title: Holt Winters initialization using backcasting like SPSS

text: I'm trying to figure out how SPSS is initializing the multiplicative Holt Winters exponential smoothing model using backcasting. Thankfully IBM roughly described their way of doing so .. very roughly. Assuming a seasonal period of 3 and 3 seasons (9 observations) for the sake of simplicity. As far as I understand they fit lines for all corresponding observations of all seasons with intercept. Fitting observations 1,4,7 with intercept $µ_1$ and slope $β_1$ .. fitting observations 2,5,8 with intercept $µ_2$ and slope $β_2$ .. and so on. So the initial trend is the average of all slopes. But how exactly are the seasonal indices computed? And isn't the intercept mostly the first observation on a regression line, is it?? (e.g. observation 1 for the first line) EDIT: I guess it will help me a lot if I know more about how backcasting is performed in general. All I found in literature is that the time series is reversed and smoothed values are generated for the season before $t_1$. But how do I get the initial values used for forecasting when I only have compound backcasts?


title: Conversion to Standard Normal with Probability

text: Good day everyone, I am currently working on a self-study question: **Question** > The number of viewers of a television show has a mean of 29 million with a > standard deviation of 5 million. Assume this distribution follows a normal > distribution. What is the probability that at least one of the four shows > has more than 36 million viewers in the next four shows. Assume the number > of viewers in each show are independent. **My Attempt** > ![enter image description here](http://i.stack.imgur.com/McTbw.jpg) However, it seemed that the answer is not correct from that of the answer key (0.2861). Appreciate some guidance on this. I know I'm missing out on something, but I just can't figure out what.


title: How to conduct Johanson's cointegration test with variables stationary at I(0)?

text: How do you interpret the Johanson's cointegration test if you are doing Granger causality and both variables are stationary at I(0)?


title: About robust regression

text: I'm working on a robust regression analysis with `SAS`. I'm running a full model first.(using `PROC ROBUSTREG`). Then I'm excluding the variable with maximum Pr>ChiSq value. And then I'm running the model again. I'm repeating the process until Pr>ChiSq value is under 0.05 for all variables in the model. This is similar with backward elimination selection.   Is this approach correct?


title: Is computer network throughput an independent variable?

text: The central limit theorem assumes independent variables and I want to apply it to study the performances of my network speed (computer networks). What I have a is the throughput measurements of my network per second. I don't know how to classify my variable, do you think it's a independent variable? I have a doubt because this measurements can be influenced by other factors like how many people are using the network and with what type of application etc.


title: Converting PDF to text format

text: I am searching for a tool/software which can convert PDF to text format for text mining. The PDF files I use, contains many tables. I tried using PDF Miner, but I found many problems while converting tables in PDF to text. Can someone recommend a tool that efficiently works? Thanks


title: Online algorithm to compute variance with a decay

text: Could somebody point me to an online algorithm that computes the variance, but gives a higher weight to more recent values?


title: Basic Probability

text: How to find the conditional/joint probabilities of two variables from their empirically estimated marginal ones and an assumed correlation?


title: Fuzzy C-means and its stages of clustering

text: when dealing with soft and hard clustering techniques such as K-means and fuzzy C-means I run into abit of difficulty on the steps that FCM takes to calculate the clusters. For instance in K-means the steps are as follows:   1. Chose number of clusters (K)   2. Initialize centroids (K patterns randomly chosen from data set)   3. Assign each pattern to the cluster with closest centroid   4. Calculate means of each cluster to be its new centroid   5. Repeat step 3 until a stopping criteria is met (no pattern move to another cluster) Could some one explain the steps of fuzzy C-means?


title: Is it a reasonable rule of thumb to examine Pearson and Spearman correlations, and use the latter if they are very different?

text: The data I am using is continuous, and may be somewhat skewed, but probably will not have many outliers. It's been suggested to me that it would be a good rule of thumb to examine both Pearson and Spearman correlations, and if they are different to report the latter, on the grounds that a difference between the two might indicate that the assumptions of Pearson's _r_ are not met.


title: What is a normal distribution with 'common variance'?

text: _Data Augmentation Approach in Bayesian Modelling of Presence-only Data_ (Divino _et al._ ) states that the data ware sampled from a normal distribution with mean 2.0 and a common variance (section 3). Does anyone know the value of a 'common variance'? Thanks.


title: How do I calculate the distribution of number of events in the busiest period?

text: I've got an estimate of the number of site visitors I'll see in a 1 hour period clicking email links in a large email campaign. I need to make sure I've got the required server capacity. That means I need to estimate the peaks in usage on a smaller time scale.   1. How do I calculate the expected distribution for the number of visitors I should expect in any given smaller time period? eg for a given 1 minute period.   2. How do I calculate the expected distribution for number of visitors in the **busiest** time period of a given length? I figure it's good enough for estimation purposes to assume that each visitor arrives at a random time that's evenly distributed within the one hour period. \-- **EDIT:** (in response to rocinante) OK, so lets consider the problem redefined in terms of the assumptions of the poisson distribution. There's some minor differences, but they don't matter to me. Also, lets not worry about possible notions like sliding windows of time for the period in 2. So long as the period in 2. is short, that doesn't matter either. Supposing I expect visitors to arrive at a rate of 10K visitors per hour over a one hour period. If I divide that hour up into 720 adjacent 5 second intervals, I expect the number of visitors in each of those 5 second periods to follow a poisson distribution. Part 2 of my question can then be restated: I can take 720 samples from a poisson distribution, and find the largest value of those 720 sample values. If I repeatedly took sets of samples like that, calculating the maximum value in each set, what would the distribution of those maximums be like? Having graphed some poisson distributions, I'm comfortable that I know enough about what to expect for planning purposes, but I'd still like to how to do part 2 of my problem as posed. Also, for large lambda, I presume some other calculation of the distribution should be used? Calculating factorials becomes impractical at some point.


title: Which test should be run on this confusion table data?

text: I have several confusion tables calculated, see, e.g., table 2 and table 1, but typically I'll have 10 or so. These confusion tables are the result of different classification processes (in fact they are the result of different sensor locations). The row and column labels are the same in each table, only the values and thus the accuracy (the last column) differs. Now I want to run a test which can tell me which table performs better (i.e, which sensor location is best). I thought of just taking the average of the accuracy, but I would also like to know which rows attribute the most to the difference.


title: Active learning through attribute & instance sampling?

text: I work in a classification problem where **collecting** instances (labeled or unlabeled) is very expensive. In fact, labeling instances is cheap. I want to explore an **active learning strategy** where I do **not** have access to unlabeled instances, but I can still ask for **labeled instances** by specifying:   * Statistics on the attributes on the new instances (e.g. give me an instance with an age > 45)   * The instance label (e.g. give me a new instance of label $l$, where $l \in L$ is one of the classification labels) or the label sampling distribution itself. My questions are:   * **Does this problem have a name?** It looks like a specific case of Active Learning, but I am not sure, since **in Active Learning one starts with a set of unlabeled instances, which is not my case**.   * What types of approaches (from the most rudimentary to the more sophisticated) can I employ to identify the most informative sampling distribution from instance attributes or instance labels?


title: How to interpret marginal R squared?

text: i have obtained a marginal R squared value of 0.01599143 after fitting a tweedie model to my data set. Is this value too low? How can i interpret this?Thanks


title: Diagnostic plots for lmer

text: I am trying to produce a glmm using the `lme4` package in `R`. To validate my model I would like to produce some diagnostic plots but am having trouble doing so as many of the packages that work for `lm()` models like CAR do not work for `lmer()` models. I have tried using the `languageR` and `asuR` packages but am just getting errors. I am particularly interested in creating partial residual plots. I have been able to produce qqplots, residual vs fitted plots and look at cooks distances manually. Any suggestions of further plots I should try or packages that would help would be greatly appreciated.               fm <- lmer(Average.payoff ~ Game + Type + Others.Type + Game*Type + Game*Others.Type +       Type*Others.Type + (1|Subjects), REML=FALSE, data=Subjectsm1)          library(asuR)     inspect(fm)     Error in .setupMethodsTables(fdef, initialize = TRUE) :        no slot of name "group" for this object of class "derivedDefaultMethod"          library(languageR)     plotLMER.fnc(fm)     warning: I don't know how to handle  Game      Error in preparePredictor.fnc(predictors[i], model, m, ylabel, fun, val = NA,  :        object 'dfr' not found


title: Why use logged variables?

text: Probably, this is a very basic question but I don't seem to be able to find a solid answer for it. I hope here, I can. I'm currently reading papers as a preparation for my own master's thesis. Currently, I'm reading a paper which researches the relationship between tweets and stock market features. In one of their hypothesis, they propose that "increased tweet volume is associated with an increase in trading volume". I would expect them, in the pairwise correlations, to correlate `tweetVolume` with `tradingVolume`, but instead they report using the logged versions: `LN(tweetVolume)` and `LN(tradingVolume)`. For my thesis, I have replicated this bit of their paper. I have collected tweets about 100 companies for over 6 months (`tweetVolume`) and stock trading volume for the same timeframe. If I correlate the absolute variables, I find `r=.282, p.000` but when I use the logged verions, I find `r=.488, p=.000`. I don't understand **why** researchers sometimes use logged versions of their variables and why correlation seems so much higher if you do so. What is the reasoning here, and why is it OK to use logged variables? Your help is greatly appreciated :-)


title: Conditions for conclusion

text: A study had shown that annual income on average was significantly higher among men than women for people over 30 years old in Portugal. The results were from 312 Portuguese men and 62 Portuguese women. The samples were chosen (not randomly) by the researchers. Can we conclude that, in general, the annual income for Portuguese men is higher than the average income for Portuguese women? I think that we can _not_ conclude this because of the following two reasons:   1. The samples for both men and women were too low.   2. The samples were _not_ chosen randomly. Am I correct? Why or why not?


title: How to use ks test for 2 vectors of scores in python?

text: I am trying to figure out how to determine whether to reject or not the null hypothesis using ks-test.   1. In matlab there is a function named kstest2 decides if it should reject $H_0$. I want to use python scipy function ks_2samp but it returns the p-value and the $\alpha$, how can I determine if one should reject $H_0$ given those 2 parameters ?   2. If I have 2 vectors of scores per day, how can I use Matlab/python kstest2 to see if the are distributed the same?   3. Is there a clean way to change the scores vector to a continuous distributed vector?


title: Exponential of a standard normal random variable

text: We know that $Z\sim N(0, 1)$. How do I prove that $e^Z$ has a mean of $e^{0.5}$? I have tried integrating $e^z$ times the pdf of $Z$ but I don't know why it isn't working out. Also what is the pdf of $\exp(Z)$?


title: What does the z-score mean in plain English (in A/B split testing)?

text: I'm relatively new to A/B split-testing and can't wrap my brain around the idea of z-scores. I know that a z-score gives one an idea how statistically significant the result is that one gets from a split-test. However, I found quite a few websites that explain A/B split-testing **without talking much about sample sizes**. These articles even suggested that the z-score itself would only show statistical significance if the datasets used to calculate it is large enough. So I ran a one-tailed test. I made a modification to an element of a page and checked how the changes performed in terms of conversion rates. This was the result:               | Treatment   | Hits | Conv | %      | z-Score |     |-------------|:----:|:----:|:------:|--------:|     | Control     | 701  | 200  | 28.53% | -       |     | Treatment 1 | 699  | 228  | 32.62% | 1.66    |      Confidence (z-Score): **95.17%** , improvement of conversion rate: **14.3%** However, I figured later that this result lacks two things:   1. I didn't define a desired minimum improvement in percentage in my hypothesis, e. g. "Treatment 1 performs at least 20% better than Control"   2. I didn't define a required sample size but just let the test run until a z-score >= 1.65 was reached (which I heard is a bad thing to do). I then decided I want at least a **22% improvement** of my conversion rate and that I needed **2,046 datasets** to tell whether or not this goal was achieved (calculated based on this blog post). The hypothesis was then: > Treatment 1 improves the conversion rate by at least 22% and I can tell with > a probability of 95% that this change is not due to chance. The final result (2,046 datasets) showed that Treatment 1 didn't improve anything really and that I can reject the hypothesis:               | Treatment   | Hits | Conv | %      | z-Score |     |-------------|:----:|:----:|:------:|--------:|     | Control     | 1020 | 305  | 29,90% | -       |     | Treatment 1 | 1026 | 308  | 30,02% | 0,06    |      **What I don't understand is:** What did the z-score of the 1,400 datasets tell me in plain English? What is it good for in this case?


title: Bootstrapping robust regression

text: I am trying to get some bootstrap CIs for coefficients obtained by robust regression. I have influential values and thus switched to `rlm`. The data are clustered and within clusters, the variance of my DV ("Durchlässigkeit") = 0. Is this all sensible for my clustered data?               > dput(d)     structure(list(Porenfläche = c(4990L, 7002L, 7558L, 7352L, 7943L,     7979L, 9333L, 8209L, 8393L, 6425L, 9364L, 8624L, 10651L, 8868L,     9417L, 8874L, 10962L, 10743L, 11878L, 9867L, 7838L, 11876L, 12212L,     8233L, 6360L, 4193L, 7416L, 5246L, 6509L, 4895L, 6775L, 7894L,     5980L, 5318L, 7392L, 7894L, 3469L, 1468L, 3524L, 5267L, 5048L,     1016L, 5605L, 8793L, 3475L, 1651L, 5514L, 9718L), P.Perimeter = c(2791.9,     3892.6, 3930.66, 3869.32, 3948.54, 4010.15, 4345.75, 4344.75,     3682.04, 3098.65, 4480.05, 3986.24, 4036.54, 3518.04, 3999.37,     3629.07, 4608.66, 4787.62, 4864.22, 4479.41, 3428.74, 4353.14,     4697.65, 3518.44, 1977.39, 1379.35, 1916.24, 1585.42, 1851.21,     1239.66, 1728.14, 1461.06, 1426.76, 990.388, 1350.76, 1461.06,     1376.7, 476.322, 1189.46, 1644.96, 941.543, 308.642, 1145.69,     2280.49, 1174.11, 597.808, 1455.88, 1485.58), P.Form = c(0.0903296,     0.148622, 0.183312, 0.117063, 0.122417, 0.167045, 0.189651, 0.164127,     0.203654, 0.162394, 0.150944, 0.148141, 0.228595, 0.231623, 0.172567,     0.153481, 0.204314, 0.262727, 0.200071, 0.14481, 0.113852, 0.291029,     0.240077, 0.161865, 0.280887, 0.179455, 0.191802, 0.133083, 0.225214,     0.341273, 0.311646, 0.276016, 0.197653, 0.326635, 0.154192, 0.276016,     0.176969, 0.438712, 0.163586, 0.253832, 0.328641, 0.230081, 0.464125,     0.420477, 0.200744, 0.262651, 0.182453, 0.200447), Durchlässigkeit = c(6.3,     6.3, 6.3, 6.3, 17.1, 17.1, 17.1, 17.1, 119, 119, 119, 119, 82.4,     82.4, 82.4, 82.4, 58.6, 58.6, 58.6, 58.6, 142, 142, 142, 142,     740, 740, 740, 740, 890, 890, 890, 890, 950, 950, 950, 950, 100,     100, 100, 100, 1300, 1300, 1300, 1300, 580, 580, 580, 580), Gebiete = structure(c(1L,     1L, 1L, 1L, 2L, 2L, 2L, 2L, 6L, 6L, 6L, 6L, 4L, 4L, 4L, 4L, 3L,     3L, 3L, 3L, 7L, 7L, 7L, 7L, 9L, 9L, 9L, 9L, 10L, 10L, 10L, 10L,     11L, 11L, 11L, 11L, 5L, 5L, 5L, 5L, 12L, 12L, 12L, 12L, 8L, 8L,     8L, 8L), .Label = c("6.3", "17.1", "58.6", "82.4", "100", "119",     "142", "580", "740", "890", "950", "1300"), class = "factor")), .Names = c("Porenfläche",     "P.Perimeter", "P.Form", "Durchlässigkeit", "Gebiete"), row.names = c("1",     "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13",     "14", "15", "16", "17", "18", "19", "20", "21", "22", "23", "24",     "25", "26", "27", "28", "29", "30", "31", "32", "33", "34", "35",     "36", "37", "38", "39", "40", "41", "42", "43", "44", "45", "46",     "47", "48"), class = "data.frame")          ## do robust regression and bootstrap the coefficients, allowing for clustered data     ## by putting "Gebiet" as strata argument (?),     ## dv variation within clusters/Gebiet = 0!     bs <- function(formula, data, indices) {       d <- data[indices, ] # allows boot to select sample       fit <- rlm(formula, data = d)       return(coef(fit))     }          results <- boot(data = d, statistic = bs, strata = d$Gebiete,                     R = 199, formula = Durchlässigkeit ~ P.Perimeter + P.Form)          # get 99% confidence intervals     boot.ci(results, type="bca", index=1, conf = .99) # intercept     boot.ci(results, type="bca", index=2, conf = .99) # P.Perimeter     boot.ci(results, type="bca", index=3, conf = .99) # P.Form


title: Fisher Distance for feature selection

text: I'm currently working for EEG signal classification from 3 electrodes. I want to have a simple feature selection algorithm that is independent with the classification process. From the feature extraction step, let's say I have this kind of matrix now (not the actual numbers/data) : A CLASS : $$ Ch_1 = \begin{vmatrix} 1 & 2 & 3 \\\ 0.5 & 0.2 & 0 \\\ 1 & 0.1 & 0.8 \\\ 1.2 & 0.8 & 1 \end{vmatrix} Ch_2 = \begin{vmatrix} 1 & 1.5 & 1 \\\ 0.3 & 0.1 & 2 \\\ 1.3 & 0.1 & 3 \\\ 1.5 & 1.8 & 2 \end{vmatrix} Ch_3 = \begin{vmatrix} 2 & 2 & 3 \\\ 1.2 & 2 & 0.8 \\\ 1.3 & 1.2 & 1.5 \\\ 1.8 & 3 & 2 \end{vmatrix} $$ B CLASS : $$ Ch_1 = \begin{vmatrix} 1 & 2 & 3 \\\ 0.5 & 0.2 & 0 \\\ 0.1 & 2 & 0 \\\ 1.2 & 0.8 & 1 \end{vmatrix} Ch_2 = \begin{vmatrix} 1.2 & 1.5 & 1 \\\ 0.3 & 0.1 & 2 \\\ 0.8 & 1.1 & 0 \\\ 1.5 & 1.8 & 2 \end{vmatrix} Ch_3 = \begin{vmatrix} 2 & 2 & 3 \\\ 1.2 & 2 & 0.8 \\\ 0.2 & 1 & 0.3 \\\ 1.8 & 3 & 1 \end{vmatrix} $$ Where on the example above, the row of the channels are the numbers of trials/observations (4 trials per class) and the column are the features extracted from each sub-band (3 features).   What I want to do is selecting which feature will give me better separation of data between classes, while maintaining close relationship within its own class.   I am trying to do Fisher Distance approach : $$ FisherDis = S_B/S_w $$ Where $S_B$ is between class matrix and $S_w$ is within class matrix. From what I read, I have to score each feature and then select some features with highest scores. Now to my question:   1\. What is "the number of samples" when I want to calculate $S_w$ and $S_B$ , is it four (as in four trials) or three (as in three features) ?   2\. Should I group the channels into one matrix? Or is it better if I'm working in each channel separately?   3\. Am I working on the right path? I have doubts in myself... Thank you very much in advance. I'd appreciate every answer from everyone because I'm fairly new to statistics (I have so much to learn..) :)


title: Calculating similarity and clustering question

text: I have a dataset of about a million companies containing their names, total employees and annual sales. I want to come up with a function that when given the company returns the 5 most similar companies in terms of their distance in total employees and annual sales. I thought of doing k-means clustering on the dataset and find clusters. Then return all the companies in that cluster. The problem with this approach is that I don't know the number of clusters I should form beforehand. Also, on a separate note, if I were to obtain a list of specialties (e.g. marketing, software, etc.) for each of the companies - how can I transform this qualitative value in to a number which can later help me calculate similarity.


title: How do I compute sales forecasts?

text: What's the statistical strategy to compute 2 month's ahead sales forecasts? I'm trying to build a chart for some software I'm creating. Piecemealing this together with a lot of thought, it appears that I'm faced with:   * **Seasonal Delta (SD)** Sales often tend to be seasonal, no matter the business. If sales went up or down in Q4 last year, they often, but not always, will mimic this pattern in Q4 this year.   * **Large Sample Delta (LSD)** There's already often, but not always, a natural rate of growth or decline occurring from now and 12 months prior to now, and one can carry that out too. That then becomes a factor.   * **Recent Sample Delta (RSD)** There's already often, but not always, a natural rate of growth or decline occurring from now and the prior month, and one can carry that out too. That then becomes a factor. Therefore, I guess if I average the top 3 delta's together, I get a **Best Guess Delta (BSD)** that I can apply to this month's average to get next month's average. Do it one more time and get the month following that. So, that's a 2 month forecast. Is that the statistical strategy to do this? Or, should I weight things? For instance, wouldn't RSD have more weight in probability than LSD and SD? And wouldn't LSD have the least weight in probability? _Please note -- my math training stopped right before Calculus. So, Algebra, Trig, Geometry -- these are things I know. I'm also a PHP programmer, if that means anything._


title: Newey-West standard errors with cross-sectional OLS?

text: Consider the cross sectional: $Y_i = a + b X_i + e_i$ where I have reason to believe that $E[e_j e_k] \not= 0$ for a concerning number of $j\not= k$. What happens if I use a serial correlation robust standard error here (such as Newey West)? Is it okay to do so even though, of course, I have no serial correlation in my data, yet I have dependence? Will Newey West correct for the fact that I have non-zero non-diagonal elements of my variance-covariance matrix?


title: Comparing sample covariance matrices to a model

text: I have several sample variance/covariance matrices, say $[C_1...C_n]$. These covariance matrices correspond to blocks of correlated multivariate normal samples (i.e. the samples used to estimate $C_1$ are not completely independent, but they are completely independent from the samples used to estimate $C_2$). In reality, there is a true variance/covariance matrix $C$; the $C_i$ should be samples from some distribution with expectation $C$. I also have two nested models that lead to predictions of this true variance/covariance matrix (say, $V_0$ and $V_1$). I'm looking to test formally whether $V_1$ is a better fit to $[C_1...C_n]$ than $V_0$ (the model for $V_0$ is nested within $V_1$). I've looked for a standard way to do this, but can't seem to find anything (e.g. Wilks' lambda and other MANOVA-type tests don't seem to be exactly what I'm looking for). Is there a way to write down a likelihood for $[C_1...C_n]$ given $V_0$ (or $V_1$)? Or some other way to compare the fit of these models?


title: How to find weekly sum from daily data in SAS?

text: I found it is quite inconvenient to do row operations in SAS. Any good ideas?


title: What is a 'bagplot', or 'bivariate boxplot'?

text: I've found a paper which introduces the multidimensional (bivariate here) version of the boxplot - a bagplot. What is that bagplot exactly? I can see the series of nested polygons based on vertices, one of those polygons being declared as a bagplot. What is the idea of nested polygon building? Which of the polygons is the bagplot (central or holding the average number of points)? Do the edges of a bagplot possess some useful properties (like specifically dividing the point set)?


title: How to use Kullback–Leibler divergence to help me choose the best distribution function?

text: I am now using **Exponential distribution** to model the the **time intervals of a sequence of random events**.Since I can choose several different lamdas for this model , I want to find out which lamda of them is fitting the training data best. I find Kullback–Leibler divergence from some paper is used to do this , but :   * I don't know if this way is suitable    * I don't know how to operate either I am new to machine learning . Anyone can help me ?


title: dealing with missing values via pro-rating

text: I'm aware there are a variety of approaches of dealing with missing values. I have been asked to analyse data that contains missing values. An author of a scale has recommended pro-rating when there are less than 3 missing values in a case. I am interested in advice regarding how this might be done. For example, if there are 25 items: 3 missing, 5 values of 2 and 8 values of 1. Dividing 18 by 22 yields a value of 0.8. Does this mean each missing value would be replaced by a value of 1? Any advice is appreciated, Bob


title: Joint entropy of two random variables

text: The joint entropy is the amount of information we get when we observe X and Y at the same time, but what would happen if we don't observe them at the same time. For example, when i toss a coin, if i got tails i will only observe the variable X, but if i got heads i will only observe the variable Y. How could i find the entropy?


title: Deriving sampling distribution

text: Assume we have the following function: $$f(p) = \frac{1}{(1-p)d}\ln\left(\frac{1}{T}\sum_{t=1}^{T}\left[\frac{1+X_t}{1+Y_t} \right]^{1-p} \right)$$ where $d$ is a constant $T$ is a constant $X_t$ for $t = 1, 2, \cdots, T$ are random variables $Y_t$ for $t = 1, 2, \cdots, T$ are random variables $p$ is defined such that it is value that satisfies $f(p) = 0$ I wish to derive a sampling distribution for the statistic $p$. By sampling distribution I mean the following: The solution to $f(p) = 0$ doesn't have a closed-form solution, but it is obvious that the resulting value of $p$ depends on $X_t$ and $Y_t$, so $p$ can be treated as a random variable that depends on the random variables $X_t$ and $Y_t$. Then for every $T$ observations of $X_t$ and $Y_t$, we have a corresponding value $p$ that satisfies $f(p) = 0$, what is the sampling distribution of $p$? Also assume you are allowed the following assumptions to achieve the above: 1) You can make any distributional assumptions regarding $X_t$ and $Y_t$, e.g., $X_t$ and $Y_t$ are independent from each other, also $X_t$, $Y_t$ for $t = 1, 2, \cdots, T$ are independently and identically distributed. 2) Rather than making distributional assumptions about $X_t$ and $Y_t$, assume you can make some assumptions about the processes $\\{X_t\\}$ and $\\{Y_t\\}$, e.g., both processes are stationary (or weakly stationary) etc. * * * EDIT: Update on my progress. So if we let $Z_t = \frac {1 + X_t} {1 + Y_t}$, then we have $f(p) = \frac {1} {(1 - p)d} \ln\left(\frac {1} {T} \sum_{i=1}^T Z_t^{1-p} \right)$. Setting this equation to $0$ and rearranging, we have: $$\sum_{i=1}^T Z_t^{1-p} = T$$ Thus, for a fixed $T$, what is the distribution of $p$?


title: Why is the posterior mode a bad Bayes estimator?

text: Suppose I have that the posterior distribution of a set of iid $X_1,...,X_n$ data that is continuous, uniformly distributed with a Pareto(a,b) prior is. Suppose that the posterior mode is: $max(x_1,...,x_n,a)$, where the mode is the maximum of the data points and the parameter a in the pareto prior distribution. Why is it that it doesn't really make sense to report this mode as the bayes estimator? Thanks!


title: How to put an exogenous variable into the ARIMA model?

text: I don't know how to put my exogenous variable in ARIMA model. I use number of tourists ('number of torism' below) in an ARIMA model and 'CLI_Index' for exogenous variable My code in R:               tourist <- ts(number of torism, start=c(2540,1),end=c(2553, 12), freq=12)     cli <- ts(CLI_Index, start=c(2540,1),end=c(2553, 12), freq=12)          # sarima (2,1,0)(0,1,1) 12 with cli     tourist.fit1 <- arima(tourist , order = c(2,1,0), seasonal = list(order=c(0,1,1), period=12) , xreg=cli)     forecast.fit1 <- forecast(tourist.fit1 , h = 1 ,xreg=cli )       I get the following error message:               > "Error in if (ncol(xreg) != ncol(object$call$xreg)) stop("Number of     > regressors does not match fitted model"      My data are available here.


title: Can I do multiple linear regression when the independent variables sum to 1?

text: I have 20 independent variables that sum to one. I have 181 vectors so my X matrix is 181 by 20 so I do have an overdetermined system. However, when I run lm() in R with this data it gives me nonsensical data with equal coefficients for every single independent variable. Also, the intercept is -1*coefficient. Could someone please explain why linear regression doesn't work in this situation?


title: Econometrics: difference-in-differences basic proof

text: At the most rudimentary level, how is a fixed effects model equal to first differencing, and how is first differencing equal to the regression $$\text{Outcome}_{it} = α_1 + α_2\text{Treat}_i + α_3\text{Post}_t + β(\text{Treat}_i\times\text{Post}_t) + ε_{it}$$ Where $\text{Treat}_i$ is a dummy which equals 1 if individual $i$ received treatment, and $\text{Post}_t$ is a dummy which equals 1 if the observation refers to post treatment time. Under what circumstances would you use the simple first difference, and when would you have to resort to the regression? Also, what are the standard error implications of using one over the other?


title: Forecasting with exponential smoothing - impact of regression to the mean? Assumptions for forecasting?

text: I have a question regarding the impact of regression to the mean when using forecasting with exponential smoothing. Unfortunately, I do not have access to the book on which this analysis was based [@Rob Hyndman, Forecasting with Exponential Smoothing: The State Space Approach (New York: Springer, 2008)], nor am I familiar with forecasting with exponential smoothing (FWES). A report was submitted to me which uses FWES to evaluate the success of a school truancy program. Numbers of absences per week were collected for each week 10 weeks prior to the intervention and again for the following 5 weeks. There is no control group, and all data are from people who received an intervention for a given number of unexcused absences. The number of unexcused absences in the 10 weeks prior to the intervention were used to forecast the expected weekly number of absences for the next 5 weeks. The slope of this line was then compared to the actual number of weekly absences. My problem with this analysis is that as far as I can tell, by definition, the intervention is going to be given at the peak, that is, the point at which the student's truancy has reached a high enough level to warrant an intervention. The weeks prior to this there would be fewer absences. As you might expect, the line in the initial 10 weeks shows a steady (jagged) increase and following the intervention, the forecasted line shows a continued increase (albeit, a slightly slower increase), while the actual number of absences plummets. The student who conducted these analyses assures me that forecasting with exponential smoothing accounts for regression to the mean, but I don't see how it can. Also, I am wondering what the assumptions are for FWES - for example, is random selection necessary, and is that a valid explanation for why these 10 weeks of data can't create a valid forecast? Or is this in fact a valid analysis? Any help would be appreciated.


title: Multinomial Naive Bayes

text: I'm looking for an article, program, algorithm that can clearly explain whats going on inside a **Multinomial Naive Bayes classifier** compared to a **Gaussian Naive Bayes Classifier**.


title: How to obtain the class conditional probability when using KNN classifier?

text: Using the KNN classifier, given a test data-point $x$, how do we get the probability of membership of $x$ to each class $y_i$, that is the probabilities $P(y_i | x)$ for $i = 1, 2, .., n$ (where $n$ is the number of classes).


title: Logistic Regression\SVM implementation in Mahout

text: I am currently working on sentimental analysis of twitter data for one of telecom company data.I am loading the data into HDFS and using Mahout's Naive Bayes Classifier for predicting the sentiments as positive,negative or neutral . Here's is what i am doing I am providing training data to the machine (key :sentiment,value:text) . Using mahout library by calculating tf-idf(Inverse Document Frequency) of text it is creating feature vector. mahout seq2sparser -i /user/root/new_model/dataseq --maxDFPercent 1000000 --minSupport 4 --maxNGramSize 2 -a org.apache.lucene.analysis.WhitespaceAnalyzer -o /user/root/new_model/predicted Splitting data as training set and testing set. That feature vector I am passing to the naive Bayes algorithm to build a model. mahout trainnb -i /user/root/new_model/train-vectors -el -li /user/root/new_model/labelindex -o /user/root/new_model/model -ow -c Using this model I am predicting sentiment of new data. This is very simple implementation what I am doing , By this implementation I am getting very low accuracy even if i have good training set . So I was thinking of switching to Logistic regression/SVM because they give better results for these kind of problem . So my question how can i use these algorithm for building my model or predicting the sentiments of tweets using these two algorithms . What steps i need to follow to achieve this ?


title: Nonnegative matrix factorization on a matrix with negative values

text: I have a gene expression profile with negative values. Is there any way to do an Nonnegative matrix factorization (NMF) on a matrix which has a few negative values? (I couldn't approximate the negatives by 1 because the range of negative values are really high.) Thanks.


title: On-line detection of over-fitting in neural networks

text: As we train a neural network, we have access to the error-rate (both on training, and test patterns). What are standard techniques to use this information to stop the learning as quickly as possible when over-fitting starts?


title: What is a "polynomially bounded" function, and why is this a requirement of the The Delta Method?

text: I am reading a paper "A note on the Delta Method" by Gary Oehlert, JASA, 1992. I am trying to estimate the variance of a function of a random variable, but first I want to understand the limitations of using the Delta Method Taylor series. Oehlert states:$^1$ > when the function we are approximating is polynomially bounded in the random > variables, then the naive Taylor series approximation will yield the correct > asymptotic approximation..." How can I tell if a function is polynomially bounded? Why is this a requirement of using the naive Taylor series approximation? Is the alternative to a naive Taylor seires a higher-order Taylor series or a different approximation? $^1$ _p 28, middle of right col, last paragraph before Examples_


title: Deviate vs. Variable

text: What is the difference between a multivariate normal random deviate and a random variable? More specifically, I suppose, what exactly is a deviate?


title: I have a discrete distribution and want to know to what extent other samples differ from it, what is the right test?

text: This is kind of a basic stats question, but I want to make sure I am doing this right. I have a distribution of objects. Specifically: `array([ 6072., 112673., 126874., 44366., 5384., 14697., 20323., 68197., 98024.,39483., 103990., 18556., 32930., 23551., 6897.])` I then have a lot of samples like [1,4,0,0,0,0...] (same length) and I'd like to know how far the samples are from the distribution above. Correlation doesn't really do it. [32,0,0,0,..] should be further away than [4,0,0,0...].


title: Time persistence in panel data

text: I am using dynamic model with panel quarter data using Stata. And my sample contain 16 nations from 2000 to 2010. Is there an approximated number of observations in the panel data to be considered as a time persistent process?


title: Estimation of change and difference in changes for repeated measurements

text: In an ongoing project concerning life style changes, I can't find the way to study the difference in effects over time between the intervention and the control group. The project is set up as a an experiment where a number of individuals have been randomly assigned to treatment or control after which the treatment group has gone through a more intense attempt at lifestyle changes. There has been followup at totally 5 time points (not equally spaced) besides baseline measurement. Due to the repeated measurements, the method chosen has been a linear mixed model in SPSS with time as a continuous variable and the baseline as a covariate. Using the "estimated marginal means" has made it straightforward to estimate the difference between the treatment and control group at chosen time points which is fine. But I would also like to estimate   1. the change _within_ each group from baseline to chosen time points;   2. the difference in change _between_ treatment and control. As for the SPSS syntax, I have basically used the following structure:               MIXED bp BY group WITH bp0 time     /FIXED=group bp0 time group*time | SSTYPE(3)     /RANDOM=INTERCEPT | SUBJECT(id) COVTYPE(VC)     /REPEATED=index | SUBJECT(id) COVTYPE(DIAG)     /EMMEANS=TABLES(group) WITH (time=3) COMPARE (group) REFCAT (FIRST).     /EMMEANS=TABLES(group) WITH (time=7) COMPARE (group) REFCAT (FIRST).      My questions could be formulated as estimation of   1. the change in response (bp) within each group between time=3 and time=7;   2. the difference in those changes  I first thought this would be fairly simple, but I am stuck. Any suggestions? Robert


title: Statistical test for difference in categorization method

text: I have a set of integer that I am trying to see use different methods to categorize them into four groups, and the 2x2 table for the outcome of the 2 methods is displayed as below:                               method_B          method_A   0   1   2   3                 0 182  11   0   0                 1  41 127   2   0                 2   0  12  18   0                 3   0   0   0   4      I am thinking of using `chisq.test` for comparing whether the difference in grouping distribution produced by 2 methods are likely to be caused by chance, but I am worrying about too many zeroes in the 2x2 table. Should I use `chi.sq` test? Or is there any equivalent `fisher.test` for my case? P.S. I tried `fisher.test` but it gives me `FEXACT ERROR 501`, which I have no idea what it is.


title: Calculating confidence interval for average hospital length of stay, case-mix adjusted, in R

text: I'm looking for sample R code, or pointers to sample R code for the following. (Gentle) critique of the approach would also be appreciated. I'm not a statistician and I'm pretty new to R. I have duration of hospitalization ("length of stay"=LOS) data for 2,000 patients from my hospital and 50,000 patients from a comparison data set. Each patient has a discharge year (YEAR) and is assigned to a diagnosis group (DG). My goal is to compare (with confidence intervals) the overall average LOS at my hospital to the expected average LOS if those patients had been part of the comparison data set. Note that for a specific year and specific DG there may be anywhere from a few to 1,000s of patients, and, of course, the distribution of LOS is not normal. There could possibly even be a DG with 1 or more patients for a given year in my hospital's data but none in the comparison data set. The approach I was considering was to create a comparison group where a random patient from the comparison data set is chosen for each patient from my hospital. The random patient would be matched by YEAR and DG. I would calculate the expected average LOS for this group and then repeat the process 10,000 times to determine the 2.5th to 97.5th percentile. I would repeat for each year and plot my hospital's average LOS versus the 95% CI for the expected average LOS. To deal with the issue of there not being a match for a patient for a given DG for a given YEAR I was thinking of loosening the match criteria to pick patients from the previous or next year. I could keep broadening the match year until there were at least N patients from which to randomly pick. Thoughts?


title: Appropriate test to compare online voting data

text: The data I want to analyze consist of set of votes similar to the voting system here on stackexchange. Votes are binary, i.e, items can receive up- or down-votes. The data have been gathered in an A/B controlled experiment. I want to compare the control group to the treatment group according to some gold standard. That is, for some items I know what the "correct" vote (according to the gold standard) should have been. However, the number of votes per item varies and can be different between control and treatment group. My hypothesis is, that the treatment reduces the voting "error". But I am not sure, which statistical test to use.


title: Ridge regression not appropriate for collinearity caused by mathematical constraints on the data

text: In this paper: Use of the Bootstrap and Cross-Validation in Ridge Regression Author(s): Nancy Jo Delaney and Sangit Chatterjee Source: Journal of Business & Economic Statistics, Vol. 4, No. 2 (Apr., 1986), pp. 255-262 The authors state that: > "As a caveat to the implementation of the ridge model, it must be pointed > out that the cement data represent a mixture problem. The four independent > variables add to approximately 100 for each observation. The ordinary ridge > model is not appropriate for collinearity caused by physical or mathematical > constraints on the data. Unfortunately, the ridge technique is often > illustrated by this data set." What I'm concerned about is whether this criticism is valid in my situation? I'm using the Nelson-Siegel-Svensson model. So I have a set of observations (which theoretically at least) are generated by the 4 factor model shown in the link above. In the extreme case that $\tau_1 =\tau_2$ then the factors relating to B2 and B3 will be perfectly collinear. (Obviously I would not envisage such a situation I am more interested in when these values are close). But for the purposes of this question suppose that they are: in that case it seems like I have a situation as described by Chatterjee and Delaney does this mean that Ridge regression is not appropriate in this case? If not what methods should I be using?


title: Overcoming small dataset anomalies in genetic algorithm

text: So I am currently making my 6th version of a model designed to predict the likelihood of of a particular medical condition based on a multifactorial genetic markers and I really would appreciate some advice with regards to a method to overcome this problem. The model I have used thus far has proved very successful, but I know there is still one thing I am not accounting for statistically. The data I have for this medical condition is particularly limited. For some parameters I have multiple candidate genetic markers lets call this one GM1, the likelihood of developing this particular medical condition is now being shown to be related to the number of repeats of this genetic marker. However when modeling them, I run into issues with limited data for some number of repeats. Parameter 1   For instance for 12 repeats of this genetic marker I could have 1000 people with 100 of them developing the condition. I input this into my model by checking what the % of people displaying 12 repeats was against the actual number that showed in a new test dataset so in this case 10% into my polynomial regression yields 8.7% chance of developing the condition. However in the same polynomial regression 35 repeats of the the genetic marker is also input as a 10% being 1 confirmed case from 10 people showing the same 8.7% chance of developing the condition. The way in which my model works is it assigns different weights to over 600 different parameters in order to best predict the actual multifactorial causation within the genetic profile. For instance parameter 2 would be a different genetic marker and a polynomial regression would be performed on the % of model data patients against the number of real observed and the same as for parameter 1, a 10% chance could input into the polynomial regression as a result of multiple different numbers of repeats of this genetic marker. What I need to know is how to adjust for different size datasets. Depending on what size dataset made up the original % chance of developing the condition. I am open to suggestion with regards to a method to achieve this.


title: Self-Study Plan for Becoming Statistical Analyst?

text: Can any one please suggest me a self-study plan for becoming a Statistical Analyst? I need suggestion something similar to Self Study Plan for Quantitative Analyst


title: Regression analysis in R using text field?

text: I'm working in R. I'd like to run a regression analysis for predicting price against terms in a text field. I have a dataset of jewellery auction listings, with price paid, date, and an unstructured description of the item type:               text,date,price_usd     "Ruby necklace, Spanish",1925,45000     "Diamond ring, 0.7 carat, bezier cut",1972,24000     "Diamond necklace",1980,87000     ...      I know how to run a linear regression for price against date:               data <- read.csv('jewels.csv')     lm1 <- lm(data$price~data$date)     summary(lm1)      Now what I'd like to do is build a similar model, using the words in the description field that are most associated with higher prices. Intuitively I'd guess these include "diamond" and "necklace", while (say) "amethyst" and "ring" were associated with lower prices, but is there a way I can build a model to look at this? My sense is that I need to do the following things:   * turn the text field into a bag of words (vector)   * remove stop words   * normalize each word for overall count(?)   * run some kind of regression against price.  I'd really welcome some guidance on how to approach each step.


title: I want a book about Normal distribution and Multivariate normal distribution

text: I have a problem about normal distribution and multivariate normal distribution .so I need a whole books on multivariate normal distributions and normal distribution. Thanks for help.


title: Variation within and between groups

text: So in my data, I have parameters(~800) and each parameter has 14 different groups. I individually know the mean and standard deviation within each group for each parameter. The data within each group has been collected with 30 measurements. Eg. http://imgur.com/nFDDRK5 "Data" I am looking to find out the variation within each group and between each group for all the parameters. I know that I can use one-way ANOVA to calculate the required variances, but given the quantity of data this doesn't seem feasible. Could I use the following approach ?   1. For calculating the variance within each group - Since I know the standard deviation for each group individually, a mean of these standard deviations may be a good measure of the standard deviation within each group ?   2. For calculating the variance between groups - Using the means of each group for every parameter, use the IQR/1.34 to calculate the standard deviation ? Since the number of observations >= 30 for each group, by CLT IQR/1.34 should give a good description of the standard deviation given the normal nature of the means ?


title: Test Statistics

text: My second question is about Test Statistics. The questions and answers are on this PDF: http://www.mediafire.com/?b74e633lxdb49rb I understand that they work out the join density and I know how to calculate this. What I really don’t understand is how they figure out the test statistic. They randomly multiply by 1 in some cases and in some cases they just pick a sum. I don’t understand WHY that’s the test statistic - is it a guess? Also, how do you check a test statistic you have calculated is correct?


title: How to implement formula for a independent groups t-test in C#?

text: I am writing a program in C# that requires me to use the Ttest formula. I have to effectively interpret the Excel formula:               =TTEST(range1,range2,1,3)      I am using the formula given here and have interpreted into code as such:                double TStatistic = (mean1 - mean2) / Math.Sqrt((Math.Pow(variance1, 2) /              count1) + (Math.Pow(variance2, 2) / count2));      However, I don't fully understand t-test and the values I am getting are completely different than those calculated within Excel. I have been using the following ranges:               R1:     91.17462277,     118.3936425,     96.6746393,     102.488785,     91.26831043          R2:     17.20546254,     19.56969811,     19.2831241,     13.03360631,     13.86577314      The value I am getting using my attempt is 1.8248, however that from Excel is 1.74463E-05. Could somebody please point me in the right direction?


title: Modelling biases in selection of one choice from many using count data

text: In a biological experiment I draw somewhere between 100k and several million datapoints, where each datapoint is one of 1024 DNA sequences. If each sequence has an equal probability of coming up, I would expect the distribution of counts for each sequence to be binomially distributed with p=1/1024. Given the very large sample size, I would expect this to look roughtly normal. However, it looks more like a log normal distribution. I tested this using a multinomial test (`chisq.test` in R) and got p<1e-16. In some experiments the top 10% most common sequences make up over 50% of all sequences sampled. I would like to try to model what factors are affecting the choice of sequence, specifically features of the sequence itself (starting with the % of the sequence made up of G and C bases). The data looks something like this:                      seq seq_gc count     5119 GCGTT    0.6    21     5120 AAATG    0.2   160     5121 GCCCG    1.0   230     5122 GCCCA    0.8   222     5123 AAATC    0.2   245     5124 GCCCC    1.0   380     5125 AAATA    0.0   162     5126 CCTAT    0.4   216     5127 TTAGC    0.4    18     5128 GCCCT    0.8   305      Looking at plots of %GC vs count, it looks like there is a relationship between log(count) and %GC for some experiments, but not for others. However, its difficult to tell because there are a limited number of values that %GC can take with only 5 bases, and some values are much more common than others. Clearly a normal least sequares regression isn't suitable here, but what would be? I tried poisson, binomial and negative-binomial generalized linear models, but they all give different answers and I don't know which one is "correct". I've looked at some diagnositc plots for them, and in all cases the standardized residules are far from normally distributed (would you expect normality in these models? Does this just mean that my explanitory variable isn't very explanatory?), and have some degree of dependace of residuals on position (how much is too much?).


title: Pivot functions

text: I have the following problem: > Suppose we have a sample $Y_1,...,Y_n \sim N(\mu, \sigma^2)$. Find the > pivotal function for $\sigma^2$ when $\mu$ is known and when $\mu$ is > unknown. I know that you need to find a statistic where its distribution doesn't depend on the parameter. For $\mu$ unknown I'm assuming you use $(n-1)/S^2$ which is a chi-squared distribution ($n-1$ degrees of freedom). For $\mu$ known, what do you use? And why? In my notes they say $(\bar Y - \mu) / (\sigma\sqrt{n}) \sim N(0,1)$ but then it goes on to say something about summing the $Y_i$'s. Why?


title: What are correct values for precision and recall when the denominators equal 0?

text: Precision is defined as: > p = true positives / (true positives + false positives) What is the value of precision if (true positives + false positives) = 0? Is it just undefined? Same question for recall: > r = true positives / (true positives + false negatives) In this case, what is the value of recall if (true positives + false negatives) = 0? P.S. This question is very similar to the question What are correct values for precision and recall in edge cases?.


title: Variance within each cluster

text: I have done some clustering to a matrix with 30 random variables , each variable has 13000 observations ). i got 10 clusters and now i need to test how good the clustering is by calculating the variance in each cluster. does anyone knows how can i calculate the variance? i can easily calculate the variance of each column in my matrix (e.g the variance of each random variable) but i want to calculate the variance of the whole cluster. does anyone know how it can be done? e.g.               data <- data.frame(x=c(2,2,2,3,7),                    y=c(30,40,40,30,10),                    z=c(1,2,3,4,5),                    cluster=c('a','a','c','a','c'))          candidates <- dlply(data,.(cluster),function(data){      laply(data[,-4],var)     })      This gives variance per column for each cluster label (a,c). I don't think it's the right approach.


title: KL divergence between two univariate Gaussians

text: I need to determine the KL-divergence between two Gaussians. I am comparing my results to these, but I can't reproduce their result. My result is obviously wrong, because the KL is not 0 for KL(p, p). I wonder where I am doing a mistake and ask if anyone can spot it. Let $p(x) = N(\mu_1, \sigma_1)$ and $q(x) = N(\mu_2, \sigma_2)$. From Bishop's PRML I know that $$KL(p, q) = - \int p(x) \log q(x) dx + \int p(x) \log p(x) dx$$ where integration is done over all real line, and that $$\int p(x) \log p(x) dx = \frac{1}{2} (1 + \log 2 \pi \sigma_1^2),$$ so I restrict myself to $\int p(x) \log q(x) dx$, which I can write out as $$-\int p(x) \log \frac{1}{(2 \pi \sigma_2^2)^{(1/2)}} e^{-\frac{(x-\mu_2)^2}{2 \sigma_2^2}} dx,$$ which I can separate into $$\frac{1}{2} \log (2 \pi \sigma_2^2) - \int p(x) \log e^{-\frac{(x-\mu_2)^2}{2 \sigma_2^2}} dx.$$ Taking the log I get $$\frac{1}{2} \log (2 \pi \sigma_2^2) - \int p(x) -\frac{(x-\mu_2)^2}{2 \sigma_2^2} dx,$$ where I separate the sums and get $\sigma_2^2$ out of the integral. $$\frac{1}{2} \log (2 \pi \sigma_2^2) + \frac{\int p(x) x^2 dx - \int p(x) 2x\mu dx + \int p(x) \mu^2 dx}{2 \sigma_2^2}$$ Letting $\langle \rangle$ denote the expectation operator under $p$, I can rewrite this as $$\frac{1}{2} \log (2 \pi \sigma_2^2) + \frac{\langle x^2 \rangle - 2 \langle x \rangle \mu_2 + \mu_2^2}{2 \sigma_2^2}.$$ We know that $var(x) = \langle x^2 \rangle - \langle x \rangle ^2$. Thus $$\langle x^2 \rangle = \sigma_1^2 + \mu_1^2$$ and therefore $$\frac{1}{2} \log (2 \pi \sigma^2) + \frac{\sigma_1^2 + \mu_1^2 - 2 \mu_1 \mu_2 + \mu_2^2}{2 \sigma_2^2},$$ which I can put as $$\frac{1}{2} \log (2 \pi \sigma_2^2) + \frac{\sigma_1^2 + (\mu_1 - \mu_2)^2}{2 \sigma_2^2}.$$ Putting everything together, I get to \begin{align*} KL(p, q) &= - \int p(x) \log q(x) dx + \int p(x) \log p(x) dx\\\ &= \frac{1}{2} \log (2 \pi \sigma_2^2) + \frac{\sigma_1^2 + (\mu_1 - \mu_2)^2}{2 \sigma_2^2} - \frac{1}{2} (1 + \log 2 \pi \sigma_1^2)\\\ &= \log \frac{\sigma_2}{\sigma_1} + \frac{\sigma_1^2 + (\mu_1 - \mu_2)^2}{2 \sigma_2^2}. \end{align*} Which is wrong since it equals $1$ for two identical Gaussians. Can anyone spot my error? **Update** Thanks to mpiktas for clearing things up. The correct answer is: $KL(p, q) = \log \frac{\sigma_2}{\sigma_1} + \frac{\sigma_1^2 + (\mu_1 - \mu_2)^2}{2 \sigma_2^2} - \frac{1}{2}$


title: Conditional distribution of quadratic forms

text: Given that $Y$ follows multivariate normal distribution ,i.e, $N_n (0, \sigma^2 I_n)$, we want to find the distribution of $Y'Y$ given that $a'Y=0$ where $a$ is a non zero constant vector. I know that the distribution of $Y'Y$ would be $\sigma^2 \chi_n^2$ if the condition is not given. How to approach for the given condition?


title: Limiting distribution of Hansen Sargan J statistic

text: I was reading about instrumental variables, where I found description of the Hansen Sargen J statistic. The model is the usual $y=X\beta+\epsilon$, where $X\in \mathbb{R}^{N\times K}$, $E(X^T\epsilon)\ne 0$ and $Z\in \mathbb{R}^{N\times M}, M>K$ and $E(Z^T\epsilon)= 0$. Assume $\epsilon\sim \mathcal{N}(0,\sigma^2I)$. The Hansen Test looks at the overidentifying restrictions $Z^T(y-X\beta_{2SLS})$ and checks if they are close to zero. The J-statistic is $$J=(y-X\beta_{2SLS})^TZVar[Z^T(y-X\beta_{2SLS})]^{-1}Z^T(y-X\beta_{2SLS})$$ I find several references that claim $J \sim \chi_{M-K}^2$, but I could not find a proof (it was not clear to me right away, anyway). The closest result for standard normal $\epsilon$ is that $\epsilon^T A \epsilon\sim \chi_{K}^2$ if and only if $A$ is idempotent and of rank $K$. Could anyone help me with the proof?


title: Comparing minitab output to R

text: I am trying to recreate analysis of diamond data in R and compare it to minitab output cited in http://www.amstat.org/publications/jse/v9n2/datasets.chu.html I'm almost certain its the same data set yet my results are different. I suspect I'm modeling it differently. I also noticed that I'm missing colourD and some others. Did the the minitab model somehow pick different dummys to exclude? **R OUTPUT**               require("Ecdat")          data(Diamond)          Diamond$ln_price <- log(Diamond$price)          summary(lm(ln_price~carat+colour+clarity+certification, data=Diamond))               Call:     lm(formula = ln_price ~ carat + colour + clarity + certification,          data = Diamond)          Residuals:          Min       1Q   Median       3Q      Max      -0.31236 -0.11520  0.01613  0.10833  0.36339           Coefficients:                        Estimate Std. Error t value Pr(>|t|)         (Intercept)       6.8011915  0.0496111 137.090  < 2e-16 ***     carat             2.8550130  0.0369676  77.230  < 2e-16 ***     colourE          -0.0295093  0.0404611  -0.729  0.46638         colourF          -0.1063582  0.0379807  -2.800  0.00544 **      colourG          -0.2063494  0.0389848  -5.293 2.35e-07 ***     colourH          -0.2878756  0.0394752  -7.293 2.81e-12 ***     colourI          -0.4165565  0.0413818 -10.066  < 2e-16 ***     clarityVS1       -0.2019313  0.0310634  -6.501 3.41e-10 ***     clarityVS2       -0.2985406  0.0333027  -8.964  < 2e-16 ***     clarityVVS1      -0.0007058  0.0311121  -0.023  0.98192         clarityVVS2      -0.0966174  0.0289396  -3.339  0.00095 ***     certificationHRD -0.0088557  0.0208641  -0.424  0.67155         certificationIGI -0.1827107  0.0249516  -7.323 2.33e-12 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 0.1382 on 295 degrees of freedom     Multiple R-squared: 0.9723,     Adjusted R-squared: 0.9712      F-statistic: 863.6 on 12 and 295 DF,  p-value: < 2.2e-16       **MINITAB OUTPUT**               ln_price = 6.08 + 2.86 Carat + 0.417 D + 0.387 E + 0.310 F + 0.210 G + 0.129 H               + 0.299 IF + 0.298 VVS1 + 0.202 VVS2 + 0.0966 VS1 + 0.0089 GIA                - 0.174 IGI          Predictor        Coef       StDev          T        P            Constant      6.07724     0.04809     126.37    0.000     Carat         2.85501     0.03697      77.23    0.000           D             0.41656     0.04138      10.07    0.000            E             0.38705     0.03082      12.56    0.000            F             0.31020     0.02748      11.29    0.000            G             0.21021     0.02836       7.41    0.000            H             0.12868     0.02852       4.51    0.000            IF            0.29854     0.03330       8.96    0.000            VVS1          0.29783     0.02810      10.60    0.000            VVS2          0.20192     0.02534       7.97    0.000            VS1           0.09661     0.02492       3.88    0.000            GIA           0.00886     0.02086       0.42    0.672            IGI          -0.17385     0.02867      -6.06    0.000                 S = 0.1382      R-Sq = 97.2%     R-Sq(adj) = 97.1%          Analysis of Variance          Source            DF          SS          MS         F        P     Regression        12     197.939      16.495    863.64    0.000     Residual Error   295       5.634       0.019     Total            307     203.574      **UPDATE** Glen_b's solution produces the correct results, but the output is not desirable.                                      Estimate Std. Error t value Pr(>|t|)         (Intercept)           6.3846350  0.0411520 155.148  < 2e-16 ***     carat                 2.8550130  0.0369676  77.230  < 2e-16 ***     C(colour, base = 6)1  0.4165565  0.0413818  10.066  < 2e-16 ***     C(colour, base = 6)2  0.3870472  0.0308241  12.557  < 2e-16 ***     C(colour, base = 6)3  0.3101983  0.0274791  11.288  < 2e-16 ***     C(colour, base = 6)4  0.2102072  0.0283593   7.412 1.32e-12 ***     C(colour, base = 6)5  0.1286809  0.0285231   4.511 9.31e-06 ***      Is there a way to get the original factor character displayed instead of the level? **UPDATE** Using relevel to set the base to the same bases as the minitab produces the same results and is still readable.                 Diamond$colour <- relevel(Diamond$colour, ref="I")            Coefficients:                          Estimate Std. Error t value Pr(>|t|)           (Intercept)       6.3846350  0.0411520 155.148  < 2e-16 ***       carat             2.8550130  0.0369676  77.230  < 2e-16 ***       colourD           0.4165565  0.0413818  10.066  < 2e-16 ***       colourE           0.3870472  0.0308241  12.557  < 2e-16 ***       colourF           0.3101983  0.0274791  11.288  < 2e-16 ***       colourG           0.2102072  0.0283593   7.412 1.32e-12 ***       colourH           0.1286809  0.0285231   4.511 9.31e-06 ***


title: Is the median score greater than 5?

text: A consumer group rated 12 manufactured products import ed from a single source on a scale of 1 to 10. The data are shown below. Do the data provide sufficient evidence to indicate tha t the median score is greater than 5? Use the sign test at the 1% level of significance. 2 9 6 7 8 5 10 2 8 5 3 4


title: Root mean square vs average absolute deviation?

text: Both Root Mean Square and Average absolute deviation seem like the measures of the magnitude of variability (especially when the variates are both +ve and -ve). What are the rules of thumb to choose one of them over the other?


title: How to find log-likelihood of multiple sequences for hmm using Kevin Murphy toolkit for MATLAB

text: I have an observation sequence of TPM, EPM and prior. I want to find the log- likelihood of around 100 sequences of length 10 at a time. How can I do this using a forward algorithm?


title: Confidence Interval for predictions for Poisson regression

text: This is a follow-up question from this post, here: Confidence intervals for predictions from logistic regression The answer from @Gavin is excellent, but I have some additional questions which I think would be useful for others. I am working with a Poisson model, so basically it is the same approach described in the other post, only `family=poisson` instead of `family=binomial`. To my first question: @Gavin writes:               mod <- glm(y ~ x, data = foo, family = binomial)     preddat <- with(foo, data.frame(x = seq(min(x), max(x), length = 100))     preds <- predict(mod, newdata = preddata, type = "link", se.fit = TRUE)      What is the point of the second line there? Is it necessary to create a data.frame with minimum and maximum of the explanatory variable? Could I not, for some explanatory variable(s) `x` (stored in some data frame `data`), just go from the first line and directly to the third? To my second question: In the beginning of his answer @Gavin writes: > The usual way is to compute a confidence interval on the scale of the linear > predictor, where things will be more normal (Gaussian) and then apply the > inverse of the link function to map the confidence interval from the linear > predictor scale to the response scale. Why are "things" more normal on the scale of the linear predictor(s)? Is this also the case when I do my Poisson regression? I assume the reason for using critical value 1.96 when constructing the CI's, is because of the assumptions that "things" are normal. Can somebody explain this further? My third question: Is there a relationship between the standard deviation which we get by using `se.fit=TRUE` in predict() and the standard deviations of the coefficients of the explanatory variables, which we simply get from `summary(mod)`? (`mod` is some glm object)