title: Properties of MLE and least squares methods for estimating parameters of ar(ma) models

text: I have annual data that seem to have a bimodal density function. My explaination is that there is a distinction between wet and dry years. For my work I would like to use an ar(1)-model for this. When I simulate from a fitted ar(1)-model with normally distributed innovations I of course get normally distributed data. Since the simulation results look different from the historical data this is not satisfactory. When I generate my innovations by just sampling from the residuals which also show a bimodal density function I get better results concerning the density function. By doing so my innovations are still white noise. But not Gaussian white noise any more. This approach is called "model based resampling" and was (I guess) introduced by Freedman (1984). That approach has some disadvantages but for me it seems to be an improvement. But I am not certain if I still can use (standard) MLE to estimate the ar- parameter(s) in the model based resampling approach. Or if I might use least squares instead. I suppose that least squares is better since there is no gaussian assumption. I was really looking for publications for this but the few I found are extremely difficult to understand. I found a lot concerning fitting of distributions and structural equation modeling, saying that small deviations from normal are not a problem for MLE as long as there are no extreme outliers. But I guess time series may not be the same. Can you help me?


title: WinBUGS/OpenBUGS Gamma distribution parameters

text: According to the WinBUGS manual the gamma distribution is defined by: `dgamma(r,mu)` However, what is `r`? Is it the shape, scale or rate parameter? I'm pretty new to statistics and googling didn't really help. Above that in most explanations I've found `r`is usually not mentioned. It's either k, `theta`, `beta` or `alpha`.


title: Implementation of Univariate logistic regression

text: For Univariate Linear Regression I can calculate the parameters (And most everything else) from simple sum of squares. Is there a corresponding method for Logistic Regression? Any pointers to code (in any language) would be helpful. (Yes there are many solvers, but I want to implement as a simple Map-Reduce algorithm, as I have for Univariate Linear Regression)


title: k-core clustering algorithm

text: I am trying to cluster data. Each point in this dataset is connected to some other points. I want to define clusters "depending on how much the points are connected to each other". After some research, I read about k-core clustering (and saw it has applications in social networking for instance). I think this is the algorithm I want to apply, but, according to what I found, it can only be used to "visualize" a network. Isn't it possible to cluster data with this method?


title: Checking whether a given sample is stationary

text: Let $x[n]$ be some time series in 1D or $x[m,n]$ in 2D, of length $N$ (resp. $N^2$) How can I assess whether it is _stationary_? At least in the weak sense. I can check whether the _stdev_ remains constant, but this is only a necessary condition.


title: Memory consumption in NLP tasks

text: I am trying to apply a simple Naive Bayes or SVM (libSVM) algorithms to a large data set, which I've constructed as an .arff file. The number of features in my set is ~180k and there are ~6k examples. Also there are 8 classification classes. The data is of size ~3.2GB. I am working with Weka's Java API and Eclipse, I am increasing JVM's memory to the maximum, but I am always getting a heap space error. I am on a MacBook Pro, 2.3 GHz Intel Core i5, 4GB 1333 MHz DDR3. Do I need to find another machine to work with or is it possible that I am having an memory leaks programmatically?


title: Divergence from uniform distribution (continuous): dispersion measure?

text: I have data of a continuous random variable within the range [-1,1], which sometimes is concentrated around 0, while other times is concentrated toward -1 and 1, while zero is relatively underpopulated. What measure can I use for both these cases to measure the divergence from a uniform distribution? In other terms: I am looking for a measure of how evenly spread out the data is within the range, but standard dispersion measures (like variance) don't seem to work, since they favor distributions in which the tails are higher than the 'peak', e.g., when the region around zero is relatively underpopulated.


title: "Bayesian and Frequentist Regression Methods" by Jon Wakefield, a good introductory Bayesian textbook for frequentist economics graduates?

text: Here's a link to a good question regarding Textbooks on Bayesian statistics from some time ago. People suggested John Kruschke's "Doing Bayesian Data Analysis: A Tutorial Introduction with R and BUGS" as one of the best options to get an introduction to Bayesian statistics. Meanwhile, a potentially interesting book called "Bayesian and Frequentist Regression Methods" by Jon Wakefield was released, which also provides code for R and BUGS. Thus, they esentially both seem to cover the same topics. Question 1: If you have read the book, would you recommend it to a frequentist economics masters graduate as both an introduction to Bayesian statstics and reference book for both frequentist and bayesian approaches? Question 2: If you have read both Wakefield's and Kruschke's book, which one would you recommend better?


title: What's a distribution with an abyss instead of a peak?

text: I am looking for a (commonly used) probability density function, which would look like a normal distribution flipped upside down. It would look like a uniform distribution with a dent in the middle. Just to be clear, I am dealing with a continuous random variable within some range, say `[-1,1]`. Sometimes I have data which is concentrated around zero, but other times I have data which is concentrated toward 1 and -1, while the region around zero is relatively underpopulated. Is there a kind of pdf which (depending on some parameter) can represent these two cases?


title: REML estimation - matrix algebra internals

text: In the context of REML estimation there is the result (ignoring some constants) that (my interest is in the matrix algebra so some notation is suppressed): $l(\mathbf V_0)=\log |\mathbf V_0| + \text{tr}(\mathbf V_0^{-1}\mathbf S) \tag{1}$ Where both $\mathbf V_0$ and $\mathbf S$ are symmetric and invertible. I am told that the following expression can be obtained from (1) by differentiating with respect to a parameter $\sigma_i$ of $\mathbf V_0$: $\text{tr}\left[\mathbf V_0 \frac{\partial \mathbf V_0^{-1}}{\partial \sigma_i}\right]-\text{tr}\left[\frac{\partial \mathbf V_0^{-1}}{\partial \sigma_i}(\mathbf V_0 - \mathbf S)\right] \tag{2}$ I'm trying to see what steps got this, but am stuck. Now I'm obviously lacking some machinery for dealing with these things but I don't know where to look. The issue is that, say, both terms in (1) are scalar, and we are differentiating w.r.t a scalar, so the it would seem the answer should be scalar. But my workings end up matrix valued. E.g. $\frac{\partial \log |\mathbf V_0|}{\partial \sigma_i} = \mathbf V_0^{-1} \frac{\partial \mathbf V_0}{\partial \sigma_i} $ Clearly the trace gets involved wrapping it all up into a scalar, but I don't know how or why!


title: Can the standard deviation be calculated for harmonic mean?

text: Can the standard deviation be calculated for the harmonic mean? I understand that the standard deviation can be calculated for arithmetic mean, but if you have harmonic mean, how do you calculate the standard deviation or CV?


title: How to handle imbalanced data using SMOTE-N (nominal data) to generate multiple synthetic data?

text: I have a classification problem with two classes working on nominal data. I want to apply SMOTE-N to deal with imbalanced data. However, it is not clear to me how to use SMOTE-N for generating N synthetic data for each feature vector in the minority class. SMOTE-N uses a modified version of the value difference metric (VDM) to find the k-nearest neighbors for each feature vector in the minority class and then the new minority class feature vector is generated by creating new set feature values by taking the majority vote of the feature vector in consideration and its k nearest neighbors (k-nn). But, how is this process repeated to generate multiple synthetic feature vectors for each feature vector in the minority class? The way the algorithm is stated, it seems that one feature vector from the minority class can generate only one synthetic feature vector (using its K-nn)?


title: Can I use one conjugate prior for samples from two normal distributions with same variance but different mean?

text: I have got a set $\{Y_t\}$ of observations consisting of two subsets $\{Y_{t,1}\}$ and $\{Y_{t,2}\} \subset \{Y_t\}$ with $\{Y_{t,1}\} \sim \mathcal{N}(\mu_1,\sigma^2)$ and $\{Y_{t,2}\} \sim \mathcal{N}(\mu_2,\sigma^2)$ i.e. _different means but the same variance_ (resulting from a regime switching model). I know the means and want to draw a sample of $\sigma^2$ in a step of a MCMC estimation. In the case of $\mu_1 = \mu_2$ I would have used the conjugate prior of the Inverse Gamma distribution (see [1], "Normal with known mean"). Can I use a conjugate prior in the case of $\mu_1 \neq \mu_2$ as well? For example by setting $\beta = \beta_0 + \frac{1}{2}\sum_{i=0}^n(Y_i - \mu_{I_i})^2 $ ($I_i$ being the correct indices according to the observations)? Or will I have to use Metropolis-Hastings to get my sample of $\sigma^2$? Best, Matt [1] http://en.wikipedia.org/wiki/Conjugate_prior#Continuous_distributions


title: Problem with classifier prediction results

text: I built a classifier with 13 features ( no binary ones ) and normalized individually for each sample using scikit tool ( Normalizer().transform). When I make predictions it predicts all training sets as positives and all test sets as negatives ( irrespective of fact whether it is positive or negative ) What anomalies I should focus on in my classifier, feature or data ??? Notes: 1) I normalize test and training sets (individually for each sample) separately. 2) I tried cross validation but the performance is same 3) I used both SVM linear and RBF Kernels 4) I tried without normalizing too. But same poor results 5) I have same number of positive and negative datasets ( 400 each) and 34 samples of positive and 1000+ samples of negative test sets.


title: What is the probability that random variable $x_1$ is maximum of random vector $X=(x_i)$ from a multivariate normal distribution?

text: Given a $n$-dimensional multivariate normal distribution $X=(x_i) \sim \mathcal{N}(\mu, \Sigma)$ with mean $\mu$ and covariance matrix $\Sigma$, what is the probability that $\forall j\in {1,\ldots,n}:x_1 \geq x_j$?


title: It is possible to implement a regression algorithm using Random Forests?

text: I have to write code to implement regression using random forests (by default Weka provides random forests for classification). Is this possible to do?


title: What is the intuition behind the independence of $X_2-X_1$ and $X_1+X_2$, $X_i \sim N(0,1)$?

text: I was hoping someone could propose an argument explaining why the random variables $Y_1=X_2-X_1$ and $Y_2=X_1+X_2$, $X_i$ having the standard normal distribution, are statistically independent. The proof for that fact follows easily from the MGF technique, yet I find it extremely counter-intuitive. I would appreciate therefore the intuition here, if any. Thank you in advance. **EDIT** : The subscripts do not indicate Order Statistics but IID observations from the standard normal distrubution.


title: Comparing model fits across a set of nonlinear regression models

