title: What is the minimum viable cell size for 2x2 ANOVA?

text: I have a 2x2, between-subjects experimental design (2 independent variables (IVs) with 2 levels each) and one dependent variable (DV). My data are unbalanced and an interaction between the IVs seems likely, so I plan to use ANOVA with Type III Sums of Squares to test whether my DV is different across one of my IVs while controlling for any influence of the other IV (or the interaction, if present). My question is this: what is the minimum number of data points that I need in each of my 4 cells, to ensure that ANOVA likely won't give spurious results? I'm aware of Simmons et al.'s recommendation (2011) of having at least 20 data point in each cell, but that recommendation only serves to control the test's rate of false negatives. What I'm more worried about is that with a small enough cell size, the statistics on which ANOVA is based probably won't be very reliable, and so the results of the test would be unreliable as well (in terms of either false negatives _or_ false positives). Are there any papers or texts out there that have studied and made recommendations regarding this concern?


title: Minimum sample size for t-test

text: Can you do a t-test or ANOVA between 2 or 3 groups respectively if one of those groups has a sample size of 2? Are the results interpretable in the same way as with larger sample sizes(n=4,5,6?). If you can provide any references to read more about very small sample size statistics, that would be helpful.


title: Neural Network Process Question - Updating weights after each training set

text: When creating a neural network, do I update the weights after each run of forward then back propogation? Or do I just keep the random weights and update the Delta variables? I am looking at slide 8 on these notes:https://d396qusza40orc.cloudfront.net/ml/docs/slides/Lecture9.pdf It says:               For i = m        Set a(1) = x(i)        Perform Forward-Propogation        Compute delta      Compute DELTA      QUESTION: Do I update the Weights that I use in Forward-propogation, or do I       use random weights and just keep updating the accumulator 'DELTA'? And if I       update the weights, do I set them to DELTA?


title: Random forest like procedure for regression or other statistical models

text: I'm wondering if there exist methods similar to one used in random forest algorithm - I mean taking simultaneously bootstrap sample and random subset of features, then building statistisal model. Have anyone took this approach for building set of regression models ? Is this approach ( random subsample plus random subset of features ) somehow universal ? Edition : the question is all about the possibility of putting some other model in the place of classification tree in random forest, what is left is some kind of meta-algorithm (as bagging can be view as meta-algorithm) = bagging + random subsets of variables


title: How does integrating the Kolmogorov forward equation give $P = \exp (Qt)$?

text: If $Q$ is a generator matrix of a continuous time Markov chain (CTMC), and I need to use this matrix to solve the Kolmogorov forward equation, I would need to start by integrating it. But I haven't got a clue how to do it. Can someone show me please? I know something like, let's assume $i$ represents the current state of a CTMC. Then, the forward equation basically tells us that we can work out $X(i + 1)$ by doing $$X(i + 1) = X(i) \cdot (Id_2 + Q)$$ To look at the difference in time, we can subtract $X(i)$ from both sides and get $$X(i + 1) - X(i) = X(i) Q$$ Thinking of this in terms of functions instead of matrices, we can say that this can be written as $$ P'(t) = X(i) Q dt$$ But I don't get how you get to this bit and how you can integrate from here. I would really appreciate any help. Thank you.


title: Why do the results in cross validation changes whenever I shuffle my training data?

text: The result that I'm talking about is the mean of cross validation from the function of sklearn. I tried to shuffle my training data then applied the CV function (shuffle then CV). I did that for several times and each result was different, some are higher than the others. Is the arrangement of data from the training data relevant? If so, then how? **Additions:** Sorry I didn't define my problem well. I'm doing a supervised binary classification. My train variable consists of the feature and classification. Here is part of my code:               import numpy as np     from sklearn import cross_validation     from sklearn.linear_model import LogisticRegression          np.random.shuffle(train)     clf = LogisticRegression()     clf = clf.fit(train[0::,1::],train[0::,0])          mean = cross_validation.cross_val_score(clf, train[0::,1::],train[0::,0], cv=cross_validation.StratifiedKFold(train[:,0],5)).mean()     print mean     #sample for the train:     #train = [[1,0.5,0.3,0.6],[0,0.3,0.2,0.1],[0,0.1,0.9,0.7]]      Here I have n_fold=5. If I remove the np.random.shuffle(train) my result for the mean is approximately 66% and it stays the same even after running the program a couple of times. However, if I include the shuffle part, my mean changes (sometimes it increases and sometimes it decreases). And my question is, why does shuffling my training data changes my mean? If the order of data in my 'train' is irrelevant, then why does my mean changes after shuffling? This is just part of my code so assume that train variable is defined somewhere. My train variable is similar to the sample, but it consists of a lot more feature.


title: Spherical platykurtic random cloud

text: I'd like to generate a multivariate continuous data which is **globular** cloud, like multivariate standard normal data is, but which is more **platykurtic** than normal data. There are many ways to get platykurtic data (e.g. from beta distribution) but the more platykurtic they are the more they are rectangular shape in multivariate space. But I want globular, hyperspherical random data. Unimodal data. Preferably from not bounded distribution, if possible. And with an option to vary the amount of flat kurtosis. **Can you suggest** a distribution or a trick to generate? P.S. Saying "hyperspherical" I mean "any dimensionality" (not "high dimensionality"). That is, I imply 2D case to be just particular case; I'm interested in it as well.


title: Why is the plot of residuals against fitted values a horinzontal line when the dependent variable is linearly related to the indenpendent variable?

