title: Fitting a copula with Poisson marginals to data in R

text: First off, I know this is a question which requires an thorough answer, so I am coming here with a very humble attitude. I have limited knowledge about both copulas and R, so I will try to explain what I know and what my thoughts are: I have a data set of n observations of (X,Y), two correlated Poisson random variables, which comes in "set" - that is - I observe (X_1,Y_1), (X_2,Y_2),...,(X_n,Y_n). I have read about using discrete marginals with copulas in: http://www.actuaries.org/LIBRARY/ASTIN/vol37no2/475.pdf and I realize it can be troublesome using Poisson marginals. I have to admit that I did not understand too much of the above article, but nevertheless I hope I am able to describe the joint distribution of (X,Y) using some copula. I understand that different copulas model different type of dependence, such as tail dependence and such. But this is not what I am looking for. For my data, it seems that X follows Y, and vice verca - so it seems like we have strong dependence on the "straight line". I found the following about fitting copulas to data:                   library(copula)           ## Toy example for gumbel copula with log-normal distribution      ## (Taken from the documentation of copula::fitMvdc)           ## Specify the copula      gumbel.cop <- gumbelCopula(3, dim=2)      myMvd <- mvdc(gumbel.cop, c("lnorm","lnorm"), list(list(meanlog = 1.17),                                              list(meanlog = 0.76)))      ## Generate some random sample to test      x <- rmvdc(myMvd, 1000)           ## Fit the random sample      fit <- fitMvdc(x, myMvd, c(1,1,2))      fit       Now, first of all this uses a gumbel copula. Could there be some other copula better fitting what I am looking for? I am curious to how 'the experts' choose which copula to work with. For example, if we are looking at financial data with strong tail dependence, it seems obvious, but for other 'types' of dependence, I am not sure if I understand the copula choice. From the above code in R, could I just replace "lnorm" with "name of poisson distribution", and what would this name be? I cannot seem to find a list of names of the marginals. Also, the above code generated data. I have my own data. I know I am asking a lot here, but this is all a bit overwhelming to me... Thank you very much for any reply!


title: Efficience of Expectation-Maximization algorithm in function of learning dataset size

text: I have datasets of increasing sizes identically distributed. I have tried to fit a gaussian mixture to these datasets using Expectation-Maximization algorithm. To check the quality of this fit, I have run goodness of fit tests between the density of the known points and the estimated density. It provides very good results. Tests hypothesis are accepted in more than 90% of the cases. The acceptation rate first increases with the number of points in the learning dataset (as I was expected). And then for the last dataset size, the acceptation rate falls from 95% to about 85%... I have run tests many times and it is always the same... Do you have some hints about what is going on ? Why would the quality of the model be so low (compared to other sizes) for one particular learning dataset size ?... Do you know articles where the behavior of Expectation-Maximization is studied in function of the sample size ? **EDIT :** As I can simulate as many datasets as I want. I do not use cross- validation. To test the quality of the density estimated with a dataset $A$, I run the goodness-fo-fit tests between a test dataset $B$, simulated in the same way that dataset $A$ and a sample of points simulated from the estimated gaussian mixture. The size of the dataset $B$ is constant through my tests and thus does not depend on the size of dataset $A$.


title: How to Test Independence of Poisson Variables?

text: I have a data set of length $10,000$. Each entry, $e_i$, of my data set is associated with two values, $x_i$ and $y_i$. Let $X$ and $Y$ represent the random variables for which my data set is a realization. I have found and tested from the data that the variables $X$ and $Y$, when considered separately, have Poisson distributions and I have used the data to estimate their parameters $\lambda_X$ and $\lambda_Y$. I believe that the two variables are nearly independent, which would give me the joint distribution of the variables--which is what I'm looking for. How do I go about testing for independence when my hypothesis is that they are independent? I was thinking about computing correlation but there are cases in which uncorrelated variables aren't independent.


title: Analyzing 3D data: What can be done?

text: I am new to this kind of analysis, and I want to know what values I can look at in 3D data. The data itself is a 3D volume $(x,y,z)$ with a floating point value in every coordinate. It is a _hyperspectral image_ , meaning: the $z$-space is the same $xy$ plane at a different wavelength. Here is an example from the internet of a hyperspectral image. The goal is to find **characteristics** in the data that can be exploited to compress it. How could I approach this? [See also PCA of hyperspectral image data ]


title: Error Bars for Monte Carlo Experiment

text: Suppose we have a random variable $X$, where $\mathbb{E}(X)$ and $\text{Var}(X)$ are known. I have computed $N$ number of MC-type samples from the distribution of $X$. Let $\bar{x} = \frac{1}{N}\sum x_i$ where $x_i$ represents a single sample from our distribution. From CLT, I know that the error bars for $\bar{x}$ can be constructed in the following way, assuming 90% CI: $$ \bar{x} \pm (1.64)\sqrt{\frac{\text{Var}(X)}{N}}$$ But I'm interested in computing the error bars for the relative error $\frac{|\bar{x} - \mathbb{E}(X)|}{|\mathbb{E}(X)|}$. Since the relative error itself is a random variable, I'm having difficulty constructing error bars for it . Any suggestions? Thanks in advance! I wonder if I can use relative error in Variance and substitute that with $\text{Var}(X)$.


title: Drawing a random Sample on a Probability Distribution

text: Say I have multiple normal or beta distributions. So I have two questions.   1. How do I Draw a random sample from a distribution?   2. How do I compare determine which of my distributions has the largest draw sample? More Info: Im trying to decipher this blog on Thompson sampling and Multi armed bandits. In it the author describes each arm as having a Beta Distribution. How do I run a random sample on each distribution and what is meant by largest drawn sample.


title: GLM poisson regression and some independent variables are itself poisson distributed

text: How can we deal with in GLM regression, when some variables of the independent variables (the x) are itself poisson distributed? Suppose we have two random poisson distributed variables X,Y with latent $\lambda1$ and $\lambda2$ where $\lambda1$ = a * $\lambda2$ + b and we only know their realisations. How do we get their linear dependency?


title: Handling Missing Values During Test Phase

text: I was searching for methods for handling missing values in case of Regression task. There are already few threads but I couldn't find what I was looking for. Suppose I have 4 independent categorical variables (A,B,C,D), each of them can takes an integer value and their limit is pre-defined. For example A can take values from 1-10, B from 2-7, etc..All are categorical variables. I want to train a regression model on these variables. Suppose during the test phase, if I don't have the value of a particular variable, how should I approach the problem. I was thinking, in lines of having priors for each variable from the training data and using them during test phase. Please share any better methods. Thanks


title: Proof: Linear Regression

text: I am at a loss how might one prove $\text{Var}{(\hat{y}_h)} = \sigma^2 \left(\frac{1}{n} + \frac{(x_h-\bar{x})^2}{S_{xx}}\right)$. Note that $\hat{y}_h$ = $b_0 + b_1X_h$ which is a regression line estimate at some given $X_h$. Moreover, what is the second term actually saying? I understand the first is just the variance divided by the degrees of freedom but the second is more complex it seems.


title: Fitting missing points to a dataset

text: I've been assigned to fit some missing values into a large dataset, and I've come across a problem. I've finally got a function describing my data (for simplicity say y=2x^2). Now, for every value of x, I have multiple values of y (so there is variance in the data). I'm using R, and used the `nls(y~ax^b)` function to find the relationship (returns y=2x^2). I can use `deviance()` to find the deviance of my fitted model. My question is, how do I use the deviance and the function to assign values to missing data points? That is, for some x values, I'm missing y. I don't think I can simply put in the value of y that my function returns, since this does not incorporate the variance. For example, if I had missing values for x={1,1,2,3,3,3,3,3}, fitting y={2,2,8,18,18,18,18,18} does not seem right. So how to I incorporate the variance to fit values of y such that my fitted values will conserve the measured variance?


title: Proving $\text{Var}{(\hat{y}_h)} = \sigma^2 \left(\frac{1}{n} + \frac{(x_h-\bar{x})^2}{S_{xx}}\right)$

text: I have asked in another question how $\text{Var}{(\hat{y}_h)} = \sigma^2 \left(\frac{1}{n} + \frac{(x_h-\bar{x})^2}{S_{xx}}\right)$. Note that $\hat{y}_h$ = $b_0 + b_1X_h$ which is a regression line estimate at some given $X_h$. This question concerns why the term $Cov(b_0,b_1)$ alone yields the RHS. Substituting $b_0 = Y - b_1X$ we get that $Cov(Y,b_1) - XCov(b_1,b_1)$ = $Cov(\frac{\sum{Y_i}}{n},\sum k_iY_i) - XVar{(b_1)}$. Here X and Y without subscript are arithmetic means. We can then rearrange to obtain $\sum \frac{k_i Var(Y_i)}{n} - \frac{X\sigma^2}{S_{xx}}$ which quickly yields the desired result. My question is, why does this work? This single term does not seem like it should alone yield the RHS. Have I made an error in algebra?


title: Can OLS be considered an optimization technique?

text: Can ordinary least squares estimation be considered an optimization technique? If so, how can I explain this? Note: From an AI perspective, supervised learning involves finding a hypothesis function $h_\vec{w}(\vec{x})$ that approximates the true nature between predictor variables and the predicted variable. Let some set of functions with the same model representation define the hypothesis space $\mathbb{H}$ (That is we hypothesise the true relationship to be a linear function of inputs or a quadratic function of inputs and so forth). The objective is to find the model $h\in\mathbb{H}$ that optimally maps inputs to outputs. This is done by application of some technique to finds optimal values for the adjustable parameters $\vec{w}$ that defines the function $h_w(\vec{x})$. In AI we call this parameter optimization. A parameter optimization technique/model inducer/learning algorithm would for example be the back propagation algorithm. OLS is used to find/estimate for $\beta$ parameters that defines the linear regression line that _optimally_ maps predictor variables to output variables. This would be parameter optimization in the scenario above.


title: Solving a regression problem