text: **CONTEXT:** I am modelling the relation between time (1 to 30) and a DV for a set of 60 participants. Each participant has their own time series. For each participant I am examining the fit of 5 different theoretically plausible functions within a nonlinear regression framework. One function has one parameter; three functions have three parameters; and one function has five parameters. I want to use a decision rule to determine which function provides the most "theoretically meaningful" fit. However, I don't want to reward over-fitting. Over-fitting seems to come in two varieties. One form is the standard sense whereby an additional parameter enables slightly more of the random variance to be explained. A second sense is where there is an outlier or some other slight systematic effect, which is of minimal theoretical interest. Functions with more parameters sometimes seem capable of capturing these anomalies and get rewarded. I initially used AIC. And I have also experimented with increasing the penalty for parameters. In addition to using $2k$: [$\mathit{AIC}=2k + n[\ln(2\pi \mathit{RSS}/n) + 1]$]; I've also tried $6k$ (what I call AICPenalised). I have inspected scatter plots with fit lines imposed and corresponding recommendations based on AIC and AICPenalised. Both AIC and AICPenalised provide reasonable recommendations. About 80% of the time they agree. However, where they disagree, AICPenalised seems to make recommendations that are more theoretically meaningful. **QUESTION:** Given a set of nonlinear regression function fits:   * What is a good criterion for deciding on a best fitting function in nonlinear regression?   * What is a principled way of adjusting the penalty for number of parameters?


title: Entering variables in multivariate logistic regression and running regression across two groups

text: I am trying to do a multiple logistic regression for 2 similar groups. I have a few questions:   1. In doing a univariate analysis, do I enter each independent variable, one at a time, first into the binary regression, before going on to do the multivariate analysis? Or is the significance values from Chi-square or t-test enough to go on with?   2. I have a test group and a control group and I want to determine the effect of independent variables (e.g HIV status, maternal weight etc) on a particular dependent variable (low birth weight). Should I perform the regression on a dataset with both the test and the control cases or should I split the file? In this case I want to see the effect of HIV on birthweight and I am having a hard time knowing how to move on.


title: Combining two autocorrelation functions measured on same stationary process

