title: Rescaling exponentially distributed variables before clustering?

text: I want to cluster data that contains binary variables, exponentially distributed (power law) variables, and normally distributed variables. I'm considering preprocessing the data in the following way and wondering whether it's reasonable. 1) shift the binary variables so that they have mean zero. no rescaling. 2) standardize the normally distributed variables, but divide by twice the standard deviation rather than by once. this is based on Gelman, A. "Scaling regression inputs by dividing by two standard deviations. Statistics in Medicine 2008 27:2865-2837. [for clarification: their point is that binary variables can be interpreted directly as indicators and can be left unscaled. But in order for numerical inputs to be interpreted in the same way, they should be divided by twice the standard deviation -- that way, both numerical and binary variables have similar standard deviations (if the binary indicators are not strongly skewed; if they are, then there is not much improvement, but also no harm).] 3) take log(x+1) [or sqrt(x)] for any exponentially distributed variable x, and standardize in the above way. Mainly, I'm worried about (3). I'm told that it's standard practice to simply standardize exponentially distributed variables. Is it standard practice because it's correct or because it's a good enough approximation? Also, would it be better to apply k-medoid clustering rather than k-means?


title: Uncertainty of conditional probability evaluated from sample

text: I'm doing some some analysis of an arbitrary string of text, modelling it as a Markov chain where the state is simply the value of the previous character. Call the current character `c` and the previous character `p`; then it is trivial to calculate `P(c | p)` for the given sample text. However, if there are not many cases of a certain (p, c) pair, or in general not many samples used. I know qualitatively that the fewer the samples, the higher the uncertainty/error in this conditional probability is. However, how can I quantify this?


title: Standardizing feature vectors for regression

text: Suppose I have a data set with the following structure: Each row of the data set indexes a town. The first column/feature variable is the total population while the other feature variables include the count of people who own various items (one feature variable for cars, one for home appliances, etc.), while still others measure average income etc. Now, it is often necessary to 'transform' the feature vectors before running some sort of regression algorithm on the data, for example standardizing them. Suppose the towns have very disparate populations (call this feature $X_1$ and let town $i$ have value $X_1^i$). Consider the feature vector, say $X_2$ measuring the number of some X in each town. My question is: > Should one, in general, first transform $X_2$ in proportion to the total > population of the towns, that is $X_2^i \mapsto \frac{X_2^i}{X_1^i}$ and > then standardize the column by $X_2^i \mapsto > \frac{X_2^i-\bar{X_2}}{\hat{\sigma}_{X_2}^2}$ or, even simply scaling the > values to the interval $[\max(X_2), \min(X_2)]$? The reason I am asking the question is: I can imagine a case where despite the towns having very different population counts, there is an item which have roughly the same count in each town. In which case, if we were to simply standardize the columns, it will reduce the values in the feature column to zero (or nearby) and intuitively, there will be tremendous loss of information. Assume that I know that $X_1$ is collinear with $X_2$ and I won't be using that feature.


title: Java implementations of the lasso

text: Are there any open-source Java implementations of lasso or least angles regression? Pure Java code would be best, but clean implementations in other languages would also be of interest. I am already aware of the existence of a variety of R packages than can do lasso/LAR fits. Thanks.


title: Explaining variation in a dependent variable based on a factorial experiment

text: I have run a factorial type test in a processing plant and have run a forward and backward step regression in R. How can I use the regression results and the anova created from the regression to know what percent of the measured variation of the dependent variable was caused by the purposeful manipulation of the independent variables?


title: Lognormal posterior?

