title: Is it posible to perform the inverse of multidimensional scaling analysis

text: We have lot of 3D data and we reduced it to 2D for performing fuzzy clustering and obtaining prototypes. We used some matlab functions that were very well documented. Now we would like to see to which initial 3D points the first prototype found corresponds, so I wonder if we can perform the reverse of MDS, is that possible?


title: Boundary/threshold test for regression-type scatter plot

text: I am looking for a way to test weather a boundary threshold exists in a physiological response--a sample of the data is plotted below. My hypothesis is that the X-variable imposes a physiological constraint on Y-values, thus producing a boundary 'ceiling' for maximum Y-values that varies inversely with X-values (indicated by the red line on figure). I assume any Y-values below the boundary are limited by some other factor not included in this model. Essentially, my goal is to determine if the boundary exists and if so to derive a confidence interval for the boundary line model--similar to a linear regression model, but describing the upper bound of the y-values, rather than the center of mass. I'm sure something like this exists, but I haven't come across it before. Also, I would appreciate any suggestions on a better title or tags for this post--I assume there are more accurate terms for what I'm describing that would help folks find this post. ![threshold](http://i.stack.imgur.com/1BQBk.jpg)


title: The a priori power of test for a simple effect within a 2x2 ANOVA

text: I need to know about the power of a simple effect that follows from a 2x2 ANOVA interaction. Let's call the factors in this design GROUP(Patient, Control) and MANIPULATION(Expectation, Random), giving me 4 cells. Our very specific hypothesis is the following: controls will adapt their behavior according to the experimental manipulation, but patients won't. So, formally, the dependent variable changes only in the Control group when the manipulation is 'Expectation'. Thus, the means of three cells do not differ significantly, whereas the fourth mean differs significantly from the rest. If I were testing data, I would simply look for the _interaction_ and then the corresponding _simple effect_. However I need to a priori determine the sample size needed to test this hypothesis, and thus need to know the power of this test. This is where I need help. So far, I've done the power analysis for a GROUP X MANIPULATION interaction. But here's the problem: _I do not know whether this power analysis also covers the subsequent simple effects tests_. I'd like to ask you the following:   * Does a power analysis for a 2x2 ANOVA interaction also 'apply' to the subsequent testing of the simple effect?   * If not, is there a way how I can formalize my hypothesis 'in a single test' whose power can be determined?   * If not, how else could I estimate power and sample size for this hypothesis?


title: What sample statistics is increasing (or non-decreasing) as data grows

text: First please forgive me of using this non-mathematical terminology in the title. Assume we have two series of real numbers $X_1, X_2, \ldots, X_M$ and $Y_1, Y_2, \ldots, Y_N$ where it is possible for $X_i=X_j$ and $i \neq j$. Now if we combine series $X$ and $Y$, the resulting set of real numbers $X_1, X_2, \ldots, X_M, Y_1, Y_2, \ldots, Y_N$ is denoted as $X+Y$. Finally let us assume a statistics $T$ can be defined for $X$, $Y$, and $X+Y$ similarly. My question is: what statistics is guaranteed to have the property: $T(X+Y) >= T(X)$? For example, suppose we define $T(X)=max(X_1, X_2, \ldots, X_M)$. It follows that $T(X+Y)=max(X_1, X_2, \ldots, X_M, Y_1, Y_2, \ldots, Y_N ) \geq max(X_1, X_2, \ldots, X_M) = T(X)$. So the maximum statistics satisfy this property. Another example would be $T(X)=card(X)$. In both cases, as data grows (from $X$ to $X+Y$), the statistics is monotonically non-decreasing. On the other hand, the median statistics does not necessarily satisfy the property. My question is, in general, what types of statistics would hold the property of $T(X+Y) >= T(X)$? I am particular in range statistics.


title: How to identify if relationship between two categorical variables is significant or not?

text: I have one dependent categorical variable and other many explanatory categorical variables. The number of nominal values for any of these categorical variables may fall in range of 5-100. My objective is to explore the relationship between dependent categorical variable and any one of explanatory variable and to find out if there is any significant relationship. The tools I can think of are like contingency tables with tests like chi- square etc. Are there any other specific suggestions? Datasets that I am going to use are large with at least 50,000 datapoints (rows) and I am using spotfire (S+) for my analysis.


title: How to build a linear mixed-effects model in R?

text: I asked a question but it was a bit long and confusing so I will attempt to keep this shorter and simple (original post Mixed effects model equations) I am basically looking to see which factors (e.g. A, B, C such as habitat type, site, disturbance rate) most affect the response value (y). I am finally understanding the concept of models and how to compare AIC values to see which combination of factors best explain y (have greater effect on), but I am new to R so wonder if my basic coding is correct. Firstly I have all my data in a spreadsheet saved as a .csv file, so I read this file into R. I also open the package lme4. Then I was thinking of using the following code (although letters would be the headings fo each set of values)               m1<-lmer(y ~ A + (1|E), REML=FALSE)     m2<-lmer(y ~ B + (1|E), REML=FALSE)     m3<-lmer(y ~ A + B + (1|E), REML=FALSE)     m4<-lmer(y ~ A + C + (1|E), REML=FALSE)     m5<-lmer(y ~ B + C + (1|E), REML=FALSE)     m6<-lmer(y ~ A + B + C + (1|E), REML=FALSE)     m4<-lmer(y ~ A + D + (1|E), REML=FALSE)     m5<-lmer(y ~ B + D + (1|E), REML=FALSE)     m6<-lmer(y ~ A + B + D + (1|E), REML=FALSE)     m7<-lmer(y ~ C + (1|E), REML=FALSE)     m8<-lmer(y ~ D + (1|E), REML=FALSE)      My basic questions are:   1. Let us say that C and D are similar factors, and derived from the same data and I do not think it useful to see if they have additive effects, is it okay not to put them in the same model, i.e. mix and match as see fit, or should all combinations be incorporated?   2. I chose ML because not all sets of values are mixed and matched. Is this correct or should I use REML?   3. Some sets of values are non-continuous, e.g. brood number or habitat type (either 1 or 2). Should I be letting R know this. Someone mentioned coding each factor but I have no clue how to do this, or is the csv file enough and then perhaps code to let R know about these particular ones?   4. Lastly, is there also a way to see p-values to see if the factor has a significant effect rather than just being the best fit to explain y? I hope that this is more simple than my first post and easier to answer all at once. I really wanted to run this model today, but don't want to have to redo it all finding that I made a simple mistake.


title: Question related to binomial distribution

text: Let $X \sim Binomial (n, p)$ with both $n$ and $p$ known. Suppose for some non-increasing function $G:[0,1] \rightarrow [0,1]$, and some fixed $c_0 \in [0,1]$, we have that \begin{align} \label{ineq} P(G(X/n) \leq c_0) \leq \alpha \end{align} where $\alpha$ can be chosen as 0.1. Let $\mathcal{G}$ be set that contains all the functions satisfy the above inequality. I want to find the smallest such function $G_0$ such that $G_0 \in \mathcal{G}$ and $G_0(x) \leq G(x), \quad \forall x \in [0,1]$, where $G$ is any other function in $\mathcal{G}$. Can anyone give some hint on looking for such a function. Thank you very much. Hannah


title: Conditional probability question

text: somebody please prove or disprove $P(x|y)=P(x|z)/P(y|z)$ I have tried to plug in numbers with example in Wikipedia, for example, $A$ is outcome number of dice 1 $B$ is outcome number of dice 2 $X$: event of $A=2$ $Y$: event of $A+B \leq 5$ $Z$: event of $B=2$ $P(x|y) = 3/10$ , $P(x|z) = 1/6$ , and $P(y|z)= 3/6$ Question: Then the above equation does not Hold, am I missing anything?


title: Should statsmodels's GLM produce the same results as R's lm?

text: Should Python's               statsmodels.api.GLM(train_y, train_X, family=sm.families.Binomial()).fit().predict(test_X)      always produce the same results as R's               predict(glm(y ~ ., data=train_X, family=binomial), newdata=test)      where `train_y` is a pandas `DataFrame` containing the `y` column in the corresponding R `data.frame`, `train`; and where `test_X` and `train_X` are dataframes containing the remaining columns from the `test` and `train` dataframes respectively? If not, are there parameters that I can supply to `statsmodels`'s `GLM` to make it produce the same results as R's `glm`?


title: Naive Bayes fails with a perfect predictor

text: Let's say I have a variable that perfectly predicts one of the classes in my dataset:               set.seed(668130)     dat <- iris     dat$X <- sample(1:3, nrow(iris), replace=TRUE)         dat$X  <- ifelse(dat$Species=='setosa', 1, dat$X)     > table(dat$X, dat$Species)              setosa versicolor virginica       1     50         12        15       2      0         18        15       3      0         20        20      Why does the NaiveBayes algorithm fail on this dataset?               library(klaR)     > NaiveBayes(Species ~ ., dat)     Error in NaiveBayes.default(X, Y, ...) :        Zero variances for at least one class in variables: X      It seems to me that it would be reasonable to output a classification of 'setosa' 100% of the time, if X=1. Other algorithms (such as randomForest) do this:               library(randomForest)     > randomForest(Species ~ ., dat)          Call:      randomForest(formula = Species ~ ., data = dat)                     Type of random forest: classification                          Number of trees: 500     No. of variables tried at each split: 2                  OOB estimate of  error rate: 4.67%     Confusion matrix:                setosa versicolor virginica class.error     setosa         50          0         0        0.00     versicolor      0         47         3        0.06     virginica       0          4        46        0.08      Is the NaiveBayes algorithm mathematically undefined in this case? I know the specific dataset is a little contrived, but the problem pops up occasionally when I am cross-validating NaiveBayes models.


title: Where to go for introduction to introduction to Machine Learning?

text: I'm watching a Machine Learning course offered by Alex Smola at http://www.youtube.com/watch?v=OmJbGrQ93CA 14 minutes in and im confused. At that part he starts to explain Linear Regression. I've watched other introduction to Machine Learning courses and also feel lost. I have difficulty understanding the equations even when their explained. I feel I've missed something. I understand there is no silver bullet and im not looking for one but for someone who is starting in this area and feels that the introduction is too difficult are there recommend resources I should be reading/watching to get up to speed in this area ? I find other algorithms such as k-means and hierarchical clustering much easier to understand but maybe im not understanding at good enough level if im finding this course so difficult.


title: Generating correlated distributions with a certain mean and standard deviation?