text: In ordinary least squares regression (OLS), if the plot of the residuals against the fitted values form a horizontal line around 0, then we can say that the dependent variable is linearly related to the independent variable. I had thought that this is true because $E(y_i - \hat{y}_I)=0$ when the dependent variable is linearly related to the independent variable, see here. However, suppose: $y_i = \alpha + \sin(x_i) + \epsilon_i$. Then $E(y_i - \hat{y}_i)$ is still 0, see here but then the plot of its residuals against its fitted value is no longer a horizontal line around 0, as this R code shows:               n <- 10^3     df <- data.frame(x=runif(n, 1, 10))     df$mean.y.given.x <- sin(df$x)     df$y <- df$mean.y.given.x + rnorm(n)     model <- lm(y ~ x, data=df)     plot(predict(model, newdata=df), residuals(model))     abline(a=0,b=0,col='blue')      ![enter image description here](http://i.stack.imgur.com/kyBwt.png) So my question is, which assumption(s) of OLS that causes the plot of the residuals and the fitted value to be a horizontal line around 0 and why/how is it true?


title: How to determine whether a sales campaign is effective?

text: Assuming that I have an online ecommerce store. Recently, for one of the products, I want to run a sales campaign for the product for a period of time. The effectiveness is defined solely in terms of increasing sales. But the problem is I have no idea how to measure this. From my experience, sales volume fluctuates from month to month, and there could be other factors ( undetected) that dampen or excites my sales record. There are just too many variables, unknown and known involved and very difficult to disentangle them. Doing controlled experiment like A/B testing the webpages( showing two versions of webpages to different visitors, one with promotion, one without promotion) is simply not practical, because the visitors may not buy instantly. They might buy after a while, so they are in a rude surprise to find that the price has changed so suddenly. Even if it is practical, I don't have the resources to invest in A/B testing, my funds are very limited. How to determine, and what are the statistical valid experiments that I can do, to check whether my sales campaign is effective or not?


title: What are some papers analysing the effect of missing data in machine learning?

text: In machine learning it seems plausible that missing data effects predictive accuracies and learning times. Does anyone know of any literature analysing this? I have found plenty of literature discussing techniques to deal with missing data, but no literature analysing the effect that it has, e.g. with 25% of the original features missing, then predictive accuracies decrease by 50%, etc.


title: Squared sum of errors

text: Why does $\mbox{SST} = \mbox{SSR} + \mbox{SSE}$ hold? It is understandable that $\sqrt{\mbox{SST}} = \sqrt{\mbox{SSR}} + \sqrt{\mbox{SSE}}$ but how can $\mbox{SST} = \mbox{SSR} + \mbox{SSE}$ hold? These are squared values.


title: CDF/ ECDF plot for data with two attributes

text: I have a data in the following format: $$ \begin{array}{rr} \textbf{colm_1} & \textbf{colm_2}\\\ 3 & 1\\\ 10 & 0\\\ 3 & 0\\\ 100 & 1\\\ . & .\\\ . & . \end{array} $$ _colm_1_ are basically whole numbers and _colm_2_ denotes binary values say success or failure. What I am expecting is that, for lower values of attributes in _column_1_ the probability of success will be less, while as the values of _column_1_ increases, the probability of success escalates. Now, how can I plot the CDF (or ECDF) of the above data?. From my perspective it should basically say _P(success | colm_1 = **Value** ) <= Some Number between 0-1_. Where **Value** is any whole number in _colm_1_. I hope I am clear in my explanation, if not please rectify me. Finally, I mainly use python's numpy and scipy. How can I transform my data so that I can plot the ecdf?


title: What is the difference between fixed effects model and random effects model for a meta-analysis of sample correlations?

text: What is the difference between fixed effects model and random effects model for a meta-analysis of sample correlations ?


title: Can I difference after fitting a time series regression model?

text: Suppose that I have a time series that exhibits a notable trend, and I want to test a hypothesis that a second variable is related to that trend. I fit a linear regression model with that second variable as a predictor, and it accounts for a significant amount of variance. However, there is some part of the trend that it does not account for. I then want to include AR or MA terms to the model to account for the autocorrelated residuals, but ARIMA models assume a stationary series, and visual inspection of the series and an ADF test indicate that the residual error series is not stationary in the mean (i.e., there is some trend left unaccounted for by the regression model). Is it fine to difference _after_ the regression model has been fit to remove this source of non-stationarity, or should one always try and find a regression model that explains _all_ the trend in the data? Thanks in advance!


title: Weighing a variable based on distance to a known date

text: In a bankruptcy model, you want to assign a higher weight to a variable as a big event date approaches (such as a company's quarterly earnings announcement date) and you reverse this weighing as you go away from the last such date. As expected, a higher positive weight would contribute more towards predicting a bankruptcy. What are the techniques that can used for deciding such a weighing scheme? I'm aware of few but wanted to see what experts might have to say.


title: Comparing two variance matrices

text: I am looking for bibliographical reference for comparing two variance matrices with he following criterion: $\text{Var}[X] \geq \text{Var}[Y] \quad \text{if} \quad \text{Var}[X]-\text{Var}[Y] \succeq 0$ This means that $X$ has grater variance than $Y$ if the difference of their variance matrices is semi definite positive. Any citable reference to books or works that explain why this is a valid and common criterion to compare variability of random vectors is appreciated. If that is not possible any help on how to search for it would be ok. Thanks


title: Statistics:Power Question Update

text: I have already asked this question but I did not initially give out all the details on the question, like the question being a multiple choice question. You manufacture and sell a iron rod whose electrical conductivity is supposed to be 10.1. You plan to make 6 measurements of the conductivity of a rod you plan to sell. You know that the standard deviation of your measurements is σ=0.1. If the product meets specifications, the mean of many measurements will be 10.1. You will therefore test H0:μ=10.1 Ha:μ≠10.1 The question is: Explain why the test you plan will not adequately protect you against selling a rod with conductivity 10.15. The possible answers are: A. Because the probability of mistakenly concluding that μ = 10.15 when the product is unsuitable is too high. B. Because the probability of mistakenly concluding that μ = 10.15 when the product is suitable is too high. C. Because the probability of mistakenly concluding that μ = 10.15 when actually μ ≠ 10.15 is too high. D. Because the power of a test against the alternative μ = 10.15 is only a probability, not a certainty. I believe the answer is Because the probability of mistakenly concluding that μ=10.15 when the product is suitable is too high because the power is 0.24 and that means there is a much greater chance that the rod is not going to eliminate the null hypothesis (0.76). Also what is the difference if I increase the sample size to more than six. Am I right or wrong? Or is it one of the other ones. Thanks everyone


title: How to deal with chromosomes of different lengths genes in genetic algorithm?

text: If I have 6 products ($x$) and 11 manufacturers ($y$) some products manufactured by some of these manufacturers; ($x,y$) may equal none or any other value. For each combination of $x$ and $y$ there is and expected values for cost, time and transportation petrol according to distance between the manufacturer and the work site. I have two problems to use binary genetic algorithm (GA) to find the Pareto sets of these combinations:   1. the gene represent the product contain only 3 digits ($6<7$) and for gene used to represent the manufacturer it is 4 digits ($8<11<16$) so the chromosome is of 7 digits. the problem is that in all examples those I have reviewed the length of genes are equal. so will this be a problem for solving this issue.   2. assume product $x_1$ and $x_5$ are not manufactured by $y_3$ and $y_7$ that means the fitness function of $(x_1,y_3)$ and $(x_5,y_7) = 0$ and since the problem is minimization the solution may be trapped in these two points. the question how such cases can be solved? Thanks for your expected valuable help.


title: is Julian day a count data or a continous data

text: A stupid question but has confused me. My dataframe looks like this:               Julian day     Year      153           1951      161           1952      167           1953       .             .       .             .      161            2007      I have Julian day for each year when an event occurs and I am trying to determine whether the Julian day of that event happening is changing with time or not. Do I treat Julian day as continuous or count data. I guess if I treat it as continuous data, I will do a linear regression of Julian day against time and if it's a count data, I will do a Poisson regression against time.


title: Hierarchical models for loss reserving

text: I am trying to fit a Weibull growth curve to my data, but I get a singularity in `backsolve` at level 0 block 1 error:               DATA     head(trian)         AY DEV     CUM     1 1991   6  357848     2 1991  18 1124788     3 1991  30 1735330     4 1991  42 2182708     5 1991  54 2745596     6 1991  66 3319994      and the formula is:                w1 = nlme(CUM~ult*(1-exp(-(DEV/theta)^omega)),                fixed=list(ult~1,omega~1,theta~1),                random=ult~1|AY,                weights=varPower(fixed=0.5),                data=trian,                start=c(ult=5000,omega=1.4,theta=45))


title: Friedman test block contents

text: For comparing computational intelligence algorithms over multiple functions I have seen the Friedman test used like the following [1]:                          Algorithm1 Algorithm2 ... Algorithmk     Function1  mean(A1F1) mean(A2F1) ... mean(AkF1)     Function2  mean(A1F2) mean(A2F2) ... mean(AkF2)     ...     Functionn  mean(A1Fn) mean(A2Fn) ... mean(AkFn)      where mean(x) is the mean result over a number of runs (say M). However the wikipedia page says that there should be "a single observation at the intersection of each block and treatment". So I was thinking something like:                               Algorithm1 Algorithm2 ... AlgorithmK     Function1_run1  A1F1_run1  A2F1_run1  ... AKF1_run1     Function1_run2  A1F1_run2  A2F1_run2  ... AKF1_run2     ...     Function1_runM  A1F1_runM  A2F1_runM  ... AkF1_runM          Function2_run1  A1F2_run1  A2F2_run1  ... AKF2_run1     Function2_run2  A1F2_run2  A2F2_run2  ... AKF2_run2     ...     Function2_runM  A1F2_runM  A2F2_runM  ... AkF2_runM          FunctionN_run1  A1FN_run1  A2FN_run1  ... AKFN_run1     FunctionN_run2  A1FN_run2  A2FN_run2  ... AKFN_run2     ...     FunctionN_runM  A1FN_runM  A2FN_runM  ... AkFN_runM      Each run is completely independent from every other run. My questions are: Is the second method valid for the Friedman test? If so, can the order of the runs be changed for a given function since they are independent? [1] S. García, D. Molina, M. Lozano, F. Herrera, A Study on the Use of Non- Parametric Tests for Analyzing the Evolutionary Algorithms' Behaviour: A Case Study on the CEC'2005 Special Session on Real Parameter Optimization. Journal of Heuristics, 15 (2009) 617-644. [2] http://en.wikipedia.org/wiki/Friedman_test


title: Decision trees variable (feature) scaling and variable (feature) normalization (tuning) required in which implementations?

text: In many machine learning algorithms, feature scaling (aka variable scaling, normalization) is a common prepocessing step Wikipedia - Feature Scaling \-- this question was close Question#41704 - How and why do normalization and feature scaling work? I have two questions specifically in regards to Decision Trees:   1. Are there any decision tree implementations that would require feature scaling? I am under the impression that most algorithms' splitting criteria are indifferent to scale.   2. Consider these variables: (1) Units, (2) Hours, (3) Units per Hour -- is it best to leave these three variables "as-is" when fed into a decision tree or do we run into some type of conflict since the "normalized" variable (3) is relatable to (1) and (2)? That is, would you attack this situation by throwing all three variables into the mix, or would you typically choose some combination of the three or simply use the "normalized/standardized" feature (3)?


title: Interpretation of descriptive statistics for dummy variable

text: How can I describe descriptive statistics for a dummy variable (gender of worker in a shop)? Let's say this is the info that I have:               mean :         0.47     median :       0     max :          1     min :          0     std. dev :     0.4998     skewness :     0.101     kurtosis :     1.01     jarque bera : 85.67     probability :  0      I know that some of the information is useless since it's a dummy variable. So how do I interpret it in words?


title: How can I plot this simple graph (Refer Image) in R?

text: I'm very new to R Programming. So please excuse for such a simple doubt. ![How to plots this graph in R?](http://i.stack.imgur.com/pcmle.png) I want to plot the above graph. The x & y values are sequence from 0 to 2560. I want plot a a curve on the points where x*y=10^6. What are the line required in R Programming Language.


title: Optimal lookback period and standard deviation level for determining in which state of Hidden Markov model with Gaussians

text: I am trying to solve the following task analytically, yet can only come up with heuristical solution for special cases based on monte carlo simulations: Given a Hidden Markov model with two states and two Gaussians where all parameters are known (i.e. transition matrix and means and standard deviations of the Gaussians): Additionally given an arbitrary instance of data with many periods that is produced by this model my task is to determine in what state I am at the moment and the only information I have is the standard deviation of a rolling window of $n$ periods: How do I determine analytically the optimal (minimal) lookback period $n$ and best separating level (of the standard deviation) $s$ to determine in which state I am at the moment? Put another way what I am asking is given the model with all its parameters how to derive at $n$ and $s$ analytically to determine the state in real time, i.e. now.


title: Probability of Selection vs. Probability of Inclusion

text: I'm studying an article about sampling (Adaptive Cluster Sampling from Steven K. Thompson) and I'm confused about two concepts. In some point, he says that the **probability of selection** of unity $i$ on any of the $n_1$ (initial sample) draws is $p_i = (m_i + a_i)/N$. He also says that the **probability that unity $i$ in included** in the sample is $$ \alpha_i = 1 - \binom{N - m_i -a_i}{n_1}/\binom{N}{n_1} $$ What is that difference between **probability of selection** and **probability of inclusion**? Initially I though that probability of selection was the probability that unity $i$ is selected for the initial sample and probability of inclusion is the probability of inclusion is the probability that unity $i$ is included in the sample after the sequential selection. But since simple random sampling is used for initial selection, the formula shouldn't take into account $m_i$ nor $a_i$ which are the size of network and edge that unity $i$ is part of. It should only take into account the total number of unities and the size of initial sample. So, it's not clear to me what is each of those probabilities.


title: How can I fit a zero inflated poisson model with only offset (without coefficients)?

text: I have already got a poisson estimated lambda, and actual result `y`, and I would like to see if the model is good. To start with, I check if the dispersion is alright.               glm(y ~ 0, offset=log(lambda), data=data, family=quasipoisson)      And then I would like to see if there is zero inflation factor by using `pscl` package.               zeroinfl(y ~ 1, offset=log(lambda), data=data ,dist="poisson")      However, I am not able to do this:               zeroinfl(y ~ 0, offset=log(lambda), data=data ,dist="poisson")      By introducing an interception, it is likely to interfere with the zero estimator so it becomes unclear if there is any zero (and how big) inflation in the original model. Is it possible to have 0 estimated parameters (besides zeros)?


title: Multiple Linear Regression Models Comparison Based on R-squared And Residual Errors

text: am currently working on a problem where I have to calibrate weather parameters at a ground location using Satellite data available (1979-2012) over rectangular grid points and surface observatory data for a location (1980-2000). I have built two models - One for January months over 1980-2000 period and other taking data from November to April over 1980-2000. Now the problem arises when I compare adjusted R-squared for both the models. Nov-Apr model gives me higher R-squared (0.7 vs. 0.2) whereas residual errors tell me the opposite story, giving poor results in the case of Nov-Apr model. One option would be to calculate different components of R-squared and see if it can explain the differences. So how do I calculate and explain these components in R? or is there any better and efficient way? Thanks in advance.


title: Basic binomial question

text: The following appeared on an assignment of mine (already turned in). I contend that not enough information is given to provide an answer.... it seems pretty cut and clear to me. However, instructor insisted it's solvable in minitab. Can you help me figure out what I'm not understanding? How do you solve this without a model of distribution of weekly demand, or at least an average value to use as constant approximation. I must be missing something simple. The problem: Consider a service company. 10% of the weekly demand is for a service category named "X" [Assume service categories are mutually exclusive]. The company must revise their resource plan iff there are too few customer orders(less than one/week) or too many customer orders (more than five/week) of service category "X". For the next 12 weeks, what is the probability that the company will **not** need to revise the resource plan? Thanks


title: Two sided permutation test

text: I'm trying to write a two sided permutation test to test the alternative hypothesis there is a difference in the medians of 2 independent samples. My question is this: am I calculating the p-value correctly? Thanks SO much!               sample1 <- groundwater$West[!is.na(groundwater$West)]     sample2 <- groundwater$East     ts <- median(sample1) - median(sample2)     > ts     [1] 0.105     R <- 9999     all <- c(sample1, sample2)     k <- 1:length(all)     reps <- numeric(R)     for (i in 1:R) {         m <- sample(k, size=length(sample1), replace=FALSE)         permsample1 <- all[m]         permsample2 <- all[-m]         reps[i] <- median(permsample1) - median(permsample2)     }     pvalue <- sum(c(ts, reps) <= ts)/9999     > pvalue     [1] 0.9223922


title: Mean Squared Error and Residual Sum of Squares

text: Looking at the Wikipedia definitions of:   * Mean Squared Error (MSE)   * Residual Sum of Squares (RSS) It looks to me that $MSE = \frac{1}{N} RSS = \frac{1}{N} \sum (f_i -y_i)^2$ where $N$ is he number of samples and $f_i$ is our estimation of $y_i$. However, none of the Wikipedia articles mention this relationship. Why? Am I missing something?


title: Dealing with underdetermination in Bayesian models

text: Bayesian models are supposedly well equipped to deal with high-dimensionality problems, and can handle sparse data well, too. But suppose I've created a model that estimate more parameters than there are data points. Are there tricks to deal with this?


title: Taylor series expansion of maximum likelihood estimator, Newton-Raphson, Fisher scoring and distribution of MLE by Delta method

text: Assume $\ell\left(\theta\right)$ is the log-likelihood of parameter vector $\theta$ and $\widehat{\theta}$ is the maximum likelihood estimator of $\theta$ then the Taylor series of $\ell\left(\theta\right)$ about $\widehat{\theta}$ is \begin{align*} \ell\left(\theta\right) & \approxeq\ell\left(\widehat{\theta}\right)+\frac{\partial\ell\left(\theta\right)}{\partial\theta}\Bigr|_{\theta=\widehat{\theta}}\left(\theta-\widehat{\theta}\right)+\frac{1}{2}\left(\theta-\widehat{\theta}\right)^{\prime}\frac{\partial^{2}\ell\left(\theta\right)}{\partial\theta\partial\theta^{\prime}}\Bigr|_{\theta=\widehat{\theta}}\left(\theta-\widehat{\theta}\right)\\\ \frac{\ell\left(\theta\right)}{\partial\theta} & \approxeq\mathbf{0}+\left(\mathbf{1}-\mathbf{0}\right)\frac{\partial\ell\left(\theta\right)}{\partial\theta}\Bigr|_{\theta=\widehat{\theta}}+\frac{\partial^{2}\ell\left(\theta\right)}{\partial\theta\partial\theta^{\prime}}\Bigr|_{\theta=\widehat{\theta}}\left(\theta-\widehat{\theta}\right)\quad\overset{\textrm{set}}{=}\quad\mathbf{0}\\\ \\\ \theta-\widehat{\theta} & =-\left[\frac{\partial^{2}\ell\left(\theta\right)}{\partial\theta\partial\theta^{\prime}}\Bigr|_{\theta=\widehat{\theta}}\right]^{-}\left[\frac{\partial\ell\left(\theta\right)}{\partial\theta}\Bigr|_{\theta=\widehat{\theta}}\right]\\\ \widehat{\theta}-\theta & =\left[\frac{\partial^{2}\ell\left(\theta\right)}{\partial\theta\partial\theta^{\prime}}\Bigr|_{\theta=\widehat{\theta}}\right]^{-}\left[\frac{\partial\ell\left(\theta\right)}{\partial\theta}\Bigr|_{\theta=\widehat{\theta}}\right]\\\ \widehat{\theta}-\theta & =\left[\mathbb{H}\left(\theta\right)\Bigr|_{\theta=\widehat{\theta}}\right]^{-}\left[\mathbb{S}\left(\theta\right)\Bigr|_{\theta=\widehat{\theta}}\right] \end{align*} As $ \theta=\widehat{\theta}-\left[\mathbb{H}\left(\theta\right)\Bigr|_{\theta=\widehat{\theta}}\right]^{-}\left[\mathbb{S}\left(\theta\right)\Bigr|_{\theta=\widehat{\theta}}\right] $ So \begin{align*} \theta^{\left(m+1\right)} & =\theta^{\left(m\right)}-\left[\mathbb{H}\left(\theta^{\left(m\right)}\right)\right]^{-}\mathbb{S}\left(\theta^{\left(m\right)}\right)\quad\quad\left({\textrm{Newton- Raphson}}\right)\\\ \\\ \\\ \theta^{\left(m+1\right)} & =\theta^{\left(m\right)}-\left[\mathbb{E}\left\\{ \mathbb{H}\left(\theta^{\left(m\right)}\right)\right\\} \right]^{-}\mathbb{S}\left(\theta^{\left(m\right)}\right)\\\ \theta^{\left(m+1\right)} & =\theta^{\left(m\right)}+\left[\mathbb{I}\left(\theta^{\left(m\right)}\right)\right]^{-}\mathbb{S}\left(\theta^{\left(m\right)}\right)\quad\quad\left(\textrm{Fisher Scoring}\right) \end{align*} **Questions**   1. I'm not sure that my derivation is correct or not. I have seen at least two different versions of derivation.   2. What would be the mean, variance and distribution of $\widehat{\theta}$ (Might by Delta method)?


title: Can constrained optimization techniques be applied to unconstrained problems?

text: I am looking at using Interior Point method for optimizing a convex function. The convex function is basically the log-likelihood of a binary logistic regression model. Can I use this technique? In generally, is there anything that prevents applying a constrained optimization technique to an unconstrained problem? From what I think, an unconstrained problem is just a constrained problem without the constraints and thus should be solvable using these techniques.


title: How to add periodic component to linear regression model?

text: I have some cumulative frequency data. A line $y=ax+b$ looks like it fits the data extremely well, but there is cyclic/periodic wiggle in the line. I would like to estimate when the cumulative frequency will reach a certain value $c$. When I plot the residuals vs. fitted values, I get a beautiful sinusoidal behavior. Now, to add another complication, note that in the residuals plots ![alt text](http://i.stack.imgur.com/3GIg5.png) there are two cycles that have lower values than the others, which represents a weekend effect that also must be taken into account. So, where do I go from here? How can I combine some cosine, sine, or cyclic term into a regression model to approx. estimate when the cumulative frequency will equal $c$?


title: fisher's test vs pairwise fisher's test

text: I have a fold enrichment table for different days as given below. I need to find the significance of every fold enrichment.               fold_enrichment <- matrix(c(0.3,0.43,0.5,0, 0.23,1.3,0.5,1, 1,0,0,6, 0.3,0.5,0.2,2), 4, 4,                   dimnames = list(Day = c("Day1", "Day7", "Day14", "Day28"),                                   Group = c("A", "B", "C", "D")))      In R, I tried using `fisher.test(fold_enrichment)` and it gives me a single p-value > Warning in fisher.test(fold_enrichment) : 'x' has been rounded to integer: > Mean relative difference: 0.7652582 > > Fisher's Exact Test for Count Data > > data: fold_enrichment p-value = 0.6182 alternative hypothesis: two.sided I have also seen another option               pairwise.fisher.test(Day,Group, p.adjust.method="bonferroni")      Can anyone suggest me the difference between fisher's test vs pairwise fisher's test and what would be the right solution for my problem?


title: Regression for really small data with high degree of multicollinearity and outliers

text: I'm working on a promotional response analysis. I have a **really small** real world dataset with 25 observations and 15 variables. The variables have a **high degree of multicollinearity** and some have **outliers**. Also, I cannot use machine learning methods because I need interpretable coefficients. Here is what I have tried so far:   1. I used GLM, but all variables seemed insignificant. I know for a fact that the dependent variable is responsive to many independent variables.   2. I used robust regression, but basically the program fails because the software showed an error that stated, "The number of observations must be at least twice the number of coefficients." I would really appreciate if you could suggest some methods/techniques or strategies to solve this problem. Thanks and regards.


title: Why are my prediction intervals a little too narrow?

text: I'm trying to make a little function to make population-level predictions from a model that includes random effects, which in turn was fit using a package that doesn't support `interval="prediction"`. (As far as I can tell, the reason why the package doesn't do so is that prediction is fraught when trying to predict within one of the model's levels -- but I just want population- level predictions). So anyway, testing my code and comparing against the in-built code in `predict.lm`, my intervals are a little too narrow. Anybody able to help me figure out why?               set.seed(13)     n=100     x1 = runif(n)     x2 = runif(n)      X = cbind(x1,x2)          y = x1 + x2 + rnorm(n)     m = lm(y~x1+x2)     summary(m)          d = predict(m,interval = "prediction")     auto = d[,3]-d[,1]          getpred = function(vals){         qt(.975,n-(2+1))*sqrt(var(m$resid)* (1+ t(vals) %*% solve(t(X)%*%X) %*%vals))         }     man = apply(X,1,getpred)          plot(density(auto),xlim=c(2.1,2.3),ylim = c(0,30))     lines(density(man),col="red")      Thanks in advance!


title: Why is a symmetric distribution sufficient for the sample mean and variance to be uncorrelated?

text: While reading, I came across the puzzling statement that the sample mean and variance _are uncorrelated only in symmetric distributions_ and there is strong correlation if the distribution is heavily skewed. First of all, is it true? I already know this holds for the case of a _normal_ population, but can that proved for _every_ symmetric distribution? Although I think the result is counter-intuitive, I would appreciate simple arguments. Thank you


title: Approach for credit scoring for an agricultural products/chemicals company

text: I am currently working on a project for a large agri-business company. We currently have a credit policy that gives scores and classifies the debtors of the company into 5 segments - VLR (Very low risk), LR, MR, HR, VHR. The customers are mostly sole proprietors who do not have audited financial details. Hence most of them get classified into HR and VHR categories. e.g. 97% of the customers do not have financial data and 98% of the total customers are getting HR/VHR rating. We are working on identifying new variables which could help in further segmentation of HR/VHR customers. As a default indicator we will be using the incremental cost to company due to default/delay. My query is can I use a supervised learning algorithm like K Nearest Neighbor for classification? Because current classification data is that has 98% HR/VHR is not correct because of missing financial data. I have seen some examples of logistic regression, discriminant analytics, k means clustering etc. So how to decide which to use? Further before running any of these tests I suppose data needs to be cleaned. Are there any clear pointers, guidelines for data cleanup? I am a management student, so I have been introduced to statistical techniques but I have only basic preliminary knowledge of advanced techniques I mentioned above. So any guidance would be appreciated.


title: Is there an R function to do a survival analysis with right censoring + nested + crossed factors

text: I have this dataset to model, but I'm not sure how to do it. I want to model the **surviving probability** of different populations of two **species** depending on a **treatment** applied. Populations should be **nested within species** , as they belong to either species 1 or species 2. And, as I have replicated measures for each population x treatment combination, I need to include this as a random effect as well. What R function can I use to do that? Is there one? Thanks.


title: PLS Regression and collinearity

text: From what i know PLS regression is used when there is more variables than observations and when there exist multicollinearity between the independent variables. I have data for a regression model that does not show signs of multicollinearity between independent variables, but the relationship between independent and dependent variables is not very strong although there is a clear positive relationship with correlation values of 0.1 to 0.5 for different independent variables compared to the dependent variable. I would like to know whether the concept of PLS regression can extend to collinearity in general, meaning if the relationship between the independent variables and dependent variables is not strong then PLS regression is a suggested regression model or are there other models that are more suitable?


title: Test for differences between parameters of models estimated from partially overlaping samples

text: Suppose I have $n$ observations $(\pmb{x}_i,y_i)$ and two subsets of indexes of $\\{1:n\\}$, $S_1$ and $S_2$, with $S_1\neq S_2$, $\\#\\{S_1\\}\neq\\#\\{S_2\\}$, and $\\{S_1\cap S_2\\}\neq\emptyset$. Suppose I call $\hat{\beta}_1$ the vector of OLS-estimated coefficients of the regression $y_i\tilde{}\pmb{x}_i|i\in S_1$ and likewise for $\hat{\beta}_2$. My question is what is the right/best test for $$H_0:\beta_1-\beta_2=0$$ Some thoughts:   * My first inclination was for a Wald test, but how do we combine the two covariance matrices from the two regressions? Perhaps by pooling them?   * Another idea I thought about was to use the Bhattacharyya distance with $\beta_1$ against $\beta_2$ (together with their respective covariance matrices).   * A third possibility would be to merge the two samples and add an interaction effect _on every coefficient_ with a dummy taking value 1 for all observations that are, say, members of $S_2$ only, but because of the overlap between, $S_1$ and $S_2$ i'm not sure this is really measuring what I want.  Any thoughts? # EDIT1: Here is a simple `R` code to implement the solution proposed by Bill below:               library("MASS")     n<-200     p<-5          b<-rnorm(p)     x<-mvrnorm(n,rep(0,p-1),diag(p-1))     y<-cbind(1,x)%*%b+rnorm(n)     L<-sample(1:200,100,replace=TRUE)     L<-list(mod1=L[1:60],mod2=L[61:100])          M1<-lm(y[L[[1]]]~x[L[[1]],])     M2<-lm(y[L[[2]]]~x[L[[2]],])          fx01<-function(M1,M2){         X1<-cbind(1,M1$model[,-1])          X2<-cbind(1,M2$model[,-1])         nt<-nrow(X1)+nrow(X2)         D1<-diag(nt)[1:nrow(X1),]            D2<-diag(nt)[1:nrow(X2),]            st<-crossprod(c(M1$resid,M2$resid))/(nt-2*ncol(X1))         P1<-solve(crossprod(X1))%*%crossprod(X1,D1)-solve(crossprod(X2))%*%crossprod(X2,D2)         P1<-tcrossprod(P1)*drop(st)         mahalanobis(M1$coef,M2$coef,P1)     }


title: What measurements or tests should I take for testing stability of patterns?

text: I have a dataset containing behaviours of users. Let's say we have 1000 users and 40 similar tasks. The tasks are similar but the difficulty level vary.To do a task, people can follow 10 patterns. Since the 40 tasks are similar, the 10 patterns apply to all of them. Each user do at least 5 tasks out of 40, following one of the 10 patterns. My goal is to see whether or not users tend to always follow the same pattern when performing these tasks. We define:   1. P1, P2.... P10 as 10 possible patterns for the tasks.   2. U1, U2, ..... U1000 as 1000 users   3. T1, T2,..... T40 as 40 tasks Since we don't care about Task-differences for now, so my data can be illustrated with P and U only:                user, patterns      U1, P1 P1 P1 P2 P2      U2, P2 P3 P5 P4 P5      U3, P1 P1 P1 P1 P1 P1 P1 P1 P1 P1 P1 P1      U4, P2 P2 P2 P2 P2 P2 P2 P1 P1 P1 P1 P1 P1       ....      U1000, P3 P2 P3 P2 P3 P2      Let's come back to my question: _Do users usually tend to stick to the same pattern?_ I would like to know how to measure this kind of stability. Now what I do is the following: For each user, I compute the number of repetitions of a certain pattern, then pick up the most frequent one. Calculate the proportion of these repetitions with regard to the number of patterns for this user. Then I get something like this:               user, proportionOfMostFrequentPattern      U1, 0.9      U2, 1      U3, 0.4      ....      With the above approach, I can capture how much proportion of time a user followed the same most frequent pattern. However, if we take a look at U4. His patterns have 7 P1 and 6 P2. If I use the approach above, the proportion is 7/13. This does not consider the repetitions of P2. Is there a standard metric that is able to calculate the stability of the patterns? In addition, is there a significance test for this kind of issue?


title: How to generate a bounded random correlation matrix?

text: I'd like to know if anyone can suggest a method for generating random correlation matrix with 90% of the entries (except 1) between $[-0.3, 0.3]$. The other 10% should be larger than 0.3 or smaller than -0.3. Thanks!


title: How to correlate attitudes and behavior?

text: I am about to turn in a paper on research methods and need some guidance on how to analyze my data statistically. The paper is about a hypnotized study where participants fill in a Likert type attitude questionnaire concerning certain types of behavior (cooperation, communication, leadership etc.). The participants will then solve a team task whilst being observed and scored by an observer on those same types of behavior. I wish to come out of this with a correlation between attitudes and behavior, but I am not sure how to do this. Can I plug it straight into a Pearson's r? or must I use Spearman rank order? What if the results does not allow for a parametric test? Also any recommendations on how to analyze inter-rater validity would be great! would Cronbach's alpha be good? Please be gentle in your suggestions, as I am a total novice when it comes to statistics:-) Thanks! Thanks for replies! Good karma to you! For simplicity I will presume there are no outliers. The Likert scale will be 1-5 for both attitudes and behavioral scores. There will be several attitudes questions and behavioral categories. Does that mean the data is on interval level? If so does that restrict me to Spearman's rank order?


title: probability distribution for heavy zero, right skewed data

text: My question is almost identical to this one. Im searching for a probability distribution for skewed data that allows for zeros. The purpose is to fit a GLM model. My data on species-area distributions are heavy right skewed. Additional variables are "partitioned area" according to latitude- so area in tropical and temperate zones, thus the data now includes zeros. So my question is; is there a probability distribution, e.g. similar to the log-normal, but which allows for zeros? Any pointers would be highly appreciated, thanks! **UPDATE** Histograms of my area data ![enter image description here](http://i.stack.imgur.com/eNLv0.jpg)


title: Plausible interpretation of PCA charts

text: ![enter image description here](http://i.stack.imgur.com/0qA0m.png) Hi, there.. what could be a plausible explanation to this plot, where the x and y axes are the first and second principal components respectively?


title: Why the sample mean is not a good (sufficient) statistic?

text: Consider a simple example of $X_{i}$ be i.i.d uniform distribution on the interval $[\theta,\theta+1]$. By strong low of large numbers, I may conclude that $$\overline{X}\rightarrow_{P} \theta+\frac{1}{2} $$ However, it is not so clear what is the pdf or cdf of $\overline{X}$ as I would have to do an $n$-dimensional integral to get the cdf as $F_{\overline{X}}(t)=P(\sum X_{i}<nt,\theta\le X_{i}<\theta+1)$. Nevertheless it seems "intuitive" that $\overline{X}$ should be a good statistic. But practice one soon learns that the real sufficient statistic is $(X_{(1)},X_{(n)})$, and proving this via the factorization theorem is not difficult. I decided to ask the professor about this after the class. He told me that $\overline{X}$ is not a good statistic because it is not close enough to $\theta$, and as statisticans one has to consider real life applications. But this explanation is not persuasive, since I can use $\overline{X}-\frac{1}{2}$ for the same purpose. I want to ask if $\overline{X}$ is really a bad statistic for this example, and if yes for what reason. Also I want to know if $\overline{X}$ is a sufficient statistic.


title: Sample size to compare ' change from baseline' proportions between two groups?

text: I have frequencies for variable with answers Yes and NO for two time points (at baseline and after treatment). The frequencies are recorded for two treatments. I would like to compare the change from baseline between these two treatments. How can I calculate the sample size? I have the following frequencies for “Yes”: treatment 1 : 28% at baseline   treatment 1 : 41% after treatment, and treatment 2 : 21% at baseline   treatment 2 : 39% after treatment.


title: Use logarithmic or linear values when performing correlation and ANOVA?

text: The data was collected using a logarithmic scale: scores (explanatory variable) are collected with values from 0-5. A score of 3 is 10x the value of a score of 2. I wasn't sure if using either the log scale (0-5) or a "linearalized" scale (0-100000), when performing corr. or ANOVA, was correct, so I tried both. Unsurprisingly, the scale I choose affects whether correlations are found or not, and whether variances are found or not. Which scale should I use in my analysis? Or what questions do I need to answer to help me decide which scale I should use in my analysis? Thanks very much!


title: Is this a normal distribution?

text: I have done some processing to my original data. Before the processing I had a nice histogram with a bell shape. Now it still looks like a bell shape to me if I would slice off the first half and stick it to the end. But I'm not quite sure what happened here. Is this still a normal distribution? Or else, what kind of distribution is it? For my task (motion recognition) I am supposed to use normally distributed data and I'm not sure whether I can use the data if the distribution looks like this. Also, I think that if I have a normal distribution at the beginning, the data should still be normally distributed after I did my processing steps (because these only included normalization and subtracting the mean from each data point and then normalizing again). ![](http://i.stack.imgur.com/34CUK.jpg) I don't know much about distributions, so if you have any tips for me where I could read about that I'd be very happy, I couldn't find any information elsewhere.


title: To log or not to log?

text: I am working on industrial output data and running a lot of tests for it, before I proceed to make an ARIMA model for it. Before I do that, I need to decide whether or not to "log" the data. I am having trouble spotting whether the data is growing exponentially or not. And even if it's not - it's logical to assume that industrial growth has constant growth and its variance goes up with time. ![Original Data](http://imgur.com/QyTgcnH.jpg) Bottom line - would you advise to log the data in general (not just for ARIMA) when I do things like seasonally adjust it and smooth it, etc.


title: Fligner test in R, with several variables indicating the grouping

text: I have to run an multiple way Anova and want to test the assumption of homoscedasticity. How can I run a Fligner.test with several independent variables (variables indicating the grouping) with R? Does it make sense to do such a thing?


title: Double integration diverges due to the non-Gaussian noise

text: I added Gaussian noise to my input data and then I integrated it twice (I used trapezoidal rule). I was wondering if the the integration itself transform the Gaussian noise to something else. I looked at the probability density function of which is $$f(x)= \frac {1} {\sigma \sqrt{2 \pi}} e^{- \frac{1}{2}(\frac{x- \mu}{\sigma})^2}$$ When we integrate the probability density function it's equal to 1 (property of Gaussian distribution). $$\int f(x)=1$$ Thus the distribution of the noise is Gaussian. When we integrated again (double integration) the probability density function is not equal to 1 anymore and it depends on the domain which we integrate it. I was wondering if we can conclude that after the second integration the Gaussian noise is not Gaussian anymore from the fact that the double integration of the probability density function is not 1.


title: Fast way to calculate difference in normal CDFs

text: I'm running a computationally intensive method where I have to calculate the difference in Normal CDF's millions of times, such as               pnorm(y)-pnorm(x)      I have not looked into the details of how the CDF is calculated in R but I was thinking there may be a way to speed this up, perhaps into one pnorm call. Or maybe there a mathematical simplification I can take advantage of. Any suggestions (maybe it's not possible)?


title: Time Series Analysis for Non-Uniform Data

text: I have wait-time data for a coffee shop for 4 weeks. Since the data is crowdsourced, it is sparse and non-uniform along the time. So my question is: how should I deal with this non-uniform data? What are the some methods in time series forecasting that can handle non-uniform data? Or is there any way to make the data uniform?


title: Repeated measures ANOVA with lme in R for two within-subject factors

text: I'm trying to use `lme` from the `nlme` package to replicate results from `aov` for repeated measures ANOVAs. I've done this for a single-factor repeated measures experiment and for a two-factor experiment with one between- subjects factor and one within-subjects factor, but I'm having trouble doing it for a two-factor experiment with two within-subjects factors. An example is shown below. `A` and `B` are fixed-effect factors and `subject` is a random-effect factor.               set.seed(1)     d <- data.frame(         Y = rnorm(48),         subject = factor(rep(1:12, 4)),         A = factor(rep(1:2, each=24)),         B = factor(rep(rep(1:2, each=12), 2)))          summary(aov(Y ~ A*B + Error(subject/(A*B)), data=d))  # Standard repeated measures ANOVA          library(nlme)     # Attempts:     anova(lme(Y ~ A*B, data=d, random = ~ 1 | subject))  # not same as above     anova(lme(Y ~ A*B, data=d, random = ~ 1 | subject/(A+B)))  # gives error      I could not see an explanation of this in the Pinheiro and Bates book, but I may have overlooked it. Could anyone tell me how to do this? Thank you.


title: The ACF of a unit root process

text: If the correlogram of $Y_t$ displays slow decay (i.e., non-stationarity), is this indicative of a unit root? I reasoned that for the model $Y_t = Y_{t-1} + u_t$, the ACF is equal to 1 for all lags, and therefore, the answer to my question is yes. But I just want to be sure.


title: How to use rejection sampling to generate draws from Unit Exponential

text: I'm working on some practice test problems, and one of them says to design a rejection sampling algorithm to produce draws from a unit exponential using draws from a Gamma(2,1). I don't understand how this is possible, because I am under the impression that the "envelope function" g(x) needs to be scalable in such a manner that for some constant $M$, $Mg(x)\geq f(x)\; \forall x$. I can't see any way to do this, as the Gamma(2,1) is going to have little mass around 0, while the exponential function has most of its mass around 0. What kind of transformation do I need to do to the Gamma function to allow it to function as an envelope? Using R, I tried flipping it to make it an inverse gamma, but that won't adequately capture the probability mass close to 0 and a K-S test confirmed that the points I generated did not arise from a unit exponential. edit: I will include my code in which I tried to use an inverse-gamma(2,1) as an envelope:               x <- c()     for(i in 1:100000)     {       g <- runif(1, 0, 1)        h <- rigamma(1, 2, 1)       M <- densigamma(h, 2, 1)       crit <- dexp(h, 1)/M       if(g < crit)         x[i] = h     }


What would be $\mathrm{Var}(X^2)$, if $\mathrm{Var}(X)=\sigma^2$?


title: Bootstrapping data envelopment analysis efficiency score using R

text: I want to perform bootstrapping for calculation of efficiency score from data envelopment analysis (DEA) using R.   * Are there any examples of data and results for this type of analysis in R to enable me to check my results?   * Are there any online or other resources that might assist my task?


title: Putting variables into the residual category

text: What does it mean to put control variables together into the residual category? Does that mean leaving them out of the model?


title: Statistical Analysis on Comments and Feedback

text: I have access to 10,000 comments for a mobile app, and I want to run some interesting statistical analysis on them. What I have done so far: Look at the frequency of each word in all the comments. Then look at a subset of these words that are relevant for performance feedback (ex: "crash", "freeze", "slow") and simply note the frequency of appearance. Would you guys suggest any other methods to extract relevant information from these comments? I have access to excel and SQL.


title: Finding which customers I am losing and gaining

text: I work in an oil analysis lab and I want to devise a statistical test to determine which of our customers are buying more oil sample tests and which customers are buying less than usual. The data I have is a list of the number of purchases on a given day. Whats a good statistical test to determine an increase or decrease in the number of samples purchased?


title: Locomotive problem with various size companies

text: I'm working through Think Bayes (free here: http://www.greenteapress.com/thinkbayes/) and I'm on exercise 3.1. Here's a summary of the problem: "A railroad numbers its locomotives in order 1..N. One day you see a locomotive with the number 60. Estimate how many locomotives the railroad has." This solution is found with the likelihood function and exponential prior like so:               class Train(Suite):       def __init__(self, hypos, alpha=1.0):         # Create an exponential prior         Pmf.__init__(self)         for hypo in hypos:           self.Set(hypo, hypo**(-alpha))         self.Normalize()       def Likelihood(self, data, hypo):         if hypo < data:           return 0         else:           return (1.0/hypo)      Conceptually this is saying, if we see a train number larger than one of our hypotheses (1...1000) then every hypothesis that's smaller has a zero chance of being correct. The rest of the hypotheses have a 1/number_of_trains chance of showing us a train with this number. In the exercise I'm working on the author then adds on a little extra. This assumes there's only one company. In real life however you'd have a mixture of big and small companies and bigger companies (both equally likely). However, this would mean that you're more likely to see a train from a bigger company since they'd have more trains. Now the question is how to reflect this in the likelihood function? This isn't Stack Overflow so I'm not really asking for coding help, but instead perhaps just help about how I might think about this problem in terms of a likelihood function.


title: Estimating the minimum of a finite probability mass function

text: Suppose we are given a discrete r.v. $X$, distributed according to some _unknown_ , finite probability mass function $p(x)$. We can assume that $p(x)>0$ for every $x$ in its domain. We can sample $X$ i.i.d. and obtain a sample $S= (X_1,...,X_m)$. Problem: how can we _formally_ lower-bound $p_*=\min_x p(x)$, within a certain confidence? To be more specific, I'd like an estimator $\hat p_*(S)$ for which it is possible to prove something like this, for any $1\geq \lambda\geq 0$: $\Pr(|p_*- \hat p_*(S)|>\lambda)\leq f(\lambda,m)$ for some rapidly decreasing function $f$ of $\lambda$ and $m$. For example, $f(\lambda,m)=\exp(-c\lambda^2m)$, for some constant $c>0$. Any reference to the relevant literature would be greatly appreciated.


title: Good introductions to time series (with R)

text: I am currently collecting data for an experiment into psychosocial characteristics associated with the experience of pain. As part of this, I am collecting GSR and BP measurements electronically from my participants, along with various self-report and implicit measures. I have a psychological background and am comfortable with factor analysis, linear models and experimental analysis. My question is what are good (preferably free) resources available for learning about time series analysis. I am a total newb when it comes to this area, so any help would be greatly appreciated. I have some pilot data to practice on, but would like to have my analysis plan worked out in detail before I finish collected data. If the provided references were also R related, that would be wonderful. Edited: to change grammar and to add 'self report and implicit measures'


title: Steps of data analysis using logistic regression

text: I am a very beginner of statistic. Recently a project require me to analyse data using logistic regression & SPSS within a specific time frame. Although I have read few books, but still very blur on how to start off. Can someone guide me through? What is the 1st ste and what next? Anyway, I have started some. Once entered the data into SPSS, I have done crosstab (categorical IV), descriptive (continuous IV) and spearman correlation. Then, I proceed to test for nonlinearity by transforming into Ln which give me some problems. I have re-coded all zero cells to a small value (0.0001) to enable the Ln transformation. Then, I re-test the nonlinearity. Question: 1) The only solution for violation is to transform the variable from continuous to categorical? I got one violation. 2) One Exp(B) is extremely large (15203.835). What does this means? Why? 3) There is one interaction has Exp(B) = 0.00. Why? Many thanks.


title: Manifold regularization using laplacian graph in SVM

text: I'm trying implement Manifold Regularization in Support Vector Machines (SVMs) in Matlab. I'm following the instructions in the paper by Belkin et al.(2006), there's the equation in it: $f^{*} = argmin_{f \in H_k}\sum_{i=1}^{l}V(x_i,y_i,f)+\gamma_{A}\parallel f \parallel_{A}^{2}+\gamma_{I}\parallel f \parallel_{I}^{2}$ where V is some loss function and $\gamma_A$ is the weight of the norm of the function in the RHKS (or ambient norm), the enforces a smoothness condition on the possible solutions, and $\gamma_I$ is the weight of the norm of the function in the low dimensional manifold (or intrinsic norm), that enforces smoothless along the sampled M. The ambient regularizer makes the problem well-posed, and its presence can be really helpful from a practical point of view when the manifold assumption holds at a lesser degree. It has been shown in Belkin et al. (2006) that $f^*$ admits an expansion in terms of $n$ points of S, $f^*(x)=\sum_{i=1}^{n}\alpha_i^*k(x_i,x)$ The decision function that discriminates between class +1 and -1 is $y(x)=sign(f^*(x))$. The problem here is, I'm trying to train SVM using LIBSVM in MATLAB but I don't want to modify the original code, so I have found the precomputed version of the LIBSVM which instead of taking input data, and output groups as parametes, gets Kernal matrix computed and the output groups and trains the SVM model. I'm trying the feed it with the regularized Kernel matrix (Gram Matrix) and let it do the rest. I tried to find the formula which regularizes the Kernal and came to this: Defining $I$ as the identity matrix with the same dimension as the Kernel Matrix, $K$ $G=\frac{2\gamma_AI + 2\gamma_ILK}{I}$ $Gram = KG$ In which $L$ is the Laplacian Graph Matrix, $K$ is the Kernel Matrix and $I$ is the identity matrix. And $Gram$ is computed using inner multiplication of two matrices $K$ and $G$. Is there anyone who can help me figure out how this is computed?


title: Is it possible for  $R^2$ of a regression on two variables be higher than the sum of $R^2$ for two regressions on the individual variables?

text: In OLS, is it possible for the $R^2$ of a regression on two variables be higher than the sum of $R^2$ for two regressions on the individual variables. $R^2(Y \sim A + B) > R^2(Y \sim A) + R^2(Y \sim B) $ Edit: Ugh, this is trivial; that's what I get for trying to problems issues that I thought of while at the gym. Sorry for wasting time again. The answer is clearly yes. $ Y \sim N(0,1)$ $ A \sim N(0,1)$ $ B = Y - A $ $R^2(Y \sim A + B) = 1$, clearly. But $R^2(Y \sim A)$ should be 0 in the limit and $R^2 (Y \sim B)$ should be 0.5 in the limit.


title: Fisher's Exact Test with weights?

text: Does anyone know of a variation of Fisher's Exact Test which takes weights into account? For instance sampling weights.   So instead of the usual 2x2 cross table, every data point has a "mass" or "size" value weighing the point. Example data:               A B weight     N N 1     N N 3     Y N 1     Y N 2     N Y 6     N Y 7     Y Y 1     Y Y 2     Y Y 3     Y Y 4      Fisher's Exact Test then uses this 2x2 cross table:               A\B  N  Y All      N   2  2   4      Y   2  4   6     All  4  6  10      If we would take the weight as an 'actual' number of data points, this would result in:               A\B  N  Y All      N   4 13  17      Y   3 10  13     All  7 23  30      But that would result in much too high a confidence. One data point changing from N/Y to N/N would make a very large difference in the statistic.   Plus, it wouldn't work if any weight contained fractions.


title: Elasticities Using GLM

text: The coefficient on a logged explanatory variable when the dependent variable also is in log form is an elasticity (or the percentage change in the dependent variable if the explanatory variable changes by one percent). Suppose I estimate a regression without logging the dependent variable but I use a log link in a General Linear Model (and family gaussian) while the explanatory variable remains in log form. Is the coefficient on that explanatory variable still an elasticity?


title: ANOVA with non-independent observations

text: Sorry for the verbose background to this question: Occasionally in investigations of animal behaviour, an experimenter is interested in the amount of time that a subject spends in different, pre- defined zones in a test apparatus. I've often seen this sort of data analyzed using ANOVA; however, I have never been entirely convinced of the validity of such analyses, given that ANOVA assumes the observations are independent, and they never actually are independent in these analyses (since more time spent in one zone means that less is spent in other zones!). For example, > D. R. Smith, C. D. Striplin, A. M. Geller, R. B. Mailman, J. Drago, C. P. > Lawler, M. Gallagher, Behavioural assessment of mice lacking D1A dopamine > receptors, Neuroscience, Volume 86, Issue 1, 21 May 1998, Pages 135-146 In the above article, they reduce the degrees of freedom by 1 in order to compensate for the non-independence. However, I am not sure how such a manipulation can actually ameliorate this violation of ANOVA assumptions. Perhaps a chi-squared procedure might be more appropriate? What would you do to analyze data like this (preference for zones, based on time spent in zones)? Thanks!


title: Machine Learning techniques for classifying response from several continuous factors

text: I have some data with rows corresponding to male-female pairs, with 8-10 continuous factors (that have been unfortunately anonymized) that apply to both the male and female, and whether they became friends in a week or so. An example row would be like:               m_factor1    m_factor2    f_factor1    f_factor2    Became_friends?     .5           .9           -.4          1            1      (Is this what is considered 'dyadic' data?) I went about building a logistic model based on this factors by using the standard procedure of running the regression and finding significant factors, keeping those, and then mixing and matching the factors, running it again, etc. I find when doing cross- validation my model predicts outcomes decently well, agreeing with the actual result about 60% of the time. However, given the problem context, I feel like there should be significant interaction between certain factors, specifically between the same factors ('m_factor1' vs. 'f_factor1') that I have not included in my model. What techniques could I use to include\find these interaction effects? I don't know much about 'machine learning' techniques beyond classical regression analysis.


title: Using mixed variables in two-step cluster analysis

text: I need to perform cluster analysis on my data set which includes both continuous and categorical type of variables. Having read around, I think K-Means is not a suitable technique for mixed data. Can I use Two-Step Clustering for the analysis? If yes, how are the categorical and continuous variables treated, technically? Thanks a lot for the help!


title: Best fit k1,k2,k3 for w[i]-w[i-1]=k1*(c[i]-(k2*w[i-1]+k3)) with c[i], w[i] known

text: At the end of each day "i", I record my weight "w[i]", and the number of calories I've consumed that day "c[i]". I believe my weight change "w[i]-w[i-1]" is proportional to my net calorie intake: "w[i]-w[i-1] = k1*(net calorie intake)". The number of calories I burn is a linear function of my prior day's weight: "k2*w[i-1]+k3" My net calorie intake is thus "c[i]-(k2*w[i-1]+k3)". My weight change is thus: "w[i]-w[i-1] = k1*(c[i]-(k2*w[i-1]+k3))" Given that I know "w[i]" and "c[i]" for each i, how do I find the best fit values of k1, k2, and k3? Note: I'm trying to solve this as a data-fitting problem. I realize it doesn't map reality perfectly.


title: What is the global model?

text: I'm trying to calculate c-hat, the overdispersion parameter for a QAIC model set. According to Burnham and Anderson, you're supposed to calculate c-hat on the global model. Is the global model the model from the a priori set that has the most parameters, or is it a separate model not in the candidate set that includes a combination of all parameters used in the candidate set? For example, with this candidate set of models:               1 grass+trees     2 grass*trees     3 grass*trees^2     4 shrubs     5 river      Is the global model simply model #3, or is i: grass*trees^2+shrubs+river ?


title: Simulated Microarray data

text: I would like to generate synthetic microarray data sets for simulation purposes. The web and literature search so far returned the following:   1. Microarray Simulator \- It provides Matlab skripts to simulates all the pipeline before the actual values. For my needs it is overpowered - and I also would like to have some methods in R.   2. Simulated two-sample microarray data from OCplus \- is actually exactly what I'm looking for. It provides methods for producing synthetic two class data and it is in R! But I'm interested in further methods for creating synthetic/simulated microarray data. Especially for more than two classes. And maybe more skrews and bolts to tinker with. I tried to find some methods in `limma`, but I wasn't succesful. Does anybody know or can recommend an R package/library for creating simulated/synthetic microarray data set ? Thanks a lot!


title: how to interpret reading decision tree result from ctree() in r?

text: After running ctree(model) My result look like this below :               1) iv == (a,b,c) ; criterion = 1 , statistic = 200          2) iv2 == (x,y) ; criterion = 0.99 , statistic = 32              3) iv3 < 30 ; criterion = 0.91 , statistic = 35                  4) weight = 100      I am new to this , but can anyone please explain to me say my total volume is 500 I know that that node (4 ) given me volume of 100 , but what exactly is criterion = 1 and statistic = 200 means ? Thank you


