title: Random seed choice changes qualitative results of elastic net parameter tuning

text: Below is a function I wrote to try and tune the $\lambda$ and $\alpha$ elastic net GLM implemented with `cv.glmnet`. I've noticed that the qualitative outcome (in terms of the alpha that yields the lowest prediction error score) changes depending on the random seed, and I suspect this is mainly because the random seed determines how the folds are partitioned. In one of my datasets, I have only 598 samples and am fitting a Cox regression. That is where this effect is strongest (and, indeed, where the error bars about the prediction errors on the $\lambda$ search grid are largest). When I fit a binomial regression to a larger data set with about 30K observations, I get less jumping around. The number of folds doesn't seem to matter much. I'm thinking of just choosing $\alpha$ = 0.5 in such cases, because the absolute difference in the cross-validation error scores are often quite small, and elastic net has inferential properties that make it more desirable that LASSO or ridge on their own. Is this the right way of thinking? Thanks.               cvglmnet2a <- function(x, y, family, nfolds, seed) {       set.seed(seed)       cvfit0 <- cv.glmnet(x, y, family = family, alpha = 0, nfolds = nfolds)       cvfitp5 <- cv.glmnet(x, y, family = family, alpha = 0.5, nfolds = nfolds)       cvfit1 <- cv.glmnet(x, y, family = family, alpha = 1, nfolds = nfolds)       plot(cvfit1); plot(cvfitp5); plot(cvfit0)       plot(log(cvfit1$lambda), cvfit1$cvm,            pch = 19, col = "red", xlab = "log(Lambda)", ylab = cvfit1$name)           points(log(cvfitp5$lambda), cvfitp5$cvm, pch = 19, col = "grey")           points(log(cvfit0$lambda), cvfit0$cvm, pch = 19, col = "blue")           legend("topleft", legend = c("alpha = 1", "alpha = 0.5", "alpha = 0"),                  pch = 19, col = c("red", "grey", "blue"))           mins <- c(alpha1 = min(cvfit1$cvm),                 alphap5 = min(cvfitp5$cvm),                     alpha0 = min(cvfit0$cvm))       mins     }


title: Question about using a multiplicative dummy variable

text: In many econometrics model, the changes in the response variables in certain intervals are more difficult than other intervals. But I believe this is often not considered when estimating the model. For example, suppose $Y_{st}$ represents the proportion of students in a certain school $s$, passing a standardized test in year $t$. Let $R_{st}$ be the academic resources students (ex. books in library), and $I_{st}$ represent average parental income of the students. In this case $Y_{st} \in [0,1],$ and we would like to estimate effect of $R_{st}$ on $Y_{st}.$ We could model this is as follows, $Y_{st} = \alpha_{0} +\alpha_{1}R_{st} + \alpha_{2}I_{st} + \delta_{t}+ u_{st}$, where $u_{st}$ is additive error term, and $\delta_{t}$ are time dummies. In this context of pass rates, intuitively it is more difficult for a school to increase the pass rates from 95% to 100%, then it is for them to go from 45% student passing, to 50% student passing. Consequently, the effect of $R_{st}$ on $Y_{st}$ should be given less weight on the latter situation (45% to 50%), than the former (95% to 100%). Suppose we were comparing two schools in which the same $R_{st}$ increase lead to these results, clearly the 95% to 100% school invested more efficiently. My idea is to use a multiplicative dummy variable with $R_{st}$, $\beta_{t}$, where $\beta_{t}$ takes on different values depending on the initial value of $Y_{st}.$ Is there a standard way to take this into consideration in the model? Are there other additional factors that could improve this model?


title: Grouping sequential values in R

text: I have a list of numbers I need to group by similarity (differences being 1 between each). For example, in a list of [198, 202, 207, 218, 219, 220], 190 would be put into a list, 202 would be put into another list, 207, into a list, and 218 219 220 into another list.


title: How to treat this OLS based on residual diagnostics

