title: How can I sample from a multinomial distribution with a fixed expected number of remaining categories?

text: I have a population of entities associated with different categories, say               blue    50000     red     300     green   80     yellow  10     pink    6     orange  3     white   2      The distribution is known up the exact counts. It is very homogeneous, i.e. if ${p_1,...,p_k}$ denotes the probabilities for each category, then $max_i(p_i)\geq0.95$. Now I would like to choose a sample size, such that the expected number of categories $N_k$ present in the sample is $n_k$, i.e. the expected number of categories whose count is larger than zero is $n_k$. How, given a count vector as shown above and $n_k$, do I choose the sample size?


title: Can a percentage be estimated as reliably with a small group?

text: I realise that my question is quite simplistic for this site but I will ask anyway. I am wondering if percentage results become more unreliable the smaller the group is (I am not referring to a sample of the group rather questioning the whole group). If yes, is there a way to determine that? For example 5% of the 30 people that received a direct-sales call hang up the phone. Does that mean that 5% of 8 people called would be expected to hang up? Or is the percentage unreliable the smaller the number which is my main question? I hope that this was clear enough. Thank you very much


title: Discontinuity Regression

text: I'm writing a political economy Master's thesis regarding the effects of the Euro Crisis on the rise of separatism in Spain (Catalan independence for example). The idea is to show that the Crisis has acted as an exogenous shock, causing a sudden rise in separatism. Is a discontinuity regression the correct method to prove this and if so, how does it work, please? The variables I was thinking of using would be time v. separatism, unemployment v. separatism, GDP growth v. separatism, etc. Can it be done on SPSS? To note: I'm not very experienced with statistics.


title: Clustering of items based on their category belonging

text: I am trying to find a clustering algorithm, but I'm working with already classified items. Basically, items belongs to one or more category, which are already known. Categories are absolutely not linked each other. I need to create clusters where the number of distinct categories among all items is minimum. The same category could be splitted in several clusters. Cluster count and cluster size are known at the beginning. **Sample:**               | Item   | Category |     | ------ | -------- |     | Item 1 | C1, C2   |                             | Item 2 | C1       |  cluster count = 2          | Item 3 | C2       |  cluster size = 3     | Item 4 | C2       |                             | Item 5 | C3       |     | Item 6 | C1       |          Result:      - Cluster 1: [Item 1, Item 3, Item 4] // 2 categories (C1, C2)      - Cluster 2: [Item 2, Item 5, Item 6] // 2 categories (C1, C3)      I have already tried custom k-means implementation, but it seems to produce results far from optimal. Am I missing something? Does a well-known algorithm correspond to this kind of problem? * * * **EDIT** _let's explain a bit more_ I need to classify items based on the only attribute I know: their category belonging. I have used a custom distance function, which can be expressed as:               d(I1, I2) = number of category that I1 and I2 do not have in common      I have only tried to use random starting point, and it gives me different results, but far from optimal.


title: Random effects models / Integrate over the random effect

text: I am trying to do maximum likelihood estimation and trying to see if the problem can be formulated using a random effect model. Here is the problem description: > There are $100$ pairs $(N_i, D_i)$ where each $N_i$ is the number of trials > and $D_i$ the number of successes in the $i$-th Binomial experiment (with > $100$ trials). The probability $p$ of the binomial depends on the NormCDF > through a parameter $\rho$ (which needs to be estimated using MLE). The > probability $p$ is really a CONDITIONAL probability, conditional on $y$ = > some constant, where $y$ can range over $(0,1)$. Thus, I have a likelihood > function $L$ which is being integrated over a variable $y$, $0 < y < 1$. > That is, $L$ = int over $(0,1)$ of Binomial(n,d,p)dy,   or equivalently,   $L$ = int over $(0,1)$ of $\\{ k1 * p^d * (1-p)^ (n-d) dy \\},0<y<1$   where $p$ is a function of $y$, and $k1$,$n$,$d$ are constants, $k1 = {n \choose d}$. Also   `p(y) = NORMCDF( (gamma - {sqrt(rho)*INVNORMCDF(y)} ) / sqrt(1-rho) )`, where gamma is a constant (over all 100 trials). I am trying to maximize $L$ to obtain an estimate for $\rho$. I would like to know if I can think of the $y$ as a random effect, and also, if specifying $y$ as a random effect in `SAS proc nlmixed` implies that the likelihood $L$ is automatically integrated over $y$. Currently I am trying to maximize the likelihood in `SAS / IML` and the problem is the optimization for $\rho$ does not work (`IML` just keeps running, and I think it is because the likelihood function is flat enough that convergence takes too long). I am trying to understand what it means for a variable $y$ to be a random effect: does it mean that the likelihood is integrated over all possible values of $y$? If that is true, can I estimate $\rho$ using `PROC NLmixed` as follows: ` gamma = 0.5; PROC NLMIXED DATA=bindata MAXITER=200000 QPOINTS=10 GCONV=1e-15; PARMS rho = 0.5; * Define likelihood; w1 = INVNORMCDF(y); p = NORMCDF( (gamma - sqrt(rho)*w1 ) / sqrt(1-rho) ); L = (n CHOOSE d)*((1-p)**(n-d)) + log(p**d); MODEL d ~ GENERAL(L); RANDOM y~ NORMAL(0,1) SUBJECT=ID; RUN; ` In this way, I don't need to integrate the binomial density over $y$ in $(0,1)$ myself in `IML` using a `QUAD` function as I am doing now.


title: What's the probability of an event happening exactly once after two independent trials?

text: **The problem:** > If event $A$ happens with probability $p$, what is the probability $P[K=1]$ > that, after two independent trials, $A$ happened exactly once? Is the answer:   * $P[K=1] = p(1-p)$, since we need $A$ to happen and then we need $A$ _not_ to happen, or   * $P[K=1] = 2p(1-p) = p(1-p) + (1-p)p$, since there are two cases $A\overline{A}$ and $\overline{A}A$ for which $A$ happen exactly once?


title: How to calculate margin of error for a binomial quality control experiment where only successes are observed (including FPC)?

text: Wikipedia's _Margin of Error_ entry says that > a random sample of size 400 will give a margin of error, at a 95% confidence > level, of 0.98/20 or 0.049 - just under 5% given an infinite population size. This means that if I polled 400 US citizens (randomly selected) and asked them "Obama or Romney?", the resulting proportion would be accurate with a 95% confidence level, to a margin of error below 5%. However, can I use this same calculation in testing that a software program will be able to deal with roughly all inputs? For example: I have an infinite population of users, I need to be sure (95% confident, given 5% error) that my software will be able to come up with a nickname for all of them based on their name and a simple algorithm. If I randomly select 400 users, and the nickname algorithm works perfectly for all 400 users, can I assume (with 95% confidence, given 5% error) that my algorithm holds for the entire population? Is this the incorrect way to calculate margin of error for this type of problem?


title: Tuning paramaters SVM, DT, k-NN, NN

text: I'm trying to compare the predictive strenght of four different algorithms:   * support vector machines   * k-NN   * decision trees    * neural networks I've got a few questions concerning the parameter tuning:   1. Some papers like [1]state that k-nearest neighbour does not need parameter tuning, don't you need to decide what k will be?   2. If you want to tune the parameters of SVM (gridsearch) and NN (different hidden neurons), on what data set do you do that if you plan on using tenfold cross-validation? If you're using your data set with tenfold, isn't there any data left to use to validate your parameters?  [1] Xiao, W., Qian, Z., & Fei, Q. (2006). A comparative study of data mining methods in consumer loans credit scoring management. Journal of Systems Science and Systems Engineering, 419-435.


title: What is the point of non-informative priors?

text: Why even have non-informative priors? They don't provide information about $\theta$. So why use them? Why not only use informative priors? For example, suppose $ \theta \in [0,1]$. Then $\theta \sim \mathcal{U}(0,1)$ is a non- informative prior for $\theta$.


title: Can't find p-values in the output from lmer() in the lm4 package in R

text: I'm a new user of mixed models and through the material I've been reading there are always probability values (p>t) or (p>z) that estimate the importance of a level of a factor in the model. However, when using the `lmer()` function in R, which supposedly gives you those probabilities, I simply don't find them. Here is the output:               Linear mixed model fit by REML      Formula: Temp ~ depth + (1 | locality)         Data: qminmatrix         AIC   BIC logLik deviance REMLdev      561.3 581.3 -273.7    551.5   547.3     Random effects:      Groups   Name        Variance Std.Dev.      locality (Intercept) 4.7998   2.1909        Residual             4.0433   2.0108       Number of obs: 128, groups: locality, 4          Fixed effects:                 Estimate Std. Error t value     (Intercept)  22.0103     1.1500  19.140     depth1        1.9564     0.6832   2.864     depth10       2.6624     0.5756   4.625     depth5        3.0209     0.4932   6.125     depthWS      -2.2585     0.5444  -4.149          Correlation of Fixed Effects:             (Intr) depth1 dpth10 depth5     depth1  -0.157                          depth10 -0.175  0.189                   depth5  -0.213  0.313  0.458            depthWS -0.191  0.334  0.373  0.441


title: Durbin-Watson  test for residuals of SARIMA(1,1,1) in R

text: I'm modelling a time series data using ARIMA. Now, I'm trying to test for the serial correlation of my model SARIMA(1,1,1) using the durbin watson test. My problem is that I don't know what linear model I would put on the formula of the `dwtest` function in R. Here's the usage of the function,               dwtest(formula, order.by = NULL, alternative = c("greater", "two.sided", "less"),            iterations = 15, exact = NULL, tol = 1e-10, data = list())      Here's my code below, Data: http://iitstat.weebly.com/uploads/7/3/4/0/7340846/chickenprod.rdata To download the data just right click the link and click "Save Link As..."               library(forecast)     library(lmtest)     ChickenProd <- ts(ChickenProd, start = 1980, frequency = 4)     SARIMA111 <- Arima(ChickenProd, seasonal = list(order = c(1,1,1), period = 4))      The residuals of my model SARIMA111 is obtain by               SARIMA111[["residuals"]]      Now, I want to test the serial correlation of it using the Durbin-Watson test, but I don't know what linear model I would use in the `formula` argument of `dwtest` function in R. Is it the SARIMA(1,1,1) model? If so, how will I extract the coefficients of the SARIMA(1,1,1) model, and make a linear model formula in R? Thank you in Advance!


title: Probability of winning a drawing/lottery

text: I can't seem to get my head around this: A random drawing will be held for which there are 900 tickets sold, for which there will be one winning ticket drawn. If you purchased 25 tickets, what is the probability that you will win. I am pretty sure this follows a binomial distribution, but I don't know how to adopt the models to solve the problem. Any help / directin would be appreciated. I am looking more to undertsand the problem domain than to get an answer. Thank you in advance,


title: Is it possible for a bayesian model to "forget" about some points?