title: Family-wise error of dependent tests?

text: Say you have a drug of which you want to test wether it increases the number of immune cells in your blood. You divide your sample in two groups where one group receives the treatment and the other group doesn't. In each group, you measure the effect both before treatment initialization and after 2 months of treatment. One can now test the effect of the drug in two ways:   1. Test the number of immune cells at two months in both groups   2. Test the increase in immune cells (=effect after 2 months minus baseline) in each group. If there would be no effect of the drug at all, but you would do both tests at the 5% level and pick the smallest p-value, what would then be the family-wise error rate (FWE)? As you can see, both tests are not independent, however, they do generate different p-values, so I already reasoned the FWE should be somewhere between 5 and 9.75%. Thanks a lot!


title: Expectation of reciprocal of a variable

text: I am confused in applying expectation in denominator. E(1/X)=? can It be 1/E(X)?


title: Frequentist confidence intervals = constant trapping probability?

text: In the case of estimating an unknown mean of a normal distribution with known variance, if I'm not mistaken, the confidence interval contains $\theta$ with probability $1 - \alpha$, regardless of the actual value of $\theta$. In other cases (e.g. when the variance is not necessarily constant), is it still the case that the confidence interval contains the actual value with probability $1 - \alpha$? Or is it only the case that for any value of $\theta$, the confidence interval contains the actual value with probability **at least** $1 - \alpha$? (which if I'm not mistaken is how confidence intervals are defined). Thanks!