text: I am struggling already a couple of days with this simple OLS, can you help? Outcome years in function of predictor score, very simple linear model. The residual plot does absolutely not look good though. Is it correct that based on the residual plot versus the outcome variable, to say that "if I predict the outcome to be 12, I most of the time over estimate the years of eduction?" In blue my fitted OLS, in red a LOESS curve. What am I doing wrong or how to improve OLS? ![enter image description here](http://i.stack.imgur.com/YNXxI.jpg) I have tried to both transform the predictor and the outcome, no luck. Is there something else you can suggest? The specific problem I have is that I am making wrong predictions when the outcome age is 12 years (ref table below, vertical are my predictions, horizontal the true values). How to solve this issue? ![enter image description here](http://i.stack.imgur.com/IRIKi.png)


title: Comparing the performance of two teams measuring many variables once

text: Which statistical tools are best suited for this problem: Team A goes and measures 3 magnitudes in 200 locations (let's imagine the volume of a room, the insulation coefficient and the average temperature outside the room). For each room they calculate the energy this room will consume to keep the desired temperature for a week. So after the measurement campaign they have measured these 3 magnitudes in 200 different rooms (each room is only measured once). So the data looks like:                (Room_volume, r_coefficient, T, estimated_energy)      Team A is assumed to be very good at getting their results. Now team B measures the same 200 rooms, and calculates all the estimated energies for each room. How should we measure the performance of team B vs team A (where A is assumed to be accurate)? Forgetting the physics, as it's just an illustrative example, how do you compare the accuracy of a measurement campaign measuring 3 values from which a 4th one can be calculated deterministically? Note that none of these values are distributed in any special way among the 200 samples (they don't look gaussian) I started looking at the distribution of errors, so:                {Estimated_energy_Ai - Estimated_energy_Bi} with i = 1..200      And seeing how many results are within a standard deviation. What other tools are appropriate?


title: Using descriptive statistics as input in SPSS

text: I want to perform a repeated measures ANOVA in SPSS using the descriptive statistics. So my input is:               Descriptive Statistics                       Mean    Std. Deviation  N     M1MI    3,8000  1,03280         10     M1MA    5,3000  2,16282         10     M2MI    7,0000  1,88562         10     M2MA    2,2000  1,54919         10     M3MI    6,2000  1,03280         10     M3MA    4,2000  1,75119         10     M4MI    4,7000  ,67495          10     M4MA    4,9000  1,19722         10      Does anyone know how I have to adjust my syntax? My design is as follows:               Dependent variables: Progress1, Progress2, Progress3, Progress4           Within subject factor 2: M (I/A)


title: Conditioning on independent random variables

text: I am in a situation where I have to compute: $$E(u(x_1)|\bar{X},S^2)$$ where $X_1$ is a normally distributed random variable and $u(.)$ some function. I know that by the student's theorem the sample mean and the sample variance are independent and moreover that $\frac{(n-1)S^2}{\sigma^2}\sim \chi^2 (n-1)$. Can I simplify the expectation with the information I possess? Is perhaps the bivariate normal distribution of use here? Thanks. **EDIT** : Yes, an iid sample on $X$ is assumed here, hence the subscript $1$ on $u(x_1)$. If the expectation cannot be simplified, what is the conditional distribution of $X_1$ given $\bar{X}$ and $S^2$?


title: How can I predict a discrete feature using regression?

text: I have a data set with ~80 records, with ~8 features. I want to predict one of the features in future records. The feature is numeric and discrete. It ranges between -30 up to 140 with steps of 5. Until now I wanted to predict another feature which is boolean, so I used logistic regression. Which method should I use here? Maybe some kind of particle filter? Thanks!


title: What does the Semivariance tell me?

text: I am looking at a Semivariogram. I know it shows me the relationship between distance and semi-variance. I also know that at the end of the range the distance no longer auto correlates. What I am wondering is, what does the semi-variance tell me at the point where to distance no longer auto correlates (at the sill)? So in my case at about 900 m the distance no longer auto correlates. The sill at this value is about 4300. What does this value tell me?


title: Estimating repeat shoppers from an incomplete sampling

text: I'm trying to estimate how many people visited the farmers market once, twice, thrice, etc. in a given time period, using sampled data. We have interview data from approximately 50% of visitors as they entered the market which lets us identify them uniquely. For the purposes of this analysis, I'm assuming that the interviews randomly captured ~50% of visits; in reality the interviews were not random, but I want to start with a simpler problem first. Within this dataset, I can identify how many people visited once, twice, thrice, etc. But, intuitively, I think the nature of sampling will lead me to underestimate the number of people who visited more than once. I've tested this intuition by randomly cutting the dataset in half a few times - starting from the full dataset (50% of visits), going down to 6.25%, and find that the larger datasets have more multiple-visit people in it (see below). However, I am unsure what happens as I go from 50% of the visits to 100% of the visits. Can you help me come up with a statistical framework to do that projection? PS - I feel that the Birthday Problem is informative here, but I can't think of how to apply it!                                    % of visits from people who visited at most.     % of visits sampled  once     twice    thrice     6.25%                83%      96%      99%     12.50%               82%      96%      99%     25%                  77%      94%      98%     50%                  67%      88%      95%     100%                 ?        ?        ?      I'm trying to think about how to apply the mark-recapture framework Gael mentioned. I agree that there is a similarity if I simplify it to the percentage of all visits from individuals who visited only once -- in essence, the population size in terms of this framework. What I'm struggling with is how to think about my sampling technique. I know how many total visits there (say, 25,000). I know that out of the 12,500 visits we "marked" in a continuous 50% sampling of visits, there were about 10,000 unique individuals, with 7,000 being marked once, and the remainder being "marked" 2, 3, or more times. I can't fit this into the simple two-stage mark-recapture model, and I'm thinking I need to do some kind of a Poisson regression (based on the wikipedia entry).


title: Set seed before each code block or once per project?

text: It is standard advice to set a random seed so that results can be reproduced. However, since the seed is advanced as pseudo-random numbers are drawn, the results could change if _any_ piece of code draws an additional number. At first glance, version control looks to be a solution to this, as it would at least allow you to go back and reproduce the version extant when you wrote down the results in your notes or paper. However, since it only takes one draw to mess things up, if you update R the results could change as well. I realize that this is probably only problematic in rare cases, but I'm curious if there are any best practices here. This is something I've been struggling with in my own work.


title: Probability of winning given competitor pair analysis

text: Given an analysis of every pair of competitors in a race, how may I determine the probability of any given competitor winning the race? For example, what is the probability of competitor 2 winning the following race? P(Cx Win) means the probability of Competitor x winning.                Cx   Cy   P(Cx Win)  P(Cy Win)     ----------------------------------      1    2       0.3        0.7      1    3       0.4        0.6      1    4       0.9        0.1      2    3       0.8        0.2      2    4       0.7        0.3      3    4       0.9        0.1      I have tried to calculate a 'rating' for each competitor by adding their individual win probabilities. For example the rating for C1 would be 1.6. The rating for C2 would be 2.2 etc. I've tried different ways to use this rating to find the probability of the competitor winning, however my gut feel tells me something is wrong. Is there a mathematical solution to this problem?


title: How to perform diallel analysis in R?

text: Diallel Analysis using the Griffing and Hayman approach is so common in plant breeding and genetics. I'm wondering if someone can share R worked example on Diallel Analysis. Is there any good referenced book which covered worked examples? Thanks References: Griffing B (1956) Concept of general and specific combining ability in relation to diallel crossing systems. Aust J Biol Sci 9:463-493 [pdf] Hayman BI (1954) The analysis of variance of diallel tables. Biometrics 10:235-244 [JSTOR] Hayman BI (1954) The theory and analysis of diallel crosses. Genetics 39:789-809 [pdf]


title: Simple machine learning: bot detection

text: I've been aching to get my feet wet with a machine learning project, and I've found one that should be relatively simple, and actually has non-negligible business value for my organization. The marketing guys have to remove bot activity from our tracking data by hand for their metrics. I wanted to pull some data from GA, and have them construct a data set (bot, not-a-bot). There are probably 5-10 (numerical) categories that we have to train the algorithm, and the data set can be made as big as the marketing guys have an appetite for. I've done a bit of reading, and played with RapidMiner/Knime/Weka a bit. I plan to do everything in Python, with `scikits-learn`, possibly working in R where I have to. My questions:   1. Is this a "not actually that easy at all" problem?   2. Given the number of categories, about how large should the training set be?    3. Given the problem, what algorithms should I start with?    4. Has anyone else done any learning around bot detection? How did it work? Am I barking up the wrong tree? Thanks in advance community!


title: Parameter estimate for linear regression with regularization

text: For given cost function $S(\beta) = (Y - X \beta)^T(Y - X \beta) + \lambda \beta^T \beta$, where $\lambda$ is regularization parameter, the $\beta$ that minimizes the given cost function is $\beta = [X^T X + \lambda]^{-1} X^T Y$. is it right?


title: Example of classification dataset where SVM with linear kernel performs well

text: I'm looking for a dataset (preferably with a story, at any rate a real dataset) where a SVM with a linear kernel performs well...in other words i'm looking for a dataset where the class boundary is likely to be linear. Ideally, there should be between 10-20 continuous variables and not too many discrete ones.... It's not really important to me how many classes there are (as long as it's a classification problem). Also, it can't be the iris dataset. Any proposals? P.S.: @modo: i don't think this is a question about obtaining a particular dataset: it's more that i have not read svm/machine learning papers since i was a master student so i'm not familiar with some comon good examples...i'm just sure they exist.


title: Comparing logistic regression models with the same IVs

text: I have multiple logistic regression models with all of the same IVs/controls and a variety of DVs (all health outcomes from the same sample). The primary IV is the sum of types of childhood abuse (emotional, physical or sexual). I made dummy variables that represent any one type of experience, any two types of experiences, or all three types of experiences (so each is mutually exclusive). This is the same type of model the CDC uses for their ACEs study which is where I borrowed the method from. Question 1: Can I compare the one experience dummy to the two experience dummy within the same model? That is, talk about the odds ratios in comparison to one another without standardizing the coefficients? My sense is yes and I've seen it done all over the place but I recently was given a dissenting opinion saying that since I am only comparing each IV to the dummy referent of 0 experiences, I can't compare them to one another without standardizing first. Question 2: What's the best method to make comparisons across models (with all the same IVs)? I'm testing the dummy IVs against a variety of physical and mental health outcomes and I'd like to compare the odds ratios for each DV based on any one type of experience, two experiences or three experiences. It would be nice to say, one experience increases the odds of this outcome by 3.2 times, this outcome by 2.1 times, etc. Therefore, I can say that one type of abuse increases the risk of depression more than anxiety disorder or two types of abuse increases the risk of PTSD over depression etc (assuming no overlap in confidence intervals). I've read Menard's 2011 piece on standardized LR coefficients and that makes sense as to what mechanism to use within a single model (as I would apply in question 1 if necessary), but I can't tell if this can be applied across DV models if I'm using all the same IVs/controls from the same sample. If I standardize each IV coefficient, then are they comparable across models? It's a random sample and each model has the same number of valid cases (1073) with no missing data.


title: Trend line does not seem to fit data

text: I'm trying to do a simple scatterplot and trend line in R, but it doesn't look right. Have I messed up something blatantly obvious? Any ideas as to why the line doesn't fit the data? ![the trend line does not fit the data](http://i.stack.imgur.com/AISqZ.png) Here is the code I used.               > plot(x, y, pch=".")     > model <- lm(x ~ y)     > summary(model)          Call:     lm(formula = x ~ y)          Residuals:          Min       1Q   Median       3Q      Max      -0.23043 -0.04340 -0.00533  0.03761  0.47882           Coefficients:                 Estimate Std. Error t value Pr(>|t|)         (Intercept) 0.051154   0.001384   36.97   <2e-16 ***     y           0.462881   0.003739  123.80   <2e-16 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 0.06365 on 71514 degrees of freedom     Multiple R-squared: 0.1765, Adjusted R-squared: 0.1765      F-statistic: 1.533e+04 on 1 and 71514 DF,  p-value: < 2.2e-16           > abline(model, col="red")


title: Paired permutations for three samples?

text: I am using two-samples paired permutation tests with Matlab. Now I have three within subject steps, forcing me to use a three-samples paired permutation test. Is there anything like that or any permutation test in general, I could use for that?


title: How to generate correlated test data that has Bernoulli, categorical, and continuous vectors (in R)?

text: I'm looking to generate a set of 5 random variables and enforce a dependence structure between them and onto a dependent variable Y. I understand how to generate correlated random variables for multivariate normal, but not when mixing different types. Below is a little more than I need, but I'm hoping someone can give me a general way of solving this problem...   * X_1 and X_2 need to be highly correlated Bernoulli variables.    * X3 needs to take one of 5 categorical values, call them "A"..."E".    * X4 needs to be normal, and negatively correlated with X1, X2.    * X5 needs to approximate test scores from 0 to 100 with a high skew, so gamma probably. X5 needs to be positively correlated with X1, X2, X4. Each of these variables must impact a "success/occurrence" Bernoulli distributed variable Y. How would I begin? I would like to enforce correlation both between the values of X, and also between each X and Y. (The categorical correlations seem particularly confusing to me.)


title: How should I set up my data for a logit model?

text: Let's say I have a sample of "events" done by a certain number of subjects, and some (although not most) of these subjects have underwent more than one event. I'd like to fit a logit model to these data to find out which characteristics contribute to a subject, over the course of an arbitrary period of time, undergoing another event within an arbitrary interval after a preceding event -- let's say three days. The way I see it, I can set my data up in one of two ways:   * take all of the events underwent by each subject, and define the dependent/indicator variable to be yes/no to "did this subject have another event within 3 days of a preceding event, over the course of the last year (for example)?" In this case, each row of the data for the model would correspond to every subject and be an aggregate of their event history.   * partition each subject's event history into pairs and define the dependent variable to be yes/no to "did this subject have another event within 3 days of this particular event?" In this case, each row would correspond to every event in the data. However, there are three outcomes: 1) the subject has another event within 3 days; 2) subject has another event within greater than 3 days; and 3) subject doesn't have any further events. Could I collapse 2) and 3) into one outcome, or would I be better off considering using a multinomial logit model instead? I'm leaning towards going with the latter option, but I'm worried that those subjects who undergo relatively more events will be over-represented in the final model (how would I deal with that, if I should?). However, I do think I gain significantly more information that way. Anyway, with that said, I'd love to hear insights on the pros/cons of each approach to setting up these particular data for a logit model.


title: Issues with using glmnet package for MATLAB

text: I am trying to use the `glmnet` MATLAB package to train my elastic net model on some huge data. My features are of size 13200, and I have around 6000 samples of these. I directly tried to use `lassoglm` in MATLAB with these features and corresponding target taking cross validation to just 3 folds and alpha = 0.5. It's already 6 hours and it hasn't finished. I have to do it for several others as well. Any suggestions what I should do?


title: Bayesian updating using $n$ noisy observations of Brownian motion

text: I am very new to Bayesian inference and can't figure out what may be an elementary problem. Also, please forgive me if I am screwing up the notation -- this is my first foray into Bayesian statistics. **Set-up:** At time $i=0$ I start with a random variable $X_0\sim \mathcal{N}(\mu_0,\sigma^2_0)$. Over time I observe the Brownian motion process such that for $j>i$, $X_j-X_i\sim \mathcal{N}(0,(j-i)D)$. Starting at time 0, I collect a sequence of $n$ observations $\{Y_i\}_{i=0}^{n-1}$ of $\{X_i\}_{i=0}^{n-1}$, which are subject to Gaussian noise such that $p(y_i|x_i)=\frac{1}{\sqrt{2\pi N_0}}e^{-\frac{(y_i-x_i)^2}{2N_0}}$. Noise is independent from observation to observation. I am interested in the mean and variance (of the estimator) of $X_{n-1}$ given $\{Y_i\}_{i=0}^{n-1}$ in terms of $n$, $\mu_0$, $\sigma^2_0$, $D$, and $N_0$. **What I've done:** Using Bayes rule and the properties of Gaussian distribution, I've found that for the first observation ($n=0$): $$p(x_0|y_0)=\frac{1}{\sqrt{2\pi S}}e^{-\frac{(x_0-M)^2}{2S}}$$ where $M=\frac{y_0\sigma^2_0+\mu_0N_0}{\sigma^2_0+N_0}$ and $S=\frac{\sigma^2_0N_0}{\sigma^2_0+N_0}$. I am not completely sure if the above is correct, but I am pretty confident that it is, as it looks like the mean squared error of the estimate of $X_0$ has decreased with the observation. However, I am having trouble extending this to $n>0$ and would appreciate any guidance. **EDIT:** I've misspecified the Brownian motion process. I think it's correct now.


title: WinBUGS error with zero values in binomial distribution: value of order of binomial <expr> must be greater than zero

