title: Cramer-von Mises using pdf instead of Cdf

text: Cramer-von Mises gives a nice way to compare distributions with: $$ \int [(\hat{F}_n(x) - F_0(x)]^2dF_0(x); $$ I'd like to use something like this to compare two empirical distributions, is there a similar statistic that instead uses: $$ \int_n[\hat{f}_n(x) - f_0(x)]^2df_0(x); $$ That is, is there something special about using the CDFs as opposed to the densities?


title: understand forecasts in linear state space models

text: The Kalman Filter provides the one-step-ahead forecasts within the recursions. We start estimating the (unkown) variance of the parameters for instance through `MCMC` or `MLE`. Then we use the variance estimates (e.g. mean) to run the Kalman filter. But, since we have already used the data to estimate the variances, can we compute forecasts using the same data again? this is what I have understood the `dlmFilter` function from the `dlm` package does and would like to know why this is correct.


title: Detecting/recognizing landscape and architecture photographs

text: I have a large set of photos of various kinds, labeled with some tags but the tags cannot really be trusted (e.g., some portraits will have a tag 'portrait' while many won't, likewise for 'landscape', etc. tags). I'd like to somehow recognize which of them are of landscape and/or architecture (as opposed to humans, whether faces or full-length, cars, food, etc.). I appreciate this is a complex topic, and I don't really have the skills or the time to address it well, so I just hope to use some open-source libraries and hope for the best, as I don't need high accuracy - both false positives and false negatives are not a disaster. I also don't care about what building is in what picture or which mountain, etc. My current idea is to run some existing face detection algorithm (perhaps Haar using OpenCV) on it, and skip all the ones which return positives. Is there anything else relatively simple which I could do to skip pictures of people from the back, plates of food, etc. (without training my own classifiers)? Is my thinking in the right direction or am I missing some easier solutions?


title: Temporal trend visualizations for regions

text: I try to reproduce a plot which I found in a publication (Mayor-Fernández et al 2012) showing unemployment rates across time for several regions. It combines a heatmap (?), a time series plot (lower panel) and boxplots (right panel). Sorry for this basic question but I am new to data analysis and R. For a better impression of the plot please see page 7 in http://www2.dse.unibo.it/wp/WP835.pdf Can someone guide me in a correct direction or share some code? I appreciate your help and time very much!


title: Combining a linear Kalman Filter with additional linear constraints?

text: This question contains a relatively long prelude, since I want to explain as clearly as possible the motivation for the question. It may well be the case that I am asking the wrong question (i.e. there is a better way of solving the problem than what I am attempting) and I would welcome any good answers along those lines. I have implemented a real time map for a physical quantity over a large region (a continent) that combines many hundreds of data points and uses a Kalman Filter to jointly estimate the coefficients of a linear combination of basis functions (the map) plus a large number of instrumental biases and other nuisance parameters. A problem with this approach is that the data points are distributed in a very uneven way; some regions have a lot of closely spaced data points and other regions are much more sparse. Since the map is represented by a global set of basis functions, there is no real happy medium in choosing the best order of expansion for the generating function. I either have a too high order for the sparse regions leading to spurious artefacts, or too low order for the data rich regions, smoothing out real structure that is supported by the data. In order to be more adaptive to local data density, I am thinking of using some kind of radial basis function (RBF) approach utilising unsupervised clustering (e.g. K-means) to distribute a set number of basis function origins to match the local data density. I still need to use some kind of filtering to estimate the nuisance parameters along with the map, but this RBF approach maintains the linearity of the system so the Kalman Filter can be used. I'd 'just' need to experiment with different forms of the basis functions and their parameters, the number of basis functions to use etc. However, it would be ideal to also incorporate a smoothness constraint, as in the case of a thin plate spline (TPS). This can be expressed as an additional set of linear constraints on the map values at each basis origin and the standard TPS formulation has a closed form solution. What I would like to do is incorporate these linear constraints into the filtering mechanism, as is required due to the noisiness of some parts of the system, rather than the one off minimisation of the TPS. I can vaguely conceive of how this might be done, but my linear algebra skills are not up to the task of fully deriving this (or even demonstrating that it can be done). So my question is whether this approach is feasible and if so what the form of the resulting filter looks like. Is there a specific piece of terminology that succinctly describes what I am trying to do? Is there a standard result out there for this situation?


title: Chi-squared test for nonlinear regression

text: In the context of a statistical signal processing problem, I have a signal model of the form $x[n]=s[n;\theta]+w[n]$ where $w[n]\sim{\cal N}(0,\sigma^2)$ and $s[n;\theta]$ depends nonlinearly on $\theta$. $n=0,\ldots,N-1$ is the time sample index. I compute the Maximum Likelihood estimate of $\theta$ as $\theta_{\rm ML}=\underset{\theta}{\arg\min}\,\,p(x[0],x[1],\ldots,x[N-1]|\theta)$ by minimizing the log-likelihood. Now my question is: can I use the $\chi^2$ distribution for a goodness-of-fit test in this case? Most statistics books I have seen present the $\chi^2$ test in the context of linear regression. I am especially interested in the general multivariate case ($\boldsymbol{\theta}$ is a vector) and $\mathbf{w}[n]\sim{\cal N}(0,\mathbf{C})$. Also, can someone please recommend a good book (at the level of Casella and Berger) presenting relevant material?


title: How do I model monthly count data (the number of lumbar fusion surgeries)?

text: I have the following data:   * Month-year number of people with back/neck problem    * number of people using surgery (lumbar fusion) as a treatment The goal is to compute usage rate of surgery (`=number of surgeries / number of people with back/neck problem`) for each month and hence come up with monthly forecast. What is the best approach to model this data? I tried to figure out seasonality in the data, but with little success. Is it possible that my data does not have any seasonality? If I use Poisson regression, then do I have to test for normality? Or can I just ignore the distribution of the data. A step by step procedure is requested as I am new to this field. I have data from January 2006 to May 2012. My counts are anywhere between 1 to 70. I am using sas and sas jmp. There is a limitation on the version of sas that I use. It has no time series license. * * * There is no seasonality in number with back/neck problem. So, I guess there will be no seasonality in number of surgeries either. Since, I have count data I am using Poisson Regression to model number of surgeries. But Poisson regression has some strong assumptions. Another possible method is Negative Binomial Regression. Currently my model contains a lagged dependent variable (lagged number of surgeries) and number with back/neck pain. I hope I am on the right path. Or do you suggest a better approach?


title: Generate data from a t-distribution with specified mean and standard deviation

text: How does one randomly sample from a T-distribution in R. From what I've found, the function `rt` in R doesn't let you specify the mean and standard deviation. For a normal distribution it is simply `rnorm(x,mu,sd)`. EDIT: For the t-distribution there is a _central_ and a _non-central_ version. I want to know the difference between the two. In addition, if I want to specify the mean and standard deviation, does that automatically mean I am dealing with the non-central version of the t-distribution? I chose to use the t-distribution because the data I am using are rather fat tailed and a t-distribution with a low degrees of freedom seems like a good idea. What other distributions are there for handling fat tailed data? Also it would be great if you can specify the function in R too.


title: Missing value treatment

text: I have a data set with 18% of AGE variable missing which is an important variable for analysis.   1. Should I try regression imputation or should I drop those observations?    2. Does even regression imputation make sense in this case (for age)??   3. I also have income variable but the correlation between age and income is negative and strength is .1  What should I do?


title: How to determine whether data is slightly or extremely non-normally distributed?

text: I'm a PhD student and doing a research on regression analysis. My question is how to determine whether the data is slightly, moderately or extremely non-normally distributed? * * * TQ to all responses of my question. But, may be my question is not so clear. Ok, let say i have some different values of skewness and kurtosis (for example: skewness = 1.5, kurtosis = 2.0). So my question is, from the values of skewness and kurtosis, what is the type of data distribution? is it moderately non-normal, or extremely non-normal or what?


title: Can the predicted value vs observed value plot have a slope not equal 1 in a LASSO model?

text: I was trying to use glmnet package in R to create a lasso regression model. The details of my data are: Dependent variables $y$: 451 observations, single value for each observation. Independent variables $x$: 451 observations, 959 values for each observation. I find the optimal value of $\lambda$ by minimizing the cross-validated mean- squared error. Then I generated the parameters $\beta$ by the respective value of $\lambda$. The number of nonzero parameters is around 20. Finally I generated the predicted values $\hat{y}=\sum\beta x$. When I make a plot of $\hat{y}$ against $y$, the slope of the plot is smaller than 1. But If I plot $y$ against $\hat{y}$ , the slope is equal to 1. My questions are:   1. Does that mean my model is flawed? Why the slope of $y$ against $\hat{y}$ is not equal to the inverse of the slope of $\hat{y}$ against $y$?   2. Is that due to the nonlinearity of my data? I know the relationship between $y$ and $x$ is nonlinear.   3. Does that mean my model can be improved by doing some nonlinear transformation on $y$? I was tried to do some sigmoid transform by it doesn't changed the slope. Many thanks.


title: Multiple cross validation consistently performing worse than chance

text: I am calculating accuracy results from a single feature using LDA. I am very curious why a certain feature consistently performs at less than chance when calculating accuracy using a stratified 10-fold cross validation technique. I have an equally distributed binary class and performing multiple xvals results in an accuracy that fluctuates fairly tightly around 0.3. This means the training data is consistently setting a boundary which when tested against is basically the opposite of what was trained on. I am performing this stratified 10-fold xval 10-20 times and I would expect this to fluctuate around 0.5, but it is fairly static around the 0.3 to 0.4 area. Thanks for the help!


title: Likelihood Ratio Testing in Multinomial logistic regression using SPSS

text: I am doing multinomial logistic regression. Dependent variable is categorical with 6 categories. 6 Independent variables are measured on an interval scale (included as covariates in SPSS) Now I am watching the LR test. None of the independent variables have a $\chi^2$ value and no degrees of freedom are shown, so significance is not specified. (except for the constant, which is significant and has 5 df's). One variable is redundant in the parameter estimates. When I omit the redundant IV in the model, there are $\chi^2$values, df's and significance values for the independent variables. Was the first model over- parametrised or something? what is the explanation for this situation? And is it valid to omit one independent variable in the model?


title: How to inverse the empirical function?

text: I fitted a empirical distribution to a set of time series data (Y) by following code in R:               Ye=rank(Y)/(length(Y)+1)      How we can find the inverse of this distribution? Thanks


title: How can a t-test be statistically significant if the mean difference is almost 0?

text: I am trying to compare data from 2 populations to tell if the difference between the treatments is statistically significant. The data sets appear to be normally distributed with very little difference between the two sets. The average difference is 0.00017. I performed a paired t-test, expecting that I would fail to reject the null hypothesis of no difference between the means, however, my calculated t-value is much higher then my critical t-value.


title: How to construct a design matrix for coxph with pspline term?

text: I am wondering how to reconstruct the design matrix for a `coxph()` model with a `pspline()` term. For example, if I fit the following model,               fit <- coxph(Surv(t,delta) ~ pspline(x,df=0))      how can I calculate the linear predictor by hand for a given value of $x$? Formally, I am modeling $h(t|x) = h_0(t)\cdot e^{(X^\intercal\beta)}$, and I want to know how to calculate $X$ or $X^\intercal\beta$ by hand. I can't find detailed documentation of how the design matrix is constructed when the `pspline()` term is used. As for why I care, I am running a large set of new predictions from a `coxph()` object with a `pspline()` term. I have stored the cumulative hazard $H_0(t)$ and want to be able to calculate $H_0(t)\cdot e^{(X^\intercal\beta)}$ efficiently. I'm finding that using               predict(fit, newdata=data.frame(x=x.new), type='lp')      is very inefficient! Any tips would be greatly appreciated.


title: Subtree replacement vs subtree raising

text: As we saw in this question, the recommended strategy of building a decision tree is postpruning. The two methods for that are _subtree replacement_ and _subtree raising_. At each node, an algorithm decides whether it should perform subtree replacement, subtree raising, or leave the subtree as is.   * _Subtree replacement_ selects a subtree and replaces it with a single leaf.    * _Subtree raising_ selects a subtree and replaces it with the child one (ie, a "sub-subtree" replaces its parent) What are the pros and cons, when should we use each?


title: Difference between dynamic programming and temporal difference learning in reinforcement learning

text: In reinforcement learning, what is the difference between dynamic programming and temporal difference learning?


title: Automatic data cleansing

text: A common problem is ML is poor quality of the data: errors in feature values, misclassified instances, etc etc. One way of addressing this problem is to manually go through the data and check, but are there other techniques? (I bet there are!) Which ones are better and why?


title: Tracking and data association with Kalman filters

text: I am trying to solve tracking problem. At certain points in time I receive object location and I should make decision whether received object location belongs to existing track or not. If not, I should create new track. Every track is filtered by Kalman filter so I can calculate speed and posterior. What algorithm or family of algorithms suit best? Maybe nearest neighbour?


title: Decision trees and backward pruning

text: By building the complete tree and pruning it afterward we are adopting a strategy of postpruning (or backward pruning) rather than prepruning (or forward pruning). Prepruning would involve trying to decide during the tree-building process when to stop developing subtrees. It's attractive because allows to avoid all the work of developing subtrees only to throw them away afterward. What are the benefits of post-prunning then?


title: What's the minimum number of runs needed to cover a certain set?

text: I have a set of 40,000 numbers. If I randomly draw 2000 numbers, how many draws do I need to cover the whole set? Kindly give the theory behind the answer as well. The problem is related to my work, I have a set of 40,000 genes and i can only build a 2000 gene network from this set. I want to cover all the possibilities of the networking between these genes. so if we say that the confidence interval is 95 % sure that i cover the data set and p=0.01 , would that change the result?


title: How to deal with dirty data without decreasing the generalization power of the model?

text: I'm working on a sequence labeling task using a conditional random field. My training data has some errors in it. These are not errors in the sense of having tokens that have incorrect labels. Instead, some of the tokens have typos in them. Not surprisingly, the labeler tends to get incorrect results when testing on the data that contains typos. I can't decide if it would be better to leave the typos in there or take them out. A possible benefit of keeping the typos as they are is that the training data would then more accurately represent the real world data that the labeler will be used on. On the other hand, I fear that training on lower quality data will result in a lower quality labeler. Is there a standard practice or conventional wisdom regarding this sort of problem? Forgive me this question is considered too subjective for this forum. I just discovered this site today.


title: Bias-variance tradeoff -- bias or variance effect

text: I understand that supervised learning is associated with an error that can be split between bias and variance: $ MSE = b^2 + var $ What does bias and variance account for intuitively ? I understand that   * bias is a measure of model prediction accuracy: bias is important if model learnt on several dataset systematically fails to predict output from another input.   * variance is the sensitivity to training dataset of the models in their capability to predict an output from a new input. In other words, if training a model on different training datasets gives very different output for a given input dataset, there is high variance. If a model trained on a wide spectrum of inputs fails to predict, it is too much 'stretched towards the training data', it is biased. How is it correct until there ? Also, if my model trained on some data fails to predict (bias), doesn't it mean that when trained on different datasets it will give very different predictors (variance) ? Would you hint to a reference to clarify ?


title: Estimate the minimal values of a gaussian normal distribution

text: My first question here and i hope that its not to obvious or stupid.   My problem: I have a number of $x$ values ( $x$ can vary) that follow a normal distribution. I can calculate all normal descriptive statistics (minimum, mean, quartiles and so on), but those $x$ values only describe a small fraction of a potentially much bigger range of values (which still follows the same distribution and the average should not vary much).   So what i really want is to estimate the minimal value of the unknown full distribution. Something like an estimation of the minimal x-intercept of a known distribution with given slope and tip value. Imagine for instance the standard `trees` dataset in R (see picture below). I don't want to know what is the minimal value of the given height values (60), but instead estimate the possible true minimal value in a most accurate way (might be $height<=50$). ![enter image description here](http://i.stack.imgur.com/GJzoN.png) Is this possible or just stupid? Or am i missing a basic statistic principle which already allows me to estimate those values?


title: Defining a multiple regression model?

text: So I have an assignment and I'm not sure what exactly they mean by defining a multiple regression model. Well I know what it is but how do I define one ? I need to use CARPRICE as the response variable and the other 3 predictors. Also what are the model assumptions ? Thanks !


title: Ridge Regression in R where coefficients are penalized toward numbers other than zero

text: Is it possible to penalize coefficients toward a number other than zero in a ridge regression in R? For example, let's say I have dependent variable Y and independent variables X1,X2,X3, and X4. Because of the multicollinear nature of the ivs, ridge regression is appropriate. But say I'm fairly certain that the coefficient of X1 is near 5, X2 is near 1, X3 is near -1, and X4 is near -5. Is there a ridge package and method in R where I can implement penalties on the coefficients of the ivs that penalize them toward those numbers instead of 0? I'd love to see an example in R with my example data, if possible. Thank you.


title: Plotting error ellipsoid from 3x3 covariance matrix in R?

text: I'm hoping to be able to take a 3x3 covariance matrix and turn this into an error ellipsoid but so far I haven't been able to achieve this. I'm very new to R (in fact turned to it to attempt to solve this problem) - and after finding a working example for a 2x2 matrix I've not managed to modify this for 3x3 successfully - can anyone point me in the right direction as to which changes are required please? The working code for the 2x2 matrix (found here in another question) is below:               ctr    <- c(0, 0)                                    A      <- matrix(c(2.2, 0.4, 0.4, 2.8), nrow=2)      RR     <- chol(A)                                    angles <- seq(0, 2*pi, length.out=200)               ell    <- 1 * cbind(cos(angles), sin(angles)) %*% RR       ellCtr <- sweep(ell, 2, ctr, "+")                         eigVal  <- eigen(A)$values     eigVec  <- eigen(A)$vectors     eigScl  <- eigVec %*% diag(sqrt(eigVal))       xMat    <- rbind(ctr[1] + eigScl[1, ], ctr[1] - eigScl[1, ])     yMat    <- rbind(ctr[2] + eigScl[2, ], ctr[2] - eigScl[2, ])     ellBase <- cbind(sqrt(eigVal[1])*cos(angles), sqrt(eigVal[2])*sin(angles))      ellRot  <- eigVec %*% t(ellBase)                                               plot((ellRot+ctr)[1, ], (ellRot+ctr)[2, ], asp=1, type="l", lwd=2)     matlines(xMat, yMat, lty=1, lwd=2, col="green")     points(ctr[1], ctr[2], pch=4, col="red", lwd=3)      Many thanks in advance for any help.


title: Comparing proportions between two mega-groups

text: Say I sit in 2 different restaurants, Cheesecake Factory and Panda Express, that has nothing in common on their menu. I randomly order items from the menu and record how much sugar it has as well as the item is an appetizer or entree or dessert. At each restaurant, I stop ordering 30 minutes after when I order my first item. At the end of the day I have a data set that may look like the picture below... I am wondering how to compare the % sugar in my orders from Cheesecake Factory to that in the ones from Panda Express? I understand that the items from the two restaurants are different and may be "incomparable" in a sense, but if I wish to simply compare the average sugar content from items listed on the menu of each restaurant, what should I do? ![Menu Sugar Content Compare](http://i.stack.imgur.com/yiGv3.png)


title: Asymptotic Distribution of beta distribution parameter

text: I am trying to create a t-interval for the first shape parameter of a Beta distributed variable based on a random sample. $ f(x;\theta,\beta) = \dfrac{\Gamma(\theta+\beta)}{\Gamma(\theta)\Gamma(\beta)}x^{\theta-1}(x-1)^{\beta-1}$ Given that $E[X]=\dfrac{\theta}{\theta+\beta}$, I found $MME(\theta) = \dfrac{\beta\overline{X}}{1-\overline{X}}$ So I want to calculate a $(1-\alpha)*100\%$ t confidence interval by using $\dfrac{\beta\overline{X}}{1-\overline{X}} \pm t_{\alpha/2,n-1}*SE\left[\dfrac{\beta\overline{X}}{1-\overline{X}}\right]$ For the sake of my estimation, I will be fixing $\beta=4$. I have tried and tried to get this in generality as I feel like at it's heart all this quanity represents is a transformation of the sample mean (Delta method/CLT) but I can not for the life of me get any traction on it. After trying and failing, I ran some simulations and in a couple examples, setting $\theta=2$ , $\beta=4$, the ratio of $\dfrac{Var\left[\dfrac{\beta\overline{X}}{1-\overline{X}}\right]}{Var(\overline{X})} \approx 81.7$ Setting $\theta=3$ , $\beta=4$, ratio $\approx 151.5$ Those don't help too much, but do give me estimates to test against once I come to a theoretical result. If anybody can shed some light on my shortcomings I would appreciate it. Either way, thank you taking the time to look. R code I wrote for testing ratios of variances.               test <- function(alpha=1,beta=1){     ratio <- numeric()     for (j in 1:100){       means <- numeric()       alphas <- numeric()         for (i in 1:1000){             sample <- rbeta(100,alpha,beta)         means[i] <- mean(sample)         alphas[i] <- beta*mean(sample)/(1-mean(sample))           }         ratio[j] <- var(alphas)/var(means)       }     paste0("Ratio of Var: ",round(mean(ratio),3),". Ratio of SD: ",round(sqrt(mean(ratio)),3))     }          mapply(test,alpha=1:5,beta=4,SIMPLIFY=F)


title: Cumulative probability of 2 independant events

text: I’m trying to find the cumulative probability of 2 independent events (at least I think its cumulative probability). I will miss work on rainy days (event 1) where my car won’t start (event 2). The probability that my car won’t start on any given day is 20%. The probability of rain on any given day in a month is 10 in 30 (33.3%). The two events are independent of each other. The probability that my car won’t start on a rainy day is 0.2 * 0.33 = .067 In a 30 day month I would expect to miss work 0.67 * 30 = 2.0 days (rounded) But how can I find the probability of missing 1 days, 2 days, 3 days …all 30 days etc. during a 30 day month?


title: Gibbs sampling for Ising model

text: Homework question: Consider the 1-d Ising model. Let $x = (x_1,...x_d)$. $x_i$ is either -1 or +1 $\pi(x) \propto e^{\sum_{i=1}^{39}x_ix_{i+1}}$ Design a gibbs sampling algorithm to generate samples approximately from target distribution $\pi(x)$. **My attempt:** Randomly choose values (either -1 or 1) to fill vector $x = (x_1,...x_{40})$. So maybe $x = (-1, -1, 1, 1, 1, -1, 1, 1,...,1)$. So this is $x^0$. So now we need to move on and do the first iteration. We have to draw the 40 different x's for $x^1$ separately. So... Draw $x_1^1$ from $\pi(x_1 | x_2^0,...,x_{40}^0)$ Draw $x_2^1$ from $\pi(x_2 | x_1^1, x_3^0,...,x_{40}^0)$ Draw $x_3^1$ from $\pi(x_3 | x_1^1, x_2^1, x_4^0,...,x_{40}^0)$ Etc.. So the part that's tripping me up is how do we actually draw from the conditional distribution. How does $\pi(x) \propto e^{\sum_{i=1}^{39}x_ix_{i+1}}$ come into play? Maybe an example of one draw would clear things up.


title: what statistical test should i use for three groups with different sample sizes?

text: What statistical test should I use if I have three independent groups of interval data and would want to check whether their means are significantly different? I'd also note that the samples sizes of each group is different (N1= 50, N2= 10, N3=13). Would it be wiser to compare one with another and then do the same for the other two other pairs? (G1=G2; G2=G3; G1=G3)? I understand that one-way ANOVA and a Tukey's post-hoc would be the best option but does the difference in sample size pose a problem? Thank you!


title: Bayesian model averaging for variable selection in R

text: I am trying to use Bayesian model averaging for variable selection with a large number of variables. In R, the BMS package allows to apply the method, with the option of using MCMC sampler (Metropolis Hastings algorithm) when the number of covariates is large. Here is a sample code:                data(datafls)      fls1 = bms(datafls, burn = 50000, iter=100000, g = "BRIC", mprior = "uniform", nmodel = 2000, mcmc="bd", user.int=F)      result = coef(fls1)      However, if you run the same code twice, the results (i.e. the posterior probabilities) would be completely different. Does anyone know how to tune the code, so that the results are consistent for every run? Thank you,


title: How to create forecast data prediction interval bands

text: I have seasonal data from which I create forecasts. The steps I perform are: deseasonalizing the data, finding the linear regression for the deseasonalized points, predicting a few points from the linear regression and adding seasonality to the predicted values to get forecast data. My input is quite sinusoidal so all works well. The problem is that the more in the future you predict, the more prediction errors increase. I'd like to show that on a chart, but I am not sure how to calculate these errors. I was thinking something like prediction interval bands for forecast data (whatever they are called). These bands would increase the further you predict in the future. Here are some images that show what I'm trying to do: sample bands image1 sample bands image2 My **question** is what is the name for these bands? (then I can do a google search for it) I'd also appreciate the formulas needed for the band calculations. I'm guessing there is a standard deviation in there somewhere. I've looked at confidence interval, but that seems to be for the data already present, not for the forecast data.


title: anova, t-test, linear models: run on individual data points or per-person means?

text: _The format of my studies is testing 10 individuals in 5 conditions 100 times per condition (this is a behavioral test) - resulting in around 5000 raw data points._ I am currently using 3 statistical method "families" for my project. These are:   * t-tests   * ANOVAa   * Linear Models I am wondering for which of these methods it would be more beneficial to analyze ALL my data in raw form (this was my initial impulse for all methods, but apparently that inflates my significance at least for some) and for which I should only take the means of every individual as a data point. Can you help me out? Presently I am using t-tests and ANOVA for per- participant means and linear models on my raw data, though that is only because I have a couple of botched trials and paired t-tests and ANOVA don't deal well with missing data. But maybe it would be the right thing to do even for full data sets?


title: Outlier detection and smoothing for multi dimensional time series

text: From kinect depth images, I have collected the following time series that represent the features = 3D joint positions, quarternion angles, difference between hip joint and centroid of the right arm (composed of right shoulder, right elbow joint and right wrist). These are the feature sets at each time frame. But these are noisy and I need to eliminate those frames where there is no data, remove outlier and smooth the feature set. My data set contains frames where the values are zero since the recording does not start immediately. I have seen some methods like moving average filters and it is applied on each dimension for smoothing so the output contains the same data points. However, for outlier removal, the length of the dataset may get altered. As a result, I need to apply on the entire dataset instead of considering each time series individually. Q1 : What are the ways by which I can detect and remove outliers for multi dimensional timeseries so that I can remove bad frames which do not contribute from the good frames? Q2 : Are there any implementations or packages for achieving this?


title: Is AIC appropriate for model selection when the parameters are fitted by least-squares rather than MLE

text: I want to compare the fit of a linear model (M1) and nonlinear model (M2):   * M1: $y = b_0 + b_1x_1 + b_2x_2 + b_3x_1x_2 + \epsilon, \epsilon \sim N(0, \sigma^2)$   * M2: $y = b_0 + b_1x_1 + b_2x_2 + b_1 b_2x_1x_2 + \epsilon, \epsilon \sim N(0, \sigma^2)$ In particular I want to know whether M1 is significantly different from M2. To estimate the parameters, I am minimizing the least-squares errors rather than maximizing the likelihood through MLE procedures. In particular, I am using the R function nls() as follow:               # Creating a sample data set     n <- 50     x1 <- rnorm(n, 1:n, 0.5)     x2 <- rnorm(n, 1:n, 0.5)      b0 <- 1     b1 <- 0.5     b2 <- 0.2     y <- b0 + b1*x1 + b2*x2 + b1*b2*x1*x2 + rnorm(n,0,0.1)     # Actual model fit     M1 <- nls(y ~ b0 + b1*x1 + b2*x2 + b3*x1*x2, start=list(b0=1, b1=0.5, b2=0.5, b3=0.5))     M2 <- nls(y ~ b0 + b1*x1 + b2*x2 + b1*b2*x1*x2, start=list(b0=1, b1=0.5, b2=0.5))      I want to compare the models using a measure of relative fit such as AIC, which can be done in R as follow:               AIC(M1, M2)        df       AIC     M1  5 -88.47849     M2  4 -90.46491      Because $\Delta AIC \approx 2$ and the models differ by only one parameter, I would conclude that both of them fit the data similarly well. In addition, I want to know whether the parameter $b_3$ from M1 significantly add to the fit using a statistical test such as an F-test. This can be done in R as follow:               anova(M1, M2)     Analysis of Variance Table          Model 1: y ~ b0 + b1 * x1 + b2 * x2 + b3 * x1 * x2     Model 2: y ~ b0 + b1 * x1 + b2 * x2 + b1 * b2 * x1 * x2       Res.Df Res.Sum Sq Df      Sum Sq F value Pr(>F)     1     46    0.40843                                   2     47    0.40855 -1 -0.00011097  0.0125 0.9115      My general question is:   * Are these analyses appropriate?  More specifically:   * Can I use AIC to compare least-squares fitted models? From a few posts such as this one it looks like AIC should be appropriate. However, I've seen posts such as this one that indicates that using AIC on non-MLE fitted models might be a problem. I understand that least-squares is equivalent to MLE if the error is normally distributed, but is this true even for non-linear models?   * Can I use a F-test to test whether $b_3$ is significantly different from $b_1 b_2$?  I know such F-test makes sense if the model are nested, but I'm unsure whether it is appropriate in this case.


title: Probability as a dependent variable in a time-series regression

text: Are there any issues to run a Newey-West time-series regression on a dependent variable that is a probability? What are the biases that I am facing? I can't find anything online that can help me out with this issue.


title: Correct calculation of effect size after Friedman test in a skewed population

text: I have some data that compares behaviour in four conditions (repeated measures)(I have several different behavioural variables and each behaviour gets analysed separately, and the behaviours are currently in % of time format). The behavioural variables are very non-normal, even after transformation. So I have run Friedman tests for the omnibus tests (friedman.test in R). For posthocs I have used friedman from the Agricolae package, but I would also like to report effect sizes. I am not sure which is the correct calculation to use because the data is not only non-normal, but also very positively skewed. I have tried the rFromWilcox function from Field but I am not sure it is the right one to use. Also I am getting r results such as -0.564NULL. I am quite new to doing my own analysis, but keen to learn, so please don't be offended if I haven't phrased something quite right. @Michael: I am trying to determine if certain behaviours occur more or less often under the different conditions, and where those differences are (nothing too deep). Two similar but different treatments are applied between condition one and two, and between three and four. I am predicting that there will be certain changes (some behaviours will increase and some will decrease) between all conditions but particularly between one and two, and three and four. My posthocs test between all pair combinations and I have used a Bonferroni correction. I was planning on reporting the summary stats, for each condition for each behaviour, the test statistic and p value, and posthoc results where Friedman is significant, and effect sizes on those posthocs.


title: Motifs for prediction

text: I'm trying to find a method for using motifs for prediction. Motifs defined as similar (sub)sequences can be found using approximate or exact methods, see Exact Discovery of Time Series Motifs. My problem is: once you have motifs of labelled time series how could the motifs be used for prediction?


title: Is principal components analysis valid if the distribution(s) are Zipf like? What would be similar to PCA but suited to non gaussian data?

text: I'm analyzing people based on their twitter stream. We are using a 'word bag' model of users, which basically amounts to counting how often each word appears in a persons twitter stream (and then using that as a proxy for a more normalized 'probability they will use a given word' in a particular length of text). Due to constraints further down the pipeline, we cannot retain full data on usage of all words for all users, so we are trying to find the most 'symbolically efficient' words to retain in our analysis. That is, we're trying to retain a subset of dimensions, which, knowing their values would allow a hypothetical seer to most accurately model the probabilities of all words (including any we left out of the analysis). So a principal components analysis (PCA) type approach seems an appropriate first step. (happily ignoring for now the fact that PCA would also 'rotate' us into dimensions that don't correspond to any particular word). But I am reading that "Zipf distributions .. characterize the use of words in a natural language (like English) " and as far as I know, PCA analysis makes various assumptions about the data being normally distributed. So, I'm wondering whether the fundamental assumptions of the PCA analysis will be sufficiently far 'off' from reality to be a ral problem. That is, does PCA rely on the data being 'close to' Gaussian Normal for it to work at all well? If this is a problem as I suspect, are there any other recommendations? That is, some other approach worth investigating that is 'equivalent' to PCA in some way but more appropriate for Zipf or power law distributed data? Note that I am a programmer, not a statistician, so apologies if I messed up my terminology in the above. (Corrections of course welcomed!)


title: The probability for two people to provide identical answers on survey questions

text: Here is the problem: A survey contains 7 binary questions (Yes/No responses). If two people are answering the survey, what is the probability for their answers on 4 or more of the questions to match? In other words, if we have four or more matching answers, we can consider the overall survey response to be similar for both people.


title: PDF of sum of one continuous logistic and one Bernoulli random variables

text: I need some help people. I have a logistic r.v. $V\sim \Lambda (0, \frac {\pi^2}{3})$ and a Bernoulli $Z \sim B(p_z)$, independent. I have also $c\in R$, $\tilde Z = cZ$, and I have to find the density of $$U = V + \tilde Z$$ I know that one can use the convolution method by using the dirac delta function to write the **pmf** of $\tilde Z$ as a weighted sum, but I have trouble with the nuts and bolts, I couldn't find something clear in the books, while some related posts I found here and in math.SE didn't do the trick for me. $1)$ Is $$P(\tilde Z=\tilde z) = p_z\delta (\tilde z+c) + (1-p_z)\delta(\tilde z)$$ the correct way to write the **pmf** of $\tilde Z$? and, $2)$ in the convolution integral should I make the substitution $u =v+\tilde z \Rightarrow v=u- \tilde z$ and integrate w.r.t to $\tilde z$, or $u =v+\tilde z \Rightarrow \tilde z=u-v$ and integrate w.r.t to $v$? I have a feeling this is trivial, but this time I did go blind.


title: Construct matrix of stacked variables in VAR regression

text: I am trying to NOT use packages for the estimation of models in order to have a deeper understanding of how things work. Currently, I am trying to estimate a VAR(1) (vector autoregression of first order) and my question is about how to construct the matrix of dependent and independent variables. Theoretically, my model is: ![model](http://latex.codecogs.com/gif.latex?%5Cbinom%7Bx_t%7D%7By_t%7D%20%3D%20%5Cbinom%7Bc_1%7D%7Bc_2%7D%20+%20%5Cbegin%7Bpmatrix%7D%20a_%7B11%7D%20%26%20a_%7B12%7D%5C%5C%20a_%7B21%7D%20%26%20a_%7B22%7D%20%5Cend%7Bpmatrix%7D%5Cbinom%7Bx_%7Bt-1%7D%7D%7Bx_%7Bt-2%7D%7D%20+%20%5Cbinom%7Bu_%7B1%2Ct%7D%7D%7Bu_%7B2%2Ct%7D%7D) As for the matrix of dependent variables, it should be easy enough. I just need to stack `xt` and `zt`. However, I don't know how to construct the matrix of independent variables. I read that VAR is very similar to SUR (Seemingly Unrelated Regressions) so I think that the matrix of independent variables should be a block diagonal matrix, but what should it look like? Especially if I want to add a constant and a time trend to the model? Thank you.


title: Can I change the proposal distribution in random-walk MH MCMC without affecting Markovianity?

text: Random walk Metropolis-Hasitings with symmetric proposal $q(x|y)= g(|y-x|)$ has the property that the acceptance probability $$P(accept\ y) = \min\\{1, f(y)/f(x)\\}$$ does not depend on proposal $g(\cdot)$. Does that mean that I can change the $g(\cdot)$ as a function of previous performance of the chain, without affecting the markovianity of the chain? Of particular interest to me is the adjustment of the scaling of Normal proposal as a function of acceptance rate. Would also greatly appreciate if someone can point out to the adaptation algorithms used in practice for this type of problem. Many thanks. [edit: Starting with the references given by robertsy and wok I found the following references on MH adaptive algorithms: Andrieu, Christophe, and Éric Moulines. 2006.   On the Ergodicity Properties of Some Adaptive MCMC Algorithms. The Annals of Applied Probability 16, no. 3: 1462-1505. http://www.jstor.org/stable/25442804. Andrieu, Christophe, and Johannes Thoms.   2008\. A tutorial on adaptive MCMC. Statistics and Computing 18, no. 4 (12): 343-373. doi:10.1007/s11222-008-9110-y. http://www.springerlink.com/content/979087678366r78v/. Atchadé, Y., G. Fort, E. Moulines, and P. Priouret. 2009.   Adaptive Markov Chain Monte Carlo: Theory and Methods. Preprint. Atchadé, Yves. 2010.   Limit theorems for some adaptive MCMC algorithms with subgeometric kernels. Bernoulli 16, no. 1 (February): 116-154. doi:10.3150/09-BEJ199. http://projecteuclid.org/DPubS?verb=Display&version=1.0&service=UI&handle=euclid.bj/1265984706&page=record. Cappé, O., S. J Godsill, and E. Moulines. 2007.   An overview of existing methods and recent advances in sequential Monte Carlo. Proceedings of the IEEE 95, no. 5: 899-924. Giordani, Paolo. 2010.   Adaptive Independent Metropolis–Hastings by Fast Estimation of Mixtures of Normals. Journal of Computational and Graphical Statistics 19, no. 2 (6): 243-259. doi:10.1198/jcgs.2009.07174. http://pubs.amstat.org/doi/abs/10.1198/jcgs.2009.07174. Latuszynski, Krzysztof, Gareth O Roberts, and Jeffrey S Rosenthal. 2011.   Adaptive Gibbs samplers and related MCMC methods. 1101.5838 (January 30). http://arxiv.org/abs/1101.5838. Pasarica, C., and A. Gelman. 2009.   Adaptively scaling the Metropolis algorithm using expected squared jumped distance. Statistica Sinica. Roberts, Gareth O. 2009.   Examples of Adaptive MCMC. Journal of Computational and Graphical Statistics 18, no. 2 (6): 349-367. doi:10.1198/jcgs.2009.06134. http://pubs.amstat.org/doi/abs/10.1198/jcgs.2009.06134. ]


title: From an email address to a quasi-random number

text: ### My Aim: I'd like to a have a function that takes an email address and outputs a quasi- random number of 1, 2, 3, or 4. ### A little detail: By quasi-random number I mean that given a typical population of email addresses, the probabilities of getting a value of 1, 2, 3, or 4 are roughly equal, and that obvious systematic properties of the email address such as the domain name do not affect the probability of getting a value of 1, 2, 3, or 4. ### A little background: I have an online experiment written in inquisit where participants log in on two occasions. I want to randomly assign participants to one of four groups. While this is easy to do for one session (I can just use a random number generator), I need some way of remembering the allocation across sessions. Thus, I thought that I could extract a quasi-random group allocation from the participant email. I'm also limited in the set of functions that I have at my disposal (see here for full list). The string functions are: tolower toupper capitalize concat search replaceall contains startswith endswith substring trim trimright trimleft length format evaluate ### Initial Thoughts: I thought about trying to extract a set of features of the email address that returned a value of 1, 2, 3, or 4 with roughly equal probabilities. Then, I could sum these properties and get the mod 4 plus 1 of that. Thus, assuming something like the central limit theorem, I might get close. Possible features that came to my mind:   * length of string   * position of first "a", "b", etc.


title: Back-testing Value at Risk with a WML investment strategy

text: I'm currently taking a course in Financial Econometrics and there is a question in the lecture notes regarding back-testing of VaR which I'm have difficulty with. First of all the procedure for back-testing of VaR using a rolling scheme is described as follows: assume we have returns data $r_{t}, t=1,...,T$ for some large enough $T$. We take the first 60 observations and calculate their sample mean $\mu_{1-60}$ and variance $\sigma_{1-60}^{2}$ and then calculate: $$\mbox{VaR}_{1-60}\left(\alpha\right)=-\mu_{1-60}+\Phi^{-1}\left(\alpha\right)\cdot\sigma_{1-60}$$ Where $\alpha$ is our confidence parameter (usually $\alpha=0.05$) and $\Phi^{-1}$ is the inverse standard normal CDF. If $r_{61}<-\mbox{VaR}_{1-60}\left(\alpha\right)$ we mark $1$ and otherwise $0$. We proceed doing the same thing for observations $2...62$ compared to $63$ and so on and in the end we count the number of times in which the result was $1$. while under the null hypothesis we expect that the proportion of times we got $1$ to be at most $\alpha$. In particular it's important to notice we had an assumption here that the returns are all normally distributed, $r_{t}\sim N\left(\mu_{t},\sigma_{t}^{2}\right)$. Now for the actual question, suppose $r_{t}, t=1,...,T$ are the returns using a WML investment strategy (momentum trading strategy). Is the procedure described above suitable for testing the VaR in such a case and if not why? Usually the phrasing of the question would hint towards there being some sort of problem with using the procedure in such a case but I don't see why that would be the case... Help would very much be appreciated! **EDIT** : WML stands for Winners minus Losers. It's an investment strategy in which you buy past winners and short-sell past losers based on the returns in an evaluation period. For example if you have 4 stocks in your portfolio then: $WML=\frac{1}{2}\left(\mbox{SmallHigh + BigHigh}\right)-\frac{1}{2}\left(\mbox{SmallLow}+\mbox{BigLow}\right)$ So in general for $n$ stocks you take the average of the $n/2$ high returns and subtract the average of the $n/2$ low returns.


title: Minimizing the sum of squares of autocorrelation function of residuals instead of sum of squares of residuals

text: I am trying to fit my multi-exponential model to some experimental data and I am using a simulated annealing algorithm. My objective function has so far been the sum of squares of the residuals: $\sum_{i=1}^n (y_i - \hat{y_i})^2 = \sum_{i=1}^n e_i^2 \:\:$ ...(1) However, it came to my attention that maybe I could set my objective function to be the sum of squares of the autocorrelation function of the residuals, with lag $j$: $\sum_{j=1}^m R^2_{j} = \sum_{j=1}^m \left( \sum_{i=1}^n e_i\,\overline{e}_{i-j} \right)^2 \:\:$ ...(2) Does that make sense? It seems like forcing the autocorrelation of the residual terms to be small may lead to wrong parameter estimates. I haven't rigorously proven this, but I suspect that minimizing (2) does not mean that (1) is minimized. Any ideas?


title: How to estimate confidence interval of the sample mean of a non-stationary time series?

text: I need an estimate of the average value $\bar{X}$ of the data within the shaded area of this graph, and its corresponding confidence interval. ![ideal case](http://i.stack.imgur.com/4nIRV.png) In the ideal case, the data would reflect some minor oscillations. I'm tempted to ignore this, and use the $CLT$ so the standard error would be $SE=\frac{S}{\sqrt{N}}$, where $S$ is the sample standard deviation, and $N$ is the total number of data points in the range. But I'm uncomfortable with this approach because I know this is not completely correct. But the main problem is that sometimes the data is not stationary for the whole range, and can show multiple stationary segments, as shown below. ![worst case](http://i.stack.imgur.com/SMHEV.png) These issues invalidate any assumption of normality. I'm struggling to identify a proper method to calculate a single value for the mean and confidence interval for the whole data range. Any ideas?


title: How Many Samples Needed for Classification

text: I am currently trying to perform classification on a set of data. I have 19 category_1 observations and 15 category_2 observations. What is the best way to sample from this dataset to train and test my classifier, or are there simply not enough observations?


title: Does SurveyMonkey ignore the fact that you get a non-random sample?

text: SurveyMonkey has steps and a chart for you to figure out what sample size you need for a given margin of error or confidence interval, based on your population size. SurveyMonkey sample size Does this chart simply ignore the fact that you will not get a random sample, since you only get the people who bother to respond to the survey? I'm getting warned as I type this that the question appears subjective so maybe I'm not asking it correctly. It's not really about SurveyMonkey but is a more general question -- can you actually calculate confidence intervals from voluntary response data using some advanced techniques that I don't know about? In exit polls or national surveys, obviously they must deal with this problem. My education didn't cover survey sampling techniques in depth but I assume it involves collecting demographic data and using that to know how representative of a sample you have. But aside from that, for a simple online survey, are they just assuming that the people who bother to respond are a random sample of the population?


title: Possibly insufficient data on a shuffling problem

text: _A company controls its employees' access to a restricted area with magnetic cards and passwords. Each of the $k$ employees received an anonymous card and chose an $n$-digit numerical password to associate to their card. Suppose the cards were shuffled and randomly redistributed, and let $A$ be an employee._ _a) What is the probability that $A$ enters the restricted area using his password and the card he got after the shuffle?_ _b) If $A$ was able to enter with his password, what is the probability that he received the same card before and after the shuffle?_ * * * a) One easily notes that the probability that an employee picked a given password is $p=10^{-n}$. Now let $L$ be how many of the $k$ cards are associated to $A$'s password, that is, how many employees (including $A$) chose $A$'s password. Then the probability is $L/k$. * * * b) Let $E_1$ be the event "$A$ entered the restricted area" and $E_2$ be "$A$ got his own card after the shuffle". $E_1$ is just the event asked on _a_ , and since $E_2\subset E_1$, this probability is $$\mathbb P(E_2\,|\,E_1)=\frac{\mathbb P(E_2\cap E_1)}{\mathbb P(E_1)}=\frac{\mathbb P(E_2)}{\mathbb P(E_1)} =\frac{1/k}{L/k} =\frac1L\quad.$$ * * * Clearly, all there is left is to calculate $L$. How do I get it deterministically? All I know how it is distributed: $$\begin{align} \mathbb P(L=l\,|\,L\ge1) &=\frac{\mathbb P(L=l\wedge L\ge1)}{\mathbb P(L\ge1)} =\\\ &=\frac{\displaystyle{\binom kl}p^l(1-p)^{k-l}I_{\\{1,...,k\\}}(l)} {\mathbb P(L\neq0)} =\\\\\,\\\\\,\\\ &=\frac{\ldots} {1-\displaystyle{\binom k0}p^0(1-p)^{k-0}} =\\\\\,\\\\\,\\\ &=\frac{\ldots}{1-(1-p)^k}\quad, \end{align}$$ and so $$\mathbb P(L=l)=\frac1{1-(1-10^{-n})^k}\binom kl10^{-nl}(1-10^{-n})^{k-l}I_{\\{1,...,k\\}}(l)\quad.$$


title: Logistic regression and success probability: Bayesian statistics

text: I have a Bayesian statistics homework question that I'm not sure of: Running this in `R`, I got completely flat lines at `0` and `1` so I was thinking this meant that the success probability does not change with this particular prior. Can anyone shed some more light on this problem?               mu0    = 0      sig2_0 = 100      n      = 10000      b0     = rnorm(n, mean=mu0, sd=sqrt(sig2_0))      x      = seq(from=0, to=10, length.out=100)      res    = matrix(0, length(x), 3)      for(i in 1:length(x)){        z       = b0        z1      = exp(z) / (1+exp(z))        res[i,] = quantile(z1, c(0.025, 0.5, 0.975))      }           plot( x, res[,2], type='l', ylim=range(res), ylab='Success probability', xlab='Covariate')      lines(x, res[,1], lty=3)      lines(x, res[,3], lty=3)


title: Generative Models: Bayesian statistics hw help

text: Can anyone shed some light on this problem. Not looking for you to write out the answer for me, just some helpful hints that will hopefully lead me in the right direction. Here is the question: > We are interested in determining the proportion of occupied sites. The > problem here is that detection is not perfect. While detection of our focus > animal at a particular site implies that the site is occupied, the lack of > detection does not imply that the animal is absent. Similar to the example > in class, write down the generative model for this presence/absence data. Just to clarify I have attempted this problem and got the following solution: I used P (present), S (signs of animal), and D (animal detected) P=1 (animal present) P=0 (animal absent) S=1 (signs of animal) S=0 (no signs of animal) D=1 (animal detected) D=0 (animal not detected) So I worked this out using rolling a dice, so: P=1 (1,2), S=1 (1-5), D=1 (1-3), D=0 (6)----> animal present P=0 (3-6), S=1 (1), D=0, S=) (2-6)-----> animal absent P~ Bernoulli (1/3) S=1/P=1~ Bernoulli= (5/6) D= 1/S=1, P= 1~ Bernoulli= (1/2) S= 1/P=0~Bernoulli= (1/6) D= 0/S= 1,P=0~ Bernoulli=(0)


title: Homoscedasticity violation: repeated-measures ANOVA

text: My study is a clinical trial with _2 treatment groups_ measured over _4 time points_...with the outcomes assessed being _blood levels_. So my independent variables are "treatment" (2 groups) and "time" (4 time points - with time being a within-subjects factor)...and my dependent variable is "blood levels". So...I'm doing a **repeated measures ANOVA (mixed model)** in SPSS, but my assumption of homoscedasticity is not met. My data meets the assumption of normality and sphericity...but not of homogeneity of variances/homoscedasticity. I'm wondering what my best options are...given this limitation. I've tried transformations and they haven't really helped. I'd really appreciate any suggestions/help. **EDIT** : I understand the possibility of options in my case such as: 1\. repeated paired t-tests with a bonferroni correction (at the cost of loss of power) 2\. non-parametric equivalents (i.e. Friedman test) to the rmANOVA ...but which is a more 'robust' test in lieu of the scientific design?


title: Bootstrap confidence interval for a biased estimator

text: I have a statistic $\hat{\theta}$ for the parameter $\theta$. Which may be biased. Assume $\mathbb{E}[(\hat{\theta}-\theta)^2]=\textit{ecm}^2$ and $\mathbb{V}[\hat{\theta}]=\sigma^2$ are known, but the distribution of $\hat{\theta}$ is unknown. I want to find a confidence interval for $\theta$ using bootstrap. I found on the internet that usually one can use bootstrap to find $k_\alpha$ such that $P( \frac{(\hat{\theta}-\theta)^2}{\sigma^2} \leq k_\alpha^2) = 1-\alpha$. I understand this is reazonable because it's analog to normal confidence intervals and in that case it's a pivotal statistic. My question is, Should I use $\textit{ecm}^2$ in the denominator instead of $\sigma^2$? If the estimator is unbiased both options are the same because ($\sigma^2=\textit{ecm}^2$) But my interest is for biased estimators. Also, What happens if instead of knowing $\sigma^2$ and $\textit{ecm}^2$ I have their bootstrap estimations? Any hint or reference on the subject would be great. Thanks


title: The state-of-the-art methods for speech recognition?

text: Recently I noticed that Google and Apple have really high quality speech- recognition services. I was wondering about the state-of-the-art methods and techniques they are/might be using to achieve such quality. I already know that Hidden Markov Models can be used for speech recognition, but I was wondering if it's the used technique now. Because as I know HMM are old. It would be great if you can link to some important recent papers that could include/descripe the used new methods.


title: Selecting gene list for subsequent analysis in problematic microarray experiement

text: I have an experimental design problem and I'm not sure which would the best way to proceed. We have a micro-array experiment in which we compare gene expression profiles between 2 groups of patients. Each group has 12 independent biological replicates. Lets call these groups CT and EXP. The goal of this experiment is to identify differentially expressed genes between CT and EXP. Two sets of hybridizations where done: first and second. First hybridization had big issue: CT samples were processed at one time-point and EXP couple of months later, creating a perfectly confounded batch effect. Second hybridization did not have this problem - all of the samples were processed at the same time. I've analyzed both sets separately and came up with 2 lists of differentially expressed genes. My question is, should I discard the results of the first hybridization completely and stick with the second one since it didn't have an obvious design issue? Or should I use both sets of hybridizations and select only genes which show up as differentially expressed in every case (there are 2900 of such genes)? Appreciate your help!


title: Transformation of random variable - lognormal distribution

text: I'm answering questions from a book that I own and I'm left scratching my head because one of my answers does not match the one provided in the book. Either the book contains a typo or I've simply done something wrong. It's a good book and it's by some quality authors, so I assume it's the latter case and hence the reason why I'm asking this question. In cases like this, I normally check if the book has an errata, but it's absent from the book's website. From the outset, this is not coursework; I'm not privileged enough to be on a course! **Question in Book** Consider the following model of asset prices, $p_{t}$, that is commonly adopted in the financial econometrics literature $$ \ln p_{t} - \ln p_{t-1} = \alpha + u_{t}~, ~~~~u_{t} \sim iid ~ N(0,\sigma^2) ~, $$ in which $\theta = \\{\alpha, \sigma^2\\}$ are unknown parameters. Use the transformation of variable technique to show that the conditional distribution of $p_{t}$ is the log-normal distribution $$ f(p_{t}|p_{t-1};\theta) = \frac{1}{\sqrt{2\pi\sigma^2}p_{t}}\exp \left[-\frac{\ln p_{t}-\ln p_{t-1}-\alpha}{2\sigma^2}\right]~~. $$ **My Answer** \begin{align} f(u_{t}) &= \frac{1}{\sigma\sqrt{2\pi}}\exp \left[-\frac{(\frac{u_{t}}{\sigma})^2}{2}\right]~~~~(\text{by definition of } N(0,\sigma^2))\\\ &= \frac{1}{\sqrt{2\pi\sigma^2}}\exp \left[-\frac{u_{t}^2}{2\sigma^2}\right]~~~~~(\text{by rules of exponents})\\\ \end{align} Applying the change of variable technique, $$ f(p_{t}|p_{t-1};\theta) = f(u_{t})\left|\frac{\partial u_{t}}{\partial p_{t}}\right| ~, $$ and recognizing that $u_{t} = \ln p_{t} - \ln p_{t-1} - \alpha$, yields \begin{align} f(p_{t}|p_{t-1};\theta) &= \frac{1}{\sqrt{2\pi\sigma^2}}\exp \left[-\frac{(\ln p_{t} - \ln p_{t-1} - \alpha)^2}{2\sigma^2}\right]\left|\frac{1}{p_{t}}\right| \\\ &= \frac{1}{\sqrt{2\pi\sigma^2}p_{t}}\exp \left[-\frac{(\ln p_{t} - \ln p_{t-1} - \alpha)^2}{2\sigma^2}\right] ,\\\ \end{align} where we note that the determinant of 1x1 matrix is its sole entry. **Problem** Comparing terms, one can see that my answer differs from the book's answer with regard to the numerator in the exponent - my answer contains a squared term while the answer in the book does not. Please point out my egregious error and give my head some rest. Cheers!


title: Is a model nested within itself before collapsing categorical variables?

text: If I have a model with a categorical variable $X_1=\\{0,1,2,3\\}$ and a continuous variable $X_2$, and I have a regression model that includes an interaction between $X_1$ and $X_2$, then I decide I want to collapse $X_1$ into an indicator variable where $X_{1\text{_new}}=0$ if $X=0$ and $X_{1\text{_new}}=1$ if $X\ge1$, is my new model with an interaction term between $X_{1\text{_new}}$ and $X_2$ nested with the original model?


title: Comparing maximum likelihood estimation (MLE) and Bayes' Theorem

text: In Bayesian theorem, $$p(y|x) = \frac{p(x|y)p(y)}{p(x)}$$, and from the book I'm reading, $p(x|y)$ is called the _likelihood_ , but I assume it's just the _conditional probability_ of $x$ given $y$, right? The _maximum likelihood estimation_ tries to maximize $p(x|y)$, right? If so, I'm badly confused, because $x,y$ are both random variables, right? To maximize $p(x|y)$ is just to find out _the_ $\hat y$? One more problem, if these 2 random variables are independent, then $p(x|y)$ is just $p(x)$, right? Then maximizing $p(x|y)$ is to maximize $p(x)$. Or maybe, $p(x|y)$ is a function of some parameters $\theta$, that is $p(x|y; \theta)$, and MLE tries to find the $\theta$ which can maximize $p(x|y)$? Or even that $y$ is actually the parameters of the model, not random variable, maximizing the likelihood is to find the $\hat y$? **UPDATE** I'm a novice in machine learning, and this problem is a confusion from the stuff I read from a machine learning tutorial. Here it is, given an observed dataset $\\{x_1,x_2,...,x_n\\}$, the target values are $\\{y_1,y_2,...,y_n\\}$, and I try to fit a model over this dataset, so I assume that, given $x$, $y$ has a form of distribution named $W$ parameterized by $\theta$, that is $p(y|x; \theta)$, and I assume this is the _posterior probability_ , right? Now to estimate the value of $\theta$, I use MLE. OK, here comes my problem, I think the likelihood is $p(x|y;\theta)$, right? Maximizing the likelihood means I should pick the right $\theta$ and $y$? If my understanding of likelihood is wrong, please show me the right way.


title: Strategies for introducing advanced statistics to various audiences

text: I work mainly with non-statisticians in fields such as medicine, social sciences and education. Whether I am consulting with graduate students, helping researchers with articles or reviewing articles for journals, I often have the problem that someone (client, author, dissertation committee, journal editor) wants to use some relatively well-known technique when it is either entirely inappropriate or when better but lesser-known methods exist. Often, I will explain the alternative technique but then be told "everybody does it the other way". I'd be interested in how others deal with this sort of difficulty. ADDITIONS @MichaelChernick suggested I could share some stories, so I will Currently I am working with one person who is duplicating a previous paper and adding one independent variable to see if it helps. The previous paper is, frankly, terrible. It treats dependent data as if they were independent; it is tremendously overfit and there are other problems too. Yet he (my client) submitted an earlier version as a dissertation and not only got his degree but was widely praised for the research. _Many_ times I have tried to convince people not to dichotomize variables. This comes up very often in medicine. I patiently point out that dicohotomizing (say) birthweight into low and normal (usually at 2,500 g) means treating a 2,499 g baby as just like a 1,400 g one; but treating the 2,501 gram baby quite differently. The clinician agrees with me that this is silly. Then says to do it that way. I had a graduate student client long ago whose committee _insisted_ on a cluster analysis. The student did not understand the method, the method did not answer useful questions, but that's what the committee wanted, so that's what they got. The entire field of statistical graphics is one where, for many, "this is how grandpa did it" is enough. Then there are people who seem to just push buttons. I remember one presentation (not by someone I helped!) who had taken an entire questionnaire and factor analyzed it. One of the variables she included was ID number! Oy.


title: Matlab RandomForest prediction error calculation

text: I am using Random Forests in Matlab for regression. After educating my model on train data, I want to get MSE on test data not used in training. I do that two ways:   1. call `predict` and directly calculate MSE using predicted and actual values   2. call `error` and use built in TreeBagger functionality to do the same task. In first case I get 10 times bigger result. Why? The only explanation I have, is that built in function somehow discounts outliers in prediction, but I am not sure how exactly it is done. Can somebody, please, explain all this to me.


title: Combination of multinomial logit model and logit model

text: I want to analyze the determinants of credit constraints of a firm. I have information for both formal and informal credit. I have 6 categories of credit-constraint statuses of a firm for formal credit, while I only have 2 categories for informal credit. As credit constraint status is determined simultaneously, and formal and informal credit market are interrelated, I want to estimate the determinants of credit constraints for formal and informal credit market at the same time. Bivariate logit/probit model may be too simplistic. Is there extension of bivariate logit model or some kind of multivariate multinomial logit model in Stata? Thank you very much in advance.


title: Regressing on many outcomes in a spatial model

text: I'm measuring many microbe populations in a 'river model' ecosystem by sampling at two points. I want to test the hypothesis that key species upstream are the best predictors for any given population at a sample site downstream (and I want to describe the best predictors). A few questions. The first is a general question about regression with many outcomes: the population level of any of hundreds of species could be important. So do I need to be concerned with using linear regression on all the upstream species to predict each downstream, species response vector? I'm thinking the problem here is false discovery due to multiple comparisons. 2) With many features and potential predictors, what is the best regression method to use? I expect that closely related species will be the best predictors of downstream populations... so CARTs to intuitively combine features? PCA/feature selection to reduce dimensions to related components? Or just let a linear model (glm?) algo 'factor out' unimportant features? 3) The spatial component is pretty simple. I expect downstream species abundance to be predictable from species upstream. So the interesting thing here would probably be species abundance that we can't explain. That being said, is it necessary to deal with potential correlation between the two sites? My preference is to ignore this component for now. Thanks everyone!


title: Could anyone identify this book?

text: I am trying to find a book I once partially read without much luck. The book's main topic, if I remember right, is stochastic processes. It was, I think, a not too heavy blue/light blue book. The word "King" somehow sounds related, but I am not sure if it was the author, or maybe written at King's College. I also think it was written originally in the mid 1980s, but I am not sure. From what I remember, the first or second chapter gave a very high level overview of stochastic processes, and showed how regular distributions, such as Bernoulli, generalize to stochastic processes. It also shows that many processes are related to the Poisson process. The way it was explained is something like (not in technical terms), a stochastic process is such that it is defined in a larger space, with each finite subspace originating in a regular distribution.


title: Route to understanding and implementing statistical timing models

text: The research paper Forecasting NBA Player Performance by Douglas Hwang describes a method for predicting player performance by utilizing a statistical timing model and a 'Weibull-Gamma' distribution. I'm approaching this from a rudimentary understanding of descriptive statistics and Bayesian inference, and would like to get some background to better appreciate, and potentially implement a similar approach to predicting individual performance (outside of the sports arena).   * What branch of statistics encompasses this sort of analysis?   * What is the usual route students/practitioners follow in learning how to design and implement these analysis?   * Any suggestions on books or resources (as applied as possible)? Or open-source libraries that provide implementations of this sort of thing?


title: Optimal parameter $\alpha$ for exponential smoothing using least squares

text: I was reading about exponential smoothing on Wikipedia. I saw a sentence saying that I could find the optimal parameter $\alpha$ by using the method of least squares. How do I go about doing that? Are there other ways I can find the right $\alpha$? **Edit:** Also, if possible, I'd like a way to handle missing values.


title: MLE vs MAP vs conditional MLE with regards to logistic regression

text: We have some set of iid RV's: $(X_i, Y_i), \; i=1,\ldots n$. We believe each to be distributed as $P(X_i, Y_i | \theta)$. So that $$ P(X,Y | \theta) = \prod_i P_i(X_i, Y_i | \theta) $$ Now using Baye's rule: $$ P(\theta|X,Y) = \frac{P(X,Y|\theta)P(\theta)}{P(X,Y)} = \frac{P(\theta)\prod_i P_i(X_i, Y_i | \theta)}{P(X,Y)} $$ As I understand it, MLE, MAP, and conditional MLE all attempt to find the best parameters, $\theta$, given the data by maximizing the left hand side by maximizing a subset of terms on the right. For MLE, we maximize the likelihood term, $\prod_i P_i(X_i, Y_i | \theta)$. For MAP, we maximize all of the numerator, $P(\theta)\prod_i P_i(X_i, Y_i | \theta)$. For conditional MLE (as in logistic regression), we have $$ \frac{P(\theta)\prod_i P_i(X_i, Y_i | \theta)}{P(X,Y)} = \frac{P(\theta) \left( \prod_i P_i(Y_i | X_i, \theta) \right) \left( \prod_i P(X_i|\theta) \right) }{P(X,Y)} $$ Conditional MLE maximizes only the $\prod_i P_i(Y_i | X_i, \theta)$ term. Is this correct? I have seen regularized logistic regression amounting to maximizing the prior, $P(\theta)$. Would modeling the third distribution for conditional MLE, $\prod_i P(X_i|\theta)$, be something different entirely? I understand that logistic regression is discriminative model. Is this a result of this? Would modeling $P(X_i|\theta)$ then give us a generative model? Thanks for any pointers.


title: Who uses R with multicore, SNOW or CUDA package for resource intense computing?

text: Who of you in this forum uses ">R with the multicore, snow packages, or CUDA, so for advanced calculations that need more power than a workstation CPU? On which hardware do you compute these scripts? At home/work or do you have data center access somewhere? The background of these questions is the following: I am currently writing my M.Sc. thesis about R and High Performance Computing and need a strong knowledge about who actually uses R. I read that R had 1 million users in 2008, but that's more or less the only user statistics I could find on this topic - so I hope for your answers! Sincerely Heinrich


title: Machine Learning conferences?

text: What are the most significant annual Machine Learning conferences? Rules:   1. One conference per answer   2. Include a link to the conference


Researchers conducted two different types of test on a large group of people. After that, the researchers subjected participants to situations like Z and noted their response R. Those with high scores on test T2 had a good response. Likewise, those with low scores on test T2 had a bad response. Therefore, researchers conclude that the score of test T2 affects the response R. However, the underlying truth is:   1. Score of test T1 affects response R.   2. Fact 2: Score of test T1 affects score of test T2.   3. Fact 3: Score of test T2 does not affect the response R. Here, the researchers wrongly assumed that the score of test T2 affects the response R. However, it is not true because the real culprit is test T1. What error have the researchers made, and what is the correct statistical procedure for this scenario? I.e., **What error is it to assume a test causes a response when it in fact only correlates with another test that causes a response?**


title: Calculating R-squared (coefficient of determination) with centered vs. un-centered sums of squares

text: I am looking at some simple regression models using both R and the `statsmodels` package of Python. I've found that, when computing the coefficient of determination, `statmodels` uses the following formula for $R^2$: $$ R^2 = 1 - \frac{SSR}{TSS}~~~~~~(\text{centered}) $$ where $SSR$ is the sum of squared residuals, and $TSS$ is the total sum of squares of the model. ("Centered" means that the mean has been removed from the series.) However, the same calculation in R yields a different result for $R^2$. The reason is that R seems to be calculating $R^2$ as: $$ R^2 = 1 - \frac{SSR}{TSS}~~~(\text{uncentered}) $$ So, what gives? Presumably there's some reason to prefer one over the other in certain situations. I haven't been able to find any information online about the cases where one of the above formulae should be preferred. Can someone please explain why one is better than the other?


title: When to use an offset in a Poisson regression?

text: Does anybody know why offset in a Poisson regression is used? What do you achieve by this?


title: Forecasting time series based on a behavior of other one

text: Apologies for this vague and unclear question, I have no background in statistics. I have two vectors of time series data, covering a six month period. The data is in daily intervals (except for weekends). The first vector contains the stock volume for a company on each of the days. The second vector contains an external volume variable (the number of Google searches) taken on the same date. ![Figure](http://i.stack.imgur.com/W2rqr.png) As visible on the graph above, (which is out of scale etc.) there is some type of relationship between the two. I want to model the two time series as to predict a future stock volume values. So far I have briefly looked into:   * support vector machines for regression: but when I plotted the data on a scatter graph there isn't any clear relationship, so this made me think that the SVM approach would fail.   * time series analysis: this approach confused me because there isn't any obvious trends in the data Can anyone offer any advice on what I approach I should take?


title: Number of segments to divide a time-series

text: Suppose we have time-series $ X_t $ and it has the following decomposition $$X_t=\mu + \varepsilon_t,$$ where $\mu$ is a mean and $\varepsilon_t$ - the error term. The model complexity will increase if we divide this time-series in to some segments,say $k$, and repeat above process. As the model complexity increases the approximation accuracy also increases. So I want to introduce a regularisation term here which will help in deciding the number of segments $k$ in which we need to divide the time-series. The error in approximation can be defined as $$ \epsilon_t= \frac{1}{k} \sum_{i=1}^k (\mu_{1}-X_{i})+\frac{1}{n-k} \sum_{i=n-k}^n (\mu_{2}-X_{i}), $$ here I have divided the time-series in 2 segments and $\mu_{1}, \mu_{2}$ are their respective means. Now I want to find out the optimal number of segments in general. Please note, that here I want to introduce a "regularisation" term which will help in deciding optimal number of segments.


title: Comparing two binary distributions

text: I have a dataset with responses coded as 0 and 1. I am trying to specify 3 models of responses and compare it to the observed results. I would therefore like to make single comparisons between two binary distributions. Is there a way to do this? Basically what I want is to compare observed results for 16 question per participant with a model that has predictions for each 16 questions. I would like to see how much these predictions match the actual observed responses.


title: Linear Hypothesis for a quantile regression in r

text: I would like to test a linear hypothesis in a median regression model similar to example below.               require(AER);require(car);require(quantreg)     data("CPS1985")          #Regular linear model     model.lm <- lm(wage ~ ethnicity + age*gender,CPS1985)     summary(model.lm)     linearHypothesis(model.lm,"age + age:genderfemale = 0")          #Quantile regression     model.quant <- rq(wage ~ ethnicity + age*gender,tau=0.5,data=CPS1985)      I'm not sure, however, how to execute a linear Hypothesis for `model.quant` to test that the sum of the coefficients on age and age:gender are equal to zero.


title: Looking for a proof of this statement

text: Let the data matrix $\pmb X$ be a collection of $n$ $p$-vectors $x_i\sim\mathcal{E}_p(\mu,\pmb \varSigma)$ where $\mathcal{E}_p(\mu,\pmb \varSigma)$ is a square integrable, continuous elliptical density on $\mathbb{R}^p$ with location vector $\mu$ and scatter matrix $\pmb\varSigma$ and $n>p$. Then let $\\{S_m\\}_{m=1}^M$ be $p$ indexes drawn from $\\{1:n\\}$ with $M=\binom{n}{p}$. Then let $A_m=\\{x_i|i\in S_m\\}$ and $y_m\in\mathbb{R}^p:\pmb A_m'y_m=1_p$ ($1_p$ denotes a $p$-vector of ones). $$\underset{1\le m\le M}{\max}\;\underset{1\le i\le n}{\max}\;(x_i'y_m-1)^2||y_m||^{-2}\text{ is bounded}$$ I'm looking for a reference to a proof (presumably of a more general variant) of (or a counter-example to) this statement. # EDIT1: a '$p$ index' is a subset of $p$ elements of $\\{1:n\\}$ drawn randomly and without replacement. There are $\binom{n}{p}$ distinct such $p$-subsets of $\\{1:n\\}$. $m\in\\{1,\ldots,M\\}$. $x_i$ is indeed a row of $\pmb X$. $\pmb A_m$ is a $p$ by $p$ matrix formed of $p$ rows of $\pmb X$ with indexes in $S_m$. # EDIT2: following Mike McCoy's question. For a given $y_m$ the inner product $x_i'y_m=1$ only for the $x_i$ with index $i\in S_m$. Here is a short R code to illustrate.               n<-20     p<-3     X<-matrix(rnorm(n*p),nc=p)     S_m<-sample(1:n,p)     A_m<-X[S_m,]     y_m<-solve(A_m,rep(1,p)) #equation of the hyperplane through the members of S_m     dis<-tcrossprod(X%*%y_m-1)/as.numeric(crossprod(y_m))      stated otherwise, the claim is that the $\max$ (over all possible $M$ subsets $S_m$) of the sup norm of `dis` is bounded whenever $\pmb X$ is drawn from a square integrable and continuous density. # EDIT3: to address Mike McCoy's questions. By bounded, I mean that $$\exists U(\mu, \pmb \varSigma):\underset{1\le m\le M}{\max}\;\underset{1\le i\le n}{\max}\;(x_i'y_m-1)^2||y_m||^{-2}\le U(\mu, \pmb \varSigma)\;a.s.$$ for some positive scalar $U(\mu, \pmb \varSigma)$ depending only on $\pmb \varSigma$ and $\mu$. Again, I suspect this theorem (or even a more general version of it) has been proven. I just can't locate the reference.


title: When, if ever, to use pairwise deletion in multiple regression?

text: I received the following question by email: > I was wondering should I use tick the option for pairwise exclusion of > missing data when I carry out regression analyses (or any analyses for that > matter) rather than using [some other missing values replacement strategy]. > Julie Pallant recommends pairwise exclusion of missing data in her SPSS > textbook. I have a few thoughts, but I was interested in first hearing your thoughts.


title: $R^2$ from a regression of two trend-stationary processes, $Y_t$ and $X_t$

text: In Estimation and Inference in Econometrics, by Davidson and MacKinnon, p.671, they claim that $R^2$ from a regression of $Y_t$ on $X_t$, where both time series are trend stationary, tends to 1 as $n$ tends to infinity. Can anyone please provide a mathematical answer to this? The models for $Y_t$ and $X_t$ are: $Y_t = \delta_0 + \delta_1t + u_t$ and $X_t = \gamma_0 + \gamma_1t + v_t$, where $u_t$ and $v_t$ are stationary.


title: How do you make forest plot with difference in means and CI?

text: I have to prepare forest plot for anti-hypertensive drugs. So far, I have seen forest plots that used odds ratios. I have difference in mean change in blood pressure in 2 groups along with 95% confidence limits. Can anyone suggest me how to draw forest plot using the difference in means, standard error and confidence limits?


title: Help on Regression Analysis (ROC, multinomial)

text: I am currently looking at continuous data and seeking a correlation with the presence (and severity, ie, normal, mild, moderate, severe) of a disease entity. Now, so far, running Spearman on SPSS, I've gotten an R of 0.370, p<0.001. I need help on a couple of things. Can/Should I perform a multinomial regression analysis? Also, I'd like to generate a ROC and find out corresponding cut-off values that may pertain to the severity of the disease, i.e., Normal, Mild, Moderate, Severe. Please advise on how to go about this. Any input is highly appreciated.


title: Estimate of $\sigma^2$ in a simple linear regression when $H_0: \beta = 0$ is true

text: I was reading some lecture notes on simple linear regression where one section said that when the slope is 0 (hence, $H_0: \beta = 0$ is actually true), $\frac{(SSY - SSE)}{(DFY - DFE)}$ estimates sigma squared. This doesn't make sense to me. If beta is actually 0, then shouldn't SSY and SSE be the same value? $$SSY = ∑(Y-\bar{Y})^2$$ $$SSE = ∑(Y-\hat{Y})^2$$ However, if $\beta = 0$, that means the best estimate ($\widehat{\text{Y}}$) is essentially $\overline{\text{Y}}$, which in turn means $$SSY = SSE = ∑(Y-\bar{Y})^2$$ In other words $$\frac{SSY - SSE}{DFY - DFE} = \frac{0}{DFY - DFE}$$ Is there something that I am missing? I hope someone here can elaborate on this.


title: Quasibinomial vs negative binomial and hurdles

text: I have some over-dispersed data and am trying to decide which model would best suit the data. The data are usually counts of symptoms or number of correct items on some cognitive tasks. As an example:               set.seed(69)     g1<-rnorm(700,30,9); g2<-rnorm(100,25,7); g3<-rnorm(100,20,5)     gt<-data.frame(score=c(g1, g2, g3), fac1=factor(rep(c("a", "b", "c"), c(700, 100, 100))), fac2=ordered(rep(c(0,1,2), c(3,13,4))))     gt$score<-with(gt, ifelse(fac2 == 0, score, score-rnorm(1, 0.5, 2)))     gt$score<-with(gt, ifelse(fac2 == 2, score-rnorm(1, 0.5, 2), score))     gt$score<-round(with(gt, ifelse(score>=30, 30, score)))      gt$cov1<-with(gt, score + rnorm(900, sd=40))/40     gt$score.30<-with(gt, 30-score)      The models I'm thinking about using are:               glmnb1<-glm.nb(score.30~cov1 + fac1*fac2, data=gt)         hur1<-hurdle(score.30~cov1 + fac1*fac2, dist="negbin", data=gt)     quasi1<-glm(cbind(score, score.30)~cov1+fac1*fac2, family="quasibinomial", data=gt)        1. How to decide between the negative binomial and the quasibinomial?   2. In this example, the hurdle model is a better fit compared to the negative binomial. However, if the quasibinomial was better compared to the negative binomial (hypothetically or otherwise), how do you compare the hurdle and quasibinomial? Is there a hurdle quasibinomial?


title: Model selection in SAS surveyreg analysis of complex surveys

text: What are model selection methods are applicable in this setting? Is it possible to calculate AIC or BIC with SAS's `surveyreg` output?


title: Two binomial probability Researcher paradox

text: A researcher knows that the probability that a person will respond his letter is 10% If the researcher sends a post-paid letter the probability that the reciever will respond is 40%. The researcher sends 5 post-paid letters and 5 non post-paid letters. What is the probability that he recieves less then 3 responses? Answer should be : 0,5193 How can you calculate this?


title: Interpretation problems with hypothesis testing

text: Two things always bugged me about hypothesis testing:   1. The chance that the population mean is _exactly_ any given number (provided that the random variable in question is continuous) is always zero, is it not? Therefore, we should always reject the null hypothesis...   2. If the result of the test is whether to reject or accept the _null hypothesis_ , what difference does it make what the alternative hypothesis states? Please, anyone can shed some light?


title: lhs Unit Eigenvectors in Python for DeGroot model Convergence for Social Networking problems

text: The influence array at time 0 is:               1/3, 1/3 1/3     1/2, 1/2, 0     0, 0.25, 0.75      It is suppose to converge to the array:               3/11 4/11 4/11     3/11 4/11 4/11     3/11 4/11 4/11      You are suppose to be able to get the convergence array by taking the left hand side (lhs) unit vector of the array at time 0. My linear algebra can be rusty at times and I want to try this in python:               import numpy as np     from scipy.linalg import eig     np.set_printoptions(precision=4)     limT = np.array([(0.33,0.33,0.33),(0.5,0.5,0.0),(0.0,0.25,0.75)])     w, vl, vr = eig(limT, left=True)     lhs = vl.conj().T.dot(limT)     >>> lhs     array([[-0.0678,  0.0335,  0.0335],            [-0.4669, -0.6231, -0.6231],            [-0.2958, -0.1006,  0.3904]])     lhs/np.linalg.det(lhs)     array([[ -2.385 ,   1.1795,   1.1795],            [-16.4266, -21.9221, -21.9221],            [-10.4052,  -3.5378,  13.7349]])      The final array is obviously not:               3/11 4/11 4/11     3/11 4/11 4/11     3/11 4/11 4/11      I could have misunderstood several issues and appreciate the assistance.


title: Cell Means Model Property

text: Could anybody demonstrate or direct me to a readily available proof of the following: For the cell means model: $$ y_{ij} = \mu_{i} + \epsilon_{ij},\ \text{ for }\ i = 1, \ldots, r\ \text{ and }\ j = 1, \ldots, n_{i}.$$ Show that: $$ \sum_{i=1}^{r}\sum_{j=1}^{n_{i}}(y_{ij}-\bar{y}_{\cdot\cdot})^{2} = \sum_{i=1}^{r}\sum_{j=1}^{n_{i}}(y_{ij}-\bar{y}_{i\cdot})^{2} + \sum_{i=1}^{r}n_{i}(\bar{y}_{i\cdot}-\bar{y}_{\cdot\cdot})^{2} $$ where the first term is $\text{SS}_{\text{TO}}$, the second $\text{SS}_{\text{E}}$ and the third $\text{SS}_{\text{TR}}$. (I have an upcoming examination in linear models and this proof was required for one of the previous year's examinations, but, so far, I haven't had much success finding the proof.)


title: Why are these hierarchical linear regression results in R and SPSS different?

text: When I try to use the data and example for HLR from this example dataset taken from this post in the R Tutorial Series on hierarchical linear regression the results don't match when I try to use the same method in SPSS. Is it because SPSS is using a different type of sum of squares (III)? The F values match for the final model, but not the 2nd one and some of the sum of squares seem off.               #R method     url <- "http://dl.dropbox.com/u/10246536/Web/RTutorialSeries/dataset_hlr.csv"     datavar <- read.csv(url, header=T)          #create three linear models using lm(FORMULA, DATAVAR)     #one predictor model     onePredictorModel <- lm(ROLL ~ UNEM, datavar)     #two predictor model     twoPredictorModel <- lm(ROLL ~ UNEM + HGRAD, datavar)     #three predictor model     threePredictorModel <- lm(ROLL ~ UNEM + HGRAD + INC, datavar)          #get summary data for each model using summary(OBJECT)     summary(onePredictorModel)     summary(twoPredictorModel)     summary(threePredictorModel)          #compare successive models using anova(MODEL1, MODEL2, MODELi)     test<- anova(onePredictorModel, twoPredictorModel, threePredictorModel)      Below here is the code for SPSS.               *SPSS method     data list free /YEAR ROLL UNEM HGRAD INC.     begin data     1   5501    8.1 9552    1923     2   5945    7   9680    1961     3   6629    7.3 9731    1979     4   7556    7.5 11666   2030     5   8716    7   14675   2112     6   9369    6.4 15265   2192     7   9920    6.5 15484   2235     8   10167   6.4 15723   2351     9   11084   6.3 16501   2411     10  12504   7.7 16890   2475     11  13746   8.2 17203   2524     12  13656   7.5 17707   2674     13  13850   7.4 18108   2833     14  14145   8.2 18266   2863     15  14888   10.1    19308   2839     16  14991   9.2 18224   2898     17  14836   7.7 18997   3123     18  14478   5.7 19505   3195     19  14539   6.5 19800   3239     20  14395   7.5 19546   3129     21  14599   7.3 19117   3100     22  14969   9.2 18774   3008     23  15107   10.1    17813   2983     24  14831   7.5 17304   3069     25  15081   8.8 16756   3151     26  15127   9.1 16749   3127     27  15856   8.8 16925   3179     28  15938   7.8 17231   3207     29  16081   7   16816   3345     end data.               REGRESSION /MISSING LISTWISE      /STATISTICS COEFF OUTS R ANOVA CHANGE      /CRITERIA=PIN (.05) POUT(.10)      /NOORIGIN /DEPENDENT ROLL      /METHOD=ENTER UNEM      /METHOD=ENTER HGRAD      /METHOD=ENTER INC.      Or did I mess something up in the procedure for SPSS?


title: Use raking to impute multivariate distribution

text: ## The Problem I have a detailed dataset of 6 variabes, but for all but one year, I only have marginal distributions of 5 variables, the rest is missing. From that, I would like to obtain a full dataset. To illustrate my problem imagine the case if I had a dataset of 3 variables (Age, Level, Year) and had only 2 margins for all but one year. 2011 (the “full year")                      Lvl1  Lvl2  Lvl3  Lvl4  Lvl5  Lvl6  Lvl7  Lvl1-7     1-20   2758  3531  2759  1472   963  1025   790   13298     21-40  2421  3530  2371  1583  1119  1376   720   13120     41-60  3497  3322  1818  1276   710   562   240   11425     61-80  4122  5554  2754  1786  1180   496   282   16174     80+    2740  5116  3050  2333  1889   656   334   16118     Sum   15538 21053 12752  8450  5861  4115  2366   70135      2010 (only 2 marginals)                      Lvl1  Lvl2  Lvl3  Lvl4  Lvl5  Lvl6  Lvl7  Lvl1-7     1-20                                              13197     21-40                                             12922     41-60                                             11369     61-80                                             16353     80+                                               15774     Sum   15151 21643 12611  8273  5586  4026  2325   69615      2009 (only 2 marginals)                      Lvl1  Lvl2  Lvl3  Lvl4  Lvl5  Lvl6  Lvl7  Lvl1-7     1-20                                              12684     21-40                                             12477     41-60                                             10844     61-80                                             16022     80+                                               14911     Sum   14367 20999 12201  8032  5235  3841  2263   66938      In this example, the challenge would be to fill in the missing values of the 2010 and 2009 data, using the 2011 distribution to make all necessary assumptions on multivariate marginal distributions. _Note that my problem is somewhat more complex. I have 5 dimensions (Age, Level, Region, Gender, Year) and three 4-dimensional marginal distributions tables (Gender, Level, Region, Year), (Age, Gender, Region, Year) and (Gender, Age, Level, Year) – but the problem is basically the same (I think … ?)._ ## How can that be done? Here’s my theory. _Survey statistics_ has done a lot of reaserach on how to make any given (sample) data agree with population margins, I could use their methods. They do so by reweighting observations. **Raking** can be used to obtain a set of such weights. The `rake` function in the `survey` package allows raking with multidimensional joint distributions. My idea was to use the full data as “sample”, and let the rake function compute weights that make the data agree with the marginal distributions. **Question 1:** Is that a good idea? Do you have a better idea? As you can imagine, there are several other problems to be solved. First one is that survey methods always assume that the “sample” data is individual data, I have _counts_ , and no other characteristics of the individuals other than the stratification variables. What I tried to do is specifying the counts as pre-raking survey weights. I tried that but the raking results in unplausible weights, that is, most weights are 0 and some are Inf. Another possibility would be to duplicate every entry of the “population” data.frame depending on the number of observations, but that would bloat the data.frame significantly (it would have approximately 18 000 000 rows instead of currently 24 480 (the product of all dimension sizes) and cause memory problems. If I specify no weights (all probabilities are 1), I get the following distribution of weights after raking:               > mean(lpg.rake$prob)     [1] 0.1289419     > median(lpg.rake$prob)     [1] 0.01700375     > quantile(lpg.rake$prob,c(.075,.1,.5,.9,.95,.99,.999))             7.5%          10%          50%          90%          95%          99%        99.9%       0.001274434  0.001697273  0.017003749  0.155742479  0.322445467  1.597854446 11.711632910      Given that the mean of frequencies in the one year for which I have a full dataset is around 700, I would expect that to be the average of the weights. **Question 2:** How should I deal with the fact that survey methods assume that observations are individuals, whereas I have one row for every combination of stratification variables plus a “count” variable? ## A working example Since my problem is very complex, and I am not familiar with the `survey` package, it is difficult for me to produce a minimal example that reproduces my problem. In some way, the complexity is the main problem. There is the example data in the `survey` package, you could take a look at that (link). However, since my data is not secret, I don’t see why I could not give you the full dataset and what I have tried so far. Here it is. The password for the zip file is: xB2KzbuK Source: Ministry for Social Affairs (1995-2012): Pflegevorsorgebericht / Bericht des Arbeitskreises für Pflegevorsorge. Vienna, Austria AND own calculations. Consequently, **question 3** is: Is there an error in my code? What could I do in order to get meaningful results?


title: Identifiying performance of individual feature components in accuracy of support vector methods

text: I'm using support vector regression (not classification) for a problem and it's working well. However, in the older method that former lab members developed (a basic linear model, with weights determined just from OLS), there is a lot of code used to determine the individual efficacy of each different right-hand-side variable from their OLS. I'm not all that impressed with these efficacy stats (they are mostly summary stats of correlations, t-stats, etc.) but they do allow you to have some notion of which RHS variables are most informative. In my case, for each individual datum, I use the whole row's worth of right- hand-side variables as the feature vector for support vector regression. So I am wondering if there is any commonly used analogue, maybe some measure from information theory or something, that gives an indication of how each marginal component of the overall feature vector contributes to the end performance. One thing that does not seem to work is training N different regressions each using only a single one of the feature components as the sole predictor. This does give some matches between the univariate regression efficacy and univariate support vector efficacy, but it doesn't capture the interactions between the different predictors, and (for example) which predictor is the "second best" predictor among all the components of the feature vector. I suppose I could work backward and test the results of removing one component at a time, to see which then causes the biggest drop in performance... but I'd prefer something less ad-hoc and perhaps backed by some statistical rigor or intuition.


title: User Defined Substitution Cost Matrix

text: I have a 30x30 symmetric matrix that I would like to use as a substitution cost matrix in TraMineR to analyze sequences of length 10 with an alphabet of 30. When I try to perform OM with this matrix, however, I get an error about the triangle inequality (see below). Other matrices that I have made by hand seem to work fine (see example below). I thought the subcost matrix could be set any way I want as long as it was symmetric and the dimensions reflected the size of the alphabet. Are there limitations on how I can set the subcosts? Or perhaps something else is causing the problem? Thanks               > diss1vs <- seqdist(bhpsfup.seqw, method = "OM", indel = .1*maxsub, sm = upsubcosts)      [>] 536 sequences with 30 distinct events/states     Error in checktriangleineq(sm, warn = FALSE, indices = TRUE, tol = tol) :        REAL() can only be applied to a 'numeric', not a 'integer'      Example matrix that works               zeroblock <- diag(6)     zeroblock[zeroblock == 1] <- 0     zeroblock          twoblock <- diag(6)     twoblock[twoblock == 0] <- 2     twoblock[twoblock == 1] <- 2     twoblock          submat.up1 <- cbind(zeroblock,twoblock,twoblock,twoblock,twoblock)     submat.up2 <- cbind(twoblock,zeroblock,twoblock,twoblock,twoblock)     submat.up3 <- cbind(twoblock,twoblock,zeroblock,twoblock,twoblock)     submat.up4 <- cbind(twoblock,twoblock,twoblock,zeroblock,twoblock)     submat.up5 <- cbind(twoblock,twoblock,twoblock,twoblock,zeroblock)     submat.up <- rbind(submat.up1, submat.up2, submat.up3, submat.up4, submat.up5)     submat.up