title: kNN search using distance fitted to a training set of similar pairs

text: I want to perform k-nearest neighbor search in multidimensional space but not using for example L2 distance but I want user to specify some "similar" pairs- examples and then perform search using this information. What algorithm I can use for this?


title: I'm looking for solution manual “A first course in Bayesian statistical methods”

text: Im looking for a solution manual for Peter Hoff's A first course in Bayesian statistical methods. I cannot find it online, does anybody know whether there is a manual available? Alternatively does anybody know where I can find good exercises and solutions for a math student who tries to understand Bayesian Statistics? Especially   * Belief, probability and exchangeability   * Conjugacy in exponential family models   * Monte Carlo approximation   * Gibbs and Metropolis sampling   * Bayesian linear regression


title: How to use CART for AdaBoost?

text: I am trying to use CARTs (Classification and Regression Trees) for AdaBoost as weak learner. My question concerns the update of the weights after fitting the best weak learner. A single CART node consists of a simple threshold (and of course the selected feature descriptor) which separates the data set for the left successor and/or right successor node. Let us assume that I want to update the weight of a single data point $x_i$ by using a single CART with depth 3.                                                                root                                               /                \                                             n1(-1)              p1(+1)                                           /      \           /      \                                        n2(-1)   p2 (+1)     n3(-1)    p3(+1)      Must I use all traversed nodes for the update of the weight or do I only need the final output of the tree for the update? I saw the first approach in the GML AdaBoost MATLAB Toolbox. I hope I could describe my problem sufficiently.