text: It seems that WinBUGS has problems **if it has only zero draws from one binomial distribution** : **1. case - simple model**               for (i in 1:sites) {         N[i] ~ dpois(lambda)         for (j in 1:sample) {             y[i, j] ~ dbin(p, N[i])         }     }      If all values `y[i,]` are zero, then the following error message appears: > **value of order of binomial y[53,1] must be greater than zero** **EDIT:** note that `N[i]` was not zero in generated data upon failure. **2. case - advanced removal model**               for (i in 1:sites) {         M[i] ~ dpois(lambda)                  for (j in 1:sample) {             y[i, j, 1] ~ dbin(p, M[i])                  after_removal[i, j] <- M[i] - y[i, j, 1]             y[i, j, 2] ~ dbin(p, after_removal[i, j])         }     }      Here it seems that if some `y[i, j, 2]` is zero then the error message (same as in case 1) appears! Note that in this case, I found a workaround: I reparametrized the model to different one which **I hope** is equivalent:               for (i in 1:sites) {         M[i] ~ dpois(lambda)                  for (j in 1:sample) {             y[i, j, 1] ~ dbin(p, M[i])             obs[i, j] <- y[i, j, 1] + y[i, j, 2]             obs[i, j] ~ dbin(p_komb, M[i])         }     }          p_komb <- p + (1 - p) * p      ... anyway an explanation, clean solution and solution for case 1 wanted: ## Why this happens? Is it a WinBUGS bug? How can it be overriden? I have WinBUGS 1.4.3 (August 2007) with immortality patch installed. Below is complete **reproducible code for R and package R2WinBUGS (with data generation)** : **1. case**               require(vcd)          sites <- 120 # 60          mean_N <- 16          N <- rpois(sites, mean_N)                p <- 0.4          sample <- 3 # 3          y = matrix(nrow = sites, ncol = sample)     for (i in 1:sites) {         y[i,] = rbinom(sample, N[i], p)     }     y[20,] = 0          sink("tmp_bugs/model.txt")     cat("          model {          # likelihood     for (i in 1:sites) {         N[i] ~ dpois(lambda)         for (j in 1:sample) {             y[i, j] ~ dbin(p, N[i])         }     }          # derived parameters     Ntot <- sum(N[])          # priors          p <- 1/(1+exp(-logit_p))     tau <- 1/(4 * 4)     logit_p ~ dnorm(0, tau)          lambda ~ dunif(0, 100)          }               ")     sink()          win.data = list(y = y, sample = sample, sites = sites)     #win.data = list(y = y)          #inits = function () { list(N = apply(y, 1, max), p = mean(apply(y, 1, mean)/apply(y, 1, max))) }     inits = function () { list(N = apply(y, 1, max), logit_p = rnorm(1, 0, 4)) }               params = c("N", "p", "Ntot", "lambda")          ni <- 500     nt <- 12     nb <- 200     nc <- 3               out <- bugs(win.data, inits, params, "model.txt",         nc, ni, nb, nt, bugs.directory = bugs.dir,          working.directory = paste(getwd(), "/tmp_bugs/", sep = ""),         debug = TRUE     )          print(out, dig = 3)          par(mfrow = c(2, 2))     hist(out$sims.list$p, breaks = 100)     abline(v = out$mean$p, col = "red", lwd = 2)     abline(v = p, col = "green", lwd = 2)     #lines(quantile(out$sims.list$p, c(0.025, 0.975)), rep(sum(par("usr")[3:4]*c(0.9,0.1)), 2), lwd = 2)     lines(quantile(out$sims.list$p, c(0.025, 0.975)), rep(par("usr")[3], 2), lwd = 4)     legend("topleft", c("estimated", "real", "95% cred. int."), col = c("red", "green", "black"), lty = 1, box.lty = 0, cex = 0.7)          #hist(mean(y) / out$sims.list$p, breaks = 1000, xlim = c(0, 50))          hist(out$sims.list$Ntot, breaks = 100)     abline(v = out$mean$Ntot, col = "red", lwd = 2)     abline(v = sum(N), col = "green", lwd = 2)     #lines(quantile(out$sims.list$Ntot, c(0.025, 0.975)), rep(sum(par("usr")[3:4]*c(0.9,0.1)), 2), lwd = 2)     lines(quantile(out$sims.list$Ntot, c(0.025, 0.975)), rep(par("usr")[3], 2), lwd = 4)     legend("topright", c("estimated total N", "real total N", "95% credible int."), col = c("red", "green", "black"), lty = 1, box.lty = 0, cex = 0.7)      **2. case** (one of the `y[i, j, 2]` will likely be zero; if not, it can be assigned directly)               require(vcd)          sites <- 120 # 60          mean_M <- 16          M <- rpois(sites, mean_M)           p <- 0.4 #0.64          sample <- 2 # 3          y = rep(NA, sites * sample * 2)     dim(y) = c(sites, sample, 2)          for (i in 1:sites) {     #   obs[i,] = rbinom(sample, M[i], p)         for (j in 1:sample) {             y[i,j,1] = rbinom(1, M[i], p)             y[i,j,2] = rbinom(1, M[i] - y[i,j,1], p)         }     }          y_sample_total = apply(y, c(1, 2), sum)     obs = y_sample_total          ############################################          sink("tmp_bugs/model.txt")     cat("          model {          # likelihood     for (i in 1:sites) {         M[i] ~ dpois(lambda)                  for (j in 1:sample) {             y[i, j, 1] ~ dbin(p, M[i])                  after_removal[i, j] <- M[i] - y[i, j, 1]             y[i, j, 2] ~ dbin(p, after_removal[i, j])          }     }          # derived parameters     Mtot <- sum(M[])          # priors          tau <- 1/(4 * 4)          p <- 1/(1+exp(-logit_p))     logit_p ~ dnorm(0, tau)          lambda ~ dunif(0, 100)          }     ")     sink()               win.data = list(y = y,          #obs = y_sample_total,          sample = sample, sites = sites)          inits = function () { list(         M = apply(y_sample_total, 1, max),          logit_p = rnorm(1, 0, 4)     ) }          #params = c("M", "M2", "p", "p2", "Mtot", "M2tot", "lambda", "lambda2")     params = c("M", "p", "Mtot", "lambda")          ni <- 2500     nt <- 16     nb <- 1000     nc <- 3          date()     out <- bugs(win.data, inits, params, "model.txt",         nc, ni, nb, nt, bugs.directory = "C:/Program Files/WinBUGS14/",          working.directory = paste(getwd(), "/tmp_bugs/", sep = ""),         debug = TRUE     )     date()                    #############################          print(out, dig = 3)          par(mfrow = c(2, 2))          hist(out$sims.list$p, breaks = 100)     abline(v = out$mean$p, col = "red", lwd = 2)     abline(v = p, col = "green", lwd = 2)     #lines(quantile(out$sims.list$p, c(0.025, 0.975)), rep(sum(par("usr")[3:4]*c(0.9,0.1)), 2), lwd = 2)     lines(quantile(out$sims.list$p, c(0.025, 0.975)), rep(par("usr")[3], 2), lwd = 4)     legend("topleft", c("estimated", "real", "95% cred. int."), col = c("red", "green", "black"), lty = 1, box.lty = 0, cex = 0.7)          #hist(mean(y) / out$sims.list$p, breaks = 1000, xlim = c(0, 50))          hist(out$sims.list$Mtot, breaks = 100)     abline(v = out$mean$Mtot, col = "red", lwd = 2)     abline(v = sum(M), col = "green", lwd = 2)     #lines(quantile(out$sims.list$Mtot, c(0.025, 0.975)), rep(sum(par("usr")[3:4]*c(0.9,0.1)), 2), lwd = 2)     lines(quantile(out$sims.list$Mtot, c(0.025, 0.975)), rep(par("usr")[3], 2), lwd = 4)     legend("topright", c("estimated total M", "real total M", "95% credible int."), col = c("red", "green", "black"), lty = 1, box.lty = 0, cex = 0.7)               hist(out$sims.list$lambda)


title: All-Purpose Sample Entropy

text: When given samples of a discrete random variable, the entropy of the distribution may be estimated by $- \sum \hat{P_i} \log{\hat{P_i}}$, where $\hat{P_i}$ is the sample estimate of the frequency of the $i$th value. (this is up to a constant determined by the base of the log.) This estimate should not be applied to observations from a continuous distribution, at least naively, because it would yield a value which depends only on the sample size. Beirlant _et al_ describe a number of approaches for the continuous problem, including estimates based on empirical CDF, nearest neighbor distances and the $m$-spacing estimate, which is given by $$\frac{1}{n}\sum_{i=1}^{n-m}\log{(\frac{n}{m}(X_{(i+m)} - X_{(i)}))}$$, where $X_{(i)}$ is the $i$th order statistic of the sample, and $m$ varies in a certain way with $n$. It is not clear how this estimate is to be computed in the presence of ties, _i.e._ it does not appear to be applicable to discrete distributions. (a naive correction for ties (drop terms which have $\log{0}$) appears to give an estimator which does not depend on the relative frequency of the classes, only their values, which seems wrong.) The question: is there an 'all-purpose' estimator which can deal with both discrete and continuous (or even mixed) distributions?


title: More than one outcome (dependent) variables in ordinal logistic regression

text: I want to run ordinal logistic regression (OLR) in SPSS. My data include 6 predictor variable (two continuous and 4 categorical ) but my outcome variables are also 6 (categorical-likert scale). e.g my dependent variable is Business development and 6 likert scale questions were asked for this, increase in profit, sales, size, asset, marketing and labour. 1-If I composite these 6 variable into one by first sum all the 6 variables and then recode them into again five categories like before (In SPSS in Transform, Compute variable and then recode into different variables). I think in this way I lost my original data and may categorize them wrong. 2-If I think of using Categorical Principle component analysis for reducing the dependent variables, it will give results (object scores) in continuous form on which I have to run linear regression. I do not want to run linear regression as my original data is categorical. 3-So option left is to run OLR without combining the dependent variable on each DV. It means that there is no model for Business development and 6 models for profit, sales.... My question is that   * Is it preferable to run OLR with each DV and then summarize the results   * OR Is there any other method for reducing DV and running OLR.   * OR Is there any better method than ordinal logistic regression for this data.


title: Difference between Norm of Residuals and what is a "good" Norm of Residual

text: I am doing some basic fitting of data and exploring different fits. I understand that the residual is the difference between the sample and the estimated function value. The norm of the residuals is a measure of the deviation between the correlation and the data. A lower norm signifies a better fit. Suppose a cubic fit has a norm of residuals of 0.85655 and a linear fit has a norm of residuals of 0.89182. What if the norms are 0.17113 and 0.24916? Is the difference between these two significant? Is a norm of residual less than 1 considered good? If not what is generally regarded as "acceptable" norm of residual.


title: Null hypothesis with multiple equally important hypothesis