text: Say you are doing some online bayesian inference over observations. You want to infer $\mu, \sigma$ for $X$, with a model $X \sim \mathcal{N}(\mu, \sigma)$. Now everytime you observe a $X_i$, you update your prior with your new posterior and get a new estimate $\mu_i, \sigma_i$. That becomes the prior for your next observation, etc etc. Now after $n$ observations, your estimates $\mu_n, \sigma_n$ come from applying "recursively" bayes to your priors $\mu_0, \sigma_0$, where : $$P(\mu_n, \sigma_n) \propto P(\mu_{n-1}, \sigma_{n-1}) \times P(X_n \sim \mathcal{N}(\mu_{n-1},\sigma_{n-1}) | \mu_{n-1},\sigma_{n-1})$$ Now is it possible to make your bayesian model "forget" about old points without having to recalculate the whole chain ? In other words, if you wanted to apply this model to a stream of data (because you have lots of points), can you make it forget about old points ? I hope I'm clear enough with my question, I'll edit as comments come if needed. Thanks a lot


title: Box Cox Transforms for regression

text: I'm trying to fit a linear model on some data with just one predictor (say (x,y)). The data is such that for small values of x, the y values give a tight fit to a straight line, however as x values increase, the y values become more volatile. Here is an example of such data (R code)               y = c(3.2,3.4,3.5,3.8,4.2,5.5,4.5,6.8,7.4,5.9)     x = seq(1,10,1)      I'm curious to know if there exists any power transform (Box cox perhaps?) that lets me get a better fit for the data than simply doing a linear fit as shown below.               fit = lm(y ~ x)      Thanks


title: Learning Curves Example

text: I'm trying to find a full example of how to plot learning curves. I watched Andrew Ng's ML class on Coursera and he mentions using learning curves to diagnose variance-bias issues. My notes show the top line as the test set or cv error slopping downward as the number of training examples increases. This makes sense to me. My problem is on the lower line of the plot which is training error. I'm not sure what this means and why the error would increase with more training examples. Should the training error be multiplied by -1? Or is there something else I need to do? Any working examples would be appreciated as would an R package that helps with this.


title: Spatial clustering with the constraint that all clusters have equal count

text: I wish to perform a spatial clustering of scattered data that represents geographic locations of individuals in an urban area. Hierarchical clustering seems to work well, and I have successfully done this for the 100,000 data points in the set. However, as an additional constraint/objective, we wish that the clusters have equal count. That is, the city center would be a small area spatially, and the suburban areas would be large. Is there a ready-made solution for this problem?


title: Literature on robustness of regression assumptions

text: In my OLS regression not all assumptions are perfectly met, but I read that due to a large sample size there is a certain robustness to assumptions (my sample is 2500 people). E.g. the DV isn't perfectly normal distributed, with significant, but small skewness and kurtosis (if you test for it, e.g. using the 'gvlma' package in R). I am now looking for literature to cite when I talk about my regression, e.g. "Inspection of fitted values vs mean plots revealed minimal deviation from homoscedasticity, but within range of robustness of regression (xxx, 1990)."


title: Feature construction and normalization in machine learning

text: Lets say I want to create a Logistic Classifier for a movie M. My features would be something like age of the person, gender, occupation, location. So training set would be something like:   * Age Gender Occupation Location Like(1)/Dislike(0)   * 23 M Software US 1   * 24 F Doctor UK 0 and so on.... Now my question is how should I scale and represent my features. One way I thought: Divide age as age groups, so 18-25, 25-35, 35-above, Gender as M,F, Location as US, UK, Others. Now create a binary feature for all these values, hence age will be having 3 binary features each corresponding to an age group and so on. So, a 28 years Male from US would be represented as 010 10 100 (010-> Age Group 25-35, 10 -> Male, 100 -> US) What could be the best way to represent features here ? Also, I noticed in some e.gs. of sklearn that all the features have been scaled/normalized in some way, e.g. Gender is represented by two values, 0.0045 and -.0.0045 for Male and female. I don't have any clue on how to do scaling/mormalization like this ?


title: What are good techniques for modeling small datasets?

text: I’m working on a classification problem. However, my training dataset is very small (just 800 items in training dataset) and each data item contains a small number of features (just 5 features). Firstly, I used Logistic Regression to create a model for this dataset. Unfortunately, prediction accuracy of my model was very bad. Next, I used Neural Network model, but could not see any progress. I suspect, number of training data items and number of features in each item are not enough for training Logistic Regressing and Neural Network. So my question is what are the good techniques for modeling small datasets?


title: Cross-correlation in R – dealing with missing values

text: I have two time series. One is an environmental variable ( _n_ = 108) organized by year and month. The other is a biological variable, also organized by year and month, but I have no data for some months ( _n_ = 97). I did a cross-correlation in R between these 2 times series, and used the `na.exclude` function for the biological variable to account for the missing values. I then fit a lagged regression, again accounting for the missing values. However, I want to know:   1. Could I do a cross-correlation by deleting these missing values (no data months) for the biological variable, which will leave me with unequal lengths? R still does a cross-correlation, but I am unsure if the results are valid.   2. Could I delete the corresponding months in the env. variable to have equal length time series, or is accounting for missing values the correct approach? Thanks.


title: How to regress a time series of proportions?

text: Every month, an organization surveys some of its customers (the total number of customers is also known). The sampled customers answer a survey with a dozen or so questions; sometimes, customers don't answer every question. A simple average of the _proportion_ of answering customers who gave the top category on Question 1 and the _proportion_ who gave the highest answer to Question 2 is calculated. (The top category for some questions is "Excellent," for other questions it is "Always.") This average creates a score that is reported at the end of each month. What command(s) in SPSS could compute the probability/magnitude of improvement in that score over a given historical window? Also, is it possible to smooth this score? Predict future scores? I've already investigated an ARIMA, but couldn't figure out how to incorporate the boundaries of the dependent variable (0,1) into the command.


title: Getting started with interior point methods

text: I want to learn about the interior point methods. I started with the wiki http://en.wikipedia.org/wiki/Interior_point_method. However, I am having difficulty grasping the idea behind it. Can anyone provide me good links, video tutorials or pointers that can help me understand it?


title: VAR vs STAR for space-time autoregression in Python

text: I want to use autoregressive model to build a predictor for some sets of spatio-temporal data. For example, I have historical traffic data (speeds at various segments of freeways). similarly, I have historical weather data for different cities. Needless to mention, there is significant spatial correlation between nearby sites. I can think of 2 ways of autoregressing (assume that the order of autoregression is small (say <= 3)): a. Use Vector autoregressive (VAR) model: where each site depends on values at all other sites in the previous time instants. b. Use Space-Time AR (STAR) model: where each site depends on values at known nearby sites in previous time instants. Which of the two models is more suitable? Python's statsmodels library has an implementation for VAR, but not for STAR. Since I am using Python for my work, I am tempted to use VAR. Could someone explain the advantages of one over the other? I am perhaps slightly better than an average coder. If I want to implement STAR model in Python, how difficult would it be? Is there a text that can help me implement it?


title: When is the PMI value good or bad?

text: Pointwise mutual information is calculated by this formula $pmi(x;y) = log(p(x,y)/p(x)p(y))$ , my question now is, When is this pmi good and when is it bad. I know if the value is low it is bad, but when is it possible to assume there is some decent associations between x and y ? Or do I have to determine this threshold myself (experimentally) ?


title: Ordinary least squares - change response and explanatory variable

text: In a simple regression model using OLS, why is it not, at least in general, possible to move the $x_i$'s to the left-hand side of the model and the $y_i$'s to the right hand side when I want to switch the explanatory and response variable? Intuitively, it seems fairly clear that if ${\hat \beta ^{OLS}}\left( {Y,X} \right) = \mathop {\arg \min }\limits_{\beta \in \mathbb{R}} \sum\limits_{i = 1}^n {{{\left( {{y_i} - x_i^T\beta } \right)}^2}} $ and ${\hat \gamma ^{OLS}}\left( {X,Y} \right) = \mathop {\arg \min }\limits_{\gamma \in \mathbb{R}} \sum\limits_{i = 1}^n {{{\left( {{x_i} - y_i^T\gamma } \right)}^2}} $, the simple "switching of sides" would not work since I'm minimizing two different sums. But is there any more formal explanation why it cannot be done?


title: Looking for a specific kind of plot/graph

text: I'm looking for info about a particular kind of graphing technique. Basically, you have a list of elements, and data indicating how closely related the elements are to each other. On the plot, elements that are closely related are plotted close to each other, and ones that are less related are drawn farther away. I'm trying to find out what this kind of graph is called and if there are any related good algorithms for it. Sorry if this is a bad question, but this has been driving me nuts and it's impossible to search for. * * * Update: So it was a _network graph_ I was looking for after all; I found this neat freeware tool called Gehpi that pretty much serves my needs. Thanks for the help!


title: Question on unexpected regression results