title: 3 groups of 8 with a different treatment each, over 45 minutes - confused as to stats test needed

text: I have three experimental groups. Each will provide a 15 minute baseline before an addition of either drug A, drug B or a saline control. Each group has an N of 8 as maximum, but possibly as low as 6 if data has to be excluded due to noise. As far as I understand it, a two way anova would be best. Possibly repeat measure as each time point will have 8 replicates. However it seems to me that the two way anova tests to see if time has significant interaction which I already know it will do. As during the baseline period I expect a decline till it is flat. I am only interested in testing for a drug dependent change and I do see a large change. Please could you advise whether I am using the best test for use in this scenario.


title: First-order Markov chain, states transition probabilities at each time are enough for the model

text: Hi markov chain specialist, Hope you can give me an answer regarding this trellis diagram that i saw on a book. why in this picture of a **general** first-order markov chain of 2 states,we should know the prob. of each state at each time? A general first-order markov chain, can be Time-dependent(non stationary) so Transition prob. can change by time. That's why we have E(time) in this picture at each time. but by having the **initial** Prob. of each state (P(t=1)), and **Transition** prob. i can calculate prob. of each **state** in any given time, so The state probabilities at each time (P(t=2,t=3,...)) are **Extra** informations. am i right or i ignore something? ![enter image description here](http://i.stack.imgur.com/SUHMH.png) pic.Statistics in Volcanology (google books) page 167 MArkov chain