text: Given a distribution A with a mean of $\mu_1$ and standard deviation of $\sigma_1$, how can I generate:   * Distribution B with a mean of $\mu_2$ and standard deviation of $\sigma_2$ and a correlation of $X_1$ with distribution A   * Distribution C with a mean of $\mu_3$ and standard deviation of $\sigma_3$ and a correlation of $X_2$ with distribution B and $X_3$ with distribution A Can someone please tell me if this even makes sense? My naive approach was the following:   1. Generate A with the given parameters   2. Generate B with the given parameters and then see if the generated values have the specified correlation with A. If not, regenerate B until this correlation is achieved.   3. Generate C using the approach in Step 2. However, I am not quite sure if this approach will terminate. Is there a better way to achieve this? I'd love to see an example in R.


title: Why does independence test use the chi-squared distribution?

text: The $\chi^2$ goodness-of-fit test uses the following _statistic_ : $$ \chi_0^2=\sum_{i=1}^n\frac{(O_i-E_i)^2}{E_i} $$ In the test, granting that the conditions are met, one uses the $\chi^2$- _distribution_ to calculate the p-value that given the $H_0$ is true one would observe such a value in a representative sample of the same size. However, in order for a statistic $\chi_0^2$ to follow a $\chi^2$-distribution (with $n-1$ degrees of freedom), it must be true that: $$ \sum_{i=1}^n\frac{(O_i-E_i)^2}{E_i}=\sum_{i=1}^{n-1}Z_i^2 $$ for independent, standard normal $Z_i$ (Wikipedia). The conditions for the test are as follows (again, from Wikipedia):   1. Sample representative of population   2. Large sample size   3. Expected cell count is sufficiently large   4. Independence between each category From conditions (1,2) it is clear that we satisfy conditions for inference from the sample to the population. (3) seems to be an assumption required because the discrete count $E_i$, which is in the denominator, does not result in a near-continuous distribution for each $Z_i$ and if it is not large enough there is an error that can be corrected with Yates' correction - this seems to be from the fact that a discrete distribution is basically a "floored" continuous one, so the shift by $1/2$ for each one corrects this. The necessity of (4) seems to come in handy later, but I cannot see how. At first, I thought that $Z_i=\frac{O_i-E_i}{\sqrt{E_i}}$ is necessary for the statistic to match the distribution. This lead me to the questionable assumption that $O_i-E_i\sim \mathcal{N}(0, \sqrt{E_i})$, which was indeed wrong. In fact, it is clear from the reduction of dimension for two sides of the equality from $n$ to $n-1$ that this cannot be the case. It has become apparent, thanks to whuber's explanations, that $Z_i$ need not equal each $\frac{O_i-E_i}{\sqrt{E_i}}$ term because $\chi_0^2=\sum_{i=1}^{n-1}Z_i^2$ (note the reduction in the number of summed variables) for standard normal random variables $Z_i$ which are _functionally_ independent. **My question** , then, is how can $\chi_0^2$ follow the $\chi^2$ distribution? What kinds of combinations of each of the $\frac{(O_i- E_i)^2}{E_i}$ terms result in squared standard normals $Z_i^2$? This requires the use of the CLT, apparently (and that makes sense), but how? **In other words** , what is each $Z_i$ equal to (or approximately equal to)?


title: Accuracy and area under ROC curve (AUC)

text: If we group examples with and without class labels using clustering techniques by treating the class as an ordinary nominal attribute, the resulting clusters can then be used for classifying test instances by assigning the most frequent class in each cluster to test instances that fall into that cluster. Why is it this method is less accurate than methods that are specifically developed for classification?


title: How to construct an interaction plot