text: I am working on an assignment for a course in statistics, and a part of the assignment is to give the null hypothesis of an experiment. The experiment is as followed (short summary): A company is testing 3 interfaces for a website. They test the interfaces with 15 students, 5 for each interface, and all 5 have a different computer knowledge level. They are all given the same task and the time it takes for each of them to complete the task is the data we're using. I know need to find the null hypothesis of interest. As far as I understand null hypothesis and the alternative hypothesis should be 2 competing propositions. However, as far as I can tell I have 3 proposition: interface 1 is better, interface 2 is better and interface 3 is better. What null hypothesis should I have here and why? Thanks in advance!


title: Singularity issues when fitting an ARIMAX model

text: When I want to apply an ARIMA model like this:               arimax(dataMat, order=c(7,0,0), xreg=xreg)      it raises an error > "singular and could not compute QR decomposition of Hessian" for xreg. How can I fit proper model by using `xreg` values?


title: Post-hoc test in two-way ANOVA

text: I am doing a two-way independent ANOVA using type of organisational change (5 levels) and job level (7 levels) within the organisation as the factors (in SPSS). I followed Andy Field's procedure ( _Discovering Statistics Using SPSS_ , third edition) and I cannot manage to get the option of Post-Hoc tests for cases where equal variance is not assumed, i.e. Games-Howell. I searched the syntax help files and somewhere there's a note saying that when using more than one factor Games-Howell cannot be performed. However, Andy Field did it! (please see page 434 of his book). Any ideas how can I do the same?


title: Bayesian Expected Loss for mean F1-score loss function

text: So I have a multi-label classification problem where the exact number of labels in each test set example is unknown. The loss function is mean F1-score which is ![enter image description here](http://i.stack.imgur.com/5bMSS.png) where p is the precision and r is the recall. I also have posterior (beta) distributions for each of the possible labels. I am wondering how I would calculate the Bayesian expected loss for each possible label. And since it is a multi-label classification problem, I would also need to calculate the expected loss for all possible combinations of labels. Is this even a tractable problem?


title: Variable selection in time-series forecasting

text: I have a time-series forecasting task and would like some input on variable selection and regularisation. My problem has the following characteristics:   * 2,000,000 sample size.   * Most of the time, no change in the response. Change in the response happens 0.3% of the time and is what I'm interested in (don't talk about asymmetric loss functions and the like - I'm purely interested in variable selection/regularlisation).   * I will be predicting direction (up/down/no change).    * 100 predictor variables. Mostly hand engineered from complex non-stationary data. All carry some predictive power - the useless ones with no univariate forecasting power were eliminated early in the process.   * There will be clusters of highly correlated predictors - in general if you select two predictors at random the correlation will not be high.   * Predictors are heteroscedastic in nature.   * Some predictors are heavily serially correlated, many are not.    * Predictors all evolve continuously and are decimal processes. They are all bell shaped.    * Some could be non-stationary, but are thought to carry critical information about the evolution of the response process that would be lost upon differencing or similar $I(0)$ transforms.   * My forecasting model will be linear, but least squares loss function is not certain (perpendicular regression and quantile regression and discrete response regression will be tried).   * Standard error unbiasedness and variance is not a direct consideration. All I care about is out of sample forecast accuracy and model generalisation.  I know I can use anything from non-negative garotte, ridge, lasso, elastic nets, random subspace learning, PCA/manifold learning, least angle regression and various yucky hacks to select based on out of sample forecast performance. But is there anything specific to forecasting, or to the characteristics of my data, that would push me one way over another (time efficiency is important - trying everything and selecting the best is not practical). Note that decision boundaries on the forecasts are not intended to be part of the regularisation algorithm, but integration could be possible. I will also be using a supercomputer.


title: Deep learning with $p \gg N$ and on multi-modal and heterogeneous data

text: I work in a problem where:   * The input data is a collection of time series from different sensor modalities, where each sensor comes with different acquisition rates and dynamic ranges   * The number of potential features vastly exceeds the number of instances (i.e. $p \gg N$) My question is: **Has deep learning been applied on this type of data before?**


title: Repeated CV vs. CV

text: Suppose a data set contains the numbers $1-100$ (i.e. $\{1, \dots, 100 \}$). In $10$-fold cross validation, the data set is divided into $10$ subsets with one used as the validation data set. For example, one such division is $S_1 = \{1, \dots, 10 \}$, $S_2 = \{11, \dots, 20 \}, \dots, S_{10} = \{91, \dots, 100 \}$. In "regular" 10-fold cross validation, each subset is used as the test data set. In 10-repeated 10-fold cross validation (for example), the part that is repeated is how $S_1, \dots, S_{10}$ are constructed?


title: Multivariate model in lme() with independent random effect, similar to MCMCglmm

text: I would like to specify a multivariate model with `lme` with a random effect for group which is independent across variables. I found this post, which explains that the model specified as:               fit.multilevel <- lme( y ~ var -1,  dd,                         random = ~ var -1| school,                            correlation =  corSymm( form = ~ v |school/id),                          weights = varIdent(form = ~ 1 | v))       is equivalent to (using `MCMCglmm`):               fit.w.null <- MCMCglmm( cbind(mathach, ses) ~ trait -1 ,             random = ~us(trait):school,                rcov = ~us(trait):units,             data = hs,             family = c("gaussian","gaussian"))      However the model I would really like has independent random effects for school across the variables:               fit.w.null <- MCMCglmm( cbind(mathach, ses) ~ trait -1 ,             random = ~idh(trait):school,                rcov = ~us(trait):units,             data = hs,             family = c("gaussian","gaussian"))      (the random effect is specified with `idh()`). I cannot find a way to specify this in `lme`.   Here is some code to generate example data and run the analysis:               gr=rep(c("A","B","C"),each=100)     grp.mean=rep(rnorm(3,5,1),each=100)     var1=rnorm(300,rep(rnorm(3,5,1),each=100),1)     var2=0.8*var1+3*rnorm(300,rep(rnorm(3,5,1),each=100),2)     dat.mcmc=data.frame(var1=var1,var2=var2,group=gr)     dat.lme=data.frame(y=c(var1,var2),var=rep(c("1","2"),each=length(gr)),                  group=rep(gr,2),id=rep(1:length(gr),2))     dat.lme$v=as.integer(dat.lme$var)      Currently this is the best I can do:               fit.multilevel <- lme( y ~ var -1,  dat.lme,                            random = ~ var -1| group,                               correlation =  corSymm( form = ~ v |group/id),                             weights = varIdent(form = ~ 1 | v))       This is (I believe) equivalent to:               fit.w.null <- MCMCglmm( cbind(var1, var2) ~ trait -1 ,                             random = ~us(trait):group,                                rcov = ~us(trait):units,                             data = dat.mcmc,                             family = c("gaussian","gaussian"))      However, I would like an `lme` command equivalent to:               fit.w.null <- MCMCglmm( cbind(var1, var2) ~ trait -1 ,                             random = ~idh(trait):group,                                rcov = ~us(trait):units,                             data = hs,                             family = c("gaussian","gaussian"))      Any help you can give me will be much appreciated.


title: How can one set up a linear support vector machine in Excel?

text: Through the last year I have been working with support vector machines for a binary text classification task. Having used software such as R and Rapidminer I have not spent much time on understanding what actually goes on inside support vector machines. This I have now started looking into in the hope of getting a better understanding of this classification/regression method. I have spent a lot of time looking for calculation examples as it tends to enhance my understanding of a concept quite well if I can actually setup a problem in Excel. Therefore I hope to get guidance by asking this question here, as I have not been able to find any step-by-step calculation examples. One can easily find descriptions of the math and optimization problems one need to understand and solve computationally, but a step-by-step calculation example I have not been able to find. If the forum approves my idea of producing such an example I will do the editing and in the end produce a nice and clear Excel sheet and a guide for future use. I suggest that we use the _Iris_ dataset (even though it is a multiclass dataset) and simply try to separate _Iris setosa_ from _Iris versicolor_. I provide three links. Link one is theory of application of SVMs which I thought one could use as a scaffold. Link two provides a regression example of how I was thinking our product would look in the end. Link three will take you to the _Iris_ dataset. Theory and application of SVMs A guide for regressions Link for the Iris dataset Below I will try to formulate the problem more neatly. Problem description: > How can one apply Excel and the technique of a linear support vector machine > with soft margins in order to solve a binomial classification task given by > separating _Iris setosa_ and _Iris versicolor_ from the _Iris_ dataset using > all available features?


title: How to estimate true value and 95% bands when distribution is asymmetrical?

text: I have a set of results of independent measurements of some physical quantity. As an example I give here real expermental data on methanol refractive index at 25 degrees Celsius published in scientific literature from 1960 to 2011:               data={{1960, 1.32652}, {1961, 1.32662}, {1963, 1.32650}, {1963,        1.32750}, {1968, 1.32698}, {1968, 1.32890}, {1970, 1.32657}, {1970,        1.32660}, {1971, 1.3260}, {1971, 1.32610}, {1971, 1.32630}, {1971,        1.3350}, {1972, 1.32640}, {1972, 1.32661}, {1973, 1.32860}, {1975,        1.32515}, {1975, 1.32641}, {1976, 1.32663}, {1977, 1.32670}, {1980,        1.3250}, {1983, 1.32850}, {1985, 1.32653}, {1991, 1.32710}, {1995,        1.32621}, {1995, 1.32676}, {1996, 1.32601}, {1996, 1.32645}, {1996,        1.32715}, {1998, 1.32820}, {1999, 1.32730}, {1999, 1.32780}, {2001,        1.32634}, {2006, 1.32620}, {2011, 1.32667}};      The first number is the year of first publication, the second is the published value. The authors do not always provide their own estimate for error of the published value and even when they do this, the published estimate is often not accurate. It is easy to see that the distribution is asymmetrical (other datasets I work with are even more asymmetrical). Here is a density histogram with bins selected by hands:               bins = {1.325, 1.3259, 1.3261, 1.3263, 1.3265, 1.3266, 1.3267, 1.3269,         1.3275, 1.328};     Histogram[dataFirst[[All, 2]], {bins}, "PDF"]      ![histogram](http://i.stack.imgur.com/x4lne.png) Because of asymmetric feature of the distribution I cannot use the mean as an estimate for the true value of methanol refractive index. The fact is that the probability to get smaller value than the true value in general is not equal to the probability to get higher value. This means that the median in general is also unsuitable. It may be assumed that each invidual measurement has normal distribution for the instrumental error of the measurement. Additional sources of error are impurities in methanol. The most common impurity is water which increases the refractive index, but authors usually dry methanol before measurement. Probably this is the reason why the density is higher for low values: other impurities such as some common organic solvents (benzene) and (possibly) diluted gases (methane, carbon oxide, carbon dioxide) lower the refractive index. From the other hand alcohols, aldehydes and ketones which may present in industrial methanol and some organic solvents (CCl4) increase the refractive index. So there at least two sources of error: one (instrumental) may be assumed normally distributed and even equal for the all measurements and the second (impurities) is probably asymmetrical (and may be even bimodal) and depends on how methanol was obtained and purified. What is the best way to estimate the true value and its 95% confidence bands for the measured quantity in such cases? P.S. References for relevant scientific papers will be appreciated. P.S.2. The above code is for _Mathematica_ system.


title: Can I use 95%CI for mean for non-normal distribution if there is a natural limit?

text: I have problem with real basics. I have variable that has p=.000 in Shapiro- Wilk. As I looked on the histogram it appears that the cause is natural limit for this measurement (value 0, skewness > 1,5).I'm wondering what whould be the best way to present this data as descriptives. What parameters of descriptive & non-paramteric statistics would be the best to present data like this? Is there equivalent of 95%CI for non-normal distribution ? (for example interquartile range) Thank you in advance!


title: What do I need to consider when using the Hessian to compute S.E.'s?

text: I use optim() in R to do a lot of MLE. I've used the approach for a lot of problems, but the one I'm working on right now consists of fitting the parameters of the generalized extreme value distribution (GEV). There are three parameters (location scale shape), and I want to get the standard errors from this fit. An example of a paper doing this is Katz et al. 2005, and here is the online code for the article. I am unaffiliated with this publication. I should also point out that the GEV question is just an example of the type of approach I take -- most of my MLE fits are time series models (state-space). I could give more examples if needed. 1) Which of the arguments to optim() are of particular relevance to the statistical validity of using the optim() Hessian to computer standard errors?   * The first thing that comes to mind is the type of algorithm -- In R, optim() is a wrapper for several types of optimization algorithms (see above link for reference on optim()). I don't even think all of them return a Hessian.   * I also think that the parameters relating to the `reltol` and `abstol` parameters could be important, but I don't fully understand what they do. 2) Is there some pre-treatment of the data (e.g., standardizing data before doing the MLE) that is advisable if I want to use the hessian to compute S.E.'s of the fitted parameters? 3) In general, of what do I need to be careful. I just don't want to incorrectly use this approach. 4) Is there anything you can think of that I definitely **should not** do? Thanks.