title: How do I calculate the difference between incidences?

text: I'm making a study on the necessity of routine omentectomy in ovarian cancer patients. The plan is to first find the incidence of occult metastases of the omentum in ovarian cancer patients, and afterwards I would like to compare my study's incidence to other published incidences... In other words, I want to know if it is statistically different or not to the other published incidences. What method should I use to calculate that? Thank you v much.


title: Mixed effect and learning curves

text: My question might be too simple but I just started to do statistical analysis and use R and is not always simple! I have performance data for 5 subjects with different level of experience, repeating the same task 6 times. I would like to predict the time to perform the task with mixed effect models. I have been building a model including the time as a response, the Trial number, the previous experience and their interaction as fixed effects and the subjects as random effects. Would this be correct? I have been then showing the learning pattern plotting my data and adding a smooth line (lowess). This is the first model               > ma <- lme(prod.time~Trial+Previous.experience+Trial*Previous.experience,                  data=  Data27_04, random=~1|Student, method="ML")     > summary(ma)     Linear mixed-effects model fit by maximum likelihood     Data: Data27_04             AIC      BIC    logLik       1517.445 1541.259 -750.7226          Random effects:     Formula: ~1 | Student             (Intercept) Residual     StdDev:    7.337648 42.42332          Fixed effects: prod.time ~ Trial + Previous.experience + Trial * Previous.experience                                     Value Std.Error  DF   t-value p-value     (Intercept)                102.44173  9.561987 137 10.713435  0.0000     Trial                       -6.48494  2.252271 137 -2.879291  0.0046     Previous.experience1       -37.36173 14.786033   2 -2.526826  0.1274     Previous.experience2        47.22627 12.451072   2  3.792948  0.0630     Trial:Previous.experience1   6.55351  3.496401 137  1.874360  0.0630     Trial:Previous.experience2  -7.55163  2.940879 137 -2.567813  0.0113     Correlation:                                 (Intr) Trial  Prvs.1 Prvs.2 Tr:P.1     Trial                      -0.841                                 Previous.experience1        0.253 -0.208                          Previous.experience2       -0.234  0.199 -0.540                   Trial:Previous.experience1 -0.207  0.264 -0.835  0.447            Trial:Previous.experience2  0.199 -0.226  0.447 -0.836 -0.550          Standardized Within-Group Residuals:            Min         Q1        Med         Q3        Max      -2.3519731 -0.6903211 -0.1031114  0.6503216  4.6699702           Number of Observations: 145     Number of Groups: 5      >      Do you think this is good enough to demonstrate a learning effect? Learning curves are exponential. I have been trying to log transform the response variable but then p-values are saying that previous experience has no significance.               >  mb<-lme(log.prodtime~Trial+Previous.experience+Trial*Previous.experience,                 data=  Data27_04, random=~1|Student, method="ML")     > summary(mb)     Linear mixed-effects model fit by maximum likelihood     Data: Data27_04             AIC      BIC    logLik       225.1042 248.9181 -104.5521          Random effects:     Formula: ~1 | Student             (Intercept) Residual     StdDev:  0.04484554 0.495812          Fixed effects: log.prodtime ~ Trial + Previous.experience + Trial * Previous.experience                                     Value  Std.Error  DF  t-value p-value     (Intercept)                 4.448206 0.10593072 137 41.99165  0.0000     Trial                      -0.060150 0.02629765 137 -2.28726  0.0237     Previous.experience1       -0.333664 0.16351518   2 -2.04057  0.1781     Previous.experience2        0.368358 0.13776525   2  2.67381  0.1160     Trial:Previous.experience1  0.051714 0.04084708 137  1.26604  0.2076     Trial:Previous.experience2 -0.043036 0.03435150 137 -1.25282  0.2124     Correlation:                                 (Intr) Trial  Prvs.1 Prvs.2 Tr:P.1     Trial                      -0.886                                 Previous.experience1        0.248 -0.221                          Previous.experience2       -0.237  0.209 -0.535                   Trial:Previous.experience1 -0.220  0.266 -0.881  0.473            Trial:Previous.experience2  0.208 -0.225  0.474 -0.883 -0.551          Standardized Within-Group Residuals:            Min         Q1        Med         Q3        Max      -2.7119095 -0.8005032  0.1127388  0.8621127  2.1988560           Number of Observations: 145     Number of Groups: 5      >      The model is surely better (AIC, BIC) also the residuals are looking better but then should I reduce the model leaving only the Trial number?