text: I am debating how to construct an interaction plot with my supervisor. We have a dataset comprising 8 independent variables. We are trying to analyse the effect of 2 of the 8 independent variables on the dependent variable. My supervisor is suggesting that to draw an interaction plot, we first fit a full model using all of the 8 independent variables we have, and for every possible combination of the levels of the 2 independent variables that we are particularly interested in (let's call them `var1` and `var2`), calculate their predicted value based on the full model that we constructed earlier. However when applying this method I was running into a problem because in order for my statistical software to make a prediction, I had to assign values for the 8-2=6 variables that are left in the dataset, which are undetermined. So I suggested to my supervisor that instead of relying on the full model for calculating predicted values, I fit a model like the one below: `y = var1 + var2 + var1*var2` (i.e. instead of `y = var1 + var2+ var3 + var4 + var5 + var6 + var7 + var8 + var1*var2`) My supervisor, however, disagrees with my view and is telling me to go on by using the mean values of `var3, var4, var5, var6, var7, var8`, which I can calculate from our original dataset, to come up with predictions. Is there something wrong with my method of analysing the effect of interaction? I prefer my method because the interaction plot looks much better with my method. However, if my method is theoretically wrong then I guess I have to stick with what is said by my supervisor.


title: Differences between regression using base R and using rms package

text: OLS models can be created using the base R lm function. They can also be created using the ols function in the rms package. Unfortunately, I cannot find anywhere online (including in the ols help page) a summary of the differences between ols and lm, and what reasons there might be in practice to prefer one or the other in particular circumstances. (I assume that there are substantial differences, otherwise why should the authors of the rms package have created the ols function.) A similar question arises with the competing mixed model functions in the lme4 and nlme packages, amongst (I suspect) many others. There does not seem to be anywhere on the CRAN site where comparisons of the features of packages with overlapping coverage can be found. Are such comparisons available elsewhere online?


title: Fitting proportions in a Bayesian model and controlling the error around the value

text: I am using a Bayesian state space model constructed in the BUGS language to fit my data. One of the outputs of the "state" part of the model is proportions. I use a beta distribution in the "observation" part to fit predicted values to observed ones with the following code.               for(i in 1:n) {      y[i] ~ dbeta(alpha[i], beta[i])      alpha[i] <-    mu[i]  * phi      beta[i]  <- (1-mu[i]) * phi     }      In that case, `phi` is controlling the spread around the mean value, but this "error" is variable when `mu` is changing. Its value follows a bell curve with highest error occurring at a proportion of 50%.   * Is there a way to fit proportions with a constant error?   * What do you think of the use of a truncated [0,1] normal distribution?


title: Good way to use adaptive learning rates in neural network

text: Adaptive learning rates means using different learning rate for different weight in neural network. Except for the emperical method which updates these learning rates based on consistency in gradient, what is the best method to use these adaptive learning rates? As I don't know much about Jacobian, Hessian matrices etc., can anyone please explain the method in simple way... Thanks in advance


title: Significance of Spatial Data

text: ### Background: It's common in many aspects of science to want to determine if the trends shown in spatial data are significant. In my field of familiarity, climate science, authors often show the trends in temperature data, for example, and overlay stippling to show trend significance at 5%. For example see Rotstayn et al. (2013) (open access) Fig. 1: ![Zonal-mean temperature trends with stippling denoting trends significant at 5%](http://i.stack.imgur.com/baCYv.jpg) Rotstayn et al. mention that > statistical significance of ensemble means is assessed using a two-sided > t-test [Sect. 2] For the example above, let's assume they have three data sets: model1 --> temperature climatology; model2 --> temperatures from an alternative model and model3 --> temperatures from another alternative model. (a) and (b) are then _model2 - model1_ and _model3 - model1_ respectively. In this case, the plot is latitude-pressure (height), but plots may be lat-lon, lon-pressure etc. Ultimately, you're working with a set of 2d matrices [x,y]. ### Question: How do you obtain significance at each discrete data point rather than obtaining one t/p statistic for the comparison of the entire data set, such that you can plot of field of 5% significance. I plot and process data in Python, so bonus points for discussion how this can be implemented in Python. Also feel free to discuss the merits of t-tests and/or suggest alternative tests that may be better suited to data that is not normally distributed. ### Note: _A similar question has been asked on Stackover flow here:Test for statistically significant difference between two arrays but I don't understand how the answers help me._


title: Formula to calculate beta matrix in multivariate analysis

text: I have to implement a multivariate analysis on $n$ random variables with a sample of $m$ data points. I would like to get a matrix with the $\beta$ (as in $n$ $\beta$ vectors put together). Is there an analytical formula for that or do I have to calculate each vector individually and put them back together ? Also is there a formula that would give me directly a vector of the estimation error for each variable?


title: Can I average samples in nested factor within-subject design if samples are from the same condition?

text: This question relates to one of my other questions: Split plot with nested factors ANOVA in R . I have subjects taking an experiment, where I have different `Scenarios` and each `Scenario` contains some number of `Trials`, where subjects can interact with each other between `Trials`, but the `Scenario` stays the same. Can I consider each `Scenario` as a different treatment condition and the trials as samples so I can compute the average on the `Trial`s of each subject for each `Scenario`?


title: use multiple Fisher test for periodic components

text: i know that there is fisher test for periodic component detection,on example is there http://www.mathworks.com/help/signal/ug/significance-testing-for-periodic- component.html but let us suppose that we have several periodic components,then what should be alternative form of Fisher method or even alternative method?let say we have 4 periodic component,and we have 4 peak in spectrum estimated by periodogram,how can i apply test to these 4 maximum value?another question if i apply fisher test to several maximum value of power spectral density,is it possible to detect out of them only periodic components using p value?please help me in this topic to understand well how can i detect periodic components in noise when we have not on, but several deterministic components with frequencies not necessity fourier frequencies?these frequencies may not be harmonically related to each other?


title: maximize mean F1 score in multilabel information retrieval problem

text: I have a multilabel text classification problem where each observation will have one or more labels associated to it. The metric I want to maximize is mean F1 score. Are there standard ways to approach the optimization of which labels and how many labels to predict for each observation to maximize mean F1 score over all labels? The approach I've taken so far is to run one-vs-all classifiers to predict each label. It's not hard to maximize each classifier's F1 score by finding an optimal cutoff. But how do you go about maximizing the mean F1 score over all labels? For example, let's say I have 3 labels with priors equal to respectively 0.4, 0.02, 0.7. For observation 1 I have: pr(label1)=0.34 pr(label2)=0.45 pr(label3)=0.87 If I sort by predicted probabilities and predict any label above priors, then I would choose first 3 and then 2. Already at this point I could ask myself if adding label 2 after label 3 is optimal from a mean F1 score perspective. I could also argue that label2's probability is much greater than its prior than label3 is of its prior. So, in this case I could predict label2 first and then ask whether entering label3 will still increase the mean F1 score. This is just a simple example, but I'm really trying to understand if there are accepted best ways to deal with this type of problems.


title: What machine learning tool is best suited for taking time series data as well as descriptive data and making a binomial classification

text: I have an interesting task of utilizing log data from computer servers in a server farm and predicting if a particular server is likely to fail in the next 24 hours. My data set will be comprised of both time series data as well as some descriptive data for each servers (see examples below) Time Series: I'll have data such as number of cycles per hour, average memory usage per hour, number of error thrown per hour. Descriptive data: Descriptive data will include things like: Computer model Processor Type Amount of Memory Time since last reboot Time in Service Given this data, I need to (1) massage the data and develop an appropriated combination of time series and descriptive data and (2) apply the appropriate predictive model to the data set to determine the probability of a server failing in the next 24 hours. I have done some research and have identified some potential predictive modeling approaches:   1. Logistic Regression   2. Neural Network   3. Support Vector Machine I'm sure there are others. My questions follow:   1. What data transformations if any would you suggest.   2. What type predictive model is best suited to my problem and why?  I have found a lot of information on various approaches, but have not come across something that says this model is best for these sorts of data sets and objectives. Thanks for your input.


title: Marginal and joint distributions of linear combinations of random vectors

text: Let $X_1,...,X_4$ be independent $N_p(μ,Σ)$ random vectors. Let $V_1,V_2$ be such that $$V_1=(1/4)X_1-(1/4)X_2+(1/4)X_3-(1/4)X_4 $$ $$V_2=(1/4)X_1+(1/4)X_2-(1/4)X_3-(1/4)X_4 $$ I need to find the marginal distributions of $V_1$ and $V_2$ and the joint density. Since they are linear combinations of random vectors I do not know the theory behind it to solve this. Any answers will be much appreciated. Thanks


title: Unsupervised Random Forest for Visual Codebook generation

text: I'm trying to apply the bag of visual words approach to make scene classification. I started to use k-means to generate my codebook, but rapidly discovered its limitations. From one codebook generation to another, I get fluctuations of 10% on the classification ratio (data dimension: 8, number of clusters: 50). In addition, k-means is slow (I'm working on Matlab). So I did some research and found that I can do unsupervised learning with random forests (paper: "Unsupervised Learning With Random Forest Predictors" by T. Shi and S. Horvath). It looked nice, so I decided to give it a try. What I did: (1) I created a synthetic data set from my original data set. To generate the coordinates of my synthetic samples, I take the corresponding coordinate of one of my original samples, chosen randomly (like explained here:http://www.stat.berkeley.edu/~breiman/RandomForests/cc_home.htm#unsup). I also tried the Addcl2 method like explained in the paper that I mentioned above. (2) The original data is labeled 1, the synthetic data is labeled 2. Then I train a random forest (I use Matlab's TreeBagger class) on this 2 classes. Then I compute the proximity matrix, but only for the original data. (3) I use this proximity matrix as a distance matrix to perform unsupervised learning with hierarchical clustering. Here is what I get for 2D data: ![Example for well separated data](http://i.stack.imgur.com/Hw0yV.png) ![Example for badly separated data](http://i.stack.imgur.com/CqbRG.png) It works approximately well for good separated data points, while the clustering result depends from one time to another, sometimes I get the same result as k-means. For badly separated data, the result is awfull. What am I doing wrong? I could need some advice ;-)


title: Establishing significance levels post-hoc

text: Proper frequentist procedure says that we should establish a sampling plan, statistical analysis and significance level before we start an experiment. We then report whether or not we reject H0 for the chosen level of significance. What if we do the same and don't choose the significance level beforehand. We simply look for the lowest significance level we can get away with, and report that, stating that we can reject the null hypothesis with the given level as significance. The lower this value is, the more weight we can give to our analysis. If we are upfront about following this procedure, is this honest? If not, where is the problem? I'm having trouble figuring out what the issue is, and why people don't report like this in general.


title: Probability distribution with standard deviation proportional to mean

text: Is there a probability distribution for which the standard deviation is proportional to the mean? More generally, I have been searching the web for a nice table comparing the properties of probability distributions (including mean and std). Any good reference?


title: Predicting proportions from time with a discontinuity

text: Imagine   * having a dependent variable $Y$ that is a proportion (i.e., the proportion of observations made at the given time point that satisfy a condition, where each time point involves 50 to 250 observations)   * $Y$ is measured at a series of time time points $X$, where $X = 1, 2, 3, ...$, typically to around 400.   * At initial time points, $Y$ typically equals zero or close to zero   * After an extended period of time $Y$ typically equals one or close to one   * At some point in between a transition occurs where values of $Y$ increase   * Throughout there is considerable time point to time point variability and given that Y is a proportion, the distribution of errors is not normal. Note also that the values of zero and one are common. Properties of the data vary across studies, such as:   * the initial value of $Y$   * the time point when the value of $Y$ starts to increase   * the duration of transition from values mainly around 0 to mainly around 1 **Questions**   * What would be a good modelling approach to such data?   * How could the onset of the transition from values close to zero to values close to one be detected, especially given the non-normal errors?


title: Kolmogorov–Smirnov "with perturbation"

text: Let $F$ be known continuous CDF of a continuous R.V. and $F_n$ represent the empirical CDF for sample of size $n$, hypothesized to be drawn from $F$. The Kolmogorov–Smirnov statistic is $D_n :=sup |F_n(x) - F(x)|$. In applying the Kolmogorov–Smirnov test, the key is that $\sqrt{n}D_n$ converges in distribution to the so-called Kolmogorov distribution, and this fact is then used to compute p-values. Now I wish to consider a slight variation. Supposed my samples are hypothesized to be drawn from $F$, but $x_i$ has been independently perturbed by an error $E_i$ (each $E_i$ has different distribution) of mean 0 before being recorded. If I construct the same empirical CDF using the perturbed values and attempt to apply the KS test, first of all, is it true that $\sqrt{n}D_n$ still converges in distribution to the Kolmogorov distribution (or some modification of it)? [Note: If $E_i$ are iid~$F'$, then I can just use the original KS test on convolution of $F$ with $F'$. But what happens otherwise?] [Maybe a proof of the original convergence theorem will also help to answer the new version with perturbation.]


title: Help creating a chart to show categorical data over time

text: The scenario: Over a course of 10 minutes 2 security guards patrol a floor and have their movement tracked. They can either be on the left or the right side of the floor, and within each side of the floor are 10 zones. So in total there are 20 zones (left zone 1, right zone 1, left zone 2, right zone 2, left zone 3 right zone 3 etc etc ) As they cross into a new zone the time is recorded. I am trying to create a chart that has time along the x axis, but I am unsure how to lay out the y axis ?? I had thought about putting the guards on the y axis and having a symbol represent each zone so I could just insert the symbol at the appropriate time as the guard moved into another zone (but that would need 20 different symbols, or just 10 symbols and use different colors for right and left zones) e.g                guard 1                            guard 2                             0   1  2  3  4  5  6  7  8  9  10                                time       I'm not so sure this is the best way to illustrate it. Ultimately this would be done on a large scale to see if there are any trends are evident. I have searched google images to find a similar structured chart but could find one. I have tried to use Excel to make a similar chart but failed. any advice as to how people would approach this ? have I got my axis muddled up perhaps ?


title: What are the assumptions of linear regression (both simple and multiple)

text: I know this is a basic question, but I get slightly different answers everywhere I look, a quick example is that I learnt at Uni that the errors of the dependent variable have to be normally distributed not the variable itself, but just recently I was brushing up on my stats and I googled 'assumptions of regression' and in one of the top links to be returned it stated that the data had to be normally distributed, I put a lot more confidence in my uni teachers than an internet link but the page claims to be written by a university lecturer. Also some places often state a good explanation of some assumptions not others. So I want to post it in an online community to get a consensus. So here are the exact questions: 1 What are ALL the assumptions of simple linear and multiple linear regression 2 How do you test each one (preferably in R) 3 Can you provide examples of how to simulate each assumption being violated and deomstrate how it has erroneously affected the results. Thanks


title: Strategy to deal with rare events logistic regression

text: I would like to study rare events in a finite population. Since I am unsure about which strategy is best suited, I would appreciate tips and references related to this matter, although I am well-aware it has been largely covered. I just don't really know where to begin. My problem is a political sciences one and I have a finite population comprising 515,843 records. They are associated to a binary dependent variable with 513,334 "0"s and 2,509 "1"s. I can coin my "1"s as rare events since they account for only 0.49% of the population. I have a set of around 10 independent variables I would like to build a model with to explain the presence of "1"s. Like many of us, I read King & Zeng's 2001 article about rare events correction. Their approach was to use a case- control design to reduce the number of "0"s, then apply correction to the intercept. However, this post says that King & Zeng's argument was not necessary if I already collected my data over the whole population, which is my case. Therefore, I have to use the classical logit model. Unfortunately for me, although I obtain good significant coefficients, my model is completely useless in terms of prediction (fails to predict 99.48% of my "1"s). After reading King & Zeng's article, I wanted to try a case-control design and selected only 10% of the "0"s with all the "1"s. With almost the same coefficients, the model was able to predict almost one third of the "1"s when applied to the full population. Of course, there are a lot of false-positive. I have thus three questions I would like to ask you: 1) **If King & Zeng's approach is prejudiciable when you have full knowledge of the population, why do they use a situation where they know the population in their article to prove their point?** 2) **If I have good and siginificant coefficients in a logit regression, but very poor predictive power, does that mean that the variation explained by these variable is meaningless?** 3) **What is the best approach to deal with rare events?** I read about King's relogit model, Firth's approach, the exact logit, etc. I must confess I am a lost among all these solutions. Thank you for your help!


title: Separating and identifying long and short term effects of statistical variables

text: A quick general question: In a practical setting, what's a good way to separate out, and then comment on, long term and short effects in a model? I had thought a good way to do so would be to include lagged versions, of the same variable, of varying distances, the longer distances indicating longer term effects, but I cannot think of anything else. Any suggestions greatly appreciated.


title: Graphing the usage of legal precedents

text: I would like to hear recommendations on what tools to use to graph the usage of precedents in a legal court in Brazil. The data is searchable at www.stj.jus.br/SCON/, and in the most recent collegiate decisions ('acórdão') database, a field ('Veja') is provided listing the cited cases. My objective would be to graph those relationships, trying to find often repeated cases and such. The dataset is relatively big, for courts anyways: about 30k collegiate decisions per year, 300k by 'single judges'. But it could be sliced to be smaller. From what I've seen, Gephi would be capable of doing this, is there a better alternative? I have about 8 months to complete the project, and have some help from a java programmer, and probably will have to learn a lot of this by myself.


title: Constructing Counterfactuals and Estimating Prevalence