title: Regressing a ratio on a component of the ratio

text: Is there anything wrong with regressing a ratio on itself? For example can we regress the log(savings ratio) on log(income) or the log(debt to income) ratio on log(income)? If not, should we use the log transformation to take income to the RHS and regress log(savings) on log(income)? If using a log transformation is the way to go, how do I interpret my results, given that the savings ratio is what I am interest in, not savings?


title: Calculate the average of hourly data of three sensors

text: I am trying to calculate the average of hourly data of three sensors but the hourly timestamps of all three sensors are different. How is it possible to measure the average of hourly data of all three sensors at a regular interval of hours say- 01:00,02:00,03:00 and so on. And also is it possible to start an extensible time series with a specific time which is before the first timestamp of the given dataset so that the extra timestamps have NA values and i can use na.rm=TRUE while calculating mean.:                 library(xts)      library(forecast)      data1$V1<-as.POSIXct(paste(data1$V1, data1$V2), format="%Y-%m-%d %H:%M:%S")           sensid<-as.factor(data1$V4)           #For Temp           for(i in 1:length(levels(sensid)))       {        l[[i]]=data.frame(data1$V6[data1$V4==i])        l[[i]]<-xts(l[[i]],order.by=data1$V1[data1$V4==i])        if(length(which (is.na(l[[i]])))<length(l[[i]]))          {                 wAvgt[[i]]<-apply.weekly(l[[i]], function(x) sapply(x,mean,na.rm=TRUE)) #Weekly   Average of Temp           dAvgt[[i]]<- apply.daily(l[[i]], function(x) sapply(x,mean,na.rm=TRUE)) #Daily Average of Temp           hrAvgt[[i]] <- period.apply(l[[i]], endpoints(l[[i]], on = "hours", 1), function(x) sapply(x,mean,na.rm=TRUE)) #Hourly Average of Temp               #For forecasting future values          aritw[[i]]<-auto.arima(wAvgt[[i]]) #Weekly Forecast ARIMA model using ARIMA           foretw[[i]]<-forecast.Arima(aritw[[i]],h=5) #Forecast of Weekly Average of Temp               arith[[i]]<-auto.arima(hrAvgt[[i]]) #Hourly Forecast ARIMA model using ARIMA           foreth[[i]]<-forecast.Arima(arith[[i]],h=50)#Forecast of Hourly Average of Temp               aritd[[i]]<-auto.arima(dAvgt[[i]]) #Daily Forecast ARIMA model using ARIMA           foretd[[i]]<-forecast.Arima(aritd[[i]],h=10)#Forecast of Daily Average of Temp         }     else         {          next         }      }      I have 27 sensors data. so i m trying to caclulate hourly,weekly and daily average of each sensors but each sensor provides 4 critical parameters like - blood pressure, temperature, Pulse and WBC. The above loop provided is only for temperature of 27 sensors. I have a bed connected to some 3 sensors and each sensor has 4 critical parameters to measure. I m concerned only with temperature for now. I need to calculate hourly average of temp for that bed only. So, initially i measure hourly average of any 3 given sensors each and then i intend to take the hourly average between these 3 sensors. The problem i m facing here is each hourly average has different number of rows and different timestamps. How to normalise all the timestamps? If the starting timestamp of sensor 1 is 2004-02-28 00:59:16, sensor 2 is 2004-02-28 01:08:23, sensor 3 is 2004-02-28 01:19:34, it is not possible to normalise them by any methods i know of yet. Please help me find the hourly average of three sensors.


title: Binary classifier probability measure?

text: I've the following situation. I've a binary classifier which classifies input feature vectors into either of two classes 'y' or 'n', along with the probability of it being in either of the classes P_y and P_n. I know that the classifier is right 60% of the time. How do I estimate the true probability of it being in each of the classes? What I'm doing now is to compute 0.6*P_y / (0.6*P_y + 0.4*(1- P_y)). This is a Bayesian adjustment, but I'm not sure if I'm doing the right thing. Any pointers much appreciated.


title: Why do Anova( ) and drop1( ) provided different answers for GLMMs?

text: I have a GLMM of the form:               lmer(present? ~ factor1 + factor2 + continuous + factor1*continuous +                      (1 | factor3), family=binomial)      When I use `drop1(model, test="Chi")`, I get different results than if I use `Anova(model, type="III")` from the car package or `summary(model)`. These latter two give the same answers. Using a bunch of fabricated data, I have found that these two methods normally do not differ. They give the same answer for balanced linear models, unbalanced linear models (where unequal n in different groups), and for balanced generalised linear models, but not for balanced generalised linear mixed models. So it appears that only in cases where random factors are included does this discord manifest.   * Why is there a discrepancy between these two methods?    * When using GLMM should `Anova()` or `drop1()` be used?    * The difference between these two is rather slight, at least for my data. Does it even matter which is used?


title: Outlier detection in very small sets

text: I need to get as accurate as possible a value for the brightness of a mainly stable light source given twelve sample luminosity values. The sensor is imperfect, and the light can occasionally "flicker" brighter or darker, which can be ignored, hence my need for outlier detection (I think?). I've done some reading up on various approaches here and can't decide on which approach to go for. The number of outliers is never known in advance and will often be zero. Flicker is generally a very large deviation from the stable brightness (enough to really mess with any average taken with a big one present), but not necessarily so. Here's a sample set of 12 measurements for completeness of the question: 295.5214, 277.7749, 274.6538, 272.5897, 271.0733, 292.5856, 282.0986, 275.0419, 273.084, 273.1783, 274.0317, 290.1837 My gut feeling is there are probably no outliers in that particular set, although 292 and 295 look a little high. So, my question is, what would be the best approach here? I should mention that the values come from taking the euclidean distance of the R G and B components of the light from a zero (black) point. It would be programmatically painful, but possible, to get back to these values if required. The euclidean distance was used as a measure of "overall strength" as I'm not interested in the color, just the strength of output. However, there's a reasonable chance that the flickers I mentioned have a different RGB composition to the usual output. At the moment I am toying with some sort of function that will repeat until a stable membership of allowed measures is reached by:   1. Finding the standard deviation   2. Putting everything outside say 2 SDs into an ignore list   3. Recalculating the average and SD with the ignore list excluded   4. Re-deciding who to ignore based on the new average and SD (assess all 12)   5. Repeat until stable. Is there any value in that approach? All comments gratefully accepted!


title: ARIMAX order stock price predictions

text: I am implementing the ARIMAX model for stock price predictions. To select the order, I looped through $0\le p,d,q\le10$ and found the best order to be (1,0,10). I'm a bit confused how to interpret these results. How do I tell if its reasonable, and what does this imply about the data?


title: How to choose prior in Bayesian parameter estimation

text: I know 3 methods to do parameter estimation, ML, MAP and Bayes approach. And for MAP and Bayes approach, we need to pick priors for parameters, right? Say I have this model $p(x|\alpha,\beta)$, in which $\alpha,\beta$ are parameters, in order to do the estimation using MAP or Bayes, I read in the book that we'd better _pick a conjugate prior $p(\alpha,\beta)$_ , which is a **joint** probability of $\alpha,\beta$, right? I have 2 questions:   1. Do we have other choices picking the prior other than this conjugate one?    2. Can we pick priors for $\alpha$ and $\beta$ respectively like $p(\alpha)$ and $p(\beta)$, other than put them together in a joint one?


title: count data: SD is higher than assumed for Poisson regression

text: I have a dependent variable (in counts, i.e., a number of technicians). My independent variable is a grouping factor (control and treatment). From the literature, the average number of technician required for the procedure is 3 with SD of 1. I would like to calculate a sample size to determine significant changes of 1 SD decrease in technicians. I think that these count data should follow a Poisson distribution and I could do a sample size calculation for Poisson regression. However, my SD is lower than one assumed for Poison regression (sqrt(mean)=1.73). What alternative tool should I use to determine a 1 SD difference and a sample size? Thanks! Julia


title: Probability of discrete events

text: An experiment has two possible outcomes: A or B. In 20 experiments, A happened 2 times, and B happens 18 times. What would I expect from a larger number of experiments? With what confidence can I say that the mean number of outcome A is 10%?


title: Where can I find a good resources for making publication quality line plots

text: I have to incorporate a number of line plots in a Latex document and I want to get a very high quality. Are there any website that show some templates or default plots styles that I can use. Myself I am not good in designing this kind of plot style, so I am searching some good basic settings for the tool I want to use. About the tools, I think I would like to use Tikz or gnuplot as I want to build tex document and plots in a single run. My data is anyway stored in text files.


