title: Definiton of the distribution of estimators, sampling and simulation methods

text: I have a question regarding the definiton of estimators. In the german wikipedia it says that the distribution is determined by g($X_1,...,X_n$) where by g is the estimator function and $X_1,..X_n$ are the observed random variables, and g is evaluated for all possible samples. Now my question is if the distribution of the estimator is evaluated for a fixed n, i.e. does each sample size has its own estimator distribution? Or is the distribution constructed for variable n? Furtheron I am interested in the following: If I have a dataset of around 200 observations given, and I want to construct the distribution of an estimator for different sample sizes, is it recommended to use sampling without replacement or with replacement (i.e. bootstrapping)? The data is normally distributed, and the estimator to be calculated is the estimator for the t-test. Thanks


title: How to test if my data is discrete or continuous?

text: It seems to me that to choose the right statistical tools, I have to firstly identify if my dataset is discrete or continuous. Could you mind to teach me how can I test whether the data is discrete or continuous with R?


title: Same p-value when comparing two GLM

text: This is my first question, please should I write something wrong correct me. I have a question when comparing two GLMs after applying stepwise selection. What I've always heard is that stepwise selection is not robust enough for the selected variables. Therefore, I use ANOVA to compare the most parsimonious model (lowest AIC) and the model without the least significant variable. Accordingly, if the ANOVA outcome is nonsignificant I remove the variable from the stepwise selection. It is said that the simple models (few predictors) dominated over difficult model (a lot of predictors). My problem is that when I try to do it with R. The p-value from the ANOVA is the same as the p-value of the stepwise selection. Do you think that I'm doing something wrong? Here is my code: I first run the command `stepAIC` and obtain the most parsimonious model:               Fullmodel<-glm(Geommean~log_CPUEig+log_CPUE+max_depth+amp_temperature+TotalP+Elevation+max_temperature+lake_area)          stepwise <- stepAIC(fullmodel)                # Step:  AIC=1280.2 Geommean ~ log_CPUE + max_depth + max_temperature     #     #                   Df Deviance    AIC     # <none>                 1474.2 1280.2     # - max_depth        1   1485.6 1280.4     # - max_temperature  1   1495.8 1282.3     # - log_CPUE         1   1937.1 1355.5      Then, I look the p-value for each selected variables from the stepwise procedure:               summary(stepwise)     # Coefficients:     # Estimate Std. Error t value Pr(>|t|)      # (Intercept)      25.3323     3.5104   7.216 5.06e-12 ***     # log_CPUE         -3.5670     0.3811  -9.359  < 2e-16 ***     # max_depth        -0.5818     0.3959  -1.469   0.1428     # max_temperature  -6.1182     3.0290  -2.020   0.0443 *      I choose the `max_depth` as the least significant variable. I run ANOVA to compare the most parsimonious model and the model without the least significant variable.               # Most parsimonious model:     parsimonious <- glm(Geommean ~ log_CPUE+max_temperature+max_depth)           # Model without the least significant variable:     variable <- glm(Geommean ~ log_CPUE+max_temperature)       I compare these two models with ANOVA:               compare <- anova(parsimonious,variable, test="F")     compare     # Resid. Df Resid. Dev Df Deviance      F Pr(>F)     # 1       279     1474.2                               #         280     1485.6 -1   -11.41 2.1593 **0.1428**      After doing ANOVA, I obtain the same p-value of the least significant variable (0.1428). In that case, I removed max_depth from the most parsimonoius model because the p-values is greater than 0.05. Do you think that this outcome is normal? Thank you very much for your help.


title: Methods of moments for t distribution

text: The parameters of a t distribution can be estimated via 1) ML or 2) method of moments If we use the method of moments we have: $\mu=E(R)$ $\sigma^2=V(R)=\frac{\beta \nu}{\nu -2}$ $\kappa = \frac{6}{\nu-4}$ we can rewrite the last two equations: $\beta = \frac{\sigma^2}{(\nu/(\nu-2))}$ and $\nu=\frac{6}{\kappa}$+4 Now my question is, how do I estimate those parameters empirically? I mean, is it ok to do the following? : use r and calculate `mean(data)`, use this value as an estimate for $\mu$. use r and calculate `var(data)` and `kurtosis(data)`, insert these values into the equations and get the values $\beta$ and $\nu$. Is this correct? Are the resulting values the same as using ML?


title: DNA: The number of 'AAAAA'-s in a randomly generated DNA sequence that's 1000 base pairs long

text: Let's say I have a _randomly_ generated sequence consisting of letters A, C, T and G that's 1000 letters long. The probability of each letter occurring is 25%. What is the probability that the sequence 'AAAAA' will occur N times within the 1000-letter sequence? The problem I have solving this is that the trials are dependent, otherwise this would be modeled nicely using the binomial/Poisson distributions. But if the sequence 'AAAAA' happens to occur at position X, then the probability of it occurring at position X + 1 is 0.25 and not 0.25 ^ 5. Thank you.


title: Predictions and confidence limits for linear terms in GAM model

text: I use semi-parametric models with linear and non-linear terms in SAS. I want to plot the predictions and confidence limits of the linear terms in the model. However, the output out statement and plots only show predictions and CL's for the non-linear terms (smoothers) in the model. How can I get the predictions and CL's of the observations of the linear terms in the model?


title: Canonical correlation analysis and time series analysis

text: Is there a way to utilize Canonical Correlation Analysis when your data are time series and repeated measures (i.e. your experimental units are not independent)? How might one approach the analysis of two sets of variables when the question is what relationships, if any, are there between one set of variables and the other. I was thinking canonical correlation analysis might help me do this, but my variables are count data (not normally distributed) taken over several consecutive years at the same location. In sum, one set of variables is the abundances of various species and the other set is the abundances of a variety of potential food resources. Perhaps it's best to look at one dependent variable at a time instead of having several dependent variables. Any advice for a statistics novice?


title: Margin-of-error calculation in survey

text: I've been trying to replicate the results in this online calculator: http://www.raosoft.com/samplesize.html However, it seems that I am missing something. Exactly how do I solve the margin of error when all the other variables (sample size, confidence level, distribution and population*) are known? I've tried to use the formulas on the page but I cannot arrive at the same result. Also using the formulas on bottom of http://www.resolutions.co.nz/sample_sizes.htm gives me somewhat different results. What am I missing? *) = a small population so it cannot be assumed infinite.


title: Gradient of Gaussian log-likelihood

text: I'm trying to find the MAP estimate for a model by gradient descent. My prior is multivariate Gaussian with a known covariance matrix. On a conceptual level, I think I know how to do this, but I was hoping for some help with the details. In particular, if there is an easier way to approach the problem, then that would be especially useful. **Here's what I _think_ I need to do:**   * For each dimension, find the conditional distribution, given my current position in the other dimensions.   * This gives me a local univariate Gaussian in each dimension, with the correct mean and standard deviation.   * I think that the gradient should just be a vector of derivatives for each of these univariate distributions. **My question has two parts:**   1. Is this the best approach to take, or is there an easier way?   2. Assuming I need to go this route, what's the best way to go about finding these conditional distributions?


title: How to find the upper bound of symmetric triangle distribution given median and lower bound

text: I have been told $x$ is drawn from a symmetric triangle distribution. If the lower bound is 3 and the median is 9, how do I calculate the upper bound? I'm not sure which value to use (I have to substitute the value I find into a formula). I was thinking if perhaps I had the upper bound, lower bound and median. I could work out the mean and use the mean value?


title: Conditional FE Binomial Logit Models (Stata)

text: I have individual panel-level data and am using conditional binomial fe logit regressions to estimate the effect of minimum wage on Employment. 1) When I try to include year dummies most of them are excluded due to multicollinearity. I have tried evertying - dropping variables, changing the values of the year dummies - but nothing seems to work. Any advise? 2) When I try to include state dummies all of them are excluded due to no within-group variation. This makes sense, as the survey I use to get the data does not follow movers. Is there still a way to control for state effects?


title: Testing randomly generated data against its intended distribution

text: I have written a program which generates random data. If the program is working correctly, that data should follow a specific, known probability distribution. I would like to run the program, do some calculations on the result, and come up with a p-value. Before anybody else says it: I understand that hypothesis testing cannot detect when the program is operating correctly. It can only detect when it is operating _incorrectly_ in a specific way. (And even then, the test "should" fail X% of the time, depending on what significance level you choose...) So, I am trying to gain an understanding of what tools might be appropriate. In particular:   * I can generate as much random data as I want. All I have to do is leave the program running long enough. So I'm not limited to any specific sample size.   * I'm interested in techniques which produce a p-value. So staring at a graph and saying "yes, that looks kinda linear" is not an interesting option. Unless there's some way of putting a hard number on the "wonkyness" of a graph. ;-) What I know so far:   * I've seen three main sorts of test mentioned which sound like they might be applicable: [Pearson] chi-squared test, Kolmogorov-Smirnov test and Anderson-Darling test.   * It appears that a chi-squared test is appropriate for _discrete_ distributions, while the other two are more appropriate for _continuous_ distributions. (?)   * Various sources hint that the AD test is "better" than the KS test, but fail to go into any further detail. Ultimately, all of these tests presumably detect "different ways" of deviating from the specified null distribution. But I don't really know what the differences are yet... In summary, I'm looking for some kind of general description of where each type of test is most applicable, and what sorts of problems it detects best.


title: Statistical comparison of two signals