text: I'm a social scientist working on a research project where I try to estimate the prevalence of lying in responding to a certain sensitive question. The way I estimate it is to rely on a quasi-experimental setting where a major event happen in the middle of a public opinion survey and compare the people's responses to certain questions before (control) and after (treatment) the event. Each individual in the survey sample is only asked once, but with the assumption that the event is exogenous and the progression of the survey is uncorrelated with personal characteristics (which appeared to be the case after some checking), the attitudes of the treatment and the control groups are comparable. My hypothesis basically suggests that there are two main types of individuals, those who lie about their attitude after the event, and those who don't. I have a measure that I believe should increase if one is lying, and decrease if not. Formally, it means that $y^{i,T}_{z}>y^{i,C}_{z}$ if $i$=liar, and $y^{i,T}_{z}<y^{i,C}_{z}$ if $i$=non-liar, where $y$ is the measure for lying, ${T,C}$ indicate treatment status, and $z$ indicates actual treatment assignment. I wish to estimate the prevalence of the behavior of lying in my sample, that is, the **_proportion_** of people in the treatment group who might have lied about their preferences. The problem is that only one of the state is observed given the nature of the survey, so I need to estimate the counter-factual using some methods. More specifically, I wish to know: $\frac{\mathbb{E}(\mathbb{I}(y^{i,T}_{z}>y^{i,C}_{z})|z=1)}{\mathbb{E}(\mathbb{I}(y^{i,T}_{z}<y^{i,C}_{z})|z=1)}$. But given the nature of the survey, I only know $y^{i,T}_{z=1}$ but not $y^{i,C}_{z=1}$, for $i\in {T,C}$, so I need to estimate the counter-factual in order to calculate this proportion. My question is: what would be the best way to do this? So far I have tried two things: (1) do a simple linear prediction from the regression equation on all the individual in the treatment group but change their treatment status from 1 to 0. (2) treat this as a missing-date problem and use various multiple imputation method to create the counter-factual. I personally prefer the second method since it takes into account the uncertainty in constructing the counter-factual. But I haven't not seen any published work employing multiple imputation in such a way, so I'm not entirely sure whether this approach is appropriate. I would appreciate if someone could tell me whether these methods are valid and/or whether there is any existing method that I can use for this task. Any suggestion would be greatly appreciated!


title: Simple approximation of tail surprisal of poisson distribution

text: I want to determine (in an algorithm) the approximate surprisal of getting an outcome "as extreme as $k$" from the $Poisson(\lambda)$ distribution. My original plan was to use $-log_2(1/2-|1/2-F_{Poisson(\lambda)}(k)|)$ but calculating $F_{Poisson(\lambda)}(k)$ requires the incomplete Gamma function which uses an iterative numerical algorithm. I am hoping to find a simple approximation (of the CDF, or even better, of the surprisal itself) that will allow me to speed things up, because I'm doing this for each pixel in a video. The fact that erf has such simple rational function approximations made me think $F_{Poisson}$ might have something similar (though we have a 2D domain here so perhaps that was wishful thinking.) As a last resort, I might create a 2D lookup table with a log-log-spaced grid and use bilinear interpolation (on the surprisal directly). ## Error The surprisal only needs to be accurate to within about 1 bit for all surprisals less than about 15 bits, translating to the requirement that $F_{Poisson(\lambda)}(k)$ is approximated to within about 50% of the distance from whichever of 0 or 1 is closer until the remaining tail probability is about $2^{-15}$. I'm just making up these specifications to give a rough idea. ## Previous related questions My original question on SO involved a very skew Binomial distribution (of unknown N), but it is so skew that I lose nothing by approximating that with the Poisson distribution, thereby getting rid of one domain variable. To be honest, I don't feel comfortable asking such an asker-specific question here, but the discussion here is what encouraged me to post it anyway. Any pointers or advice would be helpful.


title: Use repeated measures to predict binary outcome

text: I had an experiment using a 2 (between-subject variable, say Drug treatment) X 3 (within-subject variable, say Condition, which is counter-balanced, so time is not a concern here). Each condition has a score (continuous), so it can be considered repeated measures in this sense. The outcome is binary, say Succeeded (Y or N). I was interested in examining if Drug, Condition, and possibly DrugXCondition could predict Succeeded. I'm thinking to use logistic regression, but I'm not sure which procedure I should use in SAS, and what parameters I need to set. Would you give me some advice and/or sample code?


title: Are best linear unbiased predictions (BLUPs) a good indicator of the mean value for that random effect member?

text: I have data on prices of houses in different districts, and would like to determine how expensive different districts are when it comes to buying a house. However, houses vary with respect to amenities. For example, some have swimming pools and others don't. In order to remove any between-house variation in price per square meter which is due to the presence of a pool, I would like to compute the following model in `lme4`, where `has_pool` is 1 if the house has a swimming pool, and 0 otherwise:               lmer(log(rent_per_sq_meter) ~ has_pool + (1|district), data)      What I would like to know is, whether the BLUPs for districts can safely be interpreted as indicators of how expensive a district is, or whether there are any caveats I overlooked. I should probably point out that the number of houses per district is not balanced.


title: Not normalizing data before PCA gives better explained variance ratio

text: I normalized my dataset then ran 3 component PCA to get small explained variance ratios ([0.50, 0.1, 0.05]). When I didn't normalize but whitened my dataset then ran 3 component PCA, I got high explained variance ratios ([0.86, 0.06,0.01]). Since I want to retain as much data into 3 components, should I NOT normalize the data? From my understanding we should always normalize before PCA. By normalizing: setting mean to 0 and having unit variance.


title: Methods for evaluating partial autocorrelation for identification of ARIMA models

text: I am trying to programmatically identify an ARIMA model for a series of data and forecast values. Currently the problem i am facing is to find a way to evaluate partial autocorrelation. I have been looking for methods to calculate PACF for quite a long time now but in vain. Please provide some online resources which can help me in this matter. Thank you.


title: Two parameters zero at the same time

text: How can I test whether two parameters of a logistic regression model are both zero at the same time?


title: Two-way clustering in R

text: I'm wondering how to implement two-way clustering, as explained in Statistica documentation in R. Any help in this regard will be highly appreciated. Thanks


title: Unbalanced test data matter?

text: I have **balanced** training data: 300 positives and 300 negatives, but the test data is **unbalanced** : I have 15 positives and 60 negatives. Will the unbalanced test data impact classification accuracy? Actually, I tried to collect balanced test data (30 pos & 30 neg). The precision on positives are better. Is the reason the balance of test data? Thanks.


title: Finding extreme values in a normal distribution

text: I want to find extreme values (anything greater or less than three times standard deviation from the mean) after generating a set of random numbers using:               num = rnorm(1000)      My code looks like the following:               extreme = function(varname) {       for(i in varname) {         count = 0         m     = mean(varname)         sd    = 3*sd(varname)         if(i<(m-sd) || i>(m+sd)) {           count = count + 1         }       }       if(count>0) {         print(paste("There are ", count, " extreme values found.", sep = ""))       } else print("There are no extreme values.")     }      I'm always getting "There are no extreme values." I'm a beginner in R, so are there truly no extreme values in any randomly generated set of numbers following normal distribution?


title: Feature selection: all features vs a subset of them

text: I am doing a binary classification. The dataset has 3000 samples, and each sample has 10 features. But I find that the performance of using all 10 features is almost the same as that of using only the first three features. So does it mean that the last 7 features are not useful? Thanks in advance.


title: Intuitive explanation of contribution to sum of two normally distributed random variables

text: If I have two normally distributed independent random variables $X$ and $Y$ with means $\mu_X$ and $\mu_Y$ and standard deviations $\sigma_X$ and $\sigma_Y$ and I discover that $X+Y=c$, then (assuming I have not made any errors) the conditional distribution of $X$ and $Y$ given $c$ are also normally distributed with means $$\mu_{X|c} = \mu_X + (c - \mu_X - \mu_Y)\frac{ \sigma_X^2}{\sigma_X^2+\sigma_Y^2}$$ $$\mu_{Y|c} = \mu_Y + (c - \mu_X - \mu_Y)\frac{ \sigma_Y^2}{\sigma_X^2+\sigma_Y^2}$$ and standard deviation $$\sigma_{X|c} = \sigma_{Y|c} = \sqrt{ \frac{\sigma_X^2 \sigma_Y^2}{\sigma_X^2 + \sigma_Y^2}}.$$ It is no surprise that the conditional standard deviations are the same as, given $c$, if one goes up the other must come down by the same amount. It is interesting that the conditional standard deviation does not depend on $c$. What I cannot get my head round are the conditional means, where they take a share of the excess $(c - \mu_X - \mu_Y)$ proportional to the original variances, not to the original standard deviations. For example, if they have zero means, $\mu_X=\mu_Y=0$, and standard deviations $\sigma_X =3$ and $\sigma_Y=1$ then conditioned on $c=4$ we would have $E[X|c=4]=3.6$ and $E[Y|c=4]=0.4$, i.e. in the ratio $9:1$ even though I would have intuitively thought that the ratio $3:1$ would be more natural. **Can anyone give an intuitive explanation for this?** This was provoked by a Math.SE question


title: Usefulness of Z-normalization in Machine Learning

text: _Z_ -normalization means rescaling the feature $X$ by subtracting the average $\mu$ and dividing by its standard deviation $\sigma$, i.e., $(X-\mu)/\sigma$. What is the usefulness of normalizing data for machine learning?


title: Almon Model for Distributed Time Series Analysis

text: Can anybody give me a good reference on Almon model for time-series analysis. Also, what's the basic difference between Koyck model & Almon model. Regards,


title: Distribution of eigenvalues given one is known

text: I'm familiar with using insights from Random Matrix Theory to determine the number of principal components from the PCA of a covariance/correlation matrix to use to form factors. If the eigenvalue associated with the first PC is large, then it means that the remaining eigenvalues must be small (since the sum of the eigenvalues must equal the trace of the correlation matrix). When the first PC is large enough, it is thus possible that all of these eigenvalues are below the lower bounds on the Marcenko-Pastur distribution. This makes sense that they are low not because of random chance, but because the first eigenvalue is very large. However, that does not mean that they contain significant information. Rather, it would make sense to instead ask the question "given the first PC is some large number, what would the distribution of the remaining eigenvalues look like if random data were responsible for them?" Is there any research that addresses this issue? If it is possible to get the Marcenko-Pastur distribution conditional on knowing one or more eigenvalues, then it would be possible to proceed iteratively to determine whether the factors reflect significant information.


title: Show there's a disjoint partition over the sample space

text: Been working on this problem - any suggestions would be greatly appreciated! Let $e \in\varepsilon$ be a simple function over $(\Omega,A,P)$. Show that there is a disjoint partition $A_i, i = 1,\ldots,m,$ of $\Omega$ such that $e = \sum\limits_{i=1}\alpha_i1_{A_i}$.


title: Dealing with non-flat, unbounded data on a single observation