text: If I have a relationship: $ y_t = a + \theta b_t \epsilon_t,$ where I observe $y_t$ and $b_t$. $a$ is a known parameter, $\theta$ is an unknown parameter with prior distribution at time $t$: $\theta \sim logN(\mu_t, \sigma_t^2)$, and $\epsilon_t$ is an IID random variable with distribution: $\epsilon_t \sim logN(0,\sigma_\epsilon^2)$. Is it possible to obtain an analytic log-normal posterior for $\theta$? If not, do there exist other possible conjugate prior distributions with non- negative support (and whatever distribution on the noise term $v_t$ that's necessary)?


title: How tho choose the number of components in a Bayesian Hidden markov model

text: I'm implementing a bayesian Hidden markov model. I now face the problem of how to choose the number of components. I have two problems: 1) which index is better to use? 2) suppose i decide to use the BIC and that I estimates the model with an MCMC procedure. Have I to compute the BIC for each posterior sample and then take the mean value? or i have to compute the BIC only once using a punctual estimates of the parameters Thanks


title: how to compute odds ratio

text: Given a generic classification model $y=f(x_1,x_2,..,x_p)$ where $y\in \left\lbrace 0,1 \right\rbrace$ is it possible to compute the odds ratio for each variable? A theoretical explanation and an example with R code would be the best answer


title: Ranking algorithm of different weighted resources

text: I have a list of items that represent cells, now i query each cell value to get a set of possible categories that this cell might belong to, these categories are weighted. For Example:               Microsoft is: {Software=83.543266, Video Game Platform=132.297455, Degree=71.624733, Organization=208.17038, Programming Language=66.998901, File Format=67.159828}          apple is: {Physicist=50.008327, Organism Classification=145.239532, Video Game Platform=86.5653, Computer=60.457172, Organization=199.53299, Record label=53.477039}          Banana is: {Musical Group=27.214769, Organism Classification=136.063583, Writer=27.566528, Musical instrument=29.285137, Organization=36.962833, TV Program=32.693085}          Blackberry is: {Software=23.69421, Organism Classification=64.572891, Operating System=90.257011, City/Town/Village=16.836142, Computer=34.484653}      The weights that i get for each result are independent at cell level, so Microsoft having a very high confidence as an organization does not relate to the confidence score of Apple organization confidence. Now i have two problems, the first is to compare two different cells and compute a similarity score, for example how similar apple is to Microsoft having these possible categories for each. Second, i want to deduce a global common category for the column that represents these cells, for example if i run a simple algorithm that will calculate the number of categories occurrences and the average of their confidences will not be accurate. I will appreciate any help.


title: Generalization of one-sample Kolmogorov-Smirnov test for non identically distributed data?

text: Is there a standard approach for testing goodness-of-fit to a probability distribution when the samples are independent but not identically distributed? In other words, given the data $\\{\\{x_1, y_1\\}, \\{x_2, y_2\\},..., \\{x_n, y_n\\}\\}$ and the candidate density function $f(x;y)$, what is an accepted approach for testing whether the $x_i$ are sampled from $X|y_i \sim f(x; y_i)$? What are the assumptions of said approach? **Edit** In the comments, whuber suggests it is more appropriate to call $y_i$ a parameter for the distribution on $x_i$ and requests more information about how $f(x; y)$ depends on $y$. NB, I'm interested in "standard methods" more than a custom solution, so the answer should not depend on the details of $f$. That said, in my application, $f$ is a mixture of $m$ Dirichlet-Multinomial distributions (so $x$ is discrete, I asked about continuous variables above in order to reference the K-S test). The vector $\\{y_{i,1},y_{i,2},...,y_{i,m}\\}$ is the mixture weighting. The parameter for each Dirichlet-Multinomial distribution is independent of $y_i$.


title: Why are my p-values so high?

text: I am getting astronomically high p-values for all the coefficients in a logistic regression model. I am not sure why they are so high:   * 0.7096    * 0.4441    * 0.9783    * 0.7826    * 0.6890    * 0.7538    * 0.8017    * 0.9332    * 0.9564 I honestly cannot find out why they are so high. My dependent variable consists of ratios between 0 and 1, and my independent variables vary between continuous data, ordinal data, and categorical data. Just adding a bit more information. I am creating a model to predict the proportion of those who evacuated to the remaining population in an area. For example, 10 percent evacuated / 90 percent left, and so on, over 37 hours. The data is sequential where each hour before effects the value of each hour after it. In the literature they call it the sequential logit model. Currently I am using Matlab's generalized linear model function using a binomial distribution linked with a logit. I am using 9 independent (predictor) variables and one dependent variable that is continuous between 0 and 1. More information! The model being developed operates only on people that have already made the decision to evacuate. So I am not testing for if they will leave or not but when they will leave. Hopefully this clarifies what I am trying to do. I downloaded R and I ran a logistic test. Here is what I got               Deviance Residuals:           Min        1Q    Median        3Q       Max       -0.41146  -0.12817  -0.00301   0.10222   0.47937            Coefficients:                        Estimate Std. Error z value Pr(>|z|)     (Intercept)       -36.87241  570.55243  -0.065    0.948     Hour               -0.27353    1.00636  -0.272    0.786     LN.TU.             -3.96075    8.57513  -0.462    0.644     as.factor(Evac1)1   3.34697   13.38238   0.250    0.803     as.factor(Evac2)1   1.42544    3.51690   0.405    0.685     Pressure            0.04122    0.59701   0.069    0.945     WindSpeed          -0.01583    0.23526  -0.067    0.946     as.factor(Time1)1  10.11888   11.22222   0.902    0.367     as.factor(Time2)1   9.89199   10.45510   0.946    0.344     as.factor(Time3)1  10.72968    9.48993   1.131    0.258          (Dispersion parameter for binomial family taken to be 1)          Null deviance: 14.5671  on 36  degrees of freedom     Residual deviance:  1.5175  on 27  degrees of freedom     AIC: 27.927          Number of Fisher Scoring iterations: 9      My p-values remain very unappealing. I am beginning to think that instead of trying to fix this, I should just believe that what I am doing is wrong and take another approach. Thanks all who answered and helped!


title: Modeling the shape of the power function

text: In cognitive psychology, it has been argued that the learning curve of a skill follows the power function, in which the practice of the skill yields progressively smaller decrements in error. I want to test this hypothesis with the data I have, and to do this I need to fit a power law function to the data. I heard that the power function can be fit in R by `lm(log(y) ~ log(x))`. An answer to this post, however, suggests using `glm(y ~ log(x), family = gaussian(link = log))` ~~, and indeed the resulting fit prefers the`glm` approach~~ ( **EDIT:** See comments. This is not true.).               # Define a function to compute R-squared in linear space      > rsq <- function(data, fit) {     + ssres <- sum((data - fit)^2)     + sstot <- sum((data - mean(data))^2)     + rsq <- 1 - (ssres / sstot)     + return(rsq)     }          # generate data and fit a power function with lm() and glm()     > set.seed(10)     > lc <- (1:10)^(-1/2) * exp(rnorm(10, sd=.25)) # Edited     > model.lm <- lm(log(lc) ~ log(1:10))     > model.glm <- glm(lc ~ log(1:10), family = gaussian(link = log))     > fit.lm <- exp(fitted(model.lm))     > fit.glm <- fitted(model.glm)          # draw a graph with fitted lines     > plot(1:10, lc, , xlab = "Practice", ylab = "Error rate", las = 1)     > lines(1:10, fit.lm, col = "red")     > lines(1:10, fit.glm, col = "blue")     > legend("topright", c("original data", "lm(log(y) ~ log(x))",      + "glm(y ~ log(x), family = gaussian(link = log))"),      + pch = c(1, NA, NA), lty = c(NA, 1, 1), col = c("black", "red", "blue"))      ![enter image description here](http://i.stack.imgur.com/F28SM.jpg)               # compute R-squared values for both models     # (EDIT: See comments. These are not good measures to use.)     > rsq(lc, fit.lm)     [1] 0.9194631     > rsq(lc, fit.glm)     [1] 0.9205448      ~~From the $R^2$ values, it seems that the`glm` model has a better fit. When the procedure was repeated 100 times, it was always the case that the `glm` model has a higher $R^2$ value than the `lm` model~~ ( **EDIT:** See comments. $R^2$ is not a good measure for this purpose). At this point I have two questions.   * What is the difference between the `lm` model and the `glm` model above? I thought a generalized linear model with the log link function is practically the same as regressing log(y) in general linear models.   * I'm trying to compare the fit of the power function against other models (e.g., linear models as in `lm(y ~ x)`). To do this, I calculate the fitted values in linear space and compare the $R^2$ values computed in linear space, just like I did to compare `lm` and `glm` models above. Is this a proper way to compare this type of models? My final question concerns the estimation of the four-parameter power law function given in this paper (p.186; I modified the notation a little). $$ E(Y_N) = A + B(N + E)^{-β} $$ where   * $E(Y_N)$ = expected value of $Y$ on practice trial $N$   * $A$ = expected value of $Y$ after learning is completed (i.e., asymptote)   * $E$ = the amount of prior learning before $N = 0$ I have $N$ and $Y$ (`1:10` and `lc` respectively in the code above), and want to estimate $A$, $B$, $E$, and $-β$ from the data. In the example above, I assumed $A$ = $E$ = 0 mostly because I wasn't sure how to estimate four parameters concurrently. While this isn't a huge problem, I would like to estimate $A$ and $E$ from the data as well if possible. So my final question:   * Is there a way to estimate the four parameters in R? If not, is there a way to estimate three parameters, excluding $A$ from the above (i.e., assuming $A$ = 0)? The paper I linked above says the authors used the simplex algorithm to estimate the parameters. But I am familiar with neither the algorithm nor the `simplex` function in R. I'm aware of Clauset et al's work and the `poweRlaw` package in R. I believe both of them are dedicated to modeling the distribution of one categorical variable, but neither models the curve of the power law. Any help is appreciated.


title: Covariance term in simple linear regression

text: I am trying to derive the expression for the variance of $\hat{\beta_0}$ in simple linear regression. I substitute $\bar{y} - \hat{\beta_1} \bar{x}$ for $\hat \beta_0$, but in the intermediate steps the covariance term $\text{Cov}(\bar{y}, \hat{\beta_1})$ comes up and I don't know how to deal with it. Any help would be appreciated!


title: Using all categorical values as regression predictors in R

text: Let's say that you're fitting a model to predict some score as a linear function of two categorical inputs. Something like the following data set:               score   x   y     ---------------------     46.12   A   B     42.61   A   F     39.14   B   A     57.72   B   C     48.65   B   D     56.29   B   E     45.45   B   F     42.43   C   B     52.68   C   E     43.03   D   B     52.54   D   C     44.32   E   B     54.43   E   C      For interpretability, I would like to treat each input category as an indicator variable. That is, each data point has 3 input dimensions: the intercept, the x category, and the y category (all with value 1). However, the default coding behavior for lm() is to fold the first value of the first variable in with the intercept... which throws away the xA predictor in my example:               Call:     lm(formula = score ~ x + y, data = d)          Residuals:              1          2          3          4          5          6          7          8          9         10         11         12         13       4.281e-01 -4.281e-01  8.882e-16  1.171e-01 -4.441e-16 -5.452e-01  4.281e-01 -5.452e-01  5.452e-01  2.086e-01 -2.086e-01 -9.143e-02  9.143e-02           Coefficients:                 Estimate Std. Error t value Pr(>|t|)         (Intercept)  37.1562     1.0316  36.019 4.71e-05 ***     xB            1.9838     0.7379   2.689 0.074504 .       xC           -2.7167     0.8324  -3.264 0.047005 *       xD           -2.8705     0.8212  -3.495 0.039610 *       xE           -1.2805     0.8212  -1.559 0.216828         yB            8.5357     0.9824   8.688 0.003209 **      yC           18.4629     0.9439  19.561 0.000292 ***     yD            9.5100     1.0195   9.328 0.002609 **      yE           17.6952     0.9569  18.492 0.000345 ***     yF            5.8819     0.9569   6.147 0.008662 **      ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Residual standard error: 0.7209 on 3 degrees of freedom     Multiple R-squared:  0.9964,    Adjusted R-squared:  0.9856      F-statistic:  92.5 on 9 and 3 DF,  p-value: 0.001657      Of course I can achieve what I want with other machine learning libraries like scikit-learn, but they don't provide standard errors or p-values out-of-the- box. Is there a way to force R to explicitly code using every categorical value in the model (in addition to the intercept)? **EDIT** : To be clear, I have multiple categorical inputs, _and_ I want to keep the intercept. I did a pretty thorough search of this site (and others) and haven't found a discussion of how to do this yet. Also, the categories are not ordinal.


title: Speeding up R Loops in an EM algorithm implementation

text: I am writing an implementation of the EM algorithm to estimate parameters in a GLMM in R. My model is $$y_{ij}|\alpha_i \sim Bin(1,p_{ij})\textrm{, where }logit(p_{ij}) = \mu + \alpha_i \textrm{ and } \alpha_i \,\,iid \sim N(0,\sigma^2).$$ In other words, my data has a binary response that can be modeled with a fixed intercept, $\mu$, and a random effect, $\alpha_i$. I am interested in performing EM for the sake of better robustness than I'd get through an ML or REML estimator like the ones used in lmer. My procedure for estimating $\mu$ and $\sigma^2$ goes as follows:   1. Choose starting values for $\mu$ and $\sigma^2$ based on REML estimates.   2. Using Metropolis-Hastings, randomly generate L $\alpha_i$ vectors.    3. Compute the expectation for the log-likelihood for $\mu$ and $\sigma^2$ using the L simulated values of $\alpha_i$ in a pretty standard Monte Carlo scheme.    4. Maximize the expectation in (3) wrt $\mu$ and $\sigma^2$.   5. Use these maximize estimates as new starting values, and return to (1).    6. Repeat k times. This produces two k-length vectors, $\bf{\mu}$ and $\bf{\sigma^2}$. My concern is the following: speed. I have largely vectorized the procedure, but my program chokes takes awhile for values of L larger than 5,000 and k$\geq$100. I have determined that the slowest step is unfortunately Step 6. This could mean that any of the steps 1-5 are slow. I suspect that because of the slowdown with high L, the slowdown occurs during the M-step. My likelihood is separable wrt $\mu$ and $\sigma^2$ (no cross-products), so instead of one multivariate optimization, I am able to get by with two single-variable optimizations, which has sped up the computation. Things I would like:   * Is there a way to compile R code? I call an intensive series of functions thousands of times. Variables are pre-allocated before assignments and functions are all vectorized, so I can't imagine many more gains in speed in those areas. Compiling the code would drastically improve my run time. My code is already written, so I'm very averse to converting it to C or C++.   * I have been using the `optimize()` function on the relevant parts of the log-likelihood to get the appropriate maximizations of $\mu$ and $\sigma^2$ at the given iteration. In MCEM problems in the mixed models framework, is there an algorithm or function that leverages the structure better than a general optimization routine that could give me better run times?


title: compare models using anova in r

text: I got the following outcome. I am confused with the interpretation. When I consider RSS then model 2 is better than model 1. What does that high p value mean? Does in influence my conclusion?? > anova(reg4.3,reg4.4, test="Chisq")               Analysis of Variance Table          Model 1: Gas ~ Flow * CODin + A + B + C     Model 2: Gas ~ Flow * CODin + A + B + C + Blue + Green            Res.Df     RSS Df Sum of Sq Pr(>Chi)     1     60 10034.5          2     58  9851.3  2    183.21   0.5831


title: Probability distribution of Fourier coefficients

text: Suppose I have $N$ independent random numbers with a uniform distribution in $[+1, -1]$. Suppose I take the discrete Fourier transform of these numbers. What would the probability distribution of each frequency coefficient be? Each frequency coefficient is the weighted sum of $N$ independent random numbers with the same distribution. If it were an unweighted sum, then the distribution would presumably tend towards a normal distribution as $N$ increases. But since it's a _weighted_ sum... The only fact I have which seems relevant is that the PDF of the sum of some random variables supposedly equals the convolution of the PDFs of the individual variables. (Assuming they are all independent, obviously.) I'm not sure what you get if you convolve several unequal uniform distributions thought...


title: What are the most useful sources of economics data?

text: When doing research in Economy, one frequently needs to verify theoretical conclusions on real data. What are reliable data sources to use and cite? I am mainly interested in sources that provide various statistical data such as GDP, population, CPI, PPI etc. **EDIT:** Here's an aggregation of the links appearing in this thread + a few more I remembered.      **Generic:**   \- Thomson Reuters Datastream (not free, very comprehensive)   \- World Bank Data   \- United Nations Data   \- IMF Data   \- ADB Data   \- WTO Stats   \- Infochimps \- massive resource of a wide variety of public and private (commercial) datasources - plus their API   \- Freebase (now owned by Google) - open data resource   \- DBpedia \- an approach to using the Wikipedia API   \- Wikipedia API \- or go direct and access Wikipedia direct   \- CIA World Factbook   \- OECD Statistics   \- Wolfram Alpha \- a knowledge search engine   \- Zanran \- a numerical & statistics search engine   \- Inter-university Consortium for Political and Social Research   **National:**   \- UK government data project   \- US government data project   \- US FRED: Federal Reserve Economic Data   \- US Bureau of Labor Statistics   \- US Bureau of Economic Analysis   \- U.S. Census   \- US Health & Retirement Study \- U.S. Centers for Disease Control \- U.S. Bureau of Transportation Statistics \- Eurostat's Databases **Other:**   \- OpenData Initiative at Wikipedia


title: Expected value for discrete (nominal) variable?

text: I'm trying to understand the concept of the expected value. Especially, what bothers me is the expected value for discrete random variables. I will try to formulate it by examples:   * The expected value when throwing a six-sided die is $E(X) = 3.5$, which as I see it, is not even part of the probability space $\Omega = \\{1,2,3,4,5,6\\}$. Why is this not a problem at all? Is here anything wrong with my assumptions about the expected value?   * Is there an expected value when throwing a coin? I cannot calculate a weighted mean on head and tail? Is there an expected value? Is it the maximum likely value, and what if there are multiple maximum likely values? I hope my questions have sufficiently revealed my (mis)conceptions of the topic. I'd appreciate any help and insights!


title: Comparison of two simulated curves for given experimental data

text: Physicist here with a naive question, I've got a data set (points $x_i,y(x_i),i\in \mathbb N$) and two different smooth curves $f(x),g(x)$ from my simulation. I can't quite judge with my eye which one fits better with my eyes only. What are some tools to make a judgement here? My initial thought is to compute the quantities $$\sum_i(f(x_i)-y(x_i))^2,$$ $$\sum_i(g(x_i)-y(x_i))^2,$$ which comes from my familarity with least square fits. Here the situation is that I'm already given all data and curves. > What are standard ways to do this? What are generally interesting quantities > which let me judge the simulated curves?


title: How to simulate random variables according to the law of a pregiven data sample

text: Say I have the following sample:               -0.38848247;0.21655804;-1.08211969;0.00369104;0.0993393;0.25531731;1.40574444;-1.80271115;-0.58780605;-0.35026458;0.67197532;-1.29654652;-0.58874467;-0.4004176;1.46242829;0.83946654;-0.24861179;-1.604154;1.53555232;1.45589014;-0.52092604;1.65504054;0.12854317;1.66771236;1.47964458;-1.55677722;-0.40070632;-0.41957953;0.37803884;0.45349904;1.18265468;1.43979945;-1.55395031;0.01054229;0.48887945;0.12949194;-0.40283111;1.2998402;-1.02205575;-0.34370088;0.22117962;0.52840463;-0.9584811;1.22482249;-1.51189671;0.372597;1.68446854;-0.74440632;1.33920212;0.18072373;-0.35813474;0.39400846;1.18971633;1.02192759;-1.90422461;0.18459334;0.18096905;-1.73870267;1.57349896;-1.05632536;0.1864611;-1.49696658;0.51070568;-2.25007651;-0.64768552;1.80404194;2.25164576;-0.07925576;-0.83550267;-1.65513631;0.25913869;0.36030077;-0.09006407;-1.64359237;-0.1312756;-0.13790883;-0.09940115;-0.02089164;-0.60924589;-0.05451811;2.11736111;-1.40329353;-0.71788744;-0.45888623;-0.75608368;0.45762458;-0.24299548;-0.29224218;-1.0488731;-0.62028903;-0.08257067;0.92297771;0.2964071;-0.02598973;-0.20439059;-0.25195469;0.20572878;-0.49343988;0.25886695;1.20595313      How can I simulate random variables that would have the **same density function** as that of the above sample? (i.e **same law** ) A quick and dirty solution is to perform a uniform random selection out of this sample, but I am interested in more sophisticated methods. Any ides? literature references?


title: Non normal data

text: My model has seven independent variables and one dependent variable. Each are subscales for a 5 point likert scale questionairre. I tested each variable data and they are all non normal. I did spearman correlation and identified several significant correlations. I tried testing the model using linear regression but a very small R square. I am considered a beginner in statistics hence i dont know what other technique i need to use in order to test my model. I have small sample size of 247.Hope anyone can help me. Thanks


title: Is there a standard method to deal with label switching problem in MCMC estimation of mixture models?

text: Label switching (i.e., the posterior distribution is invariant to switching component labels) is a problematic issue when using MCMC to estimate mixture models.   1. Is there a standard (as in widely accepted) methodology to deal with the issue?   2. If there is no standard approach then what are the pros and cons of the leading approaches to solve the label switching problem?


title: What is your favorite layman's explanation for a difficult statistical concept?

text: I really enjoy hearing simple explanations to complex problems. What is your favorite analogy or anecdote that explains a difficult statistical concept? My favorite is Murray's explanation of cointegration using a drunkard and her dog. Murray explains how two random processes (a wandering drunk and her dog, Oliver) can have unit roots but still be related (cointegrated) since their joint first differences are stationary. > The drunk sets out from the bar, about to wander aimlessly in random-walk > fashion. But periodically she intones "Oliver, where are you?", and Oliver > interrupts his aimless wandering to bark. He hears her; she hears him. He > thinks, "Oh, I can't let her get too far off; she'll lock me out." She > thinks, "Oh, I can't let him get too far off; he'll wake me up in the middle > of the night with his barking." Each assesses how far away the other is and > moves to partially close that gap.


title: How to get to a t variable from linear regression

text: I know this must be standard material, but I had difficulty in finding a proof in this form. Let $e$ be a standard white Gaussian vector of size $N$. Let all the other matrices in the following be constant. Let $v = Xy + e$, where $X$ is an $N\times L$ matrix and $y$ is an $N\times 1$ vector, and let $$\left\\{\begin{align} \bar y &= (X^TX)^{-1}X^Tv\\\ \bar e &= v - X\bar y \end{align}\right.\quad.$$ If $c$ is any constant vector, $J = N - \mathrm{rank}(X)$, and $$\left\\{\begin{align} u &= c^T\bar y\\\ s^2 &= \bar e^T\bar ec^T(X^TX)^{-1}c \end{align}\right.\quad,$$ then the random variable defined as $t = u/\sqrt{s^2/J}$ follows a normalized Student's T distribution with J degrees of freedom. I would be grateful if you could provide an outline for its proof.


title: Sampling with inclusion/exclusion criteria from known distribution

text: Is there a simple way to determine (other than by simulation) the expected standard deviation of a sample, when sampling from a population with known distribution, for example a normal distribution with mean 25 and standard deviation 5, but the samples are included/excluded if they fall within certain limits ? In particular, I am interested in the case where the inclusion criteria is that the sampled values are between 25 and 30. Thanks in advance. This is not homework.


title: Training of a Neural Network

text: I am trying to train an Artificial Neural Network for classification. In the input layers, I have 402 neurons; the first 400 are binary, and the last two are floating points in the range -1 to 1. In the hidden layer I have 400 neurons, and in the output layer I have a single node which I want to represent values between -1 and 1. I have tried to train this network using a vectorized implementation of back- propagation which I have found online (I have tried different implementations, and also implemented one myself). My problem is, that my network does not seem to learn a lot. If my learning rate is higher than around 0.0001, I get into trouble, and quickly goes into a local minimum, and with a lower learning rate the learning is (obviously) very very slow. I can as much training data as possible, so this is not a problem, but of course time is limiting, so I would like to be able to train this network in a decent time. Do you have any intuition about what might be wrong, or how much data is needed to train this network of around 160.000 weights? If it is relevant, I can upload some of the data. **Edit:** Due to the comment by Martin, here is some learning statistics for different number of hidden neurons: Google Docs Spreadsheet Another thing which I have observed, is that for my dataset, a constant output of 0.3 will result in a SSE of around 160, so i definitely want to get below this SSE.


title: Locating Streitberg Rohmel paper

text: I'm looking for the paper: _Streitberg and Rohmel (1986) **Exact calculations for permutation and rank tests: an introduction to some recently published algorithms.** Statist. Soft.Newsletter, 12, 10-17._ Google Scholar shows a citation, but not a link, and says that it's cited by 60 other papers. It appears that the Statistical Software Newsletter was discontinued a few years back, or at least moved online to the csdassn.org, but I can't seem to find an archive section. My university library only has the newsletter from 1993 onward. Does anyone have a suggestion of where I might find this paper, or better yet, actually have a copy?


title: Estimating the distribution from data

text: I have a sample of data generated in `R` by `rnorm(50,0,1)`, so the data obviously takes on a normal distribution. However, `R` doesn't "know" this distributional information about the data. Is there a method in `R` that can estimate what kind of distribution my sample comes from? If not, I will use the `shapiro.test` function and proceed that way.


title: Create a function from a data set when there is one independent variable and 3 dependent

text: I have a data set which consists of 1.Number of users, 2.Cpu usage, 3.Ram Usage, 4.Response time, The way I understand it Number of users is my independent variable. Is there a way, a method, or a program to help me create a function that expresses Response time in terms of cpu usage, ram usage and number of usage? I believe regression is the way to go but I cannot find any example using Google. Below there is a graph that show all my values ![enter image description here](http://i.stack.imgur.com/AdPNE.png) Any help would be much appreciated UPDATE: Is it better to assume that only response time is dependent and the other 3 independent? If so what method should I use?


title: Does this graphical model describe a hidden Markov model?

text: I'm facing the problem of visual tracking in computer vision. I have some observation (image blobs by background subtraction) produced by some moving object, and the task is to infer the state (position, velocity) of the object given the observation. I assume that each state depends only on the immediate predecessor (Markov assumption), ruled by a temporal model (the dynamics) that describes how a state $X_t$ evolve in $X_{t+1}$, and that the measurement model do not depends on time, so that $y_t$ do not depends on $X_{t-1}$ but only on $X_t$. Everywhere in the web I have found that this problem can be described by this graphical model: ![enter image description here](http://i.stack.imgur.com/TgxG0.png) For inferring the posterior distribution, after receiving new measurements, a Kalman filter or Particle filter is used. But anywhere I explicitly read that this is an Hidden Markov Model. Is this a HMM? Why? Why not?


title: Which confidence interval should I use to bound actual values around a predicted value?

text: I have a 'black box' model that outputs model predictions, but I don't know what actually goes on inside the model. I also have a bunch of data that I can pass through the model so that I can compare the predictions with actual results. Which interval would I use to bound most of the actual values around the predicted value? I don't think that I would use the prediction interval because that would give me a range that the next prediction is likely to be in, correct? I would like to say that for a single prediction, I am 95% confident that the actual value falls into this interval.


title: Block randomization for different strata, how to achieve 1:1 treatment allocation

text: I am trying to randomize 60 patients in two treatments (A and B) where age will be used as strata. Since we don't know how many patients of each stratum will enter the study I have created two randomization lists of 15 blocks each (4 patients at each block). Since treatment allocation should be 1:1 (thus 30 patients receiving trt A and 30 patients receiving trt B). How should I deal the imaginary scenario: at the end of the study I will have 39 patients of stratum 1 (20 patients received trt A and 19 received trt B) and 21 patients of stratum 2 (11 patients receive trt A, 10 received trt B)? Wouldn't be considered as major deviation from the study design if i will have after unblinding 31 patients receiving trt A and 29 trt B. Will i be safe if i add the phrase ± patient along with my 1:1 ratio?


title: Ridge regression results different in using lm.ridge and glmnet

text: I applied some data to find the best variables solution of regression model using ridge regression in R. I have used `lm.ridge` and `glmnet` (when `alpha=0`), but the results are very different especially when `lambda=0`. It suppose that both parameter estimators have the same values. So, what is the problem here? best regards


title: Data normalization and sufficient statistic

text: I was taught that when we feed our data to machine learning algorithm (e.g. SVM), we should first normalize our data. Suppose I have a set of data $X = \\{x_1,x_2,...,x_n\\}$, I knew two-way of normalizing them, let $\hat{\mu}$ and $\hat{\sigma}^2$ be the sample mean and sample variance of X. I can normalize each data point by $$ y_k = \frac{x_k-\hat{\mu}}{\hat{\sigma}} $$ I think I know this when I first learn PCA. I can also normalize it using the minimum and maximum of the data: $$ y_k = \frac{x-m}{M-m} $$ where $M = \max(X)$ and $m = \min(X)$. By using this normalization, I can make sure the normalized data will be in the interval $[0,1]$. I observe a fact that: ($\hat{\mu}$,$\hat{\sigma}^2$) is sufficient statistic to a normal distribution and the projection matrix of PCA is a solution to a minimization problem that minimize $l^2$-norm. On the other hand, $(M,m)$ is sufficient statistic to a uniform distribution. My questions are:   1. The observation gives me an intuition that what normalization technique I should employ is depends on my belief (or the learning algorithm believes) of underlying distribution of my data. If I believe the distribution is normal distributed, I should use z-score normalization, if I believe the distribution is uniform distributed, I should use min-max normalization. Is my thought correct?   2. If I do not know the underlying distribution, how should I do data normalization?   3. If I am going to feed my data to an online learning algorithm (like winnow algorithm), are there any online data normalization technique?


title: Hypothesis test for difference in preference pre and post treatment

text: Two large beach holiday destinations near Bangkok, Thailand are Pattaya and Hua Hin. In a random sample of 100 individuals in Bangkok who have been on a weekend away to nearby beach destinations 6 preferred Pattaya and 3 preferred Hua Hin. The other 91 did not have a preference. These individuals were then exposed to certain tourism promotional material and after such exposure 12 preferred Pattaya and 1 preferred Hua Hin. The remaining 87 had no preference. I would like to set up a hypothesis test to determine if the promotional material had an effect on preference. Your assistance or a referral to a URL for a tutorial on this kind of hypothesis testing would be greatly appreciated.


title: Error of the sum of complex numbers

text: I'm analyzing the effect of precision error in optics experiments. One of the relevant quantities is a sum of $N$ complex numbers, each with a complex relative standard error $\epsilon$. How can I calculate the standard error of the summation? I've been thinking along the lines of a Gaussian random walk but have not made much progress with the idea.


title: Fitting a multilevel multivariate model in R with `glmer`

text: **Background** I have a large dataset that contains three binary outcomes for individuals belonging to groups. I am interested in jointly modeling these binary outcomes because I have reason to believe they are positively correlated with one another. Most of my data is at the individual level, however I also have some group-level information. Because of the structure of my data, I am treating this as a 3-level logistic regression. The first level defines the multivariate structure through dummy variables that indicate for each outcome. Therefore level-1 accounts for the within-individual measurements. Level-2 provides the between-individuals variances and level-3 gives the between-group variance. _Hypothetical Example:_ Suppose I have data for students in classrooms. I want to examine whether certain student level characteristics are important predictors for passing three different pre-tests (math, history, and gym). The pre-tests are constructed such that about half of the students should pass each exam (no floor or ceiling effect). Since some students are better than others, I expect whether or not they pass their history exam to be correlated to their probability of passing their math and gym pre-tests. I also expect that students in the same classroom will perform more similarly than students across classrooms. **Here is my attempt at writing out the model** I use $h$ to index level-1, $i$ to index level-2, and $j$ to index level-3. Recall that level-1 corresponds to within-individual, level-2 corresponds to between-individuals, and level-3 corresponds to between-groups. So I have $h$ measures for the $i^\text{th}$ individual in the $j^\text{th}$ group. Let \begin{align} \alpha_{1ij} &= 1 \text{ if outcome}_1 = 1 \text{ and 0 otherwise} \\\ \alpha_{2ij} &= 1 \text{ if outcome}_2 = 1 \text{ and 0 otherwise} \\\ \alpha_{3ij} &= 1 \text{ if outcome}_3 = 1 \text{ and 0 otherwise} \end{align} \begin{align} \text{log}\left( \frac{\pi_{hij}}{1 - \pi_{hij}} \right) &= \alpha_{0hij} + \alpha_{1hij}Z_{ij} + \eta_h \\\ Z_{ij} &= \beta_{0ij} + X_{ij}\beta_{ij} + U_{j} + \epsilon_i \\\ U_{j} &= \gamma_{0j} + X_{j}\gamma_{j} \+ \rho_j \end{align} Please leave suggests about this notation; I am not positive that it is correct. **Trying to specify the model in R** I have been using R 3.0.1, but am open to solutions using other standard software (e.g. SAS, Stata, WinBUGs, etc.). Since I am modeling a binary response I am using the `glmer` function in the `lme4` package. My data is in a long format with one row per outcome per individual per group. One of my current problems is correctly specifying a 3-level model. Given 5 individual-level measures and one group level measure, I have tried to specify the model as:               > glmer(pi ~ outcome1:(x1 + x2 + x3 + x4 + u5) +                   outcome2:(x1 + x2 + x3 + x4 + u5) +                  outcome3:(x1 + x2 + x3 + x4 + u5) +                  (1 + x1 + x2 + x3 + x4 | individual) +                  (1 + u5 | group), family=binomial, data=pretest)      but this often produces warnings and does not converge. **My questions**   1. Does my approach make sense? (i.e. does it make sense to model a multilevel multivariate model by specifying the multivariate structure in the first level?)   2. I have difficult with the proper notation and appreciate suggestions for clarifying my notation.   3. Am I using the right tools for this problem? Should I be using other packages or software? Thanks in advance for your suggestions and advice.


title: How to rate successive predictions of the outcome of an event which are made while it is taking place?

text: Consider the scenario of making successive predictions about the outcome of a sporting event while it is taking place. I will use tennis as a concrete example because it has clearly defined moments at which it is sensible to make a new prediction (i.e. after every point is played) but I am also interested in cases such as soccer where predictions can be made with arbitrary frequency (e.g. every five seconds). After every point of the tennis match is played, a model makes a prediction giving the probability that each player will go on to win the match. In a match with $N$ points played, $N$ predictions are made and my question is how should the performance of these $N$ predictions be summarised to give a single value for the performance of the model on this match? An obvious starting point would be the Brier Score, i.e. the average squared error of the predictions. However, taking this approach may not adequately penalise outliers - e.g. a single terrible prediction among otherwise good ones. At the other extreme, the maximum absolute error might favour a consistently mediocre set of predictions over a set of very good predictions and one terrible prediction. Which of these possibilities is more tolerable is open to debate but I seem to have described a general problem with summary statistics. What I am more interested in learning about are approaches that are specifically tailored to my scenario, i.e. which take into account the fact that these predictions are all made on the same outcome and each prediction is made with successively more information.


title: What does function "effects" in R do?

text: I don't understand the explanation in `R`'s help file for effects(): > For a linear model fitted by `lm` or `aov`, the effects are the uncorrelated > single-degree-of-freedom values obtained by projecting the data onto the > successive orthogonal subspaces generated by the QR decomposition during the > fitting process. Can anybody explain what this means? Are the orthogonal subspaces alluded to the one dimensional subspaces spanned by the columns of the Q-part of the QR-decomposition (and hence orthogonal to each other)? Or are they supposed to be orthogonal to something else?


title: Ideal number of variables for PCA analysis

text: I working with a dataset of around 4000 variables. I decided to carry out a PCA analysis for the data, but I am not quite sure about the suitable number of variables I should include in the test. Would feeding a big number (such as 4000) variables interfere with the PCA accuracy? As far as I can understand from the definition and the way the PCA is conducted, the number of variables should not matter, but I am not 100% sure and I couldn't find any source talking about the effect of number of variables on PCA. My question in simple terms would be: should I just include all the variables I have into the PCA test or should I do other tests to reduce the number of variables a little before doing the PCA? Many thanks,


title: Interpreting interaction effect in a Tobit model

text: I have estimated the following model on Stata 13: `tobit y x1 x2 x1x2` where:   * $y$ is aid spending by a charity in a developing country   * $x_1$ is GDP per capita of developing country   * $x_2$ is proportion of official funding that charity receives (so is a fraction between 0 and 1)   * $x_1x_2$ is an interaction term between GDP per capita and proportion of official funding (both of which are continuous variables) Edit: I am interested in the expected value of $y$ (censored and uncensored), i.e. ($E(y|x)$). Therefore I calculated the marginal effect of the interaction term (evaluated at the mean of $x_1$ and $x_2$). Specifically I found ($∂E(y│x)$)/($∂x_1 ∂x_2$). Suppose I find that ($∂E(y│x)$)/($∂x_1 ∂x_2$) is negative. I am inclined to interpret this as: _As the proportion of official funding that the charity receives rises, the charity sends less aid to countries with a higher GDP per capita_. Would this be correct?


title: Guessing test question answers from scores

text: My teacher likes to give online quizzes that are about 20-30 questions long. Every student has the same questions in the same order. We are not told after taking the quiz which questions we got wrong, but the system does tell us our score. That made me curious as to inferring the answers statistically. (Note that this is not _cheating_ , as the teacher encourages heavy collaboration.) My idea is this: when a courageous person is confident of his answers, he submits the quiz. Each answer for each question has a _fitness score_. Each time someone gets the result after taking the quiz, the square of this quiz score is added to the fitness score of every answer he selected. After some tries, we will have a good prediction of the right answers for all of the questions. Is this statistically rigorous? If not, are there better procedures?


title: Conditional expectation of $(X+Y)^2 | X = x$

text: If the density $f(x,y) = c$ when $x>0$, $y>0$, $x+y < 1$ and $0$ otherwise, find $$E((X+Y)^2 | X = x)\text{ for } x \in (0,1).$$ How to approach this question? Can we approach it by letting it equal to $$E( (X+Y)^2 | X = x) = E (E( (X+Y)^2 | X = x))$$ Just not sure...


title: How to investigate properties of the trimmed mean for a Cauchy variable?

text: > Let $X_1,X_2, \dots,X_n$ be a sample from a population with distribution > function $F(x-\theta)$, where $F$ is symmetric around $0$. The $\alpha$ > trimmed mean $T_n(\alpha)=\dfrac{1}{n-2\lfloor n\alpha \rfloor > }\displaystyle\sum_{i=\lfloor n \alpha \rfloor+1}^{n-\lfloor n > \alpha\rfloor}X_{(i)}$ is asymptotically normally distributed, i.e. $$\sqrt > n(T_n (\alpha) -\theta) \sim \mathcal N(0, \sigma^2(\alpha,F))$$ where > $$\sigma^2(\alpha,F)=\dfrac{2\displaystyle\int_0^{F^{-1}(1-\alpha)}t^2~dF(t)+2\alpha(F^{-1}(1-\alpha))^2}{(1-2\alpha)^2}$$ > Propose a method of choosing $\alpha$ from the data.compare the difference > between your proposed adaptive choice of $\alpha$ and its theoretical > counterpart. Investigate your proposed methodology when $F$ is standard > Cauchy. i.e. $f(x)=\dfrac{1}{\pi(1+x^2)},x\in \mathbb{R}.$(Density > function). From this paper and this I understand that how I can find the value of $\alpha$. But then how I can investigate it for Cauchy Distribution. Please help. Thanks in advance. _EDIT:_ According to whuber's comment here I need to show $\sigma^2(\alpha,F)=\frac{1}{2\alpha-1}+\frac{2(\pi \alpha+1)\cot(\pi\alpha)}{\pi(2\alpha-1)^2}$. I understand that here I need to use standard cauchy. I want to know how to find value of $F^{-1}(1-\alpha)$ and $\int_0^{F^{-1}(1-\alpha)}t^2~dF(t)$. Thanks in advance.


title: Interpreting "good figures of merit" in xorshift* context

text: I am trying to implement an xorshift* PRNG (practically an xorshift with a multiplication step in the end) as a long term replacement of Java's Math.Random. I have been reading an article pointed out in an answer to another question, but I am confused with the magical $M_t$ values used for the multiplication. The article comments on this number in the following way: > The goodness of the multiplier, however, is established by a figure of merit > which is a normalized best distance between the hyperplanes of families > covering tuples of length $t$ given by successive outputs of the generators. > The length t is an additional parameter, and $M_{32}$ has the best figures > of merit for $t = 32$. Clearly, if an alternative multiplier provides > improvements on both t and the associated figure of merit, we have a hint > that it could be chosen instead of $M_{32}$. > > Lacking that possibility, what if we scramble xorshift64's output with a > multiplier that has better figures of merit for a lower $t$? We thus also > ran experiment with the multiplier $M_8 = 1181783497276652981$, which has a > better figure of merit for $t = 8$ [L'Ecuyer 1999], and the multiplier $M_2 > = 8372773778140471301$, which has a better figure of merit for $t = 2$ and > was kindly provided by Richard Simard. Seeing as the author has used different multipliers in different xorshift* implementations (different state size) I am not sure how to interpret the role of $t$ in $M_t$. Here's the article by L'Ecuyer that's referenced above.


title: Expected value of sample median given the sample mean

text: Let $Y$ denote the median and let $\bar{X}$ denote the mean, of a random sample of size $n=2k+1$ from a distribution that is $N(\mu,\sigma^2)$. How can I compute $E(Y|\bar{X}=\bar{x})$? Intuitively, because of the normality assumption, it makes sense to claim that $E(Y|\bar{X}=\bar{x})=\bar{x}$ and indeed that is the correct answer. Can that be shown rigorously though? My initial thought was to approach this problem using the conditional normal distribution which is generally a known result. The problem there is that since I do not know the expected value and consequently the variance of the median, I would have to compute those using the $k+1$st order statistic. But that is very complicated and I would rather not go there unless I absolutely have to. Thanks.


title: Expectation-Maximization (EM) method for parameter estimation using fuzzy logic

text: I am sorry if my question is not fit here. If so, please recommend me the correct forum. I am thinking of estimating a fuzzy model using the EM method. I have a set of observations from a nonlinear model. I would like to Takagi- Sugeno fuzzy regression model, where the parameter estimation of the consequent parts are performed using the EM method. The consequent parts consist of linear regression model. Is there a way to do that? Or is there references available in the literature for that purpose? Thanks a lot


title: Processing unlabelled data with ANN

text: It was my understanding that the training data has to be labeled to create a model and then label any unseen data. As in the case of using ANN, how is Google able to identify cats in Youtube videos without labeling the videos or images? I am new to ML and tried to read some literature, but couldn't get the correct explanation. Would appreciate any pointers/explanation in more of layman terms.


title: How to compare the fit of two tables of percentages?

text: I have two tables, you can think of them of results of two models. I would like to compare them to a benchmark to evaluate which model is a better fit. What statistical technique do I need to do this ? Right now, I'm just summing the absolute differences, with the smallest being the best fit, this doesn't seem to be very robust. Any suggestions ? I have included an image of the data, in case my description above was not very good. ![data set](https://dl.dropboxusercontent.com/u/33103477/Untitled.jpg)


title: Can we remove trees from a random forest with poor OOB error to improve generalisation?

text: My objective is to improve out of sample generalization of my random forest _while holding the number of trees constant_. Suppose that I am only allowed to use $n$ trees on the out of sample data but it is computationally feasible to estimate $z$ trees on my in sample data where $z \geq n$ (from which, the best $n$ is selected from the $z$ estimated trees based on OOB error). Is there any research on what we expect the generalisation error to be as $z/n$ increases from $1$ to larger (theoretical asymptotic result or an actual empirical result)? Will it offer better generalization than just choosing the trees randomly (equivalent to setting $n=z$)?


title: Mean length of time spent queueing in $M/E_2/1$ system?

text: **Context:** Consider a $M/E_2/1$ queueing system, where the customer arrival rate is $\lambda$ and the service time distribution has a gamma distribution with parameters $2$ and $\mu$, i.e. with p.d.f. $\mu^2te^{-\mu t}$ , $t ≥ 0$ **Question:** Assuming that the system is in equilibrium, deduce an expression for the mean length of time spent queueing in the system by an arriving customer before she starts being served (the mean "waiting time") **What I have done so far:**   * traffic intensity is: $\rho=\frac{2\lambda}{\mu}$   * Unconditionally at equilibrium, number of customers who arrive during a service time is given by the probability generating function: $K(z) = (\frac{\mu}{\mu+\lambda-2\lambda})^2$   * At equilibrium, the PGF for the number of customers left behind in the system by the $n^{th}$ customer when she leaves is: $\Pi(z)=\frac{\mu(\mu-2\lambda)}{\lambda^2z^2-\lambda(\lambda+2\mu)z+\mu^2}$ But I don't know how to use the above to get what is being asked for. * * * Note 1: this question is related to Mean service time of a $M/E_2/1$ queueing system? Notes 2: I am looking at a more interesting modification to question 5(d) for one of the 2012 RSS exams for which they do no provide a solution: http://www.rss.org.uk/uploadedfiles/userfiles/files/GD3_2012_final%20(web%20version).pdf


title: Sampling methods

text: I have the following situation:   * I'm asking SMEs (small & medium enterprises) how they have reacted to situation ABC.   * Because I don't know who has encountered situation ABC, I ask all SMEs of one country (X).   * Only Y percent respond, because the others hadn't encountered this situation _or_ didn't want to respond. What was my sampling method? I thought it would be _systematic sampling_ because "each element in the population has a known and equal probability of selection". But I'm unsure.


title: Probability that $n$ random points lie on semicircle

text: If random points are chosen from a circle, what is the probability that all of them come from same semicircle ? My reasoning is that; Let two random points be chosen first. Obviously then both of them lie in same semicircle. Now mark any diameter including those two on the same side of it. Any more point chosen will either lie within or beyond the semicircle made by diameter each with probability 1/2. So 3 points lying on same semicircle is $ \frac{1}{2} $. 4points lying in same semicircle is $\frac{1}{2}.\frac{1}{2} = \frac{1}{2^2} $   Similarly, for n such points the total probability is $\frac{1}{2^{n-2}}$. Am I right in my reasoning?


title: Nested logit program in R

text: My question is specific to transportation modelling using a nested logit (NL) model. I wonder how to make a program including a _t_ -parameter, _t_ -value, and likelihood in R. I did estimation in a multinomial logit (MNL) model, and now I am planning to advance it using NL. Can anybody help do this? Notes for my previous estimation:   We have revealed preference (RP) data about mode choice. The modes available are car, minibus, and motorbike. In the survey we only ask the respondent's current mode and alternative mode that they use for commuting. Say we have 400 samples. Because each respondent has a different combination of current mode and alternative mode, we estimate modal choice by employing a dummy mode and using a MNL model instead of binary logit. For this reason the `glm` package doesn't work in this case. In stated preference (SP) data for train, we have 5 experiments for each respondent, thus we have 2000 data. Then we combine RP data and the SP model into a "combined RPSP model" by using the additional variable _µ_ (scale parameter). Variables considered are travel cost, travel time, frequency and income. I don't think the `glm` package for NL works in this case. I'm planning to divide the modes into 2 nests: public transport (minibus and train) and private vehicle (car and motorbike), and of course with additional variable _λ_ (scale parameter) for the nest.


title: Logistic Regression - Model Does not Fit

text: My name is Abhi and I am fairly new to statistics. I found some sample exercises online & I am trying to solve them to get a better understanding of model development. ### Problem Statement Assume a forum like stack exchange or cross validated. Using parameters like number of comments, number of upvotes, reputation points of the user, etc (there are about 10 more fields) can you predict the likelihood of the question getting answered. There about 4000 records to consider ### Current Attempt It seems to me that number of up votes and comments **should** be strong indicators. However when I graph them (and do the chi square significance test) I get very poor results ![enter image description here](http://i.stack.imgur.com/BhHPA.jpg) ![enter image description here](http://i.stack.imgur.com/4O4sG.jpg) number of up votes - p value much less than 0.05   number of comments - p value much less than 0.05 ### My Question What should be my next step from here? Is there something obvious that I am missing? Are there any transformations that I should consider. Any help would be much appreciated ### Update #1 I graphed the distribution for number of comments & number of likes. They are not normally distributed ### Update #2 Based on the suggestions below, I calculated the means and standard deviation of the 2 groups - number of comments for the people whose question was answered (m-2.094,sd-4.008) & number of comments for people whose question was unanswered (m-5.22,sd-5.688). Both of them are within 1 standard deviation of each other. I ran the t test and the difference between the 2 means is 1.91 with a p value <<< 0.0001 Does this mean this feature is useless or do I need to transform this feature


title: Estimating adjusted risk ratios in binary data using Poisson regression

text: I am interested in estimating an adjusted risk ratio, analogous to how one estimates an adjusted odds ratio using logistic regression. Some literature (e.g., this) indicates that using Poisson regression with Huber-White standard errors is a model-based way to do this I have not found literature on how adjusting for continuous covariates affects this. The following simple simulation demonstrates that this issue is not so straightforward:               arr <- function(BLR,RR,p,n,nr,ce)     {        B = rep(0,nr)        for(i in 1:nr){        b <- runif(n)<p         x <- rnorm(n)        pr <- exp( log(BLR) + log(RR)*b + ce*x)        y <- runif(n)<pr        model <- glm(y ~ b + x, family=poisson)        B[i] <- coef(model)[2]        }        return( mean( exp(B), na.rm=TRUE )  )     }          set.seed(1234)     arr(.3, 2, .5, 200, 100, 0)     [1] 1.992103     arr(.3, 2, .5, 200, 100, .1)     [1] 1.980366     arr(.3, 2, .5, 200, 100, 1)     [1] 1.566326       In this case, the true risk ratio is 2, which is recovered reliably when the covariate effect is small. But, when the covariate effect is large, this gets distorted. I assume this arises because the covariate effect can push up against the upper bound (1) and this contaminates the estimation. I have looked but have not found any literature on adjusting for continuous covariates in adjusted risk ratio estimation. I am aware of the following posts on this site:   * Poisson regression to estimate relative risk for binary outcomes   * Poisson regression for binary data but they do not answer my question. Are there any papers on this? Are there any known cautions that should be exercised?


title: Is $\phi(t) = \frac{1}{2}[cos(t) + cos(t\pi)] $ absolutely continuous?

text: X is a real-valued random variable with characteristic function: $\hspace{20mm}$$\phi(t) = \frac{1}{2}[cos(t) + cos(t\pi)$, $\hspace{10mm}-\infty<t<\infty$.   1. Is the distribution of X completely continuous?   2. Is there an r such that $E(X^r)$ does not exist or is infinite?     1. The distribution is not absolutely continuous as the limit of sin as x approaches infinity does not exist.     2. ?


title: check for players undergoing the maximum/minimum change in time course data?

text: I have time series data of many variables (This data comes from high throughput biological experiments, which measure the amount of particular biological components such as RNA content of a cell etc., ) The time series looks like $x \to x_1, x_2, x_3 $ which is the value of $x$ at time points $1,2,3$ respectively $y \to y_1, y_2, y_3 $ etc.. Now, I would like to shortlist the set of variables which have seen the biggest jump or fall (or the ones that did not change much) in their trend. I am interested in understanding which biological players changed radically over the time course and that these may be important players driving biological mechanisms. What is the best measure to calculate this?


title: Reviewing classification approach in Weka on Reuters21578

text: I'm not really familiar with Weka, I learned to use it by watching some tutorials so I am not 100% sure if my approach is the correct one. I have collected the Reuters21758 dataset and I use the documents prescribed in the ModApte split. To make it easy I load training and test instances in Weka (first all training instances followed by all test instances), perform preprocessing and during classification is specify a 75% split. I was wondering if someone can have a look at my .arff file as well as the output of this classification to tell me whether or not I did something wrong. It can be found here: https://www.dropbox.com/s/k78h5dse3hkm973/dataset.zip For clarification, I represent each topic in the Reuters dataset in its binary form {0,1} and train a classifier separately for each topic. In the zip-file in the link is the .arff file and output file of the topic **acq**.


title: Detecting varying time rates

text: I was thinking about detecting the veracity of user ratings by examining the time-stamp. Basically I have a series of times $(t_1, t_2, \ldots, t_n)$. So assuming that ratings should normally come in at a constant rate, is there a good way to detect whether some time-stamps are unusually clumped up, because a set of consecutive fake ratings was injected over a short period of time?


title: Can I use the Cholesky-method for generating correlated random variables with given mean?

text: I want to generate correlated random variables with a given correlation matrix, means, and variances. Does the Cholesky decomposition only work when the initial random variables are iids with the same mean and variance?


title: What to input into Wilcoxon rank sum test when three pre and three post measurements are compared?

text: I have measured 23 patients three times pre and three times post treatment. Now I would like to see if the treatment had a significant effect. I thought about using a Wilcoxon rank sum test for this (because I don't know how the characteristic is distributed in the population), which takes two vectors of variables, for example in R:               x <- c(subject.1.pre, subject.2.pre, subject.3.pre)     y <- c(subject.1.post, subject.2.post, subject.3.post)     wilcox.test(x, y, paired = TRUE)      But unlike this textbook example, I don't have one pre and one post measurement, but three pre and three post measurements. Now I am confused: Should each vector (x and y) be the means (pre versus post) for   1. each patient,    2. each repeated measurement, or   3. the individual values Example:               Patient Pre1 Pre2 Pre3 Post1 Post2 Post3     1       a    b    c    d     e     f     2       u    v    w    p     q     r        1. x = {mean(a, b, c), mean(u, v, w)}   2. x = {mean(a, u), mean(b, v), mean(c, w)}   3. x = {a, b, c, u, v, w} Depending on the input, a Wilcoxon rank sum test returns different values:   1. V = 155, p-value = 0.6221   2. V = 6, p-value = 0.25   3. V = 1144, p-value = 0.4867 Which one is correct? There are a couple of NAs in my data. This is the original data, if you want to play with it (times 1 to 3 are pre, 24 to 26 are post treatment; 4 to 23 [missing here] are measurements during treatment; aph is the duration of a behavior in minutes):               subject time aph     1 1 2.500000     2 1 3.789474     3 1 1.514563     4 1 1.111111     5 1 5.846868     6 1 4.588235     7 1 5.600000     8 1 5.066667     9 1 11.647059     10 1 2.000000     11 1 NA     12 1 4.461538     13 1 18.000000     14 1 1.058824     15 1 9.217391     16 1 27.840000     17 1 15.375000     18 1 NA     19 1 6.000000     20 1 9.714286     21 1 12.484848     22 1 16.503497     23 1 20.666667     1 2 3.500000     2 2 4.658824     3 2 4.750000     4 2 1.894737     5 2 4.000000     6 2 2.800000     7 2 14.228571     8 2 11.000000     9 2 NA     10 2 2.666667     11 2 3.764706     12 2 4.705882     13 2 13.272727     14 2 2.000000     15 2 18.444444     16 2 17.555556     17 2 14.222222     18 2 2.000000     19 2 4.000000     20 2 8.461538     21 2 19.200000     22 2 13.902439     23 2 13.000000     1 3 3.000000     2 3 NA     3 3 7.360000     4 3 2.588235     5 3 1.611374     6 3 1.500000     7 3 3.365854     8 3 22.375000     9 3 10.838710     10 3 2.923077     11 3 3.488372     12 3 5.176471     13 3 37.666667     14 3 1.176471     15 3 7.454545     16 3 36.235294     17 3 6.823529     18 3 2.222222     19 3 6.133333     20 3 11.428571     21 3 42.705882     22 3 28.105263     23 3 18.333333     1  24 5.375000     2  24 7.434783     3  24 9.935484     4  24 4.941176     5  24 3.732057     6  24 2.250000     7  24 13.777778     8  24 4.740741     9  24 7.428571     10  24 3.619048     11  24 7.135135     12  24 5.285714     13  24 5.000000     14  24 2.588235     15  24 11.500000     16  24 27.125000     17  24 2.875000     18  24 2.810811     19  24 9.428571     20  24 7.058824     21  24 3.833333     22  24 16.869565     23  24 16.769231     1  25 10.133333     2  25 5.454545     3  25 7.875000     4  25 6.782609     5  25 1.688312     6  25 2.666667     7  25 11.785714     8  25 3.333333     9  25 6.333333     10  25 7.400000     11  25 12.588235     12  25 4.000000     13  25 6.727273     14  25 3.030303     15  25 6.000000     16  25 30.666667     17  25 3.750000     18  25 3.000000     19  25 8.000000     20  25 8.000000     21  25 2.571429     22  25 10.695652     23  25 14.266667     1  26 7.750000     2  26 2.511628     3  26 6.333333     4  26 7.225806     5  26 3.281250     6  26 1.888889     7  26 10.400000     8  26 4.875000     9  26 13.714286     10  26 8.500000     11  26 NA     12  26 4.875000     13  26 8.181818     14  26 1.250000     15  26 4.380952     16  26 27.176471     17  26 2.375000     18  26 2.941176     19  26 11.714286     20  26 6.000000     21  26 2.375000     22  26 14.763780     23  26 14.313253


title: Difference-in-difference regression using R

text: Sorry for my bad English, could you suggest R-code for the implementation of difference-in-difference regression? I don't understand how many coefficients I need. In my analysis I compare the effect of a new law on the stock exchange volume, I have 2 periods and 2 samples. Thanks a lot! Thanks Jeremy, but in the Imbens's model the regression contains a series of multiplications between variables and dummies. I used:               lr1 <- lm(VOLUME ~ DUMMYCAP + DUMMYTIME + DUMMYCAP*DUMMYTIME )       but the result of regression is not significant. This question is about code and its economic significance. I'm sorry if this post is not appropriate for the site. But if someone has made a test of this type and knows how to write the code I would be grateful


title: Metrics of datasets?

text: It is quite common to discuss what metrics to be used for a model/algorithm/method. But when it comes to metrics of datasets, it is much less discussed. I wonder, what are the general metrics for datasets? That is to say, how does one quantitively measure if a dataset is good or not for its purpose? I understand, there may not be too many general metrics, given the purpose and nature of datasets differ vastly. So, more specifically, I am interested in the follow few cases: what are the metrics for 1) text datasets, 2) classification datasets, 3) unbalanced classification datasets? Any knowledge, insights, interesting ideas are welcome and appreciated.


title: The means are statistically equal

text: What does it mean "The means are statistically equals?" I ask this question because I'm studing ANOVA (`anova()` function in R) and I know that this test checks if the means of N groups are statistically equals. My doubt concerns what "statistically equals" means because if I have two series like:               A <- c(2,4,12,14)      the mean is **8** But if I do               B <- c(-20,-10, 10, 52)      the mean is always **8** but the number are totally differents. So the question is: does it check if the means are similar and stop OR check if the means are similar AND also check if the number that "generate" this mean are similar? Thank you


title: Sample size for cluster analysis

text: We have recently conducted a cluster analysis for an open card sorting task in which 19 experts grouped 112 items each. We achieved decent silhouette values for the cluster analysis but have been asked by reviewers to provide evidence that the number of experts we used was sufficient for this task (I think they are expecting something kind of power analysis). We have used Kendall's Coefficient of Concordance to examine the concordance in individual sort patterns between participants (this came out as .25), but it seems they want something more. After looking around I cannot find anything which might suggest a suitable sample size or even that sample size is not important and I'm pretty sure power analyses don't apply for cluster analysis. Is there anything I can report to show my sample size is sufficient (or not if the case may be)? I am using R to analyse the data. BTW by sample size I mean number of participants not number of items included in the cluster analysis.


title: Plotting only highest lines of an Excel Chart

text: I'm currently building a chart with Excel out of a fair number of lines of statistics. Since there are too many lines, i cannot plot them all, otherwise the chart would become unreadable. I'm currently selecting the lines and values by hand, which is a tedious process. I'm trying to automate some of this work. **The question** : I want to only plot the 5 highest lines of statistics. At any point in time (the "columns"), I can easily extract the 5th highest value. I then know if a value is among the 5 best ones ("IF" function). For a cell value to not be plotted, it's enough to keep the cell empty. So, if a value is not among the 5 best ones, i want the cell to be considered empty, so that it is not plotted. The problem is, i don't know how to "not output anything" or "output something which will not be plotted" from an "IF" function. For example, if output value is "", nothing is displayed in the cell, but it is still considered like a zero by the chart.


title: Cross-validation or bootstrapping to evaluate classification performance?

text: What is the most appropriate sampling method to evaluate the performance of a classifier on a particular data set and compare it with other classifiers? Cross-validation seems to be standard practice, but I've read that methods such as .632 bootstrap are a better choice. As a follow-up: Does the choice of performance metric affect the answer (if I use AUC instead of accuracy)? My ultimate goal is to be able to say with some confidence that one machine learning method is superior to another for a particular dataset.


title: Testing for differences between a sample and a subsample

text: For reasons I won't go into, I want to test for group differences between an entire sample A ($n=65$) and a subsample B ($n=45$). All of the variables I want use for mean comparisons are continuous. Is there an appropriate test for this? This doesn't seem to meet the sample independence condition for the t-test or Mann-Whitney U. Is it better to just test for group differences using a t-test between B (n=45) and !B (n=20)?


title: Getting a second-order polynomial trend line from a set of data

text: Alright, so I have about a thousand datapoints that I'm plotting on a chart (scatter plot). Here's a few of the records:               x              y     2.426032708   10     9.449923509   24     4.409997771   24     -3.125392294  11     13.04820475   26      When looking at the full data on a chart, I can visually see a 2nd order polynomial trendline is the way to go. With my full set of data, how can I calculate the formula of the line? Everything I've found in google has told me how to do it in Excel, and that's not what I want. I want to get the formula so I can predict results at other points. (Please be gentle. I'm a programmer, not a stats guy. If you need anything else, don't hesitate to ask)


title: Is there any identification assumption for IV?

text: I am asking this because I know there are IA's for cross sectional estimators and D-in-D estimators, but am unsure if there are any for IV estimation. Does anyone know if there is or isn't, and if there is - what would it be?


title: Highly Collinear Independent Variables of Interest

text: Suppose I am interested in the follows: I have county-level data. For each county, I know the share of the population that was born from one parent and the share of the population that was born from two. Suppose I am interested in the impact of the black population on tax revenues. In particular, suppose I am also interested in the impact of the share of the population born from two black parents separately from the share born from one black parent. I am using OLS. The problem is that the share of the population born from one and two black parents is highly correlated, resulting in a near multicollinearity problem. Of course I can proceed by just testing the total share of the black population, and that is a result in which I am inherently interested. However, if I include the 1 and 2 parents as together as two regressors, the results seem meaningless relative to the total share because the 1 and 2 parent variables produce a near multicollinearity. Should I run two separate regression? A regression on the share with 1 parent and a regression on the share with 2 parents? Or is there nothing that can be done because of the near multicollinearity? Thank you very much for any help you can provide.


title: conditional mean

text: Unfortunately nobody seem to know the answer to my first question...does anyone know how to compute a conditional expectation on absolute value? Let $$\boldsymbol y = \begin{bmatrix} \boldsymbol y_{a}^{\top} \\\ \boldsymbol y_{b}^{\top} \end{bmatrix} $$ $$\boldsymbol y \sim \mathcal{N}({\boldsymbol 0},{\Sigma_{y}})$$ where $$ \Sigma_{y} = \begin{bmatrix} \Sigma_{aa} & \Sigma_{ab}\\\ \Sigma_{ba} & \Sigma_{bb} \end{bmatrix}$$ we know that $E[y_a|y_b]=\Sigma_{ab}\Sigma_{bb}^{-1} \boldsymbol y_b $ does anyone know how to compute the following $E[y_a|abs(y_b)]$ ? That would be very helpful for me. Thanks a lot!


title: Sample analysis

text: I have a regression with typical Stata output. The only independent variable is a dummy (=1 if person bought certain item). How can I find how many people in the sample bought the item?


title: Interaction in regression

text: ![enter image description here](http://i.stack.imgur.com/zAuqM.jpg)In my regression model of BMR on age (yrs) , gender (1=male, 0=female), height and weight. The regression equation obtained is: BMR = 1232.059+13.281 weight+3.954 age+192.214 male-3.471 age:gender. Please interpret the negative coefficient of age:gender interaction. Do the negative coefficient have anything to do with the positive coefficient of age and gender?


title: Determine the significance of factor variables

text: I have a df like this:                     Rate   Gender             7          Male             5          Female             10         Male             20         Male             9          Female      Etc I need to run an analysis to determine significance of gender on the rate. I'm really new to this, any ideas how could I start this?


title: a guide to using the phi coefficient

text: This one is straightforward. I am trying to find a good documentation of the phi correlation. There is a Chi-Square Statistic to test the significance of the the correlation coefficient, I am hoping that the document would have good documentation of this, how many degrees of freedom for the critical value and all these things. Thank you.


title: Predicting outcomes where there's no data in rpart

text: I have a data set with about ~2000 data points. Of these ~1000 actually have features/data. (All 2000 data points have an outcome) Where there's no data, there _is_ very likely a signal. In other words, if I'm not able to get data for that data point, that's significant in and of itself, and there is a big difference in the target variable when there is data and when there's not data. The problem is that rpart shows me that these missing values result in "observations deleted due to missingness," so the model doesn't make any Is there a way to set up an rpart model such that it makes a prediction on data points with no features (i.e. the only variable available is the outcome)?               my_model <- rpart(bad ~ count_found, data = mydata)


title: Variance of the sum of a Poisson-distributed random number of (normally distributed) random variables

text: Suppose I have a fruit stand where I sell oranges. Every day I have an average of $\lambda$ people (Poisson distributed) each buying a normally distributed weight quantity of oranges with mean $\mu$ and standard deviation $\sigma$. What's the mean and variance of the number of pounds of oranges I sell per day? OK, the mean is trivial: $\lambda \cdot \mu$. But how do I calculate the variance? Strangely enough I've come up with an answer through trial and error that seems to fit really well with my Monte Carlo simulations, but I can't convince myself that it's actually the correct answer in theory. If someone can point me to the "right" answer, I'll post the one I got empirically for comparison.


title: A model for multiple nested proportions?

text: I think this is more of a stats advice than code question so I've posted here rather than SO. My data looks like this (~300 more sample rows).               > head(nest)       GroupA GroupB SubA SubB Treat     1     45    100    3    6   Pos     2     43     80    2    5   Neg     3     64     70    0    4   Neg     4     75     76    2    5   Pos     5     53     65    0    3   Pos     6     54     83    2    3   Pos      **GroupA** and **GroupB** are counts of different things within each sample. **SubA** and **SubB** are a counts of a small subset of the same things (i.e. nested). **Treat** show two different types of sample _Pos_ and _Neg_. My hypothesis is that the proportion of the subset ( **SubA/SubB** ) is a greater than the proportions of the whole ( **GroupA/GroupB** ) - particularly in _Pos_ sample types. This is suggested by plotting. However I find that if I naively just treat each row separately and apply a `binom.test` for each then it's way underpowered (too few **SubA** , **SubB** ). Given this and quite a few 0 counts in **SubA** ... Is there a more statistically powerful way I can address this idea? Particularly how might I model this in R? Any pointers are very welcome. Thanks.


title: how to do conditional forecasting with cointegration model?

text: I'm confused about multistep forecasting from VECM model for 2 cointegrated series. The model is pretty simple, in error-correction form: $$ \Delta x_{t+1} = \alpha_1 (y_t - \beta x_t -\beta_0) +\varepsilon_{1,t}, \\\ \Delta y_{t+1} = \alpha_2 (y_t - \beta x_t -\beta_0) +\varepsilon_{2,t}. $$ Where $y_t - \beta x_t -\beta_0$ is the long-term cointegrating relation and $\alpha$'s are short-term error correction coefficients. Now, what I need is multistep forecast of $y_{t+k}$, _conditional_ on $x_{t+k},\dots,x_{t}$. So $x$'s are known (including future) and I need forecast $y$'s $$\hat y_{t+k} = \mathbf{E}(y_{t+k} | x_{t+k},\dots,x_t).$$ As usual $y_{t+k}$ doesn't depend on next values of $x$, only on past. First thought - Ok, just use long-run cointegration relation $$\hat y_{t+k} = \beta \hat x_{t+k} + \beta_0.$$ Here $\hat y, \hat x$ are conditional expectations. But this doesn't account for short-term error-correction dynamics and residuals $(y_t - \hat y_t)$ from these forecasts have long'ish periods of deviations from 0 - i.e. residuals have autocorrelation. Would be good to improve forecasts using short-term error-correction as well. I try to do that using VECM formulation above. I illustrate this for 1-step forecast, but of course it's easy to iterate formulas and obtain k-step forecasts. If I use only the second equation of the model, this gives: $$ \hat y_{t+1} = (1+\alpha_2)\hat y_t - \alpha_2 (\beta \hat x_t +\beta_0) $$ and doesn't look good because forecast for $y_{t+1}$ doesn't depend on $x_{t+1}$ - we don't use all available information for forecasting. If I use both equations of the model and substitute one to the other, we have: $$\Delta \hat y_{t+1} = {\alpha_2 \over \alpha_1} \Delta \hat x_{t+1}.$$ Here $\hat y_{t+1}$ forecast depends on $\hat x_{t+1}$, but now it is only determined as short-term relation between differences via error-correction coefficients! So now long-term cointegration relation doesn't comе into play. So the question is - how to derive conditional forecast for $y$ that uses both cointegration and error-correction dynamics (I guess the latter will affect only several first steps predictions) and also conditions on all information available?


title: Estimators for linear regression when multicollinearity is present

text: I have a multicollinearity problem is a linear regression model and ridge regression was suggested as a solution. So I have spent quite some time researching different ridge regressors in the literature (there's at least a dozen), however in the course of my research I have found that Principal Components Regression can be considered a special case of ridge regression and so I am including it too. Clearly Partial Least Squares regression is closely linked to Principal Component Regression so it seems it should be added too (although I haven't researched that yet). Also when deriving Ridge regression from the Bayesian viewpoint it became obvious that a Kalman filter approach could be used to implement the Bayesian approach so I'm including that too. The Bayes approach seems fairly natural in my application as I have several years worth of data some in which the collinearity effects are present some where they are not so I can build up a decent prior for the Beta parameters for the cases where little of no MC is present. I have reviewed various comparison papers and they usually compare the estimators on the basis of their Mean Square Error Performance. Generally in the Ridge v Bayes/Kalman their seems to be no clear winner with the performance depending on the orientation of the betas being estimated to the principal components of the design matrix, the level of multicollinearity and the signal to noise ratio. So it would seem that the best estimator can only be defined in the sense of the best for a given problem. When should I use lasso vs ridge? Now through my research I have stumbled into Lasso (and Least Angle regression) and from the answers to the previous thread (see above) I can add "Elastic Net" and "non-negative garrote" (and probably more). I just wonder does it make sense to add these to my comparison too? If I'm totally honest it seems that I should? The theory behind the Lasso is pretty similar to that used for the Ridge estimator with the optimization being done in the L1 rather than L2 space. Ideally I would like to compare all the different classes of shrinkage estimator using the best estimator from each class. Right now I have at least a dozen ridge regressors and no way a priori of knowing which would be best for my problem (given the inconclusive results of the comparison papers I mentioned previously). If I only have to add say the "Lasso" and "Elastic-net" and even the "Non-negative Garrotte" then probably I can do it but if they have as many variations as the Ridge literature has and similar problems identifying which one is likely to be best for a given problem then it seems that the whole thing may become rather unwieldy I can't cover every approach that's ever been suggested in the literature. So my questions are:   1. What are the different classes of estimators to deal with Multicollinearity in linear regression: eg. Ridge, PCR, PLS, Lasso...   2. Apart from Ridge regression do these other classes have one implementation that is generally regarded as best. Any paper I have seem that compares Ridge to other classes of estimators they generally use the basic Ridge regression method which is shown to be one of the poorest performing ridge methods in the papers I have seen.   3. Do you have experience on comparing Ridge to these other techniques and which was found to be best?   4. Is it realistic to compare these different classes all together? Here is a list of the lasso techniques I have found from a review article by Tibshirani (2011): Grouped Lasso; Elastic Net Lasso; Fusd Lasso; Adaptive Lasso; Graphical Lasso; Dantzig selector; Near Isotonic Regulation; Matrix Completion; Compressive Sensing; Multivariate Methods; Now my problem consists of only 4 different X values and around 30 y values. The really high multicollinearity occurs between X variables 3 & 4\. Although there can be significant multicollinearity with variable 2 as well. (Variable 1 is a constant). Most of these techniques are suited for p>>N but this seems directly at odds with the comment from Gary here: When should I use lasso vs ridge? That lasso should be used when you have high multicollinearity effects and few variables? In any case given my problem I think that the Lasso and Elastic Net are most suitable. Can anyone shed any light as to whether any of the remaining techniques may be helpful for my problem?


title: Finding correlation of discrete distribution data

text: I'm not good at Statistics but I need to use it to summarize my case study. So I consult with a person who're good at Statistics. He suggested that My data is discrete because all the values are integer and some are duplicated as shown below. I want to find correlation between x and y of the two data sets below. In data set 1, x and y are discrete. In data set 2, x is discrete but y is continuous. **Data Set 1**                x   y     -5  -19     -5  -21     -5  -23     -5  -23     -4  -17     -4  -17     -4  -18     -4  -19     -4  -19     -4  -23     -4  -19     -4  -18     -3  -13     -3  -20     -3  -17     -3  -15     -3  -17     -3  -18     -3  -18     -3  -18     -3  -18     -3  -18     -3  -18     -3  -15     -3  -15     -2  -15     -2  -15     -2  -14     -2  -14     -2  -11     -2  -12     -1  1     0   0     0   0     0   0     0   0     0   0     0   0     0   0      **Data Set 2**               1   287.92     1   340     1   219.62     1   30     1   30     1   292.56     -2  -70.32     -2  -82.04     -2  -68.53     -2  -68.53     -2  -51.89     -2  -57.36     -3  -94.01     -3  -66.61     -3  -95.91     -3  -95.91     -3  -95.91     -3  -95.91     -3  -95.91     -3  -95.91     -3  -78.14     -3  -125.34     -3  -78.14     -4  -110.45     -4  -94.01     -4  -94.01     -4  -162.85     -4  -162.85     -4  -162.85     -4  -120.93     -5  -307.67     -5  -123.19     -5  -371.38     -5  -166.91     -5  -162.85     -5  -302.21     -5  -176.42     -5  -307.67     -5  -371.38     -6  -390.71     -6  -387.74      so he suggested me to find covariance and then find the correlation from covariance **question 1** : Should I use _sample_ covariance or _population_ covariance Its difference is that in sample covariance, the average of X and Y in the sample are used; whereas in population covariance, ux and uy are used. **question 2** : some books mentioned the equation of population covariance that pop_cov = sigma(xi-ux)(yi-uy) / N others said pop_cov = sigma(xi-ux)(yi-uy) = sigma(XY)- (uxuy) Which one is correct? **question 3** : if the answer of question 1 is population covariance, how can I get the probability of occuring each value of X (p(x)), each value of Y(p(y)), and occuring a particular X and Y (p(x,y)). Can I calculate these 3 probabilities from the sample I have, for example p(5) is 4/39 (4 is a number of 5 and 39 is total number of X). **question 4** : the formula of correlation need standard deviation (corr = covariance/(sdx sdy)), will we use standard deviation of population or sample ( is it depending on the type of covariance we use, of pop or sample)? Thank you very much. I 'm sorry that I am so novice about Statictiscs


title: Transform continuous variable to discrete variable

text: I am having an input data set( 14 x 250 ) which has attached a labels set ( 1 x 250 ). The labels set is discrete( 0 vs. 1). The problem is that the each one of the 14 characteristics are continuous. For example, for characteristic 1 I am having :               10.18 ; 11.18 ; 9.89 ; 22.45     0     ; 1     ; 0    ; 0      I want to transform the continuous variables to discrete variables, and while I am doing this I want to have a Mutual Information value between the transformed characteristics and the label values, as high as possible. Also, I want to keep the number of variables in which I split each characteristic as low as possible, because I want to use them for a Bayesian Network model. Could you recommend me some transform algorithms/techniques for this? Or some matlab build in functions. I know that I can apply the mean and standard deviation for each characteristic and then build the labels based on it, but I am not sure whether the mutual information value will be high.


title: Simple linear regression -- through a specific point

text: In follow up to a question posted here: Constrained linear regression through a specified point What are some ways that we can tell if a simple linear regression model is "good" when we constrain it through a specified point?


title: Trouble understanding posterior probabilities

text: I'm having trouble understanding what this problem is asking for. Could someone clarify how to obtain posterior probabilities using backward/forward probabilities? ![Problem](http://i.stack.imgur.com/TPsfe.png)


title: Grouped data analysis in SPSS

text: I have a question about how to analyze some questionnaire data. I have three different participant groups: A,B,C and each group people did the same questionnaire. So I have a total score for that questionnaire for each participant. I am trying to find the correlation between these scores with another number (like how many medicines they take during a day) and what the differences are between the three groups of participants. What test should I do, how can I analyze it in SPSS? I am a newbie in statistics. Thanks in advance!


title: How to develop a forecast for counts where the counts per period is low?

text: I'm trying to analyze for any trends in what I consider to be low observations or counts. Let's say you have a bus service with the following commuters:               year   average count     2001   12     2002   15     2003   17     2004   13     2005   18     2006   12     2007   9     2008   12      The above numbers are fictitious but similar to what my problem is. How can I analyze the above data for any trends and perhaps attempt to forecast for 2009 to 2012 values?


title: How to define the multiplier range for variance test based outliers detection algorithm?

text: I have a variance test based outliers detection algorithm. The algorithm is exposed with a visual application where the user can configure its parameters namely the multiplier. The question is what is the range of values for the algorithm multiplier? the algorithm basically is calculating the mean and standard deviation (sigma) of the data set and then comparing the dataset elements to the upper and lower bounds to flag an element as an outlier or not.               multiplier <- 2.3;      upper_bound <- mean + multiplier * sigma;      lower_bound <- mean - multiplier * sigma;       the evolution of the two function upper_bound (red) and lower_bound(blue) with the multiplier values is as follows ![enter image description here](http://i.stack.imgur.com/Gzcm6.png) But this graph doesn't give any clue on what range to define. Do you have any idea how to define this range?


title: How do I test whether a regression performs significantly better than another on the same data?

text: I have 2 different regression algorithms trained on the same dataset. I want to know if on a test set one of them performs significantly better than the other -- or at least if they are significantly different. What test do I use ? I am looking for something similar to McNemar's test that exists for classification algorithm.


title: Is this statement about credible intervals from Wikipedia correct?

text: I read this statement on the Wikipedia article about Credible Intervals and it seems dubious to me. Although it has a reference associated with it, the book is not available to read online. > It is possible to frame the choice of a credible interval within decision > theory and, in that context, an optimal interval will always be a highest > probability density set. > >> Wikipedia: Credible Intervals If it's correct, can someone explain how? Maybe I'm just misunderstanding one of the terms.


title: the derivation of the conditional posterior for the Poisson model setting

text: When discussing the Poisson process with changing point (Carlin, Gelfand and Smith, 1992), the model is assumed as ![enter image description here](http://i.stack.imgur.com/FCGgE.png) I am not quite clear about the derivation of conditional posterior distribution based on the above model setting, especially the following part circled with red color ![enter image description here](http://i.stack.imgur.com/WzWyK.png)


title: SVM - combining binary and continous representation of the same feature?

text: How would this influence the accuracy of the SVM model? Let's suppose that I have one variable which max value is 100 and minimum is 0. Currently, I send it to SVM as a single continuous feature, for example, if the value is 8 I send it as 0.08. I am wondering though, what would happen if I were to represent this variable as a set of 100 dummy variables? For example if the variable has a value of 56, I set the feature number 56 to 1 and the rest stays at 0. Would this work worse than a one continuous feature? How about combining these two approaches? Would that improve the accuracy or lower it? Or maybe it would make no difference?


title: Simulation for power analysis for higher order repeated measures ANOVA

text: I am attempting to find a way to perform a power analysis for a higher order repeated measures ANOVA where all factors are within-subjects (i.e., there are no between-subjects factors). I have looked extensively for a way to do this using existing tools like SPSS, or G*Power and other packages but have not found anything that seems appropriate. G*Power does not seem to support power analysis for repeated measures designs where there are multiple within- subjects factors. I am wondering if R is the right path to take to perform a power analysis of this type. I am just starting with R and was wondering if anyone had a script built for this? I found a script at the link below that seems promising but wanted to confirm that it is indeed for a power analysis of a design with ONLY within-subjects factors. The design I am interested in analyzing would have 2 factors, typically one with 2 levels and another with 3 or more. However, I also need to assess a much larger design with 4 factors. If the script in the thread linked below is for a totally within-subjects design, can it be expanded to include more factors? If so how? Power analysis of repeated measures ANOVA using simulation in R?