text: I need to develop an algorithm that will compare two signals and generate some metric(s) to describe changes between them. Signal processing and analysis isn’t my strong point so I would appreciate any help! Here is an example of what the signals look like: ![Signal example](http://i.stack.imgur.com/znm4r.jpg) The signal starts at a constant amplitude, then transitions into a lower amplitude before transitioning back up to the same constant amplitude as the start. The part in the middle is the region of interest, which is what needs to be analysed. Some of the differences that I am expecting (all within the region of interest) are:   1. Mean amplitude of the area of interest, relative to the start and end amplitudes   2. Slope of fall/rise transitions and of individual peaks and troughs   3. Number of peaks and troughs The algorithm needs to output some generic metrics which can be used to quantify changes in any or all of these parameters. Any guidance on what method(s) I could use to do this would be a great help. I know a little about cross-correlation and have heard about wavelet transform analysis, but I am wondering if these are appropriate, or if there are other methods that would work better.


title: Linear regression using standard deviations as regressors?

text: I am reviewing a study in which standard deviations (SDs) of an X variable, calculated for each individual in the study (measures on each individual were replicated 4 times), are used as predictors of a Y variable. They do not observe a significant correlation between Y and the SDs of X and they conclude by saying that it is not necessary to consider the variability among replicates of the same individual and, therefore, they can work with average values. Does that make sense?


title: DF for regression SS in multiple linear regression when sigma is known

text: In the typical environment for multiple linear regression, we have that $Y = \mathbf{X}\beta + \epsilon$ where $\epsilon$ is iid $\mathcal{N}(0, \sigma^2I) $ where $\sigma^2$ is unknown. In this case, regression sum of squares (SSR) has $\text{df}= p - 1$ ($\text{df}$ = degrees of freedom) where $p$ is the number of parameters in the model. I have two questions based on this.   1. Why is the $\text{df}= p - 1$, it seems to me that all $p$ of the parameters are free to vary.   2. If $\sigma^2$ is known, is the $\text{df}= p - 2$? Edit: I am talking about sum of squares due to regression which is defined as $\sum (\hat{y}_i-\bar{y})^2$ this has degrees of freedom p-1 where p is the number of parameters in the model. SST which is defined as $\sum (y_i-\bar{y})^2$ has degrees of freedom n-1 and SSE (sum of squares due to error/residuals) is defined as $\sum (\hat{y}_i-y_i)^2$ and has degrees of freedom n-p. SST = SSE + SSR and $df_{SST} = df_{SSR}+df_{SSE}$.


title: Overdispersion parameter

text: I am modelling a zero-truncated process with a count model, and am trying to determine whether the data are overdispersed. The Poisson distribution has a variance equal to its mean, $$\newcommand{\Var}{\operatorname{Var}} \Var(y) = E(y) = \lambda $$ The negative binomial model relaxes this assumption by estimating an overdispersion parameter $\alpha$. There are two regimens for this method, `NB1` $$\Var(y) = \lambda \alpha$$ and `NB2` $$\Var(y) = \lambda (1 + \lambda/\alpha)$$ I have seen the second version more often. I have estimated three models with the same covariates but the three different link functions (`Poisson`, `NB1`, `NB2`).               Model   logLik   alpha  p(alpha)     Poisson -7942    1      NA     NB1     -6399    1.001  0     NB2     -7944    403.4  0      Clearly the overdispersion parameter is significant and I should therefore use one of the `NB` models. But how can the log-likelihood for `NB2` be lower than for `Poisson` when there is an additional parameter in the model? When I estimate all three models, I get a warning that my variance-covariance matrix may not be positive semi-definite. But I don't have highly collinear variables, and all of the relevant matrices have strictly positive eigenvalues, so I don't see why I'm getting this warning. Could this be related to the likelihood problem? I am using the `glmmadmb` function from the `R` package of the same name.


title: How is the intercept calculated in a regression model with two independent factor variables?

text: I am struggling with a linear regression model of the shape $y = a + b_1\text{month} + b_2\text{year}$. I have 12 months for each year and 10 years. My dependent variable is a log transformed ratio. I have understood that much that when setting such a model up in R, R automatically picks a level for each variable to go into the intercept. March and 2005 goes into the intercept, in order to provide a baseline for comparison for the other factors. As stated my problem is that i cannot really figure out what the intercept represents. Is it simply an average of the ratio of March and 2005 or what is it?               Residuals:                       Min 1Q  Median  3Q  Max      -0.90339    -0.16789    -0.00373    0.15472 0.88338           Coefficients:                            Estimate    Std. Error  t   value   Pr(>|t|)     (Intercept) -3.586154   0.131642    -27.242 <2.00E-16   ***     MONTHJan    0.381735    0.140731    2.713   0.007875    **     MONTHFeb    0.256457    0.140731    1.822   0.071426    .     MONTHApr    0.072824    0.140731    0.517   0.605981         MONTHMay    0.207984    0.140731    1.478   0.142613         MONTHJun    -0.008194   0.140731    -0.058  0.953686         MONTHJul    0.363693    0.140731    2.584   0.011217    *     MONTHAug    0.195791    0.140731    1.391   0.16727      MONTHSep    0.212562    0.140731    1.51    0.134124         MONTHOct    0.124234    0.140731    0.883   0.379495         MONTHNov    0.204009    0.140731    1.45    0.15032      MONTHDec    0.175348    0.140731    1.246   0.215711         YEAR1999    0.477663    0.128469    3.718   0.000333    ***     YEAR2000    -0.027343   0.128469    -0.213  0.83189      YEAR2001    -0.166637   0.128469    -1.297  0.197612         YEAR2002    -0.060508   0.128469    -0.471  0.638684         YEAR2003    -0.173492   0.128469    -1.35   0.179948         YEAR2004    0.003592    0.128469    0.028   0.977753         YEAR2006    -0.283261   0.128469    -2.205  0.029776    *     YEAR2007    -0.267752   0.128469    -2.084  0.03972 *     YEAR2008    -0.240654   0.128469    -1.873  0.063985    .     ---                      Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1                            Residual standard error: 0.3147 on 99 degrees of freedom                         Multiple R-squared: 0.4167,                      Adjusted R-squared: 0.2988                       F-statistic: 3.536 on 20 and 99 DF,  p-value: 1.491e-05


title: Sample size in genome-wide studies

text: How to calculate sample size needed for GWAS for a given MAF, power, $p$-value and frequency of the disease ?


title: How can a vector of variables represent a hyperplane?

text: I am reading Elements of Statistical Learning and on page 12 (section 2.3) a linear model is notated as: $$\widehat{Y} = X^{T} \widehat{\beta}$$ ...where $X^{T}$ is the transpose of a column vector of the predictors / independent variables / inputs. (It states earlier "all vectors are assumed to be column vectors" so wouldn't this make $X^{T}$ a row vector and $\widehat{\beta}$ a column vector?) Included in $X$ is a "$1$" to be multiplied with the corresponding coefficient giving the (constant) intercept. It goes on to say: > In the $(p + 1)$-dimensional input–output space, $(X,\ \widehat{Y})$ > represents a hyperplane. If the constant is included in $X$, then the > hyperplane includes the origin and is a subspace; if not, it is an affine > set cutting the $Y$-axis at the point $(0,\ \widehat{\beta_0})$. Does "$(X,\ \widehat{Y})$" describe a vector formed by the concatenation of the predictors, the intercept's "$1$" and $\widehat{Y}$? And why does including a "$1$" in $X$ force the hyperplane to pass through the origin, surely that "$1$" is to be multiplied by $\widehat{\beta_0}$? I am failing to understand the book; any help / advice / links to resources would be very much appreciated.


title: Two-way robust ANOVA

text: My struggle with non-parametric methods continues... I'd like to apply a median polish instead of two-way ANOVA (normality and homoscedascity assumptions are violated, and $ n_{ij} $ are small, so I can't use CLT as an excuse). I've never used median polish so far, and our course in statistics taught us to worship ANOVA and forget about robust methods if basic assumptions are not met. I saw this post and it seems that median polish can be applied for two-way factorial design. Which technique do you find appropriate in case of violation of ANOVA assumptions? Now, what are the basic data considerations for median polish (or any other technique you find appropriate in this case)? Same shape, homoscedascity? Any resource (link/reference) is appreciated. * * * P.S. Note that I'm aware of `medpolish` function in R.


title: What is the best machine learning approach for enhancing the weights of features based on common features between classes?

text: I currently have a representation of an object as a graph, and I have some similarity measure on that graph. I currently determine the class of a query object by trying to match its graph representation against every instance in the training set, and I assign the query's class as the class of the best graph matching. Some classes have many similar features that would be good matching nodes, but certain nodes distinguish them from specific classes. For example, nodes `1` and `3` might match well in both classes `class A` and `class B`, but a match of node `2` distinguishes the object as belonging to `class A`. In that case, I would want `node 2` to have a higher weight, as it is more important to discriminating between the two classes. At the same time, nodes `1` and `2` might match well for classes `A` and `C`, but a matching of `node 3` distinguishes the object as `class A`, so `node 3` would be more important and should be weighted as such. This kind of discrimination between classes seems difficult to achieve. What is the best machine learning approach to perform this kind of weighted discrimination? **edit** Adding details. Although two nodes in a graph may match, they are almost never identical. Imagine the graph being different parts of an image, and you have some similarity measure to define how similar two image parts are. Since two images are almost never identical, they will pretty well never be exact matches.


title: Minimum sample size for unpaired t-test

text: Is there a "rule" to determine the minimum sample size required for a t-test to be valid? For example, a comparison needs to be performed between the means of 2 populations. There are 7 data points from one population and only 2 data points from the other. Unfortunately, the experiment is very expensive and time consuming, and obtaining more data is not feasible. Can a t-test be used? Why or why not? Please provide details (the population variances and distributions are not known). If a t-test can not be used, can a non parametric test (Mann Whitney) be used? Why or why not? Thanks a lot for your help! I know very little about stats, and have been searching around the internet for a detailed answer for some time.


title: How to turn a confidence set into a confidence interval

text: Element of Statistical Learning equation 3.15 states that the confidence set is $$ \{\beta: (\hat{\beta}-\beta)^T X^T X (\hat{\beta} - \beta) \leq \sigma^2 \chi^2 \} \>. $$ How do I turn this into a confidence interval?


title: PCA and SVD, R prcomp

text: I am trying to replicate Shalizi's NY Times PCA example found in his Advanced Data Analysis with an Elementary Point of View book. I found sample code and data here http://www.stat.cmu.edu/~cshalizi/490/pca/ I loaded the R workspace pca-examples.Rdata, and was able to extract a CSV out of it which I placed in https://github.com/burakbayramli/classnotes/blob/master/app-math- tr/pca/nytimes.csv In order for PCA using SVD, I did               from pandas import *     import numpy.linalg as lin     nyt = read_csv ("nytimes.csv")     labels = nyt['class.labels']     nyt = nyt.drop(["class.labels"],axis=1)     nyt_demeaned = nyt - nyt.mean(0)     u,s,v = lin.svd(nyt_demeaned.T,full_matrices=False)      Then I tried to project onto the space defined by the eigenvectors               nyt2 = np.dot(nyt_demeaned, u[:,0:2])      Then plot               plot(nyt2[:,0],nyt2[:,1],'.')      and               arts = nyt2[labels == 'art']     music = nyt2[labels == 'music']     plot(arts[:,0],arts[:,1],'r.')     plot(music[:,0],music[:,1],'b.')      I get ![enter image description here](http://i.stack.imgur.com/xiQOm.png) ![enter image description here](http://i.stack.imgur.com/sE7lz.png) This image is different from Shalizi's picture from his book (using the following R code)               load("pca-examples.Rdata")     nyt.pca = prcomp(nyt.frame[,-1])     nyt.latent.sem = nyt.pca$rotation         plot(nyt.pca$x[,1:2],type="n")     points(nyt.pca$x[nyt.frame[,"class.labels"]=="art",1:2],                pch="A",col="red")         points(nyt.pca$x[nyt.frame[,"class.labels"]=="music",1:2],            pch="M",col="blue")      ![enter image description here](http://i.stack.imgur.com/3Q7vK.png) Maybe there was a 90 degree error, flip it one way or another, etc. then you could perhaps have the two images to be somewhat close, but even then the distribution isnt exactly right, and the data points for arts are not as clearly seperated from music data points as they are in Shalizi's picture. Any ideas?


title: Determining if a set of experimental results is statistically significantly different from zero?

text: We have a dataset resulting from an experiment with a genotype A and a genotype B, which underwent either treatment X or no treatment Y, for four total sample groups. Each sample group had only three replicates. The measurements are cell counts. In a two-way ANOVA, we found a significant difference between treatment and non-treatment, and did not find a significant difference between the genotypes. That's all good. The issue is that apparently we expected the non-treatment groups to have values of zero- and they're not. A t-test shows that the difference between groups (A,Y) and (B,Y) to be not significant. I've been requested to determine if the non-zero values of those groups are statistically significantly different from zero or not. While it doesn't seem to make a whole lot of sense to me to do this, I need to provide something, and I'm not sure of an appropriate test that can be done in this case. If I pool all the values for Y, I get a mean of 36, SD of 24 and 95% CI of 19, so the CI doesn't even overlap zero. Is there an acceptable way of answering this statistically, or should I let them know that it's not exactly an appropriate question?


title: Merge covariance matrices

text: I have two 2x2 covariance matrices, stemming from bivariate datasets that are approximately normally distributed. I want to create a mixture distribution and for that I need to merge the covariance matrices. To find the variances, can I apply the following formula? Var(X) =$\Sigma_{i=1}^k\pi_i \sigma_i^2+\Sigma_{i=1}^k \pi_i(\mu_i-\mu)^2$ with $\pi_i$ being the mixture weights. But how do I mix the covariances?


title: Whitening and unwhitening for sparse coding

text: Is this procedure for whitening and unwhitening correct? Given an image $i$: decompose the image in patches:               patch=im2col(i,[8 8],'sliding');      **Whitening step** : 1) subtract the mean from each patch (i.e., from each column):               meanP=mean(patch);     patch_m=patch-repmat(mean(patch),[size(patch,1) 1]);      2) divide by the standard deviation:               devP=sqrt(sum(patch.^2));     devpatch_m=sqrt(sum(patch_m.^2));     patch_m_s=patch_m./repmat(sqrt(sum(patch_m.^2)),[size(patch_m,1) 1]);      **Unwhitening step** : 1) multiply by the standard deviation:               rec_patch= rec_patch .* repmat(devpatch_m, [size(rec_patch,1) 1]);      2) add the mean:               rec_patch=rec_patch+repmat(meanP,[size(rec_patch,1) 1]);      I have also read that before starting the unwhitening step, it is possible to subtract the mean from the "new patches" obtained from the sparse coding step. Which is the correct procedure?