text: I'm currently trying to generate a number of (rudimentary) models to perform classifications (or get probability distributions through classifications). Let's say that I have a particular video from YouTube: UMVC3 Grand Final AGE ChrisG vs TA Marvisto - Next Level Battle Circuit #3 Tournament Going through version 3 of the YouTube API I'm able to get the following Freebase topic ids about the video:               "topicDetails": {         "topicIds": [             "/m/0h1gqk6",             "/m/0ffvx_"         ]     }        * /m/0h1gqk6 - Ultimate Marvel vs. Capcom 3 (correct indicator)   * /m/0ffvx_ - Battle Circuit (incorrect indicator, different topic) As you can see, `topicIds` is an array. There's no upper bound. Now when trying to create an observation for either evaluation into a statistical model, or generating a statistical model, I don't know how to include the `topicIds` in the observations. I had the following ideas: **`Topic1`, `Topic2`, `Topic3` ... `TopicN`** I've only ever seen four topics at most on a tag, so I could have say, ten variables on an observation. However, the variables `Topic1` ... `TopicN` have an implication on the positioning. It's possible that the `/m/0ffvx_` topic id would be moved into the third slot, affecting the processing of the model (since a particular topic id may be spread out among multiple columns). **Binary classification variables** I could have a separate variable for whether _every_ tag is observed. Something like:               Has_Freebase_TopicId_/m/0h1gqk6     Has_Freebase_TopicId_/m/0ffvx_      Where I would set it to true if that topic ID exists, false otherwise. This solves the problem above of not properly capturing the observation for that topic ID, but it poses the following problems:   * If a new topic ID is observed, then there's nothing I can do about it. I don't think it's an issue anyways, as a model can't classify variables it's never observed, I can just quantify if it has something observed or not and use that.   * Currently, I have distinct 1,571 topic IDs, and there's a potential for many, many more. The current model I'm using (Maximum Entropy Model) can handle this many variables, but this is for one unbounded variable in the data I'm observing. Ultimately, this will be used to generate a maximum entropy model which will take the topic IDs (whatever the form) as well as other data points (some of it that requires dimension reduction, some that does not) and produce a probability distribution for a classification. Is there a better way to observe this kind of data which captures all of the data points?


title: Determining the statistical significance of differences in a small dataset

text: What would be the best way of determining the statistical significance of the differences in concentration between trials of a dataset such as this one:               +-------+---------------+     | Trial | Concentration |     +-------+---------------+     |     1 | 2.41          |     |     2 | 8.43          |     |     3 | 3.48          |     |     4 | 6.22          |     |     5 | 4.66          |     |     6 | 3.58          |     +-------+---------------+      Would it be a t-test?


title: Definitions of lift contradictory?

text: I understand that lift is a the ratio of the observed support to the expected one if X and Y were independent. Looking at independence it would make sense to define it as follows: $$\frac{\mathrm{P}(X \cap Y)}{ \mathrm{P}(X)\mathrm{P}(Y)}$$ Exactly this is done in the Wikipedia article on lift: Here Yet in another Wikipedia article on association rule learning (and in fact in many other publications on data mining and association rules, like e.g. here) I find this definition:$$\frac{\mathrm{P}(X \cup Y)}{ \mathrm{P}(X)\mathrm{P}(Y)}$$ **My question**   Which definition is correct or are both equivalent (if yes, how)?


title: Differentially "expressed" features

text: When any person gets a dataset of counts for two groups (metagenomics or digital gene expression), the first thing they do is to go and find features that differ significantly between the two groups, which seems obvious. However, isn't it better to first test if the two groups are in fact different, and then find those differences. While the effect of the phenomena (disease, non-disease,etc) might point that the groups are different, the obtained data might say otherwise. Or is it that the test could tell me the two groups are not significantly different, but still go and look for the differential features? Doesn't finding that they are not significantly different saves all the hassle of looking for differentially "expressed" features?


title: Binary input to ROC analysis

text: Im working on assessment of algorithm sensitivity and specificity. I've developed a simulation in order to detect true and false positives and negatives. My intersest is to know if my algorithm is capable of detect a event, or not. So, in the end, I got several yes or no, capable or not capable of detection. Is that possible to run ROC analysis/curves with binary input? Made some tests with ROCR in R without sucess. I.e. 1400 events true positives 600 events false negatives 1900 events true negatives 100 events false positives Thanks in advance


title: Interpret GLM Output (Gaussian - identity) with B0 close to zero

text: I am doing a GLM (guassian - identity) using R. I am modeling different variables (8) to try to understand what creates brand value in an organization.   Some extracts from the output are:               Null deviance: 74.343  on 81  degrees of freedom       Residual deviance: 32.870  on 69  degrees of freedom       (19 observations deleted due to missingness)       AIC: 185.75        I get that one of the coefficients is not significant (B0):               Intercept: Estimate   Std.Err. t value  Pr(>|t|)                0.005959   0.109103   0.055  0.956602       However, some of the other variables are significant. For example:               Consumer_attitude_1: Estimate   StdErr. t value Pr(>|t|)                          0.512708  0.115864   4.425 3.52e-05            Consumer_attitute_2: Estimate   StdErr. t value Pr(>|t|)                           -0.385588  0.152224  -2.533 0.013585       Based on these results, Would it be valid to say that the variables shown as significant in the model may have an impact on brand value even when B0 is close to zero? Thanks very much!


title: Likelihood and estimates for mixed effects Logistic regression

text: First let's simulate some data for a logistic regression with fixed and random parts:               set.seed(1)     n    <- 100     x    <- runif(n)     z    <- sample(c(0,1), n, replace=TRUE)     b    <- rnorm(2)     beta <- c(0.4, 0.8)     X    <- model.matrix(~x)     Z    <- cbind(z, 1-z)     eta  <- X%*%beta + Z%*%b     pr   <- 1/(1+exp(-eta))     y    <- rbinom(n, 1, pr)      If we just wanted to fit a Logistic regression with no random parts, we could use the `glm` function:               glm(y~x, family="binomial")          glm(y~x, family="binomial")$coefficients     # (Intercept)           x      #  -0.2992785   2.1429825       Or constructing our own function of the log-likelihood $$ \log\text{L}(\boldsymbol{\beta}) = \sum_{i=1}^n y_i \log \Lambda(\eta_i) + (1-y_i)\log(1-\Lambda(\eta_i)) $$ where $\Lambda(\eta_i)=\frac{1}{1+\exp(-\eta_i)}$ and $\eta_i=\boldsymbol{X_i' \beta}$ and use `optim()` to estimate the parameters that maximize it, as in the following example code:               ll.no.random <- function(theta,X,y){       beta <- theta[1:ncol(X)]       eta  <- X%*%beta       p    <- 1/(1+exp(-eta))       ll   <- sum( y*log(p) + (1-y)*log(1-p) )       -ll     }          optim(c(0,1), ll.no.random, X=X, y=y)          optim(c(0,1), ll.no.random, X=X, y=y)$par     # -0.2992456  2.1427484      which of course gives the same estimates and maximizes the log-likelihood for the same value. For mixed effects, we would want something like               library(lme4)     glmer(y~x + (1|z), family="binomial")      But how can we do the same with our own function? Since the likelihood is $$ L = \prod_{j=1}^J \int \text{Pr}(y_{1j},...,y_{n_jj}|\boldsymbol{x}, b_j) f(b_j)db_j $$ and the integral has no closed-form expression, we need to use numerical integration like Gaussian Quadrature. We can use the package `statmod` to get some quadratures, say 10               library(statmod)         gq <- gauss.quad(10)     w  <- gq$weights     g  <- gq$nodes      **UPDATE:** Using these quadrature locations $g_r$ and weights $w_r$ for the $r=1,...,R$ ($R=10$ here), we can approximate the integral over $b_j$ by a sum of the $R$ terms with $g_r$ substituted for $b_j$ and the whole term multiplied by the respective weights $w_r$. Thus, our likelihood function should be now $$ L = \prod_{j=1}^J \sum_{r=1}^{R} Pr (y_{1j},...,y_{n_jj} | \boldsymbol{x}, b_j=g_r ) w_r $$ Also, we need to account for the variance of the random part, I read that this can be achieved by replacing the $b_j \sim N(0,\sigma_b^2)$ in our $\eta$ function with $\sigma_j \theta_j$ where $\theta_j \sim N(0,1)$, so in the likelihood function above we actually replace $\theta$'s with $g$'s and not $\beta$'s. One computational issue I don't get is how to substitute the terms since the vectors won't be of the same length. But probably I don't understand that, because I'm missing something crucial here, or misunderstood how this method works.


title: Help me calculate the Precision and Recall and Correctly Classified Instances

text: How to calculate the Precision and Recall for yes, no class and Correctly Classified Instances is 35.7143 %               Data: weather.arff          @relation weather          @attribute outlook {sunny, overcast, rainy}     @attribute temperature real      @attribute humidity real     @attribute windy {TRUE, FALSE}     @attribute play {yes, no}            No.  outlook temperature humidity windy play         1   sunny       85.0    85.0    FALSE   no         2   sunny       80.0    90.0    TRUE    no         3   overcast    83.0    86.0    FALSE   yes         4   rainy       70.0    96.0    FALSE   yes         5   rainy       68.0    80.0    FALSE   yes         6   rainy       65.0    70.0    TRUE    no         7   overcast    64.0    65.0    TRUE    yes         8   sunny       72.0    95.0    FALSE   no         9   sunny       69.0    70.0    FALSE   yes         10  rainy       75.0    80.0    FALSE   yes         11  sunny       75.0    70.0    TRUE    yes         12  overcast    72.0    90.0    TRUE    yes         13  overcast    81.0    75.0    FALSE   yes         14  rainy       71.0    91.0    TRUE    no          === Run information with Weka ===          Scheme:weka.classifiers.rules.PART -M 2 -C 0.25 -Q 1     Relation:     weather     Instances:    14     Attributes:   5                   outlook                   temperature                   humidity                   windy                   play     Test mode:10-fold cross-validation          === Classifier model (full training set) ===          PART decision list     ------------------          outlook = overcast: yes (4.0)          windy = TRUE: no (4.0/1.0)          outlook = sunny: no (3.0/1.0)          : yes (3.0)          Number of Rules  :  4               Time taken to build model: 0 seconds          === Stratified cross-validation ===     === Summary ===          **Correctly Classified Instances           5               35.7143 %**     Incorrectly Classified Instances         9               64.2857 %     Kappa statistic                         -0.3404     Mean absolute error                      0.5518     Root mean squared error                  0.6935     Relative absolute error                115.875  %     Root relative squared error            140.5649 %     Total Number of Instances               14               === Detailed Accuracy By Class ===                         TP Rate   FP Rate   Precision   Recall  F-Measure   ROC Area  Class                      0.444     0.8        **0.5       0.444**     0.471      0.522    yes                      0.2       0.556      **0.167     0.2**       0.182      0.522    no     Weighted Avg.    0.357     0.713      0.381     0.357     0.367      0.522          === Confusion Matrix ===           a b   <-- classified as      4 5 | a = yes      4 1 | b = no      * * * Thanks and best regards