title: performing logistic regression with imputed variables

text: I am trying to to run a logistic regression (case-control) and the variable of interest is categorical, taking the values 0 to 6. For a subset of individuals, I do not have the exact value (0, 1 .. or 6) but a probability distribution : P(V=0)=0.03 P(V=1)=0.2 .. P(V=6) = 0.3 I need an unbiased estimate of the OR and confidence interval both when I test this variable as an ordinal one or when I am testing all the levels as dummy variables (ORs of levels 1, 2, 3 ,4 ,5 and 6 vs level 0). This looks like I need to introduce the probabilities in the likelihood maximization and I was wondering whether this is already implemented in a R script or any other program.


title: Interpretation of Scree plots and Boruta Outcomes

text: I have 37 features in my dataset. I used Boruta package in R and according to its analysis, all the features are "important" and should be retained. I examined this result of Boruta and found that if I exclude any feature from my dataset, then I am losing accuracy of my Random Forest classifier. However, when I am plotting the scree plot of my dataset in R, then I am getting the "elbow" right after the first component. Does it mean only the first component is enough? I am a newbie in data analysis and statistics. So, I would really appreciate if anyone can explain what is going on here.


title: How two distributions of same data but of different time spans can be related?

text: Suppose that I have a dynamic variable X that changes with time, and suppose that I know it's statistical distribution over a period of the last 10 years, and this dis. telling me that the most probable value for X is for example 90. Now suppose that I have another statistical dis. of the same variable, but only for the last year (that is the it includes the previous dis.), and it tells me that most probable value for X is 100. Obviously having two distributions means I have more information, and clearly the Dis. of last year shows that there is kind of a rising trend, that is why it is 100 not 90, but first Dis. also telling us that it is little bit far from the most probable value, so common sense telling me that most probable value should be somewhere in between 100 and 90, and my question is how to calculate it? P.S I'm interested in the "joined" dis. as a whole, not just the most probable value. If the answer is long, pls just refer to some link or book or how this concept of joining is called in statistics.


title: Residual, linear least squares and quality engineering

text: I wanted to ask this question on stackoverflow but I think it is more suitable here. If I am wrong, please tell me. My question concerns the use in statistics to analyse physical/engineering data. As statistical tools to be used depend on the context and are rarely universal, it is mandatory that I describe my problem which is rather long. ## My question In my case, how should I valuate the residual to consider how good are my tests? Should I simply compare their residual mean squares (RMS)? Or should I compare the square root of the RMSs? Or should I compare the root square of the sum of the square of the residual devided by the residual degrees of freedom? Or anything else? summary: $$ \frac{\sum_i(Y_i-\widehat{Y})^2}{n-2} or \sqrt{\frac{\sum_i(Y_i-\widehat{Y})^2}{n-2}} or \frac{\sqrt{\sum_i(Y_i-\widehat{Y})^2}}{n-2} or \; ELSE \;? $$ With the context you will understand why I don't want to use statistical tools using the mean $\overline{Y}$ and only residuals. ## Settings My real problem is rather complex and concerns analysis of some physics- related test data. So for the sake of the example and to avoid to go into physics (which is not our topic here): let's say I am a quality engineer and I want to check the quality of a set of _M_ samples by tensile test. (I picked tensile test because it is taught early in High school and rather simple test) My test is simple: I pull each of my sample only and strictly in the elastic regime which has **normally** a linear correlation between stress (load _L_ , Newtons) and strain to the _m_ (displacement _D_ , Meters). _m_ is the power factor of the displacement and is unknown. So I have the following law (1): $L=a*D^m $ (1) $1\leq m \leq 2$ note: for _m=1_ , a is the spring stiffness "k", related to its Young's modulus "E", and I have a typical tensile test with a typical sample. It is not important if it is the case here but it could happen. _m_ should be about the same for all my samples but is never strictly the same (you know, experiments...) and is unknown. My instrument saves the displacement and load data points in a $2*n$ matrix, _n_ being the number of data points. I have a total of _M_ samples, hence _M_ tests. ## Objective My goal is to check if each test follow correctly my law (1). If yes, my sample is good. If not, my sample is bad and I can trash it (in the recycling bin of course). When plotted, tests on "bad" samples would show steps or several regimes when $L$ vs $D^m$ for $\forall m$ (it is clearly visible for both cases). "bad" samples simply mean there are structural defaults. For good samples, $L$ vs $D^m$ (with best $m$) is a simple line (minus errors induced by the measuring instrument). ## Algorithm To find _m_ , I use a simple iteration process with a linear least square method with Matlab. For $L=0$, $D=0$ of course (no force, no strain). Here is the simple algorithm (not in perfect matlab language for clarity!) for one test:               [RMS,m,a] = function(D,L) %inputs: Displacement D and load L                               %outputs: RMS, m and a     Vm=1*101 Vector    %m will have a 0.01 precision.     Vr=1*101 Vector     Va=1*101 Vector     for 1<Vm<2         X_m = D.^Vm         Va=X_m\L     %left matrix divide built-in function of matlab:                       %Solution by least square method used.         Vr=sum((L-a*X_m).^2)./(length(L)-2)  %Residual Mean Square calculation                                              %Note: length(L)=n: number of data points     end for     [RMS, RMSPos]=min(Vr) %RMS = Residual Mean Square = min(r):least residual.                            %RMSPos: position in the vector of the least residual     m=Vm(RMSPos)         %Saves the m with the least residual     a=Va(RMSPos)         %Saves the a with the least residual          end function      And I do that for each of my _M_ tests. ## Issue Now I want to define a method to say "these samples are good" or "these samples are bad". I would compare the $M$ RMS I got from the $M$ tests: the ones close to 0 would tell me about the good tests, the ones bigger than 0 would tell me about the bad ones. I could decide of a RMS limit $RMS_{lim}$ where $RMS\leq RMS_{lim}$ would give me the good samples and $RMS > RMS_{lim}$ would give me the bad samples. In some ways, I can also sort the "goodness" of all of my $M$ samples by sorting from the smallest residuals to the biggest ones. This way, I don't need to check the plot of all of my $M$ tests manually and see if I have a beautiful line or an ugly random thingy. While I can definitely compare the RMSs within one test for different '$m$'s, I am wondering if it is both physically and mathematically acceptable to compare the RMSs between the different tests. Why? Because on the one hand my different test solutions have similar but NOT strictly equal max $L$ (or $L_n$), max $D$ (or $D_n$), '$m$'s, '$a$'s and '$n$'s ($n$: number of data points in one test). In top of that RMS is a squared value, so I should probably compare the root square of the RMS between the tests. Particularly if it is below 1... **What to compare then?** * * * ## Side-notes I cannot use the tools using the mean, like for example the coefficient of determination, because it would take into account the distance between each data point and the mean of all of them. The issue? Imagine I decide to rerun one test and "forget" one $(D_i,L_i)$ couple close to the previously calculated mean and pull a bit more my sample, still in the elastic regime, to get a higher maximum Load $L_n$. I keep the same amount of data (same $n$). What happens? Well the RMS will still be about the same as my sample will still be intact (elastic test, remember? :) ). But the coefficient of determination would get closer to one... Which has no sense as my sample didn't become suddenly "better". I tried it experimentally to make it sure. Then: in no way this is a method really used for filtering good devices from bad devices. Don't worry, no airplane will crash or building will self-destroy because a bad engineer (me) used a bad method to test the quality of its components :) . I am just a naive student with too little mathematical/statistical background working on a -rather- complex topic. ## Acknowledgement Thank you! To whomever could read this long text. Wli