text: I have two autocorrelation functions measured from the same stationary process, but the two measurements were taken with different instruments that measure different lag times. I would like to combine these two autcorrelation functions together into one single curve that spans the entire lag time of both instruments. It seems that these could be combined in the frequency domain with Wiener-Khinchin theorem, but I have not gotten anything reasonable to come out from this. Thanks! ![Two autocorrelation functions](http://i.stack.imgur.com/soTIx.jpg) The autocorrelation data: Short lag times: {{0.0166638, 1.2427}, {0.0333277, 1.16926}, {0.0499915, 1.18007}, {0.0666553, 1.17344}, {0.0833192, 1.21829}, {0.099983, 1.19867}, {0.116647, 1.14627}, {0.133311, 1.17827}, {0.166638, 1.19614}, {0.199966, 1.19341}, {0.233294, 1.18352}, {0.266621, 1.18402}, {0.333277, 1.18672}, {0.399932, 1.18258}, {0.466587, 1.15333}, {0.533243, 1.17556}, {0.666553, 1.17179}, {0.799864, 1.17035}, {0.933175, 1.18451}, {1.06649, 1.16379}, {1.33311, 1.14078}, {1.59973, 1.13816}, {1.86635, 1.13299}, {2.13297, 1.13205}, {2.66621, 1.09647}, {3.19946, 1.09374}, {3.7327, 1.07922}, {4.26594, 1.06006}, {5.33243, 1.04623}, {6.39891, 1.03004}, {7.4654, 1.02369}, {8.53188, 1.01798}, {10.6649, 1.0145}, {12.7978, 1.00935}, {14.9308, 1.00495}, {17.0638, 1.00267}, {21.3297, 1.00247}, {25.5956, 1.00451}, {29.8616, 1.00155}, {34.1275, 1.00223}, {42.6594, 0.99827}, {51.1913, 1.00015}, {59.7232, 0.996647}} Longer lag times: {{1.483, 1.2196}, {2.966, 1.1595}, {4.4489, 1.1353}, {5.9319, 1.1261}, {7.4149, 1.1126}, {8.8979, 1.0907}, {10.381, 1.0796}, {11.864, 1.0697}, {14.83, 1.0471}, {17.796, 1.0397}, {20.762, 1.0229}, {23.728, 1.0199}, {29.66, 1.0097}, {35.592, 1.0077}, {41.524, 1.0048}, {47.455, 1.004}, {59.319, 1.0023}, {71.183, 1.0008}, {83.047, 1.0011}, {94.911, 0.99988}, {118.64, 1.0013}, {142.37, 1.001}, {166.09, 0.99821}, {189.82, 0.99895}, {237.28, 1.0015}, {284.73, 1.0006}, {332.19, 0.99739}}


title: Check that state space model implemented correctly

text: I want to simulate data from the following model: $\textbf{z}_k=\textbf{H}\textbf{x}_k+\textbf{v}_k$ $\textbf{v}_k \sim N(\textbf{0},\textbf{R})$ $\textbf{H}$ does not change over time   $\textbf{x}$ is a vector of loadings   $\textbf{R}$ is a diagonal of constants $\textbf{x}_k=\textbf{F}\textbf{x}_{k-1}+(\textbf{I}-\textbf{F}){\mu} + \textbf{w}_k $ $\textbf{w}_k \sim N(\textbf{0},\textbf{Q})$ $\textbf{I}$ is the identity matrix   $\mu$ is the vector of mean values of $\textbf{x}$   $\textbf{F}$ is diagonal with the AR(1) params which do not change over time   $\textbf{Q}$ is diagonal with the innovation processes for $\textbf{x}$ I have the following code in Matlab               nDates=20000; %number of dates          mats=[1 2 3 4 5 6 7 8 9 10 12 15 20 25 30]'; %maturities      nY=length(mats); %#number of yields          z=zeros(nY,nDates);  %declare vector for yields     x=zeros(3,nDates); %declare vector for factors          R=0.00001;   %standard deviation     I=eye(3); %3*3 identity matrix          v=normrnd(0,R,nY,nDates); %generate residuals     F=[0.9963 0 0; 0 0.9478 0; 0 0 0.774]; %AR(1) matrix     mu=[0.0501; -0.0251;-.0116]; %mean of X     lambda=0.5536;           q = [0.0026^0.5 0 0;0 0.0027^0.5 0; 0 0 0.0035^0.5];     Q=q*q';     rng('default');  % For reproducibility     r = randn(nDates,3);     w= (r*Q)';                         B= [ones(nY,1),((1-exp(-lambda*mats))./(lambda*mats)),((1-exp(-lambda*mats))./(lambda*mats))-exp(-lambda*mats)];          X(:,1)=mu;          for t=2:nDates     x(:,t)=F*(x(:,t-1))+(I-F)*mu+w(:,t);     z(:,t)=B*x(:,t)+v(:,t);     end          z(:,1)=[];      It all seems straightforward enough but what test can I so to ensure that it has been implemented correctly? Ones I have thought of:   Check the correlation of the factors x on their 1 period lags match the values given in matrix F   Check that the variances of v and w are correct   Check that the mean of the simulated variables are correct I would like to check that the empirical variance of the parameters matches their theoretical equivalent, but I don't know what the theoretical equivalent should be? Please feel free to suggest further tests that will allow me to know for sure if the implementation is correct.


title: Variance of data from small versus large countries

text: In considering (for example) homicide rates reported annually for each country, it occurs to me that the U.S., being larger in population, might thereby be expected to have a rate closer to the worldwide mean. However, I don't know how to evaluate that idea quantitatively. What kind of additional data (if any) and calculations would be required to support that idea (not sure I should call it a hypothesis, because it relates to the model rather than to the world). Another way of thinking about this: with respect to statistics reported per country, does the U.S. look like a group of smaller countries combined?


title: Efficient way to classify with SVM

text: I'm doing a binary classification using SVM classfier, libsvm, where roughly 95% belongs to one class. The parameters C and gamma are to be set before the actual training takes place. I followed the tutorial but still can't get any good results. There is a script that comes with the library that is supposed to help with choosing the right values for parameters but what this script is doing is basically maximizing the accuracy metric (TP+TN)/ALL, so in my case it chooses the parameters to label all data with prevailing class label. I would like to choose parameters with recall and precision based metrics. How could I approach this problem. Accuracy is a meaningless metric for what I'm doing. Also I'm keen on changing the library libsvm to any other one that can help me with this problem as long as it takes data in the same format. 1 1:0.3 2:0.4 ... -1 1.0.4 2:0.23 and so on Can anybody help? UPDATE: yes I did try both grid.py and easy.py but even though grid search uses logarithmic scale it is extremely slow.I mean even if I run it on just small chunk of my data it takes tens of hours to finish. Is this the most efficient way to use SVM?? Have also tried svmlight but it does exactly the same it labels all data with one label. UPDATE2: I reformed my question the better reflect what sort of issues I am facing


title: Taking random items out of a container with replacement

text: I take vitamins in the morning, but one of them I only take a half tablet. So, I have an initial container with 100 full tablets, and every morning I take out a random tablet. If it's a full tablet, I break it in half, put half back, and take half. If it's a half tablet, I just take it. Given that, how many days must I do this before I get a $> 50\%$ chance of getting a half tablet? (Or, what's the % chance I get a half tablet after $X$ days?) **NB** : This isn't homework; it just occurred to me as I was doing my morning routine, and I really have no idea where to start on trying to solve this. Pointers welcome. It feels like something with an infinite series, but I'm not sure.


title: Combining dummy and continuous items in a single scale

text: In a predictive linear regression model, does it make sense to combine binary and continuous data within a single scale?   Does it make statistical sense to convert binary data into a z score? I am wondering whether I could create a scale using a combination of binary and continuous items that I could convert into a Z score.


title: Ideas for visualising and analysing grocery store data

text: I have gotten a student job in the management department of a chain of 50 grocery stores. The job includes gathering daily statistics on the economy of the stores. Every day a **per-store revenue statistic** is made, comparing the revenue to matching day last year (sunday/sunday, monday/monday) and an accumulated statistic for the month. Also a **per-store gross margin percentage statistic** is made. Now these are simple measures and I feel that there could be a lot more to be gained from the data, the data is specific down to the (cost price)/(sale price)/(number of sales)-level on every type of item. Furthermore no or very few graphics are made and the data from the previous years are not taken into use. Have any of you seen a similar problem? Do you have any ideas on how to proceed? Any easy-to-read, informative types of graphs you would want to share regarding these types of data?


title: Standard error of measurement versus minimum detectable change

text: I am currently calculating reliability estimates for test-retest data. My question is regarding the difference between standard error of measurement (SEM) versus minimum detectable change (MDC) when seeking to determine if there is a 'real' difference between two measurements. Here is my thinking thus far: Each measurement has an error band about it. For two measurements, if error bands overlap then there is no 'real' difference between the measurements.   1. For example, at 95% confidence, each measurement has an error band of $\pm 1.96 \times SEM$. So, two measurements would need to be more than $2 \times 1.96 \times SEM =3.92 \times SEM$ apart to avoid each measurement's confidence interval overlapping and for their to be a real difference between the two measurements.   2. Another method for determining if two measurements are 'different' is to use MDC where  $$MDC = 1.96 \times \sqrt{2} \times SEM =2.77 \times SEM$$ [EDIT: for second formula see e.g. p. 238 of Weir, J. P. (2005). Quantifying test-retest reliability using the intraclass correlation coefficient and the SEM. Journal of strength and conditioning research / National Strength & Conditioning Association, 19(1), 231–240. doi:10.1519/15184.1] If the difference between the two measurements is greater than MDC then there is a real difference between the measurements. Obviously the two formulas are different and would produce different results. So which formula is correct?


title: Significance test for differences in baseline and remeasurement rates between two groups

text: I need help in knowing which test to perform with this analysis. I have two groups (lets say group A and B). A study was performed with initial rates of heart failures in group A and B (i.e, baseline rates) were measured. After 6 months of follow up the heart failure rates in group A and B (i.e, re measurement rates) were measured again. By simple mathematics, I know that the improvement rate is better in group A as compared to group B.   * How do I test whether the difference in rates between the groups is statistically significant? I am using t-test for this but not sure how to go about it.


title: Asymptotic distribution of sample variance of non-normal sample

text: This is a more general treatment of the issue posed by this question. After deriving the asymptotic distribution of the sample variance, we can apply the Delta method to arrive at the corresponding distribution for the standard deviation. Let a sample of size $n$ of i.i.d. **non-normal** random variables $\{X_i\},\;\; i=1,...,n$, with mean $\mu$ and variance $\sigma^2$. Set the sample mean and the sample variance as $$\bar x = \frac 1n \sum_{i=1}^nX_i,\;\;\; s^2 = \frac 1{n-1} \sum_{i=1}^n(X_i-\bar x)^2$$ We know that $$E(s^2) = \sigma^2, \;\;\; \operatorname {Var}(s^2) = \frac{1}{n} \left(\mu_4 - \frac{n-3}{n-1}\sigma^4\right)$$ where $\mu_4 = E(X_i -\mu)^4$, and we restrict our attention to distributions for which what moments need to exist and be finite, do exist and are finite. Does it hold that $$\sqrt n(s^2 - \sigma^2) \rightarrow_d N\left(0,\mu_4 - \sigma^4\right)\;\; ?$$


title: Describing the difference between 2 types of variance

text: In my research I'm comparing the variance of a method and I would like to describe the overall variance between individuals and the variance of the replicates of these individuals. Things like 'comparing the intra-individual variance and between-individual variance' seems to get people confused. I would like to make a short brief notice of this without having to go to much in details about the experiment. What would be a way of describing this setting more clearly but still within if possible one sentence? To clarify: I have 10.000 measurements for 60 individuals. For each measurement I could calculate for example the standard deviation as a method of variance. I also have 5 replicate measurements per individual. I could calculate the standard deviation for each of the 10.000 measurement in the replicates. So now I have the variance of the measurement when looking in a population AND I have the variance when looking at replicates. When you would now have to describe these 2 types of variance in a single sentence how would you do that without going into to much details?


title: If $Z=X+Y$ ($X$ and $Y$ being random samples), what is the relationships between the respective order statistics?

text: Suppose $X$ and $Y$ are two random samples (not necessarily iid, but one can make this assumption) and that $Z=X+Y$. If one computes the order statistics of $X$ and $Y$, what can be said about the relative order statistic of $Z$? To be more clear, let $\tilde{X}_{0.99}$ and $\tilde{Y}_{0.99}$ be the 0.99th quantiles of $X$ and $Y$, respectively, does it exist a relationship $f(\cdot)$ with $\tilde{Z}_{0.99}$ (i.e., the 0.99th quantile of $Z$) such that $\tilde{Z}_{0.99}=f(\tilde{X}_{0.99},\tilde{Y}_{0.99})$? Sorry for the possibly ill-posed question... I'm not a statistician.


title: How is L2 Boosting Different from a Big Regression Tree?

text: I'm learning about boosting. I think I understand how adaptive boosting works for classification. I'm trying to get some intuition for regression boosting. At each iteration, adaptive boosting forces the next weak learner to focus more on incorrectly classified points. Intuitively, I can see why that should lead to a good classifier. I could be wrong, but L2 boosting doesn't seem to do anything like that. In L2 boosting, at each iteration, you're fitting a weak learner to the previous iteration's residuals. In a regression tree, when you're splitting a node, aren't you "fitting" the residuals from that node and that node's parents? In both cases, you're "fitting" unweighted residuals, so I don't understand why they're different. Maybe the the main advantage of L2 boosting over a single tree is that the former has many more regularization/bootstrap-ish options (e.g., randomly choosing subsets of features, the learning rate, the number of trees, individual tree depth, etc.)?


title: How to compare GCV diagnostic plots of a spline object

text: I am trying to analyse diagnostic plots of a spline object, I am using the package **Fields** from **R**. I am lost with the GCV function plot, since I can't find any guidelines of what to look for when comparing different models. Besides the value/degrees of freedom of the GCV function minimum for each model, can we also compare the GCV curve between models? This is an example of the different models and their correspondent diagnostic plots: (1) **_________________________** (2) **__________________________** (3) **____________________** ![enter image description here](http://i.stack.imgur.com/M8v7L.png) The summary of each fit:                (1)      Number of Observations:                745      Number of unique points:               745      Number of parameters in the null space 21      Parameters for fixed spatial drift     21      Effective degrees of freedom:          116.4      Residual degrees of freedom:           628.6      MLE sigma                              0.4153      GCV sigma                              0.4269      MLE rho                                0.02785      Smoothing parameter lambda             6.193           DETAILS ON SMOOTHING PARAMETER:      Method used:   GCV    Cost:  1        lambda       trA       GCV   GCV.one GCV.model      shat        6.1932  116.4326    0.2160    0.2160        NA    0.4269           (2)      Number of Observations:                734      Number of unique points:               734      Number of parameters in the null space 21      Parameters for fixed spatial drift     21      Effective degrees of freedom:          134.4      Residual degrees of freedom:           599.6      MLE sigma                              0.2788      GCV sigma                              0.2912      MLE rho                                0.02166      Smoothing parameter lambda             3.589           DETAILS ON SMOOTHING PARAMETER:      Method used:   GCV    Cost:  1        lambda       trA       GCV   GCV.one GCV.model      shat        3.5890  134.3993    0.1038    0.1038        NA    0.2912           (3)      Number of Observations:                716      Number of unique points:               716      Number of parameters in the null space 21      Parameters for fixed spatial drift     21      Effective degrees of freedom:          680.2      Residual degrees of freedom:           35.8      MLE sigma                              0.04486      GCV sigma                              0.0766      MLE rho                                11.6      Smoothing parameter lambda             0.0001734           DETAILS ON SMOOTHING PARAMETER:      Method used:   GCV    Cost:  1        lambda       trA       GCV   GCV.one GCV.model      shat     1.734e-04 6.802e+02 1.173e-01 1.173e-01        NA 7.660e-02      The lower values of lambda the smoother the surface (3), and the number of degrees of freedom is higher(3). Althought I know from the data that the observed values vary little for model (3), while for (1) and (2) I have a higher variation of observed values.


title: ANCOVA with repeated measures IV, violation of slopes

text: I am running a 2x2x5 mixed factorial ANCOVA, with one within-groups variable, and two between groups variable. I want to check for violation of homogeneity of regression slopes. My covatiate does not interact with my between groups variables. It does however interact with my within-groups variable. Is this still a violation of the assumption?


title: MLE and ordinal probit regression

text: Suppose we have a data set $X$. This data set consists of ordinal data (4 levels). To get estimates of the threshold coefficients and probit slope ($\beta_1, \dots, \beta_3$ and $\beta_4$ respectively), do most computational packages use maximum likelihood estimation? That is, given the data, MLE chooses the parameters that maximizes the probability of observing the data?


title: How can I calculate a confidence interval for the mean based on non-independent data?

text: In Chemical industrial, samples are often analyzed multiple times, e.g. 5 samples and each was analyzed 2 times with 10 data points in total. What would be the confidence interval for the mean estimated from these data? The degrees of freedom is somewhere from 4 to 10 depending on the correlation between repeats.


title: Predicting curve registration parameters in functional data analysis with noisy data

text: I have some data I'm studying where functional data analysis seems like a promising approach. But having never tackled FDA before, I'm having trouble wrapping my head around it. For background, I have Ramsay and Silverman's "Functional Data Analysis" and Ramsay, Hooker, and Graves "Functional Data Analysis with R and Matlab" but only got them recently, and am still very much the newbie with FDA. I am using R with package fda as the analysis software. My data are a sample of hundreds of young adult subjects, each of whom were measured annually by performing a task with a binary outcome over some number of trials. The outcome of interest is the success rate on the task. I'll include a short sample of data at the end of my question. Several parameters were collected on each subject, including height, weight, education level, amount of prior training on the task, and the results of an aptitude pre-test. The year-to-year change in success rate over time is the curve I'd like to model with FDA. There is significant variation across subjects in overall performance, but I am mostly interested in how a subject's success rate evolves over time rather than a subject's actual performance. The subjects generally follow a pattern of starting off with a low success rate, improving yearly until reaching a personal maximum, then declining with age. The typical subject enters the study at age 20 with a low success rate, improves until age 25, then declines until age 30, at which point no further data is collected. One of the goals is to characterize the "canonical" development curve for the population. FDA seems like a good fit for this aspect of the problem. Some of my challenges are:   * Subjects may reach their maximum success rate at different ages, and may improve/decline faster than others, so both phase variation and amplitude variation are potentially present. So it seems like curve registration will be an important part of the analysis.    * The number of trials in the test is relatively small each year, so the outcome data is noisy, and registering the data properly seems daunting since the noise in one subject can obscure the landmarks in the data.    * I am also interested in whether the other variables predict the registration characteristics of the subject's curve. i.e. Do college-educated subjects reach their maximum performance at an earlier age than others (phase), or do particularly tall subjects sustain their peak performance longer than short subjects do (amplitude)?   * We would like to be able to produce a customized curve for a subject. e.g. If we observed the tests through age 23, we would be able to predict how much improvement to expect, at what age the maximum performance would be reached, and how sharp the decline phase would be through age 30.  My questions include:   1. Is FDA an appropriate method for this problem?   2. Are there techniques in FDA for predicting a subject's phase and amplitude variation using other variables? What I've read so far seems to treat registration as a correction for noise rather than as containing features worth examination on its own.   3. What are some good ways to handle registration where the individual curves have a lot of noise?   4. Should I be able to produce individualized forecast curves that include subject-specific predicted phase/amplitude variation? What difficulties should I anticipate? Thank you for any and all suggestions. Sample data               Subject  Education    Ht   Wt   Training   Aptitude          A        College      72  200   0          Medium     B        High School  77  250   100        High     C        High School  68  160   50         Low          Subject Age Trials  Success Success%     A   20  15  3   20%     A   21  18  5   28%     A   22  30  7   23%     A   23  28  8   29%     A   24  32  13  41%     A   25  8   2   25%     A   26  20  8   40%     A   27  40  11  28%     A   28  33  10  30%     A   29  18  5   28%     A   30  10  2   20%          B   20  24  4   17%     B   21  27  5   19%     B   22  30  8   27%     B   23  33  2   6%     B   24  41  8   20%     B   25  39  5   13%     B   26  39  5   13%          C   24  13  4   31%     C   25  19  6   32%     C   26  18  5   28%     C   27  23  6   26%     C   28  16  6   38%     C   29  9   3   33%


title: how to derive the distribution of variable x conditional on x+y

text: The question is from a typical example for E-M algorithm. Let's say $(y_1,y_2,y_3)$ $\sim$ $\text{multinomial}(n;p_1,p_2,p_3)$, where $p_1+p_2+p_3=1$. How can we derive the conditional distribution of $y_2$ given $y_2+y_3=n$? The answer is $y_2|y_2+y_3 \sim \text{binomial}(n, p_2/(p_2+p_3)$). Any idea on how to derive this rigorously?


title: Why is homogeneity of variance so important?

text: I do not know why testing homogeneity of variance is so important. What are the examples that require homogeneity of variance?


title: Structural panel VAR estimation

text: I am running a structural Panel vector autoregressive model on a panel of 13 countries over the period 1970-2012. I'm having problems in implementing the model. Does anyone have estimation programs (from Matlab, Eviews, RATS or any time series package tool) so that I can use them to estimate my model? Thanks in advance.


title: Suitable distributions for the SDs of a number of groups (or participants)?

text: When analyzing data, using MLE or Bayesian methods, one needs to assume a distribution for the data. For continuous data the are a number of distributions that are often considered, for example, the normal distribution, the t-distribution, the log-normal, etc. **When analyzing the distribution of the SDs of a number of groups (or participants) what distributions could be appropriate?** Of course this is dependent on the data, and varies from case to case, but what could be some reasonable distributions to try?


title: Clustered binary data

text: I have a binary outcome (success \ failure). I have 20 subjects, half of them gave 2 samples, the others just one, so I have 30 data points. All data points without any exceptions were success (1). So my point estimate of the success rate is clearly 100%. Now I want to calculate a CI for the success rate, mainly to see what the lower limit is, I want to be able to say that with 95% confidence the success rate is higher than....(80%, 85%, whatever comes up). The problem is, as you can see, is the clustering, I can't use n = 30 because I then ignore the correlation, I feel it is a waste to use n = 20. Is there a way to calculate the variance of a single sample proportion while taking into account the cluster ? Thanks !


title: Quantile function in Stata

text: I'm new to Stata and didn't find the answer in the help, navigating the menus, or online so far. I would just like to know how to access the quantile function of a distribution. For example in R, if I want to know the 0.95 point of the $\chi^2(1)$ distribution, I do:               > qchisq(0.95,1)     [1] 3.841459      Is there an equivalent command in Stata ?


title: Voting system that uses confidence for each voter

text: Let's say, we have simple "yes/no" question that we want to know answer to. And there are N people "voting" for correct answer. Every voter has a history - list of 1's and 0's, showing whether they were right or wrong about this kind of questions in the past. If we assume history as a binomial distribution, we can find voters' mean performance on such questions, their variation, CI and any other kind of confidence metrics. Basically, my question is: how to incorporate **confidence information** into **voting system**? For example, if we consider only mean performance of each voter, then we can construct simple weighted voting system: $$result = sign(\sum_{v \in voters}\mu_v \times (-1)^{1-vote})$$ That is, we can just sum voters' weights multiplied either by $+1$ (for "yes") or by $-1$ (for "no"). It makes sense: if voter 1 has average of correct answers equal to $.9$, and voter 2 has only $.8$, than, probably, 1st person's vote should be considered as more important. On other hand, if 1st person have answered only 10 questions of this kind, and 2nd person have answered 1000 such questions, we are much more confident about 2nd person's skill level than about those of the 1st - it's just possible that 1st person was lucky, and after 10 relatively successful answers he will continue with much worse results. So, more precise question may sound like this: is there statistical metric that incorporates both - **strength** and **confidence** about some parameter?


title: Distribution of the convolution of squared normal and chi-squared variables?

text: the following problem came up recently while analyzing data. If the random variable X follows a normal distribution and Y follows a $\chi^2_n$ distribution (with n dof), how is $Z = X^2 + Y^2$ distributed? Up to now I came up with the pdf of $Y^2$: \begin{eqnarray} \psi^2_n(x) &=& \frac{\partial F(\sqrt{x})}{\partial x} \\ &=& \left( \int_0^{\sqrt{x}} \frac{t^{n/2-1}\cdot e^{-t/2}}{2^{n/2}\Gamma(n/2)} \mathrm{d}t \right)^\prime_x \\ &=& \frac{1}{2^{n/2}\Gamma(n/2)} \cdot \left( \sqrt{x} \right)^{n/2-1} \cdot e^{-\sqrt{x}/2} \cdot \left( \sqrt{x} \right)^\prime_x \\ &=& \frac{1}{2^{n/2-1}\Gamma(n/2)} \cdot x^{n/4-1} \cdot e^{-\sqrt{x}/2} \end{eqnarray} as well as some simplifications for the convolution integral ($X^2$ has the pdf $\chi^2_m$ with m dof): \begin{eqnarray} K_{mn}(t) &:=& ( \chi^2_m \ast \psi^2_n )(t) \\ &=& \int_0^t \chi^2_m(x) \cdot \psi^2_n(t-x) \mathrm{d}x \\ &=& \left( 2^{\frac{(n+m)}{2}+1} \Gamma(\frac{m}{2}) \Gamma(\frac{n}{2}) \right)^{-1} \cdot \int_0^t (t-x)^{\frac{n}{4}-1} \cdot x^{\frac{m}{2}-1} \cdot \exp(-(\sqrt{t-x}+x)/2) \mathrm{d}x \end{eqnarray} Does someone see a good way of calculating this integral for any real t or does it have to be computed numerically? Or am I missing a much simpler solution? Thanks!


title: Error exponent in hypothesis testing

text: In hypothesis testing, one must decide between two probability distributions $P_1(x)$ and $P_2(x)$ on a finite set $X$, after observing $n$ i.i.d. samples $x_1,...,x_n$ drawn from the unknown distribution. Let $A_n\subseteq X^n$ denote the chosen acceptance region for $P_1$. The error probabilities of type I and II can be expressed thus $$ \alpha_n = P^n_1(A^c_n)$$ $$ \beta_n = P^n_2(A_n)$$ (Cover & Thomas, Ch. 11 is an excellent reference for the definitions and facts mentioned in this post). Assume we have chosen the acceptance regions $A_n$'s ($n\geq 1$), so that both error probabilities approach zero as the number of observations grows: $\alpha_n\rightarrow 0$ and $\beta_n\rightarrow 0$ as $n\rightarrow \infty$. Stein's Lemma tells us that the maximum rate of deacrease of both error probabilities is determined, to the first order of the exponent, by the the KL-distance between the given distributions. More precisely $$ -\frac 1 n \log \alpha_n \rightarrow D(P_2||P_1)\tag{1}$$ $$ -\frac 1 n \log \beta_n \rightarrow D(P_1||P_2)\tag{2}$$ Now, consider the Bayesian version of the hypothesis testing problem. In this case, $P_1$ and $P_2$ are given prior probabilities $\pi_1$ and $\pi_2$, respectively, and the error probability is obtained by weighting $\alpha_n$ and $\beta_n$: $$ e_n = \pi_1\alpha_n + \pi_2\beta_n.\tag{3}$$ In this case, the optimal exponent for $e_n$ is given by Chernoff distance between the given distributions: $$ -\frac 1 n \log e_n \rightarrow C(P_1,P_2).$$ **Question** : what is wrong in the reasoning below? (Disclaimer: I'm _not_ trying to be fully formal/detailed here). By (3), the decrease rate of $e_n$ is the minimum deacrease rate of $\alpha_n$ and $\beta_n$: $$ \lim -\frac 1 n \log e_n = \min\{\lim -\frac 1 n \log \alpha_n, \lim -\frac 1 n \log \beta_n\}$$. Since $e_n\rightarrow 0$, one must have both $\alpha_n\rightarrow 0$ and $\beta_n\rightarrow 0$ as $n\rightarrow \infty$. So, by the previous considerations on Stein's Lemma, and (1) and (2), one would get $$ \lim -\frac 1 n \log e_n = \min\{D(P_1||P_2), \,\,D(P_2||P_1)\}$$ which is quite different from $C(P_1,P_2)$. **EDIT** : I realize that now that (1) and (2) cannot hold simultaneously, for the same regions $A_n$'s, so this must be the bug in the reasoning. What one can infer through a similar reasoning is just, I think, $$C(P_1,P_2)\leq \min\{D(P_1||P_2), \,\,D(P_2||P_1)\}.$$


title: Growth rate vs. Absolute correlations

text: Conceptually, what is the reason why the correlation between growth rates of A and B would be different from correlation between actual A and B? Under what circumstances would the growth rate correlation be higher? How about lower? How about equal? Thanks.


title: Classification Tree Analysis - Assessment of tree explanatory power (R-square?) using the party package in R

text: I have produced a model using the `ctree` function in R, and want to know whether this tree is actually explaining my data well. I am trying to explain the presence or absence of landscape disturbances using environmental factors. `rpart` gives an $R^2$ value, but since my data are not well-balanced (PRESENCE contains a relatively small number of 1s), I believe the `ctree` model will work better. A (simplified) version of my model is:               fit <- ctree (PRESENCE ~ A + B + C + D + E + F, data=dat);      where PRESENCE is a factor containing values of (0,1), five of the independent variables (denoted by the letters A-F) are numerical values (environmental data), and the remaining variable is categorical. I believe there is a way to make an $R^2$ value using the `predict` function (as in Getting R square value from ctree), but I'm fairly new to R and haven't had any success trying to code this. Any help would be much appreciated.


title: A statistical model for a sample of independent networks

text: Based on the lack of responses to my previous network question, perhaps this is not quite the place to ask this question, but I'll give it a try. I am planning a series of studies that involve small groups of people. In a typical study, the participants will meet in groups of 4 to have a discussion about a particular topic. As a group, the (student) participants will reach a decision about the topic (e.g., whether to support or oppose a proposed change in university policy). The discussion will be videotaped, and we will use transcripts of the discussion to obtain information about who is talking to whom about what topics. In addition, the participants will complete ratings of themselves and each other (e.g., What is your opinion about the proposed change? What do you think is this person's opinion of the proposed change?). I will also have a variety of data about the individual characteristics of all the participants. Thus, each group from each study will give me a 4-node network with directed, valued edges and a variety of attributes about each node. In a typical study, I might have 20 - 30 of these small groups (i.e., small, 4-node networks). In all the studies, I will be interested in group-level, relationship-level, and individual-level outcomes. I have provided below a sampling of the kinds of questions I would want to ask from these data:   1. Which groups express the most agreement in their discussions? Which groups express the most disagreement? Is group composition related to agreement / disagreement and to the final decision outcome?   2. Which pairs of people like each other during the discussion? Which pairs of people dislike each other? How are the pairwise patterns of liking related to the final decision outcome?   3. Who has the most influence on the course of the discussion? Who has influence over the final decision of each group? Are there any individual difference characteristics that are related to influence over the discussion and / or final decision? And so on. What I am looking is the proper statistical model to use to investigate samples of networks, rather than individual networks. I have read a bit about single-network methods (like exponential random graph models, for which there is a nice R package available), but these methods do not seem appropriate, since I am dealing with a set of independent networks, rather than a single network. In addition, a plain linear mixed model does not seem appropriate because I am collecting explicitly relational data. Finally, although a method like the social relations model seems appropriate for this situation at first blush, the social relations model seems to only partition the variance in round-robin ratings into perceiver, target, and relationship sources rather than allowing me to, for example, relate the patterns within a set of networks to overall network outcomes. Could anyone offer some advice about which statistical model might be appropriate for my situation? Any readings suggestions and / or software recommendations would be greatly appreciated. (FYI, I am a proficient R user, so R package recommendations would be especially appreciated). * * * **Edit:** At the request of Alex Williams, I have posted some example data here in csv format to give everyone a concrete example of the type of data I'm working with. In the sample data, participants, identified by `participant_id`, are each in 4-person groups, identified by `group_id`. The participants discuss a proposal and decide as a group whether the support or oppose the proposal. `group_decision` is the result of the group discussion. The participants also rate their own attitude toward the proposal (`pers_att`), their perceptions of the attitudes of the other group members (`p1_att` through `p4_att`; `person_id` tracks who `p1` through `p4` are within each group), and their own enjoyment of the discussion. In the sample data, I might be interested in the following sorts of questions:   1. Are personal attitudes related to group-level decisions?   2. Are personal attitudes related to personal enjoyment of the group discussion?   3. Were people accurate in their ratings of other peoples' attitudes?   4. Were the ratings of others related to others' enjoyment of the discussion?   5. Does disagreement within a group (quantified, for example, by the difference between the minimum and maximum attitude ratings) relate to group-level decisions?   6. Does disagreement within a group cause people to enjoy the discussion less?   7. Did people who enjoyed the discussion more have more influence over the outcome?


title: Causes of sparseness in data

text: Can data sparseness appear due to either high sample size or high dimension? How different are the situations in the two cases? Thanks!


title: Have I applied the delta method correctly?

text: I have UMVUE $$\tilde\theta = \frac{(n-1)(U-n)}{(U-1)(U-2)}$$ for $P(Y=2)=\theta(1-\theta)$ where $U=\sum_{i=1}^n Y_i$ $Y_i \sim \text{geometric} (\theta)$ I am using delta method to find the variance of $\tilde\theta$ So far, I have defined: $$g(y)=\frac{(n-1)(y-n)}{(y-1)(y-2)}$$ $$g'(y) = \frac{(n-1)[n(2y-3)-y^2+2]}{(y-2)^2(y-1)}$$ Replacing $y$ with $\theta$ for $g'(\theta)$ Since $\tilde \theta$ is unbiased estimator, $E[\tilde \theta] = \theta$ $Var(\sum Y_i)=n^2(1-\theta)/\theta^2 $ [ This is the step where I am confused, do I have to take the variance of one of the $Y_i$'s or sum of $Y_i'$s So, $Var(\tilde \theta) = [g'(\theta)]^2Var (\sum Y_i)$ Is this approach is correct? Any pointers would be helpful.


title: Converting an SPSS file to a SAS file?

text: (I'm asking this question for a friend, honest...) > Is there an easy way to convert from an SPSS file to a SAS file, which > preserves the formats AND labels? Saving as a POR file gets me the labels (I > think) but not the POR file. I tried to save to a SAS7dat file but it didn't > work. Thanks,


title: How to solve a case of unbalanced repeated measures?

text: I have animals, that could be virgin or mated (reproductive state is the fixed factor), which I've stimulated sequentially with 4 different doses of an odour (doses are the repeated measures, the same animal was blown with 4 increasing doses of the same odorant). Then, I measure the neuronal response (variable: number of spikes) of each animal to each dose of the odorant. This might be a typical case or repeated measures, however I have some missing values for the doses, I have not all the doses completed for some animals. For example, for animal 1, I missed recording 1 out of the 4 doses. what can I do? I have two statistical packages: SPSS 16 or Statistical. thanks for your help!


title: Calculate percentage difference for two sets of points

text: As the title says, I'd like to calculate the percentage difference for two sets of points. For example, suppose I have $S_{1}=\{(1,x_{1}),(2,x_{2}),(3,x_{3})\}$ and $S_{2}=\{(1,y_{1}),(2,y_{2}),(3,y_{3})\}$. How can I know the difference in percentage between both sets of data. What is the correct way to do that? Is that kind of assessment meaningful to establish to which degree of precision a set of data is preferred over the other? In my particular case, $S_{1}$ is simply a set of numerical results obtained by DSMC and $S_{2}$ was obtained by a theoretical result. I'd like to quantify how much difference exist between each other in order to establish when it is convenient to use one or the other. By "difference in percentage" I mean percent difference. Hopefully that clarifies a bit the question. **UPDATE:** Another way to formulate my question would be: How can I arrive to conclusions such as "The results from experiment A are inaccurate by 10% with respect to experiment B", when experiment A and B are a set of values.


title: Problem with GBM model

text: Alright, hopefully third time is the charm. I'm basically trying to build a predictive model with R and gbm. For various reasons, I can't explicitly state exactly what I'm doing. Basically I have a response variable that averages around zero. It has high volatility, is right skewed, and has very high kurtosis. I have a bunch of predictors I think contain information about the response. Some of them are numerical, others factors. I have partitioned available data into three sets 1, 2, 3. What I've done is trained gbm on set 1 with parameters as follows: interaction depth of 10, shrinkage of .001, 3 fold cross validation, about 200000 iterations (maybe overkill), out of bag estimation 50%, and training % at 80%. I've then tested how well the model predicts the responses in set 2 and it's significantly better than chance. The predictors are distributed similarly to the response (similar averages), but the responses have a higher stdev and kurtosis. I then proceed to update the model by training on set 1 + set 2 with the same parameters. I then test on set 3. The predictors now are extremely different from the actual responses (average of 25 for the predictor versus close to 0 for responses). I'm not sure what's causing this. I haven't changed anything so my initial thought was the data in set 2, where the range of the response variable is larger than set 1. How do I narrow down what might be causing this discrepancy with the predictions for set 3? ## I realize this is somewhat open-ended, so please let me know what other information would be helpful in answering this question. I'm not sure what exactly you guys might find useful. answers to comments the data is a time series so just took three time period of about 2 years, 1 year, and 3 months.


title: Number of observations needed to perform pca

text: I read in this paper (page 3) comparing pca to factor analysis that both methods need a number of observations about 5 times the number of variables. Why? and how would you reduce the number of variables if you have few observations only?


title: Bayesian statistics tutorial

text: I am trying to get upto speed in Bayesian Statistics. I have a little bit of stats background (STAT 101) but not too much - I think I can understand prior, posterior, and likelihood :D. I don't want to read a Bayesian textbook just yet. I'd prefer to read from a source (website preferred) that will ramp me up quickly. Something like this, but that has more details. Any advice?


title: Redefining and merging two tables in R or SAS

text: I have two tables (matrix) of same dimensions, one contains correlation coefficients and other with p values. I want to combine them into one table. For example let's say I have correlation coefficient between variable A1 and A2 of 0.75 in table 1 and p value of 0.045 in table 2. Now in my combined table 3, I want to use: condition1 for table 1: if a coefficient value in a cell of table 1 is less than 0.4 then "+", 0.4 <= coefficient <0.7 then "++" else "+++", condition2 for table 2: if a pvalue in a cell of table 2 is less than 0.01 then "+++", 0.01 <= pvalue < .05 then "++" else "+". Thus corresponding cell value for A1 and A2 in table 3 should look like: +++/++ where "+++" correspond to table 1 value of 0.75 and ++ correspond to table 2 p value of 0.045 and "/" is just a separator. I would like to do this either is SAS or R. Thanks in advance.


title: Manually calculated $R^2$ doesn't match up with randomForest() $R^2$ for testing new data

text: I know this is a fairly specific `R` question, but I may be thinking about proportion variance explained, $R^2$, incorrectly. Here goes. I'm trying to use the `R` package `randomForest`. I have some training data and testing data. When I fit a random forest model, the `randomForest` function allows you to input new testing data to test. It then tells you the percentage of variance explained in this new data. When I look at this, I get one number. When I use the `predict()` function to predict the outcome value of the testing data based on the model fit from the training data, and I take the squared correlation coefficient between these values and the _actual_ outcome values for the testing data, I get a different number. _These values don't match up_. Here's some `R` code to demonstrate the problem.               # use the built in iris data     data(iris)          #load the randomForest library     library(randomForest)          # split the data into training and testing sets     index <- 1:nrow(iris)     trainindex <- sample(index, trunc(length(index)/2))     trainset <- iris[trainindex, ]     testset <- iris[-trainindex, ]          # fit a model to the training set (column 1, Sepal.Length, will be the outcome)     set.seed(42)     model <- randomForest(x=trainset[ ,-1],y=trainset[ ,1])          # predict values for the testing set (the first column is the outcome, leave it out)     predicted <- predict(model, testset[ ,-1])          # what's the squared correlation coefficient between predicted and actual values?     cor(predicted, testset[, 1])^2          # now, refit the model using built-in x.test and y.test     set.seed(42)     randomForest(x=trainset[ ,-1], y=trainset[ ,1], xtest=testset[ ,-1], ytest=testset[ ,1])      Thanks for any help you might be willing to lend.


title: Exogenous variables in VECM

text: I found the following posts interesting and I was wondering if any of you guys know of good academic papers that describe methods/relationships of exogenous variables in VECM models. If so could you kindly point them out to me as I am very interested in learning. Thank you. Finding coefficients for VECM + exogenous variables Lagged Exogenous Variables in VECM with R


title: C++ libraries for statistical computing

text: I've got a particular MCMC algorithm which I would like to port to C/C++. Much of the expensive computation is in C already via Cython, but I want to have the whole sampler written in a compiled language so that I can just write wrappers for Python/R/Matlab/whatever. After poking around I'm leaning towards C++. A couple of relevant libraries I know of are Armadillo (http://arma.sourceforge.net/) and Scythe (http://scythe.wustl.edu/). Both try to emulate some aspects of R/Matlab to ease the learning curve, which I like a lot. Scythe squares a little better with what I want to do I think. In particular, its RNG includes a lot of distributions where Armadillo only has uniform/normal, which is inconvenient. Armadillo seems to be under pretty active development while Scythe saw its last release in 2007. So what I'm wondering is if anyone has experience with these libraries -- or others I have almost surely missed -- and if so, whether there is anything to recommend one over the others for a statistician very familiar with Python/R/Matlab but less so with compiled languages (not completely ignorant, but not exactly proficient...).


title: How to differentiate two subgroups from a histogram?

text: I have a set of samples in which I assume there are 2 definite subsets in it. I plotted their values in a histogram and found that there are two distinct modes as shown in the figure below. My question is how do I differentiate two groups. i.e how do I choose a value that differentiates the two subsets? ![enter image description here](http://i.stack.imgur.com/wbIvt.png)


title: Classifying spectra

text: I am currently trying to do the following in R: I have thousands of **measured spectra** (x,y; see below). Each spectra has one or two peaks. Also I have sets of **"training" spectra** obtained in more controlled conditions and I would like to know which of my training spectra has the closest match to the measured spectra!? I was thinking that some sort of pattern recognition would be useful but I know too little to make an informed choice as this is a bit outside of my usual work-area   * What is the most promising way/function in R to do this kind of pattern recognition I want?   * In case pattern recognition (like PCA) is not the most promising way, what other options are there? I am looking for sample bits of code or literature dealing with this kind of data analysis. ![enter image description here](http://i.stack.imgur.com/5mpqU.png) **EDIT** The peak position will most probably always be the same, however the laser used to record the spectra is temperature controlled and slight variations are possible. The intensity will change depending on experimental conditions. The two peaks sould be treated as independet peaks.


title: Standard error clustering in R (either manually or in plm)

text: I am trying to understand standard error "clustering" and how to execute in R (it is trivial in Stata). In R I have been unsuccessful using either `plm` or writing my own function. I'll use the `diamonds` data from the `ggplot2` package. I can do fixed effects with either dummy variables               > library(plyr)     > library(ggplot2)     > library(lmtest)     > library(sandwich)     > # with dummies to create fixed effects     > fe.lsdv <- lm(price ~ carat + factor(cut) + 0, data = diamonds)     > ct.lsdv <- coeftest(fe.lsdv, vcov. = vcovHC)     > ct.lsdv          t test of coefficients:                                Estimate Std. Error  t value  Pr(>|t|)         carat                 7871.082     24.892  316.207 < 2.2e-16 ***     factor(cut)Fair      -3875.470     51.190  -75.707 < 2.2e-16 ***     factor(cut)Good      -2755.138     26.570 -103.692 < 2.2e-16 ***     factor(cut)Very Good -2365.334     20.548 -115.111 < 2.2e-16 ***     factor(cut)Premium   -2436.393     21.172 -115.075 < 2.2e-16 ***     factor(cut)Ideal     -2074.546     16.092 -128.920 < 2.2e-16 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1       or by de-meaning both left- and right-hand sides (no time invariant regressors here) and correcting degrees of freedom.               > # by demeaning with degrees of freedom correction     > diamonds <- ddply(diamonds, .(cut), transform, price.dm = price - mean(price), carat.dm = carat  .... [TRUNCATED]      > fe.dm <- lm(price.dm ~ carat.dm + 0, data = diamonds)     > ct.dm <- coeftest(fe.dm, vcov. = vcovHC, df = nrow(diamonds) - 1 - 5)     > ct.dm          t test of coefficients:                   Estimate Std. Error t value  Pr(>|t|)         carat.dm 7871.082     24.888  316.26 < 2.2e-16 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1       I can't replicate these results with `plm`, because I don't have a "time" index (i.e., this isn't really a panel, just clusters that could have a common bias in their error terms).               > plm.temp <- plm(price ~ carat, data = diamonds, index = "cut")     duplicate couples (time-id)     Error in pdim.default(index[[1]], index[[2]]) :       I also tried to code my own covariance matrix with clustered standard error using Stata's explanation of their `cluster` option (explained here), which is to solve $$\hat V_{cluster} = (X'X)^{-1} \left( \sum_{j=1}^{n_c} u_j'u_j \right) (X'X)^{-1}$$ where $u_j = \sum_{cluster~j} e_i * x_i$, $n_c$ si the number of clusters, $e_i$ is the residual for the $i^{th}$ observation and $x_i$ is the row vector of predictors, including the constant (this also appears as equation (7.22) in Wooldridge's _Cross Section and Panel Data_ ). But the following code gives very large covariance matrices. Are these very large values given the small number of clusters I have? Given that I can't get `plm` to do clusters on one factor, I'm not sure how to benchmark my code.               > # with cluster robust se     > lm.temp <- lm(price ~ carat + factor(cut) + 0, data = diamonds)     >      > # using the model that Stata uses     > stata.clustering <- function(x, clu, res) {     +     x <- as.matrix(x)     +     clu <- as.vector(clu)     +     res <- as.vector(res)     +     fac <- unique(clu)     +     num.fac <- length(fac)     +     num.reg <- ncol(x)     +     u <- matrix(NA, nrow = num.fac, ncol = num.reg)     +     meat <- matrix(NA, nrow = num.reg, ncol = num.reg)     +          +     # outer terms (X'X)^-1     +     outer <- solve(t(x) %*% x)     +      +     # inner term sum_j u_j'u_j where u_j = sum_i e_i * x_i     +     for (i in seq(num.fac)) {     +         index.loop <- clu == fac[i]     +         res.loop <- res[index.loop]     +         x.loop <- x[clu == fac[i], ]     +         u[i, ] <- as.vector(colSums(res.loop * x.loop))     +     }     +     inner <- t(u) %*% u     +      +     #      +     V <- outer %*% inner %*% outer     +     return(V)     + }     > x.temp <- data.frame(const = 1, diamonds[, "carat"])     > summary(lm.temp)          Call:     lm(formula = price ~ carat + factor(cut) + 0, data = diamonds)          Residuals:          Min       1Q   Median       3Q      Max      -17540.7   -791.6    -37.6    522.1  12721.4           Coefficients:                          Estimate Std. Error t value Pr(>|t|)         carat                 7871.08      13.98   563.0   <2e-16 ***     factor(cut)Fair      -3875.47      40.41   -95.9   <2e-16 ***     factor(cut)Good      -2755.14      24.63  -111.9   <2e-16 ***     factor(cut)Very Good -2365.33      17.78  -133.0   <2e-16 ***     factor(cut)Premium   -2436.39      17.92  -136.0   <2e-16 ***     factor(cut)Ideal     -2074.55      14.23  -145.8   <2e-16 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 1511 on 53934 degrees of freedom     Multiple R-squared: 0.9272, Adjusted R-squared: 0.9272      F-statistic: 1.145e+05 on 6 and 53934 DF,  p-value: < 2.2e-16           > stata.clustering(x = x.temp, clu = diamonds$cut, res = lm.temp$residuals)                             const diamonds....carat..     const                11352.64           -14227.44     diamonds....carat.. -14227.44            17830.22      Can this be done in R? It is a fairly common technique in econometrics (there's a brief tutorial in this lecture), but I can't figure it out in R. Thanks!


title: How to express a Poisson regression as an equation

text: I have a Poisson regression where I regress a counts variable (called $H(x)$) on one covariate (called $x$) and one factor (categorical) variable (called $s$), which can take one of three values recoded to integers 1 through 3. I am including an offset variable called $P(x)$ and an intercept called $B_0$. I am trying to write the regression equation. I believe I have correctly expressed the dependent, offset, intercept, and $B_1$ terms, but am not sure about how to express the categorical variable (the $s$) term. Since it's a categorical variable that takes an integer value in the range 1 through 3, I would like to know if can I write it like this (in bold below): $$\log H(x) = \log P(x) + B0 + B1 \cdot x + B2\cdot s$$ or if instead I should express each possible value in the equation (since the regression results give one beta parameter estimate for each of the three values) like this: $$\log H(x) = \log P(x) + B0 + B1 \cdot x + B2 \cdot 1 + B3 \cdot 2 + B4 \cdot 3$$ The value for this categorical factor coded as "3" is the reference variable and so the parameter estimate is set by the regression procedure in SPSS to 0.


title: Understanding MAD Results

text: I'm using R to calculate the median absolute deviation for a few distributions, but some of the values I'm calculating do not seem realistic at all. I have the following distribution:               x <-   [1]     NA     NA     NA -0.003 -0.009  0.004 -0.001 -0.001 -0.003  0.001 -0.002  0.000 -0.003  0.000  0.006 -0.011 -0.003      [18]  0.002 -0.007 -0.002  0.006 -0.005  0.000  0.008  0.001  0.009 -0.002  0.001  0.001  0.002  0.003     NA     NA  0.001      [35]     NA  0.005 -0.002  0.003  0.016  0.007 -0.003 -0.017  0.000 -0.013  0.000  0.002  0.002  0.000     NA  0.000  0.000      [52]  0.000  0.000  0.004 -0.001  0.000 -0.002 -0.003 -0.007 -0.001 -0.001  0.000 -0.002  0.001  0.003  0.000 -0.011 -0.002      [69] -0.003  0.004 -0.007     NA -0.009  0.005 -0.001  0.001 -0.001  0.001 -0.001  0.006  0.002 -0.006  0.002 -0.002  0.004      [86]  0.006  0.001  0.000  0.002 -0.002  0.007  0.004  0.003  0.004  0.005 -0.005  0.003 -0.003  0.002  0.004  0.003 -0.002     [103] -0.002  0.001  0.002  0.000  0.000  0.003 -0.001  0.004  0.001  0.001  0.005 -0.001     NA -0.005  0.000 -0.002 -0.004     [120]  0.004     NA  0.007  0.000  0.002  0.003 -0.006 -0.002  0.000 -0.002 -0.001 -0.001 -0.001 -0.006 -0.001 -0.001 -0.008     [137]  0.000  0.003  0.001  0.001 -0.001  0.000  0.011 -0.017     NA     NA     NA      Then I used the following code to generate my MAD value:               MADx <- mad(x, center = median(x, na.rm = TRUE), constant = (1/(quantile(x, probs=0.75, na.rm = TRUE, names = FALSE, type = 1))), na.rm = TRUE, low = FALSE, high = FALSE)      I get a value of 1 when doing this, which seems unrealistic because the values I have are much less than 1. I used the quantile function to get the 75th quantile of the distribution.


title: Weights argument in lm and lme very different in R- am I using them correctly?

text: So, it seems to me that the weights function in lm gives observations more weight the larger the associated observation's 'weight' value, while the lme function in lme does precisely the opposite. This can be verified with a simple simulation.               #make 3 vectors- c is used as an uninformative random effect for the lme model     a<-c(1:10)     b<-c(2,4,6,8,10,100,14,16,18,20)     c<-c(1,1,1,1,1,1,1,1,1,1)      If you were now to run a model where you weight the observations based on the inverse of the dependent variable in lm, you can only generate the exact same result in nlme if you weight by just the dependent variable, without taking the inverse.               summary(lm(b~a,weights=1/b))     summary(lme(b~a,random=~1|c,weights=~b))      You can flip this and see the converse is true- specifying weights=b in lm requires weights=1/b to get a matching lme result. So, I understand this much, I just want validation on one thing and to ask a question about another.   1. If I want to weight my data based on the inverse of the dependent variable, is it fine to just code weights=~(dependent variable) within lme?   2. Why is lme written to handle weights completely differently than lm? What is the purpose of this other than to generate confusion?  Any insight would be appreciated!


title: Estimating sample size

text: I am working with an investigator planning a study to validate a device to measure blood pressure (BP). The investigator states a difference within +/- 5 mmHg, when the two devices are used by the same participant, as "similar". I obtained the following summary statistics (slightly modified) from a similar study which used the ANSI/AAMI SP10A-2002 protocol: n = 100 participants   mean +/- SD systolic BP = -0.70 +/- 5.00   mean difference +/- diastolic BP = 0.20 +/- 5.20 How many participants would the investigator need to recruit to validate the device considering the typical 80% power and alpha level of 0.05? I am using R (specifically power.t.test). I assume I want to focus on diastolic BP due to the greater variation. I defined the null and alternative hypotheses as follows: Null Hypothesis: The average difference between the two measurements is 5 mmHg. Alternative Hypothesis: The average difference between the two measurements are no greater than 5 mmHg (one-sided hypothesis). But apparently my thinking is off... is there a piece of information I am lacking? Thank you in advance!


title: Interpretation of feature selection task

text: So I am given the following question > Data set sample5.txt has a 20-dimensional input $x$ in $\mathbb{R}^{20}$ but > we suspect that many of these are actually irrelevant. Could you model the > function $y = f(x)$ while - at the same time - figuring which dimensions > contribute to the output? So it is a feature selection task - I understand that. But I'm sort of confused by the > at the same time part. I know many feature selection algorithms but they do not actually produce models for the data, they just produce decisions regarding which features are important and which are not. Conversely, a model (alone) doesn't really give much information regarding which features are important and which are not. Perhaps you could do simple linear regression and then select features based on the weights (but I have never heard of anyone doing this). Or do you think that I am over-analyzing the question and what I should do is simply do feature selection first, then create the model?


title: Missing data: what's the best approach?

text: In a representative sample of country population I have very few missing data, around 3%. But when I checked the missing data among communities, I found that one of them has almost 30% of data lost. That's consistent in all ages. Should I try to impute the data, remove that community from the analysis, or keep all the communities but warn the reader about this possible bias?


title: Comparing F-values of covariates in R

text: I'm new here and have a question regarding ANOVA in R. I have an ANOVA table like this from running `anova(model)` in R, where `model` is a multiple linear regression model built with the `lm` command:               Analysis of Variance Table          Response: log(price)                    Df Sum Sq Mean Sq   F value    Pr(>F)         mileage.residual    1  10.79   10.79   733.329 < 2.2e-16 ***     model               5 138.31   27.66  1880.427 < 2.2e-16 ***     body_type           4 108.87   27.22  1850.269 < 2.2e-16 ***     age                 1 654.08  654.08 44463.584 < 2.2e-16 ***     transmission        1   3.87    3.87   263.416 < 2.2e-16 ***     four_wd             1   0.86    0.86    58.718 2.353e-14 ***     nav                 1   1.41    1.41    96.168 < 2.2e-16 ***     fuel                2   2.35    1.17    79.809 < 2.2e-16 ***     age:transmission    1   0.79    0.79    53.827 2.719e-13 ***     Residuals        3429  50.44    0.01                             ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1      From reading this I have three questions:     1. How is the the model for MS(Model) chosen when calculating the different F-values for the independent variables? Is it the model with only that variable? This leads to the second question.     2. Are the different F-values comparable? Does it make sense to draw the conclusion that `age` is by far the strongest predictor among the considered covariates and that `model` and `body_type` are of about the same strength?      3. If I train a new model, but with different data (in this case representing a different car model), does it make sense to compare the Mean Sq(Residuals) as a measure across different R-models of how well the model describes the different car models. That is, if this R-model is first trained with car model A's data and has a MSE of 0.01 (as in the table) and the same R-model is then trained with car model B's data and then has a MSE of 0.1, can I conclude that the R-model is a much better fit for model A than B? This all seems like reasonable things, but I just want to be sure to not jump to any unjustified conclusions. For completeness a summary of the model is shown below:               Call:     lm(formula = log(price) ~ mileage.residual + model + body_type +          age * transmission + four_wd + nav + fuel, data = sample.data)          Residuals:          Min       1Q   Median       3Q      Max      -0.74556 -0.07250  0.00126  0.07135  0.46688           Coefficients:                                 Estimate Std. Error  t value Pr(>|t|)         (Intercept)               12.5718774  0.0104270 1205.710  < 2e-16 ***     mileage.residual          -0.0318929  0.0006863  -46.471  < 2e-16 ***     model320                   0.0827985  0.0063347   13.071  < 2e-16 ***     model325                   0.2139785  0.0088962   24.053  < 2e-16 ***     model328                   0.2689928  0.0137284   19.594  < 2e-16 ***     model330                   0.2956210  0.0105510   28.018  < 2e-16 ***     model335                   0.4724337  0.0114792   41.156  < 2e-16 ***     body_typecab               0.3254419  0.0102776   31.665  < 2e-16 ***     body_typecoupe             0.1200356  0.0087123   13.778  < 2e-16 ***     body_typehatchback         0.1082495  0.0096499   11.218  < 2e-16 ***     body_typestation_wagon     0.0165039  0.0054116    3.050  0.00231 **      age                       -0.1588759  0.0013296 -119.491  < 2e-16 ***     transmissionautomatic      0.1061707  0.0076932   13.801  < 2e-16 ***     four_wdTRUE                0.0478496  0.0067139    7.127 1.25e-12 ***     navTRUE                    0.0567863  0.0061501    9.233  < 2e-16 ***     fueldiesel                 0.0910481  0.0067196   13.550  < 2e-16 ***     fuelelectric               0.0309334  0.0471583    0.656  0.51190         age:transmissionautomatic -0.0113273  0.0015439   -7.337 2.72e-13 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Residual standard error: 0.1213 on 3429 degrees of freedom     Multiple R-squared:  0.9481,    Adjusted R-squared:  0.9478      F-statistic:  3684 on 17 and 3429 DF,  p-value: < 2.2e-16      Thank you very much!


title: Financial time series model

text: I have an interesting question that I think has not been asked yet here. I am building an AI that has as goal to predict how wrong a standard based-on- history model is. This is done based on Natural Language Processing(NLP), so from an external source (yes I have thousands of articles timed by the minute to do this!!). I come to point where I need to select my original model. The most simple way is linear regression from bars 0-40 (where 40 is 'now') and compare it to bars 5-45. Compare the difference of value 40 and 45 and voila we have the 'difference' in our model we can attempt predict/classify using news articles. Ofcourse, using linear regression is not suitable for this kind of task but what is? And even better are there libraries in .NET(C#) that can do that, since my math-to-code skills are not that great. I have been looking at polynomial regression but find that the difference can be out of proportion big. Also been looking at GARCH but I cannot seem to find a library(so not in R) that can do this without paying a part of my body. What also seems interesting is AR(i)MA but there I stopped and decided to get some advice on my problem. What would be suited for such a 'supervisor'? Would I need very complex regression or would linear even be good enough? Or in this case just use a neural network(already made a pretty good one for this once)? Greets


title: What is the probability of rejection in rejection sampling?

text: **Context:** I'm reading up on sampling in MCMC for Machine Learning. On page 5, it mentions rejection sampling:               repeat       sample xi ~ q(x) and u ~ U(0,1)       if uMq(xi) < p(xi):         accept xi     until enough samples      (where $M$ is chosen such that $p(x) \le Mq(x)$ for all $x$) **Question:** In the analysis, the paper says that it may work bad for high- dimensional settings. While I can see that, the paper gives a reason that I don't understand: $P(x \text{ accepted}) = P \left(u < \frac{p(x)}{Mq(x)} \right) = \frac{1}{M}$. This doesn't make sense to me. If the probability were constant, why even bother to evaluate $p(x)$? Should this be a "at most $\frac{1}{M}$? Or am I just misinterpreting the statement?


title: Variables which sum to a constant - ANOVA, MANOVA or none of these?

text: I asked a group of subjects to make a series of 12 binary choices regarding preferences. Let's say for arguments sake, these were between ugly (`ug`), attractive (`att`), and neutral (`neut`) faces. Hence, we have 4 `ug` vs `att`, 4 `ug` vs `neut` and 4 `att` vs `neut` choices. For each subject I summed the number of times each face was chosen. Hence, I have a 3 column table comprising a score (max 8) for `Att`, `Ug` and `Neut` for each subject. Each row sums to 12 hence the variables are negatively correlated. My questions:   * Are attractive faces preferred to ugly and if so:   * Is this driven by an attraction to `att` or an aversion to `ug` or both? - this is why we have choices with the neutral faces. I originally thought to do a repeated measures ANOVA followed by _post hoc_ tests to look for differences in ratings but i'm wondering if the fact that the DVs all sum to a constant is problematic because in essence the third variable - say $neut = 12-(ug+att)$. If so, is MANOVA the way to go, or how about chi-square?


title: Estimated distribution of eigenvalues for i.i.d. (uniform or normal) data

text: Assuming I have a data set with $d$ dimensions (e.g. $d=20$) so that each dimension is i.i.d. $X_i \sim U[0;1]$ (alternatively, each dimension $X_i \sim \mathcal N[0;1]$) and independent of each other. Now I draw a random object from this dataset and take the $k=3\cdot d$ nearest neighbors and compute PCA on this set. In contrast to what one might expect, the eigenvalues aren't all the same. In 20 dimensions uniform, a typical result looks like this:               0.11952316626613427, 0.1151758808663646, 0.11170020254046743, 0.1019390988585198,     0.0924502502204256, 0.08716272453538032, 0.0782945015348525, 0.06965903935713605,      0.06346159593226684, 0.054527131148532824, 0.05346303562884964, 0.04348400728546128,      0.042304834600062985, 0.03229641081461124, 0.031532033468325706, 0.0266801529298156,      0.020332085835946957, 0.01825531821510237, 0.01483790669963606, 0.0068195084468626625      For normal distributed data, the results appear to be very similar, at least when rescaling them to a total sum of $1$ (the $\mathcal N[0;1]^d$ distribution clearly has a higher variance in the first place). I wonder if there is any result that predicts this behavior? I'm looking for a test if the series of eigenvalues is somewhat regular, and how many of the eigenvalues are as expected and which ones significantly differ from the expected values. For a given (small) sample size $k$, is there a result if a correlation coefficient for two variables is significant? Even i.i.d. variables will have a non-0 result occasionally for low $k$.


title: Estimating population size from relative appearances

text: I know that widget X has a population of N_x. However, widget Y has an unknown population that I'd like to estimate. Both widgets appear with a differential frequency over time as I sample them. Time ... X counts .... Y counts   1 .......... 120 ............... 2   2 .......... 212 ............... 3   3 .......... 321 ............... 5   4 .......... 149 ............... 0   5 .......... 321 ............... 1 What is the total population of Y? With what confidence? I'm particularly interested in any special considerations that I many need to make when Widget Y is much rarer than the reference Widget X. Best, Paul


title: Critical region for difference in proportions

text: A new virus breaks out on a cruise ship. I want to test the hypothesis that males and females are equally likely to contract the virus. I am going to test 100 men and 100 women. Presumably if I find 87 women infected and 89 men I cannot safely reject the null. If on the other hand I find 11 men are affected and 20 women this would seem a reasonable basis on which to reject the null. So before running the test I want to define a rejection region. The rejection region is to be defined so that the probability of rejecting the null given that it is true is _at most_ 5%. Since this is a composite null – the probability of getting the disease can range from 0 to 1 – this condition must hold for each and every member of the null. How do I construct such a region?


title: Does it make sense to do CV-error-weighted model averaging?

text: We often average models together to create an aggregate prediction model. Some recent research suggests that simple model averages perform as well or better than model averages weighted by functions of information criterion scores. Weighted or simple, model averaging often (but not always) performs better than choosing the model with the best score. Many information criterion scores can be seen as approximations to different cross-validation schemes. Even though one model performs worse in cross-validation than another, it is still informative. I argue that an average model with weights as functions of cross- validation scores makes sense. Yet most people who do cross-validation to tune hyperparameters just choose the value of the hyperparameter that minimizes cross-validation error. In what circumstances might we expect a cross- validation-error-weighted average model to perform better than a single "best" model? I would love it if someone could point me to the most relevant research on this. Thank you.


title: Backpropagation algorithm

text: I got a slight confusion on the backpropagation algorithm used in multilayer perceptron (MLP). The error is adjusted by the cost function. In backpropagation, we are trying to adjust the weight of the hidden layers. The output error I can understand, that is, `e = d - y` [Without the subscripts]. The questions are:   1. How does one get the error of hidden layer? How does one calculate it?   2. If I backpropagate it, should I use it as a cost function of an adaptive filter or should I use a pointer (in C/C++) programming sense, to update the weight?


title: Stepwise logistic regression and sampling

text: I am fitting a stepwise logistic regression on a set of data in SPSS. In the procedure, I am fitting my model to a random subset that is approx. 60% of the total sample, which is about 330 cases. What I find interesting is that every time I re-sample my data, I am getting different variables popping in and out in the final model. A few predictors are always present in the final model, but others pop in and out depending on the sample. My question is this. What is the best way to handle this? I was hoping to see the convergence of predictor variables, but that isn't the case. Some models make much more intuitive sense from an operational view (and would be easier to explain to the decision makers), and others fit the data slightly better. In short, since variables are shuffling around, how would you recommend dealing with my situation? Many thanks in advance.


title: What is the interpretation/meaning of confidence intervals in misspecified models?

text: Consider the following model $Y_i = f(X_i) + e_i$ from which we observe n iid data points $\left( X_i, Y_i \right)_{i=1}^n$. Suppose that $X_i \in \mathbb{R}^d$ is a $d$ dimensional feature vector. And suppose that a ordinary least squares estimate is fit to data, that is, $\hat \beta = {\rm arg} \min_{\beta \in \mathbb{R}^d} \sum_i (Y_i - \sum_j X_{ij} \beta_j)^2$ Since a wrong model is estimated, what is the interpretation for the confidence interval around estimated coefficients? More generally, does it make sense to estimate confidence intervals around parameters in a misspecified model? And what does the confidence interval tell us in such a case?


title: Clustering (k-means, or otherwise) with a minimum cluster size constraint

text: I need to cluster units into $k$ clusters to minimize within-group sum of squares (WSS), but I need to ensure that the clusters each contain at least $m$ units. Any idea if any of R's clustering functions allow for clustering into $k$ clusters subject to a minimum cluster size constraint? kmeans() does not seem to offer a size constraint option.


title: How to compute the confidence intervals on regression coefficients in PLS?

text: The underlying model of PLS is that a given $n \times m$ matrix $X$ and $n$ vector $y$ are related by $$X = T P' + E,$$ $$y = T q' + f,$$ where $T$ is a latent $n \times k$ matrix, and $E, f$ are noise terms (sssuming $X, y$ are centered). PLS produces estimates of $T, P, q$, and a 'shortcut' vector of regression coefficients, $\hat{\beta}$ such that $y \sim X \hat{\beta}$. I would like to find the distribution of $\hat{\beta}$ under some simplifying assumptions, which should probably include the following:   1. The model is correct, _i.e._ $X = T P' + E,y = T q' + f$ for unknown $T, P, q$;   2. The number of latent factors, $k$, is known, and used in the PLS algorithm;   3. The actual error terms are i.i.d. zero-mean normal with known variances; This question is somewhat underdefined because there are scores of variants of 'the' PLS algorithm, but I would accept results for any of them. I would also accept guidance on how to estimate the distribution of $\hat{\beta}$ via _e.g._ a bootstrap, but perhaps that is a separate question.


title: Why would it not make sense to apply factor analysis after PCA?

text: I think that after apply PCA the covariance between components is reduced since the aim of PCA is to maximize variance. Covariance is needed by Factor analysis and thus it does not make sense to operate on the data in PCA space. Does this make any sense? Does someone have any other explanation?


title: probability that a variable is ONE OF the top k out of n when ordered

text: Suppose ($h_1,h_2,...,h_n$) is an $n\times 1$ vector. Let $h_i=g_iX_i$, where $g_i$ is a non-random variable which can vary across $i$ and $X_i$ is a random variable with Pareto Type I distribution. $$f_X(x) = \begin{cases} {\alpha x_m^\alpha \over x^{\alpha +1}}, & \text{if $x\ge x_m$} \\[2ex] 0, & \text{if $x\lt x_m$} \\ \end{cases}$$ What is the probability, as a function of $g$'s, that any given $h_i$ will be among the top $k$ largest of the $n$ $h$'s?


title: Difference of means when population mean is known?

text: I’ve approached a unique problem, and I want to make sure that I am going about the test that I am running in the correct way. Here’s a brief background: The goal of this project is to automate an advertisement system, and the goal of my question is to determine the p-value that one type of advertisement has a lower cost-per-user than another type of ad. **Data / Definitions:**   1. _Arrivals_ - Number of people that clicked on ad   2. _Signups -_ Number of people that signed up   3. _Spend_ - Total amount spent on each ad   4. *Cost_Per_User = Spend / Signups I have daily data on all of these variables for two different ads. This allows me to get the aggregate population cost_per_user, as well as the variance of the aggregate cost_per_user. Now, I want to test the probability of two different ads having different performance today, so that I might be able to adjust the budgeting in an informed manner. What would be the best way to use this data to run a proper test? Thanks in advance! PS: I think that this test Hypothesis Test for Difference of Means is probably the correct one, but I haven't been able to find the mechanics of it.


title: Sufficient statistics

text: I have independent observations $Y_1, Y_2, ...., Y_n$ that take values 1 or 0 with $$P(Y_i = 1) = \frac{1}{1 + \exp\{-(\alpha + \beta x_i)\}}$$ where $\alpha$ and $\beta$ are unknown constant and $x_1, x_2, ....,x_n$ are known, I have to show that $(\sum Y_i,\sum x_i Y_i)$is a minimal sufficient statistics for $(\alpha,\beta)$. We know that any statistics $T(Y)$ is minimal sufficient if $\frac{\prod_{i=1}^nf(y_i;\alpha, \beta)}{\prod_{i=1}^nf(z_i;\alpha, \beta)}$ is independent of $\alpha$ and $\beta$ if and only if $T(Y) = T(Z)$. But here the probability density of $Y$ is independent of $y$. Is there any other method to show the minimal statistics?


title: MATLAB - sig figs in MATLAB's normal approximation of p-value?

text: I am using MATLAB's ranksum function with the computation `method` of the p-value set to `approximate`. Is there a way to determine with what precision the approximation is being made?


title: Learning statistical concepts through data analysis exercises

text: I find that simple data analysis exercises can often help to illustrate and clarify statistical concepts. What data analysis exercises do you use to teach statistical concepts?


title: If you run OLS regression on cross sectional data, should you test for autocorrelation in residuals?

text: I have a set of observations, independent of time. I am wondering whether I should run any autocorrelation tests? It seems to me that it makes no sense, since there's no time component in my data. However, I actually tried serial correlation LM test, and it indicates strong autocorrelation of residuals. Does it make any sense? What I'm thinking is that I can actually rearrange observations in my dataset in any possible order, and this would change the autocorrelation in residuals. So the question is - should I care at all about autocorrelation in this case? And should I use Newey-West to adjust SE for it in case test indicates so? Thanks!


title: Training multiple models for classification using the same dataset

text: For my classification problem, I am trying to classify an object as Good or Bad. I have been able to create a good first classification step that separates the data into 2 groups using SVM. After tuning the parameters for the SVM using a training/holdout set (75% training, 25% holdout), I obtained the following results from the holdout set: Group 1 (model classified as Bad) consisted of 99% Bad objects, and Group 2 (model classified as Good) consisted of about 45% Good objects and 55% Bad objects. I verified the performance of the model using k-fold CV (k=5) and found the model to be stable and perform relatively consistently in terms of misclassification rates. Now, I want to pass these objects through another round of classification by training another model (may or may not be SVM) on my group 2 of maybe good/maybe bad objects to try and correctly classify this second group now that I have gotten rid of the obviously bad objects. I had a couple of thoughts, but am unsure of how to proceed. (1) My first idea was to use the data from the classified objects from the Holdout set to train another model. I was able to train another classification model from the results of the holdout set. The problem is I am using less than 25% of the original data, and I am worried of overfitting on a very small subset of my data. (2) My second idea was to gather the results of the 5-fold CV to create another dataset. My reasoning is that since the data is partitioned into 5 parts, and each part is classified into two groups from a model trained by the other 4 parts, I thought that I could aggregate the predicted results of the 5 parts to obtain a classified version of my original dataset and continue from there. The only problem is, I have a sinking feeling that both methods are no good. Could CV shed some light on some possible next steps? Thank you! **EDIT** Sorry, my question was badly worded. Let me try to clarify what I am trying to do. It can be thought of like a tree...   * Let me call the original dataset Node 0.    * I used classification method 1 to split Node 0 into Node 1 and Node 2.      * Node 1 has low misclassification rate (Mostly consists of bad objects)     * Node 2 has high misclassification rate (Roughly even mix of good and bad objects)   * I now want to use classification method 2 to split Node 2 into Node 3 and 4 The "classification method" can be anything (LDA, QDA, SVM, CART, Random Forest, etc). So I guess what I am trying to achieve here is a "classification" tree ( _not_ CART), where each node is subjected to a _different_ classification method to obtain an overall high "class purity". Basically, I want to use a mix of different classification methods to obtain reasonable results. My problem lies in the loss of training data after the first split. I run out of usable data after I run it through "classification method 1", which was SVM in my case.