title: Word ranking using Bayes theorem

text: I have problem with find ranking of words while computing multiple contents. I'm splitting each words from the contents except stop words. I have individual word counting with all contents, i.e total contents count=30, word$_1$ count=some $x$, word$_2$ count=some $y$ and so on. We cannot consider word occurrence in separate content. Now we can find probability of word$_1$ with word$_2$, word$_3$, ..., word$_m$ using Bayes theorem. How can we find it?


title: Comparing models of discrete and continuous data

text: I have a data-analysis problem with two closely-related timeseries. Both are binned in time at the same rate, but one is discrete (takes values in $\mathbb{Z}$), and the other is continuous. I am interested in modeling each of these processes' dependence on a common set of predictors, $X$. Given $X$, the question I would like to answer is, > "How well can I predict the discrete series from $X$, compared to my ability > to predict the continuous timeseries from $X$ (and visa versa)?" I'd be interested in knowing if there's any standard way of addressing this problem, or if I should just avoid comparing prediction performance. What I've looked at so far is log-likelihood - a Gaussian log-likelihood for the continuous data, and Poisson for the discrete data. This leads me to another question, the answer to which I suspect to be false (but it's not obvious to me), > "Is it meaningful to compare log-likelihoods from different models? If so, > in what cases?" Hopefully I've left enough details for a meaningful response. Thanks for any input!


title: Two DV's in a 2x2 factorial experiment

text: I'm trying to investigate the effect on (a) tear resistance and (b) gloss of a piece of film due to additive and extrusion method. I can see that the best results for tear resistance is low additive high extrusion, and the best results for gloss is low additive, low extrusion. Although is there a way to say the best overall method - i.e. for best gloss AND tear resistance? Or a way to see if there is an effect between the two DV's - i.e. will the glossier a piece of film is mean that there is less tear resistance? Would a scatter plot with no correlation confirm that there is no effect between gloss and tear resistance? I want to be able to have an overall recommendation stating that the best overall results occurr when additive is at 'x' level and extrusion at 'y' level? I also cannot verify normality (as there are only 3 replications) or homogeneity of variance as SPSS states that there are fewer than two non singular cells. Any help would be much appreciated...


title: Validating previous analysis with new data

text: I am interested in trying to validate some results from a previous experiment. The experiment involves selecting probes of interest that separate our two populations of interest, e.g. smoker vs. non-smoker, high protein diet vs. low protein diet. I have called some probes on an initial data set as significant but what I am interested in is how to statistically validate them in a second set of samples done at a separate time. Is there a methodology or statistical test to say that there is a difference in a probe on the second set of samples given that they have been called different initially. One way (I thought of) is to call differences on this new data and see if there are overlapping results. This however, does not take into account the fact the existing information we have. I am sorry if this is a simple question but I have not been able to find an answer as of yet.


title: Any reference to 20 Newsgroups multi-label dataset?

text: I usually refer this webpage for multi-label dataset. These dataset are in arff format, which can be loaded easily by Weka/MATLAB. Currently, I'm looking for similar ready to use multi-label dataset for 20 Newsgroups data. The original webpage for 20 Newsgroups dataset provides a link to download matlab format file for the file 20news-bydate.tar.gz, which is the single label dataset (it doesn't includes the cross-posts). This page provides the matlab format dataset for the original file 20news-19997.tar.gz, but they are also single label dataset. Can anyone point me to the location where i can get the multilabel matlab/arff format for original 20 Newsgroups dataset? P.S: This is a repost from metaoptimize forum, as i didnt get any response.


title: Recursive (online) regularised least squares algorithm

text: Can anyone point me in the direction of an online (recursive) algorithm for Tikhonov Regularisation (regularised least squares)? In an offline setting, I would calculate $\hat\beta=(X^TX+λI)^{−1}X^TY$ using my original data set where $λ$ is found using n-fold cross validation. A new $y$ value can be predicted for a given $x$ using $y=x^T\hat\beta$. In an online setting I continually draw new data points. How can I update $\hat\beta$ when I draw new additional data samples without doing a full recalculation on the whole data set (original + new)?


title: How to compare two samples without knowledge of the distribution

text: lets say I have the mean of two different measurements of something, a and b. I also have the standard deviation on a, and b, but I do not have access to all the individual measurements of a and b. Is it possible to compare the two values to see if the means are statistically different (i.e., reject the null hypothesis that they are the same)? Any help would be greatly appreciated


title: Using exploratory data analysis

text: I am trying to use exploratory data analysis to decide which model to use with my data for prediction either linear regression/neural networks etc, basically I am focusing on linear and non linear models. What kind of exploratory analysis can I do to decide whether a linear model is good enough for me or a non linear model would be good?


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