title: Where does the summation sign come from when working out expected value of time

text: > The lifetime of a machine is modelled by an exponential random variable $X$ > with $P(X>x) = e^{-\lambda x}, \lambda, x > 0$. This machine cannot be > repaired. A maintenance crew checks this machine at times $T, 2T, 3T, ...., > $ wjere $T$ is a specified length of time > >   * Find the probability that the machine fails between the $(k-1)$th and > the $k$th inspection, $k = 1, 2, ...$. >   * Find the expected value of the first time when the maintenance crew > finds that the machine is down. [HINT: $\sum_{j = 1}^{\infty} jq^j = (1 - > q)^{-2}q$ if $|q| < 1$.] > I did the first bit and correctly got the probability to be $e^{-\lambda kT}(e^{-\lambda T} - 1)$. But I'm stuck as to what to do from here. I don't see how the summation would come into calculating the expected value? EDIT: > Consider the events $A_k=\{(k-1)T < X \leq kT\}$ and denote $U$ the duration > of the time when the machine is down before it is discovered to be down. > Then use the formulas > > $$EU = \sum_{k=1}^{\infty} E(U |A_k)P(A_k)$$ > > and > > $$E(U|A_k) = \int_0^{T} P(U > u |A_k) du$$ > > to find the expected duration of the time when the machine is down before it > is discovered to be down. I'm still working through this bit, but my first question is, in the first equation, should that read $E(U)$? If not, what does $EU$ mean? EDIT 2: In the answers, it also says $EU$, so what does this mean?


title: Estimating covariance of the difference of directional distributions derived from Gaussian mixtures

