title: Test if probabilities are statistically different?

text: I completed a Monte Carlo simulation that consisted of one million ($10^6$) individual simulations. The simulation returns a variable, $p$, that can be either 1 or 0. I then weight the simulations based on predefined criteria and calculate the probability of $p$. I also calculate a risk ratio using $p$: $$\text{Risk ratio} = P(p|\text{test case}) / P(p|\text{control case})$$ I had eight Monte Carlo runs, which consist of one control case and seven test cases. I need to know if the probabilities of $p$ are statistically different compared to the other cases. I know I can use a multiple comparison test or nonparametric ANOVA to test individual variables, but how do I do this for probabilities? * * * For example are these two probabilities statistically different?: _Probabilities_ : $P(p|\text{test #3}) = 4.08 \times 10^{-5}$ $P(p|\text{test #4}) = 6.10 \times 10^{-5}$ _Risk Ratios_ : $\text{Risk Ratio}(\text{test #3}) = 0.089$ $\text{Risk Ratio}(\text{test #4}) = 0.119$


title: Regressing, analysing data with points rather than polynomial?

text: I am looking into making a regression of a bunch of data that is contained on some range of real numbers. In my case, x is between 0 and 1 and y is between 0 and 10. If I have 150 data points on this plane, and want to best model the original data with N data points (not polynomial, just points), how can I do this? How can I then analyse how good the fit is? Any help is appreciated. Sorry if my notation is not perfect, I am from a physics background so please bear with me. Cheers. Data example:               x.vector <- cumsum(rnorm(150)^2) y.vector <- 10*rnorm(150)^2      I want to fit 2 points, 3, 4, ... and check fit at each one.


title: Determine if three is statistically different than ten for a very large number of observations (1,000,000)

text: For 1,000,000 observations, I observed a discrete event, X, 3 times for the control group and 10 times for the test group. How do I determine for a large number of observations (1,000,000), if three is statistically different than ten?


title: How are these negative values understood?

text: I can't find literature to understand this table. What I know doing my research is that an effect size should be between 0 and 1. when 0.2 is slow, 0,5 medium and 0.8 and higher , high. But In an article of Cochrane database I found these results and I can't really make a sense of this. I need this really badly so anyone who can help me to understand these values (including those ranges in brackets) I'd really appreciate it. ![Effect size](http://i.stack.imgur.com/k4uR3.png)


title: lm() producing many NAs for coefficients

text: I am trying to run a regression using about 80 independent variables. The problem is that the last 20+ coefficients return NA. If I condense the range of data to within 60, I get coefficients for everything just fine. Why am I getting these NAs and how do I resolve it? I need to reproduce this code using all of these variables.               composite <- read.csv("file.csv", header = T, stringsAsFactors = FALSE)     composite <- subset(composite, select = -ncol(composite))     composite <- subset(composite, select = -Date)     model1 <- lm(indepvariable ~., data = composite, na.action = na.exclude)      composite is a data frame with 82 variables. UPDATE: What I have done is found a way to create an object that contains only the significantly correlated variables, to narrow the number of independent variables down. I have a variable now: sigvars, which is the names of an object that sorted a correlation matrix and picked out only the variables with correlation coefficients >0.5 and <-0.5. Here is the code:               sortedcor <- sort(cor(composite)[,1])     regvar = NULL          k = 1     for(i in 1:length(sortedcor)){       if(sortedcor[i] > .5 | sortedcor[i] < -.5){         regvar[k] = i       k = k+1      }     }     regvar          sigvars <- names(sortedcor[regvar])      However, it is not working in my lm() function:               model1 <- lm(data.matrix(composite[1]) ~ sigvars, data = composite)      Error: Error in model.frame.default(formula = data.matrix(composite[1]) ~ sigvars, : variable lengths differ (found for 'sigvars')


title: Statistical podcasts

text: What are some podcasts related to statistical analysis? I've found some audio recordings of college lectures on ITunes U, but I'm not aware of any statistical podcasts. The closest thing I'm aware of is an operations research podcast The Science of Better. It touches on statistical issues, but it's not specifically a statistical show.


title: What statistical test to use for A/B test?

text: We have two cohorts of 1000 samples each. We measure 2 quantities on each cohort. The first one is a binary variable. The second is a real number that follows a heavy tail distribution. We want to assess which cohort performs best for each metric. There are plenty of statistical tests to choose from: people suggest z-test, others use t-test, and others Mann–Whitney U.   * Which test or tests should we choose for each metric for our case?    * What happens if one test suggests significant difference between cohorts and some other test suggests non-significant difference?


title: ROC curve and its function beginner

text: I have 3 features of a signal (example: amplitude, frequency, energy). I want to check which feature is the best to represent that particular signal. That signal is classified into two categories (good or bad). I used ROC to find the AUC. From my understanding: the bigger the AUC is the better it can discriminate between two classes (good:bad). Is this true? If this is true, then I pick the feature that has the biggest AUC to be used in a classifier. Am I on the right path? I appreciate any answer.


title: Performance of hierarchical clustering for binary data in R

text: I am trying to use Hierarchical clustering to see how well it performs in classifying a dataset which I previously know its true classification. I am new to clustering in general. I was able to draw the dendogram of the cluster, but I am not sure how I can measure the performance of this method. Here's a reproducible example of what I have:               y <- matrix(rnorm(40), 10, 4, dimnames=list(paste("item", 1:10, sep=""), paste("variable", 1:4, sep="")))  ## suitable data from clustering     class <- c(rep("Yes",5),rep("No",5))      x <- cbind(y,class) ## the complete dataset     plot(as.dendrogram(hr), edgePar=list(col="blue", lwd=2), horiz=F) ## a nice plot!      How can I decide whether this method was actually able to classify the data correctly? can I also split the data according to the branch it belongs to? I would like to mention that I am working on a big dataset, so looking by the naked eye on the dendogram isn't the best thing to do. Note: It might not be the best method, but I was trying to find a way that takes into account 10 variables to cluster the data. Many thanks,


title: Neural nets method? Am i doing something wrong?

text: To set the scene - i am using neural nets in SAS program i have ~ 1000 points variables have been selected prior to this on a combination of human knowledge and variable recombination (PCA etc) I am using 3 way 60-20-20 validation and 70-15-15 I have come up with a method that im sure has a flaw somewhere as it seems to work quite well! Previously i found it difficult to run two independant neural nets which are identical in every way and get comparable results. Also with no real fixed rules around neural nets except (from reading) \- 1 layer is usually ok \- usually you only need a number of neurons per layer equal to number of inputs Sometimes i get a good one, sometimes a bad one for a specified layer/neuron layout. I know this is a product of 'random' learning but if 1 or 2 out of 5 work it doesnt sit well with me. Also looking at a graphical representation of these one off good fits usually ( **but not always** ) produces very archaic and non real world fits when looking at graphical representations. So i came up with an idea where i can run X number of models of for any specified layer/neuron number arrangement. vary the number of neurons OR layers by 1 and rerun X times ..repeat. Effectively modelling all reasonable layouts This allowed me to collect many many (100s and 1000s) models with multiple iterations and look for consistency and well fitted models where by R^2s for validation, test training etc and other figures are auto tabulated this SEEMS to work. If i get a good and consistent bunch (say 3 sigmoidal neurons in 1 layer and 1 in the second) perform well in and consistently across 10 repeats i then look at those models graphical output for realistic curves As another option i can run 2 sets of these models for 2 different validation column layouts. so double the number of models from above. This in as one randomly generated 3 way validation result may by chance collect too many low results in the test section. I do see differences in differing 3 way validation One fear with this is am i crossing the boundary of refining after i have used a TEST set? or as i am running the models blind and then looking at test sets is this ok? The test, training and validation are all fixed for all models Also, does it even matter if 1 out of 10 models of same spec is good and the others are all bad? Are tehre any other blinding errors i am making? Appreciate any responses as it is a lengthy read


title: Correlation of likert scales for different questions

text: Can you correlate 2 Likert scales when the question for each one was different? I have 1 question which asks how concerned someone is about pollution. They rank this on a Likert Scale of 1-5 where 1 is most concerned and 5 not at all. Then my other question asks whether people intend to take any extra measures to reduce their own pollution in future, again on a scale of 1-5 but where 1 is very likely down to 5 not at all. Can I run a correlation test on these two questions? Naturally I would expect the more concerned one is, the greater their intention to act. Obviously both scales are 1-5 but as the meaning behind what the number means is slightly different (I.e in Q1, 1= Most Concerned and in Q2, 1= Very Likely), can they still be correlated?


title: RMSE for Model Validation and Calibration

text: Based on UK's Framework for Marine and Esuarine Model Specification,the % RMSE for current speed is calculated by dividing RMS by maximum peak current. $$\text{RMSE(%)} = \left(\frac{\text{RMSE}}{\text{Max current speed}}\right) \cdot100$$ Why maximum peak current is used and not mean average? Any logical reasons behind it?


title: Surrogate loss functions that are not consistent in the binary-class case

text: I'm looking for examples of surrogate loss functions that are useful for binary-class classification but are not consistent with respect to the zero- one loss, where I'm using consistency in the sense of surrogate loss functions converging to Bayes predictor [1] (also known as classification calibration, Bayes consistency and infinite-sample consistency). A classical result from [1], states that if $\psi$ is a convex surrogate loss function then $\psi$ is consistent with respect to the zero-one loss if and only if $\psi$ is differentiable at zero and $\psi'(0) < 0$. This seems to include all widely-used loss functions (i.e., hinge loss, logistic loss, exponential, etc.). I'm looking for examples of surrogate loss functions that are useful but which are not consistent in this setting. Any ideas? [1] P. L. Bartlett, M. I. Jordan, and J. D. McAuliffe, “Convexity , Classification , and Risk Bounds,” J. Am. Stat. Assoc., pp. 1–36, 2003. ↩


title: Testing residuals from a cox model with time dependent covariates

text: I'm doing survival analysis with time dependent covariates, using the counting process style. I already have a set of models and I want to test de residuals. I'm having trouble with the lack of information about this subject, **it makes sense to test Schoenfield residuals for time-dependent covariates**? I'm using R.


title: Decision tree for output prediction

text: I have satellite data that provides radiance which I use to compute the Flux (using surface and cloud info). Now using a regression method, I can develop a mathematical model directly relating radiance and flux and can be used to predict the flux for new radiance values. Is it possible to do same using decision trees or regression trees? In a regression there is mathematical equation connecting dependent and independent variable? Using decision trees, how could you develop such a model?


title: Can the sample space contain outcomes that do not occur?

text: The definition of the probability density function of the continuous uniform distribution on Wikipedia is: $$ f(x) = \left\\{ \begin{array}{} \frac{1}{b-a} & \text{for } a \leq x \leq b \\\ 0 & \text{for } x < a \text{ or } x > b \end{array} \right. $$ As usual, $\mathbb{P}(X=x) = 0$ for any $x \in [a,b]$. But this just mean that the it "almost surely" does not occurs, rather than it is "impossible". However, when $x < a$ or $x > b$, $\mathbb{P}(X= x) = 0$ and in this case, it actually is impossible. Values in intervals such as $[c,d]$ where $c > b$ never ever occurs either, rather than just "almost surely" not occur. I thought a nice way to resolve this would be to simply restrict the sample space to exclude outcomes that is impossible or never occurs. So the pdf would instead be: $$f(x) = \frac{1}{b-a} \text{ for } a \leq x \leq b.$$ This way, a probability of $0$ always mean that the event almost surely does not occur, which (to me) feels intuitively satisfying. I was wondering if this is already being done in measure theoretical probability theory? The definition of the sample space $\Omega$ that I have read always define it as the set of all _possible_ outcomes.


title: Best method of calculating line of best fit / extrapolate to compensate for delays

text: Let's suppose there is a project which is expected to take a certain amount of time to complete. As certain jobs are done, we can quantifiably measure how much of the project has been completed at any time. Consider a progress bar - it starts at 0%, and is updated at some discrete time values as parts of the project are completed, until it reaches 100%. It is expected that this completion percentage will increase linearly over time. While the project has yet to be completed, it seems logical that the best way of estimating the final completion date based on this data is to plot the (time,%) values, calculate a line of best fit, and find where that intersects the 100% value. As time goes on, these approximations should get closer and closer to the true value, oscillating around it on both sides. However, let's say some small part of the project (for example, a 5% range) gets delayed (or perhaps gets completed faster than expected). This period is going to drag down (or up) the line of best fit every single time it is calculated - meaning that our estimated completion date is going to be permanently too low, and will increase each time we calculate it, rather than oscillating about the true value like you'd expect. Is there a better method of estimating the completion time that can take into account local periods of faster or slower increase? Some sort of moving average may smooth the data, but doesn't seem logical with the variation not being cyclic. It is more like you'd want a piece-wise model - but how could you calculate the actual pieces and use this information appropriately without just looking at the graph and eyeballing where it 'seems' different?


title: (Why) do overfitted models tend to have large coefficients?

text: I imagine that the larger a coefficient on a variable is, the more ability the model has to "swing" in that dimension, providing an increased opportunity to fit noise. Although I think I've got a reasonable sense of the relationship between the variance in the model and large coefficients, I don't have as good a sense as to _why_ they occur in overfit models. Is it incorrect to say that they are a symptom of overfitting and coefficient shrinkage is more a technique for reducing the variance in the model? Regularization via coefficient shrinkage seems to operate on the principle that large coefficients are the result of an overfitted model, but perhaps I'm misinterpreting the motivation behind the technique. My intuition that large coefficients are generally a symptom of overfitting comes from the following example: Let's say we wanted to fit $n$ points that all sit on the x-axis. We can easily construct a polynomial whose solutions are these points: $f(x) = (x-x_1)(x-x_2)....(x-x_{n-1})(x-x_n)$. Let's say our points are at $x=1,2,3,4$. This technique gives all coefficients >= 10 (except for one coefficient). As we add more points (and thereby increase the degree of the polynomial) the magnitude of these coefficients will increase quickly. This example is how I'm currently connecting the size of the model coefficients with the "complexity" of the generated models, but I'm concerned that this case is to sterile to really be indicative of real-world behavior. I deliberately built an overfitted model (a 10th degree polynomial OLS fit on data generated from a quadratic sampling model) and was surprised to see mostly small coefficients in my model:               set.seed(123)     xv = seq(-5,15,length.out=1e4)     x=sample(xv,20)     gen=function(v){v^2 + 7*rnorm(length(v))}     y=gen(x)     df = data.frame(x,y)          model = lm(y~poly(x,10,raw=T), data=df)     summary(abs(model$coefficients))     #     Min.  1st Qu.   Median     Mean  3rd Qu.     Max.      # 0.000001 0.003666 0.172400 1.469000 1.776000 5.957000               data.frame(sort(abs(model$coefficients)))     #                                   model.coefficients     # poly(x, 10, raw = T)10                  7.118668e-07     # poly(x, 10, raw = T)9                   3.816941e-05     # poly(x, 10, raw = T)8                   7.675023e-04     # poly(x, 10, raw = T)7                   6.565424e-03     # poly(x, 10, raw = T)6                   1.070573e-02     # poly(x, 10, raw = T)5                   1.723969e-01     # poly(x, 10, raw = T)3                   6.341401e-01     # poly(x, 10, raw = T)4                   8.007111e-01     # poly(x, 10, raw = T)1                   2.751109e+00     # poly(x, 10, raw = T)2                   5.830923e+00     # (Intercept)                             5.956870e+00      Maybe the take-away from this example is that two thirds of the coefficients are less than 1, and _relative to the other coefficients_ , there are three coefficients that are unusually large (and the variables associated with these coefficients also happen to be those most closely related to the true sampling model). Is (L2) regularization just a mechanism to diminish the variance in a model and thereby "smooth" the curve to better fit future data, or is it taking advantage of a heuristic derived from the observation that overfiited models tend to exhibit large coefficients? Is it an accurate statement that overfitted models tend to exhibit large coefficients? If so, can anyone perhaps explain the mechanism behind the phenomenon a little and/or direct me to some literature? Thanks!


title: Limitations of the Nataf Transformation

text: I'm working with the Nataf Transformation (i.e a Gaussian copula) and am looking at bivariate joint distributions. However I can't find any publications discussing the limitations of the transform. For example it seems that some joint distributions cannot be completely describe by the transform and I'm trying to find the types of input distributions it can handle. Does anyone have any ideas? Or does anyone know of some good references that do clearly talk about the limitations? Many thanks.


title: Standard deviation of percentage changes in a time course

text: I ran into a problem analysing some data I produced. I measured concentrations in biological samples in a time course manner. I measured three replicates for each of the seven time points. I can display the data as means +- SD, that's easy. Now I want to express the data as percent changes relative to the first time point. I am completely clueless on how to apply SD or any other measurement of variability to the percentages. Sorry if this is a double post, but the other answers somehow don't fit to my problem. Another problem I will run into soon is the following: I will analyze a different genotype in the same manner and want to test for significant differences between the genotypes. Any recommendation for a test? Preferably something that can be done in Excel, I am just beginning to learn R.


title: How to fill in missing data series?

text: I have monthly water quality (Nitrate) measurements from July, 2005 till Oct., 2013 (84 measurements). However, there is no measurements from Oct., 2007 till June, 2009 (21 measurements). Could you please let me know how do I fill in the missing data using other dataset? What is the best software do you recommend to be used to fill in the missing data? ![enter image description here](http://i.stack.imgur.com/yjEhf.jpg) Thanks in advance!


title: Interpreting spline results

text: I'm trying to fit a spline for a GLM using R. Once I fit the spline, I want to be able to take my resulting model and create a modeling file in an Excel workbook. For example, let's say I have a data set where y is a random function of x and the slope changes abruptly at a specific point (in this case @ x=500).               set.seed(1066)     x<- 1:1000     y<- rep(0,1000)          y[1:500]<- pmax(x[1:500]+(runif(500)-.5)*67*500/pmax(x[1:500],100),0.01)     y[501:1000]<-500+x[501:1000]^1.05*(runif(500)-.5)/7.5          df<-as.data.frame(cbind(x,y))          plot(df)      I now fit this using               library(splines)     spline1 <- glm(y~ns(x,knots=c(500)),data=df,family=Gamma(link="log"))      and my results show               summary(spline1)          Call:     glm(formula = y ~ ns(x, knots = c(500)), family = Gamma(link = "log"),          data = df)          Deviance Residuals:           Min       1Q   Median       3Q      Max       -4.0849  -0.1124  -0.0111   0.0988   1.1346            Coefficients:                            Estimate Std. Error t value Pr(>|t|)         (Intercept)             4.17460    0.02994  139.43   <2e-16 ***     ns(x, knots = c(500))1  3.83042    0.06700   57.17   <2e-16 ***     ns(x, knots = c(500))2  0.71388    0.03644   19.59   <2e-16 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          (Dispersion parameter for Gamma family taken to be 0.1108924)              Null deviance: 916.12  on 999  degrees of freedom     Residual deviance: 621.29  on 997  degrees of freedom     AIC: 13423          Number of Fisher Scoring iterations: 9      At this point, I can use predict function within r and get perfectly acceptable answers. The problem is that I want to use the model results to build a workbook in Excel. My understanding of the predict function is that given a new "x" value, r plugs that new x into the appropriate spline function (either the function for values above 500 or the one for values below 500), then it takes that result and multiplies it by the appropriate coefficient and from that point treats it like any other model term. How do I get these spline functions? (Note: I realize that a log-linked gamma GLM may not be appropriate for the data set provided. I am not asking about how or when to fit GLMs. I am providing that set as an example for reproducibility purposes.)


title: ROC graph shape

text: Could you explain to me how the shape of a ROC curve is determined? From the following illustration, it seems that for every time the actual class (C) is positive, it goes up and when it's negative, it goes to the right. Why does it only seem to be based on the true class rather than the actual right or wrong classification of the object? I assume it's to do with the mentioned score, but I don't know how this would be calculated and ranked. Please explain. ![](http://habrastorage.org/files/267/36b/ff1/26736bff158a4d82893ff85b2022cc5b.gif) Thank you.


title: Stratified Sampling and the Central Limit Theorem

text: What can be said about the convergence rate of stratified sample means to a normal distribution, given different allocation schemes? Obviously, under very poor allocation, this convergence can fail (for instance, keeping the sample size in one stratum fixed); but for any reasonable scheme, and under reasonable hypotheses, a CLT applies (Bickel and Freedman 1984, etc). However, is anything known about the sample sizes required to obtain close approximations to a normal distribution, given, say, a reasonable allocation scheme, and information about the (sample size needed for) normal approximation under simple random sampling from the distribution? Or at least for proportional and/or optimal allocation?


title: Using the Lagrange multiplier statistic in regression

text: I don't understand how the LM statistic works because if we run the regression and then record the residual, and then run the regression using the residual as the dependent variable and then using the unrestricted model, how can we run a regression of a constant on 5 different variables? For example, I'm using Stata and here's my commands:               reg bwght cigs parity faminc      where `bwght` is my dependent variable and `cigs`, etc. are the dependent variables               scalar epsilon=e(rss)      here is where I saved my error term (I think)                reg epsilon cigs parity faminc motheduc fatheduc      I tried to run this as a regression but I was told that there was an error because epsilon isn't a variable. Which I can understand why. Which leads to my question. How can we conduct a Lagrange Multiplier test?


title: Mixed effect logistic regression: defining random effect for repeated measure design with replications

text: I measured a binary response for each subject in 5 different conditions. For each subject and condition, I replicated the experiment 36 times. I thus have 36 binary values per condition per subject. I am trying to build a model for those data. I suppose a logistic regression is what I'm looking for, and I am working with the `lmer` package. My aim is to check whether the conditions significantly influence the observed values, so I would have two models:               lmH1<-lmer(value~condition, (random effects), data=dataset, family=binomial)      and               lmH0<-lmer(value~1, (random effects), data=dataset, family=binomial)       By looking at the output from `anova(lmH0, lmH1)`, I would be able to determine the significance of the effect of my condition. I am just not sure what to specify as random effect; the models I defined so far are:               lmH1 <- lmer( value ~ condition + ( 1 | subject ), data = dataSet, family = binomial )      and               lmH2 <- lmer( value ~ condition + ( 1 | subject/condition ), data = dataSet, family = binomial )      However I am not sure about how lmer handles the replicates, so I don't know whether I should include those replicates in my random effects or not. I could modify the proposed models so that the grouping defined by the random effects refers to a specific binary values instead of a group of binaries values. My new models would then be               lmH1a <- lmer( value ~ condition + ( 1 | subject/(condition:replicate) ), data = dataSet, family = binomial )      and               lmH2a <- lmer( value ~ condition + ( 1 | subject/condition/replicate ), data = dataSet, family = binomial )      With those models R returns the warning message `Number of levels of a grouping factor for the random effects is equal to n, the number of observations`. But the model is still computed. All 4 models return very similar values for the fixed effects and for the random effects that they have in common (e.g. the subject random effects are very similar for all 4 models and the condition within subject random effects are very similar for `lmH2` and `lmH2a`). How can I check which random effect structure is the most appropriate for my design and collected data?


title: What is acceptable mean value for a Likert scale?

text: I have a question about mean value in a Likert scale. My objective is to identify important factors for selecting short or long term relationship with subcontractor. I used Likert-type scale ranging from "highly agree" to "highly disagree" on 22 factors. Then I have asked 41 main contractors and I used that result to analyse mean and standard deviation of all factors. After I got mean and SD I would like to ask whether and how I can cut some factors of my 22 factors by mean value. If there is an accepted method, please tell me too. And if you can provide reference on that method it would be appreciated. Sorry that make u confuse. Actually, my purpose of study is to identify important factors for selecting short or long term relationship with subcontractor. the meaning of " how can i cut some factors of my 22 factors by mean value" means i want to reduce the number of factor to less than 22 factor by using result of mean value. for example, below result is the result of mean and SD in 22 factors, so after the result, i want to find important influencing factors by mean value. for me, i want to take (mean value > and = 4) so all the factors that mean value > and = 4 i will take as important influencing factor. anyway am i make u clear? if so , is this method accepted (mean > and= 4) is ok? or if you have any method pls share me. one more thing, if u a agree with my thinking, do u have any reference to support cos now i am trying to find reference to support too. thanks you Factors Mean SD Time control in planning 4.537 0.596 Work quality 4.463 0.552 Cooperation 4.268 0.549 Experience 4.220 0.725 Resources 4.171 0.587 Commitment 4.146 0.527 Honesty 4.122 0.557 Monitoring 4.122 0.510 Coordination 4.049 0.631 Trust 4.024 0.724 Joint problem solving 4.024 0.689 Clear understanding 4.000 0.632 Safety training for employees 3.976 0.851 Profit base 3.976 0.612 Communication 3.976 0.651 Employee skill training 3.951 0.545 Flexibility to change 3.927 0.608 Innovation 3.878 0.678 Safety control system 3.854 0.691 Wastage disposal control 3.732 0.633 Financial status 3.683 0.820 Knowledge 3.683 0.567


title: Why would someone use a Bayesian approach with a 'noninformative' improper prior instead of the classical approach?

text: If the interest is merely estimating the parameters of a model (pointwise and/or interval estimation) and the prior information is not reliable, weak, (I know this is a bit vague but I am trying to establish an scenario where the choice of a prior is difficult) ... Why would someone choose to use the Bayesian approach with 'noninformative' improper priors instead of the classical approach?


title: Problem understanding what type of test to use and how to proceed with the questions provided

text: I have a case study where a company made new golf balls with stronger coatings than current golf balls. One of the technicians is concerned that the new coating may decrease the driving distance. An experiment was conducted in which 40 types of each ball was hit by a machine to determine mean driving distance of the two balls. The distances for each ball were given in a spreadsheet. Here are the samples:               Current: 264.00 261.00 267.00 272.00 258.00 283.00 258.00 266.00 259.00 270.00 263.00               264.00 284.00 263.00 260.00 283.00 255.00 272.00 266.00 268.00 270.00 287.00               289.00 280.00 272.00 275.00 265.00 260.00 278.00 275.00 281.00 274.00 273.00               263.00 275.00 267.00 279.00 274.00 276.00 262.00     New;     277.00 269.00 263.00 266.00 262.00 251.00 262.00 289.00 286.00 264.00 274.00               266.00 262.00 271.00 260.00 281.00 250.00 263.00 278.00 264.00 272.00 259.00               264.00 280.00 274.00 281.00 276.00 269.00 268.00 262.00 283.00 250.00 253.00               260.00 270.00 263.00 261.00 255.00 263.00 279.00      These are the questions:   1. Formulate and present the rationale for a hypothesis test that Par could use to compare the driving distances of the current and new golf balls   2. The CEO tells you (in contrast to what Bill said) that he would like to sell the new ball unless there is overwhelming evidence that the cut‐resistant ball is slower than the old ball so you should frame your hypothesis accordingly. Hint: you are more willing to make a Type II error. Analyze the data to provide the hypothesis testing conclusion and statements (at .05). Also include narrative versions of what ever of the above forms you use such as: Ho: The mean is ... What is the p‐value for your test? What is your recommendation for Par, Inc.?   3. Provide the following descriptive summaries of the data for each model: ONLY the means and variances and the z values.   4. What is the 95% confidence interval for the difference between the means of the two populations?   5. Do you see a need for larger sample sizes and more testing with the golf balls? Discuss. I've looked for help elsewhere and everyone seems to use the t-test to conduct the analysis. Shouldn't the z-test be used since the sample sizes are both over 30? I would also like some direction in understanding how to answer the other questions but my main question is the most important.


title: One class classifier Cross validation

text: I am working on a problem which requires one-class classifier. I am using LIBSVM. I know there are tons of material out there but still I could not find the answer to my query.   1. How do I estimate the optimum parameters for the RBF kernel?   2. Using `svm-train -s 2 -t 2 -v 5 train.scale`, I am getting 49% accuracy and my training set has no outlier data. So, does this actually mean 98% accuracy in real world scenario?


title: Spacings between discrete uniform random variables

text: Let $U_1, \ldots, U_n$ be $n$ i.i.d discrete uniform random variables on (0,1) and their order statistics be $U_{(1)}, \ldots, U_{(n)}$. Define $D_i=U_{(i)}-U_{(i-1)}$ for $i=1, \ldots, n$ with $U_0=0$. I am trying to figure out the joint distribution of $U_i$'s and their marginal distribution and possibly their first few moments. Can anyone give some hint on this. Also can you please recommend a book on order statistics? Thank you.


title: Tool to work with dataflow easier/more visual

text: I've written some Python code which does some (network) data analysis. It's mostly large data coming from SQL where I play around with statistics and algorithms. However configuring the interplay between modules and tracking/debugging what's going on with the data isn't very convenient. Do you know a tool which could nicely display the dataflow modules, with maybe some tooltip information, some table data browser etc.? Something where I could incorporate my Python code and display table-like data.


title: Latent Dirichlet analysis document comparison

text: I have been using Latent Dirichlet Analysis for a while but I am a bit confounded as to it's practical ability to compare two documents. It is of course ideal for classification when you want to see in what category a certain document or word belongs, but in content comparison I am at my wits end. An inferred document yields a distribution over n topics, summing to one. Using Kolmogorov-Smirnov was recommended somewhere, but it has the annoying property of giving two identical distributions a very high score, even if they do not make sense. Two words not in the dictionary will give a perfect match, which is of course absurd. The normalized dot product actually works very well because it punishes flat distributions (trivial documents), but I thought there would be a better one. Any suggestion is appreciated.


title: Spline df selection in a general additive Poisson model problem

text: I've been fitting some time series data using a Poisson general additive model using SAS's `PROC GAM`. Generally speaking, I've been having it's built-in generalized cross-validation procedure generate at least a decent "starting point" for my single spline, which is a non-linear function of time along with a single parametric term (the one I'm actually interested in). So far, it's worked rather swimmingly, with the exception of one of my data sets. There are 132 observations in that data set, and GCV suggests a spline of 128 degrees of freedom. That seems...wrong. Very wrong. More importantly, it's also not at all stable. I tried a second approach, using something like a "Change in Estimate" criteria to stop adding degrees of freedom when the estimate of the parametric term stops changing because why continue to add control if nothing's different? The problem is that the estimate isn't at all stable. I tried the following degrees of freedom, and as you can see, the parametric term bounces around wildly:               DF: Parametric Estimate:     1   -0.76903     2   -0.56308     3   -0.47103     4   -0.43631     5   -0.33108     6   -0.1495     7    0.0743     8    0.33459     9    0.62413     10   0.92161     15   1.88763     20   1.98869     30   2.5223     40-60 had convergence issues     70   7.5497     80   7.22267     90   6.71618     100  5.83808     110  4.61436     128  1.32347      I've no intuition at all about what I should be using in terms of df for this particular bit of data. Any other ideas for how to choose a df? Should I be looking at the significance of the spline? Doing some more looking between df = 10 and df = 15, it looks like df = 12 is the closest you can come to the estimate generated by 128 and still be in the "reasonable degrees of freedom" range. Along with the linear term, the intercept and the single parametric term, that feels like a pretty heavily saturated model. Is it justifiable to just go with 12? As a second update, changing the smoothing from `spline(t)` to `loess(t)` is resulting in much more well-behaved df estimates - should I just switch to loess smoothing?


title: " random effects must be less than the number of observations" error in lmer package

text: I'm trying to implement a regression model with both fixed and random effects. The package I use is the `lme4`. I want to find the relationship between the continuous variables `X1` and `X2`, using the categorical variable `F1` as a fixed effect. When I try to implement the model `model <- lmer(X1 ~ F1 + (1|X2))` I get the following error: > `Number of levels of a grouping factor for the random effects must be less > than the number of observations` My understanding for the error is that the `X2` varies less than the `X1` (which has some similar values). Why is this a problem?


title: What are statistical methods for comparing different brackets

text: I am not particularly familiar with statistics and am looking at methods for analysing numbers that have been broken into different "brackets" or "groups" for different entities. Consider three companies.   * A announces that they have 10 factories generating between 10,000 and 20,000 widgets per year and 15 factories generating between 20,000 to 40,000    * B announces that they have 7 factories generating between 15,000 and 25,000 widgets and 5 factories generating between 25,000 and 35,000   * C announces they have 5 factories generating between 10,000 and 23,000 widgets and 14 between 23,000 and 40,000 To be clear all widget production is included in my simplified example above. What statistical methods can be used (and do any exist) that would allow you to aggregate this data together and make generalised predictions (all companies have on average x factories that generate 15,000 widgets per year and y factories that generate 40,000)?


title: How can I find tendencies in specific cases between two groups?

text: I'm doing this project where we are investigating whether social networks can influence peoples distinctiveness. Briefly summarized, we have primed a group with a social media and another is the control group. Afterwards, they were given a questionnarie where they have...   1. Written 12 identity elements 2. Rated these elements on how central they are to their identity on a likert-scale from 1-7 where seven is very central to their identity (this scale is used throughout) 3. Rated these elements on how a) different they make them b) how it constitutes a special social role and c) how seperated these elements make them. We are not using the words further, just their ratings. So far so good. We have alot of variables in our database: 1: Age 2: Sex 3: Priming or control-group 4: Centrality of the elements (the mean of 2 questions about the centrality asked in different ways) (12 variables) 5: How different it makes them (12 variables) 6\. How "roly" it makes them (12 variables) 7: How separated it makes them (12 variables) The way i have tackled this so far: Found the mean of centrality for each element, multiplied it by variable 5-6-7 (look up) casewise and squarerooted the number (gives a 1-7 number). Then the mean of these three new variables (CentralityXDifference, CentralityXRole, CentralityXSeperateness) gives us a number of "Total distinctiveness". Unfortunately these means don't differ depending on priming or controlgroup, but that's okay. We can probably argue why. Here comes some questions   1. What do i do if I want to check for tendencies that the priming or controlgroup accordingly have high scores in variable 5-6-7 when they have high in variable 4; Do the priming group tend to rate elements of high centrality as highly distinctive? I cannot figure out a way to do this in a smart way, so if you have an idea i'll be grateful.   2. Have I done something wrong the way i've done it already? Is there room for improvement in my method?


title: Interaction of 2 categorical variables with an empty cell

text: If I have a logistic regression equation in the form of : Logit($\pi$($X_1$, $X_2$)) = $\beta_0$ + $\beta_1$$X_1$ + $\beta_2$$X_2$ + $\beta_3$$X_1$$X_2$ where $X_1$ and $X_2$ are both binary categorical variables and there is an empty cell in their contingency table, what happens to $\beta_3$? I am assuming that there will be no convergence when trying to calculate the MLE for $\beta_3$ similar to when there is an empty cell in the contingency table of a dependent and independent variable. Could someone elaborate on this?


title: How to get scored combination of features

text: My data looks like this (F=Features)                       F1     F2     F3     F4     F5     F6     F7     F8....     ID1     0.67   0.76   0.3    0.54   0.21   0.88   0.97   0.45....     ID2     0.76   0.68   0.10   0.45   0.12   0.44   0.79   0.54....     ID3     0.67   0.76   0.3    0.54   0.21   0.88   0.68   0.76....     ID4     0.67   0.10   0.3    0.45   0.3    0.88   0.97   0.45....     ...     ...     ...      I have about 40 features (I have just put 8 here). If I set a threshold for, say 4 features, what I am looking for is a combination of 4 features which, together, are most relevant and significant in the dataset. I need some kind of score that measures how good a combination of features is. This is what I meant by the confidence score (or whatever score we may call it). So instead of selecting 1 feature I want to select a combination of 4 features. For example.               F1-F4-F9-F12 = 0.92     F2-F3-F7-F6 = 0.85     F5-F3-F4-F8 = 0.667      Here, F1-F4 is not subtraction. I am just attaching the two features together. The scores above, I do not know how to get it. That is my question? So I do not know what kind of a test to use or can be used. How can I go about it? Thanks


title: What's the probability that a bookmaker is mispricing odds on soccer games?

text: An English soccer team plays a series of matches against different opponents, of varying ability. A bookmaker offers odds for each match as to whether it will be a home win, away win, or draw. Part-way through the season, the team has played n matches, and has drawn $k$ of them, which is more than might be expected from the odds. What is the probability that the bookmaker is mis- pricing the odds on these matches, rather than just being unlucky? If the bookmaker continues to price the team's remaining matches in a similar way, and I bet $1 that each one will be a draw, what is my expected return?


title: Independence of real and imaginary parts of the scalar product of multivariate complex normal vectors

text: Reading the paper _"Distribution of Inner Product of Two Complex Gaussian Vectors and its Application to MPSK Performance" by Mallik R.K._ I can not get one of the steps.   The author assumes two complex Gaussian vectors $$\vec{X}\sim\mathcal{CN}(\vec{\mu_X},\mathrm{I}), \\\ \vec{Y}\sim\mathcal{CN}(\vec{\mu_Y},\mathrm{I}),\\\ \mathrm{E}\left[\Re(x_i)\Re(x_j)\right]=\mathrm{E}\left[\Re(y_i)\Re(y_j)\right]=\delta_{ij},\\\ \mathrm{E}\left[\Im(x_i)\Im(x_j)\right]=\mathrm{E}\left[\Im(y_i)\Im(y_j)\right]=\delta_{ij},\\\ \mathrm{E}\left[\Re(x_i)\Re(y_j)\right]=\mathrm{E}\left[\Im(x_i)\Im(y_j)\right]=0. $$ and is after the joint distribution of real and imaginary parts of their inner product:   $$Z=\vec{Y}^H\\!.\vec{X}, \\\\\ Z_1=\Re(Z),\\\ \ Z_2=\Im(Z)$$ It is said in the paper that > ... it is clear that $Z_1$ and $Z_2$ conditioned on $Y$, are independent But then I went on and found the covariance between $Z_1, Z_2$ and got a nonzero result, so they can not be independent. Or am I wrong somewhere in my derivations?   Here they are.   Starting from the definition: $$\require{cancel} \begin{eqnarray} &&\mathrm{cov}(z_1,z_2)= \mathrm{cov}\left(\sum_i(\Re(y_i)\Re(x_i)+\Im(y_i)\Im(x_i)),\sum_i(\Re(y_i)\Im(x_i)-\Im(y_i)\Re(x_i))\right)=\\\ &=&\mathrm{E}\left[\left(\sum_i(\Re(y_i)\Re(x_i)+\Im(y_i)\Im(x_i))\right)\\!\\!\left(\sum_i(\Re(y_i)\Im(x_i)-\Im(y_i)\Re(x_i))\right)\right]-\\\ &-&\mathrm{E}\left[\sum_i(\Re(y_i)\Re(x_i)+\Im(y_i)\Im(x_i))\right]\\!\\!\mathrm{E}\left[\sum_i(\Re(y_i)\Im(x_i)-\Im(y_i)\Re(x_i))\right]=\\\ &=&\mathrm{E}\left[ \sum_i\sum_j\Re(y_i)\Re(y_j)\Re(x_i)\Im(x_j)\right]+\mathrm{E}\left[\sum_i\sum_j\Im(y_i)\Re(y_j)\Im(x_i)\Im(x_j)\right]-\\\ &-&\mathrm{E}\left[\sum_i\sum_j\Re(y_i)\Im(y_j)\Re(x_i)\Re(x_j)\right]-\mathrm{E}\left[\sum_i\sum_j\Im(y_i)\Im(y_j)\Im(x_i)\Re(x_j) \right]-\\\ &-&\left(\sum_i (\Re(y_i)\mathrm{E}\left[\Re(x_i)\right]+\Im(y_i)\mathrm{E}\left[\Im(x_i)\right])\right)\left(\sum_i (\Re(y_i)\mathrm{E}\left[\Im(x_i)\right]-\Im(y_i)\mathrm{E}\left[\Re(x_i)\right])\right)=\\\ &=& \sum_i\sum_j\Re(y_i)\Re(y_j)\cancelto{0}{\mathrm{E}\left[\Re(x_i)\Im(x_j)\right]}+\sum_i\sum_j\Im(y_i)\Re(y_j)\cancelto{\delta_{ij}}{\mathrm{E}\left[\Im(x_i)\Im(x_j)\right]}-\\\ &-&\sum_i\sum_j\Re(y_i)\Im(y_j)\cancelto{\delta_{ij}}{\mathrm{E}\left[\Re(x_i)\Re(x_j)\right]}-\sum_i\sum_j\Im(y_i)\Im(y_j)\cancelto{0}{\mathrm{E}\left[\Im(x_i)\Re(x_j)\right]}-\\\ &-&\left(\sum_i (\mu_{\Re(x)_i}\Re(y_i)+\Im(y_i)\mu_{\Im(x)_i})\right)\left(\sum_i (\Re(y_i)\mu_{\Im(x)_i}-\Im(y_i)\mu_{\Re(x)_i})\right)=\\\ &=&\\!\sum_i\\!\sum_j\left(\\! \Re(y_i)\\!\Im(y_j)\\!(\mu_{\Re(x)_i}\mu_{\Re(x)_j}\\!\\!-\\!\mu_{\Im(x)_i}\mu_{\Im(x)_j})\\!\\!+\\!\mu_{\Re(x)_i}\mu_{\Im(x)_j}(\Im(y_i)\Im(y_j)\\!-\\!\Re(y_i)\Re(y_j)\\!)\\!\\!\right) \end{eqnarray} $$ And the last term is clearly not zero (unless the mean is zero).   So where's the mistake?


title: Generalization of Heckman's two step procedure

text: Following my question Heckman sample selection vs. OLS about the meaning of Mills ratio I am wondering why some researchers estimate a generalization of Heckman instead of the actual Heckman's procedure. (based on Das et. al 2003 http://www.jstor.org/stable/3648610 (gated version). Instead of using Mills ratio only, the generalized version also includes predicted probabilities, Mills ratio, their squares and interaction terms in the 2 step. I do no understand the advantages / disadvantages of this approach compared to "basic" Heckman sample selection? Which model should I prefer for testing sample selection?


title: Are there useful distributions for ternary variables (e.g. $-1,0,1$ data)?

text: The title says it. If one wishes to analyze ternary outcomes—that is categorical outcomes with specifically three values (-1,0,1)—are $\chi^{2}$ / contingency table tests and multinomial-logistic regression type models the only choices for making inferences, or are there well-articulated and useful ternary distributions and related inferential methods that we might consider as an alternative? Such data result in the differences of two 1/0-coded nominal variables, in the qualitative (sign) representation of causal effects (e.g. in loop analysis), in ternary logic, and perhaps elsewhere.


title: How to handle skewed binary target variables?

text: > **Possible Duplicate:**   >  Supervised learning with "rare" events, when rarity is due to the large > number of counter-factual events I am trying to predict diabetes using the BRFSS dataset by using a supervised learning classification model. But I see that the target variable which is having diabetes or not is skewed. That is 90% of the records are non-diabetic and only 10% of the records are diabetic. How do I handle the skewness in the target variable?


title: Linear vs. nonlinear regression

text: I have a set of values $x$ and $y$ which are theoretically related exponentially: $y = ax^b$ One way to obtain the coefficients is by applying natural logarithms in both sides and fitting a linear model:               > fit <- lm(log(y)~log(x))     > a <- exp(fit$coefficients[1])     > b <- fit$coefficients[2]      Another way to obtain this is using a nonlinear regression, given a theoretical set of start values:               > fit <- nls(y~a*x^b, start=c(a=50, b=1.3))      My tests show better and more theory-related results if I apply the second algorithm. However, I would like to know the statistical meaning and implications of each method. Which of them is better?


title: Determining distribution of a set of ordered numbers

text: I've always wondered this. Let's say someone gives you a spreadsheet with a large set of numbers, say 10,000, and says "What distribution does this set of numbers follow (Gaussian, exponential, Poisson, etc.)"? How do you go about determining the answer?


title: Help with computing the Expected Opportunity Loss

text: I'm trying to compute the expected opportunity loss (EOL) for a simple investment problem, where the cost **_C_** is set to be exactly 5 million and the value **_V_** is estimated to be normally distributed between 2.5 million and 25 million (with 90% CI). I understand that the EOL is simply _the cost of making the wrong decision_ x _the probability of making the wrong decision_. Here the default decision is to invest. In such case, I assume the EOL can be computed by:   1. chopping the value distribution into small segments   2. for each such segment from the interval (0, _C_ ] computing its mean and its probability   3. multiplying the mean of each segment by its probability    4. summing up all the products from step 3 Now for this particular case, I computed the EOL to be 246,000. However, this is different from the EOL that was computed for the same problem in the literature that uses a mysterious formula for the computation. I am simply looking for a confirmation that either:   1. the algorithm is wrong   2. the algorithm is correct, but the EOL result is wrong, in which case I need to revise my code   3. both the algorithm and my EOL result are correct, which means that the literature is wrong


title: Confidence interval of the expected value

text: Is this possible to compute? From a normal distribution we took a sample $( 11, 13, 13, 14, 14, 14, 15, 15, 17, 18).$ Compute the $98 \%$ confidence interval of the expected value. I'm studying statistics on my own but I don't know what formula should I use here.


title: Adjust p-values for variance inflation factors

text: I'm fitting linear mixed models and I have slightly inflated p-values due to multicollinearity. I deleted the factors with the highest VIFs until none of them was larger than 3. The VIFs tell me, how much the variances in my model are inflated due to collinearity, so just out of curiosity, I then calculated, what p-values I would have gotten with a VIF of 1. From the SE in the model summary I calculated the variance, divided it by the VIF and then re-calculated the new SE and t-value from which I then calculated "adjusted" p-values. Would you regard this as a valid approach? I did not find any publications on that topic and I'm afraid I'm making things a little too easy?!? I added some self-containing Code below:               # create data.frame with "inflated" t- and p-values              summary = data.frame(Factor = c("(Intercept)", "temperature", "color", "amount of beer"), Value = c(-3.2,2.1,0.5,3), Std.Error = c(3,1,0.8,0.9), DF = c(20,20,20,20))         summary$t.value = NA             summary$t.value = summary$Value/summary$Std.Error         summary$p.value = NA              for (i in 1:4)         {           p.value = 2*pt(-abs(summary$t.value[i]), df = 20)               summary$p.value[i] = p.value         }            # some moderate variance inflation factors          vif = c(1.4,3,1.03,2.4)          # calculate new t- and p-values taking VIFs into account         summary$var_new = sqrt(summary$Std.Error)/vif # variance / VIF = "adjusted" variance         summary$t.value_new = summary$Value/(summary$var_new)^2 # "adjusted" t-values             summary$p.value_new = NA              for (i in 1:4)  #calculating p-values from adjusted t-values         {           p.value_new = 2*pt(-abs(summary$t.value_new[i]), df = 20)               summary$p.value_new[i] = p.value_new         }


title: Modeling prior probability as a delta function

text: I'm using approximate Bayesian computation to find the true value of a parameter. My prior distribution is uniform over $(0, 1)$. I was watching this video on Bayesian learning and the lecturer states (around 36:00) that this is making a huge assumption. That is, we are assuming the value of the unknown variable has the same mean as the random variable (in my case $0.5$) and also has a variance that we can compute (I assume based off the number of observations?). He goes on to say that if you want to really model a prior like this, you should use a delta function centered around a value, $a$, which is unknown. My questions are these:   1. What is wrong with assuming that the mean and variance are the same as the random variable?   2. How can I use a delta function to generate probabilities for my prior?


title: Independence between residual sum of squares and mean?

text: Let $Y_1, \cdots, Y_n$ be independent random variables each with the distribution $N(\mu, \sigma^2)$, then \begin{align*} S^2 &= \sum(Y_i - \overline{Y})^2 \\\ &= \sum((Y_i - \mu) - (\overline{Y} - \mu))^2 \\\ &= \sum((Y_i - \mu)^2 + (\overline{Y} - \mu)^2 - 2(Y_i - \mu)(\overline{Y} - \mu)) \\\ &= \sum(Y_i - \mu)^2 + \sum(\overline{Y} - \mu)^2 - 2(\overline{Y} - \mu)(n\overline{Y} - n\mu) \\\ &= \sum(Y_i - \mu)^2 - n(\overline{Y} - \mu)^2 \end{align*} Dividing by $\sigma^2$ on both sides we can see that $S^2 \sim \chi^2(n-1)$, some books say from this you can declare $S^2$ and $\overline{Y}$ independent, why?


title: Should I use weighted data when analysing a population survey subsample

text: I have a question related to using data weightings when looking at a subgroup of a population sample. I am using survey data which was provided with weights, but I only want to include the taller participants in my analysis. When I have selected my subsample,the data weights obviously won't sum to 1. Can I just rescale the weights to sum to 1? Or do I need to recalculate the weights in some other way. The data is in SPSS. Should I be using the complex sample module?


title: Comparison of horizon forecasted values for two models (ARIMA and NAR)?

text: I fitted ARIMA and NAR models to 10 different univariate time series. These time series are water characteristics at a water treatment plant (CO$_2$, TDS, TH, pH, ...). After fitting suitable models I forecasted seven lags ahead with models for each series. I have seven observed at the end of each series that should be plotted with the corresponding seven forecasted. After plotting I found that the forecasted values by NAR model are visually near the observed but the values by ARIMA were almost constant across observed values. For more comparison I decided to calculate MSE for the forecasted values for two models. My expectation was that MSE for NAR should be smaller than ARIMA but this was completely different. For example, for one of my series forecasted values by NAR were visually close in comparison with those from ARIMA but its MSE was more than for ARIMA. Why? Please let me how can I compare them statistically.


title: Stack of HMMs. Does this have a name?

text: Consider a sequence of tokens. One could fit a Markov model to this sequence by computing the empirical stochastic transition matrix. A refinement to this model would be to build a HMM on top of it, with a small tweak: we let the emission probability of a token depend both on the hidden state and on the previous token. This is still a hmm but with a structure constraint. We then take the Viterbi path of the hidden states, and fit an HMM on top of it. And so on, with many HMM on top of HMMs. Of course, the end result is still a HMM, but hopefully one with a lot of structure and easy to fit. One example would be a natural language text. The first layers of HMM would describe common phonemes, higher levels common suffixes and prefixes, then words, parts of speech and sentence structure. This is _not_ a hierarchical HMM. Does this have a name?


title: How does the test error pattern depend on the regularizer function?

text: This question is regarding the role of regularizer in an objective function. Given a loss function $f(x)$, a regularizer function $r(x)$, and $\lambda$ being a trade-off function, our aim is to $\min f(x) + \lambda*r(x)$. I'm able to understand that as the value of $\lambda$ increases our training error increases. I wanted to know how the test error behaves as $\lambda$ increases. From "Pattern Recognition and Machine Learning" by Bishop, I came to know that "the test error initially decreases and then starts increasing". Consequently, we choose the $\lambda$ which has the lowest validation error. Does the same pattern (decrease and then increase) repeat for any given regularizer function? Or, does this pattern depends on the regularizer function which we use?


title: Extension of conditional probability result for non independent random vectors

text: Let $X$ & $Y$ be random vectors. Let $Z$=$f(X,Y)$ be a random variable. Then if $X$ & $Y$ are independent, the following is a well known result: $P(f(X,Y)\in B\mid Y=y)$ = $P(f(X,y)\in B)$ _**I need some help to 'formalize' what appears to be an intuitive generalization when $X$ & $Y$ are not assumed to be independent, namely that_** $P(f(X,Y)\in B\mid Y=y)$ = $P(f(X,y)\in B\mid Y=y)$ without assuming $Y$ is discrete. One 'immediate' problem for example is that, I am not sure that the RHS of the above equation is well defined. _Breiman's text 'Probability'_ , defines the conditional probability $P(C\mid Y=y)$ ( _Definition 4.7_ ) as a measurable function in $y$ satisfying $P(C, Y\in A )$=$\int_{A} P( C\mid Y=y) \, P_{Y}(dy)$ where $C$ is a measurable subset of the sample space and $A$ is an arbitrary Borel subset in the state space of $Y$. _**Using the concept of Regular Conditional Probability, I think the statement above can be formalized as follows._** Let $X$ & $Y$ be random vectors. Let $Z$=$f(X,Y)$ be a random variable, and $P_{Y}$ be the marginal distribution of $Y$. Let $Q_{X\mid Y}(\cdot\mid y)$ be a regular conditional distribution for $X$ given $Y=y$. Then $P(f(X,Y)\in B\mid Y=y)$ = $P(f(X,y)\in B\mid Y=y)$, is I believe an **_intuitive interpretation_** of the following 'unproven' statement $Q_{Z\mid Y}( B\mid y)$= $Q_{X\mid Y}(B_y \mid y)$ where $B$ is any arbitrary Borel subset of $R$, and $B_y$ = $[x \mid f(x,y) \in B]$ and $Q_{Z\mid Y}(\cdot\mid y)$ the regular conditional distribution for $Z=f(X,Y)$ given $Y=y$. Written alternatively we need to prove that $P(Z\in B, Y\in A )$=$\int_{A} Q_{Z\mid Y}( B\mid y) \, P_{Y}(dy)$=$\int_{A} Q_{X\mid Y}( B_y\mid y) \, P_{Y}(dy)$ So to conclude my question is how to prove the above statement about Regular Conditional Probabilities, and if the generalization to non independent vectors is correct. _**FYI I posted a similar question on MathStackExchange and have not got a response._**


title: Understanding the mean shift algorithm with Gaussian kernel

text: I am referring here to the Wikipedia article regarding the mean shift algorithm. Here the Gaussian kernel is denoted as $e^{-c\|x_i - x\|^2 }$ My questions are the following:   * Is $c$ always positive? I.e. are hence the "weights" given by the kernel always between 0 and 1?   * Looking at a picture explaining the algorithm, I have seen that in the two-dimensional case, the kernel can be understood as a circle..how do I see it from the formula above that it is indeed a circle, i.e. only points in this circle are considered at each iteration? Thanks


title: Confirmatory factor analysis for a newly constructed questionnaire

text: I am about to develop a new questionnaire with 142 questions I have performed exploratory factor analysis and 30 latent factors were extracted. I was wondering if it is **necessary** to perform **confirmatory factor analysis** for my newly constructed questionnaire too.


title: Relation between best fit line and eigenvector of maximum eigen value of an estimated covariance matrix

text: ( _This question is from my pattern recognition course._ ) There is this exercise: Imagine we have $N$ samples with $n$ dimensions. First it asks to find a point $m$ where the summation of Euclidean distances from $m$ is minimum. Then imagine another vector $e$ with size of 1. $e$ passes from $m$. Every point on this line crosses from: $x=m+\alpha*e$. The $\alpha_k$ value is the distance of a point on the line where the distance from that point and the $x_k$ is minimum. Then the exercise asks me to find values of $\alpha_k$ where the distance is minimum (i.e., the dashed line). The last part wants me to prove that the desired values of $\alpha_k$ are actually the eigenvector with the maximum eigenvalue of the below estimation of covariance matrix: $\Sigma=1/N\sum_{k=0}^{k=N} (x-m)(x-m)^t $ ![enter image description here](http://i.stack.imgur.com/rljb0.jpg)


title: algorithm for finding similar items

text: I have a group of items S, say there are 100 items in S, and I know some features of these items, i.e., color, size, ...Now, I have another group of items P, say there are 10000 items in P. What can I do to find a sub group of items in P that are similar to items in S? I know the corresponding features of all items in P, i.e., color, size... My original thinking is that I can use K nearest neighbor. But it seems that as KNN uses a majority vote, in order to use KNN, at least two different labels are in need in training data set while in my case I only have one label - all items in S would have the same label. Any ideas on how to find items in P that are similar to S??


title: Statistical notation in plain English

text: I have a very limited understanding of statistics, and am having trouble translating an answer to this question from quant.SE, "Separating the wheat from the chaff", the answer being: > Larry Harris has a chapter on performance evaluation in Trading and > Exchanges. He states that over a long period of time, a skilled asset > manager will consistently have excess returns whereas a lucky one will be > expected to have random and unpredictable returns. Thus, we start with the > portfolio's market-adjusted return standard deviation: > > \begin{equation} \sigma_{adj} = \sqrt{\sigma^2_{port} + \sigma^2_{mk} - > 2\rho\sigma_{port}\sigma_{mk}} \end{equation} > > where $\rho$ is the correlation between the market and portfolio returns. > > For a sample size $n$ (generally number of years), the average excess > returns, and the adjusted standard deviation from above, we have a > t-statistic: > > \begin{equation} t = \frac{\overline{R_{port}} - > \overline{R_{mk}}}{\frac{\sigma_{adj}}{\sqrt{n}}} \end{equation} > > Now we can simply determine the probability that the manager's excess > returns were luck by plugging this t-statistic into the t-distribution's PDF > with degrees-of-freedom $n - 1$. The lower the probability, the more we can > believe the manager's excess returns were from skill. At this point, I now understand the top formula, though this was based on the finding the book's the section on topic referenced in the answer above on Google Books. What I don't get now is this formula: > \begin{equation} t = \frac{\overline{R_{port}} - > \overline{R_{mk}}}{\frac{\sigma_{adj}}{\sqrt{n}}} \end{equation} I understand it's somehow related to the t-distribution's "Probability Density Function" but I have no idea what the R or "overline" notation means, or for that matter, where the values are R-port and R-mk come from; "port" being portfolio, and "mk" being market. If it matters, since it took me a bit to figure this out "adj" just flags the sigma (standard deviation) as the "portfolio's market-adjusted return standard deviation." Any suggestions?


title: MCMC Geweke diagnostic

text: I'm running a Metropolis sampler (C++) and want to use the previous samples to estimate the convergence rate. One easy to implement diagnostic I found is the Geweke diagnostic, which computes the difference between the two sample means divided by itsestimated standard error. The standard error is estimated from the spectral density at zero. $$Z_n=\frac{\bar{\theta}_A-\bar{\theta}_B}{\sqrt{\frac{1}{n_A}\hat{S_{\theta}^A}(0)+\frac{1}{n_B}\hat{S_{\theta}^B}(0)}},$$ where $A$, $B$ are two windows within the Markov chain. I did some research on what are $\hat{S_{\theta}^A}(0)$ and $\hat{S_{\theta}^B}(0)$ but get into a mess of literature on energy spectral density and power spectral density. (http://en.wikipedia.org/wiki/Spectral_density) but I'm not an expert on these topics; I just need a quick answer: are these quantities the same as sample variance? If not, what's the formula for computing them? Another doubt on this Geweke diagnostic is how to pick $\theta$? The above literature said that it is some functional $\theta(X)$ and should imply an existence of a spectral density $\hat{S_{\theta}^A}(0)$, but for convenience I guess the simplest way is to use the identity function (use samples themselves). Is this correct? The R coda package has a description but neither does it specify how to compute the $S$ values.


title: Transforming t-distributed residuals to normal

text: I have an AR process with t-distributed noise at each step. Does anyone know a transformation i can apply to get approximately normally distributed residuals? Thanks!


title: Nonlinear PCA ( R package homals)

text: I'm using the NLPCA to reduce the dimensionality of nine variables (4 nominal / 3 ordinal /2 numeric) to obtain the object-scores to be used as dependent variable in a regression model. I'm using the package homals (http://www.jstatsoft.org/v31/i04/paper). The output is:               X <- homals(sintomi[,2:10], ndim = 1, active = TRUE, level = c("numerical", rep("ordinal",3),"numerical",rep("ordinal",4)))          X     Call: homals(data = sintomi[, 2:10], ndim = 1, level = c("numerical",          rep("ordinal", 3), "numerical", rep("ordinal", 4)), active = TRUE)          Loss: 0.0002077596           Eigenvalues:         D1      0.0208           Variable Loadings:                D1     V1  0.2096132     V2  0.2420691     V3  0.2171551     V4  0.1370880     V5  0.1508221     V6  0.2537438     V7  0.1959783     V8  0.1769447     V12 0.2154044               predict(X)          Classification rate:       Variable Cl. Rate %Cl. Rate          1       V1   0.5149     51.49      2       V2   0.7292     72.92      3       V3   0.7589     75.89      4       V4   0.7768     77.68      5       V5   0.6336     63.36      6       V6   0.9048     90.48      7       V7   0.8869     88.69      8       V8   0.8036     80.36      9      V12   0.8661     86.61      I'm having the following questions:   1. Is it best to consider Ndim = rank = 1 or Ndim = rank = max (rank) to reduce the dimensionality of data in one component?   2. Is there a command to automatically calculate the proportion of variance explained by the first component? Otherwise, how can I calculate it by hand?    3. Is it necessary to standardize numeric variables before perfoming "homals"? If anyone has any thoughts for this, responses would be greatly appreciated. Thanks.


title: Differences in kurtosis definition and their interpretation

text: I have recently realised that there are differences in the kurtosis values provided by SPSS and Stata. See http://www.ats.ucla.edu/stat/mult_pkg/faq/general/kurtosis.htm My understanding is that the interpretation of same would therefore be different. Any advice on how to deal with this?


title: Finding the distribution when the observations are dependent

text: How do we find information about the distribution of a variable in presence of dependency among our observations? This dependency is coming from measuring the variable on the same group of subjects multiple times. These are two examples of what I have in mind:   1. Suppose we have a website for customers to purchase goods. We are interested in the distribution of the amount of time the customers spend on this site per visit. A customer might come back to the site many times. These visits are considered distinct observations. A customer's second visit might be because of an earlier promotion that they saw on the site. Therefore, our observations are not i.i.d. (or a simple random sample).    2. In a clinical trial, where, say, blood glucose levels are measured for subjects in 3 minute intervals and we're interested in their "typical value" in the 20 minute interval after exercise. Again, the values are dependent for each subject.  My naive answer is to reduce all observations from the same subject to one value, and to look at the distribution of these values, but I'm certain there's a better way to do this. * * * Update: After reading the responses from @PeterFlom and @NickCox, I realized that my question was somewhat vague and I needed to add some more detail. Let's assume that I'd like to find the mean of this distribution. The weak law of large numbers requires i.i.d. variables (probably there is a more general version, say for exchangeable variables, but that's still not a valid condition in the situations I mentioned above.) How does one go about finding "average time a user spends on the website per visit"? I welcome any suggestions, but I'm most interested in answers that have some theoretical justification. * * * Update 2: This is the mathematical formulation of my question. Hopefully, it'll clarify the ambiguities in the wording: Suppose for $N$ subjects, we've made observations:   * For subject 1, $X_{1,1},\ldots,X_{1,i_1}$.   * For subject 2, $X_{2,1},\ldots,X_{2,i_2}$.   * ...   * For subject N, $X_{N,1},\ldots,X{N,i_N}$. (So, the number of observations for different subjects is not the same.) All the $X_{i,j}$ have the same distribution, but they are not necessarily independent. My question is, if my goal is to estimate $E(X_{i,j})$, what would be the best estimator? I don't believe the answer $$\frac{\sum_{i,j}X_{i,j}}{i_1+\cdots+i_N}$$ is correct, because the law of large numbers doesn't apply here. My answer is $$\frac{\sum_i\frac{\sum_j X_{i,j}}{i_j}}{N}$$ but as I mentioned earlier, I feel there must be something better.


title: Repeated Measures in R, SAS or SPSS: Whale Body Temperature Data

text: I have a question regarding a repeated measures experiment. I have 5 whales (randomly selected) that I recorded body temperature on. Data was collected via satellite, programmed to transmit records at 2-minute intervals, but data for each whale are not a complete time-series data set (i.e. unbalanced records/missing data points). A state-space model was creating using the whale's GPS location to determine what state each whale was in while tracked at sea. The states are the following: resident foraging, resident transiting, and post-resident transiting. I've partitioned the body temperature data into the 3 states. I have also divided up the body temperature data into 2 diel periods: daytime and nighttime. So I have two subplot factors: Diel (2 levels: day and night) and State (3 levels: res. forage, res. transiting, post-res. transiting). I have more than 500 body temperature values for each state/condition for all 5 whales, and I have more than 1000 body temperature values for daytime and nighttime periods for all whales. But again, I don't have the same number of data points for all whales, due to lack of satellite transmissions and the fact that there were differences among individual whales (e.g. some foraged for longer periods). Can anyone offer non-parametric statistical advice (data transformations were all unsuccessful) for setting up a model analyze the effect of daytime/nighttime and state (and the interaction between these factors) on whale body temperature in R, SAS, or SPSS? Would be happy to explain more about the data or experimental design. Thanks, Jim


title: neuralnet error

text: I am trying to train a `neuralnet` model by using `caret` package by this command:               model.neuralnet <- train(X, Y, method='neuralnet', trControl=myControl,          stepmax=700, rep=1, learningrate=0.05, algorithm = 'backprob',          preProcess=PP, tuneGrid = expand.grid(.layer1=4, .layer2=1, .layer3=1),          act.fct = 'tanh', linear.output = TRUE)      `X`, `Y` are correct inputs, `myControl` should be proper settings, because I build in same session successfully other models with them. After every second output on console it writes something like this (I have repeated CV but folds must be OK, because other models had no problems...):               model fit failed for Fold02.Rep1: layer1=4, layer2=1, layer3=1      In the end of the computations with above command I get:               Aggregating results     Selecting tuning parameters     Fitting layer1 = numeric(0), layer2 = numeric(0), layer3 = numeric(0) on full training set     Error in if (tuneValue$.layer1 == 0) stop("the first layer must have at least one hidden unit") :      argument is of length zero     In addition: There were 50 or more warnings (use warnings() to see the first 50)     Timing stopped at: 0.16 0 0.15       But I think I do not have set `layer1` to zero (`tuneValue$.layer1 == 0`). Please - were I am doing mistake? Milan


title: Probability of finding a particular sequence of base pairs

text: Thinking about probability always makes me realize how bad I am at counting... Consider a sequence of $n$ base letters $A,\; T, \; C, \text{ and } G$, each equally likely to appear. What is the probability that this sequence contains a particular sequence of base pairs of interest of length $r\leq n$? There are $4^n$ different (equally likely) sequences possible. Start with the sequence of interest at the beginning of the full sequence; $4^{n-r}$ sequences like this are possible. We can start our sequence of interest in $n+1 -r$ different locations. Hence, my answer is $(n+1-r)/4^r$. This probability is increasing in $n$, which makes sense to me. But this probability exceeds 1 when $n>4^r +r-1$. But that can't be. The probability should approach 1 in the limit (seems to me), but not exceed it. I assume that I'm double counting something. What am I missing? Thanks. (FYI, not homework, just a toy example in preparation for exams. A question posed by my molecular biologist friend.)


title: Comparing rates of events between an exposed and unexposed group

text: I'm looking at 2 papers which are addressing the same issue using a retrospective cohort study design but use different methodology to analyse data and I would like to know why the methods used differ? When and why would you use standardized rate ratios? When and why would you use Poission regression? Thanks in advance.


title: Identifying a number of significance from a histogram

text: I'm creating a program that identifies trending music artists that uses Twitter metrics. I have data in a histogram format that represents the frequency of twitter @mentions of an artist for the last 30 days. I need my program to recognize a "significant" change (number that is quite greater than the rest in the sample) in the frequency of @mentions. Here is a data sample: **histogram: [ 43, 17, 70, 137, 198, 113, 126, 96, 100, 130, 107, 112, 438, 215, 76, 30, 119, 78, 38, 27, 27, 35, 29, 38, 59, 35, 164, 113, 29, 111 ],** In this scenario, `438` is the number of significance. It's fairly obvious in this data set b/c it is 2x greater than any other value. I need to build an equation that recognizes this number of significance across various data sets. In my first attempt at writing this equation, I will calculate the average of the 30 values and compare the average to the highest value in the data set. After doing some sampling, I will come up with a percentage difference that justifies significance. This is a simple approach, but I'm worried its too subjective. Are there any other ways of doing this?


title: ABC model selection

text: It has been shown that ABC model choice using Bayes factors is not to be recommended due to the presence of an error coming from the use of summary statistics. The conclusion in this paper relies on the study of the behaviour of a popular method for approximating the Bayes factor (Algorithm 2). It is well known that Bayes factors is not the only way for conducting model choice. There are other features, such as predictive performance of a model, that might be of interest (e.g. scoring rules). **My question is** : is there a method, analogous to Algorithm 2, for approximating some scoring rule(s) or other quantities that can be used for conducting model choice in terms of the predictive performance in contexts with complicated likelihoods?


title: Is there a sequential version of probabilistic latent semantic analysis?

text: Does someone know if it exists some way to do online learning with pLSA? The model training is really time consuming, so it is not feasible to rebuild it after every changes in the data.


title: Suggestions for computing correlations between evaluators

text: I have a data frame consisting of evaluations made by students of other student fellows. An excerpt of the evaluations data frame looks like this (columns are evaluators, rows correspond to the evaluated student):                      s1          s2    s3       s4     1       7           8    NA        8     2      NA           7     8       NA     3       9           6     6        8     4       9           9     7        7     5       9           8    NA        7            ...      I want to know how evaluators are related to each other in terms of similarity of their non-NA assessment values. Would it be OK If I use `corrgram(evaluations)` or is it advisable to use some other method? Thanks!


title: spss covariates variables in regression model

text: in the regression model in the spss statistic software can you explain me what are the covariates variables? if i want express a model like                y = a0 + a1*x1 + ....... + an * xn       which are the covariates?


title: Likelihood Function for Censored Structural Model

text: I am estimating a structural model with the following data generating process: $$x=\begin{cases}\begin{array}{c}y-\theta\\\0\end{array} & \begin{array}{c}if\ y>\bar{y}(\lambda_1)\\\if\ y\leq\bar{y}(\lambda_1)\end{array}\end{cases}$$ Further, y is a random variable with distribution $F(\cdot|\lambda_2)$, and $\bar{y}(\cdot)$ is an increasing function. The parameters I want to estimate are $\theta$, $\lambda_1$, and $\lambda_2$. I thought that this was a standard Type 1 censored model, and so wrote the likelihood as follows: $$L=\prod f(x+\theta|\lambda_2)^{I(x>0))}*F(\bar{y}(\lambda_1)|\lambda_2)^{1-I(x>0)}$$ However, after estimating the model, I do not think this is correct. You can see that the parameter $\lambda_1$ only contributes to the likelihood when the data is censored (ie $x=0$). Hence, since $\bar{y}(\cdot)$ is an increasing function, the value of $\lambda_1$ that maximizes the likelihood is positive infinity, which is essentially what my estimate for $\lambda_1$ was. However, $\lambda_1$ is clearly identfied, since it controls the selection equation. I then thought that instead the likelihood should be this: $$L=\prod \left(\frac{f(x+\theta|\lambda_2)}{1-F(\bar{y}(\lambda_1)|\lambda_2)}\right) ^{I(x>0))}*F(\bar{y}(\lambda_1)|\lambda_2)^{1-I(x>0)}$$ The reasoning is pretty straight forward; when $x>0$, y is truncated from below at $\bar{y}(\lambda_1)$. However, this turns out also not to work, for the same reason as before: the likelihood is strictly increasing in the parameter $\lambda_1$. I am almost positive that $\lambda_1$ is identified; if I was simply estimating whether or not x=0, then the likelihood would be: $$L=\prod \left(1-F(\bar{y}(\lambda_1)|\lambda_2)\right)^{I(x>0))}*F(\bar{y}(\lambda_1)|\lambda_2)^{1-I(x>0)}$$ which would have a unique maximizer in $\lambda_1$. So what should the correct likelihood be? Thanks!


title: How to print table from packet in R?

text: In package MASS there is data Insurance packet. And I need to: Print table, which row will be Group and column will be Age, in the table should be the sum of Claims meaning. Can someone help me to do this?


title: Measuring variation of one data set with respect to another data set

text: Each of the points in the plot below is a mean of 1000 samples. The standard deviation of a point on the RED line is "small" in that it is smaller than the difference between the mean of the RED line and the mean of the GREEN line at that x value. What is a good clear way to state that the variation of the data in both of these data sets is small? Should I state the standard deviation is less than the half the difference between the set means, or is there a more formal metric to describe the same? ![Example of two data sets](http://i.stack.imgur.com/PCsy1.png) PS: This plot is intended to show that GREEN is greater than RED with high likelihood.


title: In hypothesis testing, does a p-value less than alpha always mean you reject NH?

text: I have learned that a p-value that is less than alpha means you reject the null hypothesis. But this seems to only make sense in the case of lower tail tests. What do you do for upper tail tests?


title: On properly setting up a control chart

text: Let me present to you the way I have being 'doing' control charts up until now. I have a series of $k$ samples of size $n$. I calculate the the sample means $\bar{x}_i$ and the grand mean $\bar{\bar{x}}$ --- the mean of the sample means (or the mean of all the $k\times n$ data points). I also calculate the sample ranges $w_i$ as well as the grand range $\bar{w}$. Now I use this data to produce what is an approximate 95% confidence interval with AKA inner control limits: $$\bar{\bar{x}}\pm A'_{0.025}\bar{w},$$ where $A'_{0.025}$ is read from a table and depends on the sample size $n$. Similarly, I produce a 99.8% confidence interval with AKA outer control limits: $$\bar{\bar{x}}\pm A_{0.001}'\bar{w},$$ where again $A'$ depends on $n$ and is read from a table. Now I set up my control chart with the grand mean and both pairs of inner and outer limits and I plot the sample means and depending on where the sample means lie we make decisions. The rules we have declared are 1\. one sample outside the outer control limits and we have a problem 2\. two consecutive samples outside the inner control limits and we have a problem Now I have read some answers to other questions and I want to say that I have no problem making assumptions about my data (strictly speaking this is a math class). However I have a number of questions and problems with the way I have being doing things.   1. _Do I understand out of (statistical) control? If not is there a difference between being out of statistical control and out of 'quality' control._ \--- my understanding is that a process is out of control when a change affects the statistics of the process. This change in turn could be something we can affect or not although I am not concerned: all I will say is that the process should be stopped and examined if it is out of control. This leads onto my next question.   2. _We use the possibly-out-of-control data to decide what is normal_ \--- --- my understanding is that **under the assumption** that the population mean is equal to the grand mean we should see our sample means within the control limits with the relevant probability. However, this process could be 'out of control' as I understand it. This could mean two things (a) the population mean is changing or (b) the population mean has changed. In both scenarios I see the out of control data feeding back and essentially saying what is normal both in scenario (a) where the range $\bar{w}$ is potentially big or worse (b) where the grand mean is not estimating the population mean. I would prefer that the the control chart parameters are based on historical or aspirational data.   3. _Is it not backwards?_ \--- surely what should happen is that there is _a priori_ a control chart and a worker does the sample, reacts in 'real time' to adverse findings and then stops the process --- or is this more quality control/process capability.   4. _In the field?_ \--- how are control charts actually done in the field (at a low level, neglecting trend, stratification, etc). The people learning about this are training to be civil engineers --- although I also do it with Industrial Measurement Control classes (who I am sure see more of it in their other classes beyond maths).   5. _In other learning centres?_ \--- is this still taught? What approach is typically taken? Again I am not going to be going into trends, cycles, etc.   6. _Should I just teach quality control --- or does nobody do this?_ \--- would it not be a control chart if I worked off historical or aspirational data. Thank you very much for your help on this, I appreciate it. Not happy teaching it just by plug and chug --- it has to make sense the way we are 'doing' it and at the moment it doesn't.


title: Why are cumulative residuals from regression on stock and index returns mean reverting

text: In the paper > M. Avellaneda and J. H. Lee, Statistical arbitrage in the U.S. equities > market, July 2008, in the Appendix on page 44, I have some questions. First he runs the regression of stock-return ($R_n^S$) with index/ETF-return ($R_n^I$). $R_n^S = \beta_0 + \beta R_n^I + \epsilon_n, ~~~~~n=1,2,...,60$ Then he defines an auxiliary process $X_n$ as sum of the residuals from regression of $R_n$, and estimates it to be a mean reverting process. $X_n = \sum_{j=1}^{n} \epsilon_j ~~~~~n=1,2,...,60$ I have two related questions? What is the significance of taking cumulative sum of residuals as opposed to taking just residuals as a mean reverting process? I have a basic intuition but lack a good understanding. Second, he says that the regression on stock-returns "forces" the residuals to have mean zero. Why is this? How does this imply that the sum of all residuals, $X_{60}=0$?


title: How can I model random graphs (trees) with deterministic constraints?

text: I have been trying to represent a tree structure as a random graph, but so far my research has led to variations of random graph models and I can't see studies which may help me leverage some of the benefits of random graph models in a more specialized family of graphs. My data is in the form of tree structures, and I want a probability distribution representation of this tree structure so that I can assign probabilities to various instances of trees. There are rules that govern the structure of the tree: every node other than the root is dependent on its parent to simply exist. If the parent is not there, the child can't exist. If the child exists, its attributes have their own distributions, i.e normal, beta etc.. This almost sounds like a DAG with conditional probability tables. But I can't see how I could include the existence/non-existence of a child node in a CPT. My goal is to get a distribution for these tree structures and assign probabilities to them, but I could not find which modelling approach could represent them. Where should I be looking?


title: Comparing two means with same sample size

text: I am trying to do some web statistics (AdWords). I have a cost, a revenue, and a number of clicks. How do I determine if the difference between cost/click and revenue/click is statistically significant. I do not have access to the variance nor the standard deviation, of neither the cost nor the revenue.


title: Population stability index - division by zero

text: Population stability index quantifies the change of a distribution of a variable by comparing data samples in two time periods. It is very commonly used to measure shifts in scores. It is calculated as follows:   1) The sample from base period is discretized. Usually it is partitioned into deciles   2) The sample from target period is discretized using the same intervals as in first step   $PSI = \sum_{i} (A_{i} - B_{i}) \cdot ln(\frac{A_{i}}{B_{i}})$   Where:   $A_{i}$ - share of i-th bin in base period. $B_{i}$ - share of i-th bin in target period. **Question** : What should be done when one of the bins from target sample is empty?


title: How can I demonstrate non-linearity without categorising a predictor?

text: I don't know what is the appropriate term for my question. The scenario is described as following. In the analysis there one dependent variable Y and two independent variable X1 and X2. All three variables are continuous. I converted X1 into a categorical variable which has three levels A, B, and C. It was found that Y and X2 are positively correlated in group A and B, but negatively correlated in group C. I was told that convert a continuous variable into categorical was generally a bad idea and I understand this. My question is, how can I demonstrate the above pattern without break X1 into categories? I was suggested to go with multiple regression, but I still don't know how to demonstrate this kind of relationship in the three variables with multiple regression. Thanks.


title: How would you test for high order autocorrolation

text: The question is the following; ![enter image description here](http://i.stack.imgur.com/zuJye.jpg) And I determine the intervals ![enter image description here](http://i.stack.imgur.com/Kpkcs.jpg) If I Will do breusch-godfred test, I guess accourding to these interval, I can find high order correlation. But I dont know how to do this test. I dont understand it and dont understand the question properly. I have just started econometrics, so I am not good. Please help me solving the question. ThanK you:)


title: Test for statistical difference between two curves on same plot?

text: I've got data plotted for A/Ci curves (plant physiology, A=photosynthesis rate, Ci=internal CO2 concentration--essentially it's a dose-response relationship). Generally, when this type of data is presented, people just add SE error bars (like the provided image) and if they don't overlap call it significant, or they run an ANOVA on individual pairs of points (such as the two points, black and white, that sit above Ci=120). ![A/Ci curve](http://ejournal.sinica.edu.tw/bbas/content/2004/3/1747-Fi0.JPG) However, I would like to find a single test for the overall dose-response relationship. So, essentially I can assert that the dose-response relationship is significantly different between the two groups (in my case, the groups are all for the same species--it's a fertilizer experiment and I want to be able to claim that the addition of fertilizer changes the kinetics of the A/Ci relationship). I would also like to be able to test more than two curves on the same plot in the same manner (but I only have 3 curves, so I could do permutations if necessary). Also, one other issue (not represented in this image) is that there is also an X-axis error term, so the error bars should actually be crosses, not just error on the y axis--which makes using an ANOVA approach difficult, since you have to test points in both dimensions. Any suggestions for an appropriate test? I work mostly in R, so any link to existing scripts would be ideal--but I also have access to Graphpad Prism and JMP if anyone knows how to run this type of test in either of them.


title: Adding individual variance for each observation in DLM-model

text: I have a problem with my analysis made with R-package dlm. I have an averaged data vector y, where every yi has variance s2i. On the other words, my data is monthly averages and the standard errors for the means. What I want to do is take the monthly variance into account in the analysis. In my current model there has level and trend component and a seasonal component. The model is given by: dlmModPoly ( order=2, dV=(exp(par[1])+kk.s2), dW=c(0,exp(par[2])) ) \+ dlmModSeas ( 12, dW=c(exp(par[3]), rep(0, 12 - 2) ) ) So, how can I add the individual variances to this? Best regards Harry


title: Using GPML in Matlab for MultiClass Classification

text: I am using Rasmussen's GPML code in Matlab R2011a_student. I have training data (2560x29707) w/ labels (6 classes), and test data (640x29707). To prep the data I have 1. converted from sparse to full, 2. binarized the classes (ie. All classes that = 1 are 1, everything else is -1). I planned on running this 6 times to accommodate all the classes. I ran the following code (taken right from the documentation, but replaced the x,y,t values with my data): meanfunc = @meanConst; hyp.mean = 0; covfunc = @covSEard; ell = 1.0; sf = 1.0; hyp.cov = log([ell ell sf]); likfunc = @likErf; hyp = minimize(hyp, @gp, -40, @infEP, meanfunc, covfunc, likfunc, x, y); [a b c d lp] = gp(hyp, @infEP, meanfunc, covfunc, likfunc, x, y, t, ones(n, 1)); * * * I get the following error, and I'm not sure what it means. Any help would be greatly appreciated. ??? Error using ==> gp at 76 Number of cov function hyperparameters disagree with cov function Error in ==> minimize at 75 [f0 df0] = feval(f, X, varargin{:}); % get function value and gradient Error in ==> ds1pleasework at 7 hyp = minimize(hyp, @gp, -40, @infEP, meanfunc, covfunc, likfunc, full_TrainSet_feature, L_train); Thanks so much.


title: Predicted probabilities of square term in Stata

text: I am estimating a multivariate logit model. I have included a square term on my main independent variable, IDV and IDV:IDV. And of course a set of control variables. The results indicate an inverted U shaped relationship; IDV is positive, while IDV:IDV is negative. I have graphed the probabilities using `-postgr3-`:               postgr3 IDV, asis(IDV IDV*IDV)      However, I am unsure about how to get the probabilities of different values of IDV. To be clear, I am looking for something similar to the information `-clarify-` gives, i.e. the probabilities for different percentiles of my IDV, while taking into account the squared term. Any ideas? It would be very much appreciated. Let me know if I my question is not clear. I recently started using Stata.


title: Explained variance of non-negative matrix factorisation

text: How can I calculate R2 or another measure of explained variance for non- negative matrix factorisation. Currently I calculate the total sum of squares (TSS) as the sum of the squared distance of every data point to the mean and residual sum of squares (RSS) as the squared sum of the distance of every residuum to the origin. Then R2=1-RSS/TSS. So basically I did the same as one would for a regression. Now sometimes it can happen that one component/factor yields a R2 of .15 whereas two (with the same data) results in an R2 of .90. Can this be correct? After checking, all data points are roughly the same "size" (distance to origin) and the first component is similar to the mean. Can it be that TSS is wrongly calculated? Thanks. * * * edit: I think what I am really looking for is a method to judge how many components/factors are needed to explain my data. For example a method can be to want to have an R2 of >.9. Is there a better way than to use R2, especially regarding the fact that one factor (under special circumstances might explain very little variance, whereas 2 already describe a lot. * * * edit2: Would it be meaningful to describe the relate the root mean square residual to Dr=sqrt(norm(A-W*H,'fro')/(N*M)) to the root mean square of the original matrix Dtot=sqrt(norm(A,'fro')/(N*M)). A is the original matrix of size N-by-M and, W and H the factors as determined by the NNMF. 'fro' is the matlab code for the Frobenius norm. Thus, I would use 1-(Dr/Dtot) and a cutoff value to choose the number of 'factors' that describe enough of my data.


title: Dynamic calculation of number of samples required to estimate the mean

text: I am trying to estimate the mean of a more-or-less Gaussian distribution via sampling. I have no prior knowledge about its mean or its variance. Each sample is expensive to obtain. How do I dynamically decide how many samples I need to get a certain level of confidence/accuracy? Alternatively, how do I know when I can stop taking samples? All the answers to questions like this that I can find seem to presume some knowledge of the variance, but I need to discover that along the way as well. Other are geared towards taking polls, and it's not clear to me (beginner that I am) how that generalizes -- my mean isn't w/in [0,1], etc. I think this is probably a simple question with a well known answer, but my Google-fu is failing me. Even just telling me what to search for would be helpful.


title: Generating visually appealing density heat maps in R

text: While I know that there are a series of functions for generating heat maps in R, the problem is that I'm unable to produce visually appealing maps. For example, the images below are good examples of heat maps I want to avoid. The first clearly lacks detail, while the other one (based on the same points) is too detailed to be useful. Both plots have been generated by the density() function in the _spatstat_ R package. How can I get more "flow" into my plots? What I'm aiming for is more of the look the results of the commercial SpatialKey (screenshot) software is able to produce. Any hints, algorithms, packages or lines of code that could take me in this direction? ![enter image description here](http://i.stack.imgur.com/Wbhei.png)