text: Sorry this is a bit long, but the explanation of what I want to do needs to be clear to avoid issues such as this quote..."It is impossible to speak in such a way that you cannot be misunderstood." Karl Popper. I am running linear regression models, but I am getting expected results. I wonder what else I might try to derive an estimated value of bat echolocation parameters based on forearm measurements. It is known that the size of the bat is negatively related to the characteristic frequency (Fc) of their echolocation calls (decades of my field work) . So in general larger guys have lower frequency calls and smaller guys have higher frequency calls. I have run the regressions based on the FA (valid forearm measurements) and the known and valid Fc ranges for a dozen species or so and using the lm models to "predict" Fc values for a few species that have FA values but have not yet been recorded. Hence there are no valid echolocation call parameters. R Code used is below discussion. I have valid ranges for the known species FA (forearm measurements) and Fc(minimum) and Fc (maximum). So I do two separate runs with the data using the lm model one with FA~Fcmin and one FA~Fcmax. The goal is to provide the predicted (estimated) values for the species with known FA values but w/o verified Fc value ranges. My concern is that the predicted values returned are much lower than the true values for the verified species. Therefore I am not confident the predicted values for those w/o verified Fc ranges are useful. One very helpful person looked at one simple data set I sent and showed that the statistical differences between the true values and predicted were not significant. However Krebs' admonishment to students eons ago "Do not confuse statistical significance with ecological significance" is true here. The values of the predicted ranges are far lower than reality so the few species that do not have field recorded Fc values are suspect. These differences in predicted values from a true range will will make a difference for potentially IDing the unknown calls. A difference of 10kHz Fc generally suggests a different species, albeit some are much closer and may only have a 5 kHz difference. I am looking at acoustic data sets of calls from South America and there are many "sonospecies." These are clearly separate species based on echolocation call parameters that have yet had "faces & voices" matched. We know that call parameters are diagnostic for families and genera even when the species is unknown. It is then the Fc values that assist in identifying the species within a cluster of calls from the same genus. Sample of R code used: Bats <\- dget('C:/=Bat data working/Acoustic Parameters/_Working/=Vespertilionidae/Bats.robj') model.lm <\- lm(formula=Fc ~ as.factor(FA),data=Bats,na.action=na.omit) > Anova(model.lm,type='II') Error in solve.default(L % _% V %_ % t(L)) : > system is computationally singular: reciprocal condition number = 0 > summary(model.lm) Call: lm(formula = Fc ~ as.factor(FA), data = Bats, na.action = na.omit) Residuals: ALL 5 residuals are 0: no residual degrees of freedom! Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 53.3 NA NA NA as.factor(FA)34.4 -4.6 NA NA NA as.factor(FA)35.4 2.3 NA NA NA as.factor(FA)35.5 9.0 NA NA NA as.factor(FA)40.5 -7.3 NA NA NA Residual standard error: NaN on 0 degrees of freedom (2 observations deleted due to missingness) Multiple R-squared: 1, Adjusted R-squared: NaN F-statistic: NaN on 4 and 0 DF, p-value: NA > tmp<-predict(model.lm) Bats[names(tmp),"predicted"]<-tmp rm('tmp') > rm('model.lm') > > model.lm <\- lm(formula= Fc~ FA,data=Bats,na.action=na.omit) > > Anova(model.lm,type='II') > > summary(model.lm) > > tmp<-predict(model.lm,Bats) > > Bats[names(tmp),"Predicted.Fc"]<-tmp > > rm('tmp') rm('model.lm') With the results it can be seen that the predicted Fc values on right are not close to the true Fc values on left and then make me hesitant to accept the 2 with NA predicted values. FYI Species are simple 6 letter coded for genus and species. Species FA Fcmin Fcmax FcMinpredic FcMaxpredic Myoalb 35.3 45.7 48.7 51.73 55.26 Myoata 37 NA NA 49.52 52.59 Myokea 33.7 57.8 61.3 53.80 57.77 Myonig 34.5 51.6 55.7 52.77 56.52 Myooxy 40.5 45.7 47.6 44.98 47.09 Myorip 36 53.3 57.5 50.82 54.16 Myosim 38 NA NA 48.23 51.02 Perhaps simple linear regression is not the method to use? Thanks for any additional suggestions. Bruce


title: Factors contained in the unobserved effect?

text: I try to study differences in health care cost between counties over a several years period, using a fixed effects model. I use a cost measure as dependent variable and have some demographic, geographical, social and size-related factors as explanatory variables. What are some of the factors you think could be contained in the unobserved effect, $a_i$?


title: Bayes Rule. Given I flipped 10 heads in a row, what is likelihood I own a double headed coin?

text: $$\Pr(A | B) = \frac{\Pr(B | A) \Pr(A)}{\Pr(B)}$$ So my $\Pr(A)$ is probability that a coin is double headed. For the sake of argument let's say 1 in 10,000 coins are double headed. My $\Pr(B)$ is the probability of flipping 10 heads, which is 1 in $2^{10}$. So naively I thought $\Pr(B|A)$ well the probability that I flip 10 heads given that I have a double headed coin is 1? That leaves the probability that I have a double headed coin at $\frac{10^{-4}}{2^{-10}}$. However this means that the probability will be > 1 for larger flips. The step in which I have made a mistake is evaluating $\Pr(B|A)$. Can anyone help?


title: applying random forest in place of decision tree method in R

text: I am relatively new to random forest technique. I applied decision tree successfully to a data predicting 1/0 outcome. The predictors are mixture of continuous and categorical. I was trying similar to tree method: train_tree <\- randomForest(Response_variable ~ predictor1 + predictor2+...predictor10 ,trainData,ntree=50) But its throwing me error. Can you please help me?


title: Transition rates in continuous time markov chain

text: > A house has 2 rooms of similar sizes with identical air conditioners > equipped with thermostats which turn on and off as needed to maintain the > temperature in each room to a desired level of 22 degrees. Suppose that a > thermostat remains on or off for exponential amounts of time with means > $1/\mu$ and $1/\lambda$, respectively, independently of other thermostats. > Consider the Markov process $\\{X(t), t \ge 0\\}$ whose state space is the > number of active air conditioners. Write down the matrix of transition > rates. I'm not sure how to exactly approach this type of question. My working is as follows but if someone could clarify my confusion that would be good. Working: So clearly there are 3 states, 0 for no air conditioners are on, 1 for one air conditioner is on (active), and 2 for two air conditioners are on (active). Now to work out $q_{01}$, i.e., the transition rate from state 0 to 1, assume currently no air conditioners are on. Consider the two independent poisson processes both with parameter $\lambda$, where the interarrival time is the duration of "off", then merging these two processes gives a poisson process with parameter $2\lambda$, so the transition rate from state 0 to state 1 is $2\lambda$. Now what about the transition rate from state 0 to 2? I am told that it's 0, but why? Isn't it possible for both air conditioners to both go from "off" to "on"? What is the argument that $q_{02} = 0$?


title: glmm unbalanced nested random factors

text: imagine you have the following data structure                                                    data                                            |                                            |                        env1    env2    env3    env4      env5                         |       |       |        |         |         replicates   1,2,3    4,5,6    7,8       9         10      At each replicate you apply your nominal treatment A & B and measure your response for each replicate at under seqeuntial A and B treatments. So each replicate has a data point for A and B treatments. Using `lme4` a GLMM of this data might look like:               glmer(response ~ treatment + (1|env/replicate), family=poisson, data=data)      I would like to know what are the effects (repercussions for my model) of having only one replicate per environment in some of the cases?


title: Sample size calculation for comparison of 3 proportions

text: How do I calculate sample size for estimating difference in proportions between 3 groups? Most of the free sample size calculators have options for only 2 groups.


title: Question about Dynkin Lehmann Scheffe Theorem