text: Given Gaussian mixtures $X_1, X_2 \in \mathbb{R}^p$ defined as $$P(X_i = x) = \sum_s \omega^{(s)}_i \mathcal{N}(x; \mu^{(s)}_i, \Sigma_i)$$ where the superscript $(s)$ indexes the $s$-th component of the mixture distribution $X_i$, the covariance of all mixture components are identical, and $\sum_s \omega^{(s)}_i = 1$ to keep the resulting density normalized. Define new random variables $Z_i$ having support only on the unit sphere in $\mathbb{R}^p$ by marginalizing $X_i$ over the distance from the origin, \begin{equation} P(Z_i = z) = \begin{cases} \int^\infty_0 P(X_i = t z)\,dt & \|z\| = 1 \\ 0 & \|z\| \neq 1 . \end{cases} \end{equation} I am looking to find computationally efficient methods for estimating $\Gamma = \mathbf{E}[(Z_1 - Z_2)(Z_1 - Z_2)^T]$. I have thought of two approaches so far:   1. Crudest possible approach: estimate a Gaussian $Y_i$ that best approximates $Z_i$ (by performing importance sampling with $X_i$ as the proposal distribution to compute $Z_i$'s mean and covariance), then use $\mathbf{E}[(Y_1 - Y_2)(Y_1 - Y_2)^T]$ as the estimate of $\Gamma$.   2. Perform importance sampling directly on $Z_1 - Z_2$ by drawing samples $x_1, x_2$ from $X_1, X_2$, computing the corresponding values of $z_1, z_2$, and the importance weight as $P(Z_1 = z_1, Z_2 = z_2)$.


title: Is there standard nomenclature for Confidence Distributions?

text: I am wondering if there are standard naming conventions in general for confidence distributions or for the particular cases of the non-centrality parameters of (non-central) $t$- and $F$-distributions. For the unitiated, let random variable $X$ have CDF $F(x;\theta)$. The confidence distribution CDF for $\theta$ is just $G(\theta;x) = F(x;\theta)$. This is somehow 'dual' to the original CDF. For creating confidence intervals one typically uses the quantile function of the confidence distribution, here $G^{-1}(q;x)$. So, given data $x$, the interval $\left[G^{-1}(1-\alpha/2;x),G^{-1}(\alpha/2;x)\right]$ is a $1-\alpha$ CI for the unknown $\theta$. Are there standard names for these things? 'dual-t'? 'dual-F'? Also, are there guidelines in R for this kind of thing? For example, `dt, pt, qt, rt` are the density, distribution, quantile and generator for the (non- central) $t$-distribution. What should I call the same functions for the confidence distribution? `dcdt, pcdt, qcdt, rcdt`? If anyone has more knowledge of the history of the $t$- and $F$-distributions, I would gladly name them for pioneering researchers. For example, perhaps the $t$-distribution CD should be name for Johnson and Welch? ( **edited** again for clarity, probably still lacking same. Just looking for the _names_ for these things.)


title: estimating beta parameter

text: I have data on 500 examinees with responses to 20 questions. Because the response is dichotomous, I use Beta(1,1) as the conjugate prior. Now I’m interested in using Beta(alpha, alpha) as the prior. How can I estimate alpha? I hope to do it in R.


title: Testing null hypothesis with dependent samples

text: Let's say I conduct a survey of randomly selected 100 people, asking each of them two questions:   * "Do you like NY?"   * "Do you like LA?". All 100 people responded. Among them, 60 answered "yes" about NY, and 65 answered "yes" about LA. My null hypothesis is that that the fraction of US population who like New York is the same as the fraction who like Los Angeles. How can I test it? I was hoping to use two-sample t-test, but I guess it's not right since my two samples are not independent: the same 100 people are answering both questions.


title: Correlating volume timeseries

text: Consider the following graph: ![twitter and trading volume](http://i.stack.imgur.com/2OrfU.gif) The red line (left axis) describes the trading volume of a certain stock. The blue line (right axis) describes the twitter message volume for that stock. For instance, on May 9 (05-09) about 1.100 million trades and 4.000 tweets were made. I would like to calculate whether there is a correlation between the timeseries, either on the same day or with a lag - for instance: tweet volume correlates with trading volume one day later. I'm reading many articles who have made such analysis, for instance Correlating Financial Time Series with Micro-Blogging Activity, but they do not describe how such an analysis is made in practical terms. The following is stated in the article: ![enter image description here](http://i.stack.imgur.com/q5Qy2.gif) However, I have very little experience with statistical analysis and don't know how to execute this on the series that I have. I use SPSS (also known as PASW) and my question is: what are the steps to take to make such an analysis from the point where I have a datafile underlying the above image? Is such a test a default feature (and what is it called) and/or how could I else execute it? Any help would be greatly appreciated :-)


title: How do I use vector auto regression using the statsmodels library in python?

text: Sorry about the rookie question but I have been at the documentation for three days and couldn't figure much out. (Link To documentation page) First, how do I load my own data? Must I store it in a particular format or can I load it from a csv, mat or xls file. If so how do I do that? And second, I keep getting the following error while using the `model.plot()` function.               model.plot()     AttributeError: 'VAR' object has no attribute 'plot'      (I used the same example as given in the documentation. I couldn't figure out what is wrong.) Did I miss something? Any help would be appreciated!


title: Are variables, which linear combination results in a endogenous variable, endogenous?

text: I'm a little bit confused.... Lets say I assume $x$ to be endogenous and I observe that $x = a v_1 + b v_2 + c v_3$ so that $x$ can be expressed as the perfect linear combination of $v_1, v_2, v_3$. By this I mean that if I regress $x$ onto $v_1,v_2,v_3$ the $R^2$ will be 1 i.e. the model is deterministic. In this particular case is it now wrong to assume that $v_1, v_2, v_3$ are all exogenous? I mean how can a perfect fit of exogenous variables be endogenous. One of the $v$'s has to be endogenous to begin with? In terms of Regression- Analysis I would guess this is not true since we do this kind of thing all the time. The only difference is that we won't get a perfect fit hence my confusion. If this would not be a problem I could always use one of the $v$'s as an instrument for $x$.


title: Bootstrap Confidence Intervals for Weir & Cockerham's Fst

text: I'm working on calculating bootstrap confidence intervals for Weir & Cokerham's Fst. I want to use the percentile-t method as described in this paper. I'm calculating the $F_{st}$ value between two sub-populations in the single allele, mulitple loci case as: $$ \hat \theta_W = \frac {\sum \limits_l a_l}{\sum \limits_l (a_l + b_l + c_l)} $$ as described in Weir & Cockerham's 1984 paper. My current plan is to estimate the variance of $\hat \theta$ by jackknifing over the loci: $$ var(\hat \theta) \mathrel{\hat{=}} \frac {m - 1} {m} \sum \limits_{L = 1}^{m} (\hat \theta_{(L)} - \frac {1} {m} \sum \limits_{L = 1}^{m} \hat \theta_{(L)})^2 $$ where $\hat \theta_{(L)}$ is the estimate of $\hat \theta$ obtained by omitting locus $L$ and $m$ is the number of loci. To get the confidence intervals, my steps are:   1. To perform around 1,000 bootstrap replicates (I haven't worked out how many will be appropriate, I was going to try this as a base, then go up.)   2. In each round, I would gather a simple random sample with replacement of the loci and calculate the $F_{st}$ value ($\hat \theta^*$).   3. I would then using jackknifing to estimate the variance of $\hat \theta^*$. (I would use the same method as above)   4. Then I would calculate the t-statistic: $t^* = \frac {\hat \theta^* - \hat \theta} {\hat se_{\hat \theta^*}}$   5. Once the for loop is done. I would obtain the confidence interval with: $$ (\hat \theta - t^*_{(1 - \alpha / 2)} \hat se_{\hat \theta}; \hat \theta - t^*_{(\alpha / 2)} \hat se_{\hat \theta}) $$ (The confidence interval is from the Wikipedia entry on bootstrapping) My main questions are:   1. Is this the correct way to calculate the studentized bootstrap confidence intervals?   2. Weir and Cockerham suggest that one could also bootstrap over samples instead of loci, the first paper I cited and Weir & Cockerham bootstrap over loci, is this the best approach for comparing two populations? I was going with because my understanding was that we can assume independence between locus (assuming little dependence by linkage) and might not be able between samples in a sub-population as they might be related. Let me know if there are any errors, or if more info is needed. Thank you! Edit: Here's some R code I've written to do this. First I find the a, b, and c arrays and calculate $\hat \theta$, then I do this to find it's variance:               findVar <- function(a, b, c) {       theta.L <- c(rep(0, 10))       for (i in 1:ncol(data)) {         # remove locus i         theta.L[i] <- sum(a[-i]) / sum(a[-i] + b[-i] + c[-i])       }       m <- ncol(data)       mean.theta <- sum(theta.L) / m       var.theta <- ((m - 1) / m) * sum((theta.L - mean.theta)^2)       return (var.theta)     }     var.theta <- findVar(a, b, c)      Then I did this to find the bootstrap confidence intervals:               t.i <- c(rep(0, 10))     m <- ncol(data)     aTemp <- c(rep(0, 10))     bTemp <- c(rep(0, 10))     cTemp <- c(rep(0, 10))     for (i in 1:1000) {      list <- sample(1:m, m, replace = TRUE)       for (j in 1:m) {        aTemp[j] <- a[list[j]]        bTemp[j] <- b[list[j]]        cTemp[j] <- c[list[j]]      }      theta.i <- sum(aTemp) / sum(aTemp + bTemp + cTemp)           # find the variance      var.theta.i <- findVar(aTemp, bTemp, cTemp)           t.i[i] <- (theta.i - theta) / var.theta.i     }          alpha <- 0.05     t.i <- sort(t.i)     firstT <- t.i[ceiling(length(t.i) * (1 - alpha / 2))]     secondT <- t.i[ceiling(length(t.i) * (alpha/2))]          lower <- theta - (firstT * var.theta)     upper <- theta - (secondT * var.theta)      Here's a sample result:               Theta:  -0.1748     Confidence Interval: ( -0.2244 ,  0.6153 )      Note: I didn't have very large data set, just wanted to see if it worked, maybe. Edited again, to fix a mistake with getting the $\alpha / 2$ and $1 - \alpha / 2$ percentiles, and changed the example results


title: Why do we use a one-tailed test [F-test] in analysis of variance (ANOVA)?

text: Can you give the reason for using a one tailed test in the analysis of variance test? Why do we use a one-tail test - the F-test - in ANOVA?


title: Bayesian random walk: updating on samples from posterior

text: Suppose that, at first, I am trying to estimate the mean and standard deviation of some data that I assume to be normally distributed. My prior is gaussian with mean $\mu_0$ and variance $\sigma^2_0$. Now I observe some data $D_0$ and update my prior according to Bayes' rule, leading to a gaussian posterior with mean $\mu_1$ and variance $\sigma^2_1$. Next (and bear with me, this is weird), suppose that I sample from _the posterior_ to create data $D_1$. Again, I update by Bayes' rule, leading to a new posterior ($\mu_3,\sigma^2_3$). Rinse and repeat. This defines a random walk over $u$ and $\sigma^2$. Has this sort of thing been studied? Can anyone point me towards something relevant?


title: When to use gamma GLMs?

text: The gamma distribution can take on a pretty wide range of shapes, and given the link between the mean and the variance through its two parameters, it seems suited to dealing with heteroskedasticity in non-negative data, in a way that log-transformed OLS can't do without either WLS or some sort of heteroskedasticity-consistent VCV estimator. I would use it more for routine non-negative data modeling, but I don't know anyone else that uses it, I haven't learned it in a formal classroom setting, and the literature that I read never uses it. Whenever I Google something like "practical uses of gamma GLM", I come up with advice to use it for waiting times between Poisson events. OK. But that seems restrictive and can't be its only use. Naively, it seems like the gamma GLM is a relatively assumption-light means of modeling non-negative data, given gamma's flexibility. Of course you need to check Q-Q plots and residual plots like any model. But are there any serious drawbacks that I am missing? Beyond communication to people who "just run OLS"?


title: Parametric or non-parametric test?

text: In outline shape analysis after getting the results of chain code (in shape software) and principal component values should I use the parametric (ANOVA) or non parametric (Kruskal-Wallis) test to evaluate the variation between populations ($n=40$)?


title: R - Calculate trendline and extend it for 365 days automatically

text: I have daily data with logarithmic decreasing trend (around 30-days-data). ![enter image description here](http://i.stack.imgur.com/kw13r.png) In excel I'm calculating these data's trendline and with that trendline equation, I'm projecting these values to 365-days data. Then find the sumation of these 365 days' values. Is there any ways to automize that in R by writing a function or using any package ? Thanks.


title: Is there anything similar to Central Limit Theorem for Variance for Hypothesis Testing?

text: Central Limit Theorem says, "given certain conditions, the arithmetic mean of a sufficiently large number of iterates of independent random variables, each with a well-defined expected value and well-defined variance, will be approximately normally distributed". So this can be used in Hypothesis Testing, that even a distribution $M(\mu, \sigma)$ is not a normal distribution, if $x \sim M(\mu, \sigma)$, the mean of $n$ samples, $\bar x$, shall follow normal distribution $N(\mu, \sigma / \sqrt{n})$. The parameter $\mu$ has a way of hypothesis testing. However, how about the Variance? For normal distribution $N(\mu, \sigma)$, the $\sigma$ of $n$ samples follows $\chi^2(n)$ distribution, based on this one can set up the hypothesis testing. If the distribution $M(\mu, \sigma)$ is not normal, is there any similar way to set up a hypothesis testing for the Var (or standard deviation)? What I can think of now is find the 5% and 95% percentile from the distribution, and change this to a binomial test. But this sounds a bit loose?


title: Continuous dependent variable between 0 and 1 fitted with sigmoidal function

text: I've tried to find an answer to this to no avail. It is similar to this question, but not quite identical. I have a continuous dependent variable that occurs between zero and one (with some data points landing precisely on zero and one). I want to fit a $S$-shaped regression. My data is on coral reefs and the $Y$ variable is percentage dead per m$^2$. I don't have frequency counts in terms of successes and failures, just percentage mortality by cover. Is logit a terrible idea here? Since I only have two sig-figs would turning 80% mortality into 80 out of 100 observations be a bad idea? Is there some other sigmoidal function I might explore?


title: Autocovariance of an ARMA(2,1) process - derivation of analytical model for $\gamma( k)$

text: I need to derive analytic expressions for the autocovariance function $\gamma\left(k\right)$ of an ARMA(2,1) process denoted by: $y_t=\phi_1y_{t-1}+\phi_2y_{t-2}+\theta_1\epsilon_{t-1}+\epsilon_t$ So, I know that: $\gamma\left(k\right) = \mathrm{E}\left[y_t,y_{t-k}\right]$ so I can write: $\gamma\left(k\right) = \phi_1 \mathrm{E}\left[y_{t-1}y_{t-k}\right]+\phi_2 \mathrm{E}\left[y_{t-2}y_{t-k}\right]+\theta_1 \mathrm{E}\left[\epsilon_{t-1}y_{t-k}\right]+\mathrm{E}\left[\epsilon_{t}y_{t-k}\right]$ then, to derive the analytic version of the autocovariance function, I need to substitute values of $k$ - 0, 1, 2 ... until I get a recursion that is valid for all $k$ greater than some integer. Therefore, I substitute $k=0$ and work this through to get: $$ \gamma \left(0\right) = \mathrm{E}\left[y_t,y_t\right] = \phi_1 \mathrm{E}\left[y_{t-1}y_t\right] + \phi_2 \mathrm{E}\left[y_{t-2}y_t\right]+\theta_1 \mathrm{E}\left[\epsilon_{t-1}y_t\right]+\mathrm{E}\left[\epsilon_ty_t\right]\\ $$ now I can simplify the first two of these terms, and then substitute for $y_t$ as before: $$ \gamma\left(0\right) = \phi_1 \gamma\left(1\right) + \phi_2 \gamma\left(2\right)\\ + \theta_1 \mathrm{E}\left[\epsilon_{t-1} \left(\phi_1 y_{t-1} +\phi_2 y_{t-2} +\theta_1 \epsilon_{t-1} + \epsilon_t \right)\right]\\ + \mathrm{E}\left[\epsilon_t \left(\phi_1 y_{t-1} +\phi_2 y_{t-2} +\theta_1 \epsilon_{t-1} + \epsilon_t \right)\right] $$ then I multiply out the eight terms, which are: $$ +\theta_1\phi_1\mathrm{E}\left[\epsilon_{t-1}y_{t-1}\right]\\ +\theta_1\phi_2\mathrm{E}\left[\epsilon_{t-1}y_{t-2}\right]\\ +\theta_1^2\mathrm{E}\left[\left(\epsilon_{t-1}\right)^2\right]=\theta_1^2\sigma_{\epsilon}^2\\ +\theta_1\mathrm{E}\left[\epsilon_{t-1}\epsilon_{t}\right]=\theta_1\mathrm{E}\left[\epsilon_{t-1}\right]\mathrm{E}\left[\epsilon_{t}\right]=0\\ +\phi_1\mathrm{E}\left[\epsilon_{t}y_{t-1}\right]\\ +\phi_2\mathrm{E}\left[\epsilon_{t}y_{t-2}\right]\\ +\theta_1\mathrm{E}\left[\epsilon_t\epsilon_{t-1}\right] = \theta_1\mathrm{E}\left[\epsilon_{t}\right]\mathrm{E}\left[\epsilon_{t-1}\right]=0 \\ +\mathrm{E}\left[\left(\epsilon_t\right)^2\right] = \sigma_{\epsilon}^2 $$ So, I am left needing to resolve the four remaining terms. I want to use the same logic for lines 1, 2, 5 and 6 as I used on lines 4 and 7 - for example for line 1: $\theta_1\phi_1\mathrm{E}\left[\epsilon_{t-1}y_{t-1}\right] = \theta_1\phi_1\mathrm{E}\left[\epsilon_{t-1}\right]\mathrm{E}\left[y_{t-1}\right] = 0$ because $\mathrm{E}\left[\epsilon_{t-1}\right]=0$. Similarly for lines 2, 5 and 6. But I have a model solution that suggests the expression for $\gamma\left(0\right)$ simplifies to: $\gamma\left(0\right) = \phi_1\gamma\left(1\right)+\phi_2\gamma\left(2\right) +\theta_1\left(\phi_1+\theta_1\right)\sigma_{\epsilon}^2+\sigma_{\epsilon}^2$ This suggests my simplification as described above would miss the term with the coefficient $\phi_1$ - which under my logic should be 0. Is my logic at fault, or is the model solution I found incorrect? The worked solution also suggest that "analogously" $\gamma\left(1\right)$ can be found as: $\gamma\left(1\right) = \phi_1\gamma\left(0\right)+\phi_2\gamma\left(1\right) + \theta_1\sigma_{\epsilon}^2$ and for $k>1$: $\gamma\left(k\right) = \phi_1\gamma\left(k-1\right)+\phi_2\left(k-2\right)$ I hope the question is clear. Any assistance will be much appreciated. Thank you in advance. This is a question related to my research, and is not in preparation for any exam or coursework.


title: How to predict when the next failure of a machine will be

text: I need to predict when the next failure of a machine will be. I have a long history of its past failures, repairs, and preventive maintenance. By common sense, I can guess there must be some positive or negative influence of the past incidents to the next failure. But I cannot find the right way to shape up a model counting all these past transactions.


title: Machine learning classifiers big-O or complexity

text: To evaluate the performance a new classifier algorithm, I'm trying to compare the accuracy and the complexity (big-O in training and classifying). From Machine Learning: a review I get a complete supervised classifiers list, also a accuracy table between the algorithms, and 44 test problems from UCI data repositoy. However, I can't find a review, paper or web-site with the big-O for common classifiers like:   * C4.5   * RIPPER (I think this might not be possible, but who knows)   * ANN with Back Propagation   * Naive Bayesian   * K-NN   * SVM If anyone has any expression for these classifiers, it will be very useful, thank you.


title: p-value zero in hypothesis testing for survival curves

text: I have done survival analysis. I used Kaplan-Meir to do the survival analysis. Description of data: My data set is large and data table has close 120,000 records of survival information belong to 6 groups. Sample:                  user_id   time_in_days   event total_likes total_para_length group     1:       2          4657     1       38867        431117212   AA     2:       2          3056     1       31392        948984460   BB     3:       2            49     1          15            67770   CC     4:       3          4181     1       15778        379211806   BB     5:       3            17     1           3            19032   CC     6:       3          2885     1       12001        106259666   EE      After fitting the survival curves and plotting it, I see they are similar but yet at any given point in time their survival proportions don't seem to look like identical. Here is the plot: ![Survival Curves](http://i.stack.imgur.com/9wLYH.png) I ran a hypothesis test where my H0: There is not difference between the survival curves and here is the results that I got.               > survdiff(formula= Surv(time, event) ~ group, rh=0)     Call:     survdiff(formula = Surv(time, event) ~ group, rho = 0)                       N Observed Expected (O-E)^2/E (O-E)^2/V     group=FF 28310    27993    28632      14.3      19.0     group=AA 64732    63984    67853     220.6     460.1     group=BB 19017    18690    16839     203.4     245.6     group=CC  9687     9536     8699      80.6      91.0     group=DD 13438    13187    11891     141.3     164.2     group=EE  3910     3847     3324      82.4      89.7           Chisq= 788  on 5 degrees of freedom, p= 0        I am little confuse by trying to figure out what it means, specially since I got `p-value=0`. I am fairly new to survival analysis so after reading and digging through I realized that this is a non-parametric as I understand which means that it doesn't make any assumptions of the underline distributions of the time. After reading about cox-proportional hazard function and going over c-cran pdf I performed a cox regression test and here is what I got from that:               > cox_model <- coxph(Surv(time, event) ~ X)     > summary(cox_model)     Call:     coxph(formula = Surv(time, event) ~ X)            n= 139094, number of events= 137237                    coef  exp(coef)   se(coef)       z Pr(>|z|)         X1 -7.655e-05  9.999e-01  1.504e-06 -50.897   <2e-16 ***     X2 -1.649e-10  1.000e+00  5.715e-11  -2.886   0.0039 **      ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1             exp(coef) exp(-coef) lower .95 upper .95     X1    0.9999          1    0.9999    0.9999     X2    1.0000          1    1.0000    1.0000          Concordance= 0.847  (se = 0.001 )     Rsquare= 0.111   (max possible= 1 )     Likelihood ratio test= 16307  on 2 df,   p=0     Wald test            = 7379  on 2 df,   p=0     Score (logrank) test = 4628  on 2 df,   p=0      My big X is generated by doing rbind on total_like and total_para_length. Looking at Rsquare and P-Values I am not sure what really is going on here. If I can't throw away the Null-Hypothesis I should give a larger p-value.


title: Running a Pearson's correlation calculation on binary survey data

text: I have the following raw survey data output from a survey app: > ![enter image description here](http://i.stack.imgur.com/QvBv0.png) I was asked to run a Pearson's R on the results (1s and 0s). The format that is returned is as follows: > ![enter image description here](http://i.stack.imgur.com/hYzwO.png) I was further asked to get the results in the following format: > ![enter image description here](http://i.stack.imgur.com/mcmUt.png) Which is basically the correlation between the gender the other questions. I'm using SPSS to run these. First question: Does it make sense to run a Pearson's R on such data? Second Question: How do I achieve the second format of results? Does it make sense to sum then average out the individual correlation values for each question? Thanks! P.S. All figures and questions are purely made-up.


title: Regression and transforming variable with square root

text: I have a variable that can range from -100 to +100. When the number is -100, the output from this transformation should be something like 10. For +100 it should be around -10, and for 0 it should be 0. However, the points in between may be skewed so that the value for 1 is not necessarily .1, it may be 2 or .05. What formula or transform should I use when putting this variable into a linear model for R? I feel like there's a sqrt or log involved but negative values are no good. There are many options, what's the best? I'm looking for something like: $$y = B_1(x_1) + B_2(x_2) + B_3(transform(x_3))\\$$ I'd like the B3 to be able to tell me the skewness of the variable.


title: Normalizing a 2D-histogram and getting the marginals

text: I have a 2D-histogram for two vectors, _s1_ and _s2_ , generated using the _hist3_ function in Matlab:               [hist2D, binC] = hist3([s1' s2']);      I am normalizing it by making its total volume equal to unity, in the following manner:               L = binC{1}(2) - binC{1}(1);     B = binC{2}(2) - binC{2}(1);     totalVolume = sum(sum(hist2D.*L*B));     prob2D = hist2D/totalVolume;      **Question:** Is this correct way to normalize a 2D-histogram? I have also normalized the 1D-histograms for _s1_ and _s2_ , as shown below.               [hist1, binCentres1] = hist(s1);     binWidth1 = binCentres1(2) - binCentres1(1);     prob1 = hist1 / (sum(hist1) * binWidth1);     %same for s2      **Question:** How do I get the marginal (1D-)histograms from the normalized 2D-histogram? I have tried doing this in the following manner:               prob1M = sum(prob2D, 2); %extract marginal for s1     prob2M = sum(prob2D, 1); %extract marginal for s2      If I were doing this correctly, I'd expect _prob1_ to be equal to _prob1M_. I _seem_ to be on the right track because the bar graphs (below) look similar, but scaled on the vertical axis. Maybe I'm doing the normalizations wrong? ![Bar Graphs fro prob1 and prob1M](http://i.stack.imgur.com/3yMps.png) I tried normalizing _prob1M_ as well, using:               prob1M = prob1M / (sum(prob1M) * binWidth1);      **Question:** Is normalization still necessary if you get the marginals from a normalized 2D-histogram? Why / Why not? After normalization, the area of _prob1_ and _prob1M_ is equal to 1, where _area_ is calculated as:               area = sum(binWidth1.*prob1)     %same for prob1M      However, _prob1_ and _prob1M_ (and _prob2_ and _prob2M_ ) are still slightly different:               prob1 = 0.9412 0.4412 0.3235 0.3235 0.2941 0.3235 0.3235 0.4118 0.4706 1.1471     prob1M = 0.9706 0.4706 0.3824 0.3235 0.3235 0.3235 0.3530 0.3824 0.4412 1.0295      Thank you for your suggestions.


title: How can I convert annual standard deviation to a longer period?

text: Quicken provides annual standard deviation of returns for a given portfolio using analysis done by the Newport Group. I'd like to convert this to a longer term number--say 10, 20, or 30 years. According to this article, you can use the equation: ![equation for converting standard deviation](http://i.stack.imgur.com/dR1vq.gif) (Substituting T-years and annually for monthly and daily respectively.) Is this approximation accurate for such long T?


title: Derivation of the pdf of the Beta distribution

text: $$ f(x)= \frac{x^{(\alpha−1)} * (1−x)^{(\beta−1)}}{\mathcal B(\alpha,\beta)} $$ Questions:   1. How can we derive this famous pdf?    2. What is the intuitive meaning of the Beta distribution? (Please don't give the baseball example).    3. Please give some numerical problems with solutions where it is not explicitly mentioned that we have to use the Beta distribution (from our intuition we have to decide which distribution to use).


title: What statistical test should be used to accomplish the objectives of the experiment

text: I'd like to have your feedback to analyze my data from the following experiment. **Experiment Information**   1. A total of 30 (15 from one city and other 15 from another city) slaughtered animals were sampled.   2. Three tissues (Muscle, Kidney and Liver) were collected from each sampled slaughtered animal.   3. Each sample was analyzed for the presence of three drugs and each drug amount was measured.  **Objectives** The objectives of the experiment are to compare the level of drugs between   1. Types of drugs   2. Cities   3. Tissues.  Thanks in advance for your time and help.


title: Incorrect reasoning in obtaining average family income

text: This is a relatively easy question but I am new to statistics so I hope I could get a better idea. Me and my friend were discussing the following: My friend said that if someone wants to find the average monthly family income of a country and goes about finding this by saying: if the average family size is 4, I just multiply this by the average monthly income per person. However, I think think reasoning is incorrect but I can't really explain why.


title: Why law of large numbers does not apply in the case of Apple share price?

text: Here is the article in NY times called "Apple confronts the law of large numbers". It tries to explain Apple share price rise using law of large numbers. What statistical (or mathematical) errors does this article make?


title: How to conduct a factor analysis on questionnaire data based on 30 items in three blocks?

text: I have a questionnaire with 30 questions divided into 3 blocks of 10 questions each. Each question is answered with a number 1 to 10 and then an average of them all is calculated. I need to know how much each block has contributed to the final score (and if possible - how much each question has contributed). How do I do this? I have access to R.


title: Problem with p-values from A/B testing

text: I'm pretty confused right now and the more I look at it, the more confused I get. Overthinking! I'm doing some A/B tests for my website. I'm using the significance calculator that's in this article: http://visualwebsiteoptimizer.com/split-testing- blog/ab-testing-significance-calculator-spreadsheet-in-excel/ It tells me if I've reached significance levels of 90%, 95%, or 99%, based on my views and conversions. However, the p-value differs for different results. What I mean is... Say I have the following data:                         views   signups  z-score  p-score     control   13329   217     variant   13328   245      -1.3149  0.09426      The calculator tells me I have significance at 90%. The p-value here is .09426, putting me at 9.4%. The only way I see to make that actually match with over 90% is by subtracting it from 100. So 100 - 9.4 = 90.6. 90.6 > 90 = significance of 90% But then I have this situation:                         views   signups  z-score  p-score     control   12643   237           variant   12643   205      1.5356   0.93769      The calculator again tells me that I have significance of 90%. However here the p-value is .93769, putting me at 93.7%. 93.7 > 90 = significance of 90% But why is the p-value .09 for one and .93 for the other, when they both reach significance of 90%? For the first, I have to subtract the 9.4% from 100% to get the 90.6% significant rate. The other is .93, already 93%. No subtraction required. I don't understand. :( Plain English preferred if possible.


title: Inferring hidden categorical variables from multimodal response data

text: Lets say I have a response variable that is bimodally distributed. Below is a simple example with 1 predictor variable: ![](http://i.stack.imgur.com/bgXmI.png) This can be dealt with by adding a categorical variable that splits the data into two groups. For example, in the toy data above you could imagine y=child’s height (full grown) and x=parent’s height. The data above the best fit line would be males and the data below the best fit line would be females. Let’s say, for arguments sake, that we didn’t originally record the gender of each child. **First (small) question:** I suppose one could use a clustering algorithm to recover the two groups in the data, and then create a linear regression model with the clustering output as an additional categorical predictor variable. Is this ever done? Are there problems with this logic? **Second (more important) question:** Let’s say we want to predict the height of a new child, but we don’t know his/her gender. Rather than giving the expected value (red line), I want to return two possible values (i.e. height if male, or height if female). What is the best way to do this? **Preliminary answer:** Use the clustering approach I mentioned above, and fit an expanded linear regression model with an additional categorical variable which nominally refers to gender. Then, given a new observation x, predict y in the expanded model for both possible cases. **Difficulties:** I am interested in how to solve this problem for much higher dimensional data sets, and specifically how to automate it. In general, it will not be known how many clusters there are in the data, or which variables to cluster the data along. The data I am interested in might also be highly non-linear. Consider the example below: ![](http://i.stack.imgur.com/TRAAH.png) I think this is an interesting question for dynamical systems that switch between multiple attractor states. For example a neuronal network with a noisy (and uncontrolled) input that can “bump” or switch the activity pattern between multiple (unknown) states. This might lead to different clusters in the response variable as in the above example. See: http://www.scholarpedia.org/article/Attractor_network Matlab code for generating figure 2:               N = 2e2;     x = rand(N,1);     y1 = nan(N,1); y2 = nan(N,1);     y1(x<0.4) = 1.2*x(x<0.4)+0.2+0.03*randn(sum(x<0.4),1);     y2(x<0.4) = x(x<0.4)+0.6+0.03*randn(sum(x<0.4),1);     y1(x>=0.4 & x<0.6) = 0.5*x(x>=0.4 & x<0.6)+0.01*randn(sum(x>=0.4 & x<0.6),1);     y2(x>=0.4 & x<0.6) = 0.6*x(x>=0.4 & x<0.6)+0.01*randn(sum(x>=0.4 & x<0.6),1);     y1(x>=0.6) = 0.5+0.02*randn(sum(x>=0.6),1);     y2(x>=0.6) = 0.5+0.02*randn(sum(x>=0.6),1);     y = [y1; y2];     x = repmat(x,2,1)+0.02*randn(2*N,1);     plot(x,y,'ob')     xlabel('x','FontSize',20)     ylabel('y','FontSize',20)


title: Variational Posterior Dirichlets in LDA

text: I am running the c code for LDA provided on David Blei's website. The code outputs several files. The output file final.gamma is supposed to include the "Variational Posterior Dirichlets". If I understand this correctly (which it is very likely I do not) it is supposed to be for every document the probability that that document was generated by each of the k topics. If my understanding is correct then the sum of the posteriors for each document should add up to 1. i.e. the sum over the k topics of the probability that that topic generated the document should be 1 (every article is assumed to be generated by some mixture of the topics). I ran LDA with 2 topics and 20 documents and got the following final.gamma file:               0.0163467211 439.0163467211     476.0163467211 0.0163467211     0.0163467211 141.0163467211     151.0163467211 0.0163467211     0.0163467211 526.0163467211     0.0163467211 78.0163467211     52.3257109467 314.7069824954     0.0163467211 356.0163467211     2448.0163467211 0.0163467211     0.0163467211 356.0163467211     0.0163467211 554.0163467211     0.0163467211 632.0163467211     373.0163467211 0.0163467211     0.0163467211 533.0163467211     726.0163467211 0.0163467211     0.0163467211 501.0163467211     38.8105963538 347.2220970883     470.0163467211 0.0163467211     0.0163467211 437.0163467211     723.0163467211 0.0163467211      As you can see for each document there are two numbers representing the two topics that generated that article. These numbers do not sum up to 1, what do they mean? For this output I ran LDA on twenty documents the first ten relating to baseball and the second ten relating to football. I would expect LDA to assign the first ten documents similar posteriors and the second ten documents similar posteriors but that is not the case. The words that the topics generated are clearly split on the football/baseball axis i.e. topic 0 generates only baseball words and topic 1 only football words so I would expect the posteriors to work. Is this a bug in the LDA implementation or am I misunderstanding the meaning of the posteriors?


title: Error in scipy Wilcoxon signed-rank test

text: I was digging in the scipy code for this test Wilcoxon signed-rank test (stats.Wilcoxon) and I found that in scipy they compute the sum of the ranks for the differences that are positive and separately for the ones that are negative. Then they picked the smaller one and use that as W. That is substantially different from the test explanation in Wikipedia, or other sites where `W = sum(all_differences)`. **Is this approach valid?**


title: Notation for the sample total variation

text: Given a sample $x=(x_1, \ldots, x_n)$, what is your favorite notation for the sum of squares $\sum (x_i -\bar x)^2$ ? I like to keep $s^2_x$ for the sample variance, and I wonder whether $ss^2_x$ is a nice notation for the sum of squares. Now if there are two paired samples $x$ and $y$, how to denote $\sum(x_i-\bar x)(y_i - \bar y)$ (and by the way how do we call it) ? If we use $ss^2_x$ for the total variation, should we use $ss_{xy}$ or $ss^2_{xy}$ ? Actually I'm currently writing a doc about Deming regression and here is an example of an expression using the above notations : $$\frac{ss^2_x ss^2_y - ss^4_{xy}}{ss^4_{xy}}. $$ Strange for the eyes, no ? Moreover the first time I typed this expression I did a mistake by confusing $ss_{xy}$ and $ss^2_{xy}$. What do you think about $v^2_x$ ?


title: Is Mahalanobis distance dependent from the vector dimensions?

text: I would like to know if Mahalanobis distance is a normalized measure and the feature dimensions (or vector dimensions) do not affect the measure? In other words if I have a vector $ v \in \Re^n $ with an associated Mahalanobis distance $ D_v $ and a vector $ x \in \Re^m $ where $ m >> n$ and with a distance $ D_x$ associated, could I directly compare $ D_v $ with $ D_x$ value or there is a dimension effect? If is there a dimension effect how it is possible to compensate for it? Thank you.


title: MLE of multinomial distribution with missing values

text: I am trying to solve a problem and the results I get seem counter-intuitive. We randomly throw $n$ balls into an area partitioned into 3 bins $b_1,b_2,b_3$. The size of each bins is proportional to the probability the ball will fall in it. Let's call these probabilities $p_1,p_2,p_3$. This can be described by a multinomial distribution. Now let's say I throw 12 balls, and I know how many landed in each bin ($x_1=3,x_2=6,x_3=3$). I would like to estimate the size of the bin from the observations. For this I use Maximum Likelihood. It can be shown that the MLE will be $p_1=3/12,p_2=6/12,p_3=3/12$. This is pretty intuitive. It turns out that the actual likelihood at this point is: $L(p_1=0.25,p_2=0.5,p_3=0.25|x_1=3,x_2=6,x_3=3)=$ $=\frac{12!}{3!6!3!}0.25^30.6^60.25^6=0.07050$ Now, let's assume I knew in advance that $p_1=p_3$. How would that change my result? It would not - I would still get the same parameter values $p_1=0.25,p_2=0.5,p_3=0.25$. The twist comes now: let's assume I cannot observe balls that landed in $b_3$. If I know that 12 balls were thrown I am fine, since I can calculate $b_3=n-b_1-b_2=12-3-6=3$. but what happens if I don't know $n$? I figure that in this case, I would need to estimate $x_3$ (or equivalently $n$) as well. However, if I use MLE, the results start looking weird. Intuitively, I would expect that if I observe $x_1=3,x_2=6$ and I know that $p_1=p_3$, then the MLE will probably be $p_1=0.25,p_2=0.5,p_3=0.25,x_3=3$. However, it is clearly not the maximum, since for example: $L(p_1=0.24,p_2=0.52,p_3=0.24|x_1=3,x_2=6,x_3=2)=$ $=\frac{11!}{3!6!2!}0.24^30.52^60.24^2=0.07273$ So from this it seems that $x_1=3,x_2=6,x_3=2$ is more likely than $x_1=3,x_2=6,x_3=3$ even if I know that $p_1=p_3$, which seems very counter- intuitive. **My questions are whether my logic is sound, whether my intuition is misleading me and whether this is the correct way to estimate the parameters and missing data.** **EDIT:** For future reference, I found this highly relevant paper, which addresses this exact problem. It also discusses the slight skew mentioned in whuber's answer.


title: Help me plan/design a survey

text: We are doing a survey regarding what members of a club think about the club. Say we have roughly 1000 members, who will be divided into five groups of 200 in each. Then each group will be given a survey at different times of the year (one group March, one April, and so on). Besides the fact that people may think differently about the club over the year, and that we make sure that the groups are fairly similar. Is there anything else statistically (or otherwise) that we need to think of?


title: Mixed, repeated measure model specification and results interpretation using LMER in R

text: I am working with data from a computer task which has 288 total trials, each of which can be categorically classified according to **Trial Type** , **Number of Stimuli** , and **Probe Location**. Because I want to also examine a continuous variable, the total Cartesian **Distance** between stimuli per trial (divided by number of stimuli to control for varying numbers), I have opted to use a mixed linear model with repeated measures. In addition to each of these task variables, I am also interested in whether folks in various diagnostic groups perform differently on the task, as well as whether or not there is a **Dx** interaction with any of the above variables. Thus (if I'm not mistaken), I have the following effects in my model: **Trial Type** , a fixed effect **Number of Stimuli** , a fixed effect **Probe Location** , a fixed effect **Dist** (ance), a fixed effect **Dx** , a fixed effect **Dx*Trial Type** , a fixed effect **Dx*Number of Stimuli** , a fixed effect **Dx*Probe Location** , a fixed effect **Dx*Dist** , a fixed effect **Trial** , a random effect, nested within **SubID** , a random effect Based on my examination of documentation, it seems that the nesting of random effects does not seem to be important to lme4, and so I specify my model as follows: `tab.lmer <\- lmer(Correct ~ Dx+No_of_Stim+Trial_Type+Probe_Loc+Dist+Dx*No_of_Stim+Dx*Trial_Type+Dx*Probe_Loc+Dx*Dist+(1|Trial)+(1|SubID),data=bigdf)` This would be my first question: **1) Is the above model specification correct?** Assuming so, I am a bit troubled by my results, but as I read and recall my instruction on such models, I am wondering if interpretation of particular coefficients is bad practice in this case:               Linear mixed model fit by REML ['merModLmerTest']     Formula: Correct ~ Dx + No_of_Stim + Trial_Type + Probe_Loc + Dist + Dx *           No_of_Stim + Dx * Trial_Type + Dx * Probe_Loc + Dx * Dist +           (1 | Trial) + (1 | SubID)        Data: bigdf          REML criterion at convergence: 13600.4          Scaled residuals:           Min       1Q   Median       3Q      Max      -2.89810 -0.03306  0.27004  0.55363  2.81656           Random effects:      Groups   Name        Variance Std.Dev.      Trial    (Intercept) 0.013256 0.11513       SubID    (Intercept) 0.006299 0.07937       Residual             0.131522 0.36266      Number of obs: 15840, groups:  Trial, 288; SubID, 55          Fixed effects:                              Estimate Std. Error         df t value Pr(>|t|)         (Intercept)             4.196e-01  4.229e-02  4.570e+02   9.922  < 2e-16 ***     DxPROBAND               8.662e-02  4.330e-02  2.920e+02   2.000  0.04640 *       DxRELATIVE              9.917e-02  4.009e-02  2.920e+02   2.474  0.01394 *       No_of_Stim3            -9.281e-02  1.999e-02  4.520e+02  -4.642 4.53e-06 ***     Trial_Type1             3.656e-02  2.020e-02  4.520e+02   1.810  0.07097 .       Probe_Loc1              3.502e-01  2.266e-02  4.520e+02  15.456  < 2e-16 ***     Probe_Loc2              3.535e-01  3.110e-02  4.520e+02  11.369  < 2e-16 ***     Dist                    1.817e-01  2.794e-02  4.520e+02   6.505 2.06e-10 ***     DxPROBAND:No_of_Stim3  -1.744e-02  1.759e-02  1.548e+04  -0.992  0.32144         DxRELATIVE:No_of_Stim3 -2.886e-02  1.628e-02  1.548e+04  -1.773  0.07628 .       DxPROBAND:Trial_Type1  -9.250e-03  1.777e-02  1.548e+04  -0.521  0.60267         DxRELATIVE:Trial_Type1  1.336e-02  1.645e-02  1.548e+04   0.812  0.41682         DxPROBAND:Probe_Loc1   -8.696e-02  1.993e-02  1.548e+04  -4.363 1.29e-05 ***     DxRELATIVE:Probe_Loc1  -4.287e-02  1.845e-02  1.548e+04  -2.323  0.02018 *       DxPROBAND:Probe_Loc2   -1.389e-01  2.735e-02  1.548e+04  -5.079 3.83e-07 ***     DxRELATIVE:Probe_Loc2  -8.036e-02  2.532e-02  1.548e+04  -3.173  0.00151 **      DxPROBAND:Dist         -3.920e-02  2.457e-02  1.548e+04  -1.595  0.11066         DxRELATIVE:Dist        -1.485e-02  2.275e-02  1.548e+04  -0.653  0.51390         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1      In general, these results make sense to me. The troubling portion, however, comes in the positive, significant (yes, I am using LmerTest) p-value for DxProband, particularly in light of the fact that in terms of performance means, Probands are performing worse than Controls. So, this mismatch concerns me. Examining the corresponding ANOVA:               > anova(tab.lmer)     Analysis of Variance Table of type 3  with  Satterthwaite      approximation for degrees of freedom                    Sum Sq Mean Sq NumDF   DenDF F.value    Pr(>F)         Dx             0.8615  0.4308     2   159.0   1.412   0.24662         No_of_Stim     0.6984  0.6984     1   283.5  37.043 3.741e-09 ***     Trial_Type     8.3413  8.3413     1   283.5   4.456   0.03565 *       Probe_Loc     25.7223 12.8612     2   283.5 116.405 < 2.2e-16 ***     Dist           5.8596  5.8596     1   283.5  43.399 2.166e-10 ***     Dx:No_of_Stim  1.4103  0.7051     2 15483.7   1.590   0.20395         Dx:Trial_Type  2.0323  1.0162     2 15483.7   0.841   0.43128         Dx:Probe_Loc   3.5740  0.8935     4 15483.7   7.299 7.224e-06 ***     Dx:Dist        0.3360  0.1680     2 15483.7   1.277   0.27885         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1      ...the results seem to more or less line up with the regression, with the exception of the **Dx** variable. So, my second question is **2) Can anyone clarify what is going on with the Dx variable? Is interpreting individual coefficients from the regression model bad practice in this case?** Finally, as a basic (and somewhat embarrassing) afterthought, **3) If I attempt to reduce the model, I should favor the model with the _lower_ REML, yes?** In summation, **1) Is the above model specification correct?** **2) Can anyone clarify what is going on with the Dx variable? Is interpreting individual coefficients from the regression model bad practice in this case?** **3) If I attempt to reduce the model, I should favor the model with the _lower_ REML, yes?** **ADDENDUM** By request, I'll describe the task and data a little further. The data come from a computer task in which participants are presented a number of stimuli, either two or three, in various locations about the screen. These stimuli can either be "targets" or "distractors". After these stimuli, a probe stimulus is presented; if it appears in a position where a previous target has appeared, participants should respond "yes"; if it appears in the position of a previous distractor or elsewhere, the correct answer is "no." There are 288 trials of this nature; some have two stimuli and some have three, and some lack distractors entirely. The variables in my model, then, can be elaborated as follows: **Number of Stimuli:** 2 or 3 (2 levels) **Trial Type:** No Distractor (0) or Distractor (1) (2 levels) **Probe Location:** Probe at Target (1), Probe at Distractor (2), or Probe Elsewhere (0) (3 levels) **Distance:** Total Cartesian distance between stimuli, divided by number of stimuli per trial (Continuous) **Dx:** Participant's clinical categorization **Sub ID:** Unique subject identifier (random effect) **Trial:** Trial number (1:288) (random effect) **Correct:** Response classification, either incorrect (0) or correct (1) per trial Note that the task design makes it inherently imbalanced, as trials without distractors cannot have Probe Location "Probe at Distractor"; this makes R mad when I try to run RM ANOVAs, and it is another reason I opted for a regression. Below is a sample of my data (with SubID altered, just in case anyone might get mad):                    SubID      Dx Correct No_of_Stim Trial_Type Probe_Loc      Dist Trial     1 99999999 PROBAND       1          3          0         1 0.9217487     1     2 99999999 PROBAND       0          3          0         0 1.2808184     2     3 99999999 PROBAND       1          3          0         0 1.0645292     3     4 99999999 PROBAND       1          3          1         2 0.7838786     4     5 99999999 PROBAND       0          3          0         0 1.0968788     5     6 99999999 PROBAND       1          3          1         1 1.3076598     6      Hopefully, with the above variable descriptions these data should be self- explanatory. Any assistance that people can provide in this matter is very much appreciated. Sincerely, peteralynn


title: R CCH CI Univariable Analysis

text: Using the survival package in R I am trying to calculate 95% confidence interval (CI) of the proportional hazards model coefficient. Normally this is pretty simple to do with the survival package.               > confint(cch(Surv(time, event) ~ a + b + c, ...))                   2.5 %    97.5 %     a   -0.4216482 0.1225872     b    0.2624281 1.3339170     c    0.3344311 1.2224415      However when I try to do this in a univariable setting, the confint function fails to return anything, except a matrix with no rows.               > confint(cch(Surv(time, event) ~ a, ...))                   2.5 %    97.5 %      Is this a bug in R or am I missing something else? If this is a bug how can I calculate these CIs 'by hand'?


title: Normalize data in unnormalized data after normalization

text: i have some data in numeric. I want to do some classifications method. So i decided to check the normalization of the data. I have normalize my data using weka tools. And i think weka had normalized my data correctly. And then i check the normality of my data. But it still unnormalized data.. So, what should i do to this condition? Should i do normalization again until the normality is reached? Solution will be appreciated.. :)


title: Experimental design & questions on use of generalized linear models

text: I have an ecological experiment for which I need to analyze bird count data. Here is the set up: 2 treatments (open/control), 3 regions. Not quite a full 3x2 factorial because in 2 regions there are 3 plots (250m x 250m) for both open & control (6 plots). For “reference” region there are just 3 plots (more similar to open plots in the upstream region). Birds were counted weekly (can be considered independent) in each plot over 4 years. Bird type & densities are highly dependent on water depth (have depth for each survey date) – would like to use as covariate (note: the response isn’t linear, it’s unimodal, and it differs by species). The goal is to examine treatment vs control differences in numbers of types of birds, & examine differences/interactions of region & year (as open plots mature & vegetation grows in). The closed plots contain many zeros (I may need zero-inflated model). The literature is pointing me towards negative binomial or zero-inflated Poisson GLM. If I choose to do this, though, I have following questions:   * Can I nest plots within region for this type of model? I wanted unit of replication to be the week for each treatment, without sacrificing replication for each treatment x region plot.    * Also, any suggestions on what to do with a covariate like depth? It could really dominate the terms of a GLM. I'm still learning R, but assuming based on this forum & speaking with others that I will need to do this in R or SAS. Advice is welcome on the software front, including what R packages I'll need to load.


title: Estimating parameter using curve-fitting and model comprised of uncorrelated product of two functions

text: Suppose that I have a time-series dataset represented by the function $f[t,\omega,q]$, where $t$ is the time as an independent variable, $\omega$ is the frequency as an independent variable and $q_k$ is an unknown that I want to determine for each discrete $t_k$ timestep. The $q_k$ is constant over all $\omega$ for a certain $t_k$. In other words, each timestep $t_k$ has a number of $f[t_k,\omega,q_k]$ corresponding to different $\omega$, but $q_k$ is constant for each $t_k$. The observations are modeled by the product of two functions: $f[t,\omega,q] = A[t,\omega]\exp(g[\omega,t,q])$ In the above, $A[t,\omega]$ is unknown, and changes with $t$ and $\omega$. The mathematical formula for $g[\omega,t,q]$ is known, but I do not know $q_k$, the parameter to be determined at every timestep $t_k$. I can make the assumption that $A[t,\omega]$ is uncorrelated with $g[\omega,t,q]$, and I will not be able to measure $A[t,\omega]$, so I cannot make very many assumptions about it. How would I set up a curve-fitting problem to estimate $q_k$ at each timestep $t_k$? Is there a good reference on this type of problem? What will affect the accuracy of the $q_k$ estimation, even if $A[t,\omega]$ is uncorrelated with $g[\omega,t,q]$? Can Bayesian methods be useful in such a problem? In addition, $A[t,\omega]$ and $g[\omega,t,q]$ and $f[t,\omega,q]$ can be complex numbers. How would I deal with this in curve-fitting?


title: How many data points for my new product do I need to tell if my new product performs differently to the old?

text: We've been using the old formula to produce a product for years. We take a sample every quarter or so and run it through 8 different quality tests (every quarter because it takes that long for the tests to complete). Currently 53 samples of the old formulation have been put through these tests. One of the components of the old formula became unavailable so we're using a new formula that replaces the unavailable component with what we believe is a suitable substitute. Again, every quarter or so, we're going to take a sample and run it through the same tests as we did for the old formula. How many samples of the new formula are we going to have to test to be confident we can tell whether the new formula is "significantly different" or "not significantly different" to the old formula? (Whether it is significantly different or not is beside the point, I just want to know when we can be 95% confident that we are reaching the right conclusion - different or not different - rather than having to wait for 53 samples and 13 years to find out). If the answer to this question is "it depends", could someone please explain what it depends on? Much gratitude.


title: Conditional density plots with two continuous variables in R

text: I am very new to `R`. I have to create conditional density plots with two continuous variables. I have found a nice example with which to start, however this sample has `y` as factor variable. I have tried to change the command to two continuous variables, however unsuccessfully so far. Can someone help me to re-write this example command for the situation where both variables are continuous and `y` variable, e.g., takes following values: `c(55, 58, 60, 66, 66, 67, 65, 64, 70, 69, 80)`? Here is the example I’m following: ## NASA space shuttle o-ring failures               fail <- factor(c(2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 2),                    levels = 1:2, labels = c("no", "yes"))     temperature <- c(53, 57, 58, 63, 66, 67, 67, 68, 69, 70, 81)      ## CD plot               cdplot(fail ~ temperature)     cdplot(fail ~ temperature, bw = 2)     cdplot(fail ~ temperature, bw = "SJ")      ## compare with spinogram               (spineplot(fail ~ temperature, breaks = 3))      ## highlighting for failures               cdplot(fail ~ temperature, ylevels = 2:1)      ## scatter plot with conditional density               cdens <- cdplot(fail ~ temperature, plot = FALSE)     plot(I(as.numeric(fail) - 1) ~ jitter(temperature, factor = 2),          xlab = "Temperature", ylab = "Conditional failure probability")          lines(53:81, 1 - cdens[[1]](53:81), col = 2)      ## NASA space shuttle o-ring failures               fail <- c(2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 2)     temperature <- c(53, 57, 58, 63, 66, 67, 67, 68, 69, 70, 81)      ## CD plot               cdplot(fail ~ temperature)     cdplot(fail ~ temperature, bw = 2)     cdplot(fail ~ temperature, bw = "SJ")      ## compare with spinogram               (spineplot(fail ~ temperature, breaks = 3))      ## highlighting for failures               cdplot(fail ~ temperature)      ## scatter plot with conditional density               cdens <- cdplot(fail ~ temperature, plot = FALSE)     plot(I(as.numeric(fail) - 1) ~ jitter(temperature),          xlab = "Temperature", ylab = "Conditional failure probability")     lines(53:81, 1 - cdens[[1]](53:81), col = 2)


title: How does Fisher LDA work?

text: Intuitively, how does Fisher LDA work? From this Linear discriminant analysis and Bayes rule I completely understood the Bayesian approach but I'm not able to relate it to the Fisher's one described there. What's the relation between the linear combination of predictor and posterior probability in Bayesian approach? Why do the coefficients' estimate have that form? I know that this might be a trivial questions but I'm studying this on my own and my mathematical background is weak.


title: Avoiding matrix inversions with the multivariate Gaussian distribution?

text: Let $(X, Y)$ be a Gaussian random vector with mean $\mu$, variance $\Sigma$ and both $X$ and $Y$ are vectors. No structure is assumed on $\Sigma$ (aside from being symmetric positive definite). Is there a way of calculating $(X - \mu_x)^T \Sigma_{xx}^{-1} (X - \mu_x)$ and getting the parameters of the conditional distribution of $Y|X$, $\mu_{Y|X} = \mu_y + \Sigma_{yx}\Sigma_{xx}^{-1}(X - \mu_x)$ and $\Sigma_{Y|X} = \Sigma_{yy} - \Sigma_{yx}\Sigma_{yy}^{-1}\Sigma_{xy}$ without having to invert $\Sigma$? Asking because the "golden rule of numerical linear algebra" suggests that matrix inversion can usually be avoided if you are smart enough, but I am not very smart. Or is this a situation where you really do need to do a matrix inversion? In retrospect this seems more likely to be the case since people often impose structure on $\Sigma$ for computational reasons.


title: What type of statistics test do I use for the following

text: I have three sets of means. For instance               mean.a = [3,2,7,3,8,2,5,7,3,2,5]     mean.b = [9,7,9,10,7,13,30,17,16,9,9]     mean.c = [23,24,25,23,22,22,21,9,25,34,12]      And for each mean, I have standard deviations.               sd.a = [.1,.18,.12,.24,.12,.16,.25,.1,.1,.23,.11]     sd.b = [.3,.13,.1,.24,.2,.1,.24,.11,.1,.23,.1]     sd.c = [.3,.1,.33,.21,.32,.35,.31,.3,.33,.42,.42]      I want to statistically show that `a` is less than `b` and `c`. How can I do this, and find a p value for this? Note that the values are paired. So, for instance, the 5th value of `a` should be less than the 5th value of `b` and 5th value of `c`. It was suggested to me to use a Wilcoxon signed rank test to show, firstly, that c > b then that b > a. I found the R function `wilcox.test(a, b, paired=T)`. However, I have means and SDs for both `a` and `b`, so how would I do this?


title: Filter out information in Kalman Filter

text: I have a disagreement with friend over if an observable that depend on lagged observable, should be called a state variable. I have the following system:               xt=u+z(1,t-1)+x(t-1)+d·[y(t-1)-x(t-1)]+e(x,t)     yt=u+z(2,t-1)+y(t-1)-d·[y(t-1)-x(t-1)]+e(y,t)          z1 and z2 are AR(1) processes     z(1,t)=p1z(1,t-1)+e(1,t)     z(2,t)=p2z(2,t-1)+e(2,t)      The observations are xt and yt I say the state variables are z1 and z2, but my friend say I am wrong in that the state variable must also include x(t-1) and y(t-1) Since I think that x(t-1) and y(t-1) are observed and known at time t, can I claim that they are known and thus can be filtered out And as far as the Kalman filter gain and the variance of estimate P , as well as their steady state values are concerned, I can derive K and P and their steady state values by pretending the system is the following, which will yield the exact same K and P:               xt=u+z(1,t-1)+x(t-1)+e(x,t)     yt=u+z(2,t-1)+y(t-1)+e(y,t)          z1 and z2 are AR(1) processes     z(1,t)=p1z(1,t-1)+e(1,t)     z(2,t)=p2z(2,t-1)+e(2,t)      Am I right that the 2nd system will yield the same K and P and their respective steady state value as the first system? But my friend say I am wrong because x(t-1) and y(t-1) are not known at time 0, they must be state variables, and thus xt and yt are both observables variables and state variables.


title: Does the standard deviation of the folds of LOO cross-validation have any practical meaning in comparison/evaluation of classifiers?

text: There are dozens of questions regarding LOO and variance. Most of the answers are purely theoretical or too general. I have also read many papers like this paper. Specifically: I have two not too stable classifiers/learning algorithms (suppose C4.5 and RIPPER for simplicity), the dataset is small, thus LOO is desirable. I need to compare the accuracy of models generated from both algorithms (it is not about model selection nor accurate absolute estimates of the prediction error nor how to setup experiments). **Is there any value calculating, e.g. the standard deviation in this case? ps. it is a single run, the the standard deviation is over the hits obtained by the models during LOO** _It seems to me that LOO's accuracy estimate and LOO's StDev are function of each other, almost redundant information. If you take from one, immediately you have given to the other._ * * * Another way to think about this is to realize that the "internal" variance of LOO is due to _scarce testing data_ and not due to _lack of accuracy of the model_. However, this brings another question, which is how much the scarceness of testing data affects the the "internal" variance of other CVs besides LOO.


title: Use of Deviance Residuals for Leave-One-Out Cross Validation

text: I am a newbie to stats and having some difficulties understanding how to use deviance residuals for leave-one-out cross validation for a logistic regression model. The problem that I am trying to solve is unique in that the structure of the data is not particularly suitable for applying one of the procedure in software such as SAS. The problem is similar to this one (http://support.sas.com/kb/24/992.html) but assume that the data is based on opinions from 13 individuals. To validate the model, I wanted to leave data out for one individual at a time and generate deviance residuals for each case. I now have the deviance residuals for the full case with data from 13 individuals, and additional 13 cases where data from one individual is removed. My questions are (1) how do I use these deviance residuals to perform the cross validation (leave-one-out) (2) what conclusions about my model can I be able to draw from the results? I would like to do this manually because of the nature of my data (calculated the residuals separately from each case). Also, I am also open to other suggestions for accessing the predictive capability of such model.


title: Metropolis Algorithm for a 3 variated function

text: I need to implement a program that generate a sample from a 3 variated $pdf$ using the Metropolis algorithm (and its variations). I was thinking to use a 3-variated normal distribution as my proposal distribution, but I am having some difficult and doubts doing it. For example, I've seen this simulation in R: http://blog.abhranil.net/2014/02/08/r-code-for-multivariate-random-walk- metropolis-hastings-sampling/ 1) He uses metropolis random walk with mean $X^t$ and a given covariance matrix. It seens he uses the fact that this distribution is symmetric. I know it is a normal, but is it always symmetric? 2) What would be better: Draw from a trivariated normal directly or draw each parameter from univariated normal? In the last case, we can't adjust covariance, so it looks less efficient and I don't even know it is valid. Thanks in advance!


title: Beginning anova in R: Phantom significance when testing for interaction terms

text: Suppose I create a dummy scenario as such:               > A <- rnorm(10000)      > B <- rnorm(10000)      > C <- rnorm(10000)      > Y <- A*B + rnorm(10000,sd=0.1)      Doing a simple ANOVA correctly identifies that none of the variables are significantly predictive of the outcome:               > anova(lm(Y~A+B+C))     Analysis of Variance Table          Response: Y                 Df  Sum Sq Mean Sq F value Pr(>F)     A            1     1.5 1.54411  1.4209 0.2333     B            1     0.3 0.28909  0.2660 0.6060     C            1     1.6 1.62425  1.4946 0.2215     Residuals 9996 10862.8 1.08672          But not let's say I decide to include the interaction terms:               > anova(lm(Y~A*B*C))     Analysis of Variance Table          Response: Y                Df  Sum Sq Mean Sq    F value    Pr(>F)         A            1     1.5     1.5 1.5281e+02 < 2.2e-16 ***     B            1     0.3     0.3 2.8610e+01  9.05e-08 ***     C            1     1.6     1.6 1.6074e+02 < 2.2e-16 ***     A:B          1 10761.8 10761.8 1.0650e+06 < 2.2e-16 ***     A:C          1     0.0     0.0 9.8700e-02    0.7534         B:C          1     0.0     0.0 1.5062e+00    0.2197         A:B:C        1     0.0     0.0 1.6790e-01    0.6820         Residuals 9992   101.0     0.0                              ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1      It has correctly identified the interaction between A and B as being the most significant, but now for some reason the individual terms A and B have also gained significance... and C which had nothing at all to do with creating the model is significant as well? Either I have not written the test correctly or I am completely misunderstanding how a Two-Way ANOVA with interaction terms works Using a simple linear model gives expected results:               > summary(lm(Y~A*B*C))          Call:     lm(formula = Y ~ A * B * C)          Residuals:          Min       1Q   Median       3Q      Max      -0.29566 -0.06667 -0.00092  0.06665  0.33620           Coefficients:                   Estimate Std. Error  t value Pr(>|t|)         (Intercept) -0.0003212  0.0009707   -0.331    0.741         A            0.0003483  0.0009613    0.362    0.717         B            0.0003184  0.0009619    0.331    0.741         C           -0.0003213  0.0009702   -0.331    0.741         A:B          1.0008711  0.0009370 1068.214   <2e-16 ***     A:C         -0.0014855  0.0009588   -1.549    0.121         B:C          0.0008860  0.0009561    0.927    0.354         A:B:C       -0.0002489  0.0009085   -0.274    0.784         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Residual standard error: 0.09705 on 9992 degrees of freedom     Multiple R-squared:  0.9913,    Adjusted R-squared:  0.9913      F-statistic: 1.634e+05 on 7 and 9992 DF,  p-value: < 2.2e-16


title: Equal number of training instances of each classification label?

text: I am using Naive Bayes to perform binary classification. In my training set, the two class labels occur with probability Pr(label A) = 0.95 and Pr(label B) = 0.05. Should I prune the training set so that there is an equal number of training instances of each label? Does the answer apply to any classifier, not just Naive Bayes?


title: What statistical technique would be most appropriate for comparing two instructional approaches?

text: I work in a school and would like to use a statistical technique for comparing the performance of students using two different instructional programs. A first grade teacher uses a new, specialized spelling program not used by any other teacher in the school. The following year, a second grade teacher would like to find out if students from the first grade teacher's class are better spellers than students from other classes. The second grade teacher has fairly small samples. There are 3-5 students in the classroom from the previous first grade teacher's class and 11-13 students from other first grade teachers who didn't use the specialized program. The teacher would like to evaluate spelling by the number of words correct on a well-recognized 100 word spelling inventory. Which statistical techniques would be appropriate for helping the teacher determine if the students who did or didn't receive the specialized spelling program are achieving at the same level? Due to the small sample sizes, would it be useful to use a Mann-Whitney U Test? Are the sample sizes too small to look at the relationship through a Point-Biserial Correlation Coefficient? Thank you very much for any assistance.