title: Alternative ways to define markov chain states?

text: I have a time series of daily measurement data, if you plot the histogram of the data points of the time series, you see long tail and high positive scewness. Also zero values is an important state. My goal is to use Markov Chain to model the time series, but the measurement data are integers without upper bound, so in order to model it as markov chain, I need to define the a mapping from the measurement data to state space. There's no obvious definition of different states of the time series data. My first attempt is to use 25 percentile, 50 percentile and 75 percentile as threshold to define the states(0s are a separate state), say if I'm aiming at 5 states. Is there other ways to define meaningful states?


title: Fractional polynomial model not converging in Stata

text: I am having trouble getting a fractional polynomial logistic model to converge using Stata 12.               fracpoly, logistic: y x      I have tried   * centering the dependent variable   * transforming the variable    * adding a fixed value to every observation so that there aren't any 'zeros' I am trying to fit a 2nd order fractional polynomial using the default powers (-2, -1, -.5, 0, .5, 1, 2, 3). I have experimented, and found that the problem seems to occur with particular powers (-2, 2, 3) and so had been trying with just the remainder. However, I recently added another year's worth of data, and with this larger dataset, the model fails to converge even with this reduced set. The variable looks like this               Variable: x     Obs:147778     Mean: 1.608095     SD: 1.572316     Min: .0792616     Max: 51.26371         I have left the model running overnight, and it still fails to converge. If I `set trace on` then Stata is 'stuck' at               version 11: mata: Mopt_maxmin()      Here is a running line smoother of the data if that helps. ![running_line_smoother](http://i.stack.imgur.com/PfS8a.png)


title: Poisson Distribution of a radioactive decay

text: You have received a radioactive mass that is claimed to have a mean decay rate of at least 1 particle per second. If the mean decay rate is less than 1 per second, you may return the product for a refund. Let X be the number of decay events counted in 10 seconds. Part a. If the mean decay rate is exactly 1 per second, what is P(X <= 8)? I got 0.33282. Part b. Based on the answer to part a, if the mean decay rate is 1 particle per second, would eight events in 10 seconds be an unusually small number? I'm a little confused on part b because isn't it really obvious that yeah we should return it because the mean for particles per second is 8/10 =0.8 when it should equal 1.0?


title: Linear regression on a sample spanning many orders of magnitude

text: Beer's law from chemistry says that the absorbance of a liquid $A$ is proportional to the concentration $C$, so: $$A = kC$$ The standard thing to do then is to prepare a set of solutions with known concentrations, measure the absorbance to form a 'standard curve' (a calibration curve basically), and do a simple linear regression on that data to get the proportionality (which can then be used to predict the concentrations of unknown solutions). One easy way to do this is to start with a known concentration and perform a serial dilution, which would get you 2x dilution, 4x, 8x, 16x....etc. I.e. if you start with a solution of $100\mu \mathrm g/\mathrm {mL}$ you'd get solutions with $50\mu \mathrm g/\mathrm {mL}$, $25\mu \mathrm g/\mathrm {mL}$, $12.5\mu \mathrm g/\mathrm {mL}$ etc... Now when you do the linear regression, you have a data set with lots of data points at low concentrations and very few at higher concentrations. It seems much more natural to represent this problem on a log scale. My question is then, should I be doing a linear regression of $A$ vs. $C$ or $\log A$ vs. $\log C$? When I compare the models, they seem to give answers that are on the same order of magnitude but on the order of 30% different. ![data plotted on linear scale](http://i.stack.imgur.com/MuHNw.png) ![data plotted on log log scale](http://i.stack.imgur.com/sqTg4.png)


title: Marginal model versus random-effects model – how to choose between them? An advice for a layman

text: In searching for any info about **marginal model** and **random-effects model** , and how to choose between them, I have found some info but it was more-or-less mathematical abstract explanation (like for example here: http://stats.stackexchange.com/a/68753/38080). Somewhere I have found that there were observed substantial differences between a parameter estimates between these two methods/models (http://www.biomedcentral.com/1471-2288/2/15/), however the opposite was wrote by Zuur et al. (2009, p. 116; http://link.springer.com/book/10.1007%2F978-0-387-87458-6). Marginal model (generalized estimating equation approach) brings population-averaged parameters, while outputs from random-effects model (generalized linear mixed model) take into account random effect – subject (Verbeke et al. 2010, pp. 49–52; http://link.springer.com/chapter/10.1007/0-387-28980-1_16). I would like to see some layman-like explanation of these models illustrated on some model (real-life) examples in language familiar to non-statistician and non-mathematician. In detail, I would like to know: **When should be used marginal model and when should be used random-effects model? For which scientific questions are these models suitable?** **How should be outputs from these models interpreted?**


title: Inter-rater reliability using Intra-class correlation with ratings for multiple objects on multiple properties

text: ### Initial Question: I am trying to calculate inter-rater reliability. Previous researchers in this area have used intraclass correlation. SPSS has options for two-way random, mixed random and one-way random models. SPSS help says choose the right one based on whether the 'people effects are random' and 'item effects are random', can anybody please explain to me what these terms mean? Also I wanted to know how to structure my data file - I think I need to put each raters' rating of each item in a separate column (rather than each rater on a separate row), is this right? ### UPDATED question: Thanks to the help I have received, I have visited the websites, which were useful. I have decided that I need to use the mixed model, and I have transposed my files so that each rater is a column. ### Additivity: The websites also mentioned some assumptions, one of which is _additivity_ , tested with Tukey's test of non-additivity. Apparently, this tests the hypothesis that there is no multiplicative interaction between the cases and the items. I think maybe I have too many variables in my calculation, because the examples others use have a lot less, e.g. each rater making only one rating, as described below, my raters have made many more ratings than that.   * What does the additivity assumption mean?   * What should be done when the test of non-additivity is statistically significant? ### Inter-rater reliability with multiple objects and properties I will explain in more detail what I have done. I have five raters and they have rated facial and vocal expressions of participants watching different emotional films. Raters make 18 ratings per film, and there are 16 films, so there is a total of 288 variables per rater, per participant rated. Each rater has rated the same 4 participants each, and I have placed these participants into different files. I thought about making seperate files for each of the 16 films, but that would be so many files.   * How should I calculate ICC and organise my files in this specific situation?   * Should I calculate a mean of individual ICCs or perhaps select a random sample?


title: $R^2$ for mixed models with multiple fixed and random effects

text: I have a mixed model with three random effects and many fixed effects. I set aside some of my data to validate my final model and now I would like to compute the $R^2$ on the validation set, in `R`. I found some code for a mixed model with one random effect but not for a model with more than one random. Does anyone know a package or some code to compute $R^2$ for mixed models with more than one random effect? If this isn't possible, does anyone have a good suggestion for comparing model fits?


title: Capture-recapture sampling valid in literary analysis?

text: So I've been compiling a list of fictional (anime/manga) characters which meet a certain criteria (http://www.gwern.net/hafu#list) from a universe of all anime/manga characters since 1963 (which is of unknown total size - but very large!), and I've been wondering how I could estimate how complete my list is at any point. My final goal is to slice the character data and look for trends by year or decade; getting an idea of how large my sample actually is may help me estimate systematic bias. If I've gotten a large fraction of estimated characters, then I can hope that any decade trends may be real and not just a case of 'looking under the lamppost'. What sort of techniques would be useful here? I've looked, and "capture- recapture" seems like an answer (especially reading Predicting total number of bugs based on number of bugs revealed by each tester ) But I'm not clear whether it really applies here. What would a "recapture" be, in my scenario? Would I have to keep track of every Google search or Google Alerts result or list of _hafu_ characters I run into, and write down whether each entry was already present in my master list? Are there assumptions in capture-recapture that aren't satisfied here? UPDATE: since no one replied yea or nay, I went ahead and did it:   * results: http://www.gwern.net/hafu#corpus-estimation-result   * data & source code: http://www.gwern.net/hafu#capture-recapture-code


title: How to understand optimal Scaling in R: The Package homals for novices

text: Does anyone know of a step-by-step guide for the practical implementation of Gifi Methods for Optimal Scaling in R: The Package homals? Although I have an OK theoretical understanding (thanks chl for directing me to articles), I'm a tech novice and some of the language/technical information is a bit too complicated for me. On that note, can I do the same thing within SPSS??


title: Stratified Random Sampling of Ratio (X/Y) where X itself is an estimate

text: I am trying work out the variance of the population estimate var (R) where R = X/Y (X = sum(x) and Y = sum(y)). For each member of the population I know y and I have stratified the population and taken a random sample of each strata. For each member of the sample I need to estimate x. x is not known with certainty and the estimate of x has its own standard error (which I can estimate from a pre sample calibration exercise). However, which I know enough to know that this error should be taken into account my knowledge of sample statistics is limited to the situation where x is a quantity that can be exactly determined. If anyone can help with an answer or point me to anappropriate reference source it would be appreciated.


title: Does an exact test always yield a higher P value then an approximated test?

text: I ran a simulation on that for mcnemar test, and the answer seemed to be yes. I was wondering if this can always be said to be the case that the exact P value is higher (or not smaller) then the p value which is arrived to through an approximation. some code for example:               set.seed(234)     n <- 100 # number of total subjects     P <- numeric(100)     P_exact <- numeric(100)     for(i in 1:100)     {     x = table(sample(1:2, n, T), sample(1:2, n, T))     P[i] <- mcnemar.test(x, correct = F)$p.v     P_exact[i] <- binom.test(x[2,1],x[1,2]+x[2,1])$p.valu     }          #for different n - the level of problem is worse     #plot(jitter(P,0,.01), P_exact )     plot(P, P_exact )     abline(0,1)


title: Over enthusiastic MANOVA in R

text: In recent research I pursue a MANOVA model on a dataset comprised of 50 measured characters and several grouping factors that I need to test. When I use something like this               summary(manova(malesM ~ popMales*manageMales*biomeMales), test = "Wilks")       I get 2.2e-16 significance for every possible output (also with other tests beside Wilks). The same thing I noticed when I use dummy, randomly generated factors. My data comes from a real-life study and I really can not belive that almost any theoretical, grouping factor will be significant. Any advice on how to refine my model, and to test this further?


title: Scaling data that are on different scales for a graph

text: Looking at the following dataset:                Date        Visits   Carts      carts       Orders                           Created   converted    Created     2011-11-11    12277     161        9          36       2011-11-12    11871     93         5          19         2011-11-13    13072     107        8          8          2011-11-14    13594     112        4          34         2011-11-15    12741     129        8          43         2011-11-16    15491     261        16         57      2011-11-17    13418     186        17         42          I've been asked to plot this on a graph, using the Date has the X-Axis and the rest of the data on the Y-Axis. The problem is that the scale of the data is dramatically different. where Visits are in the thousands and Orders Created are in the low tens, the data doesn't plot well on a graph. I was wondering what a statistician would do in this scenario, I could divide the the Visits by a 1000 and then put in the description (Visits (K)), but then I start to have the same problem with Carts Created, as they are in the hundreds and everything else is in the low tens. What is the "done" thing in this scenario? Kind regards, Mike


title: Is it possible to calculate 95% confidence limits of a predicted value of a regression equation when we don't have the raw data?

text: Reading a paper, I realized that a polynomial regression equation was published with its associated standard error and sample size. That regression relates enamel thickness ($y$) and enamel formation time ($x$). I only know the regression equation, not the raw data. This is the regression equation: $$y=8.7+0.37x-0.00005x^2\\\ (R^2=0.97, \text{standard error}= 0.01, n=20)$$ I have my own data from my sample, and I would like to use this regression and know the confidence interval. For example, I have one value ($x=1210 \text{ microns}$), and its formation time following the regression is 383 days. I read in several places that the raw data are necessary to calculate the confidence interval for a predicted value. Is there any way to know it when we don't have the raw data? If so, how?


title: Evaluation measure of clustering (without having truth labels)

text: I'm clustering a set of data but I don't have truth document that allow me to evaluate the result of clustering (I have unlabelled data), so I can not use an external evaluation measure. In this case, is there any efficient evaluation measures, that will allow me to incrementally evaluate the clustering result each time I cluster some new data to the current clustering ?


title: How to do significance testing on averages of regression coefficients across multiple groups?

text: I have over 1000 respondents each with their own regression coefficients. I am analyzing data by different segments by averaging the constant and coefficient scores for each group. Sometimes I will have 2 groups (gender, yes/no questions, etc.) and sometimes I will have multiple groups (Ethnicity, age groups, etc.) What significance test should be used to see if the average coefficient scores are significantly different from each other?


title: FactorGraph approximations: Binary variables to unary

text: I am reading the Expectation Propagation chapter in Chris Bishop's book and there is a bit on approximating factor graphs. So, the original factor graph and the approximating factor graph looks as follows: ![enter image description here](http://i.stack.imgur.com/wkAVG.png) So, the original joint distribution can be written as: $$ p(x) = f_a(x_1, x_2) \; f_b(x_2, x_3) \; f_c(x_2, x_4) $$ Now, the assumption is that each of these factors also factorise according to the variables they depend on, so we have an approximating distribution as: $$ q(x) = \widetilde{f}_{a1}(x_1) \; \widetilde{f}_{a2}(x_2) \; \widetilde{f}_{b2}(x_2) \; \widetilde{f}_{c2}(x_2) \; \widetilde{f}_{b3}(x_3) \; \widetilde{f}_{c4}(x_4) $$ My question is when we do this approximation, how do we initialise these new intermediate factors? So, for example imagine that in the original graph, $f_a$ which depends on $x_1$ and $x_2$ represents some sort of probability as $P(x_2|x_1) P(x_1)$. Now assume $P(x_1)$ is modelled as a gamma distribution and $P(x_2|x_1)$ is modelled as a multivariate normal. When I approximate the graph and introduce the factors $\widetilde{f}_{a1}$ and $\widetilde{f}_{a2}$, what is the nature of these factors? Is $\widetilde{f}_{a1}$ going to be a gamma distribution and $\widetilde{f}_{a2}$ be a MVN? I know that this should be the case because it makes it easy to multiply messages flowing into a variable node as long as they are of the same type but cannot understand whether we are free to give these messages whatever form we want? Also, why is this approximation more easier to solve than the previous case? I can see each factor is much simpler but what does it mean in terms of problem complexity?


title: The right way to report random effects in a Cox survival model

text: I have a data set (patients (15000) nested in hospitals (50)), with a number of covariates. I built a proportional hazards model in Stata, and entered site as a frailty term. I would like to report how survival varies between hospitals. I can extract the frailties and report these, but this has no straightforward, face-value interpretation. I have experimented with the following:   * Extract the (estimated) baseline hazard, and multiply this for each site and plot the different curves. There is no way to quantify the certainty of the frailty with this, and it seems to use a lot of ink to describe the variation in one parameter. ![enter image description here](http://i.stack.imgur.com/gM6L3.jpg)   * Report a median (or other percentile) survival for each site, and do so ideally using something like a funnel plot. I have done this by extracting the (estimated) baseline hazard after centering the covariates, and then multiplying this by the random effect, and then finding the greatest reported survival time below my chosen centile. For now, I have plotted this against the sample size (for that hospital). There are clearly problems with bounding on the y-axis, and I do not know how I might construct control limits. ![enter image description here](http://i.stack.imgur.com/PBVm0.png)   * This reference recommends Silcocks P. Hazard ratio funnel plots for survival comparisons. J Epidemiol Community Health. 2009;63:856-861., and uses a fixed-effects model, eventually reporting relative, centred hazard ratios. Does anyone have any recommendations or comments? NB: I am using Stata if this is relevant.


title: Predictive model on categorical variable with high cardinality

text: So I have a categorical variable with ~ 8000 types (or levels, or cardinality). The dataset is rather large (a couple million records). Most machine learning techniques are failing do to memory limits on my computer, so I think an approach may be to run a loop and do individual regressions for each category step by step. Is there an easy way to do this in R? So for example, if the category was States in America, I would like 1 regression model for Virginia, 1 regression model for Texas, 1 regression model for Vermont, ect. Any other ideas for dealing with a very large categorical variable with a very large cardinality?


title: An incremental Gaussian mixture model

text: **Question 1:** Suppose that data is modelled by a mixture of K probability distributions which are actually Gaussians. $P(x_i|\theta_j)$ is the probability density of the j'th cluster, for which the parameters that we want to learn are $\theta_j = \\{\mu_j, \sigma_j\\}$. The probability that data-point $x_i$ belongs to cluster $c_j$ (i.e. to the j'th mixture) is $P(x_i \in c_j | x_i, \theta_j) = \frac{P(x_i|\theta_j)}{\sum_l P(x_i|\theta_l)}$. Thus the algorithm consists of take some initial values for the parameters of each mixture $\theta_j$ (how do we choose this initial parameter values ?). Then we repeat the flowing two steps many times:   * Compute $P(x_i \in c_j | x_i, \theta_j)$ for each data-point i = 0,..,N, in a mixture j = 1,..,K. Thus we will get a matrix M of N x K.   * Optimize the parameters $\theta_j$ of each cluster j. (how do we optimize/update them ? does this mean just trying with other parameter values ?) I'm confused on how to choose the parameters of each cluster and how to optimize/update them (the second step). So if I understand, we should repeat the following code algorithm until convergence:               X = {x1, ..., xn} set of data-points     M = {m1, ..., mk} set of k means choosen randomly from X     S = {s1, ..., sk} set of k covariance matrix, each si is choosen to be the identity matrix          for i from 1 to n     {         for j from 1 to k         {             compute Pij = pdf(xi | mj, sj) / SUM_l pfd(xi | ml, sl)         }     }          for i from 1 to n     {         for j from 1 to k         {             mj = mj + (Pij * xi) / SUM_l Plj // update the mean         }     }      **Question 2:** Is there an incremental version of a Gaussian mixture model which is incremental (the data-points are processed sequentially)? That is, without "repeating" the 2 steps of the EM algorithm "until convergence".


title: How to plot results from text mining (e.g. classification or clustering)?

text: In text classification and clustering, the number of features are normally big, e.g. I currently get are around 5,000 features which is already really small compared to many other text mining tasks. Given I am totally new to visualization, I have no clue about how should I plot the results from text classification and clustering. For example, in one clustering task, I have three categories, hundreds of documents, and some thousands of features in these texts. After I apply a cluster, e.g. the KMeans cluster from scikit-learn, I have no idea about how to plot these results onto a 2D figure. It is the same for classification, say, I get texts sorted into three categories, but I get no idea how to plot them. I've tried to learn from some examples, but most ones I found are based on pure data with very few features rather than text. **Question:** Could you refer me to any tutorial or paper on either drawing text clustering/classification results or explaining the math part of such visualization. Or any help would be appreciated.


title: How can I sample from a log transformed distribution using uniform distribution?

text: I am transforming an unscaled density function to log scale to avoid underflow issues. BI was performing integration on this function on a grid of values before I used the log transormation, to build a grid based cumulative distribution function. Then, using a draw from uniform[0,1], I was choosing the the largest point on the grid which had a cumulative probability value smaller than the draw from the uniform. This worked fine as long as the univariate density function could be integrated. With the transform to log-scale, I can't really get my head around this mechanism. The joint likelihood is so small that I can't back transform the density, so I have to perform the same uniform distribution based sampling on the log scale. Is this an established practice? Feedback and pointers would be appreciated.


title: Machine learning approach to a selection problem

text: I'm thinking how to tackle such a problem. Let's say we have a set of products, each have some numeric characteristics `(x1, x2, x3)`. A customer is given a choice of `n` such products and buys one. Only a few of products are shown and the choice is made only between them based solely on the characteristics. The number of products for selection varies. Now assuming that I have a good amount of data describing the sets given for choice and the choice made, how do I build a model to predict the outcome of similar experiment in the future? That is, we give `n` products and want to find which will be the choice based on their characteristics. My first thought was neural network, but that would work only if I have a constant set of inputs and here it varies. Any suggestions are welcome.


title: Off-diagonal range for guaranteed positive definiteness?

text: Consider the problem of randomly generating a $p \times p$ covariance matrix $\Sigma$ with the diagonal entries constrained to be 1, and the off-diagonal entries $\Sigma_{ij}=\Sigma_{ji} \sim Unif[-a,a]$. What is the maximum value of $a$, as a function of $p$, such that $\Sigma$ is positive definite with probability 1?


title: Fitting piecewise linear curves with lmer

text: (cross-posted from the R-sig ME list - maybe this is too basic for it) Hello I would like to fit piecewise linear growth curves to a number of subjects (each with a number of observations at different abscissa), where all the slopes for a common interval (e.g, 2 <= x < 3) presumable arise from the same normal distribution, and therefore should be "shrunk" towards each other. What is the best way to about doing this in lme4? I have a toy, cleaned up (although I am not very good with R) example of the type of data I am looking at. Thank you for your time               set.seed(1)     # our growth functions are piecewise linear with these breaks     breaks <- seq(0,5)     mean.slopes <- seq(5,2,-.6)          # generate the true slope parameters for 20 subjects (variance isn't         really constant or known, just simplifying)...     true.slopes <- sapply(mean.slopes,function(x) rnorm(n=20,sd=3,mean=x))          # ...and the intercepts - they really don't come from any particular         distribution     true.intercepts <- runif(n=20,min=30,max=50)          # generate sample data - we don't necessarily have the same number of     observations per subject     mat <- matrix(nr=20*10,nc=3)     for (i in 1:20)     {        xs <- runif(n=10,min=0,max=6)        slopes <- true.slopes[i,]        ys <- true.intercepts[i] + rnorm(n=10,mean=0,sd=0.9) +        sapply(xs, function(x) sum(slopes[1:floor(x)]) + (x - floor(x)) *        slopes[ceiling(x)])        id <- rep(i,10)        mat[(1+10*(i-1)):(10*i),] <- cbind(xs,ys,id)     }     mat <- as.data.frame(mat)     names(mat) <- c("x","y","id")     # qplot(x=x,y=y,colour=as.factor(id),data=mat)      Eventually I might like to shrink the slopes for each individual towards each other, in the sense that if for a given subject the slope in the interval [2,3) is very small compared to the slopes of other subjects in the same interval, then the same would hold for the slope in [1,2), but I would like to get the basics first.


title: Suggestions for improving the fit of probability density

text: I am comparing the error between wind turbine production and predictions from 24 hours ahead. The graphs below are using the normal distribution. But can I do better? ![enter image description here](http://i.stack.imgur.com/oNqHC.jpg) ![enter image description here](http://i.stack.imgur.com/tvsoC.jpg)


title: Probability in random selection - Exposing corruption in a Sweepstake

text: We are running a sweepstake for the football world cup and the first three people to choose "drew" the past three winners while in a closed room. Just for fun, I'd like to show the probability of this happening by chance? As an extension, what is the probability that each of these 3 got the country that they got.


title: Wald test and Likelihood ratio test, where do the confidence intervals on the regression coefficients come from?

text: So I'm trying to build my own Wald test and likelihood ratio test code within a machine learning pipeline. I can get the final fitted logistic regression coefficients from liblinear. I'm coding in MATLAB. How would I get the variance of the regression coefficients and also the confidence intervals of the regression coefficients? Clearly for a variance and confidence interval, you need a sample of multiple sets of coefficients. But I thought you only get a single set of coefficients at the end of the log likelihood optimization. Basically trying to replicate the results in the following link http://www.ats.ucla.edu/stat/mult_pkg/faq/general/nested_tests.htm


title: Find correlation between 2 non-ordered variables

text: I have two variables A={some values} and B={some values). Each value in A is not related to the corresponding value in B. How I can find the correlation between these two variables? or how can I know if these two variables have any kind of association "relationship" between them? Just an example to make the question clear: Assume A refers to the number of daily deaths in a very big hospital (30 values within a month) and B refers to the daily number of absent doctors in that hospital (30 values within a month). Values in A and B are not sorted by the date (i.e. A[i] is not related to B[i]), so they are in different order and we want to know if they are correlated or have any kind of relationship. I assume correlation coefficients will not work hear as each value in A does not match the corresponding value in B.


title: In a binary classification soft-margin SVM problem, is there a case where there are no bounded support vectors for one class?

text: For two class labels +1 and -1, is there a case in conventional soft-margin SVM where one class does not have a margin support vector and just margin error support vectors or correct points?


title: Evaluating Expected Shortfall

text: I am writing my thesis on VaR and ES risk measurements and have encountered some issues with how to best test the accuracy of ES statistics. My understanding of the topic is that backtesting ES adequately is extremely difficult or even impossible, as it is not an elicitable risk measure. And I also assume the number of approaches in the literature are scarcity because of exactly this property. However very recently D.Tasche et.al. (http://arxiv.org/pdf/1312.1645v2.pdf) uploaded a paper where it was argued that ES is not directly elicitable, but indirectly elicitable because it can be approximated by several VaR estimates. Thus, they state ES can be backtested reasonably by backtesting several VaRs (at different confidence levels) related to the ES estimate. Their specific proposition is:$E{S_\gamma }(L) = {1 \over {1 - \gamma }}\int_\gamma ^1 q u(L)du \approx {1 \over 4}\left( {{{\rm{q}}_\gamma }{\rm{(L) + }}{{\rm{q}}_{0.75\gamma + 0.25}}(L) + {{\rm{q}}_{0.5\gamma + 0.5}}(L) + {{\rm{q}}_{0.25\gamma + 0.75}}(L)} \right)$ where L is the loss distribution and gamma the chosen confidence level. From my simple calculations, proxying financial data with a student-t(6) the accumulated VaRs with this approach seems to deviate much from the analytical ES. I wonder if anyone find the Tasche et.al. approach appealing/adequate for backtesting ES? Further I am also interested in knowing what the "state of the art" approach is in the industry (if any particular) for backtesting ES. Any help would be greatly appreciated


title: CDF for uncorrelated bivariate normal

text: I have an uncorrelated bivariate normal process: $X \sim N(0, \sigma_1), Y \sim N(0, \sigma_2), \rho = 0$, and I'm interested in the distribution of $Z = \sqrt{X^2 + Y^2}$. I know that if I'm willing to normalize _X_ and _Y_ by their variance then _Z_ has the chi distribution with 2 degrees of freedom. But I am looking for closed-form expressions of the variance and CDF of Z that incorporate the difference in variance in the two dimensions. I.e., without factoring out the sigmas, and knowing $\sigma_1 \neq \sigma_2$,   * What is the variance of Z?   * What is the CDF -- i.e., Pr($Z \leq \alpha$)? (I know this is also an elliptic distribution, but the general form for those is far too complicated.)


title: What have to be normally distributed: groups or whole sample?

text: It is a general and maybe stupid question, but it will help me to avoid a mistake. To test the dependence of continuous variable on nominal with parametric test (t-test, ANOVA), "the data has to be normally distributed" (distribution have to be not much violated from normal). But I don't clearly understand: I have to test the normality of the continuous data together (whole sample) or the normality of data in every single group? And if my groups have small sizes (n=4-6), and therefore I fail to reject null, does this indicate a permission to use parametric methods? Can anyone provide me a link to some _good_ examples of statistical (biomedical) hypothesis testing with R (google did not helped)?


title: Entropy-based methods in R

text: I was wondering if anybody knows of the existence of an R package which implements entropy-based methods (maximum empirical likelihood, maximum exponential empirical likelihood, minimum discrepancy estimators, etc.)? Thanks


title: Nelson-Aalen estimator in discrete time

text: I was wondering whether the Nelson-Aalen estimator can be applied to discrete- time multi-state data (in months). I ask this because it seems that literature tends to present it only in the continuous-time framework (in particular within the so called "counting processes" approach). Since in discrete time data are incomplete (in the sense that they are interval-censored), I guess that it is necessary to assume non-informative interval censoring, in addition to the random right-censoring assumption. Is this rigth? Can you suggest some literature about it? Thanks for any help


title: Getting the good prediction standard error with R and/or SAS

text: Consider for instance the balanced one-way random effects ANOVA model: $$(y_{ij} \mid \mu_i) \sim_{\text{iid}} {\cal N}(\mu_i, \sigma^2_w), \quad j=1,\ldots,J, \qquad \mu_i \sim {\cal N}(\mu, \sigma^2_b), \quad i=1,\ldots,I. $$ In this post I asked how to get a prediction interval and I found an answer in a published paper by Lin & Liao. The approach is similar to the case of the simple Gaussian sample model and it runs as follows. Denote by $\sigma_{\textrm{tot}}^2=\sigma^2_b+\sigma^2_w$ the total variance and by $\tau^2=\sigma^2_b+\frac{\sigma^2_w}{J}$ the variance of the sample mean $\bar{y}_{i\bullet}$. Then, considering a new observation $y^{\textrm{new}}$ one has $$ y^{\textrm{new}} - \bar{y}_{\bullet\bullet} \sim {\cal N}\left(0, \sigma^2_{\textrm{tot}} + \frac{\tau^2}{I}\right).$$ The variance can also be written as $$ \left(1+\frac{1}{I}\right)\tau^2 + \left(1-\frac{1}{J}\right)\sigma^2_w$$ and it is estimated without bias by $\frac{I+1}{I(I-1)}\frac{SS_b}{J} + \frac{SS_w}{IJ}$ where $SS_b$ is the between sum of squares and $SS_w$ is the within sum of squares. Moreover $SS_w$, $SS_b$, $y^{\textrm{new}}$ and $\bar{y}_{\bullet\bullet}$ are independent. Then we have the distributional approximation $$\frac{y^{\textrm{new}} - \bar{y}_{\bullet\bullet}}{\sqrt{\frac{I+1}{I(I-1)}\frac{SS_b}{J} + \frac{SS_w}{IJ}}} \approx \mathrm{t}_\nu,$$ where $\nu$ are the Satterthwaite degrees of freedom. Now I wonder whether it is possible to get the prediction standard error (or alternatively an approximation) $se=\sqrt{\frac{I+1}{I(I-1)}\frac{SS_b}{J} + \frac{SS_w}{IJ}}$ with R and/or SAS. I am not only interested in this example, but also in more general mixed models. Here I asked whether it is possible with Frank Harrell's `Gls` package but the author himself replies that this package does not handle predictions for a model "without predictor". Below is a numerical example for the balanced one-way random effect ANOVA:               ranovapred <- function(y, group, conf=0.95){         group <- factor(group)         means <- aggregate(y~group, FUN=mean)$y  # groups means         I <- length(levels(group))         J <- length(y)/I         sizes <- table(group) # groups sizes         if(!all(as.numeric(sizes)==J)){ stop("balanced only!") }         ssw <- crossprod(y-rep(means, times=sizes))  # within sum of squares         ssb <- J*crossprod(means-mean(y)) # beween sum of squares         a <- (1/J*(1+1/I))/(I-1)         b <- 1/I/J          v <- a*ssb+b*ssw # estimates the variance of (Ynew-Ybar)         nu <- v^2/((a*ssb)^2/(I-1)+(b*ssw)^2/I/(J-1)) # Satterthwaite degrees of freedom         alpha <- 1-conf         bounds <- mean(y) + c(-1,1)*sqrt(v)*qt(1-alpha/2,nu)     return(list(bounds=bounds, std.error=sqrt(v), ss=c(ssw=ssw,ssb=ssb)) )     }     # fictive data     set.seed(421)     I <- 3     J <- 4      dd <- data.frame(y=rpois(I*J,10), group=gl(I,J))          > dd         y group     1  12     1     2  11     1     3  13     1     4  14     1     5  10     2     6   8     2     7   8     2     8  12     2     9   8     3     10  6     3     11  6     3     12  9     3     > ranovapred(y=dd$y,dd$group)     $bounds     [1] -1.136421 20.636421          $std.error              [,1]     [1,] 3.338538          $ss       ssw   ssb      22.75 55.50


title: Visualize survival analysis with time dependent covariates

text: As a follow-up to Model suggestion for a Cox regression with time dependent covariates is the Kaplan Meier plot accounting for the time dependent nature of pregnancies. In other words the dataset is now broken down into a long dataset with multiple rows according to number of pregnancies. The KM graph, and also the extended cox model, seems to hint at a beneficial effect of pregnancy on outcome. However, looking at the KM graph I'm wondering-should the line for first birth start at 1.0? Wouldnt it be more intuitive to start this line at the y-value for 0 births at x equal to when the first birth is given? EDIT: After researching this closer I see that regular KM is not good. Rather I should use the method of Simon and Makuch which is used in Stata (Simon R, Makuch RW. A non-parametric graphical representation of the relationship between survival and the occurrence of an event: application to responder versus non-responder bias. Statistics in Medicine 1984; 3: 35-44. ) Has anyone seen this implemented for R? ![enter image description here](http://i.stack.imgur.com/1z1Ct.png)


title: Difference between unrestricted VECM and restricted VECM?

text: What's the difference between an **unrestricted** and a **restricted** VECM? I believe a hint lies within the `cajorls()`[1] function of R language's `urca`[2] package.


title: Sample write up for a repeated measures MANOVA using one IV and two DVs?

text: I have run a repeated measures MANOVA for my dissertation study in clinical psychology but am struggling a little bit with the proper write up. I have found several examples on the internet and through my school library for written results, but none specifically for a RM MANOVA with one IV and only two DV's (for pretest and post test results). An ideal example (I can dream, right?) would show the output and the corresponding write up. I don't believe I have written my results correctly. This is my draft so far: > In looking at within-subject effects, we look to see if five weeks of > training compared to five weeks without training has an influence on post > test scores, and find that there was a statistically significant effect > within groups F(2,15)=31.57, p< .001. > > In looking at a pairwise comparison of the interaction between time and > group (meaning that a change over time differs between the two groups from > pretest to post-test measures), there was a significant interaction effect, > with the immediate start group showing significant improvement, > F(2,15)=24.48, p<. 001 and the wait list control group showing no change p < > .001. Therefore, the first hypothesis was accepted, working memory training > did have a significant effect on working memory for the immediate start > group. > > In looking at between-subject effects, we first look to see if there was a > difference in post-test scores between the immediate start group and the > wait list control group, and find that there was a non significant > difference between groups F(2,15)=.514, p= .608. Therefore, the second > hypothesis that working memory training had an effect on attention and > impulsivity is not confirmed, and this hypothesis is rejected.


title: How to compute the CI of an Inverse Gaussian distribution?

text: I have come across some formula for calculating the confidence interval (CI) of data from an Inverse Gaussian distribution. But, being really new in stats, I am unable to figure out what to do. On page 51 of this document, there is a formula. But I am not sure how to compute the different elements of the formula. Can anyone kindly help me out ? Thanks. (In the picture below, while I know that $\hat{\mu}$ means the mean and $n$ means the number of observation. I am not sure what to make of $k$ or $\lambda$... .) ![enter image description here](http://i.stack.imgur.com/rO0Gp.png)


title: A 5pt likert scale and a 6 stage algorithm analysis correlate will ordinal regression analysis work if not than what

text: the two likert scales will determine disease perception A.severity B. susceptibility. The alogorithm will determine stage the of intention to be screened for disease which the participants can fall into 1 of 6 stages. I am not sure of the best way to conduct my analysis I was leading toward ordinal regressional analysis .????? Tiffiny


title: How to generate normal variates subject to mixed constraints?

text: I want to randomly generate 1000 normal variates (using `rnorm`, _e.g._ ) that have mean 100. 25% of the 1000 numbers should be over 110. How can I do this in `R`? I only got this far:               x <- rnorm(1000,100,1)


title: When regression is run on population

text: I have a small population and it is possible to estimate Population regression function from the data. The size of the population is very small. It is only eight. I think since the regression relates to poulation, there is no need to conduct any testing. Am I correct?


title: Unable to figure out right transformation

text: I have obtained some data about how complexity of Java open source projects varies with time. I want to fit a curve to the data, however I am unable to figure out the right kind of transformation. I tried some that came immediately to my head: $\ 1/y \sim x,\ \log(y) \sim x,\ y \sim 1 + x + x^2$, etc. but I never got a good $R^2$ value. Can you please suggest how do I proceed? Below is the plot: ![enter image description here](http://i.stack.imgur.com/Gs4C8.jpg) I have also shared the data here.


title: How to make adjustment for correlation coefficient?

text: Suppose that there are two two-dimensional maps. For simplification, let's say one map is the temperatures in the 48 continental U. S. states, and the other map is the corresponding humidities. The data of temperature and humidity are collected at one specific moment at some equally-spaced (e.g., 1000) locations. If I directly calculate the correlation between the two maps, the coefficient would be inflated because the maps are spatially correlated? Suppose that I know the full width at half maximum (FWHM) that characterizes the spatial correlation, is there a way that I can adjust the correlation coefficient or its degrees of freedom?


title: Modeling Future Logistic Regression Covariates

text: I have a logistic regression model which implicitly estimates the probability that a business will default in each of 8 future quarters conditional upon survival to that quarter. In other words, I have 8 sets of coefficients to give me 8 different conditional probabilities based on the covariates at t=0. (cumulative for each borrower calculated as 1-(1-Pr_conditional_t1)*(1-Pr_conditional_t2)*...(1-Pr_conditional_t8)) Now, I am building a different yet related model, and I want to explicitly estimate the future covariates for the borrowers (same covariates as the first model) as opposed to using the same covariates at t=0 for all 8 quarters. I want to use macroeconomic factors to model these future covariates. I know that Logistic Regression is more intended for use with covariates that have been observed or empirically estimated, but I was wondering if anyone had any suggestions on ways to go about modeling this without creating biases in future observations based on the convex nature of the logit curve and different risk categories, etc... I was thinking OLS would suffice but don't know if I am thinking too simply. I have re-fit my coefficients explicitly where the defaults (0,1) at t=k are regressed on the covariates at t=k for all businesses who have survived t>=k quarters, as opposed to defaults at t=k being regressed on covariates at t=0 in the original model (same number of observatsions, different time alignment). Thanks much for any help!


title: How to create a random variables in a simulation using skewness and kurtosis as well as average and standard deviation input?

text: I am curious to learn whether there are any best practises in creating random variables for a Monte Carlo simulation using input such as skewness and kurtosis information of a particular distribution. Of course those need to be used in addition to the average and standard deviation. My objective is to create a number series that better captures the shape of a distribution than just using the average and standard deviation. A next question would be: what functions in Excel and/or R to use generate the numbers on this basis?


title: How to implement dummy variable using n-1 variables?

text: If I have a variable with 4 levels, in theory I need to use 3 dummy variables. In practice, how is this actually carried out? Do I use 0-3, do I use 1-3 and leave the 4's blank? Any suggestions? NOTE: I'm going to be working in R. UPDATE: What would happen if I just use one column that uses 1-4 corresponding to A-D? Will that work or introduce problems?


title: Recurrence definition for a Markov chain

text: We define a state i to be recurrent if $\sum\limits_{n=0}^\infty P(X_n=i,X_k \neq i$ for $1\leq k < n | X_0=i)$=1. Why do we take infinite series over the probability? Why don't we define recurrent just as: $P(X_n=i,X_k \neq i$ for $1\leq k < n | X_0=i)$=1 Doesn't this already mean, that if we start in s, we will with probability 1 (almost surely) return to s, after n steps? Edit: Also, i is said to be recurrent if: $\sum\limits_{n=0}^\infty p_{ii}^{(n)}= \infty$. How does this fit logically with the other definition? Why would the series of the transition probabilities be unbounded?


title: Is it valid to analyze signal detection data without employing metrics derived from signal detection theory?

text: A signal detection experiment typically presents the observer (or diagnostic system) with either a signal or a non-signal, and the observer is asked to report whether they think the presented item is a signal or non-signal. Such experiments yield data that fill a 2x2 matrix: ![alt text](http://i.stack.imgur.com/JktFQ.png) Signal detection theory represents such data as representing a scenario where the "signal/non-signal" decision is based on a continuum of signal-ness on which the signal trials generally have a higher value than the non-signal trials, and the observer simply chooses a criterion value above which they will report "signal": ![alt text](http://i.stack.imgur.com/Vb2UM.png) In the diagram above, the green and red distributions represents the "signal" and "non-signal" distributions, respectively and the grey line represents a given observer's chosen criterion. To the right of the grey line, the area under the green curve represents the hits and the area under the red curve represents the false alarms; to the left of the grey line, the area under the green curve represens misses and the area under the red curve represents correct rejections. As may be imagined, according to this model, the proportion of responses that fall into each cell of the 2x2 table above is determined by:   1. The relative proportion of trials sampled from the green and red distributions (base rate)   2. The criterion chosen by the observer   3. The separation between the distributions   4. The variance of each distribution   5. Any departure from equality of variance between distributions (equality of variance is depicted above)   6. The shape of each distribution (both are Gaussian above) Often the influence of #5 and #6 can only be assessed by getting the observer to make decisions across a number of different criteria levels, so we'll ignore that for now. Additionally, #3 and #4 only make sense relative to one another (eg. how big is the separation relative to the variability of the distributions?), summarized by a measure of "discriminability" (also known as d'). Thus, signal detection theory proscribes estimation of two properties from signal detection data: criterion & discriminability. However, I have often noticed that research reports (particularly from the medical field) fail to apply the signal detection framework and instead attempt to analyze quantities such as "Positive predictive value", "Negative predictive value", "Sensitivity", and "Specificity", all of which represent different marginal values from the 2x2 table above (see here for elaboration). What utility do these marginal properties provide? My inclination is to disregard them completely because they confound the theoretically independent influences of criterion and discriminability, but possibly I simply lack the imagination to consider their benefits.


title: When is the shrinkage applied in Friedman's stochastic gradient boosting machine?

text: In boosting, each additional tree is fitted to the unexplained variation in the response that is currently un-modelled. If we are using squared-error loss, this amounts to fitting on the residuals from the aggregation of the trees fitted up to this point. I am not clear on whether it is at this point that the shrinkage (learning rate) is applied? E.g. we fit the first tree, and compute fitted values. Do we now shrink (down weight) these fitted values by the learning rate yielding a new set of fitted values, from which we compute the residuals that are used as the responses for the second tree in the ensemble, and iterate through many trees?


title: When is "Nearest Neighbor" meaningful, today?

text: In 1999, Beyer et al. asked, When is "Nearest Neighbor" meaningful? Are there better ways of analyzing and visualizing the effect of distance flatness on NN search since 1999? > Does [a given] data set provide meaningful answers to the 1-NN problem? The > 10-NN problem? The 100-NN problem? How would you experts approach this question today? * * * Edits Monday 24 Jan: How about "distance whiteout" as a shorter name for "distance flatness with increasing dimension" ? An easy way to look at "distance whiteout" is to run 2-NN, and plot distances to the nearest neighbor and second-nearest neighbors. The plot below shows dist1 and dist2 for a range of nclusters and dimensions, by Monte Carlo. This example shows pretty good distance contrast for the scaled absolute difference |dist2 \- dist1|. (The relative differences |dist2 / dist1| -> 1 as dimension -> ∞, so become useless.) Whether absolute errors or relative errors should be used in a given context depends of course on the "real" noise present: difficult. Suggestion: always run 2-NN; 2 neighbors are useful when they're close, and useful when not. ![enter image description here](http://i.stack.imgur.com/N23Vs.png)


title: Do likelihood ratios and Bayesian model comparison provide superior & sufficient alternatives to null-hypothesis testing?

text: In response to a growing body of statisticians and researchers that criticize the utility of null-hypothesis testing (NHT) for science as a cumulative endeavour, the American Psychological Association Task Force on Statistical Inference avoided an outright ban on NHT, but instead suggested that researchers report effect sizes in addition to p-values derived from NHT. However, effect sizes are not easily accumulated across studies. Meta-analytic approaches can accumulate distributions of effect sizes, but effect sizes are typically computed as a ratio between raw effect magnitude and unexplained "noise" in the data of a given experiment, meaning that the distribution of effect sizes is affected not only by the variability in the raw magnitude of the effect across studies, but also variability in the manifestation of noise across studies. In contrast, an alternative measure of effect strength, likelihood ratios, permit both intuitive interpretation on a study-by-study basis, and can be easily aggregated across studies for meta-analysis. Within each study, the likelihood represents the weight of evidence for a model containing a given effect relative to a model that does not contain the effect, and could typically be reported as, for example, "Computation of a likelihood ratio for the effect of X revealed 8 times more evidence for the effect than for its respective null". Furthermore, the likelihood ratio also permits intuitive representation of the strength of null findings insofar as likelihood ratios below 1 represent scenarios where the null is favoured and taking the reciprocal of this value represents the weight of evidence for the null over the effect. Notably, the likelihood ratio is represented mathematically as the ratio of unexplained variances of the two models, which differ only in the variance explained by the effect and thus is not a huge conceptual departure from an effect size. On the other hand, computation of a meta-analytic likelihood ratio, representing the weight of evidence for an effect across studies, is simply a matter of taking the product of likelihood ratios across studies. Thus, I argue that for science seeking to establish the degree of gross evidence in favour of a effect/model, likelihood ratios are the way to go. There are more nuanced cases where models are differentiable only in the specific size of an effect, in which case some sort of representation of the interval over which we believe the data are consistent with effect parameter values might be preferred. Indeed, the APA task force also recommends reporting confidence intervals, which can be used to this end, but I suspect that this is also an ill-considered approach. Confidence intervals are lamentably often misinterpreted (by students and researchers alike). I also fear that their ability for use in NHT (by assessment of inclusion of zero within the CI) will only serve to further delay the extinction of NHT as an inferential practice. Instead, when theories are differentiable only by the size of effects, I suggest that Bayesian approach would be more appropriate, where the prior distribution of each effect is defined by each model separately, and the resulting posterior distributions are compared. Does this approach, replacing p-values, effect sizes and confidence intervals with likelihood ratios and, if necessary, Bayesian model comparison, seem sufficient? Does it miss out on some necessary inferential feature that the here-maligned alternatives provide?


title: Generating random numbers from a t-distribution

text: How can I generate random numbers that follow a student- _t_ distribution? From several sources I understand that this can be done using a random sample of size $n$ drawn from a normally distributed population, as follows: $$ t = \frac {(x - m)} {(s/\sqrt n)} $$ Where $x$ is the sample mean, $m$ is the mean of the normal distribution (I assume you can just use the standard normal distribution, so $m=0$?), and $s$ is the sample standard deviation. The degrees of freedom of the student- _t_ distribution will then be $n-1$. Do I understand correctly that in order to generate a random student- _t_ value with $f$ degrees of freedom, I should first generate $f+1$ normally distributed values (i.e. standard normal), and then calculate the mean ($x$) and standard deviation ($s$) of these and apply the above formula? And if I do this repeatedly many times, the resulting random values will approach a student-t distribution with $f$ degrees of freedom? I tried this in Excel using a macro that uses the above formula and another macro that generates random Gaussians (which works, I tested it) but the resulting random values do not seem to be completely student- _t_ distributed. For instance with 6 degrees of freedom, the variance of 10,000 random values is about $1.7$ while it should be $6/(6-2) = 1.5$.


title: Identifying argmax result in information theory

text: I'm reading through a thesis at the moment which uses a "well known result from information theory", namely $\mathrm{argmax}_\mathbf{g}\,\mathbf{f}\,\mathrm{log}\,\mathbf{g} = \mathbf{f}$, for two probability vectors $\mathbf{f}$ and $\mathbf{g}$. I'm having trouble actually finding the source or a name for this result, though and my textbooks and Google seem to be drawing a blank. Any ideas? Thanks in advance!


title: How to calculate the shared variance between multiple IVs in predicting DV?

text: I would like to know how to calculate the degree to which multiple IVs share variances with each other in predicting for the DV. For example, A, B, C, D and E predicts for Z. How do I calcuate the degree of shared variance for A, B, C, D and E in predicting for Z? Thank you!


title: Establish minimum sample size for linear regression with historical data

text: I am working with a set of historical data that I did not collect myself. I cannot add to this data set in any way. For a 350-year period I have many thousands of data points, each associated with one of ~50 time points. Some of the time points have hundreds of data points associated with them, so I'm sure that those time points can be used in the linear regression. However, others have very few, even just one. So my question is, how can I determine a cut-off for the time points that have enough data (and should be included in my regression) and the time points that do not have enough data (and should be excluded)? I'd like to do this in R. Also, perhaps this is relevant - I'm interested in the minimum for each time point, not the mean. Thank you!


title: What is the meaning of McFaddens Axiom: Irrelevance of Alternative Set Effect?

text: On page 110 of _McFadden,1973 - Conditional logit analysis of Qualitative Choice Behavior, Frontiers in Economics, ed Zarembka, New York: Academic Press, pp. 105-142_ the following three Axioms are made: > **Axiom 1: Independence of Irrelevant Alternatives (IIA):** > > For all possible sets $B$, measured attributes $s$, and members $x$ and $y$ > of $B$: > > $P(x|s,{x,y})P(y|s,B) = P(y|s,{x,y})P(x|s,B)$ This is essentially the first part of the IIA - Axiome made by Luce, _Luce, R. D. 1959/2005. Individual Choice Behavior: A Theoretical Analysis. New York: Wiley. Reprinted by Dover Publications_ > **Axiom 2: Positivity:** > > $P(x|s,B)>0$ for all possbile alternative sets $B$, vectors of measured > attributes $s$, and $x\in B$. This is essentially the second part of the IIA - Axiome made by Luce. Now McFadden does a additional Axiom which says: > **Axiom 3: Irrelevance of Alternative Set Effekt:** > > The function $V(s,x,z)$ determining the selection probabilities in > > $P(x|s,B) =\ \frac{e^{V(s,x,z)}}{\sum_{y\in B}e^{V(s,y,z)}}$ > > has the additively separable form: $V(s,x,z) = v(s,x) - v(s,z)$ Further, he states that Axiom 3 follows from Axiom 1 and 2 if there exsists some "universial benchmark" alternative $z$ such that if $B$ is a possible alternative set, then $B\cup\\{z\\}$ is also. This follows by noting that Equation (9) holds for $z\not\in B$, provided Axiom 1 and 2 holds for $B\cup\\{z\\}$. Then, taking $z$ to be the universial benchmark in Equation (10) and defining $v(x,s)=V(s,x,z)$ for all alternative sets yields the result. Equation (9) $\frac{ p_{yx} }{ p_{xy} } = \frac{ p_{yz}/p_{zy} }{ p_{xz}/p_{zx} }$ follows from Axiom 1 and two where $p_{yx} = P(y|s,{x,y})$ and equation (10) $P(x|s,B) =\ \frac{e^{V(s,x,z)}}{\sum_{y\in B}e^{V(s,y,z)}}$ follows from Axiom 1 and 2 if we set $V(s,x,z) = log(\frac{p_{xz}}{p_{zx}})$ **My question in this context is:** Why do we need this additional assumption, i.e., Axiom 3? I intuitively understand, that alternative set effects cannot be measured if no replications on the same individual are measured. What I do not understand is, how this _"universal benchmark"_ or _"benchmark"_ idea is related to this problem? Will every choice from $B$ be related to some item $z\in B$, which will then be called benchmark? What does the " _irrelevance"_ in this context mean?


title: A function like lmer?

text: Can you suggest me any random effect linear model in R? At the moment I am using **lmer** , but I would like to know if there is any other similar function. Thank you


title: Shapiro test on replicated values

text: I have a distribution of frequencies (1000 data points) expressed as %. The list of data looks like this: > 3.10%   >  1.80%   >  1.70%   >  1.70%   >  1.60%   >  1.60%   >  1.50% I would like to check for normality distribution of such data using for example the Shapiro test. I would like to know if I have to remove replicated or if they will stay there, because I performed the test and with the replicated numbers the output is: > Shapiro-Wilk normality test   >  W = 0.9422, p-value < 2.2e-16 and after removing replicates the output is: > Shapiro-Wilk normality test   >  W = 0.8911, p-value = 0.05804 So the output change. Is it correct to maintain replicated values since they correspond to real measurements?


title: Scikit-learn reports memory error when fitting Gaussian process model

text: I want to fit a Gaussian Process with about 50,000 training examples and 130 features using Scikit-learn. Right now, I only have 1 theta hyperparameters as I run the process with all defaults. But I get a MemoryError when it tries to compute the covariance matrix. This is where I get the error https://github.com/scikit-learn/scikit- learn/blob/master/sklearn/gaussian_process/gaussian_process.py#L52 What can I do?


title: Requirements for evaluation function in evolutionary algorithms

text: For a publication I'm looking for references for generic requirements of evaluation/fitness functions in evolutionary (more specifically genetic) algorithms. I could come up with some requirements myself (efficient implementation, intuitive results, etc.) but I would like some references to existing work, since it should not be new. Any reference, preferably to papers but if necessary to books, is appreciated.


title: Visualization and navigation through a Pareto front

text: I've recently extended my genetic algorithm in such a way that it constructs a Pareto front, instead of weighing the different quality dimensions. The quality dimensions are not related to 'spacial' information, so no natural x,y,z candidates. They are real quality dimensions (for instance recall, precision, cost, etc.). Now I'm faced with the challenge of providing visualization and navigation options to the user to investigate the Pareto front of candidates. Currently I'm using simple navigators per dimension and I've implemented a simple dot plot like below: Example of used chart I've seen research that uses a 3D version of the dot plot as visualization aid. However, I'm keen to know of other visualization solutions or best practices regarding Pareto fronts. Any help is appreciated.


title: Looking for sparse and high-dimensional clustering implementation

text: I'm looking for a clustering implementation with the following features:   * Support for high-dimensional data. Now I have approximately 160.000 dimensions/features.   * Be able to manage sparse matrix. That is, not only to read sparse matrices, but also capable of making operations in this format.   * Properly shows the centroid for each cluster. I've tested some packages:   * Rapidminer, which seems to be a memory eater, I suppose because although capable of reading a sparse matrix, it is not capable of working with them as they are.    * Cluto, which is very fast and low-memory consumption, but it is not able of show properly the centroid elements (source code not available). It shows descriptive features together with a percentage of how that feature contributes to the average similarity, but there is no clear info (here is a question about that, with no clear answer) about how is calculated that, and also I have clusters where there is 0.0% but it is not clear for me if this means the program is unable to show an upper precision or if that feature has nothing to do tho the average similarity. I appreciate any comment or answer about it.


title: Generate random correlation matrix that has approximately normal entries

text: I would like to generate a random correlation matrix such that the distribution of the entries (except 1) look approximately like normal. The motivation is this. For a set of n time series data, the correlation distribution often look quite close to normal. I would like to generate many "normal" correlation matrices to represent the general siutation and use them to calculate risk number. Do you think I answer your question?


title: Interpretation of odds ratio when outcome is a percentage

text: I ran a GEE model, with a dependent variable of "percent of total students with an unexcused absence," using a binomial family. My dependent variable is basically a proportion, with range from 0 to 1 to obey the rules of the binomial family. I understand that the difference between the non-event/event 0-1 is "zero percent of students with an unexcused absence" and "100% of students with an unexcused absence." Questions:   * How can I interpret the Odds Ratio for this in a way that makes policy-level, real-world sense? An interpretation of 1-unit change in independent variable X would be e.g. "a school with a 1-unit increase in X has 0.9 times the odds of having 100% of students unexcused from class at least once, compared to a school without the increase in X." Seems clunky! You may ask why I used a binomial family rather than a Gaussian family, but according to "Comparison of Logistic Regression and Linear Regression in Modeling Percentage Data" by Zhao, Chen and Schaffner (2001), logistic regression can and should be used for any models with a dependent variable that's modeled as a percentage for various reasons. If there's no other accurate way to interpret this coefficient through logistic regression and ORs, does anyone have any suggestions about how to model the data to be better interpretable?


title: Is there any way I can manually set the intercept in SPSS?

text: As part of my dissertation I am looking into the relationship between temperature and methane emissions. A large proportion of this work relies on statistical interpretation and analysis. When plotting the link between flux and temperature SPSS outputted a 'quadratic' line of best fit with very little variation from the original linear plot. This doesn't fit the expected trend, however when I select the option to suppress the intercept the line of best fit follows the expected trend far better. A previous thread here pointed out some issues with this technique though. If I chose to suppress the intercept it makes the assumption at temperature 0, flux is also 0 which is not the case. I was wondering if there was a way to 'suppress' the intercept as it were in order to maintain the improved line of best fit, but manually set the intercept to -5 for example. Any pointers would be greatly appreciated.


title: How should one interpret the beta weights for overlapping regressors in a finite impulse response model?

text: I am trying to fit a general linear model in which each regressor consists of 1s and 0s only. Suppose my design matrix is as follows:               1 0 0 0 0 0 0 0 0      0 1 0 0 0 0 0 0 0     0 0 1 0 0 0 0 0 0     0 0 0 0 0 0 1 0 0     0 0 0 0 0 0 0 1 0     0 0 0 0 0 0 0 0 1     0 0 0 1 0 0 0 0 0     0 0 0 0 1 0 0 0 0     0 0 0 0 0 1 0 0 0     1 0 0 0 0 0 0 0 0     0 1 0 0 0 0 0 0 0     0 0 1 0 0 0 0 0 0     0 0 0 1 0 0 0 0 0     0 0 0 0 1 0 0 0 0     0 0 0 0 0 1 0 0 0     0 0 0 0 0 0 1 0 0     0 0 0 0 0 0 0 1 0     0 0 0 0 0 0 0 0 1      This matrix is 18 X 9: It includes 9 regressors for 18 data points. The regressors can be divided into three groups: The second regressor is the same as the first, but shifted down one row. The third is the same as the second, but shifted down one more row. Within the middle three and last three regressors the same relationship is seen. I know how to interpret the beta weight for each regressor in this design matrix. The weight is simply the mean of the data points modeled by the regressor. This is true because each row in the design matrix contains a single 1; there is no overlap between any of the regressors. In the case of the first regressor (the leftmost column), the beta weight is the mean of data points 1 and 10. Suppose I were to replace this design matrix with the one below, for a new dataset in which there are only 13 data points:               1 0 0 0 0 0 0 0 0      0 1 0 0 0 0 0 0 0     0 0 1 0 0 0 1 0 0     0 0 0 0 0 0 0 1 0     0 0 0 1 0 0 0 0 1     0 0 0 0 1 0 0 0 0     1 0 0 0 0 1 0 0 0     0 1 0 0 0 0 0 0 0     0 0 1 1 0 0 0 0 0     0 0 0 0 1 0 0 0 0     0 0 0 0 0 1 1 0 0     0 0 0 0 0 0 0 1 0     0 0 0 0 0 0 0 0 1      Notice that the only difference between this design matrix and the one above it is that the rows have been squashed: Every third row in the original design matrix has been merged with the row after it, such that some rows in the new design matrix have two 1s. My question is, how do I interpret the beta weights for regressors that overlap other regressors in exactly this way? What I mean is, if I wanted to compute the beta for each regressor by hand, how would I do it? In the case of this particular design matrix, is there only one way to assign beta weights? Do other kinds of matrices with overlapping regressors not produce unique weights? How can I tell whether a given matrix will produce unique weights?


title: How to properly test for difference in number of user actions?

text: Suppose you divide users into two groups by flipping a coin for each user. Now you measure the number of actions in each group (say, number of total clicks). How do you test if the two groups are different? For the number of users in each group, I would use a chi-squared test between expected # users (total/2) and observed in each group. What test can I use for the number of actions? The counts of actions are correlated, because the same user contributes multiple actions (clicks in our example). Moreover, it may be that some users contribute far more than others. I could use a permutation test or a bootstrap, but that does not scale well with the number of users. (A chi-squared test that computes group sums on the database or Hive side is much easier than a bootstrap where maybe millions of user counts have to be transferred into R.) Thoughts? Maybe I can use a chi-squared test anyway because the sum of # of actions should converge to normal (central limit theorem)?


title: Do edges in directed acyclic graph represent causality?

text: I am studying _Probabilistic Graphical Models_ , a book for self-study. Do edges in a directed acyclic graph (DAG) represent causal relations? What if I want to construct a Bayesian network, but I am not sure about the direction of arrows in it? All the data will tell me is the correlations observed, not the inter-linking between them. I know I am asking too much, as I am sure following chapters will address these issues, but it's just that I can't stop thinking about it.


title: Given probability x, how likely is it that probability y is due to chance?

text: I wrote a Python script searching text for words and separating unique and non unique words. Running the pages I have through my program, I find 1041 unique words out of 3742 total, so roughly 27% of all the words are unique. Yet when I look at the first word of each page, 33 of the 44 words starting each page are unique, for a 75% ratio. How do I measure how likely this 75% result is due to chance, given the 'normal' 27% unique word probability? I have some experience in R if that helps at all in explaining things to me.


title: Negative Hessian matrix in R optim

text: I used R optim() to find min log likelihood, however the diagonal elements of the inverse of Hessian matrix turns to be negative..               > round(diag(solve(KFopt$hessian)),2)     [1]    0.00   -0.08   -0.03    0.00   -0.47   -0.18 -167.32      Does that mean the optimization is wrong? So what should I do if I want to extract standard deviation of each estimates? Thanks in advance!


title: Appropriate normality tests for small samples

text: Title says enough... so far, I've been using Shapiro-Wilk's statistic in order to test normality assumptions in small samples. Could you, please, recommend another technique?


title: Intraclass correlation and aggregation

text: Imagine that:   * You have a sample of 1000 teams each with 10 members.   * You measured team functioning by asking each team member how well they think their team is functioning using a reliable multi-item numeric scale.   * You want to describe the extent to which the measure of team effectiveness is a property of the team member's idiosyncratic belief or a property of a shared belief about the team. In this and related situations (e.g., aggregating to organisations), many researchers report the intraclass correlation (e.g., Table 1 in Campion & Medsker, 1993). Thus, my questions are:   1. What descriptive labels would you attach to different values of the intra-class correlation? I.e., the aim is to actually relate the values of the intra-class correlation to qualitative language such as: "When the intraclass correlation is greater than x, it suggests that the attitudes are modestly/moderately/strongly shared across team members."   2. Do you think the intraclass correlation is the appropriate statistic or would you use a different strategy?


title: Using mixed effects modelling to estimate and compare variability

text: Say I observe two groups of 10 people, measuring some quantity 100 times in each person. There will presumably be some variability across these 100 measures in each person. Can I use mixed effects analysis to assess whether this within-person variability is, on average, different between the two groups? For example, using traditional statistics, I could compute the standard deviation (SD) within each person then submit these SDs to an anova comparing the groups, but I wonder if this two-stage process can be replaced by a single mixed effects model, consequently obtaining the various advantages of mixed effects modelling (shrinkage, accounting for different numbers of observations per person, etc) as well. To be clear, here is R code depicting the scenario described above and the SD/anova-based approach:               set.seed(1)          group_A_base_sd = 1     group_B_base_sd = 2     within_group_sd_of_sds = .1          n_per_group = 10     obs_per_id = 100          temp = data.frame(         id = 1:(n_per_group*2)         , group = rep(c('A','B'))     )          #generate example data     library(plyr) #to avoid loops (for coding convenience only)     obs_data = ddply(         .data = temp         , .variables = .(id,group)         , .fun = function(x){              #generate a unique sd for this individual              # based on their group's sd plus some              # within-group variability              id_sd = ifelse(                   x$group=='A'                   , rnorm(                        1                        , group_A_base_sd                        , within_group_sd_of_sds                   )                   , rnorm(                        1                        , group_B_base_sd                        , within_group_sd_of_sds                   )              )              #generate data points with the above generated              # variability              to_return = data.frame(                   obs_num = 1:obs_per_id                   , measurement = rnorm(obs_per_id,0,id_sd)              )              return(to_return)         }     )          #first step of an anova-based approach:     # compute SDs within each Ss     obs_sds = ddply(         .data = obs_data         , .variables = .(id,group)         , .fun = function(x){              to_return = data.frame(                   obs_sd = sd(x$measurement)              )         }     )          #second step of an anova-based approach:     # compute the anova on the SDs     summary(         aov(              formula = obs_sd~group              , data = obs_sds         )     )


title: Normalization for dynamic Bayesian network inference?

text: I would like to know whether preprocessing of the dataset is required for dynamic Bayesian network inference or not?


title: Linear regression across $n$ samples

text: Consider a random variable $Y$ with mean: $y = \alpha + \beta x$ and variance $\sigma^{2}$. Assume $n$ independent values of $Y$ are observed and are split into two samples such that the samples can be labelled: $n_{1}$ and $n_{2}$ ($n_{2} = n-n_{1}$); furthermore, assume that the observations in both groups share a slope/gradient ($\beta$), but have differing intercepts ($\alpha$). How can it be shown that: $$\hat{\beta} = \frac{w_{1}\hat{\beta_{1}} + w_{2}\hat{\beta_{2}}}{w_{1}+w_{2}}$$ where $\hat{\beta_{1}}$ and $\hat{\beta_{2}}$ are the estimates for $\beta$ derive from each sample and $w_{1}$ and $w_{2}$ are the sums of squares of the $x$ values (independent variables) about their mean for each sample. In addition, what form does $Var(\hat{\beta})$ have? Finally, for clarity, following on from the above model, is it correct to assume that, the intercept ($\alpha$) for each sample will be given by: $$\alpha_{j} = \bar{Y}^{(j)}-\beta_{j}\bar{X}^{(j)} \text{, for } j=1,2$$ And the gradient/slope ($\beta$) is given by: $$\hat{\beta_{j}} = \frac {\sum_{i=1}^{n} X_{i}^{(j)}Y_{i}^{(j)}-n\bar{X}^{(j)}\bar{Y}^{(j)}} {\sum_{i=1}^{n} X_{i}^{(j)2}-n\bar{X}^{(j)2}} = \frac{S_{XY}}{S_{XX}} \text{, for } j=1,2$$ (Note: The $(j)$ notation in the above ($X^{(j)}$, $Y^{(j)}$, etc.) is intended to denote sample to which each data point belongs -- It is _not_ exponential notation.) Edit: I should probably mention that, originally, I had assumed that I could use the method of least squares to sum the $SS_{E}$ terms from each sample, set the partial derivative w.r.t. $\beta$ equal to zero and solve for $\beta$; such that: $$ \sum_{j=1}^{2} \frac{\partial SS_{E}^{(j)}}{\partial \beta} = \frac{\partial SS_{E}^{(1)}}{\partial \beta} + \frac{\partial SS_{E}^{(2)}}{\partial \beta} = 0$$ $$ = \frac{\partial} {\partial \beta} [\sum_{i=1}^{n - n_{1}} (y_{i}^{(1)}-\hat{\alpha_{1}}-\hat{\beta}x_{i}^{(1)})^{2} + \sum _{i=n+1-n_{1}}^{n} (y_{i}^{(2)}-\hat{\alpha_{2}}-\hat{\beta}x_{i}^{(2)})^{2}] = 0$$ $$ = -2 \sum_{i=1}^{n - n_{1}} [(y_{i}^{(2)}-\hat{\alpha_{1}}-\hat{\beta}x_{i}^{(2)})(x_{i}^{(2)})] -2 \sum_{i=n+1-n_{1}}^{n} [(y_{i}^{(2)}-\hat{\alpha_{2}}-\hat{\beta}x_{i}^{(2)})(x_{i}^{(2)})] = 0$$ However, this does not appear to be equivalent to the form given above, so I'm assuming this is not a valid assumption.


title: Better machine learning methods for medical imaging segmentation

text: I am wondering if there are some better machine learning methods for medical imaging segmentation? Currently, I have some relatively low-resolution MRI images and I tried to use histogram and k-means algorithm, which is an old method to do a general clustering based segmentation. However, pixel may help estimate more than one cluster in the k-means procedure. Is there any novel direction for clustering pixel values in order to do the segmentation? I think this can be done by two ways:   1. Change the clustering similarity measurements but still use k-means, e.g., change histogram based clustering into something like using mutual information, but it might be well done by other people (if you know any, please give me the reference).   2. Change k-means to other clustering methods but keep the similarity measurements, i.e., histogram. Any ideas would be appreciated. Thanks a lot. A.


title: Multiple Linear Regression with Combined Variables

text: Question: > Suppose we have an ordinary household scale such as might be used in a > kitchen. When an object is placed on the scale, the reading is a combination > of the true weight plus a random error. You have two coins of unknown values > B1 and B2. To estimate the weights of the coins, you take four observations: > >   * Put coin 1 on the scale and observe y1   > >   * Put coin 2 on the scale and observe y2   > >   * Put both coins on the scale and y3   > >   * Put both coins on the scale and observe y4   > > > > (a) Write a linear model in matrix form and find the least-squares estimate > of the coins using the usual formula (X'X)^-1 X'y.   >   1. Regarding the question from letter (a), my understanding is that for either B1 or B2, they can be estimated by averaging the weights (for example, for B1 it would be adding the weights from y1, y3 and y4). Is this assumption valid?      I thought it was, and so I tried to set up the formula for B1:              B1 = 1/3*K1*Y1 + 0*K2*Y2 + 1/6*K3*Y3 + 1/6*K4*Y4 + e1 + e2 + e3 + e4, and       B2 = 0*K1*Y1 + 1/3*K2*Y2 + 1/6*K3*Y3 + 1/6*K4*Y4 _ e1 + e2 + e3 + e4,       where `Ki` is the slope parameter estimate for the `ith` coin and `ei` is the residual for the `ith` coin. The numbers 1/3 and 1/6 (1/3 * 1/2) to find the average weight of `B1`. Am I totally off base with this reasoning and subsequent setup of my linear model?   2. Can I also assume that `y3` is found by simply adding `y1` and `y2`? I didn't think I could, because the readings would probably be different and thus have different residuals.   3. Does the `y` in the formula correspond with the estimated column vector containing predicted `B1` and `B2`?


title: Classifier predicts only one class

text: I was trying myself in kaggle CIFAR competition, I trained lots of classifiers but get the same result/fail (don't know how to treat them), maybe someone could help me figure what i'm doing wrong.   1. So in the beginning I converted images to vector (included all colours) and tried to train classifier on a sample of training set. There was no preprocessing and I used random forest method. It scored about 70% percent on a training set, however when I try to apply my classifier on test set it classifies all points as one class. (in sample size was about 10k points in 3k dimensions with 500,1000 trees)   2. Then I figured that I have to do dimensionality reduction, I figured out the amount of features I'm needed (`nFactors` library in R) and tried to train random forest classifier again on ICA projections. Prediction on test set again spits out only one class.   3. Next thing I figured that maybe my method is not calibrated enough, so I sampled small amount of data and tuned svm classfier on it. Then I trained it on test sample, which resulted in 80% error on training set, but when I predict class on test set I still get one class prediction. Then I also tried create LogitBoost classifier, it classifies 50% of data quite correctly (leaving rest 50% NA's) but when i try it on test sample it results in all NA's.  What am I doing wrong? It looks like my classifier can't catch any underlying structure, though I tuned parameters, increased amount of trees to like 5k and did feature extraction, but still I get only one class prediction. (usually the one that have biggest error) Since other teams are actually doing quite well it means that something gone really wrong.


title: Computing conditional expectations

text: > **Possible Duplicate:**   >  Conditional Expectation / Estimator Confusion Let $X_1, X_2, X_3 \sim N(0, d^2)$ and $T = X_1^2 + X_2^2 + X_3^2.$ I have an estimator for $d$, $$\hat{d} = \frac{\sqrt{T\ 2\pi}}{4},$$ and another estimator for $d$, $$\tilde{d} = \frac{1}{3} \sqrt{\pi / 2}\ \left(|X_1| + |X_2| + |X_3|\right).$$ I need to show that $\mathbb{E}(\tilde{d} | T) = \hat{d}$ and that $\mathbb{E}(|X_1| | T) = \frac{1}{2} \sqrt{T}$. I'm mostly confused about how to proceed with the algebra to compute the expectations. I have never had to compute E[estimator | something] before


title: Does a positive Mantel-Haenszel chi square statistic in a trend test indicate a positive trend?

text: When running my analysis on a dataset, I find a positive Mantel-Haenszel chi square trend statistic. Does this mean that the trend in the data is positive or simply that there is a trend? The data deals with the improvement of patients' symptoms after taking a drug or a placebo.


title: How to select best parameter for polynomial kernel?

text: I am using LibSVM library for classification. For my problem I am using polynomial kernel and I need to select best parameters (`d` = degree of polynomial kernel, and `C` = soft margin constant). The LibSVM guide suggests for grid search for this job. In this library there exists a parameter selection tool (`grid.py`), but it seems to me this is implemented to tune parameters of RBF kernel (gamma). My questions are:   1. Is there any other good solution for selecting suitable parameters rather than Grid search?   2. Can anybody give me some hints / sample code for selecting parameters for polynomial kernel?


title: Calculating the probability of the Nth sample in the birthday problem w/o knowing the previous samples

text: If we draw two random variates from a discrete uniform distribution $[1, D]$, the probability that the samples are distinct is $(D-1)/D$. Explanations of the birthday problem state that if we sample a third time, the probability that this third sample is distinct from the first two is $(D-1)(D-2)/D^2$. Doesn't this assume that the first two samples were distinct? If they are not distinct, then the probability that the third sample is distinct is $(D-1)/D$. The solution to the birthday problem requires us to know the probability that the $n^\textrm{th}$ random variate is distinct from the previous samples. How can we calculate the probability that the $n^\textrm{th}$ sample is distinct from the previous samples without knowing the values of these previous samples?


title: Which weighting factor to use for text categorization

text: I am working on a text categorization task, and I possess 21,000 documents for training, and (for the time being), 7000 documents for testing. I construct the doc-term matrix for both training corpus and testing corpus, with two different weighting factors, i.e. **_TF_** (term frequency) or **_TF-IDF_** (term frequency–inverse document frequency). Then I used SVM with Gaussian Radial kernel, to classify the documents. The F1-measure for tf-idf weighting is nearly 0.8, while for tf weighting, the performance is less better, around 0.7. So, logically, We'll be prone to use tf-idf weighting. However, one problem will come out in an **_INCREMENTAL CONTEXT_**. That is, when we have to categorize a single or a few new documents from time to time (with a pre-trained model). It will not be suitable to use tf-idf weighting for one single document, since tf-idf is often used to measure a word importance in a collection of documents. Should I compromise to use tf weighting or some other tricks exist? I will be very grateful for any hints or advices to solve this problem! Thanks a lot !


title: Creating variogram for a 10,000 set data

text: I am trying to create a variogram for a data set with 10,000 points. However, if I try to actually calculate the distance of each point with the other then I will have 10,000*9999/2 pairs. Out of these pairs of distances I can round off the distance values to lets say 2 decimal places. Then I can find the distances having equal values. Then I can take the average of the variances having equal distances to find the variance for that particular distance. Then I can get the experimental variogram. This process will definitely be very slow. Is there any efficient way ? I mean instead of creating the variogram from all the observation point, I can only take a subset of it surrounding the point where I want to interpolate the value. I can create a variogram out of this subset. Then I can further take k neighbors of the destination point and use this subset variogram to interpolate. Will this be more efficient and correct?


title: Multiple comparisons

text: I would like to compare a handful of percentage values - whether they are significantly different from each other. Since the t-test only compares means, it doesn't seem appropriate. Multiple comparisons such as Scheffe's gave me a warning in SAS that only CLASS variables were accepted. What kind of test should I be looking for?


title: Detecting and preventing collusion in ratings

text: I am not sure whether or not this problem can be solved through statistical methods, but here is the problem: I have 4 professional reference providers for a job candidate rating the candidate on a scale of 0-100 on say 10 items. Behind the scenes we check for the consistency of ratings for each item from the 4 reference providers and arrive at a total score for each of the 10 items. The problem is what is the 4 reference proviers collude with the candidate and give a rating of 100/100 on every item? Such scores cannot be prevented because there may really be an excellent candidate who can score 100/100 on every item from each of the 4 reference providers. So, is it possible to design a statistical method to differentiate between the fake candidate colluding with reference providers and the excellent candidate getting all 100 scores?


title: Streaming k means clustering Mahout

text: First of all, excuse me if this is not a good place to ask this question. Can anyone explain to me how streaming kmeans algorithm implemented in Mahout works? And can it be used for anomaly detection?


title: Regression problem with "included classification"

text: I have basic knowledge on machine learning techniques and I'm working on a real-world regression problem which is basically predicting the time consumption of a certain task based on a number of predictor values. Unfortunately (as usual) the underlying relationship is not known to me. However, I do assume that there are different types of tasks in my data and therefore I'm searching for a technique to handle this. Breaking it down to a simpler regression problem: I have one predictor p1, a categorial predictor p2, and the resulting time t, and my data looks looked like this: ![enter image description here](http://i.stack.imgur.com/bCT50.jpg) The regression becomes rather easy, as long as p2 is used by the regression model. My Questions are:   1. Is there any type of machine learning technique for regression that does this kind of "included classification"? What I means is that it automatically builds regression models for the seperate classes, if that improves the model.   2. Would that technique work if the second predictor would not be categorial but "Type A" had p2 between 10 and 20, and "Type B" had p2 values between 30 and 40?   3. What are good approaches to handling categorial predictors in regression?   4. Any other advice for my described real-world problem? Thanks for your help //Edit 17.12.: Removed Example without p2


title: In a census-like situation, does it make sense to use Bayes Rule?

text: I'm involved in a census at our condo association. We may, for logistics reasons, survey 1/4 of the population per quarter so that we have surveyed everyone by the end of the year, but let's assume we survey everyone at once. A year later, we survey again. Since we're attempting to cover 100% of the population each year, the results are what the results are each year, right? It wouldn't make sense to run last year's results (prior) and this year's results through Bayes Rule would it? EDIT: In response to the excellent questions, I'm thinking of a single question on the survey that is asked both years. For example: what percentage of residents believe that the Management Office does a good job of handling complaints? At first glance, this seems textbook-like: you have a prior distribution based on the previous poll, you gather data, use Bayes, and now you have a posterior which in principle better reflects reality than the prior or current data. EDIT2: For this question, I'm assuming that there is no other data: we simply have (anonymous) polling results (and uncertainty) from one question asked in two different polls and wonder if the data for that question can be combined in some manner. EDIT: As the answers point out, the problem is that the two samples are not independent, regardless of how far apart they are, because many residents will be asked both years. (Which gets me thinking: how independent is "independent"? Do we say that we need a method of choosing respondents that has less than 95% (99%, etc) probability of picking the same respondent twice?) If we were taking samples from a much-larger population, I can see how Bayes Rule would be helping us to solidify our answers, but in this case... Actually, I'm not sure how using Bayes Rule on data that is evolving over time (polls, say) would work out, since the prior is "out of date" in some sense. Questions:   1. Does it make any sense in a census (attempting to survey the entire population) situation to run the prior census' results and the current census' results through Bayes Rule?   2. If we only get a 50%, say, response rate would that change the answer to #1?   3. Does it make sense to use, say, polling data as a prior for polling if we know that the data (opinions) are changing over time? (i.e. aren't priors a statement of prior belief about the current situation, not just arbitrary beliefs at some time in the past?) EDIT: It appears that once the population's answers are subject to change, we then have to actually model this change (at least accounting for time), at which point we've moved beyond a simple application of Bayes Rule. IS that correct?


title: How would you filter this?

text: This image show a histogram (200 bins) of accumulated distances from a radar distance meter (very noisy). The peak around 7 meters is an object. At thought this looked kind of like a normal distribution, at least if you ignore values <4m (which for this application is reasonable). I have also achieved a reasonable fit with a log-normal distribution. What I am trying to do is to filter out true distances based on the probability distribution. ![alt text](http://i.stack.imgur.com/tYgBf.png) **EDIT:** log-normal distribution: ![alt text](http://i.stack.imgur.com/9UMoY.png) **Note:** this question can also be found at math.stackexchange.com


title: Statistical significance with sample size of 5 or less

text: Should small sample sizes of 5 and less be included in reports? This is an area of medical training that is currently being reported only to the course level trainers and statistically not included in reports to upper management. All of the data is mined from course-level surveys.


title: No valid coefficients for NegBin regression

text: I am doing multiple regression with some data (5 predictors, 1 response). Since the response is discrete and non-negative, I thought I would try Poisson regression. However, the data are significantly overdispersed (variance > mean), so I am now trying negative binomial regression. I was able to fit the model with this code.               library(MASS)     model.nb <- glm.nb(Response ~ Pred1 + Pred2 + Pred3 + Pred4 + Pred5 - 1, data=d)      Now I would like to see if I can get a better fit by including interactions between the predictors. However, when I try to do so, I get the following error.               > model.nb.intr <- glm.nb(Response ~ Pred1 * Pred2 * Pred3 * Pred4 * Pred5 - 1, data=d)     Error: no valid set of coefficients has been found: please supply starting values      Any ideas what may be causing this?


title: Help on SVM for road image processing

text: I am new to SVM. I would like to use SVM to segment/cluster/classify a road image into two distinct regions i.e. drivable region and non-drivable region. Unfortunately, I do not have any images where the regions are already labeled. I really don't know how to go about this i.e. specifying the training set. Will weka be a good tool to use for the SVM? Thanks in advance


title: Better method to estimate gender gap: difference between medians or median gap by occupation?

text: I'm doing some analysis of Bureau of Labor Statistics data on wages earned by men and women, and I'm trying to determine which of the two methods I've used is the better estimate of the gender wage gap. The first compares the median weekly earnings of all women employed full-time (669 USD) against the median weekly earnings of all men employed full-time (824 USD) and determines that the wage gap is about 81:100. The second takes all of the occupations for which a wage gap can be computed and weights them according to the number of women in each occupation. The median wage gap, based on this method, comes out to about 86.5:100. Are both methods fairly sound? If so, is one better than the other?


title: What kind of mean to use for a mix of negative and positive percent values?

text: I've got a question from a colleague about what sort of mean to use in calculating the average of the deviations of actuals from a forecast, which are in percent. The data look(s) like this:               -0.26     -0.21     -0.9     1.56     1.88      For example, the last unit above beat the forecast by 1.88%. My intuition was geometric (based on this question), but I am not sure how to to deal with the negative values. Edit: I think I've figured it out once I had some more coffee. I think you want to transform the data to (1+percentage). For example, enter +4% as 1.04. Similarly, a percentage of -2% would be .98. Would that be correct?


title: Weighting variables for an index

text: I have been tasked with trying to modify our current "index" which basically takes 4 observations per person and calculates a score based on what they achieve. Here is how the score is created (all variables are ordinal in nature, higher scores representing more desirable outcomes):   * Variable A - gets a score of 1 to 3    * Variable B - gets a score of 1 to 3   * Variable C - gets a score of 1 to 5   * Variable D - gets a score of 1 to 8 Basically what I am looking to do is revise our index but have little experience doing so using a statistical framework. What techniques exist? My "problem" is that I will most likely be required to keep the same max score (19, maybe 20) and use the same 4 variables. Lastly, I will need to be able to explain the score in a way that our MIS department can code the index back into our database. Any help or insight you can provide to get me started will be very much appreciated.