text: I'm self-studying for an examination, and I would like to understand how to use the Dynkin Lehmann Scheffe theorem for an applied question. I am using Bickel and Doksum's "Mathematical Statistics" (2007 edition), and there is only a sentence describing the process: Let P = {$P_{\theta}$ : $\theta$ $\in$ $\Theta$} where $P_{\theta}$ is discrete concentrated on X = {$x_1$,$x_2$,....}. Let $p(x,\theta)$ $\equiv$ $P_{\theta}$$[X = x]$ $\equiv$ $L_x(\theta)$ > 0 on X. It can be shown that $\frac{L_x(.)}{L_x(\theta_0)}$ is minimally sufficient. I found a question with the solution, but I would like to see how the process works. Here is the question: ![enter image description here](http://i.stack.imgur.com/zc9Ys.jpg) Here is the solution: ![enter image description here](http://i.stack.imgur.com/24AKX.jpg) I apologize about not having shown any work, but there is a dearth of information available on this theorem, so I just wanted to see how the process for this type of problem works.


title: Population specific attributes to pairwise distances: How to choose a distance?

text: I've collected a lot of data for different populations and I need to convert them into a pairwise distance measure. There are several distances such as , Bray-Curtis and Manhattan distances. My question is How to choose a distance? Is there a book or paper has a clear criteria on how to choose? thanks


title: Access tables created in SAS Enterprise Guide Client into SAS Enterprise Miner Client?

text: How can I access tables created in SAS Enterprise Guide Client into SAS Enterprise Miner Client?


title: Cross Validation on Training and Testing

text: I have 500 instances of student dataset and 10 attribute. If I am dividing 170 for training, 330 for testing. These ration is correct for crossvalidation?


title: Weak stationarity and ARMA-ARCH/GARCH models?

text: I am slightly irritated about weak stationarity in connection to ARCH/GARCH models. I do not know the answer and I am not sure about it: The basic question is: > Do we have to test weak stationarity before applying an ARMA-GARCH model? Further on it can be said: > ADF and others test the mean equation, but this is not for the volatility > equation, so what test do we have to use for the autocovariance- > stationarity? My thoughts: Standard ARMA models assume the unconditional mean and unconditional variance to be constant. For ARMA-GARCH models this is also the case: The unconditional mean and unconditional variance need to be constant, whereas in case of the ARMA-GARCH models the conditional variance does not need to be constant. It is correct that for the mean equation we may have to think about using a trend-stationary or difference-stationary model. But this is only concerning the mean equation, yes. Conditional variance can be tested by testing for ARCH effects (Box-Ljung, Lagrange Multiplier). So for ARMA-GARCH models we still need weak stationarity, since the unconditional mean and unconditional variance need to be constant. So I am not sure, but we have to test for weak stationarity before applying an ARMA-GARCH model? And especially with financial returns, do we also have to test for it? And which test do we use (and which command is implemented in R, so what command can you suggest?) I know that if the unconditional variance is nonstationary (not finite and then it is also not constant) an integrated GARCH may be appropriate. But just because it is not constant I cannot say I use an integrated GARCH model? I also know that for ARMA-GARCH processes all the "characteristical" roots lie outside the unit circle. So in case of a ARMA-GARCH(1,1) $\alpha_1+\beta_1<1$ is necessary. But I only see this after estimation? This is not a test for covariance stationarity? **EDIT** : It basically pins down to (see the comments): How can I test the unconditional variance to be constant? I mean in order to apply a GARCH model I have to make sure that I have constant unconditional mean (ADF test and so) and I have to test for constant unconditional variance (how?). I know that I have to further make sure that I have nonconstant conditional variance for GARCH processes, otherwise having a constant conditional variance ARMA is sufficient (test for ARCH effects). EDIT 2: There is a Wavelet Spectrum test in the locits package, what about this test?


title: scikit-learn score metric

text: I am using `scikit-learn` in Python and they define a quantity called _score_. It's defined in the middle of the documentation page. Reproduced here: > Returns the coefficient of determination R^2 of the prediction. The > coefficient R^2 is defined as (1 - u/v), where u is the regression sum of > squares ((y_true - y_pred) ** 2).sum() and v is the residual sum of squares > ((y_true - y_true.mean()) ** 2).sum(). Best possible score is 1.0, lower > values are worse. A few questions regarding this:   1. What's the intuition behind this metric?   2. What is considered a good score? What is considered bad?


title: What to do when the means of two samples are significantly different but the difference seems too small to matter

text: I have two samples ($n \approx 70$ in both cases). The means differ by about twice the pooled std. dev. The resulting $T$ value is approximately 10. Whilst it's great to know that I have conclusively shown that the means are not the same, this seems to me to be driven by the large n. Looking at histograms of the data I certainly do not feel that such as small p-value is really representative of the data and to be honest don't really feel comfortable quoting it. I'm probably asking the wrong question. What I'm thinking is: ok, the means are different but does that really matter as the distributions share a significant overlap? Is this where Bayesian testing is useful? If so where is a good place to start, a bit of googling hasn't yielded anything useful but I may not by asking the right question. If this is the wrong thing does anyone have any suggestions? Or is this simply a point for discussion as opposed to quantitative analysis?


title: How defensible is it to choose $\lambda$ in a LASSO model so that it yields the number of nonzero predictors one desires?

text: When I determine my lambda through cross-validation, all coefficients become zero. But I have some hints from the literature that some of the predictors should definitely affect the outcome. Is it rubbish to arbitrarily choose lambda so that there is just as much sparsity as one desires? I want to select the top 10 or so predictors out of 135 for a cox model and effect sizes unfortunately are small.


title: Decision tree unclassified data

text: not sure if this is the place for this, please let me know if not! I'm trying to complete a piece of homework and one of the questions is : Given this tree derived using the ID3 algorithm:               Gender = Male         TV  = Yes             Income = 40 - 50                 Age = 30 – 40 -> Take up = Yes                 Age = 40 – 50 -> Take up = No             Income = 30 – 40 -> Take up = Yes             Income = 50 - 60 -> Take up = Yes             Income = 20 – 30 -> Take up = Yes         TV = No -> Take up = no      and this unseen data case :               Gender -> Male     TV Promotion -> Yes     DVD Promotion -> No     Income -> 40 – 50     Age - > 20 – 30      Take up -> ?      What should be the take up value of this case. I have a feeling it's yes because it falls under the 30 - 40 age bracket but I'm trying to find out if there's actually a factual way of working it out but I'm not sure what I'm searching for. I seem to have searched the first 10 odd google pages of "unseen data on ID3 decision tree" but I'm not getting anything remotely useful. Relevant links would be appreciated! Thanks.


title: An event-log data-set for evaluating a filtering system?

text: For my research, I am working on an event filtering/matching system that considers each event to be a set of key-value pairs (e.g. "action=logon", "user=john.doe"). I would like to evaluate its validity, but I have been unable to find a public data-set with my requirements:   * It has to be a time series of samples with a consistent internal structure.   * It should be sizable - at least several millions of records/events.   * It has to be pre-classified to a number of categories/use-cases. A mere classification into e.g. normal/abnormal is not enough. Alternatively, I could also use a set with a domain-specific tool that can perform this classification accurately. There are several data sets (e.g. web server logs) around, but I have yet to find one that satisfies this requirement; I need to be able both to train my system and evaluate its performance using _an accurate golden standard_.   * I would like the categories to be conceptual equals e.g. various bank transaction types. I am currently using the KDD99 data-set, but it is too focused on anomaly detection; while it does have separate anomaly classes, all "normal" events are lumped into the same category... Is there a data-set somewhere that I could use? Preferably something that has already been peer-reviewed to a degree? As an alternative, is there a standard tool or other way to _generate_ such a data-set? I have considered setting up e.g. a Web application and hammering at it with a load tester, or even using a generic simulation tool, but I believe that none of these approaches would really stand to peer review, which would compromise the evaluation itself.


title: Testing for statistical significance in survey results for different countries

text: I have just run the same survey on websites in different countries and I want to be able to tell whether there are significant differences between countries. For example if 20% of Germans donate to charity, 30% of British and 40% of Dutch, can you confidently say that the Dutch donate twice as much to charity as the Germans. Sorry if this is a stupid question, I have not studied statistics and have tried reading about t-testing and anova testing but it was all very confusing :( Thanks.


title: Is this a bug in auto arima or am I doing something wrong?

text: I must be doing something very wrong here, as auto.arima in R is completely dying, but I can't see what it is. I have the latest version of forecast and R and I think this happens on both Windows and Unix. It works for some/most time series I tried (equities) but fails for others. It seems to fail more often when I diff the high/low with the previous day's close as opposed to just diffing the closes as below. Is this a bug or am I somehow giving arima bad data? (and causing it to die with a horrible error message) I tried searching for this error but didn't come up with much. Thanks a lot.               library(tseries)     library(forecast)     dwa <- get.hist.quote(instrument="DWA", start="2010-01-01", end="2013-10-31")     logreturns <- diff(log(dwa$Close))*100      fit <- auto.arima(logreturns, trace=TRUE)      Output:               This is forecast 4.8           trying URL 'http://chart.yahoo.com/table.csv?s=DWA&a=0&b=01&c=2010&d=9&e=31&f=2013&g=d&q=q&y=0&z=DWA&x=.csv'     Content type 'text/csv' length unknown     opened URL     .......... .......... .......... .......... ....     downloaded 44 Kb          time series starts 2010-01-04     time series ends   2013-10-07           ARIMA(2,1,2) with drift         : 1e+20 *      ARIMA(0,1,0) with drift         : 1e+20 *      ARIMA(1,1,0) with drift         : 1e+20 *      ARIMA(0,1,1) with drift         : 1e+20 *      ARIMA(1,1,2) with drift         : 1e+20 *      ARIMA(3,1,2) with drift         : 1e+20 *      ARIMA(2,1,1) with drift         : 1e+20 *      ARIMA(2,1,3) with drift         : 1e+20 *      ARIMA(1,1,1) with drift         : 1e+20 *      ARIMA(3,1,3) with drift         : 1e+20 *      ARIMA(2,1,2)          : 1e+20 *Error in if (diffs == 1 & constant) { : argument is of length zero     Calls: auto.arima -> myarima     In addition: Warning messages:     1: In if (is.constant(x)) { :       the condition has length > 1 and only the first element will be used     2: In if (is.constant(x)) return(d) :       the condition has length > 1 and only the first element will be used     3: In if (is.constant(dx)) { :       the condition has length > 1 and only the first element will be used     Execution halted


title: Graphical Probabilistic Representation of Naive Bayes

text: Given the Naive Bayes graphical representation below, I want to calculate $P(X|Y_1,Y_2)$. Are the calculations below correct? ![enter image description here](http://i.stack.imgur.com/xWOqG.png) The factored joint distribution regarding the system is: $$P(X,Y_1,Y_2)=P(X) \cdot P(Y_1|X) \cdot P(Y_2|X)$$ $$\frac{P(X,Y_1,Y_2)}{P(X)}=P(Y_1|X) \cdot P(Y_2|X)$$ $$P(Y_1,Y_2|X)= P(Y_1|X) \cdot P(Y_2|X)$$ Using Bayes rule to flip the left hand side: $$\frac{P(X|Y_1,Y_2) \cdot P(Y_1,Y_2)}{P(X)}= P(Y_1|X) \cdot P(Y_2|X)$$ $$P(X|Y_1,Y_2) = \frac{P(Y_1|X) \cdot P(Y_2|X) \cdot P(X)}{P(Y_1,Y_2)}$$ $$P(X|Y_1,Y_2) = \frac{P(Y_1|X) \cdot P(Y_2|X) \cdot P(X)}{P(Y_1) \cdot P(Y_2)}$$


title: Bootstrapping - do I need to remove outliers first?

text: We've run a split test of a new product feature and want to measure if the uplift on revenue is significant. Our observations are definitely not normally distributed (most of our users don't spend, and within those that do, it is heavily skewed towards lots of small spenders and a few very big spenders). We've decided on using bootstrapping to compare the means, to get round the issue of the data not being normally distributed (side-question: is this a legitimate use of bootstrapping?) My question is, do I need to trim outliers from the data set (e.g. the few very big spenders) before I run the bootstrapping, or does that not matter?


title: Holt's Linear and Holt-Winters in R

text: With the below code, I have run Holt's linear and Holt-Winters forecasts using Excel / Solver. I wanted to replicate this using R (Excel can be a pain) but I am getting the below error with `hw()`. It is my first time running this with R so I am trying to go by the docs. I assume running the function will pick the initialization and also optimize the parameters? Or am I way off? Either way, I thought `hw()` would work, but clearly I am doing something wrong here.                  x<-c(12638.8,11583.3,13024.1,12594.2,13068.1,12765.4,13125.6,13316.3,13054.7,13879.3,14436.6,17861,13923.8,12418.4,13854.4,13558.8,14335,13639.9,14132.2,14457.2,14045.4,14941.2,15556,18984.2,14825.5,13138.3,14473.6,14507.6,15055.9,14384.2,15199.7,15169.7,15223.5,16412.9,16662.9,20638.9,16154.8,14696,15793.8,15804.2,15895,15868.6,16435.6,15945.9,16381.4,17084.2,17261,21652.6,16330.6,14840.5,16419.8,16104.2,16326.2,16453.4,16750.3,16678,16816.4,17431.7,17920.6,22491.3,17016.6,15534.8,17339.2,17065.6,17313.3,17448.2,17617.5,17814.2,17747.6,18592.4,19223.3,23636.1)         x.ts <- ts(x, start=2001, frequency=12)          plot(x.ts)          ####Holt's Linear     holt(x, h=12, damped=FALSE, level=c(80,95), fan=FALSE,          initial=c("optimal"), exponential=FALSE,         alpha=NULL, beta=NULL)          #Holt-Winters     hw(x, h=2*frequency(x), seasonal="multiplicative", damped=FALSE,          level=c(80,95), fan=FALSE, initial=c("optimal"),          exponential=TRUE, alpha=NULL, beta=NULL, gamma=NULL)          # hw gives the following error     #Error in ets(x, "AAA", alpha = alpha, beta = beta, gamma = gamma, damped = damped,  : Nonseasonal data


title: $E[e^{cX}]$ where $c < 0$ and $X$ is lognormally distributed

text: I am trying to calculate the expectation $$E[e^{cX}]$$ for arbitrary $c<0$ (for $c>0$ the expectation is infinite) if $X$ is lognormally distributed, i.e. $\log(X) \sim N(\mu, \sigma)$. My idea was to write the expectation as an integral, but I did not see how to proceed: $$E[e^{cX}] = \frac{1}{\sqrt{2\sigma\pi}}\int_0^\infty \frac{1}{x}\exp\left(cx - \frac{(\log x - \mu)^2}{2\sigma^2}\right)dx$$ I also tried the Itô formula (the actual task is to find $E[e^{cX_T} \mid X_t = x]$ where $X$ is a geometric Brownian motion, but it reduces to the problem above because we are looking at a Markov process), but that did not look very promising either. Can anybody help me?


title: Pre and Post test analysis utilizing a likert scale and SPSS

text: I will openly admit that statistics and statistical programming is not really in my wheelhouse. That said, my job requires it. I am trying to analyze the data gathered from a pre and post test survey where responses were measured with a 1-5 likert scale. 1-strongly agree and 5 is strongly disagree. Essentially we would like to prove that our intervention changes attitudes closer to strongly agree for the majority. I would like to be able to show the percentage of growth for each of our 14 questions from pre to post test. Essentially I need to know what is the change? Again I am using SPSS and in way over my head.


title: help on how to include term $\exp(β_t)/(1+\exp(β_t))$ in AR(2) model

text: I am trying to include a term in an AR(2) model: $$Y_t=\left( a_0+a_1 \frac{\exp(\beta_t)}{1+\exp(\beta_t)}\right)Y_{t-1}+bY_{t-2}+\delta\epsilon_t$$ Can anyone please help me with this? I don't seem to be able include this in the R code.               data=read.table("water.txt",header=T)     n=ncol(data)     koef=matrix(0,n-1,4)      rownames(koef)=names(data)[-1]      colnames(koef)=c("a1","a2","var.a1","var.a2")     for(i in 2:n) {           ser=log(data[,i])           ar2=arima(ser,order=c(2,0,0),xreg=data$year)               koef[i-1,1]=ar2$coef[1]           koef[i-1,2]=ar2$coef[2]               koef[i-1,3]=ar2$var.coef[1,1]           koef[i-1,4]=ar2$var.coef[2,2]     }          data:     yr  water     1986    0.01     1987    -0.63     1988    -0.14     1989    2.52     1990    1.96     1991    0.6     1992    1.82     1993    1.82     1994    0.89     1995    1.32     1996    -1.43     1997    0.75     1998    0.02     1999    0.69     2000    1.65     2001    -1     2002    0.66     2003    0.09     2004    -0.08     2005    0.4     2006    -0.87     2007    1.2     2008    1.1     2009    0.07     2010    -2.57


title: beanplots: plus sign?

text: I would like to draw **bean plots** using the statsmodel package for Python. In the example provided on the documentation, I see a **red plus sign** in each beanplot: ![enter image description here](http://i.stack.imgur.com/7RsJ2.png) What does it represent?


title: Combining errors in a histogram (binned data)

text: I'm processing some data that requires binning before it goes through a regression algorithm. The script is in Python and uses the Numpy `histogram` function, but the code should be self explanatory. For reference, `histogram` outputs either an array containing the integer number of points in each bin, or you can weight by the value of the points in the bin (e.g. a sum). The y errors are standard devs. So the data is binned like:               #Number of items in each bin     (binned, dump) = np.histogram(x,bins)     #Binned data, weighted by value of each point     (xbinned, dump) = np.histogram(x,bins, weight=x)     #etc for y and y error ...      The errors on the y values (this assumes no x errors) are weighted by the variance so the errors are added in quadrature:               (varbinned, dump) = np.histogram(x,bins,weight=error**2)      To finish off the function returns an averaged array of points, e.g.               xbinned = xbinned/binned      What I'm not sure about is how to handle the errors. The result should be square rooted (to convert back from the summed variances), but how is the number of bins factored into the answer? Is it as straightforward as doing:               errorbinned = np.sqrt(varbinned/binned)      Thanks!


title: Chi-square test for checking if values are close to zero?

text: I have a number of set of data points. For each set I have done linear regression to find the lines of best fit for the data. I hypothesise that the gradients of all the lines of best fit should be zero. I took the average and standard deviation of all the gradients and these are the values I calculated: mean = $5.72 \times 10^{-5}$, standard deviation = $5.02 \times 10^{-4}$ To me, a layman, I would say that this shows that the gradients should be zero given the very small mean value and the tight standard distribution. However I've been told that I need to be a bit more rigorous - apparently I should do a chi-squared test to prove my hypothesis. Could someone please give me a set of steps detailing how I would go about this? Thanks.


title: How to translate R to SQL for a Cox Proportional Hazards model?

text: I have built a cox model in R using the coxph function in the survival package, and now I need to replicate the model in SQL for scoring. From my understanding, the model has the form described on the bottom of page 2 of this document, http://cran.r-project.org/doc/contrib/Fox-Companion/appendix- cox-regression.pdf, which gives it semi-parametric flexibility. Since there is an unspecified alpha term, I cannot just take the coefficients and use the model like a typical linear model or generalized linear model (and exponentiate). There are ways to estimate this alpha term, and I believe this added term to the hazard is needed to specify the complete model. If this is the case, how do I get my hands on this alpha term?


title: Using the gap statistic to compare algorithms

text: I want to compare the performances of two clustering algorithms that give me different numbers of clusters. I recently learned about the gap statistic. However, from what I have learned, this statistic is used to find the optimal number of clusters for one algorithm (On that page, for example, it is used to find the best number of clusters for k-means). Is it possible to use it to compare which algorithm clusters gives the best performances? (finds the clustering that minimises the distance in clusters and maximized the distance between them)


title: Conditional distribution of Inverse Wishart

text: Suppose $\begin{bmatrix} K_{11} K_{12}\\\K_{12}^T K_{22} \end{bmatrix}\sim\mathcal{IW}\left(\eta,\begin{bmatrix} \Sigma_{11} \Sigma_{12}\\\\\Sigma_{12}^T \Sigma_{22} \end{bmatrix}\right)$.   1. What is the conditional distribution of $K_{11}|K_{22}$ and $K_{12}|K_{22}$   2. What is the expectation of $K_{12}^TK_{22}^{-1}$


title: Applying Delta Method to generate prediction bands for custom distributions step by step

text: I think that this is an interesting issue for all those who whan to not only adjust one prediction function, but also to obtain the prediction bands to asses the graphical or tabular error I'm trying to apply step by step the Delta Method to my data, however, as I would like to share with you my advances to get to the final. Sumarizing and reviewing the standard method using in comercial software above all, we find: ![enter image description here](http://i.stack.imgur.com/h23aM.jpg) Where the right square root involves the inverse of the Hessian matrix I think, although they call it design matrix. As well as the gradient matrix function. So I share here my code to try to fit a custom non-linear function and asseseing its prediction intervals (confidence and prediction, or int Mathematica language, mean and simple prediction intervals). CODE WITH MY EXAMPLE DATA:               (*Example data*)     In[1]:= data = {31, 46, 70, 87, 87, 93, 114, 128, 133, 134, 143, 155, 161,      161, 163, 177, 181, 207, 207, 226, 302, 315, 319, 347, 347, 362,      375, 377, 413, 440, 447, 461, 464, 511, 524, 556, 800, 860, 880,      954, 5200, 12000};          In[2]:=dist = ProbabilityDistribution[{"CDF",      Exp[-a Exp[-x b] - c Exp[-x d]]}, {x, 0, Infinity},      Assumptions -> {{0 < a < 10}, {0 < b < 1}, {0 < c < 10}, {0 < d <      1}}];          In[3]:= res = FindDistributionParameters[data,           dist, {{a, 3.8}, {b, 0.006}, {c, 0.08}, {d, 0.0002}}];                      (* ACOV is:*)     In[4]:= acov[data_, dist_, paramlist_, mleRule_] :=           Block[{len, infmat, cov}, len = Length[data];          infmat = -D[LogLikelihood[dist, data], {paramlist, 2}]/len /.           mleRule;     cov = Inverse[infmat]];          In[5]:= acov[data, dist, {a, b, c, d}, res];                     (* Gradient function G' *)     In[6]:=D[Exp[-a Exp[-x b] - c Exp[-x d]], {{a, b, c, d}}];          In[7]:=Partition [%,4];     In[8]:= Transpose[%7];                      (* To multiply Transpose-Gradient by acov by Gradient*)          In[9]:=%7 %5 %8                      (*T-student critical value to alpha 0.05 and 42 data*)     Quantile[StudentTDistribution[41], 0.95]      The problem now is to validate this code and to do it working with a numerical example to get the predicted y and the their confidence and prediction intervals. That is, the final matrix have to give a simple number and then we have to calculate the MSE as well as the t-student queantile. Reference posts: http : //mathematica.stackexchange.com/questions/6498/standard - errors - for - maximum - likelihood - estimates - in - finddistributionparameters; http : //stats.stackexchange.com/questions/15423/how - to - compute - prediction - bands - for - non - linear - regression?answertab = active #tab - top http : //stats.stackexchange.com/questions/68080/basic - question - about - fisher - information - matrix - and - relationship - to - hessian - and - s I look forward to your help, thank you.


title: How to compare observed and expected outcomes for continuous data

text: I am working on some data, more specifically some predictions of some outcomes. The predictions vary on the continuous scale, between $-3$ and $3$. They can for example be: $x_1=-2.4, x_2=-2.1, x_3=1.4, x_4=0.4, \cdots, x_n=-1.2$ Now, for each prediction, I also have the corresponding results, $y_1, y_2, \cdots , y_n$, which also are continuous. How can I check if my predictions are good or not? I do not have any knowledge whatsoever on how the predictions are made, thus cannot create any confidence interval around the predictions.


title: Random walk: kings on a chessboard

text: I have a question about the random walk of two kings in a 3×3 chessboard. Each king is moving randomly with equal probability on this chessboard - vertically, horizontally and diagonally. Τhe two kings are moving independently from each other in the same chessboard. Both of them start in the same square, and then they move independently. How could we find the probability in time $n$ both of them are in the same square, as $n$ goes to infinity?


title: Why bother with the dual problem when fitting SVM?

text: Given the data points $x_1, \ldots, x_n \in \mathbb{R}^d$ and labels $y_1, \ldots, y_n \in \left \\{-1, 1 \right\\}$, the hard margin SVM primal problem is $$ \text{minimize}_{w, w_0} \quad \frac{1}{2} w^T w $$ $$ \text{s.t.} \quad \forall i: y_i (w^T x_i + w_0) \ge 1$$ which is a quadratic program with $d+1$ variables to be optimized for and $i$ constraints. The dual $$ \text{maximize}_{\alpha} \quad \sum_{i=1}^{n}{\alpha_i} - \frac{1}{2}\sum_{i=1}^{n}{\sum_{j=1}^{n}{y_i y_j \alpha_i \alpha_j x_i^T x_j}}$$ $$ \text{s.t.} \quad \forall i: \alpha_i \ge 0 \land \sum_{i=1}^{n}{y_i \alpha_i} = 0$$ is a quadratic program with $n + 1$ variables to be optimized for and $n$ inequality and $n$ equality constraints. When implementing a hard margin SVM, **why would I solve the dual problem instead of the primal problem?** The primal problem looks more 'intuitive' to me, and I don't need to concern myself with the duality gap, the Kuhn-Tucker condition etc. It would make sense to me to solve the dual problem if $d \gg n$, but I suspect there are better reasons. Is this the case?


title: What is the default contrast coding for the interaction test in a 3x2 ANOVA?

text: As far as I understand it, the interaction term in a traditional 2x2 ANOVA is great for detecting a disordinal interaction, because the test assumes that the interaction is of form (-2, -2, 2, 2) (Buckless and Ravenscroft, 1990). In cases where I would expect an ordinal interaction, I would then want to use contrast coding for appropriate power. However, my question is how this changes when there is an added level of one factor. What contrast coding does the interaction term in a 3x2 ANOVA assume?


title: Confused about confidence interval

text: I am confused about the concept of confidence interval. Specifically, assume there is a Gaussian variable $X \sim N(\mu, \sigma)$ with $\sigma$ known, and I am interested in the lower-bound $\mu_L$ of the mean with $95\%$ confidence level. I will do the experiment for $5$ times, and observe $X_1$, $X_2$, $X_3$, $X_4$, $X_5$. **Option 1:** I treat each sample separately, and I can compute $\mu_L = X_i - \sigma z$ for each $X_i$. And then I guess there is some way (I don't know how) to compute the actual lower bound from these 5 $\mu_L$'s. **Option 2:** On the other hand, if I take $T = (X_1+X_2+X_3+X_4+X_5)/5$, I can compute $\mu_L = T - \sigma/\sqrt{5}z$. (assuming $T$ is normal, we can use t-stat too.) My question is, is there any method other than option 2 to compute a lower- bound based on the $5$ samples? And for option 1, is there a way to compute the lower-bound based on the 5 lower-bounds computed?


title: Constructing a discrete r.v. having as support all the rationals in $[0,1]$

text: This is the constructivist sequel of this question. If we can't have a discrete uniform random variable having as support all the rationals in the interval $[0,1]$, then the next best thing is: Construct a random variable $Q$ that has this support, $Q\in \mathbb{Q}\cap[0,1]$, and that it follows _some_ distribution. And the craftsman in me requires that this random variable is _constructed_ from existing distributions, rather than created by abstractly defining what we desire to obtain. So I came up with the following: Let $X$ be a discrete random variable following the Geometric Distribution- Variant II with parameter $0<p<1$, namely $$ X \in \\{0,1,2,...\\},\;\;\;\; P(X=k) = (1-p)^kp,\;\;\; F_X(X) = 1-(1-p)^{k+1}$$ Let also $Y$ be a discrete random variable following the Geometric Distribution-Variant I with identical parameter $p$, namely $$ Y \in \\{1,2,...\\},\;\;\;\; P(Y=k) = (1-p)^{k-1}p,\;\;\; F_Y(Y) = 1-(1-p)^k$$ $X$ and $Y$ are independent. Define now the random variable $$Q = \frac {X}{Y}$$ and consider the _conditional_ distribution $$P(Q\leq q \mid \\{X\leq Y\\})$$ In loose words "conditional $Q$ is the ratio of $X$ over $Y$ conditional on $X$ being smaller or equal than $Y$." The support of this conditional distribution is $\\{0,1,1/2,1/3,...,1/k,1/(k+1),...,2/3,2/4,...\\} = \mathbb{Q}\cap[0,1]$. **The "question" is: Can somebody please provide the associated conditional probability mass function?** A comment asked "should it be closed-form"? Since what constitutes a closed- form nowadays is not so clear cut, let me put it this way: we are searching for a functional form into which we can input a rational number from $[0,1]$, and obtain the probability (for some specified value of the parameter $p$ of course), leading to an indicative _graph_ of the pmf. And then vary $p$ to see how the graph changes. If it helps, then we can make the one or both bounds of the support open, although these variants will deprive us of the ability to definitely graph the upper and/or lower values of the **pmf**. Also, if we make open the upper bound, then we should consider the conditioning event $\\{X<Y\\}$. Alternatively, **I welcome also other r.v.'s that have this support(s), _as long as they come together with their pmf_**. I used the Geometric distribution because it has readily available two variants with the one not including zero in the support (so that division by zero is avoided). Obviously, one can use other discrete r.v.'s, using some truncation. I most certainly will put a bounty on this question, but the system does not immediately permit this.


title: Imputing missing observation in multivariate time series

text: Suppose I have a dataframe consisting of six time series. In this dataframe, some observations are missing, meaning at some timepoints all time series contain a NA-value. In R, one possible imputation package that can be used to impute time series data is Amelia. However, this package does not work for observations that are completely missing. Are there other ways to impute my data? For what it's worth, the amount of missing observations is less than 20% of all observations.


title: Interpretation and understanding results from a chi square test for independence

text: I'm making a chi square test on data that has information about students. I want to find out whether there is a relationship between how well the students did on a particular test and the level of dropout from education. I have a 2×2 matrix with the variables `Level in test` which takes the values `level 1` and `level 2`, and the variable `dropout` which has the values `not active` and `active`. I have performed a chi square test for independence on a whole faculty ( _n_ = 1688) and the p value = 0.0042, indicating that there is a difference between students in `level 1` and `2` in relation to their `dropout`. But when I perform the same analysis on each of the 5 educations in the faculty, I get the p values (0.5275, 0.6499, 0.1190, 0.3298, 0.2660). These individual conclusions are that there is no difference in students in `level 1` and `level 2` in relation to `dropout`. Are these conflicting conclusions? What should I trust? (See attached data.) Thank in advance for the help! > ![enter image description here](http://i.stack.imgur.com/5uG1b.png)


title: Reference for $\text{Var}(X)\le (b-a)^2/4$

text: I am not a statistician, but am working a proof for the upper bound of an expression which contains the variance of a variable which obtains its values from a closed interval, [0,1]. I have seen in several places (mostly forums and web pages) the result that $\text{Var}(X)\le \frac{(b-a)^2}{4}$ for $a\le X\le b$. I would like to use, and cite, this result but am having an unusually difficult time finding a solid reference. Could someone here recommend one for me?


title: Maximum of standard deviation, $\sigma$

text: I would like to know how we can determine the maximum value of standard deviation, $\sigma$. Given a range for a population or sample, $[a, b]$ where $b$ is the maximum value and $a$ is the minimum value, I suppose there is also a range of $\sigma$ which I denote as $[0, m]$, i.e. $0 \leq \sigma \leq S$. I think the minimum $\sigma$ has to be 0 because if all samples are the same value then there is no deviation from the mean. But what about the maximum $\sigma$ from the mean? Essentially I want to know how to calculate the lower and upper bounds for $\sigma$, given the range of the population or sample. It would be great if I could have some reference to a paper or book to validate the answers I get.


title: $\sigma^2 \le (\mu-a)(b-\mu)$ for all probability distributions bounded on $[a,b]$?

text: Let $\mu$ be the mean and $\sigma$ the standard deviation of a probability distribution defined on the bounded interval $[a,b]$ (that is, the probability that the random variable lies outside $[a,b]$ is zero). Does the following inequality hold generally for any such probability distribution? $$\sigma^2 \le (\mu-a)(b-\mu)$$ **Motivation:** This inequality holds for the Beta distribution (as is obvious from the formulae here).


title: Variance of a bounded random variable

text: Suppose that a random variable has a lower and an upper bound [0,1]. How to compute the variance of such a variable?


title: Which machine learning algorithm is the slowest BUT surest?

text: CrossPost: https://stackoverflow.com/questions/24301743/which-machine- learning-algorithm-is-the-slowest-but- surest?noredirect=1#comment37556042_24301743 Perhaps my perception of time is augmented by the faster machine speeds these days, but I was wondering if there was a form of machine learning that takes longer but will yield drastically better results on large datasets with lots of noise. I am assuming here that faster convergence somehow has a relationship with the likelihood of becoming stuck in a local extrema. I notice little fluctuation after 100-500 epochs at which point I simply have to restart. I am currently using a feed-forward neural network for both regression and classification. I suppose genetic algorithms seem to be the most time intensive types of "brute force" machine learning. I was also thinking that other types of neural networks could be modified (such as their momentum or learning rate to increase its range over the function). Obviously, I have tried adjusting both of these but this has not solved my problem.


title: How Sensitive Are Neural Networks?

text: CrossPost: https://stackoverflow.com/questions/24301472/how-sensitive-are-ff- neural-networks I am aware of pruning, and am not sure if it removes the actual neuron or makes its weight zero, but I am asking this question as if a pruning process were not being used. On variously sized feedforward neural networks on large datasets with lots of noise:   1. Is it possible one (or some trivial amount) extra OR missing hidden neurons OR hidden layers make or break a network? Or will its synapse weights simply degrade to zero if it is not necessary and compensate with the other neurons if it is missing one or two?   2. When experimenting, should input neurons be added one at a time or in groups of X? What is X? Increments of 5?   3. Lastly, should each hidden layer contain the same number of neurons? This is usually what I see in example. If not, how and why would you adjust their sizes if not relying on using pure experimentation? I would prefer to overdo it and wait longer for a convergence than if larger networks will adapt itself to the solution. I have tried numerous configurations, but it is still difficult to gauge an optimum one.


title: AIC, BIC and GCV: what is best for making decision in penalized regression methods?

text: My general understanding is AIC deals with the trade-off between the goodness of fit of the model and the complexity of the model. $AIC =2k -2ln(L)$ $k$ = number of parameters in the model $L$ = likelihood Bayesian information criterion BIC is closely related with AIC.The AIC penalizes the number of parameters less strongly than does the BIC. I can see these two are used everywhere historically. But generalized cross validation (GCV) is new to me. How GCV can relate to BIC or AIC? How these criteria, together or separate used in selection of penalty term in panelized regression like ridge ? **Edit:** Here is an example to think and discuss:                   require(lasso2)         data(Prostate)         require(rms)              ridgefits = ols(lpsa~lcavol+lweight+age+lbph+svi+lcp+gleason+pgg45,                method="qr", data=Prostate,se.fit = TRUE, x=TRUE, y=TRUE)         p <- pentrace(ridgefits, seq(0,1,by=.01))         effective.df(ridgefits,p)         out <- p$results.all             par(mfrow=c(3,2))             plot(out$df, out$aic, col = "blue", type = "l", ylab = "AIC", xlab = "df"  )             plot(out$df, out$bic, col = "green4", type = "l", ylab = "BIC",  xlab = "df" )             plot(out$penalty, out$df,  type = "l", col = "red",               xlab = expression(paste(lambda)), ylab = "df" )             plot(out$penalty, out$aic, col = "blue", type = "l",                 ylab = "AIC", xlab = expression(paste(lambda))  )             plot(out$penalty, out$bic, col = "green4", type = "l", ylab = "BIC",            xlab= expression(paste(lambda))          require(glmnet)     y <- matrix(Prostate$lpsa, ncol = 1)         x <- as.matrix (Prostate[,- length(Prostate)])         cv <- cv.glmnet(x,y,alpha=1,nfolds=10)         plot(cv$lambda, cv$cvm, col = "red", type = "l",            ylab = "CVM",   xlab= expression(paste(lambda))      ![enter image description here](http://i.stack.imgur.com/HMAXf.jpg)


title: Optimal weights for a linear combination

text: I am trying to find optimal weights for a linear combination. I already know the signs of the weights (+/-). Specifically I have gene expression data where 100 genes belong to a certain pathway. I know which genes are regulators or not (signs) for those 100 genes. I am trying to compute a pathway score for each sample using appropriate weights. Is there any Bayesian type (so that I can use prior information +/-) approach that can be used to obtain weights?


title: What is a good way to present data in a bar-chart when the spacing between bars isn't equal?

text: I have some MATLAB data; the first column is some message length, and the second column is how often that length occurs. Something like this:               15   1     40  12     200  5     215  2      My first instinct to visualize this data was to put it in a bar graph, as so: ![bar graph](http://i.stack.imgur.com/Rcqak.png) The problem with bar graphs, though, is that there is an implied equal spacing. In the case of the picture above, there is an interval of 10 units (presumably years in this case) from bar to bar. My data does not have equal spacing everywhere, as shown in the example above. And the range of the message lengths is large; one chunk of data I just processed has a range of 44 to 7756 with very irregular spacing (anywhere between 10 and 1000 apart). I don't think a histogram is the right choice here either, since the data I'm looking at are so irregularly spaced, there's no good choice of bin size. How should I go about visualizing this data?


title: Principal Component Analysis among matrices

text: I have the following setting. I have _n_ Hermitian Positive Semidefinite (HPSD) matrices, and a metric induced by a matrix norm. I am primarily interested in the Frobenius norm and the operator norm. I want to extract the principal "principal component" for this set of observations, i.e., a 1-dimensional subspace in HPSD such that that the sum of squares of the minimum distances between the matrices with this subspace is minimized. I am not interested in second, third etc components, since they are not even well defined in this space, since I have not defined a scalar product for simplicity. In the case of the Frobenius norm, the problem can be reduced to traditional PCA, by using as input vectors the stacked versions of the input matrices. But in the case of the operator norms, I can't find a strategy to attack the problem. Questions:   1. Has anyone seen this specific problem before? Recommendations and references are highly appreciated.   2. Has anyone dealt with computation of PCA in the case of non-euclidean distances?


title: Intuitive explanation of the bias-variance tradeoff?

text: I am looking for an intuitive explanation of the bias-variance tradeoff, both in general and specifically in the context of linear regression.


title: Calculate z-score from odds-ratio

text: I want to calculate the z-score for a distribution for which I know the odds- ratio, and the p-value. I also know the upper and lower 95% confidence limits for the odds-ratio. An example with two-data points is shown here.               OR       OR_95    OR_95U   p-value     0.997804 0.970573 1.025798 0.876215     1.039562 1.010116 1.069866 0.00815


title: Aggregation of propensity scores with varying reliability

text: When trying to estimate the number of sampling units with an attribute, is there a good algebraic way to aggregate over propensity scores for that attribute which each have their own error? For example, when the propensity scores may be calculated with varying amounts of information from each sampling unit and each has its own standard error on that propensity. Or they could be expected values of beta distributions. In the latter case I know I can simulate beta-bernoulli outcomes from each sampling unit and add up the results many times; but is there a consistent estimator of the result of this difficult to scale process? In short, how do people aggregate propensity scores of varying reliabilities? Edit: I suppose I worded it poorly; the data I have is all either binary or categorical and each observation is accompanied by the probability it was observed correctly. So suppose I have 5 persons; 3 of which had a value of 1 for an attribute, 2 of which had a value of 0 for that attribute, each of which had a probability .8,.81,.82.,.83,.84 respectively of being observed correctly. What is the expected value of p(having that attribute)?


title: Residuals in R using auto.arima and forecast package

text: I am fitting a model using the `auto.arima` function in package `forecast`. I get a model that is AR(1), for example. I then extract residuals from this model. How does this generate the same number of residuals as the original vector? If this is an AR(1) model then the number of residuals should be 1 less than the dimensionality of the original time series. What am I missing? Example:               require(forecast)     arprocess = as.numeric(arima.sim(model = list(ar=.5), n=100))     #auto.arima(arprocess, d=0, D=0, ic="bic", stationary=T)     #  Series: arprocess      #  ARIMA(1,0,0) with zero mean               #  Coefficients:     #          ar1     #       0.5198     # s.e.  0.0867          # sigma^2 estimated as 1.403:  log likelihood=-158.99     # AIC=321.97   AICc=322.1   BIC=327.18     r = resid(auto.arima(arprocess, d=0, D=0, ic="bic", stationary=T))     > length(r)       [1] 100      Update: Digging into the code of auto.arima, I see that it uses Arima which in turn uses `stats:::arima`. Therefore the question is really how does `stats:::arima` compute residuals for the very first observation?


title: Pooled OLS with time-variant regressors

text: Could someone please explain how a pooled OLS is set-up when the regressors are aggregate variables that are time-varying only (i.e. do not vary for the cross-section at each time _t_ )? If the pooled regression is expressed as: $$\;Y_{i,t} = \;X_t \;\beta_i + \;\epsilon_{i,t} $$ The response vector $Y$ is stacked: in other words, we stack the time-series of each $Y_i$ to form only one response vector of dimension (TN x 1). The regressor X should be stacked similarly to obtain a vector of dimension (TN x 1). The issue here is that X is time-varying only, which means that the stacked vector would be `rep(X, N)`: i.e. the same vector X is repeated N times. Would it be valid to run a pooled OLS with such a stacked regressor? Thank you,


title: How do I write the mathematical equation for ARIMA(0,0,5)(1,0,1) period 12?

text: I would appreciate if someone could help me write the mathematical equation for the seasonal ARIMA (0,0,5) x (1,0,1) period 12. I'm a little confused with how to go about this.


title: chi square test or Z test?

text: SO I have an outcome from a biological experiment, where I have counted the number of dead cells upon a mutation. For example, in once case - out of 120 cells we saw 16% of the dead cells, but our expectations were around 10%. So to check the difference is significant I have performed 'Z-test for changes in proportions'. But my question is, whether is it possible to do ch-squared test on a single experiment like the above ? I have googled a lot, but all examples I see use a contingency table or involve multiple variables.


title: SVM failing entirely with when test set is varied

text: I am experiencing a strange problem with varying the test set size. This is mildly confusing to explain, but I'll do my best. I'm using octave to train an SVM on timeseries-like data and it's generally yielding good results. However, occasionally when we've introduced "bad" features, the SVM will "give up" almost instantly during training and produce a model that effectively pins predictions on a single feature. In these cases, the train completes in ~1 second vs the normal 4-5 minutes. I've never really known what causes it, but I just presumed it couldn't find a useful model, which seemed reasonable. The problem I'm facing now is that it appears something as simple as normalization can be the difference between a very effective model and a "give up instantly" result. The specific case is this: I have a dataset of say 5000 training examples. They are split into a 90/10 split for train/test and work well. Similarly, if I split the data set 95/5, that also works fine (better, in fact for test performance, simply because it's smaller). However, weirdly, if I split the dataset 90/5/5 and throw out the last 5%, the machine refuses to train. So all I'm doing is _removing_ some of the test set, which is causing the entire process to fail. To make things even more confusing, I found that normalization is the cause: If I normalize (zero-mean + scaling) then slice, everything works fine, however, if I slice, _then_ normalize, the machine fails. Basically, my conclusion is that the exclusion of certain training examples in the test set is causing normalization to (naturally) come up with differently values scaled that are for some reason problematic. I have eyeballed the min/max/mean of the scaled values and none of them vary by more than a few percent, so I'm somewhat at a loss of what to pursue next. Any guidance of how to go about diagnosing this or what may be happening would be most appreciated. Cheers! **Update:** After some further work on this, I've realised that the machine "gives up" instantly if there is a feature with all values being zero (this can happen either on certain subsets of training data or just bad luck). Although I admit a feature of all zeroes is of little value, is there a specific reason why this would cause the algorithm to behave as described? Incidentally, we're using a Gaussian kernel if it makes a difference.


title: Question re permutation tests for a relatively complex data set

text: I regret that I cannot be overly explicit due to confidentiality agreements. I will try to state my problem by analogy. Suppose that you have two groups of people, one with a condition and one normal. You make roughly 80 measurements of circumference of different body parts on each person. Clearly, these measurements will all be correlated, some more highly than others (measurements that are close to each other will tend to be more correlated than others). You have an _a priori_ hypothesis that the largest differences between groups will occur at certain places (e.g. "patients with the condition will have larger arms than patients without"). As a first step, we looked at t-tests comparing the measurements between the groups; there were quite a few significant differences and most of the largest differences were in the hypothesized region. My idea was to use some form of permutation test to account for multiple testing while accounting for the correlations. This would not only be to test the specific hypothesis but to correct for multiple testing on all 80 measurements (instead of using Bonferroni). However, I need help: 1) Is this reasonable? 2) If reasonable, any guidance on how to do it in either `R` or `SAS`? (any hints or references welcome) ADDING A REFERENCE One paper that does something similar to what I want is Nichols and Holmes But that is for FMRI data, which is somewhat different; and they give a program in Matlab, which I don't have or know. Also, they don't consider exactly the sort of problem I have.


title: How to explain difference of importance between feature selection and model quality?

text: I have a data collection with a mixed feature set consisting of both numerical features and text features. The number of numerical features is quite small, i.e., 6, comparing to the number of text features, i.e., 2000. ( The data set is a csv file, where there are 6 columns consisting of digit values only; other columns are just plain texts). We found that removing these numerical features (columns) in fact improves the performance. However, when we conduct the feature selection based on the information gain, looks like those numerical features gets more weight. How to explain this scenario?


title: Advice for mixed model with multilevel effects in R

text: Im migrating from SAS to R with a few difficulties, not to mention some large gaps in statistical knowledge. Im investigating the effect of ethnicity on blood pressure, in a longitudinal study. For all ethnic groups, it is evident that their blood pressure declines the first 3-4 years and then it increases steadily thereafter (no other trend observed). I have 100.000 observations, from 10.000 unique individuals. Data is heavily unbalanced; some individuals have 1 observation while others have 20. Observations are gathered at different time points. Fixed covariates: age, sex, treatment, BMI, ethnicity. Random covariates: ethnicity Repeated measure unit: individual (ID). How would You model this? I'm interested in the ethnic differences and must therefore have ethnicity as a fixed covariate. But ethnicity could be appreciated as a level in terms of multilevel models, and thus modeled as random in mixed models. Så my subjects are, if im not wrong, nested within ethnic groups. I have read instruction manuals for both nlme package and lme4 package. I decided to go with lme4, despite non-linear trend in blood pressure but tried to adjust for this by taking the second polynomial of duration (time*time). How would you model this? E.g:               lmer(hba1c_up ~ age + sex + (1|Ethnicityx/ID)) ?          lmer(hba1c_up ~ age + sex + (1|ID) + (1|Ethnicityx))      Any suggestions? Help would be immensely appreciated! /Adam


title: How to adjust measured values depending on control values?

text: I have taken 20 photographs. Each photograph contained a museum specimen of a bird and a colour chart. In each photograph, I measured the brightness of a specific part of the bird's plumage and the brightness of a specific part of the colour chart. In each photograph, the value of the colour chart should be 50. However, due to light conditions varying from photograph to photograph, the brightness values of the colour chart are sometimes below 50 and sometimes above 50. Therefore, the brightness of the bird's plumage should be adjusted accordingly. The dataframe is a representation of the brightness values obtained for the colour chart and the brightness values obtained for the bird plumage:               set.seed(50)     data.frame(brightness.colour.chart = c(sample(40:55, 20, replace=T)),                brightness.bird.plumage = sample(seq(1, 10, 0.1), 20))             brightness.colour.chart brightness.bird.plumage     1                       51                     6.6     2                       47                     2.7     3                       43                     7.3     4                       52                     8.4     5                       48                     3.6     6                       40                     6.5     7                       51                     5.9     8                       50                     3.1     9                       40                     9.6     10                      41                     4.0     11                      46                     3.9     12                      44                     4.2     13                      50                     4.4     14                      41                     5.8     15                      44                     9.9     16                      50                     5.7     17                      53                     1.9     18                      45                     6.1     19                      41                     3.7     20                      42                     4.3      What is the best way to adjust the brightness values for the bird plumage so they are constant across each photograph?


title: Expectation of a generalization of Dirichlet distribution

text: For the standard Dirichlet, the expectation of $X_i$ is $\alpha_i/\alpha_0$, where $\alpha_0 = \sum_i \alpha_i$ (http://en.wikipedia.org/wiki/Dirichlet_distribution). I am considering the following generalization. Suppose we are playing a simple poker game as follows. We are player 1 and observing player 2's plays. Player 2 can be dealt one of two hands with probability 1/2 (or more generally probability $p$) -- K or Q (player 2 sees this and player 1 does not). Then player 2 selects one of $i$ actions. Let $X_i$ denote the rv for the probability player 2 plays action $i$ with a K, and let $Y_i$ denote the rv for Q. Suppose that he is following a static probability distribution for all rounds, and that the $X_i$'s and $Y_i$'s are independent (his strategy for a K is independent of his strategy for a Q). Player 1 only observes the action $i$ of P2, and not his card. I am trying to find a closed form for $E[X_i]$ and $E[Y_i]$. This generalization differs from the "Generalized Dirichlet distribution," which has a relatively simple closed-form solution for the expectations (http://en.wikipedia.org/wiki/Generalized_Dirichlet_distribution). I think the pdf of this distribution is the following, where $p$ denotes the probability P2 is dealt a K, $\alpha_i$ is the observed number of times he has taken action $i$ so far, and $N(\alpha)$ is some normalization constant I'm not sure how to compute: $$f(x,y;\alpha) = N(\alpha) * \prod_i [p x_i + (1-p)y_i]^{\alpha_i}$$ where $x_i, y_i \geq 0, \sum_i x_i = 1, \sum_i y_i = 1.$ One can apply the binomial theorem here, but I am not sure how to proceed and if that will even help. Somehow the gamma/beta function will need to come in to play. Thanks a lot.


title: Is cross-validation still valid when the sample size is small?

text: I have what seems to be a very basic confusion about cross-validation. Let us say I am building a linear binary classifier, and I want to use cross- validation to estimate the classification accuracy. Imagine now that my sample size $N$ is small, but the number $k$ of features is large. Even when features and classes are randomly generated (i.e. "actual" classification accuracy should be 50%), it can so happen that one of the features will _perfectly_ predict the binary class. If $N$ is small and $k >> N$, such a situation is not unlikely. ~~In this scenario I will get 100% classification accuracy with any amount of cross-validation folds, which obviously does not represent the actual power of my classifier, in a sense that probability to classify a new sample correctly is still only 50%.~~ [Update: this is **wrong**. See my answer below for the demonstration of why it is wrong.] ~~Are there any common methods of dealing with such a situation?~~ For example, if I want to assess statistical difference between my two classes, I could run MANOVA which in case of two groups reduces to computing Hotelling's T. Even if some of the features yield significant univariate differences ("false positives"), I should get an overall non-significant multivariate difference. However, I do not see anything in the cross- validation procedure that would account for such false positives ("false discriminants"?). What am I missing? One thing that I can think of myself, would be to cross-validate _over features_ , e.g. to select random subset of features (in addition to randomly selecting a test set) on each cross-validation fold. But I do not think such an approach is often (ever?) used. **Update:** Section 7.10.3 of "The Elements of Statistical Learning" entitled "Does Cross-Validation Really Work?" asks exactly the same question and claims that such a situation can never arise (cross-validation accuracy will be 50%, not 100%). So far I am not convinced, I will run some simulations myself. [Update: they are right; see below.]


title: Cumulative standard deviation calculation

text: I have a calculated daily survival probability, and associated SE. From this I can obtain the standard deviation. I have to report the cumulative survival over a specific period. So for example, for a 38 day period, where daily survival rate (DSR) = 0.97, and SD = 0.16, the probability of survival over the period is simply 0.97^38 = 0.314. Is there a way to calculate the cumulative SD as well? Is it appropriate to do so? Note too that in this example DSR is constant, but DSR and SD can, in some situations, vary over the period. This may be pretty straight forward, but I can't seem to find a comprehensible answer. Thanks for any assistance.


title: Proposal distributions for covariance matrices in MCMC implementation of hierarchical models

text: In a MCMC implementation of hierarchical models, with normal random effects and a Wishart prior for their covariance matrix, Gibbs sampling is typically used. However, if we change the distribution of the random effects (e.g., to Student's-t or another one), the conjugacy is lost. In this case, what would be a suitable (i.e., easily tunable) proposal distribution for the covariance matrix of the random effects in a Metropolis-Hastings algorithm, and what should be the target acceptance rate, again 0.234? Thanks in advance for any pointers.


title: Project feature vectors into 2D plane for visualization

text: I would like to project my observations which consist of more than 2 variables into a scatter plot. Some time ago I saw an R package that does this by reducing the dimensions (possibly using PCA and other algorithms) and plots similar observations close to each other. I have been googling to find how to accomplish that but haven´t found anything.


title: What programming language for statistical inference?

text: just for curiosity... What language is used most here? R? MATLAB? Python? Java? What for prototype or for production? For example I think MATLAB is mostly used for prototyping, python for both prot. and production...


title: Integrability of a sequence of iid random variables

text: I'd really appreciate some hints on the first part of the following question: Let $f_n, n\in \mathbb{N}$ be a sequence of iid random variables over $(\Omega, A,P)$. That is, $P(\\{f_1 \in A\\})=P(\\{f_2 \in A\\})= ···$ for all $A \in \mathbb{B}$. Show that, if there is a $c \in \Re$ such that $P\Bigg(\Bigg\\{\lim_{n \to \infty}\frac{1}{n}\sum\limits_{i=1}^n f_i=c\Bigg\\}\Bigg) >0,$ then $P(\\{|f_n|>n$ for infinitely many $n\\})=0.$ Then conclude that $f_1$ (and therefore also $f_2,f_3,..$ because of the identical distribution) is integrable.


title: A good summary measure on analyzing repeated measures of longitudinal data?

text: I have an interesting question on longitudinal data and I'm looking for a proper summary measure that enables me to answer my question of interest without knowing/using anything about correlated data. Suppose we have two different samples - one treatment and one control. At each trial, subjects should solve a puzzle and we count number of attempts to solve the puzzle. Number of attempts could be 1, 2, 3, or if the puzzle is not solved after 3 attempts the subject is failed in this trial. We repeat the trial 10 times so that each subject gets to solve 10 puzzles. Obviously, we have correlated data since we have 10 measurements on each subject. My hypothesis of interest is whether there is any learning effect. Would you think of any summary measure on each subject (each subject has 10 measurements) that summarizes the measurements on that subject into a one measurement so that we can avoid (or at least attempt to avoid) getting correlated data? Thank you,


title: k-fold cross validation vs k times hold-out validation

text: I am facing the evaluation of a genetic programming algorithm. I am using the Proben1 cancer1 dataset to evaluate the models created by this algorithm. This dataset contains 699 samples, which is currently divided into 50% training, 25% validation and 25% test data. Many academic articles use k-fold validation for evaluation of the resulting models. I do understand the creation of k models to reduce the variance in the error rating percentage. However, I do not understand why it would not be preferable to do x times the hold-out (k=2) method, where every time the data is randomly partitioned into training and test data. The reason for my lack of understanding is that I think that models evaluated with a higher order k generalize less due to simple fact that they are trained with a higher percentage of the data. Given my data set of 699 samples which of the following two methods would be preferable and why;   * the division in training/validation/test and perhaps this repeated x times randomly assigning sample to each set with a test set of 25%    * 10-fold cross validation but thus with a test set of just 10%.


title: Derivation of closed form lasso solution: KKT optimality conditions

text: As defined in this thread Derivation of closed form lasso solution I would like to see a derivation of the closed form soft-thresholding solution for the LASSO using the KKT optimality conditions. Does anybody know how to do it?


title: How can one extract meaningful factors from a sparse matrix?

text: I am interested in finding some practical (and reasonably well accepted) techniques for finding the underlying factors of a sparse matrix. Specifically, I have a very large sparse matrix whose cells appear to be populated from an approximately geometric distribution. In its natural form the matrix is square. The values in the cells represent item x item co- occurrences under case 1 over the diagonal and under case 2 under the diagonal. If necessary I can subset the matrix to particularly interesting items in order to make it rectangular. I believe that there are meaningful factors underlying this structure. However my understanding is that because the matrix is sparse factor analysis is not an appropriate approach. What approach can I take that will make it most likely that I can find interpretable patterns in the data? I saw that there was another question asking for references on sparse variants of PCA, but I think I'm looking for something more akin to an obliquely rotated factor solution. I'm willing to dig into suggested readings somewhat, but my prior experience with factor analysis (and related techniques) is limited, and I prefer a relatively straightforward answer (one with R code is even better).


title: Level of uncertainty when measuring a period of time?

text: I need to create a raw data table for my science experiment, but I'm unsure what the level of uncertainty would be when measuring days. How would I determine the level of uncertainty when measuring a period of time?


title: Is it needed to normalize data before rule model extraction algorithms like ID3?

text: I will use naive Bayes or decision tree that gives rule model both. Is is necessary to normalize data before working with such algorithms.


title: Violation of normality assumption

text: I am trying to run an econometric panel data (random effects) model with about 950 observations (so not a small dataset). My data consists of different European public companies and of a few financial statistics of theirs during the last 10 years. However, I am having enormous amount of trouble with the following aspect: According to Jarque-Bera test (p value of 0) I am violating the assumption of normality. While my skewness seems small enough (under 3), the kurtosis has a value of 14. However, I am very puzzled as to what to do now. As one of the options I read that you could take logarithms, but one of my variables (market value added) is often negative and Eviews, obviously, doesn't let me opt for that option. Is there any other solutions to this problem or can anyone offer some advice to an econometric novice? Thank you in advance.