text: I have a set of results of independent measurements of some physical quantity. As an example I give here real expermental data on methanol refractive index at 25 degrees Celsius published in scientific literature from 1960 to 2011:               data={{1960, 1.32652}, {1961, 1.32662}, {1963, 1.32650}, {1963,        1.32750}, {1968, 1.32698}, {1968, 1.32890}, {1970, 1.32657}, {1970,        1.32660}, {1971, 1.3260}, {1971, 1.32610}, {1971, 1.32630}, {1971,        1.3350}, {1972, 1.32640}, {1972, 1.32661}, {1973, 1.32860}, {1975,        1.32515}, {1975, 1.32641}, {1976, 1.32663}, {1977, 1.32670}, {1980,        1.3250}, {1983, 1.32850}, {1985, 1.32653}, {1991, 1.32710}, {1995,        1.32621}, {1995, 1.32676}, {1996, 1.32601}, {1996, 1.32645}, {1996,        1.32715}, {1998, 1.32820}, {1999, 1.32730}, {1999, 1.32780}, {2001,        1.32634}, {2006, 1.32620}, {2011, 1.32667}};      The first number is the year of first publication, the second is the published value. The authors do not always provide their own estimate for error of the published value and even when they do this, the published estimate is often not accurate. It is easy to see that the distribution is asymmetrical (other datasets I work with are even more asymmetrical). Here is a density histogram with bins selected by hands:               bins = {1.325, 1.3259, 1.3261, 1.3263, 1.3265, 1.3266, 1.3267, 1.3269,         1.3275, 1.328};     Histogram[dataFirst[[All, 2]], {bins}, "PDF"]      ![histogram](http://i.stack.imgur.com/x4lne.png) Because of asymmetric feature of the distribution I cannot use the mean as an estimate for the true value of methanol refractive index. The fact is that the probability to get smaller value than the true value in general is not equal to the probability to get higher value (the reasons for this are subtle and I will not explain them here). This means that the median is also unsuitable. It may be assumed that each invidual measurement has normal distribution. What is the best way to estimate the true value and its 95% confidence bands for the measured quantity in such cases? P.S. References for relevant scientific papers will be appreciated. P.S.2. The above code is for _Mathematica_ system.


title: Reducing the dimension of an embedding

text: Let $O \in \mathbb R^{p\times m}$ be a data matrix of observations. Suppose we are given a model $\mu : \mathbb R^n \rightarrow \mathbb R^m$ which is able to approximately fit the observations. Fix $l<n$. We wish to find a sub-model $\varepsilon:\mathbb R^l \rightarrow \mathbb R^n$ which will come close to minimizing the following expression: \begin{equation} \sum_{i=1}^p \left[\min_{y \in \mathbb R^l} \Vert \mu(\varepsilon(y)) - O_i\Vert ^2\right]\end{equation} Are there standard approaches (algorithms) to solve this problem? What if $\varepsilon$ is restricted to being linear?


title: SV model estimation in R using tsbugs

text: I have been trying to estimate the basic stochastic volatility model using `OpenBUGS` via R and at an stage of the following command. Please can you comment for the command that can give me the estimated parameter values of the model?               library(tsbugs)     y <- svpdx$pdx          # create and write bugs script     sv0 <- sv.bugs(y)     writeLines(sv0$bug, "sv0.txt")          # estimate model parameters in R2OpenBUGS     library(R2OpenBUGS)     sv0.bug <- bugs(data=sv0$data,                      inits=list(inits(sv0)),                     param=c(nodes(sv0, "prior")$name, "h"),                     model="sv0.txt",                     n.iter=11000, n.burnin=10000, n.chains=1)


title: Assessing correlations, adjusted for time

text: I have time series data for two signals, $X$ and $Y$. At $\text{time}=0$, the both $X$ and $Y$ are roughly 50. As time progresses, $X$ increases to settle at 100 and $Y$ decreases to settle at 0. Hence, the time-series go from (50,50) to (0,100). My question is, if I plot this data, clearly $X$ any $Y$ will appear negatively correlated. Is there a way to _control_ for the correlation that develops just by virtue of the time? Or am I thinking about this all wrong?


title: Ratios of means - statistical comparison test using Fieller's theorem?

text: I would really appreciate any suggestions with the following data analysis issue. Please read till the end as the problem at first may appear trivial, but after much researching, I assure you it is not. The situation is a little complicated because I want to compare the ratios of means: For example, in one experiment I have collected data (Electric current levels, a continuous variable) from 7-8 cells (replicates) which express a particular ion channel gene (Gene1). This gives me a mean current level for Gene1. I have then measured current levels from a different set of 7-8 replicate cells which also express Gene1 PLUS an activator gene (the "treatment"). Now I have a mean current value for Gene1+activator. I repeat the two measurements for Gene2 (a different ion channel gene), again recording mean current values from 7-8 cell replicates for Gene2 alone and then for Gene2+activator. The quantity I am interested in comparing is the %activation or fold activation caused by the activator for Gene1 versus Gene2. So, I would obtain a ratio by dividing the mean current for Gene1+activator by the mean current for Gene1 alone. This would give me the fold activation for Gene1. I would compare this to a similar ratio obtained for Gene2. I have done some research on this and using Fieller's intervals to compute the error bars or confidence intervals seems promising. However, I don't know how to convert that to hypothesis test and get an appropriate p-value for the comparison of means being same/different. Furthermore, the best solution would also allow multiple comparisons and allow me to compare "fold-activations" for Gene1 and Gene2 and Gene3 and Gene4 at the same time. Fieller's intervals seem like the perfect tool to compute 95% confidence intervals etc around each fold-activation but as we are submitting to a journal, they will insist on p-values for our comparisons. As of now, I am reduced to insisting the lack of CI overlap clearly signals a significant difference but I know that lack of 95%CI overlap is an overly stringent comparison test which represents an alpha<0.05. I would truly appreciate any suggestions to the appropriate comparison test (single or multiple comparison, either at this point will be helpful). Thanks in advance.


title: Problem with proof of Conditional expectation as best predictor

text: I have an issue with the proof of > > $E(Y|X) \in \arg \min_{g(X)} E\Big[\big(Y - g(X)\big)^2\Big]$ which very likely reveal a deeper misunderstanding of expectations and conditional expectations. The proof I know goes as follows ( another version of this proof can be found here) \begin{align*} &\arg \min_{g(X)} E\Big[\big(Y - g(x)\big)^2\Big]\\ = &\arg \min_{g(X)} E \Big[ \big(Y - E(Y|X) + E(Y|X) - g(X)\big)^2\Big]\\ =&\arg \min_{g(x)} E \Big[ \big(Y - E(Y|X)\big)^2 + 2 \big(Y - E(Y|X)\big) \big(E(Y|X) - g(X)\big) + \big(E(Y|X) - g(X)\big)^2\Big]\\ =&\arg \min_{g(x)} E \Big[ 2 \big(Y - E(Y|X)\big) \big(E(Y|X) - g(X)\big) + \big(E(Y|X) - g(X)\big)^2\Big]\\ \end{align*} The proof then typically continues with an argument showing that $2 E\Big[ \big(Y - E(Y|X)\big) \big(E(Y|X) - g(X)\big)\Big] = 0$, and hence \begin{align*} \arg \min_{g(x)} E\Big[\big(Y - g(x)\big)^2\Big] = \arg \min_{g(x)} E \Big[\big(E(Y|X) - g(X)\big)^2\Big] \end{align*} which can be seen to be minimized when $g(X) = E(Y|X)$. My puzzles about the proof are the following:   1. Consider > > $E \Big[ 2 \big(Y - E(Y|X)\big) \big(E(Y|X) - g(X)\big) + \big(E(Y|X) - > g(X)\big)^2\Big]$. It seems to me that, independently of any argument showing that the first term is always equal to zero, one can see that setting $g(X) = E(Y|X)$ minimizes the expression as it implies $\big(E(Y|X) - g(X)\big) =0$ and hence > > $E \Big[ 2 \big(Y - E(Y|X)\big) \big(E(Y|X) - g(X)\big) + \big(E(Y|X) - > g(X)\big)^2\Big] = E( 0 + 0)$ = 0. But if this is true, then one might repeat the proof replacing $E(Y|X)$ by any other function of $X$, say $h(X)$, and get to the conclusion that it is $h(X)$ that minimizes the expression. So there must be something I misunderstand (right?).   1. I have some doubts about the meaning of $E[(Y−g(X))^2]$ in the statement of the problem. How should the notation be interpreted? Does it mean > > $E_X[(Y−g(X))^2]$, $E_Y[(Y−g(X))^2]$ or $E_{XY}[(Y−g(X))^2]$?


title: What is a density function?

text: I know about histograms and also know that if we connect the mid-points on the top of bars in a histogram we will get a frequency polygon. This polygon could then be 'smoothed' in a way that it generates a frequency curve like in normal distribution, Poisson distribution, etc. What is smoothing? How is it done? I have heard about the 'filtering techniques' and 'removing noise' but I don't understand the terms. Googling yields things like Kalman filter, moving average filter, low-pass filter, kriging, convolution, etc but they are very abstract concepts for me. Is there any good resource to start learning about these things from scratch?


title: How to plot the likelihood for the location parameter of the Cauchy distribution with a given dataset

text: I am trying to plot the likelihood function of $\theta$ for the $\text{Cauchy}(\theta,1)$ distribution given a specific data set. I am using the code below and it doesn΄t work.               data    <- read.table("mydata", head=F)     theta   <- seq(1, 3, l=100)     density <- function(theta) { prod( dcauchy(data, theta, 1) ) }     plot(theta, density(theta), type="l", xlab="theta",                 main="Likelihood function", col="blue")      Does anyone know why?


title: Is it possible to statistically test relationships between counts and percentages?

text: I have counted adult butterfly numbers over 3 areas and would like to compare these counts with the percentage heather cover on each of the 3 areas. Is this possible? My data looks like this (the cover is the average for each area whilst count is the total number of adults observed over 10 days – let me know if this should be an average instead): $$\begin{array}{c|cc}\rm Area&\rm Cover&\rm Count\\\hline \rm Field\ 1& 53.7& 216\\ \rm Field\ 2&19.5& 2\\ \rm Field\ 3&39.8& 6106\end{array}$$ I suspect from looking at these figures that there probably won't be a relationship, but I would still like to test the relationship nonetheless. I know that there are specific statistical tests for proportion data and other tests for count data because of the way both count and percentage data behave, but I wonder if there is one that will analyse both? I am using the R software environment.


title: Logistic regression with random and with fixed intercept

text: I have made logistic regression with fixed intercept and with random intercept on the same data, and the results (the odds-ratios) are very different. Could somebody tell me what can be the reasons of these big differences between the two regressions?


title: Plotting decay in clustering?

text: I've got a method that produces clusters of elements and singletons. The singletons usually comprise about half of the initial element set, and there usually is a single biggest cluster of a several million elements, followed by a second cluster of a few hundred thousand elements, then the rest of the clusters with thousands to hundreds of elements. I would like to plot this in a way that lets me compare different sets that have initial element sizes. For example, set 1 can have 30 million elements, but set 2 can range from 5 million elements to 50 million elements. How can I plot these characteristics of decay of clustering in a way that is normalized for comparison of multiple sets of different element size?


title: K-fold or hold-out cross validation for ridge regression using R

text: I am working on cross-validation of prediction of my data with 200 subjects and 1000 variables. I am interested ridge regression as number of variables (I want to use) is greater than number of sample. So I want to use shrinkage estimators. The following is made up example data:                #random population of 200 subjects with 1000 variables          M <- matrix(rep(0,200*100),200,1000)         for (i in 1:200) {         set.seed(i)           M[i,] <- ifelse(runif(1000)<0.5,-1,1)         }         rownames(M) <- 1:200              #random yvars          set.seed(1234)         u <- rnorm(1000)         g <- as.vector(crossprod(t(M),u))         h2 <- 0.5          set.seed(234)         y <- g + rnorm(200,mean=0,sd=sqrt((1-h2)/h2*var(g)))              myd <- data.frame(y=y, M)     myd[1:10,1:10]          y X1 X2 X3 X4 X5 X6 X7 X8 X9     1   -7.443403 -1 -1  1  1 -1  1  1  1  1     2  -63.731438 -1  1  1 -1  1  1 -1  1 -1     3  -48.705165 -1  1 -1 -1  1  1 -1 -1  1     4   15.883502  1 -1 -1 -1  1 -1  1  1  1     5   19.087484 -1  1  1 -1 -1  1  1  1  1     6   44.066119  1  1 -1 -1  1  1  1  1  1     7  -26.871182  1 -1 -1 -1 -1  1 -1  1 -1     8  -63.120595 -1 -1  1  1 -1  1 -1  1  1     9   48.330940 -1 -1 -1 -1 -1 -1 -1 -1  1     10 -18.433047  1 -1 -1  1 -1 -1 -1 -1  1      I would like to do following for cross validation - (1) split data into two halts - use first half as training and second half as test (2) K-fold cross validation (say 10 fold or suggestion on any other appropriate fold for my case are welcome) I can simply sample the data into two (gaining and test) and use them:               # using holdout (50% of the data) cross validation      training.id <- sample(1:nrow(myd), round(nrow(myd)/2,0), replace = FALSE)     test.id <- setdiff(1:nrow(myd), training.id)           myd_train <- myd[training.id,]      myd_test  <- myd[test.id,]         I am using `lm.ridge` from `MASS` R package.               library(MASS)     out.ridge=lm.ridge(y~., data=myd_train, lambda=seq(0, 100,0.001))     plot(out.ridge)     select(out.ridge)          lam=0.001     abline(v=lam)          out.ridge1 =lm.ridge(y~., data=myd_train, lambda=lam)     hist(out.ridge1$coef)         out.ridge1$ym     hist(out.ridge1$xm)      I have two questions - (1) How can I predict the test set and calculate accuracy (as correlation of predicted vs actual)? (2) How can I perform K-fold validation? say 10-fold?


title: Problem with Sales Regression Residuals

text: I'm trying to build a model for the ticket sales for different sporting events over a period of 30 days before the game to the day of the game. The problem that I'm having is that I can't seem to fit the right model to it. I've tried negative binomial, linear, poisson, but nothing seems to fit right. A problem that tends to occur is I get a residual plot that looks like the one below. Does anybody have any suggestions for things that I can try? ![Linear Model](http://i.stack.imgur.com/zpkvz.png) Residuals for the negative binomial: ![enter image description here](http://i.stack.imgur.com/G0QyN.png)


title: Quantitative cluster evaluation

text: I have a clustering results by different unsupervised algorithms in the form of coordinates (x,y) and class (it's currently binary). I would like evaluate the quality of clustering by some value. Most of the used techniques assuming that the closely located entries in the plot are related, and if they are from different classes that's not good. I have checked the R packages like clusterSim and clv (cluster evaluatin) by all of them seems to work with original data, but not with the already produced coordinates. Can you please suggest me the quantitative evaluation technique for described situation? PS: My data is imbalanced with ratio around 1:10


title: Why do means&error bars in an ANOVA graph depend on the factors and covariates that define it?

text: I have a mixed-model design with between-subjects factors Group (treatment, control) and Difficulty (easy, hard), and within-subjects factor Timepoint (pre, post). I want to use my statistical package (Statistica) to plot the dependent variable (memory test score) with separate categories for easy&hard, separate lines for the 2 groups, and repeat such a graph for each time point. However, I was surprised to see that the means & error bars in the graphs produced by the ANOVA window depend on how the ANOVA is defined: 1) if I do each between-subjects ANOVA separately (once for each level of the within factor, i.e. once for pre and once for post), the means & error bars showing the interaction in each ANOVA are different from the graph showing the 3-way interaction (browen down for pre and post) in the mixed-model 3-way ANOVA 2) if I take IQ as a covariate, I get means & error bars that are again different to the no-covariate case 3) if I enter the covariate into the model crossed with the other factors (e.g. if I look for a "Group X Difficulty X IQ" interaction), I get means & error bars that are yet-again different as compared to when the covariate is added to the model separately, i.e. uncrossed with any other factor. Why is that? Shouldn't the mean and error bars simply be calculated as the average and SEM (or whatver the error bars represent) respectively? And shouldn't those two statistics be independent of what factors are included in the model (ANOVA) and _how_ they are included (see points 1-3 above)? If it is in fact legitimate that these things influence the means and error bars, then which graph should I believe? In other words, how should I plot my two between-subjects (Group x Difficulty) interactions: from within the large 3-way (Group x Difficulty x Timepoint) ANOVA or from within each of the 2-way ANOVAs (Group x Difficulty) done separately for pre and for post? Thanks and sorry for the long question! Please let me know if I need to make it more clear.