title: Is sample kurtosis hopelessly biased?

text: I am looking at the sample kurtosis of a fairly skewed random variable, and the results seem inconsistent. To simply illustrate the problem, I looked at the sample kurtosis of a log-normal RV. In R (which I am slowly learning):               library(moments);           samp_size = 2048;     n_trial = 4096;          kvals <- rep(NA,1,n_trial); #preallocate     for (iii in 1:n_trial) {         kvals[iii] <- kurtosis(exp(rnorm(samp_size)));     }     print(summary(kvals));      The summary I get is                  Min. 1st Qu.  Median    Mean 3rd Qu.    Max.        11.87   28.66   39.32   59.17   61.70 1302.00       According to Wikipedia, the kurtosis for this log-normal RV should be around 114. Clearly the sample kurtosis is biased. Doing some research I found that sample kurtosis is biased for small sample sizes. I used the 'G2' estimator as provided by the `e1071` package in CRAN, and got very similar results for this sample size. **The question** : which of the following characterize what is going on:   1. The standard error of the sample kurtosis is simply very large for this RV (even though the hand-wavey common estimate of the standard error is of order $1/\sqrt{n}$). Alternatively, I used too few samples (2048) in this study.   2. These implementations of sample kurtosis suffer from numerical problems which might be corrected by _e.g._ Terriberry's method (in much the same way that Welford's method gives better results than the naive method for sample variance).   3. I computed the population kurtosis incorrectly. (ouch)   4. Sample kurtosis is inherently biased and you can never fix it for such small sample sizes.


title: libSVM for unbalanced data

text: I'm using libSVM for binary classification and my training data is very unbalanced (-1:90%, +1:10%). According to libSVM's documentation, it's better to set different penalties for positive and negative classes. For example, the SVM problem is: $\min\limits_{w,b,\xi} \frac{1}{2}{\bf w^Tw} + C^+\sum\limits_{y_i=1} \xi_i + C^-\sum\limits_{y_i=-1} \xi_i$ My question is which penalty should be larger and why. Thanks


title: Bivariate random variable with R

text: I am trying to transform this Matlab code into R. My goal is to generate a bivariate random variable with a pre-specified correlation. The code uses the idea of t-copula. I can't figure out how to write this in R. Will appreciate any help. The code is:               n = 1000;      rho = .7;      nu = 1;      T = mvtrnd([1 rho; rho 1], nu, n);      U = tcdf(T,nu);      X = [gaminv(U(:,1),2,1) tinv(U(:,2),5)];


title: How to do meta-regression analysis with SPSS?

text: I have done a meta analysis and heterogeneity is too high. I am working with (even/Total) for experimental and control groups to calculate the Odds Ratio. I have done fixed-effect and random-effect modeling. I now need to use meta- regression via SPSS. Can anyone direct me to a good set of materials to learn how to do this? I have read some links but none has really solved this question. Thanks


title: Sample/Filter nonuniformly distributed values so that result follows a uniform distribution

text: I have a dataset with a nonuniform distribution. I want to sample it so that the result is uniformly distributed. If I know the class of the example, with what probability should I choose it to get a uniform distribution of each class in the resulting sample? Example:                        before:  after:     Class1:     40%     20%     Class2:     30%     20%     Class3:     15%     20%     Class4:     10%     20%     Class5:      5%     20%      It probably need some inverse distribution function, but I can't remember any specifics.../


title: Expectation of $\left(X-M\right)^T\left(X-M\right)\left(X-M\right)^T\left(X-M\right)$

text: If $X=[x_1,x_2,...,x_n]^T$ is an $n$-dimensional random variable and we have $E\left\\{X\right\\} = M = \left[m_1,m_2,...,m_n\right]^T$ $Cov\left\\{X\right\\} = \Sigma = diag\left(\lambda_1,\lambda_2,...,\lambda_n\right)$ how can I express the following expectation in terms of $M$, $\Sigma$, and $n$ (and maybe raw $m_i$'s and $\lambda_i$'s)? $E\left\\{ \left(X-M\right)^T\left(X-M\right)\left(X-M\right)^T\left(X-M\right)\right\\}$ Supposing $x_i$'s are i.i.d and have normal distribution would be acceptable, but are these assumptions necessary? **Update:**   1. I know that $E\left\\{ \left(X-M\right)^T\left(X-M\right)\right\\} = \sum_{i=1}^n \left(\lambda_i\right)$ but don't think this would help in this case.   2. In the section ~~_6.2.3 Cubic Forms_~~ _8.2.4 Quartic Forms_ of Matrix cookbook there is a formula for calculated quadratic expectations like this, but i don't want just a formula to solve it. I think there should be a simple question for this problem because the covariance matrix is diagonalized.


title: Image compression by ML techniques

text: I am interested in recent papers or books on image compression using Machine Learning-based techniques. I know about Jiang's survey on Image compression using neural networks, but this article seems a little outdated (1999). Do you know some more recent resources (articles, books) that treat about this subject? Thanks, Lucian


title: What statistical test is appropriate for correlations within repeated measures?

text: We have $N$ individuals measured repeatedly in time. For each individual, at each measurement, we measure the values of two variables: $A$ and $B$. The number of measurements per individual might not be the same. The $i^{th}$ measurement of the $n^{th}$ individual might be taken at a different time than the $i^{th}$ measurement of the $m^{th}$ individual. How can we compute the correlation between $A$ and $B$ and measure the significance of the correlation ? What statistical test can we employ ? Are there any package in `R` that does this ?


title: Calculating the price elasticity and income elasticity of demand

text: The demand function for air travel between the U.S. and Europe has been estimated to be $$\ln Q = 2.737 - 1.247 \ln P +1.905 \ln I$$ where $Q$ denotes number of passengers (in thousands) per year, $P$ the (average) ticket price and $I$ the U.S. national income. Determine the price elasticity and income elasticity of demand.


title: Partial correlation plot, split by groups SPSS

text: Is it possible to illustrate partial correlation scatter plots for 2 subgroups on the same graph? e.g. I want to make scatter plots of data controlled for age, differentiated by males or females. I've tried doing partial regression plots generated by linear regression analysis, but I can't split it by groups. Options to do it in excel or R would be fine too. Thanks!


title: Presenting marginal effects of logit with fixed effects

text: I have a traditional logit model with a dichotomous dependent variable and several independent variables. One of the IVs is categorical and in my R code, I treat it as a factor. In the past, I've presented marginal effects by creating a simulated dataset in which I hold all variables at their mean or median and then increment the variable of interest over its range. I then use the estimated model to get the predicted probability of a given outcome over the range of the variable of interest. This is then easily plotted to create a nice graphical interpretation of the estimated effect of changes in the variable of interest. But in the fixed effect case, where one variable in the model is a factor, how does one generate an analogous graph? Do I need a separate estimated effect graph for each value of the factor variable or is there someway to estimate an average effect of the changes in the variable of interest over all the categories in factor variable while all other variables are held at their mean?


title: Converting categorical variables to numerical ones

text: I am working on a binary classification task where one of the input variables takes values in {A, B, C}. A lot of classifiers only take numerical attributes as inputs so I was wondering what the best way to do this is. If I just map them to integers like A:1, B:2, C:3, I am worried that this would mean that A is more similar to B and less similar to C (because |1-2|=1 and |1-3|=2) but that is not necessarily true. So won't this method introduce some sort of bias?


title: Deep Learning Networks: Fundamental differences

text: I have read and heard in several places that Deep Learning Networks take considerably longer to train than, say, support vector/kernel machines, random forests or boosting methods, but they can give better performance. My question is, what is fundamentally different about deep learning networks in relation to other learning methods that explains this difference? What is it known about their training complexity? For example, are DNN discriminative (in contrast to the other methods)? It looks like they could be considered generative since Wikipedia mentions: > Once sufficiently many layers have been learned the deep architecture may be > used as a generative model by reproducing the data when sampling down the > model (an "ancestral pass") from the top level feature activations. I have also read on Wikipedia that DNNs can be pre-trained in an **unsupervised manner**. Is this a notable difference in relation to other methods. Is this pre-training step that common? But most importantly, **why do DNNs take so long to train in comparison to other methods**?


title: Doubs on binary classification

text: I am new to machine learning and have few questions regarding how to do classification using Neural Network. I have the Matlab toolbox for it.   I have motion trajectories representing some action. These are multi dimensional information containing the centroid of the pixel obtained after tracking result across successive frames. \- I need to classify it as a 2 class problem : Suppose I have the following training data set for activity A               -0.9 0.2     0.89 0.3      and so on I have data set $D_A$ containing 10 samples of this activity. So, how do I put a label across this and feed into a neaural network. I tried putting 1 across the third dimension of the data A for the entire length of the data. How do I train with 10 samples and classify based on some k fold cross validation? A code snippet will be greatful.