text: I wanted to solve such a regression problem: $$Y = Xb + e$$ where $X$ is a $m$ by $n$ matrix, resulting in: b = (X'X)-1X'Y as a solution. Since $n$ is quite large (2400), I can't use the conventional methods to calculate the inverse of $X'X$. So I thought about using LU Decomposition using Crout method. In this link an implementation of this method is available which requires an $n$ by $n$ matrix $A$ to solve the equation $Ax=B$. Does it mean that I should use $X'X$ matrix as the input ($A$)? Please note that the problem is that it requires an $n$ by $n$ matrix instead of an $m$ by $n$ matrix.


title: Expected value of non-standard quadratic form

text: Say I have two $n$-vectors $f(t)=(f_1(t), \dots, f_n(t))$ and $f(s)=(f_1(s), \dots, f_n(s))$, both with expected value 0.   Let $\operatorname{Cov}(f(t)) = a \Sigma$ and $\operatorname{Cov}(f(s)) = b \Sigma$ with scalar $a,b$. Further assume that the two vectors are not independent and that $Cov(f_i(t), f_j(s)) = c \Sigma_{ij} \quad \forall i,j$.   I know that $\operatorname{E}(f(t)' \Sigma f(t)) = a \operatorname{trace}(\Sigma^2)$, but can any of you help me derive an expression for $\operatorname{E}(f(t)' \Sigma f(s))$? Background: This is not homework (i'm too old for that...), I'm working on an estimator for the temporal covariance structure of spatially correlated functional data.


title: Comparing player percentages effectiveness with different sample sizes

text: I have a spreadsheet of player data. One of the categories of player data is pass percentage. Pass percentage is calculated as follows:               AccuratePasses/TotalPasses      In total I have around 400 players to compare, however the sample size of passes between each player is very different. For example on player may have executed 1592 passes (1426 of which were accurate), and another may have executed only 7 passes (of which all were accurate). How can I use the data I have to extrapolate a score or ranking based on the player effectiveness or likelihood of being effective. At the moment I have a formula to do this as follows:               PlayerAccuratePasses-(PlayerTotalPasses*(AverageAccuratePasses/AverageTotalPasses)      I did not pull this formula together and want to be sure that this makes sense mathematically to produce a score of who is the most effective passer.


title: How to tell the "closeness" of two variables

text: I have completed an experiment measuring the magnetic field of a solenoid. The answer I got was $0.0075 \pm 0.0011$T. The uncertainty by the way was generated from Microsoft Excel's Regression under data analysis and is the standard error of the gradient of some graph. The gradient is the magnetic field, B in Tesla. I would like to compare this with a theoretical value obtained from an approximation formula of $B = 0.0106623 \pm 0.0000007$T. How do I determine whether my experimental result is "close" to the theoretical one?? I would imagine it being easy if my theoretical result had no uncertainty. Then I could just say that the theoretical result is about three standard errors away. But now the theoretical value has an uncertainty. How do I get across this? Thanks.


title: Maximum likelihood estimator of discrete distribution with known marginals

text: I would like to compute the maximum likelihood estimator of a general two- dimensional discrete probability $p_{ij}$ with $i=1,\ldots,N$ and $j=1,\ldots,N$ given a set of counts $c_{ij}$ and the marginal probabilities $\sum_j p_{ij} = p^{(x)}_i$ and $\sum_i p_{ij} = p^{(y)}_j$. The values $c_{ij}$, $p^{(x)}_i$ and $p^{(y)}_j$ are inputs. I would like to know if there is a closed form solution to the above problem. I have worked on the problem using Lagrange multipliers. The objective function is $$ F\left(\\{p_{ij}\\}, \\{\lambda^{(x)}_i\\}, \\{\lambda^{(y)}_j\\}\right) = \sum_{ij} c_{ij} \log p_{ij} - \sum_{i} \lambda^{(x)}_i \left(\sum_j p_{ij} - p^{(x)}_i\right) - \sum_{j} \lambda^{(y)}_j \left(\sum_i p_{ij} - p^{(y)}_j\right). $$ The $p_{ij}$ derivatives result in the equations $p_{ij} = \frac{c_{ij}}{\lambda^{(x)}_i+\lambda^{(y)}_j}$ for all $i,j$. Substituting these equations into the $\lambda$ derivative conditions gives $$ \sum_j \frac{c_{ij}}{\lambda^{(x)}_i+\lambda^{(y)}_j} = p^{(x)}_i \text{ for all $i$} $$ and $$ \sum_i \frac{c_{ij}}{\lambda^{(x)}_i+\lambda^{(y)}_j} = p^{(y)}_j \text{ for all $j$}. $$ I do not know how to solve the above equations except by numerical methods. I would like a closed form if possible.


title: Unequal sample sizes: When to call it quits

text: I'm peer reviewing an academic journal article and the authors wrote the following as justification for not reporting any inferential statistics (I deidentified the nature of the two groups): > In total, 25 of the 2,349 (1.1%) respondents reported **X**. We > appropriately refrain from presenting analyses that statistically compare > group **X** to group **Y** (the other 2,324 participants) since those > results could be heavily driven by chance with an outcome this rare. My question is are the authors of this study justified in throwing in the towel with respect to comparing groups? If not, what might I recommend to them?


title: Box and whisker plot or bar plot with mean and SE?

text: I am trying to present my data in a meaningful way. To do so, I thought it is a good idea to use a box and whisker plot, because, it shows min, max, median and skewness of data. However, due to variety in my data, I am getting a lot of outliers. For example, the many little circles are outliers. Is this bad? ![enter image description here](http://i.stack.imgur.com/GIfaD.jpg) How about using a Standard Errors bar plot instead of a whisker bar plot as below? ![enter image description here](http://i.stack.imgur.com/Rj7RU.gif) In which case Standard Errors bar plot is preferable?


title: Clarification on LDA and the multivariate Gaussian

text: From my understanding, to calculate the posterior probability of a sample $x$ belonging to a class $k$ using Linear Discriminant Analysis you would first calculate the eigenvector matrix $W$ required to transform any given sample into the $p$-dimensional discriminant space. Then, using $W$, you would project $x$ and the mean for the class $k$ training data ($\mu_k$) into the discriminant space to get $x'$ and $\mu_k'$. You would then use the results of these projections to calculate $P(x|k)$ via the multivariate Gaussian: $$P(x|k)=\frac{1}{\sqrt{(2\pi)^p|\boldsymbol\Sigma|}} \exp\left(-\frac{1}{2}({x'}-{\mu_k'})^T{\boldsymbol\Sigma}^{-1}({x'}-{\mu_k'}) \right)$$ Next, you would use Bayes Theorem to estimate the proportionate probability $\widehat{P(k|x)} = P(k)P(x|k)$ given $P(k) = N_k / N$ with $N_k$ being to total number of training samples in class $k$ and $N$ being the total number of training samples. Finally you would normalize to get an actual estimate of $P(k|x)$ via: $$P(k|x) = {\widehat{P(k|x)}} / {\sum_{i\epsilon N} \widehat{P(i|x)}}$$ In the answer to Linear discriminant analysis and Bayes rule it is noted that if the pooled within-class covariance matrix is used for $\boldsymbol\Sigma$ then $|\boldsymbol\Sigma| = 1$ and the squared Mahalanobis distance $d = ({x'}-{\mu_k'})^T{\boldsymbol\Sigma}^{-1}({x}-{\mu_k'})$ simplifies to become the Euclidean distance $d = ({x'}-{\mu_k'})^T({x}-{\mu_k'})$. **I'm a little confused as to how this simplification works and have been unable to find any good references that demonstrate it.** Also, would that mean $\boldsymbol\Sigma = \boldsymbol{S_w}$ used to calculate the discriminant space or is it something else altogether? Any help would be greatly appreciated.


title: Lilliefors-Test and Jarque-Bera-Test results interpretation

text: I have conducted both of the above tests via the commands:               lb <- lillie.test(x)     jb <- jarque.bera.test(x)      I think `library(nortest)` is required to do this. I want to test `x` for normality. `x` is a simple time series from t = 1 to t = 380 (I cannot hand out the data, but I can provide the results). The results are: lillie.test: `0.06***` (highly significant) jarque.bera.test: `28.97***` (highly significant) So can anyone help me with the interpretation? I did some searching and I also read documentations, but this only confused me more :( What is the Hypothesis of these tests and what do these results tell me? I thought the hypothesis was normally distributed data, but I am not sure anymore.               kurtosis(x)     [1] 3.28458     skewness(x)     [1] -0.6423145      Thanks.


title: Why is my coefficient so large?

text: I have a GLM modeling basketball data where the dependent variable is the proportion of wins for a season. The independent variables are either percentages, or numeric. The family is binomial, and the link is logistic. Why is the "fgp" coefficient so large? ![enter image description here](http://i.stack.imgur.com/uVp4d.png) Especially compared to "ftp," since the both have similar relationships with the dependent variable, winp: ![enter image description here](http://i.stack.imgur.com/PsAOY.png) Here are some diagnostic plots for reference: ![enter image description here](http://i.stack.imgur.com/Mx62Z.png) ![enter image description here](http://i.stack.imgur.com/0Mg5g.png) Also, I know I'm supposed to model the DV as a proportion of successes over total trials. I read the R command for this is: DV <\- cbind(myData$wins, myData$losses). Shouldn't it be: DV <\- cbind(myData$wins, myData$totalTrials)


title: Name of algorithm (or paper) that scikit-learn cluster.estimate_bandwidth() function implements for meanshift bandwidth selection

text: Can someone tell me the name of the algorithm (or paper) that sklearn.cluster.estimate_bandwidth implements and is used by the meanshift algorithm implemented in Scikit-Learn to automatically select its bandwidth parameter? Thanks. **EDIT** : From the link to the code given by Franck, it seems that the bandwidth sigma is estimated for a set of feature vectors X using the average of the furthest distance between vectors in X to their k nearest neighbors. These vectors (whose pairwise distance we are considering) are either sampled randomly or computed using all the vectors in X in batches, depending on the parameters given.


title: When was the k-means clustering algorithm first used?

text: K-means is probably one of the most used algorithms for clustering. I was looking for bibliography for its first use, but it has been around a lot, so what's the first one? Also, when was the algorithm first named "k-means"?


title: MCMC and data augmentation

text: I have been looking at an MCMC data augmentation question; the general form of the question is as follows: Suppose data gathered on a process suggests $X_{i} \sim \text{Pois}(\lambda)$ and a prior for the rate parameter is suggested as $\lambda \sim \text{Exp}(\lambda_{0})$. The data is recorded and presented in a typical form (i.e. the number of occurrences of each value for $X_{i}$ from $0$ to $n$), however, the data gathered does not discriminate in instances where $X_{i} \leq 1$ (i.e. all occurrences where $X_{i} = 0$ and $X_{i} = 1$ are grouped into one category). Given the data, the likelihood and the prior described above, the question asks for:   * The posterior form of $\lambda$,   * The number of occurrences where $X_{i} = 0$. I'm not really sure how to answer this question, but I am aware that Gibbs Sampling can be used in data augmentation. Does anybody have any information on how this could be done? **EDIT:** I should specify that it's primarily the second part (the number of occurrences where $X_{i} = 0$) that I'm unsure about. For the first part (the posterior form of $\lambda$), given the likelihood and the prior suggested, I have reasoned (although I'm happy to be corrected): Given: $$ \pi(\lambda|\vec{x}) \propto p(\vec{x}|\lambda) \times p(\lambda) $$ So, for the model given above: $$ \pi(\lambda|\vec{x}) = \frac{\lambda^{\sum_{i=1}^{n}x_{i}}}{\sum_{i=1}^{n}x_{i}!}e^{-n\lambda} \times \lambda_{0}e^{-\lambda \lambda_{0}} $$ Simplifying yields: $$ \pi(\lambda|\vec{x}) = \frac{\lambda^{\sum_{i=1}^{n}x_{i}}}{\sum_{i=1}^{n}x_{i}!}e^{-\lambda(n+\lambda_{0})}\lambda_{0} $$ which is proportional to (and hence the posterior form is given by): $$ \pi(\lambda|\vec{x}) \propto \lambda^{\sum_{i=1}^{n}x_{i}}e^{-\lambda(n+\lambda_{0})}\lambda_{0} $$


title: Reference on examples with R codes for Bayesian simulation based methods of posterior approximation

text: I have been trying to learn the Bayesian simulation based methods of posterior approximation. Although the theories are now quite clear but I am seeking for some examples with R codes so that I can practice and implement to my own work. I am particularly looking for the examples with R codes of the non-iterative and iterative Monte Carlo methods. Mainly-   1. Importance sampling    2. Rejection sampling   3. Metropolis algorithm   4. Metropolis-Hastings algorithm    5. Gibbs sampler So, is it possible to find examples with R codes of these methods of simulation somewhere collectively? I mean it is a kind of heuristic learning, so basically I need some tutorial or some blog or website where the examples are properly described with corresponding R codes. **Edit:** I am sorry I forgot to mention, I was actually looking for creating functions in R for these sampling methods to know the basic of these algorithms. Not looking for any packaged software.


title: Confidence estimation for data points in a recommender system

text: I have a 100 by 100 matrix. Each cell is either 0, 1, or missing (denoted NA). Rows denote 'users' and columns denote 'items'. My goal is to impute the missing values, and provide a confidence level for each imputed value. To impute the missing values I have implemented a set of recommender systems including naive Bayes and k-Nearest Neighbours. I have calculated the mean accuracy (miss-classification rate), and standard deviation of this rate, for each recommender system, using 5-fold cross validation. I would next like to calculate a level of confidence for each data point. However, I have only one data point for each user-item pair, so more well- researched methods are difficult to apply. My solution is the following. I have run the cross validation process 20 times, each with a different splitting of the data set. This has provided 20 predictions for each data point, each based on a different subsample of the data. I have then calculated an estimated confidence for each data point as the following: Number of correct predictions/Total number of predictions. I interpret this as the probability of a correct prediction for a particular pair. Is this a principled metric to estimate confidence in this situation, and are there any other approaches that you recommend I look into?


title: How to perform exploratory factor analysis on associative network?

text: In an article by Teichert and Schontag ("Exploring Consumer Knowledge Structure Using Associative Network Analysis", 2010), the authors perform (page 387) an exploratory factor analysis (EFA) on an associative network, but they do not provide any practical detail in order to produce a similar analysis. Does anyone know how to perform EFA of a network? I guess the input should be the adjacency matrix (which is here called "associative matrix" because the network is a network of associations between words/concepts). It would be great if you could point me to some documentation, maybe in R (that would be perfect!). Thank you very much.


title: Same (poor) results with linear and non-linear classifiers

text: I am doing classification of "text quality" using four classes and using 30 features with 1300 samples. I am using the following classifiers:   1. LDAC based on linear discriminant analysis from mlpy.    2. svm with rbf and gridsearch    3. svm with polynomial kernel degree 5 and gridsearch and    4. random forest.  The svm and random forest classifiers are based on scikit-learn. I am getting similar (poor) accuracy results with all the classifiers. I get approx 57.1% with LDAC, svm and random forests vary from 57.8- 58.1%. The accuracy is based on cross validation and the corresponding confusion matrices (sum of diagonal/total). Why am I getting similar results with all the classifiers? Any suggestions about the reasons for having similar results when using linear and non-linear classfiers? Is this agreement between linear/non-linear classifiers at low- level a strong hint of flaws in some specific area? I appreciate any comments here in order to try to improve this situation!


title: Transformation to increase kurtosis and skewness of normal r.v

text: I'm working on an algorithm that relies on the fact that observations $Y$s are normally distributed, and I would like to test the robustness of the algorithm to this assumption empirically. To do this, I was looking for a sequence of transformations $T_1(), \dots, T_n()$ that would progressively disrupt the normality of $Y$. For example if the $Y$s are normal they have skewness $= 0$ and kurtosis $= 3$, and it would be nice to find a sequence of transformation that progressively increase both. My idea was to simulate some normally approximately distributed data $Y$ and test the algorithm on that. Than test algorithm on each transformed dataset $T_1(Y), \dots, T_n(y)$, to see how much the output is changing. Notice that I don't control the distribution of the simulated $Y$s, so I cannot simulate them using a distribution that generalizes the Normal (such as the Skewed Generalized Error Distribution).


title: How to model my data into a probit

text: I am beginning my thesis, and I need some advice. I am trying to estimate a probit model. The binary dependent variable is employment status, and the independent variables include network size, age, years since migration, and network quality. There are 18,000 migrants in the sample with IDs, and there is a year variable (2001–11). Network size refers to the number of migrants living in the same area as that migrant. I have annual information (2001–11) for network size, age, and years since migration, but only 2011 information for employment status and network quality, so it a very unbalanced panel. Given this, can I estimate this model using an ordinary probit command?               probit empstatus nwsize age years nwquality if migrant == 1      I do not want to use random effects, etc. by using `xtprobit`. Will I have to create 11 network size variables (`nwsize01`, `nwsize02`, ...., `nwsize11`) and regress:               probit empstatus nwsize01 nwsize02, ...., nwsize11 age years nwquality if migrant == 1      Taking away the panel aspect? Might both be possible? Any help would be appreciated.


title: Best way to find non-randomness regions in these or similar count data?

text: Let say I have data in a shape:   [0,0,0,0,1,0,0,0,1,1,0,0,0,0,0,0,0,0,1,0,0,2,0,0,0,....] - so mainly zeros....   However I know how long is my 'signal' and how many counts are they.   * Is it possible to find **regions of non-randomness** in it?  I assume that in my signal counts are spread in both random and non-random fashion, and I would like to find only the last one.   * Can I do it with HMM?


title: Goniometer data analysis

text: I am designing a health science experiment for uni and am focussing the study on measuring range of motion in a person's knee before and after they wear a copper bracelet. To measure the ROM, we are using a goniometer, which is a giant protractor with 2 arms that align with the angle of the joint to provide the degree of movement in that joint. This might sound basic, but I'm having a fight with myself as to whether this data is ordinal or ratio. I think it might be ratio because measurements are taken when the person is lying down (ie flat at 0 degrees), and you cant have a -ve degree, but it could be considered ordinal if you use the baseline measurement as an arbitrary '0 degree' point and then have + or - results from that point. Which way is the 'normal' way of using this data, because ultimately I just want to have 'before' and 'after' data to compare - please help!


title: Bootstrapping estimates of out-of-sample error

text: I know how to use bootstrap re-sampling to find confidence intervals for in- sample error or R2:               # Bootstrap 95% CI for R-Squared     library(boot)     # function to obtain R-Squared from the data      rsq <- function(formula, data, indices) {       d <- data[indices,] # allows boot to select sample        fit <- lm(formula, data=d)       return(summary(fit)$r.square)     }      # bootstrapping with 1000 replications      results <- boot(data=mtcars, statistic=rsq,           R=1000, formula=mpg~wt+disp)          # view results     results      plot(results)          # get 95% confidence interval      boot.ci(results, type="bca")      But what if I want to estimate out-of-sample error (somewhat akin to cross- validation)? Could I fit a model to each boostrap sample, and then use that model to predict for each other bootstrap sample, and then average the RMSE of those predictions?


title: Weighting datapoints by uncertainty - asymmetric CIs

text: I have several estimates obtained using different softwares, these programs give me an estimate as well as a 95% confidence interval. However, these CIs are not symmetric around my point estimates, the upper confidence interval is generally much larger. I want to weight these estimates by their uncertainty, but I am not sure how to do this when dealing with asymmetric CIs.


title: Standard error of regression coefficients without an assumption of homoscedastic normal noise

text: I have a time series that is affected by two (or more) kinds of events. When event $A$ happens, some signal is linearly added to the time series (the signal lasts, for example, for 100 time points). When event $B$ happens, another signal (also of 100 time points length) is linearly added to the time series. The noise is not completely normal - there are some very strong outliers here and there. My goal is to model the two signals. If the responses to the two events didn't overlap in time, I could simply average all the time points delayed t time points after each event kind. The standard deviation of the values in this subset of time points divided by $\sqrt{n}$ (n is the number of event occurrences) performs well enough in estimating the error: if in one of event $A$'s occurrences there was a strong outlier at time $x$, it will increase the local estimate of error for that particular delay. My problem is that the two events do overlap in time. Therefor, I have to use regression in order to deconvolve the two responses. This is done by creating a set of pulse predictors of 100 delays from the occurrences of each event type (a Finite Impulse Response model). By OLS I can get a good estimation of the shapes of the two signals. However, I can't get a good estimate of the error: the standard OLS error estimator is assuming homoscedastic error. Hence, a time point at time $t_1$ that was distorted by a strong outlier has the same error estimation as time point $t_2$, which wasn't noised so much. This makes the error estimator useless in telling a true increase in the signal from a momentary noise. As a crude approximation, I thought about estimating the error for each delay $t$ (and each event kind) by the standard deviation of the particular residuals at that specific delay, divided by $\sqrt{n}$ (n is the number of event occurrences, not total time points as in the usual OLS SE). I'd appreciate any help with this problem. And in particular, answering these two questions:   1. Is the latter solution biased? And if it is, in what way?    2. Is there a better way to estimate the noise in this case? Due to the limitations of my scientific field, it should not involve rich parametric assumptions.


title: Software for making semi variograms and analyses?

text: Our lab used to have a program called GS+ that let us make semi variograms on our data and analyse them. Unfortunately the licence has expired. Is there another piece of software that lets you do this? If so, how specifically? We prefer a Free and open source licence, but all suggestions are welcome. Thank you.


title: How to combine three different scores

text: 1. I have a population of 300 cases. It's split in three sub-populations 50, 50 and 200 in size.   2. I have developed three (different) models resulting in a score variable which rank orders each of the sub-populations separately. The score-ranges for all three segments are different... 40 to 240, 20-410 and 4 to 600. Within each sub-population, higher score means higher level of toxicity i.e. level of risk. Score are not comparable at the moment (that's what I'm trying to achieve).   3. I know (a priori) that the absolute level of toxicity for first sub population is highest, average for middle sub-population and lowest for the third sub-population. What I'm trying to achieve is single rating scale. Ideally, I should be rescale each score variables to reflect my a priori knowledge.... first population should have highest scores, then second and then third... something like this               range A      from 40 to 240     to       600 to 1000     range B      from 20 to 410     to       400 to 700     range C      from 4 to 600      to       0   to 500      NOTE: I intentionally allowed for overlaps Im currently stuck... I know its more art than science, and I will appreciate any ideas you may have. Thanks


title: Deep learning vs. Decision trees and boosting methods

text: I am looking for papers or texts that compare and discuss (either empirically or theoretically):   * Boosting and Decision trees algorithms such as Random Forests or AdaBoost, and GentleBoost applied to decision trees. with   * Deep learning methods such as Restricted Boltzmann Machines, Hierarchical Temporal Memory, Convolutional Neural Networks, etc. More specifically, does anybody know of a text that discusses or compares these two blocks of ML methods in terms of speed, accuracy or convergence? Also, I am looking for texts that explain or summarize the differences (e.g. pros and cons) between the models or methods in the second block. Any pointers or answers addressing such comparisons directly would be greatly appreciated.


title: forecasting: problem with accuracy

text: A few weeks back I used the following command: **accuracy(train,test)** where train and test are training and test data respectively. Last night I updated R and the forecast package and used the same command and I got error. After trying a little I used the following command **accuracy(train,test[1:30])** and it worked. (I was checking the accuracy of 30 forecasted values). Is there some change in the forecast package or did I do something wrong? regards Leo


title: Predictive model for error of another model

text: Does it make any sense to build a regression model for a certain target variable on a certain training set. Then build a regression model for the errors of the previous model ( real values vs predicted). And then add the results of both models on a test set. Has anyone heard of anything like this? I tried to do some prior research but dont know how to even begin searching for this. To explain a little bit more, I'm trying to do an interpolation of maximum tree heights per squared KM over a certain territory. I tried with a random forest regression wich does quite well but I wanted to improve it bit if there is a chance. I have big errors (this is all on a 20% test set) when I try to predict very tall trees or very short trees. The model is biased and doesnt have enough information to "see" these big and small heights. I was looking for a way to correct this if theres is one. I uploaded an image of what is happening, I ordered the error from smallest (overestimation) to largest (subestimation). I hope it helps. ![enter image description here](http://i.stack.imgur.com/Hb7NB.png)


title: Resources for learning R

text: I'm interested in learning R on the cheap. What's the best free resource/book/tutorial for learning R?


title: Best tutorials for R beginner?

text: Heyo- I just started fooling around with R and want to find a good tutorial. Suggestions? Also, are there other things I should be doing to beef up my R chops? Context: very little coding experience, intermediate knowledge of SPSS Much obliged.


title: Where can I find useful R tutorials with various implementations?

text: I'm using **R** and the manuals on the R site are really informative. However, I'd like to see some more examples and implementations with R which can help me develop my knowledge faster. Any suggestions?


title: Regression analysis of different ranges in percent?

text: I have two variables, both percentages but have been converted to decimal numbers, and one variable ranges from roughly 15-30 percent and the other ranges from 55-70 percent. The lower range is my response, dependent variable. The residuals are really high but should have some correlation. Would scaling one of the variables help and what type of scaling? Logarithmic? I hope I explained my question well. I am not an expert statistician but trying to learn. Any help would be appreciated. BTW, I am using R. EDIT Would scaling one of the variables help lower the residuals?


title: Where should N/A be placed and how should it be coded on a 7-point Likert scale?

text: **Background:** I'd like to use a 7-point Likert scale for an upcoming survey (focused primarily on opinions). For the most part, the questions will be centered around the "To what extent...?" stem. That is, the answers will range from, e.g., "To a very small extent" to "To a very great extent". Using the 7-point Likert scale, it is most often recommended to include the "neutral" as the center/middle score. Based on my research thus far, however, I'd like to discourage people from not taking a _stance_ on the subject. That is still debatable though. Now, to ensure I'm in no violation of "knowledge liability" (i.e., respondents must know the answer), it was recommended to include the "N/A" in the question/scale. In my view, "N/A" does _not_ equal the neutral. So, I don't think the center position would be the right location on the spectrum. ### My questions:   * **If I were to include the neutral in the survey questions, what's its location in the scale?**   * **How should N/A answers be numerically coded?** If it would be the most-outer left point (code value = 1), it would take the place of the "To a very small extent"... which then has been shifted by 1 score to the right (code value = 2) and thereby isn't the total opposite of the "To a very great extent" (code value = 7).


title: Pros and cons of clustering algorithms?

text: So I am somewhat new to the realm of applied statistics and machine learning and am currently trying to figure out how to approach a problem I'm working on. Let me first describe the data that I have. I am collecting data from a device during user interaction. There are a discrete number of possible events that can happen on this system, and the data will be labeled with (1) what event occurred and (2) what user was using the device. For each event, I will extract features from data gathered on the device during that event. For instance, let's say that the device is a machine running Linux, the user is Bob, and the event is issuing the 'cd' command via command line. When this even occurs I log features such as CPU usage and system time. So, I have a number of events, each of which are labeled with a user and event descriptor. Each of these events also has N feature values associated with it. I want to cluster per user and event (meaning that I want to run a clustering algorithm for Bob across all the times that he has used the 'cd' command using the extracted feature values). I need the output of this clustering to be a number of descriptors that represent the clusters in N-dimensional space (I was thinking perhaps a centroid and radius per cluster). So, hopefully that was a sufficient explanation of what I'm working with. Given that I'm a bit of a n00b with this stuff, I'd really appreciate any advice with the following questions:   1. What clustering algorithms would be well-suited to this type of data and what are their pros/cons?   2. What are (generally) the variables in clustering algorithms that can be used to fine-tune your model?   3. Are there a number of standardized ways to succinctly describe clusters other than centroid + radius?   4. Are there tools outside of Weka that you all would recommend for running this sort of analysis?   5. Are there any inherent flaws to my approach? Lastly, if any of you know of any good resources that could be considered "crash courses" in clustering algorithms, I'd love to hear about them! Thanks in advance!


title: Derivation of the posterior over topics in LDA

text: When studying the latent Dirichlet allocation, I am not very clear about some procedures in their deriving equations. Please refer to the attached figure, how to understand those two steps, marked as 1 and 2 in the figure. ![enter image description here](http://i.stack.imgur.com/03JPU.jpg)


title: QQ plot does not match histogram

text: I have a histogram, kernel density and a fitted normal distribution of financial log returns, which are transformed into losses (signs are changed), and a normal QQ plot of these data: ![http://tinypic.com/r/34ocwvr/6](http://i.stack.imgur.com/8joMW.png) The QQ plot shows clearly that the tails are not fitted correctly. But if I have a look at the histogram and the fitted normal distribution (blue), even the values around 0.0 are not fitted correctly. So the QQ plot shows that only the tails are not fitted appropriately, but clearly the whole distribution is not fitted correctly. Why does this not show up in the QQ plot?


title: Example of discontinous effect of x on y dataset (for paper)

text: For a paper, I need an example of a dataset $(x_i,y_i)$ where the residuals are $iid$ (the $x$ do not represent time) with a discontinuity on the effect of $x$ on $y$.   * I have already found a dataset in Berger and Pope 2010 but the jump there is not clearly visible   * The 'LifeCycleSavings' dataset has been suggested, but it has only 50 observations (ideally, I would be looking at 100). Does anyone know of other datasets?


title: What skills are required to perform large scale statistical analyses?

text: Many statistical jobs ask for experience with large scale data. What are the sorts of statistical and computational skills that would be need for working with large data sets. For example, how about building regression models given a data set with 10 million samples?


title: Adjusting sample weights in AdaBoost

text: I am trying to read up about AdaBoost from Tibshirani (page 337 onwards), and would appreciate some help in understanding it better. The book says that "For each successive iteration m = 2, 3, . . . , M the observation weights are individually modiﬁed and the classiﬁcation algorithm is reapplied to the weighted observations. At step m, those observations that were misclassiﬁed by the classiﬁer Gm−1 (x) induced at the previous step have their weights increased, whereas the weights are decreased for those that were classiﬁed correctly. Thus as iterations proceed, observations that are diﬃcult to classify correctly receive ever-increasing inﬂuence. Each successive classiﬁer is thereby forced to concentrate on those training observations that are missed by previous ones in the sequence." I am not able to understand what it means to "reapply the algorithm to the weighted observations". Say for example if I am doing 2-class text classification, what am I doing to my observations (document vectors)? How does this "force" the classifier to concentrate on some samples vs others?


title: What is the median of a non-central t distribution?

text: What is the median of the non-central t distribution with non-centrality parameter $\delta \ne 0$? This may be a hopeless question because the CDF appears to be expressed as an infinite sum, and I cannot find any information about the inverse CDF function.


title: Automation of standardization

text: Can anyone point towards applications, SQL algorithms or similar that automate the production of standarized rates/ratios (ie directly/indirectly age/sex etc.-standardized rates/ratios) with confidence intervals from numerator, denominator & reference data? The process would ideally prevent inappropriate calculations, eg from small counts.


title: I'm not sure how to interpret my binary logistic regression output from SPSS

text: My dependant variable is diagnosis of cancer malignant being 0 and benign being 1. And my covariate is mean radius (of the tumour). I get this:               Variables in the Equation:                                B       S.E.   Wald    df   Sig.   Exp(B)  95% C.I.for EXP(B)                                                                       Lower   Upper     Step 1a radius    1.020      .104   96.518  1   .000    2.774   2.263   3.401             Constant -15.039    1.474   104.091 1   .000    .000                  a Variable(s) entered on step 1: radius.                  How do I interpret this in words?


title: Which is the null hypothesis for testing whether I've broken my simulation?

text: The situation: I'm writing agent-based computer simulations in which there are random effects which can be biased by various parameters. I run the simulation with the same parameters many times in order to get a distribution of outcomes. All simulations run in the same underlying simulation framework. Now I've changed some source code in the framework. This code randomizes a certain sequence; the first n elements are then chosen from this sequence. This change should not make any meaningful difference in results, I believe, but I'd like to be sure. While I'm not going to test the change for every simulation I've ever run or ever will run, I can at least check the simulations I'm running now to see whether I get different outcome distributions. I ran the same simulation 300 times using the old code and then the new code. My general question is: Using frequentist statistics, what should my null hypothesis be? The most natural null, I believe, is that the two distributions are the same. However, the dangerous outcome is that I treat the code change as not affecting distributions, when in fact it did affect them. That makes it seem as if the null should be that the distributions are different. (In order to try to get answers that will be especially helpful, maybe I should say up front that I find _some_ of the traditional ways of talking about choice of null in terms of "what you want to reject" as _sometimes_ unhelpful, counterintuitive, or contrary to what scientists do--as opposed to what they say when they're being fastidious. These answers How to specify the null hypothesis in hypothesis testing fit my understanding pretty well, but I would like to some confirmation about how I'm thinking about my particular problem. In general, my feeling is that the real issue in the choice of null concerns risk. The null should be the alternative which is less costly--for you, for society, whatever--to accept. It's the safe choice. That's why we accept the alternative only if the results are extremely unlikely given the null. In my case, the safe choice is to distrust the intuition that the code change made no difference. Feel free to tell me that I'm wrong, of course! My question is also related to Which one is the null hypothesis?, but I didn't find that discussion sufficiently helpful for my case.) * * * Here are further details which don't matter, I think. You can stop reading unless the following looks relevant: I did test the "same-distribution" null for several different pairs of outcomes. I used a bootstrapped Kolmogorov-Smirnov tests to get p-values, because the simulations are too complex to know the nature of the underlying distributions, and because the outcomes are roughly discrete (outcomes can vary from -1 to 1, but always converge to near one of several fractional values in practice). In all cases the p-values were greater than .05, so I do not have good evidence the distributions are the same in each case. If I take the null to be that _all_ of these distributions are the same, then I am doing multiple testing--8 tests in all, as it happens, one for each of the numbers generated by a single simulation run--so the relevant level is much lower, .00625: Even better, I think. None of my p-values are anywhere near that cutoff! Do you see what I mean: That this seems like the wrong way to think about testing whether the code change has resulted in a difference in outcome distributions in this case? It's too easy to find that the distributions are the same, and too hard to find that they're different. Should the null be that the distributions are different?


title: I am stuck with proof of MLE for Logistic regression

text: ![MLE of Logistic regression](http://i.stack.imgur.com/6qhRY.jpg) Hi all, I am doing a self study of logistic regression and working through the proof of the MLE; please see picture from Montgomery and Peck text. Could someone please fill in the intermediate steps that take us from the right hand side of 14.7 to the final statement in the photo. I understand going from L() to L-L() and the change from the product operator to the summation operator.. but I don't understand the algebraic manipulation to get $\ln(\frac{\pi}{1-\pi})$ on the left hand side or where $1-y_i$ went when it was brought down from the exponent position. I hope that makes sense and I appreciate your help.


title: Exponential decay function

text: I have a model where I assumed a variable X would have an exponential decay over time t. I modeled it in R as $$N = N_1e^{kt},$$ where $k < 0$ and $t$ is 1:95. However, looking are some real data that I have captured, it seems that the decay is slightly slower to begin with but is much faster towards the end. The greatest gap in my projected data and the actual data is in the middle. Beginner question, but I was wondering is that type of decay was a known function, or is there some way for me to model it with the exponential decay function. I have a plot of the two functions, but am not allowed to upload it here. The plot shows a symmetry between the two sets of data, as if the actual data plot is a mirror image of the model data. I generate the model data using this               gf <- 1.35 * exp(-0.00946 * 1:95)      and the actual data gf2 is                   c(1.35, 1.33776990615813, 1.32552158334356, 1.31325480792751,      1.30096935107516, 1.28866497856582, 1.27634145060489, 1.26399852162705,      1.25163594009013, 1.23925344825924, 1.22685078198044, 1.21442767044339,      1.20198383593233, 1.18951899356456, 1.17703285101575, 1.1645251082312,      1.15199545712208, 1.13944358124587, 1.12686915546963, 1.11427184561532,      1.10165130808571, 1.08900718946957, 1.07633912612483, 1.06364674373805,      1.05092965685848, 1.03818746840501, 1.02541976914382, 1.01262613713477,      0.999806137143936, 0.986959320019858, 0.974085222030623, 0.961183364158636,      0.948253251349736, 0.935294371712896, 0.922306195666412, 0.909288175026071,      0.896239742030314, 0.883160308296893, 0.870049263704956, 0.856905975195814,      0.843729785484912, 0.830520011676707, 0.817275943773186, 0.803996843065683,      0.790681940398468, 0.777330434291129, 0.763941488905233, 0.75051423183888,      0.737047751730707, 0.723541095652442, 0.70999326626635, 0.696403218720661,      0.682769857252295, 0.669092031461845, 0.655368532220595, 0.64159808716336,      0.627779355713794, 0.613910923580398, 0.599991296651462, 0.586018894205227,      0.571992041337234, 0.557908960489589, 0.543767761945991, 0.529566433130946,      0.51530282652052, 0.500974645933667, 0.486579430925755, 0.472114538946752,      0.457577124852233, 0.44296411726136, 0.428272191135991, 0.413497735800733,      0.398636817423197, 0.383685134710617, 0.368637966230041, 0.353490107290977,      0.338235793692888, 0.322868608762547, 0.30738136887844, 0.291765980931759,      0.276013262639391, 0.260112712872577, 0.24405221348003, 0.227817635241112,      0.211392306411718, 0.19475627881785, 0.177885285875707, 0.160749213513732,      0.143309764509418, 0.125516708879546, 0.10730146998924, 0.0885651903719612,      0.0691537518211836, 0.048795218179725, 0.0268838502712772)


title: Read cointegration relations R

text: I want to figure out how to read the cointegration relation between 5 cointegrated variables. I present the results down here with r=number of cointegrated equations is 3. Variables are normalized respect to A.                   A   B   C    D     A   1   1   1    1     B   5   -3  -5  -16     C   4   20  8   -14     D   -2  6   9   -15      Please can you help by confirming me that this is the right way to write the 3 cointegrated equations: $A=-5B-4C+2D$; $B=A+3*B-20*C-6*D$; $C=A+5*B-8*C-9*D$ with lag length=1 But in this case we have $B$ in both sides and it doesn't help. Please me help me to find how to write the equations.


title: Combining features to improve F1 score in random forest

text: In a project we have to use random forests to classify data into 5 sets. The data set is large (300k data elements and 9 features) and the result is based of F1 score rather than mean square error. A suggested method is to use new features apart from the given 9 features. Can some one please suggest why would taking combinations of features like X1+X2*X3^2 where Xi is a feature improve the accuracy? I am skeptical since all the trees in random forest are already grown to full depth.


title: Split-split-plot design and lme

text: I’m working on a data set in order to evaluate the impact of drying on sediment microbial activities. The objective is to determine if the impact of drying varies among sediment types and/or depth within the sediment. The experimental design is as follows:   The first factor **Sediment** corresponds to three types of sediment (coded _Sed1_ , _Sed2_ , _Sed3_ ).   For each type of Sediment, sampling was carried out on three sites (3 sites for Sed1, 3 sites for Sed2, 3 sites for Sed3). **Site** is coded : _Site1_ , _Site2_ , ..., _Site9_.   The next factor is **Hydrology** : within each site, sampling is carried out in a dry plot and in a wet plot (coded _Dry_ / _Wet_ ).   Within each of the previous plot, sampling is carried out at two **Depths** ( _D1_ , _D2_ ) in triplicate. There are a total of n = 108 samples = 3 Sediment * 3 Sites * 2 Hydrology * 2 Depths * 3 Replicates. I use the lme function in R (lnme package) as follows :               Sediment<-as.factor(rep(c("Sed1","Sed2","Sed3"),each=36))     Site<-as.factor(rep(c("Site1","Site2","Site3","Site4","Site5","Site6","Site7","Site8","Site9"),each=12))     Hydrology<-as.factor(rep(rep(c("Dry","Wet"),each=6),9))     Depth<-as.factor(rep(rep(c("D1","D2"),each=3),18))     Variable<-rnorm(108)          mydata<-data.frame(Sediment,Site,Hydrology,Depth,Variable)          mod1<-lme(Variable~Sediment*Hydrology*Depth, data=mydata, random=~1|Site/Hydrology/Depth)      I found an example of a comparable split-split-plot design and its analysis in : http://www3.imperial.ac.uk/portal/pls/portallive/docs/1/1171923.PDF Could someone confirm that this is the right way to analyze these data?   Do you think that the random structure is correctly specified according to my experimental design? Thank you in advance for your help.


title: Factoring expected prediction error

text: I'm trying to understand the derivation of Expected Prediction Error, as described in The Elements of Statistical Learning. Specifically, it says: > $EPE(f) = E(Y - f(X))^2$ > > By conditioning$^1$ on X, we can write EPE as: > > $EPE(f) = E_XE_{Y|X}([Y - f(X)]^2 | X)$ And, a footnote for this says: > $^1$Conditioning here amounts to factoring the joint density $Pr(X,Y) = > Pr(Y|X)Pr(X)$ where $Pr(Y|X) = Pr(Y,X)/Pr(X)$ and splitting up the bivariate > integral accordingly. I'm confused on the following:   1. Is $E(Y - f(X))^2$ equivalent to $E((Y - f(X))^2)$ and is $E_XE_{Y|X}(g(X, Y))$ equivalent to $E_X (E_{Y|X}(g(X, Y)))$? I'm not used to seeing $E$ used without an argument surrounded by parentheses.   2. I don't understand why $E_{XY}(g(X,Y)) = E_XE_{Y|X}(g(X,Y) | X)$.   3. What is meant by "splitting up the bivariate integral accordingly" in the footnote?   4. Besides Wikipedia, what are some resources I should be looking at so that I don't get confused by these concepts?


title: How does bootstrapping in R actually work?

text: I've been looking into the boot package in R and while I have found a number of good primers on how to use it, I have yet to find anything that describes exactly what is happening "behind the scenes". For instance, in this example, the guide shows how to use standard regression coefficients as a starting point for a bootstrap regression but doesn't explain what the bootstrap procedure is actually doing to derive the bootstrap regression coefficients. It appears there is some sort of iterative process that is happening but I can't seem to figure out exactly what is going on.


title: How to compare two discrete series?

text: Suppose I have two series of numbers and they don't have to be time series. $S_1=$ {9,8,1 ...} $S_2=$ {7,9,2 ...} What are the methods to compare these two series and make a conclusion such that these two series are statistically similar, they have similar patterns or they are completely different? I'm looking for pattern similarity, for example do they both draw an S curve, or they are linear, and I want to ignore time shifts.


title: Maximum likelihood optimization error in R

text: **Edit: A link to the data and nlm() tracing information has been added, and the code has been changed in accordance with suggestions from the comments** Download the data from here. (Dropbox link). Download the nlm() tracing information from here. (Dropbox link) I am working on a problem to find the maximum likelihood estimators $\left( \hat\sigma_n^2,\hat l,\hat\sigma_f^2 \right)$ for some data. The log likelihood is $$\log \mathcal{L} \left( \sigma_n^2,l,\sigma_f^2 \right)= -\frac{1}{2} y^T \left(\mathbf{K}+\sigma_n^2\mathbf{I} \right)^{-1}y - \frac{1}{2} \log \left| \mathbf{K}+\sigma_n^2\mathbf{I} \right| - \frac{n}{2} \log2\pi$$ $y$ is a $n \times 1$ vector obtained from my data (response variables). $\mathbf{I}$ is the identity matrix. $\mathbf{K}$ is a $n\times n$ matrix depending on $\sigma_f^2$ and $l$; its matrix elements are: $$K_{pq} = \sigma_f^2 \exp\left( - \frac{1}{2l^2} \left| x_p -x_q \right|^2 \right)$$ $x_i$ are $3\times1$ vectors (independent variables). I tried to find the maximum likelihood estimators in R by minimizing the negative log likelihood. The functions I tried using were _nlm_ and _optim_. Both gave the same error - NA/Inf were produced. How do I proceed to calculate the estimators? I would prefer more accessible answers. (I am only a beginner) I will write my code below for reference. The matrix $\mathbf{A}$ in the code (and download) is given by: $$A_{pq} = \exp \left( - \frac{1}{2} \left| x_p -x_q \right|^2 \right) $$               load("A.Rdata") # Contains matrix A to help calculate K (in the formula above).                     # A is the exponential part of the formula. (but without the 'l squared' part)      load("y.Rdata") # Contains the response variable y          num_unique <- 786          Calculate_K_plus <- function(vect){       sn2 <- vect[1]*vect[1]       exponent <- 1/(vect[2]*vect[2])       sf2 <- vect[3]*vect[3]       B <- A^exponent       B <- sf2 * B       B <- B + sn2*diag(num_unique)       return(B)}          minus_log_likelihood <- function(vect){      K_plus <- Calculate_K_plus(vect)      K_plus_inv <- solve(K_plus)      out = 0.5 * ( t(y) %*% K_plus_inv %*% y) + 0.5 * log(det(K_plus)) + (num_unique/2)*log(2*pi)      return(out)}          nlm(minus_log_likelihood,c(1,1,1))


title: DLM with autocorrelated residuals

text: I've created a DLM (using the dlm package in R) with a linear trend (dlmModPoly(2)) and a seasonal component (dlmModSeas(seasons)). After filtering, the residuals look Gaussian according to a qqplot. The confidence band around the estimates looks very good as well. However, the acf plot of the residuals shows definite autocorrelation. I assume I need to improve my variance estimates from dlmMLE? Does anyone have any suggestions on where to start?


title: bootstrapping the density: bias from parametric bootstrap of polynomial fit of binned frequencies

text: Bootstrap examples for density estimation usually start with resampling from the empirical distribution of the "micro data", e.g. here on Normal Deviate talking about confidence intervals for nonparametric estimates of the density at some point. But if one is interested in a property of the density (e.g. excess mass at/around a prescpecified point) is it wrong to calculate frequencies by bins, fit some flexible parametric form, and bootstrap with resampling the error, where the error is the discrepancy between the frequency in the bin vs the predicted density? This seems to be what Chetty et al. describe on page 23 (p. 25 of the PDF) of an influential paper with their posted code used in other papers since. Does this "macro level", aggregate (block?) bootstrap do something different from resampling the empirical CDF? In any case, is the error in the parametric fit of the density the only noise we care about in the problem if we think that individual values (not frequencies) themselves have noise in them? Larry Wasserman made a brief comment on this on his blog that > That seems unnecessary and kind of complicated. And it would add additional > bias. The regular bootstrap as I described should work fine. How bad is the additional bias? What is its source exactly? Note the related question about whether the estimand is regular enough to bootstrap or subsample in the first place, and with what rate of convergence.


title: Linear regression simulation

text: I would really appreciate if someone could help me with this question. I want to simulate $n$ datasets in R with eight predictors where $β=(3,1.5, 0, 0, 2, 0, 0, 0)$ and the pairwise correlation between $x_i$ and $x_j$ is set to be $\text{corr}(i,~j) = 0.5$ for $i-j$ or for all $i=j$.


title: Are two Pearson correlation coefficients different?

text: I am aware of this question here, which mine was listed as a duplicate of, but it does not fully answer my question. It did however help me progress a little further so thanks, I couldn't find it before. The online calculator on the above answer also disagrees (though marginally) with the vassarstats.net calculator so I think that further backs up my reasoning for not using a black box. So I'll re-explain my problem and hope that this gets opened to answers: I have two Pearson correlation coefficients which I would like to compare. Each comes from 2 sets of 40 genetically identical lines and are correlations between male and female trait values. I have managed to do $z$-transformations on my correlation coefficients in `R` using the function `atanh()` and replicated that with a home-made function `RtoZ <- function (r) 0.5*log((1+r)/(1-r))`. The problem is what to do next: how do I actually test, in `R`, whether the two correlations are different?


title: Estimation for covariance in the context of mixed model

text: There are 4 subjects and 4 measured data of each subject and the data should be interpreted in the context of a mixed model. What is a reasonable estimation for the covariance of y?               y <- data.frame(response = c(10,19,27,28,9,13,25,29,4,10,20,18,5,6,12,17),                     treatment = factor(rep(1:4, 4)),                     subject = factor(rep(1:4, each = 4)))


title: What is the frequency of a time series for hourly data?

text: I am using R for time-series analysis and predictions, the package 'forecast' to be more precise. I am in a dilemma. I have hourly data that needs a prediction and needs to be analysed. I am using the STLF function, since I set the frequency to 24 (because it's greater than 13). But, when I make the forecasts for the next 6 hours, with a data set containing 300 points, I get the following forecast:               Point          Forecast Lo 80    Hi 80    Lo 95    Hi 95     13.50000       29.60251 21.28421 37.92081 16.88077 42.32425     13.54167       27.84124 18.89136 36.79111 14.15358 41.52889     13.58333       30.89487 21.33420 40.45554 16.27309 45.51665     13.62500       36.04991 25.89498 46.20484 20.51928 51.58053     13.66667       40.40386 29.66798 51.13975 23.98474 56.82298     13.70833       41.13250 29.82644 52.43856 23.84138 58.42362      As you can see, the next points are 13.5, 13.54, ... etc. This is like this probably because the data set is containing 300 points, and 300/24 = 12.5, 301/24 = 12.54167, ... etc, and assuming that the first point is considered to be 1 and not 0, so there you have the way the points are provided. My question is: will I get better results if I adjust the seasonality in a way that will give me forecasts for each hour? I.e. if the next point is 12, then 13, then 14, ... etc, up until 23 and then to start from 0 (24 hours span). If yes, please tell me how to adjust the seasonality to my data? Is there a way for making even more complicated seasonalities? (say, if the data is taken every 5 minutes or so). Thank you in advance for your answer.


title: Multinomial regression after PCA

text: I have extracted 7 factors from 25 variables using principal component analysis with varimax rotation. Now I want to use these factors as independent variables and study their relation with an ordinal dependent variable. For this I applied multinomial logistic regression in SPSS. But it gives me a warning that a number of cells have zero frequencies. What should be done?


title: Determining true mean from noisy observations

text: I have a large set of data points of the form (mean, stdev). I wish to reduce this to a single (better) mean, and a (hopefully) smaller standard deviation. Clearly I could simply compute $\frac{\sum data_{mean}}{N}$, however this does not take in to account the fact that some of the data points are significantly more accurate than others. To put it simply, I wish to preform a weighted average of these data points, but do not know what the weighting function should be in terms of the standard deviation.


title: Can you combined two sources with difference variance to reduce error?

text: I have two samples of data each estimates of a position `x, y` with Gaussian noise. One source has a larger variance than the other. Is this source in any way useful in providing a better estimate of the true position than the source with lower variance or should it just be discounted?


title: Estimate normal distribution mean given standard deviation and samples

text: I have samples from normal distributions with the same mean but different standard deviations. I know the corresponding standard deviation for each sample. Is there a way to more accurately estimate mean other than taking an average over all the samples?


title: How to describe/explain the shape of a distribution which has two peaks?

text: How can I describe the shape of the following distribution, which has two peaks? What are the important things in describing the distribution shape? It is the output of the following `R` code: plot(density(Boston$tax)) ![enter image description here](http://i.stack.imgur.com/MWayt.png) (The data set is in the `MASS` package) I would like to know what rug (these red points) are showing in the previous PDF? ![enter image description here](http://i.stack.imgur.com/L05Zy.png) I created it by the following code on `Boston` data from `MASS` package:               plot(density(Boston$tax))         rug(Boston$tax, col=2, lwd=3.5)


title: Odd shift of statistical significance when some variables are dropped from model

text: I am testing the statistical significance of some explanatory variables in a factor model, using the **Fama-MacBeth two-pass regression** (also known as the Cross-Sectional Approach). Assuming my dependent variables are stock returns, and that my explanatory variables are macro factors, such as Inflation, Oil prices, Industrial Production:   * I first run a rolling-time-series regression (60 months window) of each stock with respect to all my explanatory variables. This will result in a time-series of sensitivities of each stock to the macro variables.   * Then, I run a cross-sectional regression for each month of my sample period, using stock returns as dependent variables, and the estimated sensitivities as independent variable. When the cross-sections are ran, we obtain the time-series factor values. For each factor, I test whether significantly different from zero. Consider the following cross-sectional regression: $R_i,t = \alpha_i,t + \beta_{i,1}(t)F_1,t + \beta_{i,2}(t)F_{2,t} + ... + \beta_{i,n}(t)F_{n,t} + \epsilon_i,t $ where the betas are estimated in previous periods, using time-series regression. The purpose of the cross-sectional regressions is to get time- series values of $F_1$, ..., $F_n$. After I test for significance, I only collect the _k_ significant variables, and include them in a regression. I redo the same test as above (time-series and cross-sectional regressions), only using my selected _k_ variables. BUT NOW, none of my _k_ variables are significant!! I am not sure how to construct my final model based on these results. Any help would be appreciated! Thank you!


title: How to choose a data subset in RapidMiner?

text: I'm working with a CSV which contains approximately 220,000 entries. My aim is to predict one of the attributes (ATT1) using the other 3 (ATT2, ATT3, ATT4). I've been able to do this using NaiveBayes, but now I feel unsatisfied with the result. The reason is that ATT1 can be one of 6 values (VAL1-6), but these are not evenly distributed into the dataset. I'm afraid this could lead to an unprecise prediction. How do I select a given number of entries for each value of ATT1 from within RapidMiner?


title: How do I determine confidence intervals around weighted, probilized events?

text: How would I go about finding out the confidence intervals around a set of distinct bianary occurrences where each occurence has a different associated probability and each occurence is weighted? To be more specific, we typically determine a milestone budget by assigning reasonable probabilities to each milestone. Each milestone typically usually has a differnt probability of success and a different dollar value. For example if we need to determine the amount of money we need to hold out for 3 milestones, we would do this: $$P(\text{milestone}_1) \text{cost}(\text{milestone}_1)+P(\text{milestone}_2) \text{cost}(\text{milestone}_2)+\dots $$ So it could look something like this: Budget = (10% x $\$$1,000,000) + (50% x $\$$5,000,000) + (20% x $\$$20,000,000) = $\$$6,600,000. Each milestone is bianary in that it either happens or it doesn't. so for the first milestone, we are saying that there is a 10% chance that we will need to pay $\$$1M dollars and a 90% chance that we will pay nothing. Obviously, we won't ever actually be paying the probilized amount of $\$$6.6M. How would I go about determining the probability range around where the spend is most likely to fall. In other words, if I wanted to demonstrate where the actual spend might fall with a certainty of 80% how would I go about calculating this? Note that we have more events in my real world problem... Thanks in advance for any help. I am not a statistician at all so please forgive my ignorance and any lack of clarity.


title: Origin of strange formula for equilibrium standard deviation

text: In the paper > M. Avellaneda and J. H. Lee, Statistical arbitrage in the U.S. equities > market, July 2008, in the Appendix on page 46, how does he get equilibrium standard deviation as following: $$\sigma_{eq} = \sqrt{\frac{\text{Variance}(\zeta)}{1 − b^2}}$$ If anyone knows the paper, please explain. Much appreciated,


title: How to Calculate a Z-Score from Power Log Distributions

text: _Please forgive me if the question is non-sensical as I only have a elementary knowledge of statistics, and I have no idea where to start_. I want to perform an analysis of some web analytics information, and determine which of our partners we should be building stronger relationships with because they are sending us highly engaged users. ![enter image description here](http://i.stack.imgur.com/cSO2A.png) My first thought was to calculate the z-score for each of the main metrics in the data set:   * PageViews   * Visits   * Visitors   * Repeat Visits   * etc.  And then for partner sites that have a positive Z-Score for:   * Repeat Visits   * Pages / Visit   * Avg. Time on Page  with a negative Z-Score for Bounce Rates, those were the sites that would be more desirable to improve our relationships with. Additionally, we would focus on the partnersites that met that criteria and **also had a negative Z-Score for PageViews or Visits.** It would be my assumption that **these sites would be the ones that are driving highly engaged users but at a very low volume**. However, **I'm not sure if this assumption is correct**. Given the Z-Score is based on Standard Deviation, which from what I understand is really meant for normal distributions, it doesn't seem like the right path given PageViews, Visits, Visitors, and sometimes Pages/Visits fits more of the Power Law distribution than a normal distribution.


title: Tests for quality of permutation (Y-randomization)

text: Can you please suggest the metrics/tests for accessing the quality of Y-randomization (permutation)? I have some doubts in quality of Y-randomization that I am using with **sample.int** from R. For 100 Y-randomization runs I found that in average 2% of the instances are falling into themselves, and around 5% of instances are falling into the similar instances with correlation higher then "95% quantile of correlation matrix for initial data". The dataset is rather small ~80 instances and have some similar entries with correlation median for all against all around 0.55. Update: part of a code I am using for modeling (RSNSS)               # t_comb_orig - original data, input only          # original data and it's correlation     t_comb_orig <- rbind(t_train_in, t_test_in)     cor_orig <- matrix(NA,nrow(t_comb_orig),nrow(t_comb_orig))     for(l in 1:nrow(t_comb_orig)) {       for(k in 1:nrow(t_comb_orig)) {         cor_orig[l,k] <- cor(t_comb_orig[l,],t_comb_orig[k,])       }     }     cor_orig[lower.tri(cor_orig,diag=T)] <- NA     cor_orig_median <- median(c(cor_orig),na.rm=T)     cor_orig_095qutl <- quantile(c(cor_orig),0.95,na.rm=T)     # modeling     for (i in 1:N) {         # randomize         t_comb_rand <- t_comb_orig[sample.int(nrow(t_comb_orig)),]         t_train_inR <- t_comb_rand[1:nrow(t_train),]          t_test_inR <- t_comb_rand[(nrow(t_train)+1):nrow(t_comb_rand),]         # train model         t_model <- mlp(x=t_train_inR, y=t_train_out,size=c(floor(ncol(t_train_inR)/2)), learnFunc="Rprop", inputsTest = t_test_inR, targetsTest = t_test_out, maxit=50)         t_ext_prediction <- predict(t_model,t_ext_in)         write.table(ifelse(t_ext_prediction >= 0.5, 1, 0), file=paste("test__",i,".txt",sep=""),sep="", col.names = F, row.names = F)         # plot error         png(filename=paste("test__",i,".png",sep=""))         plotIterativeError(t_model)         dev.off()             # randomization test (correlation between scrambled vs real)             for (j in 1:nrow(t_comb_orig)) scrambling_test[j,i] <- cor(t_comb_orig[j,],t_comb_rand[j,])     }     #     write.table(scrambling_test,file="scrambling_vs_real_cor.txt",sep="\t", col.names = F, row.names = F)     scrambling_test_per_run <- matrix(NA,ncol(scrambling_test),1)     for (o in 1:ncol(scrambling_test)) scrambling_test_per_run[o,1] <- length(which(scrambling_test[,o] > cor_orig_095qutl))/nrow(scrambling_test)     write.table(scrambling_test_per_run,file="scrambling_vs_real_cor_per_run.txt",sep="\t", col.names = F, row.names = F)     #


title: Comparing ln(y) =bX vs. OLS and other count data models

text: My data is a balanced panel: 12 years $\cdot$ 600 ids $\cdot$ 12 months = 86,400 observations. For about 150 of my 600 ids, at a certain point in this 12 year period a treatment is applied switching them to the "treated" or "on" state. A dummy variable AFTER indicates when this is true. My outcome is a count, positive only and quite skewed, looks like a power law. I'm interested in estimating the impact of treatment on outcome. I've tried these 3 specifications in STATA and I dont seem to understand why they don't accord well with each other.   1. OLS $\rightarrow$ `reg y AFTER i.id i.year, robust cluster(id)`   2. Log $\rightarrow$ `reg ln(y) AFTER i.id i.year, robust cluster(id)`   3. Poisson $\rightarrow$ `xtpoisson y AFTER i.year, fe vce(robust)`   4. Negative Binomial $\rightarrow$ `nbreg y AFTER i.id i.year, vce(robust)` Here i.id and i.year are fixed effects for id and year respectively. Estimated coeffs with standard errors in parentheses are as follows:   1. OLS: `403922 (223946) *`   2. Log : `-0.664 (0.0252) ***`   3. Poisson : `-.0188 (0.019)` Could someone interpret these for me and tell me why one model would give me a poistive coeff while other give me a negative one? My theory tells me it should be positive, but im not sure why this should be the case that the non linear models give me a different answer. Sum stats for my outcome variable are: range: `[0,3.025e+08]`   unique values: `59,191`   mean: `1.3e+06`   std. dev: `1.1e+07` percentiles: `10% (0), 25% (13480), 50% (74690), 75% (322982), 90% (1.3e+06)`


title: Details on specifying random parameters in 3-level HLM in R

text: I am used to fitting 3-level HLM with HLM software, but now I am moving to `nlme` or `lme4` in `R`. Unfortunately, though, I am struggling a little bit with the proper specification/notation of the random parameters of the 2nd and the 3rd levels. What I am trying to do is to estimate a 3-level linear model where some of the random effects are in the 2nd level, some are in the 3rd, but I do not know how to differentiate them in the "random" part of `lme()`. To make it clear, let me give you the specification equations:               Y = b0 + b1*X1 + b2*X2 + e              b0 = a00 + a01*X3 + r0         b1 = a10 + a11*X4 + r1         b2 = a20 + r2                   a00 = c00              a01 = c10 + u1              a10 = c20 + u2              a11 = c30              a20 = c40      and the final model is               Y = c00 + c20*X1 + c40*X2 + c10*X3 + c30*X1*X4 + r0 + r1*X1 + r2*X2 + u2*X1 + u1*X3 + e      To estimate this model, so far I could only think of the following working codes:               lme(Y ~ X1 + X2 + X3 + X1:X4, random= ~ X1 + X2 + X3 | level3 / level2, data=mydataset)     lmer(Y ~ X1 + X2 + X3 + X1:X4 + (1|level2) + (X1| level3/level2) + (X2|level2) + (1|level3) + (X3|level3), data=mydataset)      The `lme` seems to be somewhat wrong, as it estimates the random parameters of X1, X2 and X3 both in level 2 and in level 3, not the random compenents of X1 and x2 in level 2 and of X1 and X3 in level 3. In `lmer` I tried to do so, but I don't know if it is right to so in this way. Please, could any of you give me some basic light here? Thanks a ton!


title: The value of adding the ROC graph if the AUC is given

text: I always see in papers that when they want to show how they classifiers performed, they provide ROC graph and the AUC score. Now as far as I know only the AUC tells how well the classifier performed, so what is the advantage of adding the ROC graph? What can one tell from it?


title: Orthogonalization for dimensionality reduction

text: There are some methods like singular value decomposition (SVD), principal component analysis (PCA), factorial analysis and many more that are used to reduce a high-dimensional dataset into fewer dimensions while retaining important information. What would happen if for instance the problem was not about finding those very informative features, but synthetizing the data in a way to assure that correlation between the resulting dimensions is minimized regardless of whether the resulting number of dimensions is low enough and capable of retaining most of the variability. With traditional dimensionality reduction techniques such as SVD, PCA or FA you don't really assure that the resulting dimensions will be uncorrelated, I wonder if there is any solution to this problem.


title: Should I use ANOVA to analyse this experiment?

text: We performed an experiment on bacterial growth whereby optical density was recorded for three flasks of _E. coli_ bacteria every 30 minutes for 4 hours. 2 flasks upon the 2 hour mark were inoculated with two different antibiotics, whilst the other served as a control. Would a one-way ANOVA be appropriate in this case? Should I, for example, select two time points: one before adding the antibiotic, one after and do two separate ANOVAs? I will be using class data so I will have 5 replicates for each time point.


title: Interpreting multiple regression coefficients with 2 continuous variables interacting and 2 categorical variables interacting

text: I find it challenging to interpet interaction effects in OLS multiple regression where the interactions are between two categorical variables and between two continuous variables. Say for instance that the interaction between sexfemale: medmobility is -10.1, where sexfemale is a categorical variable with male as the reference group and medmobility is a categorical variable where low mobility is the reference group. What is the meaning of that coefficient, -10.1? Likewise, what if you have an interaction between two continuous variables, say for instance between weight and IQ (weight:IQ) If the coefficient is, for example, 1.3 what does that mean? How does the interpretation change if both have been centered? What about the main effects of these variables? Thanks for the helping me to understand this!


title: AcF and Stationarity

text: Very often in time series literature, it is remarked that if a series is non- stationary the AcF will decrease to zero very slowly while the opposite occurs for a stationary series. What's the basis for this "rule of thumb"? I know that for a strictly stationary process the autocorrelation is independent of time, whereas for a wide-sense stationary process the autocorrelation is a function of the time lag but these don't explain the "rule of thumb".


title: Multinomial logistic predictor matrix

text: I'm reading Tuerlinkx & Wang, "Models for Polytomous Data". They write: "Apart from the intercept, the predictors that are included in the predictor matrix X can be classified into seven groups: item predictors, person predictors, logit predictors, and the four combinations of them. ... A logit predictor has the same value in all rows of the predictor matrix referring to the same component of the vector-valued link function. Interactions between a logit predictor and any of the other predictors are also possible." I don't understand what a logit predictor is. If it has the same value in all rows, then it sounds like it's just the intercept term. But from the first sentence it sounds like there's also another intercept term? And how can you have an interaction between an intercept and another predictor? It'd be helpful to see a model matrix, but table 3.1 doesn't answer my questions.


title: How to use Gaussian Mixture Models for clustering new data?

text: I applied Gaussian Mixture Model on my data and train the model in MATLAB. How I can test my model or use it to cluster new data? Thanks for any answer or comment.


title: High kurtosis and bad skewness

text: Is it necessary to have data that looks normal if you want to apply a dynamic correlation coeffient (dcc)? Should the explanatory variables in the dcc method also be normalized? My dataset has a high kurtosis and bad skewness. I used the log10 to try to make it more normal, but remarkably the kurtosis is much higher after taking the log... Very strange. What other options do I have to make my data more normal and how do I apply that in Stata?


title: Generating non-zero binomial probabilities (n,k ) with small p and large n - k

text: I am trying to generate binomial probabilities (in R) as follows: ${N \choose{k}} p^{k} (1-p)^{(n-k)}$ My problem is given $p \approx 0.03$, and $N =400$, $k>270$, I get the probability equal to $0$. I am actually evauating all the $k$'s up to $n$, so a high $n$ returns a probability matrix that has $0$'s for $k >270$. I was wondering if there is a function in R / or an approximation (binomial or normal) that you would recommend in this case? ## Edit I ve tried dbinom function in R as well - `dbinom(400,280,0.03) = 0`. So that does not help.. What I am after is a non-zero approximation..


title: Likelihood of events coinciding at random

text: I have what is probably a very naive question, but I've been unable to find a suitable explanation elsewhere. I am trying to calculate the likelihood of two errors occurring at the same position if they are randomly distributed amongst all positions available. The problem is difficult to explain in context, so I will simplify the scenario. Say I have two positions that can have a number of identities across a number of events. So I might represent this with two vectors:               X<-(1,1,0,0,0,0,0,0,0,0,1)     Y<-(1,1,0,0,0,0,0,0,0,0,0)      My assumption is that the 1's are errors and the zeros are not. What I want to calculate is the probability of the errors coinciding at corresponding positions in X and Y if the errors are distributed at random, and then how unlikely it is that I observe a certain number of corresponding errors (at the same position in X and Y) in observed data? To do this I have calculated the frequency of an error at each position (so in the e.g. above, 3/11 in X and 2/11 in Y) and then calculated the expected chance of overlap by multiplying these together (3/11*2/11=6/121). I can then calculate the probability of seeing at least two overlaping 1's (from observed data) using the binomial cumulatively -> 1-pbinom(1,11,6/121) = which is the chance of two or more observed overlaps in 11 trials given and expectation of 6/121 per trials. However, does this fully capture the probability of observing this particular scenario (two or more matching errors) if the errors are randomly distributed, or am I missing some information in the fact that the 0's also match up frequently too? In the scenario I am working with, having zeros match up is also good evidence that the positions of 0's and 1s are not random, but as they are often way more frequent than 1's I don't know if I should be considering these matches. Sorry if this question is unclear, but any help would be much appreciated!


title: Hidden Markov model library

text: I am looking for a free hidden markov model libraries in C/C++. I can google and get a list, but a it is impossible to tell which free ones are industry standard and well established. Can anyone provide guidance on this. Maybe the right way to put it would be, what hidden markov model library would you use on the speech recognition system you used to lock your house with a passphrase (ya, I expect anyone on this forum knows not to do that) -- assume that you are impatient and don't like waiting around for your house to decide it is you, i.e. performance counts for a lot.


title: Repeated measures anova or two-way anova?

text: I have carried out an experiment that measured soil chemistry parameters at different times (0, 3, 6 and 9 months) from 12 plots in the field (3 treatments arranged in a fully randomized design with four replications). In each case, soil was randomly sampled in ten squares to give a composite soil sample per plot. Therefore, I believe that each soil is different. I statistically analysed with two-way ANOVA (with treatment and sampling date as factors). However, a colleague told me that since sampling dates are not independent from each other, it doesn’t make sense to use a two-way ANOVA with treatment and sampling date as orthogonal factors. Instead, a repeated measure ANOVA should be used. Which is correct?


title: Poisson processes

text: I have two realizations of a poisson stochastic process, they are over the same space with rate $\lambda_{1}$ and $\lambda_{2}$. What is the probability that N elements in both sequences are the same, e.g., $X(t) \in [0,1000]$? Note: The two processes are independent, i.e. X(s) is independent of X(t). Lets say the number of elements within a process is finite, e.g. K.


title: Global optimization of deep learning models

text: Apparently, deep neural networks have been making an impact recently. The layer-by-layer training of these networks has made it feasible to construct complex, deep, and well-performing neural networks. Still, I feel that some applications of deep learning models might benefit from global optimization approaches that do not easily get stuck in local optima. I haven't seen any research in this direction, though. Couldn't evolutionary/biologically inspired algorithms (e.g., Particle Swarm Optimization, Differential Evolution) be used to make deep learning models more powerful? Or is the computing power necessary for this particular combination of techniques currently a limiting factor?


title: Evaluating mean-squared error

text: Hello I am running a Regression Tree experiment. I am new to Regression Trees, and I am using Mean Squared error to test my tree. I am confused because I am getting a large Mean Squared Error but I am not sure how to evaluate if it is too high. Should different successful trials expect varying MSE? Do there exist methods for checking this? Thanks!


title: How to compute efficiency?

text: Suppose instead of maximizing likelihood I maximize some other function g. Like likelihood, this function decomposes over x's (ie, g({x1,x2})=g({x1})g({x2}), and "maximum-g" estimator is consistent. How do I compute asymptotic variance of this estimator? **Update 8/24/10** : Percy Liang goes through derivation of asymptotic variance in a similar setting in An asymptotic analysis of generative, discriminative, and pseudolikelihood estimators. **Update 9/14/10** : Most useful theorem seems to be from Van der Vaart's "Asymptotic Statistics" Under some general regularity conditions, distribution of this estimator approaches normal centered around its expected value $\theta_0$ with covariance matrix $$\frac{\ddot{g}(\theta_0)^{-1} E[\dot{g}\dot{g}^T] \ddot{g}(\theta_0)^{-1}}{n}$$ Where $\ddot{g}$ is a matrix of second derivatives, $\dot{g}$ is gradient, $n$ is number of samples


title: Is there an accepted definition for the median of a sample on the plane, or higher ordered spaces?

text: If so, what? If not, why not? For a sample on the line, the median minimizes the total absolute deviation. It would seem natural to extend the definition to R2, etc., but I've never seen it. But then, I've been out in left field for a long time.


title: Calculating Kurtosis with Beta distribution in R

text: I am working on an assignment question in R. The problem I am having is to calculate the Kurtosis using Beta distribution. So far I am able to get the kurtosis value of `alpha = 3 beta = 2` using this code:               > kurtosis=(3*alpha^3*beta+3*alpha*beta^3+6*alpha^2*beta^2+alpha^3+beta^3+13*alpha^2*beta+13*alpha*beta^2+alpha^2+beta^2+14*alpha*beta)/((alpha*beta*(alpha+beta+2)*(alpha+beta+3)))      > kurtosis [1] 2.892857 The question asks me to generate 1 000 000 independent identically distributed random variables following a Beta distribution. That is easy:               n=100 (I am using 100 for now as 1 000 000 pretty much just hangs my machine up for hours)     p<-rbeta(n, shape1=alpha, shape2=beta)     p               [1] 0.6264632 0.8556988 0.8947960 0.2508657 0.2478296 0.3652293 0.4058185 0.9133942 0.2521895 0.3846150 0.7311879 0.6123567 0.5687362 0.5272337 0.7321268 0.5167320 0.2628915      [18] 0.6827633 0.4034805 0.7767011 0.9246785 0.5946597 0.7077071 0.4014387 0.5738073 0.7225824 0.2736515 0.2674327 0.3449036 0.8527233 0.9142692 0.8225190 0.8081674 0.7760414      [35] 0.5947241 0.8211220 0.6825275 0.1994827 0.3839823 0.8674370 0.4383073 0.6895321 0.6948742 0.4175201 0.4870160 0.4839224 0.5868314 0.6400326 0.5956450 0.7255520 0.4290051      [52] 0.8997959 0.8925137 0.7729661 0.6168284 0.5989350 0.6112734 0.1849604 0.5710712 0.5314022 0.8364838 0.7888875 0.8285151 0.5491049 0.9450382 0.6302862 0.6035714 0.4513216      [69] 0.1588688 0.6744496 0.6010260 0.6908575 0.5219169 0.6897020 0.1486987 0.6108051 0.5032918 0.8729911 0.8216785 0.4243647 0.6615128 0.5479400 0.6326745 0.7449108 0.7917081      [86] 0.3582359 0.7631012 0.6198699 0.4402614 0.5099704 0.3586977 0.4702083 0.5465371 0.7308831 0.3773918 0.5209375 0.2836313 0.6724058 0.7975775 0.6507358      Now I am using this formula to calculate Kurtosis:               (sum(Xi-Xbar)^4/n)/(sum(Xi-Xbar)^2/n)^2      Putting this into R               for (i in 1:n) {     value[i] <- ((p[i]-mean(p))^4/n)/(((p[i]-mean(p))^2/n)^2)}      When I call value I get this return:                value                [1] 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100      [44] 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100      [87] 100 100 100 100 100 100 100 100 100 100 100 100 100 100      I tried to do the first value of (p[1]-mean(p))^4/n)/(((p[1]-mean(p))^2/n)^2) in mathematica and i get value of 7.389056099 I found that the problem is when I try to divide numerator by denominator, is there another way to do division other than / ? Or am I approaching the problem wrong?


title: Estimation of variance

text: In a simple linear regression, the estimation of ![enter image description here](http://i.stack.imgur.com/ocOzM.gif) why is it n-2 and not n-1? thanks


title: Combined individual error

text: This is from my college project management course. Reading through an example question here it says: > When estimating in parts, the total error will be less than the sum of the > part errors. This makes sense to me. > For a 1000-hour job with estimating accuracy of ± 50%, the estimate range is > from 500 to 1500 hours. This also makes sense > If the estimate is independently made in 25 parts, each with 50% error, the > total would be 1000 hours, as before and the estimate range would be from > 900 to 1100 hours Okay, How did we get that? > To combine independently-made estimates > >   * Add the estimated values. >   * Combine the variances (squares) of the errors. > Cool, following you... > With 25 estimates for a 1000-hour job > >   * Each estimate averages 40 hours >   * The standard deviation is 50%, or 20 hours > Here's where it stops making sense. It goes on to take the square root of sum of the variances which makes sense, but where did he get the standard deviation being 50%? I could be missing something here?