title: Adjusting a Logistic Curve through R

text: I have data of certain proportion $p$ in a population. I need to approximate the behaviour of this population through a Logistic curve of the form: $p=\frac{1-e^{-\gamma t}}{1+ \alpha e^{-\beta t}}$ I have consulted several webpages and I have found out that the more simple process is to estimate the logistic curve as $p=\frac{1}{1+ \alpha e^{-\beta t}}$ And then estimate the value for $\gamma$, but I still don't know how to do it because when I use a logistic regression the range of the dependent variable is one and zero, not all the range between them. I have found a document about it, but it's implemented in MatLAB, how can I do it in R? http://home2.fvcc.edu/~dhicketh/diffeqns/activities/logistic.pdf


title: Deep learning: what is done with sampling in practice?

text: I have come against a small issue when implementing some deep learning techniques. Sometimes, if one is implementing a classic binary-binary Restricted Boltzmann Machine (RBM), the sampling functions (i.e. particularly the function that samples the hidden units given the visible units, $p(h|v)$) requires using a binomial distribution. However, I have found in my code that it's entirely possible to generate a row in matrix using the binomial distribution (using the visible means computed from a propagate down function as the p value and $n = 1$) with all 0's or all 1's. This can lead to certain numerical problems (akin to numerical analysis issues such as NaN's appearing in the computation given that later matrix operations use division or Log) in later computations. My question is: What is done in practice (obviously most deep learning papers to my knowledge, as well as tutorial sites, do not address these small technical issues) in such cases? A binomial distribution sampling could generate a row in a matrix with all 0's. What I've done is decided that I need to check for such cases and add some very slight noise to one of these 0's (let's say I write a small sub- routine to check that a matrix doesn't contain any rows of all 0's and will then randomly select a cell in any such bad-case-row and then, insert a noise value such as $10^{-2}$). I know I'm giving up some accuracy in introducing some slight noise to nudge the computation in the right direction, but once I do I get correct predictions/label constructions (for all cases) so it appears my code works fine. So while I believe I may have found one technical workaround, I would like to know what is actually done to handle this scenario? _Update_ : Is this a technical challenge no one really encounters? Is there something special about Python (and usu. Theano) that magically ensures there are never any sample rows in a sampled matrix that are all 0's or all 1's? even though a binomial sampling is used)?


title: Testing for independence and strength of association

text: I am having difficulty at formulating properly a null hypothesis regarding the following setup: Suppose that we are interested in probing a possible association between bundled products sold at supermarket chains vs online. We have P products, and only bundles of two products are available for sale. Thus, there are P /2 * (P-1) possible bundles. At the end of the month, when we look at consumer logs for **online** sales, we see that consumers purchased **O bundles** (number of bundles). Similarly, when we check the logs for sales at chains, we observe that **C bundles** were sold. Additionally, we notice the intersection between **O bundles** and **C bundles** , which is **S bundles** ( S <= O,C) that were sold both online and at the chains. Now, we would like to know two things: (1) whether observing a particular bundle of two products (Pi, Pj) sold online is associated or non-independent of observing the same bundle to be sold at chains? (2) if there is an association, how can we assess the strength of dependence or association? Any pointer / comment would be very much appreciated. Thank you.


title: If $X \sim {\rm Bernoulli}(p)$, what is the distribution of $Y=aX-b$?

text: Let $X \sim {\rm Bernoulli}(p)$ and $Y=aX-b$. I want to find the distribution of $Y$. I am assigned a homework problem with specific $a$ and $b$. My textbook covers methods for solving this, given that x is continuously distributed. However, $X$ is discretely distributed. I'm not looking to be shown exactly how to do this. I am just unsure of where to start, as this is not covered in our text. Any help you could provide to me would be greatly appreciated.


title: Understanding what Dassault iSight is doing?

text: I've been tasked with taking up some modeling (perhaps predictive) for chemical formulations and some resultant performance characteristics. I have experience with finite element analysis (mechanical and fluid) and some more basic statistics, which identified me as a candidate for taking this on. I'd like to learn more in order to better understand this request. I'm being trained by someone and it seems the software my company is using for this is Dassault iSight. My questions:   * I'm a bit disadvantaged as it was hard for me to search this site without even understanding the proper terminology. I'm hoping someone with iSight experience (or from scanning the page) could let me know what "class" of software this is -- neural network? Machine learning? Bayesian modeling? Predictive analysis?   * Might I be pointed to some resources to learn more about this type of modeling so I can understand what iSight is doing on a more theoretical level?   * If someone knows iSight, are there other packages that aim to do the same? I'm not sure how they chose this software package; if something else would be better suited (preferably in R, which I have experience with), I'd love to know about it. For reference and without getting very specific we've conducted a DOE of various chemical formulation ratios for a material and have five or so responses we've tested on each of our ~80 lots. Now we're looking to determine the interactions between the X's (formulation components) as well as how the X's affect the Y's (performance attributes). I realize this is a total noob/ignorance question... just looking to understand the proper terminology and direction for learning more as it's been hard to search on my own.


title: Predicting Naive Bayes model in R on a test data with a single record

text: I built a naivebayes model using the Housevotes84 data(discrete data) in mlbench package-               model <- naiveBayes(Class~., data=HouseVotes84)      I took one record from the data to predict the model. Irrespective of which record it is the model always provides the output- "Democrats".               housevotes_test <- read.csv("path")          predict(model, housevotes_test[1,])     [1] democrat      I then tried a test data with 2 records- 1 wrt republicans and the other wrt Democrats. I find that while predicting the model on this test data the model classifies accurately. Here are my questions-   1. Why does the model always classify the test record into "Democrats" when only a single record is given in test data?   2. Is the classification dependent on the target variable? Since it classifies the test records accurately when the data consists of atleast 1 record pertaining to "democrats". I am not sure if the code is wrong. So please point out where the error is or provide an explanation for this mis-classification when only 1 record exists in test data. Please try this experiment and please help me with finding an answer to this.


title: Reversible jump MCMC code (Matlab or R)

text: Does anyone know of some well written code (in Matlab or R) for reversible jump MCMC? Preferably a simple demo application to compliment papers on the subject, that would be useful in understanding the process.


title: Can probability be used to calculate (any?) effect size measurement?

text: For example, some studies report the probability of survival, on a scale of 0-1, with SE, and sample sizes for two treatments, rather than the mean and SD. Is there a sensible way to use the probabilities directly or convert to a useful format?


title: Interpretation of Output

text: Here is the R code that produced the output below:               library(caret)     set.seed(934)     fitControl <- trainControl(method= "repeatedcv", number=10, repeats=10)     logitfit <- train(z~ a+b+c+d+e, data=train, method="glm", trControl = fitControl)     logitfit      Consider the following output from R:               Generalized Linear Model           900 samples      20 predictors       2 classes: '0', '1'           No pre-processing     Resampling: Cross-Validated (10 fold, repeated 10 times)           Summary of sample sizes: 846, 847, 846, 845, 847, 846, ...           Resampling results            Accuracy  Kappa  Accuracy SD  Kappa SD       0.80     0.58  0.076       0.18          This is a logistic regression model using a 10-fold cross validation repeated 10 times. Does the output above indicate that the average accuracy is 0.80?


title: Ordinal model comparison: lower deviance but lower accuracy, is this possible?

text: I am performing ordinal regression (6 response categories) on a large dataset (>15000 points) and comparing two models, one with one predictor X1 (Model 1) and the other with two predictors X1 and X2 (Model 2). The models I fit are also cumulative over the categories: as in, I am looking at outcomes being _at least_ at a certain level - and the probability of being at level 1 is always 1. The likelihood ratio test reveals that Model 2 provides a significantly better fit to the data when compared to Model 1. However, I am also computing the models' respective accuracies, by summing up the number of correct estimations for each category (again cumulative) over the total number of outcomes, and find that accuracies are lower for Model2 than for Model 1. X2 has got missing values, so I tried to run the regression again on a subset of the data which had values for all X1 and X2, when I do this the accuracy remains the same. I am very puzzled by these results. How can the "better" model have a lower accuracy? I checked my code a hundred times for programming issues without success, so now I am starting to wonder whether there might be another possible reason for this. Does missing data in a predictor affect the accuracy - without affecting the deviance? Are there any issues with large datasets when it comes to computing these measures?


title: Drawing balls with two features

text: Suppose that I have a total of 1000 balls in a bag. 400 balls are labeled `a` and remaining 600 balls are labeled `b`. In addition, within the same 1000 balls, 60 balls are `white` and remaining 940 balls are `black`. Now, if I conduct an experiment of picking a ball one after other without replacement, lets say 30 times, how can I compute the probability that I would have at least 20 balls that are `white` and labeled `a`? I would like to basically know how probable is this event of picking 20 balls that are `white` and labeled `a` out of 30 events, given this setup, due to chance. And its not a homework question. Thank you!


title: Hierarchical multinomial logit with R/JAGS

text: I am working on a small project where I have to do a Choice Based Conjoint (CBC) analysis. In order to get the part-worths for the different consumers I need to estimate a hierarchical multinomial model. Does anyone know how to do this in R/JAGS, and how to save the estimated coefficients for each consumer in a .csv file?


title: Expected number of failures preceding the first success

text: I am a beginner in statistics and probability. I have a question that says > What is the expected number of failures preceding the first success in an > infinite series of independent trials with the constant probability of > success equal to $p$? I have tried of solution of this and is not quite sure of this: The probability of failure is $1-p$ The probability of $x$ failures in a row is $(1-p)^x$. Now, expectation of $x$ is $E(x) =\sum_0^\infty x (1-p)^x\\ =(1-p) + 2(1-p)^2+3(1-p)^3 +\ldots \infty \\ =\frac{1-p}{p^2} $ Thanks in advance.


title: Compute the power of a matrix in R

text: (1) I am looking for a package for computing the power of a matrix. If you have some good recommendation please let me know. (2) I searched on the internet and followed what some said to install a package called "Malmig" in R but after selecting the mirror site, it failed: > In install.packages("Malmig") : package ‘Malmig’ is not available Some idea why? Thanks!