title: How to calculate cex.axis automatically so to fit to the maximal width possible (in R)?

text: Let's say we created a plot, with labels:               plot(1, xlab = "text", ylab = "text times three") # text too small      In this plot there is still a lot of free region to enlarge the text. We could do so by using cex.lab, like this:               par(cex.lab = 2)     plot(1, xlab = "text", ylab = "text times three")      Questions:   1. Is there a way to know what is the largest cex.lab possible, without the text getting over the borders?   2. Is it possible to do that separately for each component of the axis? (so to have a different fitting of cex.axis to the xlab vs the ylab) ?


title: Use robust standard errors in confint

text: Despite this being of a program-related question I was asked to submit this here instead stackoverflow. Is it possible to have confint use the robust variance-covariance matrix obtained by vcovHC (from the sandwich package) after fitting a model? I\ve been looking around without finding it..


title: For what class of distribution is this function concave?

text: consider the following sorted sample: $x_1<x_2<\ldots<x_n$ and the kernel function: $h(j)=\log(x_j-x_1), j>1$. Now, it turns out that when I draw data from many continuous distributions, sort them and compute $h(j)$, the plot of the $h(j)$ as a function of $j$ is concave in this sense: $h((j+k)/2)>(h(j)+h(k))/2$ for $1<j<j+1<k<n$ My question is the following: will this observation (concavity of $h(j)$ as a function of $j$) be true for all continuous distributions? If not is there a class of distributions for which it will always be true? Thanks in advance, At this point, given the useful comments from @Glen_b and @Momo, I think I should give an example. this is a vector of values of $x_j-x_1,j>1$ (where the $x$'s are drawn from a continuous distribution, but have been rounded here to the third digit so that they don't take too much place.) My problem, is that when I plot $h(j)$ for this data-set, the result is clearly nor concave (in the sense I outlined above).               x=c(1.403, 1.406, 1.408, 1.416, 1.417, 1.42, 1.427, 1.441, 1.448,      1.456, 1.458, 1.465, 1.466, 1.472, 1.472, 1.477, 1.477, 1.479,      1.482, 1.491, 1.5, 1.504, 1.518, 1.52, 1.52, 1.53, 1.544, 1.545,      1.561, 1.573, 1.595, 1.595, 1.599, 1.603, 1.605, 1.612, 1.617,      1.618, 1.618, 1.628, 1.64, 1.644, 1.646, 1.653, 1.679, 1.682,      1.682, 1.693, 1.694, 1.71, 1.71, 1.711, 1.741, 1.75, 1.756, 1.773,      1.794, 1.799, 1.804, 1.808, 1.86, 1.882, 1.895, 1.955, 1.992,      1.995, 2.009, 2.009, 2.063, 2.123, 2.329, 2.356, 2.405, 2.535,      2.632, 2.635, 2.725, 2.763, 2.783)


title: Estimating risk parameters from observed choices

text: My background is in labor economics/applied micro but I am currently working with data collected from a field experiment. My experimental/micro skills are rusty to say the least. In my data, I observe individuals choosing one of the choices $A, B, C, D$. The value of the choice is realized after the choice. I know that $$E(A) <= E(B) <= E(C) <= E(D)$$ $$VAR(A) <= VAR(B) <= VAR(C) <= VAR(D)$$ If I assume that $A, B, C, D$ are continuous, independent and normally distributed and individuals are aware of expected payoffs $E(i)$ and variances $VAR(i)$.   * Is there a way to estimate a expected risk aversion parameter (e.g. CRRA) for a group of people?    * Would it make sense to compare these risk estimates from two groups of people, say men and women or high educated and low educated? I would appreciate if you guys could come up with some relevant literature.


title: Is there an advantage to rescaling a set of Likert items from 1-5 to 0-4?

text: I have a set of Likert type items asking respondents their attitude on different topics. For example, attitudes toward drugs/alcohol where 1=strongl disagree and 5=strongly agree to the statement "I think it is okay for someone my age to smoke marijuana." The 6 items in this particular scale are then scored using the mean of non-missing items. Is there some advantage to first rescaling the items so that the range of responses goes from 0-4 rather than 1-5? The only difference is in interpretation of the scale where the average "attitude" is 3.5 on a scale of 0-4 and 4.5 on a scale of 1-5. From an interpretation standpoint, what is easier for a lay audience? Or is there some standard of practice in terms of reporting?


title: What's wrong with the intuition that kernel measures similarity between observations?

text: Near the middle of page 16 of Andrew Ng's notes on SVM, he explained an intuitive view of kernel as measuring similarity between observations, but then added the caveat that > there are things wrong with this intuition, but nevermind The notes don't really explain what's wrong with this intuition. Does anyone know what he might be referring to?


title: How best to communicate uncertainty?

text: A massive problem in communicating the results of statistical calculations to the media and to the public is how we communicate uncertainty. Certainly most mass media seems to like a hard and fast number, even though except in a relatively small number of cases, numbers always have some uncertainty. So, how can we, as statisticians (or scientists describing statistical work), best communicate our results, while keeping the uncertainty in tact, and making it meaningful to our audience? I realise that this isn't actually a statistics question, rather a psychology question about statistics, but it's certainly something that most statisticians and scientists will be concerned about. I'm imagining that good answers might reference psychological research more than stats textbooks... * * * **Edit:** As per user568458's suggestion, a case-study may be useful here. If possible, please keep answers generalisable to other areas. The particular case that I'm interested in serves as a nice example: _the communication of climate science to politicians and the general public, through mass media_. In other words, as a scientist, it is your job to convey information to a journalist in such a way that they have little difficulty in accurately conveying that information to the public - that is, the truth, although not necessarily the _whole_ truth, which won't usually fit in a news- bite. Some particularly common examples might be the communication of the uncertainty in the estimate of degree of warming over the remainder of the century, or in the increased likelihood of a specific extreme weather event (i.e. in response to a "was this storm caused by climate change" type question).


title: Regression when you have both I(1) and I(2) processes

text: I want to estimate $\log Y_t=\beta_0+\beta_1\log X_{1t}+\beta_2\log X_{2t}+e_t$ Using ADF-tests, I have found that both $\{\log Y_t\}$ and $\{\log X_{1t}\}$ are I(1) (i.e. they contain a unit root), while $\{\log X_{2t}\}$ is I(2). I have checked that $\{\Delta \log Y_t\}$, $\{\Delta \log X_{1t}\}$ and $\{\Delta_2\log X_{2t}\}$ are all stationary. Now my question is what regression do I run using these differences, so that I can still say something about the original relationship above? Can I simply take second differences of $\{\log Y_t\}$ and $\{\log X_{1t}\}$ and estimate $\Delta_2\log Y_t=\beta_1\Delta_2\log X_{1t}+\Delta_2\beta_2\log X_{2t}+v_t$ or does taking the second difference of an I(1) series introduce some kind of problem?


title: How to design experiments for Market Research (with a twist)?

text: Consider a kind of auction where you are presented with, say, 1000 prospective clients. Based on information about these prospectives--age, sex, race, income, educational achievement and the like--you may 'bid' to pitch your product to some fraction of them, say 250. (Ignore the cost of bidding.) To maximize your chances of selecting the right subset, I would probably use a model of our product's 'likability' built using logistic regression. I can fumble my way through this part well enough. However, to _build_ the model of product likability, I have to do some market research, testing the product pitch on subjects that we may recruit from the population at large. This is actually fairly costly. Moreover, it should probably be tuned to the demographics of the population from which the prospectives are drawn. For example, a factorial design, say, might ask us to recruit test subjects in equal proportion across the levels of the race factor, when in reality we are likely to encounter very few Native Americans, say, in the 1000 prospectives, and can simply choose not to pitch to them at all as a rule. (Sad, but true.) How should such an experiment be designed? To be concrete, the design variables are all categorical and ordinal factors, the bidding fraction is an input parameter (1/4 in the example quoted above), as is the maximum number of subjects that can be recruited. It seems like maybe some mix of experimental design and random sampling might be appropriate, but I am open to all reasonable suggestions and pointers. I should also note that given the likely small effect sizes and the small sample recruitment pools that we can afford, it is unlikely that the market research will yield statistically significant regression coefficients. And so overoptimizing the experimental design is probably silly, and any reasonably non-insane procedure will suffice.


title: Measuring how close a given sample distribution is to an ideal distribution

text: I'm trying to fashion survey experiments for which the demographic composition of the samples needs to be as close as possible to a 'golden' composition. For concreteness, let's say the 'golden' composition is `14% middle class, 35% upper-middle class, and 51% upper class`. I have a whole bunch of samples, each with slightly different compositions. For example:               S1 = 10% middle class, 41% upper-middle class, 49% upper class     S2 = 33% middle class, 33% upper-middle class, 34% upper class      I'm trying to find a good statistical measure to describe how close the distribution in each sample is to the ideal, 'golden' composition. Is there a simple metric here that I could use?


title: Compare mean effect of one predictor based on another predictor

text: I would like to compare the effect of one IV ("sensation seeking") on a DV ("intended infidelity") based on another IV ("gender"). Actually I really just want to compare the means. So I would like to compare the mean effect of sensation seeking (a continuous variable with 5 levels) for males and females on intended infidelity (a continuous variable with 4 levels). Why can't I seem to do this in SPSS? Do I have to do a moderation analysis? Can I re-compute the variables and into what? Thanks very much for any suggestions!


title: What are the case studies in public health policy research where unreliable/confounded/invalid studies or models were misused?

text: I am drafting a literature review on a current public health issue where data are confounded: **What are common historical case-studies that are used in public health/epidemiology education where invalid or confounded relationships or inferences were intentionally or erroneously employed in public health policy and legislation?** The automobile fatality surge of the 1960s and subsequent evidence-based, government-led study which determined seatbelts and eventually airbags should be required by law is a great example of **HOW** public health policy should be driven by statistically powerful inferences and models. I am looking more for examples of cases of the opposite type (bad science to make policy hastily). However, if nothing else I would like to learn of more cases similar to the previous example of powerful studies for successful public health benefit. I want to use these as examples to demonstrate how evidence-based statistical public health research is important to the making of policy.


title: Best way to convey highly skewed data

text: I am estimating a linear probability model with a hierarchical structure. I am bootstrapping standard errors to test the significance of the coefficients in order to double check the parametric estimates of the standard errors. In the bootstrapping process I am finding that the coefficients are all positive (which I expected) but there are a few extreme positive outliers. The outliers make the standard errors huge and subsequent t-statistics small. This is misleading because an empirical p-value has a value under 0.01. My questions are   1. What is the best way to convey that this coefficient is statistically different from zero?    2. What would be the best way to report these results? Normally in a LaTex table I report coefficients plus standard errors.


title: Hypothesis test for population rather than for sample

text: I have a set of data for companies in a certain industry. The data set spans two consecutive years, 2000 and 2001. Let's say I have 100 unique companies in the dataset and I calculate average service time for each company and for each year. How do I test whether the average waiting time for company x in 2001 is less than the average waiting time for company x in 2000? The data set contains all service transactions performed by company x in 2000 and 2001, and this is what confuses me. The way I understand hypothesis test is that you test the null hypothesis using sample mean. However, with this data I have the population mean so what does it mean when your professor ask you to test whether the average service time in 2001 is less than the average service time in 2000 using "simple mean comparison test such as t-statistic"?


title: How to approach prediction of new observations with incomplete data from model built from complete data

text: I currently have a gradient boosting model that uses the `gbm` package in R that classifies observations at the end of a year. Daily behaviors are logged for each observation (with most of them zero). The model is built on observations which have a year's worth of complete data and performs well out of sample. For illustrative purposes, assume the complete data look like this:               id day1 day2 day3 day4 day5, ..., day365      1    1    1    2    0    0, ...,      0      2    3    0    1    0    0, ...,      0      3    1    1    2    0    0, ...,      0      4    1    0    0    0    0, ...,      1      However, I want to predict the probability of each class for new individuals who have not yet been in the sampling frame for a year. Thus, all of the days past the present day will be NA. This precludes the use of many algorithms, but treating them as 0s also implies having more information than I actually have. What is the best way to treat the unobserved, but observable missing data when making predictions? I will not be re-estimating the model on a daily basis but rather using the model fit from the existing training set.


title: Difference between quasi-likelihood estimating equations, IEE and GEE?

text: 1. What is the difference between quasi-likelihood estimating equations, and GEE? From a note > Quasi-score function is for independent, over-dispersed data (Poisson or > binomial), while GEE1 is for correlated data. But it seems to me that in the note, their equations both don't take into account the correlation between different responses. So how shall I understand them correctly?   2. How is independence estimating equations (IEE) different from the above two? Thanks!


title: Motivation behind random forest algorithm steps

text: The method that I'm familiar with for constructing a random forest is as follows: (from http://www.stat.berkeley.edu/~breiman/RandomForests/cc_home.htm) To build a tree in the forest we:   1. Bootstrap a sample of size N where N is the size of our training set. Use this bootstrapped sample as the training set for this tree.    2. At each node of the tree randomly select m of our M features. Select the best of these m features to split on. (where m is a parameter of our Random Forest)   3. Grow each tree to largest extent possible -- i.e. no pruning.  While this algorithm makes sense at a procedural level and certainly produces good results, I'm not clear what the theoretical motivation is behind the steps 1, 2, and 3. Could someone explain what motivated someone to come up with this procedure and why it works so well? For example: why do we need to perform step 1? It doesn't seem like we're bootstrapping for its usual purpose of variance-reduction.


title: How do I compute a prediction interval from generalized estimating equation output?

text: Suppose I have repeated observations on individuals, which I've analyzed using a generalized estimating equation (GEE). What is the procedure for computing the prediction interval for future observations given values of the predictors? I would appreciate either instructions or references. Thank you.


title: Are sampling weights necessary in logistic regression?

text: When should I use weights when performing a logistic regression? The weights I'm referring to are sampling weights from a survey? Or should I just use the unweighted data?


title: Information out of the hat matrix for logistic regression

text: It is clear to me, and well explained on multiple sites, what information the values on the diagonal of the hat matrix give for linear regression. The hat matrix of a logistic regression model is less clear to me. Is it identical to the information you get out of the hat matrix applying linear regression? This is the definition of the hat matrix I found on another topic of CV (source 1): $H=VX ( X'V X)^-1 X' V$ with X the vector of predictor variables and V is a diagonal matrix with $\sqrt{(π(1−π))}$. Is it, in other words, also true that the particular value of the hat matrix of an observation also just presents the position of the covariates in the covariate space, and has nothing to do with the outcome value of that observation? This is written in the book "Categorical data analysis" of Agresti: > The greater an observation’s lever- age, the greater its potential influence > on the fit. As in ordinary regression, the leverages fall between 0 and 1 > and sum to the number of model parameters. Unlike ordinary regression, the > hat values depend on the fit as well as the model matrix, and points that > have extreme predictor values need not have high leverage. So out of this definition, it seems we can not use it as we use it in ordinary linear regression? Source 1: How to calculate the hat matrix for logistic regression in R?


title: Modelling the tail only

text: I'm trying to model a real-world random variable that behaves approximately as a Gaussian, so a Normal distribution fit is reasonable but far from perfect. However, I only care about its tail, that is, I only care about finding $x$ such as $P(X > x) < p$, with $p$ anywhere between 1% and 5%. I'm modelling the random variable as $X\sim\mathcal{N}(\mu, \sigma^2)$, with parameters estimated from all the samples. I tried to search for specific ways to handle this and I got a lot of results regarding Extreme Value Theory, but it seems to apply mostly to even smaller probabilities than I care or to a heavy-tailed distributions, which is not my case. Any pointers?


title: How does the power of a logistic regression and a t-test compare?

text: Is the power of a logistic regression and a t-test equivalent? If so, they should be "data density equivalent" by which I mean that the same number of underlying observations yields the same power given a fixed alpha of .05. Consider two cases:   1. [The parametric t-test]: 30 draws from a binomial observation are made and the resulting values are averaged. This is done 30 times for group A (which has a binomial Pr of .70 of occurring) and 30 times for group B (which has a binomial Pr of .75 of occurring). This yields 30 means per group that represent a summary of 1,800 draws from a binomial distribution. A 58df t-test is performed to compare the means.   2. [The logistic regression]: A logistic regression is performed with a dummy coded slope representing group membership and each of the 1,800 draws. My question has two parts:   1. Given a set alpha of .05, will the power of these methodologies be the same or different? Why? How can I prove it?   2. Is the answer to question 1 is sensitive to the sample sizes going into the t-test, sample size of each group in the t-test, underlying binomial probabilities, or some other factor? If so, how can I know (without simulation) that the power is indeed different and what sort of changes will produce what sort of changes in power? Alternatively, provide worked out R code that solves the issue using simulation.


title: How to get coefficients from princomp object in R?

text: Basically, I want to manually compute what `predict(princomp(somedata))` produces. In order to do that I need to have coefficients matrix so that               newdata %*% coefficients      would produce the predicted values of each component. Can anybody give me a clue to get the coefficients?


title: Test randomness of a generated password?

text: As you know, there are many password generators out there to increase computer security. Suppose I am given such a password (say, a string of letters, numbers, symbols, etc.), is there a way for me to test how random it is? Is there some sort of index for this measurement? Thanks!


title: Combining principal component analysis and partial least squares

text: I know PCA and PLS are considered as alternative method to each other. But I am thinking about a kind of combination of the two in case of lots of predictors with little variability. In that case, when I run 1-component PLS with original predictors, it does not produce a meaningful model in terms of prediction. But if I first compute 10-20 PCA components and run 1-component PLS with those PC scores as predictors, practically, the model is quite good in terms of prediction power. But I would like to know why. Can anybody explain why this is better than 1-component PLS with original predictors?


title: Alpha Expansion algorithm over graph with vertices having different sets of labels

text: The alpha-expansion algorithm described here (Probabilistic Graphical Models: Principles and Techniques By Daphne Koller, Nir Friedman, Page 593: algorithm 13.5), is designed for graphs where each vertex has the same set of labels. I am interested in a more general problem where each vertex can have different set of labels having different cardinalities. Can the alpha-expansion algorithm be used in such cases? For example, let the graph have 3 vertices. The first vertex can assume labels (a1, a2, a3). The second vertex can assume labels (b1, b2) and the third one can have (c1, c2, c3, c4, c5). The node and edge potentials are defined accordingly. Can alpha-expansion still help me in finding an optimal (approximate) label combination?


title: How to help people understand "...per 100,000"?

text: It's well known that members of the public have difficulty comprehending numbers over a certain size. A common tricky issue for people communicating demographic trends is helping people make sense of rates out of 100,000. This is an important issue in helping people make sense of changing health and crime statistics, for example. For those of you who do statistical communication for the public, what frames of reference or other analogies and aids have you found to be successful in helping people to put rates out of 100,000 in perspective? There's another question with a discussion alluding to 100,000 people being a small to medium sized city, and it seems there's disagreement over whether this analogy is helpful or not. Answers citing evidence would be great, personal experience and ideas are helpful too. _(p.s. I'm alsostarting a discussion on meta on whether to tag/allow/encourage stats communications questions like this on this site)_


title: How can you make linear discriminant analysis reduce dimensions to the number of dimensions you are looking for?

text: Let's say I have a $m \times n$ matrix where $m$ is the number of points and $n$ is the number of dimensions. I would like to give a target dimension parameter which is let's say d. d can be a set of values like $\{2,4,6,\ldots,n\}$. I would then approach the problem using Fisher's linear discriminant analysis to give me output matrix in $m \times d$. Is this possible to do? I still don't understand how LDA reduces dimensions from let's say 10,000 to 1,000.


title: Absolute error loss minimization

text: From _Robert (The Bayesian Choice, 2001)_ , it is proposed that the Bayes Estimator associated with the prior distribution $\pi$ and the multilinear loss is a $(k_2/(k_1+k_2))$ fractile of $\pi(\theta|x)$. The proof follows that $$ E^\pi[L_{k_1,k_2}(\theta,d|x)] = k_1\int\limits_{-\infty}^{d} (d-\theta)\pi(\theta|x)d\theta + k_2\int\limits_{d}^{+\infty}(\theta-d)\pi(\theta|x)d\theta $$ Then, using the identity $$ \int\limits_{c<y} (y-c)f(y)dy = P(y>c) $$ I would get to $$ E^\pi[L_{k_1,k_2}(\theta,d|x)] = k_2P^\pi(\theta>d|x) - k_1P^\pi(\theta<d|x) $$ But he gets to $$ E^\pi[L_{k_1,k_2}(\theta,d|x)] = k_1\int\limits_{-\infty}^{d} P^\pi(\theta<y|x)dy + k_2\int\limits_{d}^{+\infty} P^\pi(\theta>y|x)dy $$ And then takes the derivative in $d$. What am I'm missing, and why he does this last step? Thanks!


title: Does a regression tree strictly dominate OLS in prediction?

text: Since OLS tries to measure E[Y|X], and regression trees try to partition the data into different branches, then take means of the response under different branches, is it reasonable to say that regression trees strictly dominate OLS in prediction? Also, since OLS is not robust to outliers, would it be reasonable to say that regression trees are a better first, off-the-shelf algorithm to try in a prediction problem? For interpretation, it definitely makes sense to me that knowing $R^2$, regression coefficients, p-values, etc. are more useful than seeing the tree structure.


title: Robust standard errors in multiple regression

text: I use Andrew F. Hayes' macro for SPSS (HCREG at http://www.afhayes.com/spss- sas-and-mplus-macros-and-code.html) to perform multiple regression analyses with robust standard errors. The information I get in the output is limited compared to the one SPSS provides when applying multiple regression. For example, there are no standardized coefficients (Beta) for predictors and there is no adjusted R-squared. Can I assume that the standardized coefficients will be the same as in the model without robust standard errors? If not, does anyone know how can I go about computing them by hand? I have the same question for the adjusted R-squared. Any help would be greatly appreciated! Thank you very much in advance.


title: What is bias correction?

text: I have seen many places where they have input/output datasets where they first create a linear regression line, correct the bias, and then only use that data for their model. I didn't get what this bias correction is?


title: coin flipping distribution and probability

text: A person is trying to claim that psychics exist by flipping a coin 10 times for each person behind a screen. n=10 p=1/2 Its created into a one tailed test with significance level at .001. The first 999 people do not get enough right to reject the null hypothesis. But the 1000th person gets all 10 right. if Y is the number of people out of 1000 who get all 10 coin flips right, would the distribution of y be binomial too? Probability that one or more of the 1000 people would correctly guess 10 coin flips?


title: How to calculate expected value and standard deviation if I have 100 values divided into 15 groups (normal distribution)?

text: These are the values:               18.2    17.6                                                     ------------------------------------------------------------------------------------------------------------     18.9    18.6                                                     ------------------------------------------------------------------------------------------------------------     20  19.7    19.8    19.6                                             ------------------------------------------------------------------------------------------------------------     21.4    20.6    21                                               ------------------------------------------------------------------------------------------------------------     21.9    22.1    22.2    22.2    21.6                                         ------------------------------------------------------------------------------------------------------------     23.1    22.6    23.2    23  22.9    23.2    23.2    23  22.5    22.5                     ------------------------------------------------------------------------------------------------------------     23.8    24.1    23.8    23.6    24.3    23.7    24  24.2    24.4    23.8    23.5                 ------------------------------------------------------------------------------------------------------------     24.7    24.9    25  25.4    25.1    25.3    25.3    25.2    25.2    24.7    24.6    24.5    24.5         ------------------------------------------------------------------------------------------------------------     25.6    26  25.6    25.6    25.6    25.7    25.6    26.3    26.3    26.2    26  25.9    26.2    26     ------------------------------------------------------------------------------------------------------------     27  27.2    27  27.2    26.7    27  26.6    27.3    26.8    26.6    26.9    27  26.5         ------------------------------------------------------------------------------------------------------------     28.2    27.7    28.1    27.9    27.6    27.7    28.1    28.2    27.5                         ------------------------------------------------------------------------------------------------------------     28.6    29.1    28.7    29.1                                             ------------------------------------------------------------------------------------------------------------     30.3    29.6    29.7    29.5                                             ------------------------------------------------------------------------------------------------------------     30.9    31  31  30.5                                             ------------------------------------------------------------------------------------------------------------     32.3    31.6     ------------------------------------------------------------------------------------------------------------         These are the groups:                - 17.5 -   18.5      - 18.5 -   19.5      - 19.5 -   20.5      - 20.5 -   21.5      - 21.5 -   22.5      - 22.5 -   23.5      - 23.5 -   24.5      - 24.5 -   25.5      - 25.5 -   26.5      - 26.5 -   27.5      - 27.5 -   28.5      - 28.5 -   29.5      - 29.5 -   30.5      - 30.5 -   31.5      - 31.5 -   32.5      I'm supposed to calculate the standard deviation and the expected value ... How is it different comparing to having only 15 values and not divided them into groups? If I only have 15 values, I know how to calculate the stand.deviation and the exp.value. I use the formulas that are well-known:               Exp.value: ............(1/n)*Sum(xi)..................where xi = all the 15 values <br>     Stand.deviation: ...(1/(n-1))*Sum( xi - u)^2...where xi = all the 15 values and u =  exp.value


title: Split plot with factorial design vs three way ANOVA

text: I have three site temperature types: cold, medium, and hot. This is a fixed effect - I am specifically am testing if temperature affects my response variable. I believe temperature type is my 'block'. I will have three sites for each temperature type. I think site is a random effect. Within each site, I will have a factorial design, where two other fixed factors, A and B, are crossed. I am interested in the interaction of all three fixed effect (tempxAxB). From reading my stats book, I think this is a split plot design, but my stats book does not have information about factorial design (the cross of A and B) within blocks. **Is what I just described a split plot design or a three way ANOVA?** And if it is a funky split plot, could someone point me to an ANOVA table (how to calculate sum of squares) with a similar set up? Here is a schematic of the setup: ![Set up](http://i.stack.imgur.com/E60hm.png)


title: Cumulative / Cumulative Plot (or "Visualizing a Lorenz Curve")

text: I don't know what such plots are called and thus I just gave this question a stupid title. Let's say I have an ordered dataset as follows               4253  4262  4270  4383  4394  4476  4635  ...      Each number corresponds to the amount of postings a certain user contributed to a website. I am empirically investigating the "participation inequality" phenomenon as defined here. In order to make it easy to grasp I would like to produce a plot which allows the reader to quickly deduce statements such as "10% of the users contribute 50% of the data". It should probably look similar to this admittedly pretty lousy paint sketch: ![enter image description here](http://i.stack.imgur.com/ysrzN.png) I have no clue how this is called thus I don't know where to look for. Also, if somebody had an implementation in `R`, that would be awesome.


title: Detecting outlier cash movements

text: If I'm watching a series of accounts for transactions going in and transactions going out, I want to notice unusually large or transactions for any particular account on any particular day. So if account A typically moves a few hundred dollars and one day moves five thousand dollars, that's a clear outlier. If account B typically moves a few million dollars in or out and one day moves 20 million dollars, that's a clear outlier. What I'd like to do is present a measure that should highlight outliers - I was thinking number of standard deviations versus a population of the rolling last 60 days, but I'm wondering if that's correct. I'm checking to see if it's a gaussian distribution, but are there better ways to hit what I'm looking for? I think this poses a different set of questions than Robust outlier detection in financial timeseries.


title: Gibbs sampler for conditionals that are exponential: Example from Casella & George paper

text: I am trying to work out Example 2 from Casella and George's paper "Explaining the Gibbs Sampler" in R. The example is:               f(x | y) = y*exp(-yx)     f(y | x) = x*exp(-xy)          where x, y are defined on [0,B]      The idea is to generate the distribution of f(x). My code (below) works, but what I do not understand is the role of the interval on which x and y are defined. In the paper this interval is [0, B]. In their example, the authors say that their simulation uses B = 5, but I am not clear where to incorporate this into my code...               X = rep(0, 500)     Y = rep(0, 500)          k = 15          for (i in 1:500) {         x = rep(1, k)         y = rep(1, k)              for (j in 2:k) {             x[j] = y[j-1]*rexp(1, y[j-1])             y[j] = x[j]*rexp(1, x[j])                }              X[i] = x[k]         Y[i] = y[k]      }     print(max(X))     print(max(Y))          hist(X, breaks=40, freq=F)


title: Time series as cross-sectional data

text: I have time series, for example, gdp and unemployment(`unemp`), freq= 4. What if I interpret it as cross-sectional data and do cross-sectional analysis instead of time series? My task is to test how unemployment affects gdp. Is it allowed to do that kind of analysis? Do the coefficients in the model `lm(gdp~unemp)` have an economic explanation?


title: Interpolation of spectra: uneven sampling to even sampling

text: I have a spectrum. Specifically, my data is relative intensity $[I_{\tilde{\nu}}]$ versus wavenumbers $[\tilde{\nu}]$. The wavenumbers are equally sampled so that ${d\tilde{\nu}} = c$, where $c$ is some constant. However, I need the data to be given in terms of wavelength $\lambda = \frac{1}{\tilde{\nu}}$. When I perform this transformation, the step size of wavelength $d \lambda$ decreases with increasing wavelength. However, I need the step size to remain constant for the entire spectrum. I am guessing this calls for some kind of interpolation or fitting. How can an even sampling of wavelengths $\lambda$ be accomplished?


title: Comparing two ordinal regression models

text: I am trying to find out if I can somehow assess if one model fits my data significantly better than another model? They are both ordinal regressions and we have introduced an additional interaction term in the second model. Now we want to know if the model has improved. Do I only look at the model fit indices? How do I know if the model has significantly improved?


title: Is MANOVA appropriate?

text: So I've been reading various online tutorials and I've not come across and example which is similar to what I'm trying to do. I've thought up this example because I was struggling to explain with words alone. It's really basic but hopefully illustrates the concept:               movie_eg=data.frame(score_children=c(1.1,2.5,3.9,1.5,4.3),                         score_adults=c(12.5,13.1,18.1,13.1,18.8),                         brad_hrs=c(0.1,0.2,0.3,0.1,0.6),                         angelina_hrs=c(0.2,0.6,0.8,0.9,0.1),                         cgi_budget=c(0.9,0.1,0.2,0.3,0.8),                         seat_cost=c(0.1,0.3,0.6,0.1,0.9))     rownames(movie_eg)=c("shrek","bourne","lockstock","scream","sharknado")      So in reality I have over 400 "movies" each with 2 continuous outcomes ("scores") and 4 predictive measurements ("brad","ang","budget","seats") which are also continuous. The two "scores" are moderately correlated across "movies". I'd like to see if there is a strong connection between any/all of the predictive values and the two outcomes - so here, what (if anything) influences both "scores"? Therefore, I believe that the "dependent variables" are the two outcomes and the "independent variables" are the four predictions. So I would do: `summary.aov(manova(cbind(brad_hrs,angelina_hrs,cgi_budget,seat_cost)~cbind(score_children,score_adults),data=movie_eg))` and then see all P values are greater than 0.05, so nothing had a strong effect on score? But that seat_cost was the best predictor (although still not significant)?               Response seat_cost :                                         Df  Sum Sq  Mean Sq F value  Pr(>F)       cbind(score_children, score_adults)  2 0.45417 0.227085  17.583 0.05381 .     Residuals                            2 0.02583 0.012915                         (Probably worth mentioning in reality the "scores" are measures of drug response and the predictions are measures of different cellular functions, so really nothing to do with movies!)