title: Which is the correct mean?

text: My colleague and I disagree on the following data, and the calculation of the mean: Given $n$ respondents, each respondent has three chances to score a hoop at basketball. A respondent has 'hooped' successfully if they score 1 or more hoops. Given the following table of scores               respondent id | try 1 | try 2 | try 3     1                1        1       1         2                0        0       1     3                0        0       0     4                1        1       0      This becomes:               respondent id | success     1                1         2                1     3                0     4                1       I argue 3/4 successes, he argues 3/(4 * 3) successes justifying this by saying they each have 3 attempts. For me this would only make sense if we were not reducing, but counting each success as well. Am I wrong? Have I missed the hoop?


title: Out of sample bootstrapping and significance

text: I am investigating the relationship between temperature fields obtained from numerical weather models and electricity demand. I am applying a PCA-based approach, i.e. I study the linear relationship between main temperature patterns/modes and main demand patterns/modes. Given that I am working on summer yearly demand, I have time-series with few samples (<20) and for this reason I've decided to apply the following bootstrapping procedure:   1. I create a temperature and elec. demand datasets with the usual sampling with replacement   2. I create my linear model between the two fields   3. I create a temperature dataset with the not-selected samples and I project them on the PCA-space I've just computed    4. I calculate out-of-sample output I do this for about 5000 times and in the end I obtain a matrix with only the out-of-sample outputs. I calculate the mean on all the out-of-sample predicted demands for each year and I use it to calculate RMSE error. I think this approach could be considered a .632 bootstrap procedure. I'd like to compute the significance of the obtained results. I was thinking about the possibility to shuffle at each bootstap iteration the temperature dataset in order to see whether I obtain similar results breaking the direct temporal link between demand and temperature. Given that I don't have a robust statistical background I'd like your opinion about any method to obtain the statistical significance of my bootstrap procedure.


title: Doane's formula for histogram binning

text: I'm implementing various algorithms to estimated the best number of bins to use for histograms. Most of the ones I am implementing are described on the wikipedia page: https://en.wikipedia.org/w/index.php?title=Histogram&oldid=548769683#Number_of_bins_and_width * I am stuck on a problem with Doane's formula:               1 + log(n) + log(1 + kurtosis(data) * sqrt(n / 6.))      where `n` is the data size. The problem is when the kurtosis is negative and `n >> 1` beacuse the argument of the `log` becomes negative. * (that page has changed since this was posted, link edited to point to the page as it was at the time of posting)


title: EM algorithm manually implemented?

text: I want to implement the EM algorithm manually and then compare it to the results of the `normalmixEM` of `mixtools` package. Of course, I would be happy, if they both lead to the same results. The main reference is Geoffrey McLachlan Finite Mixture Models 2000. I have a mixture density of two Gaussians, in general form, the log-likelihood is given by (McLachlan page 48): ![loglikelihood](http://i.stack.imgur.com/PsPRz.png) The $z_{ij}$ are 1, if the observation was from the $i$th component density, otherwise 0. The $f_i$ is the density of the normal distribution. The $\pi$ is the mixture proportion, so $\pi_1$ is the probability, that an observation is from the first Gaussian distribution and $\pi_2$ is the probability, that an observation is from the second Gaussian distribution. The E step is now, calculation of the cond. expectation: ![condexp](http://i.stack.imgur.com/J6puy.png) which leads, after a few derivations to the result (page 49): ![tau](http://i.stack.imgur.com/RGYw1.png) in case of two Gaussians (page 82): ![taugauss](http://i.stack.imgur.com/HCneA.png) The M step is now the maximization of Q (page 49): ![q](http://i.stack.imgur.com/4RG1v.png) This leads to (in case of two gaussians) (page 82): ![mugauss](http://i.stack.imgur.com/lK4Ly.png) ![sigmagauss](http://i.stack.imgur.com/PwJTn.png) and we know that (p. 50) ![pigauss](http://i.stack.imgur.com/u4G56.png) We repeat the E, M steps until ![diffchange](http://i.stack.imgur.com/K9yrV.png) is small I tried to write a R code (data can be found here)               # EM algorithm manually     # dat is the data               # initial values     pi1<-0.5     pi2<-0.5     mu1<--0.01     mu2<-0.01     sigma1<-0.01     sigma2<-0.02     loglik[1]<-0     loglik[2]<-sum(pi1*(log(pi1)+log(dnorm(dat,mu1,sigma1))))+sum(pi2*(log(pi2)+log(dnorm(dat,mu2,sigma2))))                    tau1<-0     tau2<-0     k<-1          # loop     while(abs(loglik[k+1]-loglik[k]) >= 0.00001) {          # E step     tau1<-pi1*dnorm(dat,mean=mu1,sd=sigma1)/(pi1*dnorm(x,mean=mu1,sd=sigma1)+pi2*dnorm(dat,mean=mu2,sd=sigma2))     tau2<-pi2*dnorm(dat,mean=mu2,sd=sigma2)/(pi1*dnorm(x,mean=mu1,sd=sigma1)+pi2*dnorm(dat,mean=mu2,sd=sigma2))          # M step     pi1<-sum(tau1)/length(dat)     pi2<-sum(tau2)/length(dat)          mu1<-sum(tau1*x)/sum(tau1)     mu2<-sum(tau2*x)/sum(tau2)          sigma1<-sum(tau1*(x-mu1)^2)/sum(tau1)     sigma2<-sum(tau2*(x-mu2)^2)/sum(tau2)          loglik[k]<-sum(tau1*(log(pi1)+log(dnorm(x,mu1,sigma1))))+sum(tau2*(log(pi2)+log(dnorm(x,mu2,sigma2))))     k<-k+1     }               # compare     library(mixtools)     gm<-normalmixEM(x,k=2,lambda=c(0.5,0.5),mu=c(-0.01,0.01),sigma=c(0.01,0.02))     gm$lambda         gm$mu     gm$sigma          gm$loglik      The algorithm is not working, since some observations have the likelihood of zero and the log of this is -Inf. Where is my mistake? Help would be very appreciated, since this topic is pretty advanced for me.


title: What is a complete list of the usual assumptions for linear regression?

text: What are the usual assumptions for linear regression? Do they include:   1. a linear relationship between the independent and dependent variable   2. independent errors   3. normal distribution of errors   4. homoscedasticity Are there any others?


title: Forecasting and auto-correlation

text: I'm reading this chapter forecasting principles and practise from a forecasting book. The author has explained a linear regression model. Now this linear regression model will definitely have some errors in fit. For these errors, the author has listen certain characteristics, which should be taken into account. He assumes that these errors:   1. Have mean zero; otherwise the forecasts will be systematically biased.   2. Are not autocorrelated; otherwise the forecasts will be inefficient as there is more information to be exploited in the data.   3. Are unrelated to the predictor variable; otherwise there would be more information that should be included in the systematic part of the model. The problem is that I'm not able to comprehend these and come up with an example to get a hang of these points. Could somebody point out some examples / situations from daily life so as to understand these points. For example the author wrote in point 2.] ..."there is more information to be exploited in the data"... What is this information? How does it come into being? If I don't even know how this information comes into existence, then how will I be able to understand its existence or effects for that matter? Is there a simple explanation with example for my question? I have been trying to understand this since quite some time, please help.


title: When doing ANOVA, what do we need to assume is normally distributed?

text: Do we assume that the population distributions are normal, or that the sample distributions are normal, or that the sampling distribution is normal? If the latter, what do we mean by sampling distribution in this context? Why do we need to assume this? Does the normality of these things entail or suggest the normality of others?


title: Validating data set using standard deviation

text: Suppose I have two data sets. The average of one data set is 5.5 secs and the standard deviation is 0.0435. The average of the other data set is 5.2 secs but the standard deviation is 0.5123. I know that something is wrong with the second data set because the standard deviation is too high. My question is how I should interpret the standard deviation value so that I could use the value to determine if data set is valid. Standard devation might not be the one I should use to do this. I will appreciate suggestions for any other methods to validate this data set.


title: Determining similar users from hierarchical clustering

text: I use hierarchical clustering to cluster users which are similar to each other based on a Jaccard coefficient. I have now coded a solution to extract similar users based on hierarchical clustering: To group users that are similar I compute the Jaccard coefficient between all the different users. I then check each user for similarity based on their Jaccard score. So if user1 and user2 have a similarity score of .2 then these users are very similar. I perform same for all users building a similarity table of users. So once grouping is finished I have a table structure which shows which users are close to each other based on the Jaccard coefficient. Is this the correct approach to extract similar users based on hierarchical clustering?


title: Multiple regression analysis

text: I have a data set based on 14 field sites, the dependent variable I am investigating is count data recorded at each site (on 3-5 visits). The independent variables (I have around 45) are fixed at each site (apart from wind speed and temperature). An example data set, illustrating this is below:                           Trees    water   butterflycount      Site 1      6        8       3          Site 1      6        8       12          Site 2      3        3       8          Site 2      3        3       0      What is the best way of preparing this kind of data for a multiple regression analysis in R? Would data1.glm = glm(butterflycount ~ 1, data = data1, family = poisson) and the working through adding the lowest AIC valued variables to the equation work? Any help would be greatly appreciated.


title: How to use coxph for prediction of event?

text: I fitted a cox model as: `cfit <- coxph(Surv(tr$start, tr$end, tr$event1) ~ imp+sec + cluster(tr$compid))` Now i want to get the probability of event1 for the test dataset for each start and stop time interval, given it has other variables present. How can i use the predict function for that?


title: Inferring transition matrices in continuous time Markov processes

text: Suppose I have a process **X1** , for which I do not have a generator matrix, only a transition (probability) matrix **P1** for some time interval **T** , e.g. **T** =100. Suppose I have another process **X2** , such that **X2** is identical to **X1** except that if **X2** transitions to a state below a certain threshold, it is returned to one of the states above the threshold (with equal probability of each) at 1 time period thereafter, where 1 is the finest granularity time interval under consideration. My question is, would there be any way to estimate a new transition matrix **P2** for **X2** (over time period **T** ) directly from the old one, without the need to worry about the so-called embedding problem for Markov chains (or solving a quadratic optimisation to infer the best generator)? To illustrate, if there are 3 states a-c above the threshold and another 3 states d-f below, then my known 1-period transition probabilities would be:                     a   b    c   d    e    f         a[ .    .    .    .  .  .,         b  .    .    .    .  .  .,         c  .    .    .    .  .  .,         d  1/3  1/3  1/3  0  0  0,         e  1/3  1/3  1/3  0  0  0,         f  0     0    0   0  0  1];      where ' **f** ' is an absorbing state. How would I use this information to obtain **P2** for the **T** -period case, given that I know **P1**?


title: Normalization vs. scaling

text: What is the difference between data 'Normalization' and data 'Scaling'? Till now I thought both terms refers to same process but now I realize there is something more that I don't know/understand. Also if there is a difference between Normalization and Scaling, when should we use Normalization but not Scaling and vice versa? Please elaborate with some example.


title: Linear regression with repeated measures in R

text: I was unable to figure out how to perform linear regression in R in for a repeated measure design. In a previous question (still unanswered) it was suggested to me to not use `lm` but rather to use mixed models. I used `lm` in the following way:               lm.velocity_vs_Velocity_response <- lm(Velocity_response~Velocity*Subject, data=mydata)      (more details on the dataset can be found at the link above) However I was not able to find on the internet any example with R code showing how to perform a linear regression analysis. What I want is on one hand a plot of the data with the line fitting the data, and on the other hand the $R^2$ value along with the p-value for the test of significance for the model. Is there anyone who can provide some suggestions? Any R code example could be of great help. * * * **Edit**   According to the suggestion I received so far, the solution to my analyze my data in order to understand if there is a linear relation between the two variables Velocity_response (deriving from the questionnaire) and Velocity (deriving from the performance) should be this:               library(nlme)     summary(lme(Velocity_response ~ Velocity*Subject, data=scrd, random= ~1|Subject))      The result of summary gives this:                   > summary(lme(Velocity_response ~ Velocity*Subject, data=scrd, random= ~1|Subject))         Linear mixed-effects model fit by REML          Data: scrd                 AIC      BIC   logLik           104.2542 126.1603 -30.1271              Random effects:          Formula: ~1 | Subject                 (Intercept) Residual         StdDev:    2.833804 2.125353          Fixed effects: Velocity_response ~ Velocity * Subject                                    Value Std.Error DF    t-value p-value     (Intercept)               -26.99558  25.82249 20 -1.0454288  0.3083     Velocity                   24.52675  19.28159 20  1.2720292  0.2180     SubjectSubject10           21.69377  27.18904  0  0.7978865     NaN     SubjectSubject11           11.31468  33.51749  0  0.3375754     NaN     SubjectSubject13           52.45966  53.96342  0  0.9721337     NaN     SubjectSubject2           -14.90571  34.16940  0 -0.4362299     NaN     SubjectSubject3            26.65853  29.41574  0  0.9062674     NaN     SubjectSubject6            37.28252  50.06033  0  0.7447517     NaN     SubjectSubject7            12.66581  26.58159  0  0.4764880     NaN     SubjectSubject8            14.28029  31.88142  0  0.4479188     NaN     SubjectSubject9             5.65504  34.54357  0  0.1637076     NaN     Velocity:SubjectSubject10 -11.89464  21.07070 20 -0.5645111  0.5787     Velocity:SubjectSubject11  -5.22544  27.68192 20 -0.1887672  0.8522     Velocity:SubjectSubject13 -41.06777  44.43318 20 -0.9242591  0.3664     Velocity:SubjectSubject2   11.53397  25.41780 20  0.4537754  0.6549     Velocity:SubjectSubject3  -19.47392  23.26966 20 -0.8368804  0.4125     Velocity:SubjectSubject6  -29.60138  41.47500 20 -0.7137162  0.4836     Velocity:SubjectSubject7   -6.85539  19.92271 20 -0.3440992  0.7344     Velocity:SubjectSubject8  -12.51390  22.54724 20 -0.5550080  0.5850     Velocity:SubjectSubject9   -2.22888  27.49938 20 -0.0810519  0.9362      Correlation:                                (Intr) Velcty SbjS10 SbjS11 SbjS13 SbjcS2 SbjcS3 SbjcS6 SbjcS7 SbjcS8 SbjcS9 V:SS10 V:SS11 V:SS13 Vl:SS2 Vl:SS3     Velocity                  -0.993                                                                                SubjectSubject10          -0.950  0.943                                                                                SubjectSubject11          -0.770  0.765  0.732                                                                                SubjectSubject13          -0.479  0.475  0.454  0.369                                                                                SubjectSubject2           -0.756  0.751  0.718  0.582  0.362                                                                                SubjectSubject3           -0.878  0.872  0.834  0.676  0.420  0.663                                                                           SubjectSubject6           -0.516  0.512  0.490  0.397  0.247  0.390  0.453                                                                    SubjectSubject7           -0.971  0.965  0.923  0.748  0.465  0.734  0.853  0.501                                                             SubjectSubject8           -0.810  0.804  0.769  0.624  0.388  0.612  0.711  0.418  0.787                                                      SubjectSubject9           -0.748  0.742  0.710  0.576  0.358  0.565  0.656  0.386  0.726  0.605                                               Velocity:SubjectSubject10  0.909 -0.915 -0.981 -0.700 -0.435 -0.687 -0.798 -0.469 -0.883 -0.736 -0.679                                        Velocity:SubjectSubject11  0.692 -0.697 -0.657 -0.986 -0.331 -0.523 -0.607 -0.357 -0.672 -0.560 -0.517  0.637                                 Velocity:SubjectSubject13  0.431 -0.434 -0.409 -0.332 -0.996 -0.326 -0.378 -0.222 -0.419 -0.349 -0.322  0.397  0.302                          Velocity:SubjectSubject2   0.753 -0.759 -0.715 -0.580 -0.360 -0.992 -0.661 -0.389 -0.732 -0.610 -0.563  0.694  0.528  0.329                   Velocity:SubjectSubject3   0.823 -0.829 -0.782 -0.634 -0.394 -0.622 -0.984 -0.424 -0.799 -0.667 -0.615  0.758  0.577  0.360  0.629            Velocity:SubjectSubject6   0.462 -0.465 -0.438 -0.356 -0.221 -0.349 -0.405 -0.995 -0.449 -0.374 -0.345  0.425  0.324  0.202  0.353  0.385     Velocity:SubjectSubject7   0.961 -0.968 -0.913 -0.740 -0.460 -0.726 -0.844 -0.496 -0.986 -0.778 -0.718  0.886  0.674  0.420  0.734  0.802     Velocity:SubjectSubject8   0.849 -0.855 -0.807 -0.654 -0.406 -0.642 -0.746 -0.438 -0.825 -0.988 -0.635  0.783  0.596  0.371  0.649  0.709     Velocity:SubjectSubject9   0.696 -0.701 -0.661 -0.536 -0.333 -0.526 -0.611 -0.359 -0.676 -0.564 -0.990  0.642  0.488  0.304  0.532  0.581                               Vl:SS6 Vl:SS7 Vl:SS8     Velocity                                           SubjectSubject10                                   SubjectSubject11                                   SubjectSubject13                                   SubjectSubject2                                    SubjectSubject3                                    SubjectSubject6                                    SubjectSubject7                                    SubjectSubject8                                    SubjectSubject9                                    Velocity:SubjectSubject10                          Velocity:SubjectSubject11                          Velocity:SubjectSubject13                          Velocity:SubjectSubject2                           Velocity:SubjectSubject3                           Velocity:SubjectSubject6                           Velocity:SubjectSubject7   0.450                   Velocity:SubjectSubject8   0.398  0.828            Velocity:SubjectSubject9   0.326  0.679  0.600          Standardized Within-Group Residuals:             Min          Q1         Med          Q3         Max      -1.47194581 -0.46509026 -0.05537193  0.39069634  1.89436646           Number of Observations: 40     Number of Groups: 10      Warning message:     In pt(q, df, lower.tail, log.p) : NaNs produced     >       Now, I do not understand where I can get the R^2 and the corresponding p-values indicating me wether there is a linear relationship between the two variables or not, nor I have understood how my data can be plotted with the line fitting the regression. Can anyone be so kind to enlighten me? I really need your help guys...


title: What to do when every SVM returns "not a member"?

text: We have trained about 200 Linear Support Vector Machines (hand-coded in C#) belonging to our 200 categories and we use them in text categorization. Due to the shortness of the text in our training samples and other uncontrollable factors, very often every machine returns a negative number for (classifies as negative) a given document. Returning "None" as an answer is not feasible because it happens maybe about 40% of the time, and not classifying 40% of the documents is unacceptable from a business point of view. What we do is, we don't bother about the sign of the results and return the category with the maximum value anyway. However, I feel this makes our error rate larger than need to be. So here is my question: Is there a principled middle-of-the-road way? That is it will sometimes (but not too often, maybe about 5% of the time) return "None" and it will sometimes return the maximum of the negative numbers, so that |misclassified documents| + k.|unclassified documents| will be minimum (for some k<=1)? For example, is there a way to translate the distance given by an SVM to a probability value? Or, relatedly, is there a principled way to compare the output of two Support Vector Machines, i.e. +0.2 denotes a higher probability of membership in one of them? Any suggestion, source, paper is welcome. Thank you.


title: Moderated moderation - conditional effects suggest 3-way interaction but no main effect

text: I am doing a modearated moderation using the process add on. There is no significant 3-way interaction for my 3 moderator (one is gender) but the conditional effects suggests that the signficant 2-way interaction is different for males and females. How do I interpret/report this?


title: Mean accuracy for sum of squares/variance

text: Formality says that when computing the sample variance or sum of squares one computes a sample mean with the sample data given. If for some reason I had access to a better estimate of the mean (from a bigger sample) then can I use that mean instead of the mean from the given sample? I want to compute the $R^2$ and t-statistic from a regression analysis, for which I need $\text{SSX}$ and $\text{SSY}$ as the sum of squares from data set Y and data set X. This SSX/Y are basically the variance(X)*(n-1)=SSX and variance(Y)*(n-1)=SSY. So if for my regression analysis I use a smaller sample than the one I actually have. After this I need to compute $R^2$ for example, in which I wish to sue a better estimate of the mean/variance using a bigger sample. This seems consistent with the conceptual theory with me, as a better estimate is always desired, and considering that both data sets (smaller and extended sample) are from the same population. Is there a reason not to do this?


title: Statistical method to compare time-series with different periods

text: I have 2 study areas Study area 1: 12 meteorological stations-years available 1981-2000(same data and step) Study area 2: 10 meteorological stations-years available 1985-2011(same data and step) From the 1st study area I have to use the above mentioned data series. I would like to know which statistical test I can use in order to choose the appropriate years from the 2nd study area. Thank you in advance.


title: Canonical correlation analysis with rank correlation

text: Canonical correlation analysis (CCA) aims to maximize the usual Pearson product-moment correlation (i.e. linear correlation coefficient) of the linear combinations of the two data sets. Now, consider the fact that this correlation coefficient only measures linear associations - this is the very reason why we also use, for example, Spearman-$\rho$ or Kendall-$\tau$ (rank)correlation coefficients which measure arbitrary monotone (not necessarily linear) connection between variables. Hence, I was thinking of the following: one limitation of CCA is that it only tries to capture linear association between the formed linear combinations due to its objective function. Wouldn't it be possible to extend CCA in some sense by maximizing, say, Spearman-$\rho$ instead of Pearson-$r$? Would such procedure lead to anything statistically interpretable and meaningful? (Does it make sense - for example - to perform CCA on ranks...?) I am wondering if it would help when we are dealing with non-normal data...


title: Detecting patterns of cheating and collusion in competition

text: I'm looking for ideas on how to discover cheating in sports that have judges awarding points. I can't name the sport since I participate in it at the moment. It's been known for a long time that this sport has rampant collusion behind the scenes. In each competition there are competitors who are evaluated on skill by a panel of judges. The spring is simple since the judges only have to place the first six competitors in rank order of preference. The thing is that sometimes the judges are actually teaching the competitor. The better the competitor is ranked in the competition, the better the judge looks and hence gets more business, and since judges know each other they have the opportunity to collude. Therein lies the incentive to cheat. Any suggestions in how I would go about detecting if judges are fairly scoring competitors? The main constraint is that due to the size of the competitions there is not really a lot of data. Judges and competitiors would not really meet very often, four to five times at the most.


title: Pairwise comparison of vectors with unequal sizes and unequal variances

text: I have 3 vectors which are unequal in size (lengths: 21, 33 and 7). Each vector contains ratio of bug-to-duplicate bugs in 3 different time periods. I have to find whether there is a difference between the 3 groups. I cannot use one-way ANOVA since the groups have unequal sample sizes and unequal variances. I cannot use TukeyHSD because of unequal variances. Is Dunnett's modified Tukey Kramer test the only test available for such a kind of data? If, so how to interpret the results of that test since the Dunnett's test in R does not provide a p-value? Is chi-squared test applicable to this data? Thanks.


title: Rpart and tree models -- cross validation beyond cp value

text: I normally cross-validate my tree-models (rpart) only on cp-value, e.g. by using caret or xpred or the internal rpart cross validation function. However, there are other rpart parameters -- for example minimum size and surrogate usage. Should one cross-validate over these parameters as well, or is this superfluous (and hence not included in the standard cross-validation packages?


title: What is standing in the way of pulsed neural networks being used in applications?

text: Pulsed or Spiking neural networks incorporate more of the membrane dynamics of biological neurons, where pulses carry information to the next layer. The neurons do not necessarily have to "fire" all at the same time, as they would in a backprop, for example. Yet, there seem to be barriers against employing these models for machine learning problems. What specific issues stand in the way of machine learning practitioners using models that are more biologically realistic?


title: Time Series Regression: serially correlated errors vs autocorrelation of residuals?

text: If I have a time series regression model with serially correlated errors: $Y_t=\beta_0+\beta_1X_t+\epsilon_t$ $\epsilon_t=\rho\epsilon_{t-1}+v_t$ why is the estimate of the correlation ($\rho$) in the second equation NOT the autocorrelation of the residuals from the first equation? The autocorrelation of the residuals ( _from the first equation using OLS regression_ ) is easy to get, BUT the rho is estimated with more difficulty [e.g. nonlinear least squares]. But I'm not following why.


title: Asymptotic Efficiency of Two-Stage Least Squares

text: Apparently Wooldridge, Introductory Econometrics, 2002ed is the only book showing that two-stage least squares (2SLS) is asymptotically efficient. I cannot get a copy of the proof. Is it correct to just use the generalized method of moments (GMM) proof substituting in the instrumental variables?


title: proof FGLS asymptotically efficient

text: Prove that FGLS is asymptotically efficient. Does one have to use Cramer Rao to do this?


title: Marginal of a uniform distribution

text: Given $f\left(x|\theta\right)=1/\theta, 0\leq x\leq \theta,L\left(\theta, a\right)=\left(a-\theta\right)^2,$ and $\pi\left(\theta\right)=\theta e^{-\theta},\theta\gt 0$ I've seen Problem calculating joint and marginal distribution of two uniform distributions, but it's not what I'm looking for. Problem: calculate the marginal density of $x$, i.e., $f\left(x\right)$. \begin{eqnarray} f\left(x\right) &=& \int f\left(x|\theta\right)\pi\left(\theta\right)d\theta\\\ &=& \int \frac{1}{\theta}\cdot\theta\cdot e^{-\theta}d\theta\\\ &=& \int e^{-\theta}d\theta. \end{eqnarray} Somehow the last equation should turn out to be $e^{-x}$, according to the book. I don't see how this is true, since I'm integrating over $\theta$. I don't understand why $x$ should end up in the solution since it's nowhere in the original equation. I'm also not sure what the limits of integration are. If $\theta > 0$, then I though maybe they should be $\left[0, \infty\right)$, but that gives $1$ as the solution to the integration.


title: Anova from R output interpretation

text: I have a question on how a statistician would normally interpret an anova output. Say I have anova output from R.               > summary(fitted_data)          Call:     lm(formula = V1 ~ V2)          Residuals:          Min       1Q   Median       3Q      Max      -2.74004 -0.33827  0.04062  0.44064  1.22737           Coefficients:                 Estimate Std. Error t value Pr(>|t|)         (Intercept)  2.11405    0.32089   6.588  1.3e-09 ***     V2           0.03883    0.01277   3.040  0.00292 **      ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 0.6231 on 118 degrees of freedom     Multiple R-squared: 0.07262,    Adjusted R-squared: 0.06476      F-statistic:  9.24 on 1 and 118 DF,  p-value: 0.002917           > anova(fit)     Analysis of Variance Table          Response: V1                Df Sum Sq Mean Sq F value   Pr(>F)        V2          1  3.588  3.5878  9.2402 0.002917 **     Residuals 118 45.818  0.3883                         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1       From the above, i guess the most important value is Pr(>F), right? So this Pr, is less than 0.05 (95% level). How should my "explain" this? Do I explain it in "association", ie, V2 and V1 are associated (or not) ? or in terms of "significance"? I always felt that I couldn't understand when people say "This value is significant....". So what is "significant"? Is there a more intuitive form of explanation? like "I am 95% confident that ...." . Also, is the Pr value the only important piece of information? or can i also look at residuals and the rest of the output to "explain" the result? thanks


title: Visualising generalized linear model

text: What kind of plot are normally used in Generalized linear model and what are their interpretations? Especially for Standardized deviance residual vs fitted value plot, what can we see from the plot?


title: Comparing different ways of specifying the same response variable

text: I have a question about different ways of specifying a response variable in a model and what effect this could have on my results. In my example, I wanted to model the number of times a bird visits a nest to feed its young during a set period of time. So the data would be a set of counts, let's call it `Number of Visits`. Now one thing that might affect how many visits an adult makes is the size of the brood (`Brood Size`). I have seen two ways of modeling this in the literature:   1. Have Number of Visits as the response variable and include Brood Size as an explanatory variable.   2. Have Number of Visits/ Brood Size as the response variable and include other explanatory variables as appropriate. How would these approaches differ and is either superior to the other, would each approach give a different answer? I can see that using each approach you might be modeling subtly different things. In option `1)` you seem to be modeling the overall number of visits while controlling for brood size, whereas in option `2)` you seem to be measuring the per capita number of visits.


title: Fast alternatives to the EM algorithm

text: Are there any speedy alternatives to the EM algorithm for learning models with latent variables (especially pLSA)? I'm okay with sacrificing precision in favor of speed.


title: Finding weights for variables in kNN

text: I'm using euclidean distance for kNN. I have labeled data, I have took logarithm of some variables to make them look more like normaly distributed and scaled them all. And now I would like to multiply some variables by weights, then compute euclidean distance and train kNN. But how to find those weights ? My idea is to determine centers of classes this going to be set C, and then make optimization of kNN on set C by random search, I think that I can't do it on subset of training set, because it size would by to high or too small for accurate representation/sampling of dataset Do you have any other ideas ? I don't think that changing parameters k and l going to have the same approach as mine or mayby does it ?


title: What is the use of the line produced by qqline() in R?

text: The `qqnorm()` R function produces a normal QQ-plot and `qqline()` adds a line which passes through the first and third quartiles. What is the origin of this line? Is it helpful to check normality? This is not the classical line (the diagonal $y=x$ possibly after linear scaling). Here is an example. First I compare the empirical distribution function with the theoretical distribution function of ${\cal N}(\hat\mu,\hat\sigma^2)$: ![comparison of cumulative distribution functions](http://i.stack.imgur.com/ny3n2.jpg) Now I plot the qq-plot with the line $y=\hat\mu + \hat\sigma x$; this graph roughly corresponds to a (non- linear) scaling of the previous graph: ![qqnorm along with the "good" line](http://i.stack.imgur.com/ALKb2.jpg) But here is the qq-plot with the R qqline: ![qqnorm and qqline](http://i.stack.imgur.com/0ZwkI.jpg) This last graph does not show the departure as in the first graph.


title: Is it okay to convert PCA scores to absolute value?

text: I am calculating PCA using a numpy/python approach. The principal components are in the form of GIS grids (rasters), created in a manner very similar to the GRASS GIS approach described here: http://grass.fbk.eu/gdp/html_grass62/r.covar.html Every cell in an output component grid contains a PC score. Within each component grid there is a range of negative to positive PC scores. Is it okay to convert the scores to absolute value? Is the direction of the value important, or just the relative within-component (i.e., within-grid) magnitude? Thanks.


title: When to use residual plots?

text: I have performed a simple regression analysis between one dependent variable (DV) and one explanatory variable (IV). If the p-value from the regression analysis for the IV is not significant, should I still use residual plots to verify that the regression model used was correct (and the statement of non-significant IV is correct)? Or should one use residual plots only for models that include significant IVs?


title: How to best analyse relationship between individual-level predictor and organisation level outcomes?

text: I have inherited a dataset with two types of data from an employee opinion survey project spanning multiple organizations. We have responses by individual employees on questions to do with their work situation (such as "I have the resources I need to do my job"). Then we also have metrics data from each organization, such as number of days of sickness absence. However, I only have an average value from each organization for the metrics data. So I have "average number of days of sickness absence" from each organization. Here's an example of what the data looks like: the first column has the organisation ID, the second column an employee's response to a survey question, and the third column the average number of sickness absence days for the organization               Organisation#    Employee survey response    Average #days absence for org.     1                4                                  4.2     1                2                                  4.2     1                5                                  4.2     2                1                                  2.8     2                4                                  2.8     2                3                                  2.8     3                1                                  1.7     3                5                                  1.7     3                4                                  1.7      What is the best way to analyse the relationship between employee survey responses and organizational level metric data? In this case, we are not in the first instance interested in how the relationship varies by organization - the focus is on the nature of the relationship between the survey questions and the organization level metrics. There is clearly a mix of levels going on here, so I sense that a multi-level approach may be necessary, and a straight correlation between the two variables isn't appropriate. However, I don't have much experience with multi- level designs and analyses, so would be very grateful for any suggestions and input.


title: Reference On Panel Data Analysis

text: I'm looking for good panel data analysis references. I've Baltagi & "Panel Data Econometrics in R". But I want something which is short (max 100 pgs), lucid (know this criteria differs) & have examples to work through. I'll be grateful if somebody can give me some references.


title: When can RMSEA be zero?

text: If $\text{RMSEA} = 0$, it means $\chi^2 < df$. Does it disqualify RMSEA as a criterion to evaluate the model fit, or is it just the explanation why it is zero?


title: What are subjective interestingness measures?

text: I came across the word _subjective interestingness_ measures in my book, where author says that: > Subjective interestingness measures are based on user belief in the data. > These measures find patterns interesting if they are unexpected > (contradicting user's belief) or offer strategic information on which user > can act. (Data Mining and Concepts and Techniques, by Jiawei Han and Micheline Kamber) I very much confused how can "unexpected" data be an interesting pattern? Am I understanding the concept wrongly? Could you provide me with an illustration?


title: Multiple regression and hypothesis test $H_0$:$\beta_2=0$

text: Multiple regression model   $H_0$:$\beta_2=0$, $H_1$:$\beta_2 \neq 0$   where $\beta_2$ is the vector of elements ($\beta_2, \beta_3, \dots, \beta_k$) and $\beta$ is slope of regression line. Why it is equivalent to a test based on the statistic $$\frac{R^2/(k-1)}{(1-R^2)/(T-k)}$$ where $R^2$ is the square of the multiple correlation coefficient of the equation. I don't know how to solve this. Please give me any suggestion or hint.


title: Neighborhood analysis in Matlab using a dot plot

text: I have points in a 2D graph (coordinates: X,Y property: Z). I would like to find for every point the closest, for example, 5 points and save their properties. What would be the easiest approach? Update:   Using the following code:               %Synthetic data     A = {[1,1]; 'A'};     B = {[2,2]; 'B'};     C = {[3,3]; 'C'};          %plot     D = {A{1};B{1};C{1}};     VarPlot = cell2mat(D);     plot (VarPlot,'.');          %knnsearch     [IDX,dist] = knnsearch(VarPlot(:,1),VarPlot(:,2))      I receive the following result:               IDX =           1      2      3          dist =           0      0      0      What does this mean? And how can I link the result to the properties A,B and C. I am new to this kind of questions.


title: What distribution is most commonly used to model server response time?

text: I have a servlet-based application wherein I measure the time taken to complete each request to that servlet. I already compute simple statistics like the mean and maximum; I'd like to produce some more sophisticated analysis however, and to do so I believe I need to properly model these response times. Surely, I say, response times follow some well-known distribution, and there are good reasons to believe that distribution is the right model. However, I don't know what this distribution ought to be. Log-normal and Gamma come to mind, and you can make either one sort of fit real response time data. Does anyone have a view on what distribution the response times ought to follow?


title: Issues with the Huynh-Feldt values using either Anova or ezANOVA in R?

text: Is there an issue with the output of the Huynh-Feldt correction when running an ANOVA for both between and within factors (mixed design) in R? My issue is that when I run a 'within-subject (only)' ANOVA using either the 'Anova' or 'ezANOVA' functions, the output for the Huynh-Feldt correction is the exact same I get when using SPSS. However when I run an ANOVA on a mixed design, and consider both between and within subject effects using either the 'Anova' or 'ezANOVA' functions, my Greenhouse-Geisser corrections match up to SPSS, but the Huynh-Feldt corrections are slightly off compared to the output from SPSS? Any reason why this could be? I appreciate any advice/help anyone could provide!


I need to calculate the Probabilistic Sensitivity Analysis for a function. I was given this: beta distributions assigned to represent uncertainty   And have this parameter with this data:               variable = d_progress     probability of variable : 0.1     n = 100     r = cases = 3      Now, I need to use this function or something?               p <- rbeta(n, shape1=alpha, shape2=beta)      I already have the parameters so I can find `p`. If I have:               p<-rbeta(100, 1, 99)       What I do is I create a vector now? Of 100 values? So, I need to create a loop, 100 times for the following formula:               for each of 100 p generated (       EV = p*b_par       I then submit each EV in a vector or something...       I need to Find the different values of EV for different p,         so I just find the mean of all the 100 EV in the end     )


title: Checking a panel unit root test in R done manually

text: I have spent much time looking for a special package that could run the Pesaran(2007) unit root test (which assumes cross-sectional dependence unlike most others) and I have found none. So, I decided to do it manually; however, I don't know where I'm going wrong, because my results are very different from Microsoft Excel's results (in which it is done very easily). My data frame is made up of 22 countries with 506 observations of daily price indices. Following is the model to run using the Pesaran(2007) unit root test: (i) With an intercept only $$\Delta Y_{i,t} = a_i + b_iY_{i,t-1} + c_i\overline{Y}_{t-1} + d_i\Delta\overline{Y}_{t-1}+ e_i\Delta\overline{Y}_{t-2}+ f_i\Delta\overline{Y}_{i,t-1}+ g_i\Delta\overline{Y}_{i,t-2} + \varepsilon_{i,t}$$ where $\overline{Y}$ is the cross-section average of the observations across countries at each time $t$ and $b$ is the coefficient of interest to us because it will allow us to compute the ADF test statistic and then determine whether the process is stationary or not. I constructed each of these variables in the following way: $\Delta Y_t$               dif.yt = diff(yt)      ## yt is the object containing all the observations for a specific country      ## (e.g. Australia)      $Y_{t-1}$               yt.lag.1 = lag(yt, -1)      $\overline{Y}_{t-1}$               ybar.lag.1 = lag(c(rowMeans(x)), -1)      ## x is the object containing my entire data frame      $\Delta \overline{Y}_{t-1}$               dif.ybar.lag.1 = diff(ybar.lag.1)      $\Delta \overline{Y}_{t-2}$               dif.ybar.lag.2 = diff(lag(c(rowMeans(x)), -2))      $\Delta Y_{t-1}$               dif.yt.lag.1 = diff(yt.lag.1)      $\Delta Y_{t-2}$               dif.yt.lag.2 = diff(lag(yt, -2)      After constructing each variable individually, I then run the linear regression               reg = lm(dif.yt ~ yt.lag.1[-1] + ybar.lag.1[-1] + dif.ybar.lag.1 +                        dif.ybar.lag.2 + dif.yt.lag.1 + dif.yt.lag.2)     summary(reg)      It is obvious that the explanatory variables in my regression equation differ in length, so I'd like to know whether there is a way in R to make all the variables of equal length (perhaps with a function). Also, I'd like to know whether the procedure I used was correct and if there are more optimal ways.


title: Normality test for small samples

text: I get a little bit confused by the conclusions I can take or not with these small samples. I have been measuring the degradation of a pollutant for 9 days. I measured the remaining concentration and only took triplicates each day, so that gives me a very small sample size. I want to check at which point there is no significant difference between 2 days anymore. I first thought of anova followed by Tuckey. But then when i checked for normality, I did the Shapiro test for normality and my data didnt pass the test. That doesn't allow me to say my samples dont come from a normal distributed population, does it? What does it tell me then? And would the qqplot and histogram tell me more or is that also not very relevant with small samples? And does that exclude me using ANOVA? Would a permutation test make sense? but I read the equivalent of Anova requires a lot of calculation time and that it wasnt used a lot.. I would be gratefull if someone could clarify the situation a little bit.. thank you in advance.


title: How do you do bootstrapping with time series data?

text: I recently learned about using bootstrapping techniques to calculate standard errors and confidence intervals for estimators. What I learned was that if the data is IID, you can treat the sample data as the population, and do sampling with replacement and this will allow you to get multiple simulations of a test statistic. In this case of time series, you clearly can't do so because autocorrelation is likely to exist. I have a time series and would like to calculate the mean of the data before and after a fixed date. Is there a correct way to do so using a modified version of bootstrapping? Any other suggestions is also welcomed.


title: Advanced Excel graphs

text: I should not be here, I am a designer by trade and all these numbers scare me. I am attempting to help someone display their information; they sent me an excel spreadsheet and asked me to make graphs out of certain information. Unfortunately the spreadsheet contains ~13,000 rows and 32 columns, and I need to find out how to group them so I have all 0-25 values in one bar, 26-35 in another, etc. I do not know how to even begin to calculate this information or organize it (other than by hand) can someone point me in the right direction as to what I need to do?


title: Expected value of spurious correlation

text: We draw $N$ samples, each of size $n$, independently from a Normal $(\mu,\sigma^2)$ distribution. From the $N$ samples we then choose the 2 samples which have the highest (absolute) Pearson correlation with each other. What is the expected value of this correlation ? Thanks [P.S. This is not homework]


title: Multiple testing correction alternatives

text: I have a data set with about 6,000 Fishers exact tests. If I want to do a multiple testing correction to control the FWER what are the advantages of one of these methods vs the other calculating a q-value for each test Cite Using a Benjamini-Hochberg control of the FDR Hochberg's step up method or Holm step down method Cite Sorry if this is a stupid question, I'm a not a classically trained statistician


title: SVD in movie recommendation

text: Assume that there is a 5 $\times$ 6 matrix that records the ratings of six users on five movies. I have computed the singular value decomposition (SVD) for such a matrix. Suppose I add another movie so now it becomes 6 $\times$ 6 matrix. However, I only have a rating of that movie for two people. How can I estimate the ratings the other users would give to this new movie?


title: How to calculate the standard error of the marginal effects in interactions (robust regression)?

text: what I am interested in learning is how to calculate the std error of the marginal effects of a X variable when it is part of an interaction, especially in robust regression. There are tipically two cases that interest me: when there is an interaction between two continuous variables and when there is an interaction between a continuous and a binary variable. Let's our model be a simple one: $Y= b_0 + b_1X + b_2Z + b_3XZ + e$ I know how to calculate the marginal effect of X in both cases, but not how to calculate its standard error. Please, any hints on how to do this, both theoretically or in R code, may be very helpfull.


title: The f-distribution in Working-Hotelling confidence bands

text: Okay, so I'm having a conceptual issue with using the f-distribution to build a Working-Hotelling confidence band. It makes sense to me that we use a tα/2 or a zα/2 value when building a standard confidence interval, as there are two tails to equally take into account. Now, doesn't the f-curve also have two tails? Why would I use a 1-α value of f when building the CI? What about the other tail? Shouldn't I have to account for both tails? Why does (for example) a 95% W-H band use a α-value of 0.05? What am I missing here?


title: What is the optimal distance function for individuals when attributes are nominal?

text: I do not know which distance function between individuals to use in case of nominal attributes. I was reading some textbook and they suggest Simple Matching function but some books suggest that I should change the nominal to binary attributes and use Jaccard Coefficient. However, what if the values of the nominal attribute is not 2? what if there are three or four values in that attribute? Which distance function should I use for nominal attributes?


title: MDP value iteration

text: In Markov decision processes, what is the guarantee that value iteration chooses the same policy action from a given state for every iteration? I am referring to the slides given by AWM at http://www.autonlab.org/tutorials/mdp09.pdf. Thanks Udaya


title: Markov chains #2

text: A very new European “Rapid Reaction Force for Fire” has been created today and begins operation between three Countries “A”, “B” and “C”. It’s main resource is a super aircraft “Funderbird2” with a massive water cannon that even carries a small mini-submarine for fighting fires at sea. Unfortunately, it can only be in one Country at a time. ![enter image description here](http://i.stack.imgur.com/sXDkM.jpg) What is the matrix representation of this problem? The probabilities of going to a fire in another Country, given that the force are in a given Country to begin with, are shown on the sketch above (any resemblance to any particular Country or Nation is purely coincidence and not intended). The probabilities were obtained as a weekly average using statistics for fires over many recorded weeks. When will the plane take up permanent residence in one Country? Footnote: The final time computed here is the “Laplace transform” of the rate of rescue operations for country C. In this question the matrix approach is used to represent the differential calculus. This is my homework. I need help from where to start. Would like you to give me directions and to be solved by me, so that I can learn something. I think I have solved under a) but I am not sure if matrix is correct.               [1.0 0.0 0.0]     [0.1 0.6 0.3]     [0.1 0.4 0.5]


title: Modeling multinomial problems with unknown sample size in BUGS

text: I'm trying to estimate the design effect of a series of relatively small sample size surveys ($n\sim 70$) with multiple responses. Design effects roughly correspond to how much larger actual sample variance than would be expected from naive random sampling. The simplest way to parametrize this is for Effective Sample Size to be equal to SampleSize/D, where D is an unknown parameter. Throwing it naively into BUGS is impossible due to the inability to handle multinomial models with unknown sample size. I tried to do a multinormal approximation with a structured covariance matrix, but I can't figure out/find a closed form for the precision matrix of a multinomial variable. Does anyone have any ideas about the best way to proceed? This seems like it should be a common problem.


title: Getting started with creating an index based on multiple survey items

text: I'm compiling a survey that will have several questions which would lend to the creation of an index of a main dependent variable (level of engagement in sucession planning). The questions will involve topics like:   * PROCESS: linking strategic planning to succession planning, identifying critical positions that need to be filled,identifying competencies, identifying high potential employees, coaching and mentoring, providing leadership development;    * ROLE IDENTIFICATION: identifying roles (several options); and    * RESULTS: zero people ready for each position, 1 or 2 people ready for each position, or a pool of people ready for each position.  I'm not sure how to start. Any suggestions?


title: Is MANOVA the correct way to handle multiple response variables that are additive?

text: ### Context: I am trying to analyze an experiment on plant community response to two treatments. Here’s a simplified description of the experiment, there are a few extra complications in reality. Treatments were applied to small patches of ground arranged in blocks with a mix of naturally occurring plant species present in each patch. I have measured biomass of the plants in each patch in years one and two as the response variables. Analyzing total biomass seems simple enough with a repeated measures GLM or linear mixed model and block a random effect. However I am also interested in breaking down total biomass into four response groups of interest: perennial grasses, perennial forbs, annual grasses and annual forbs (forbs are plants that aren't grasses). ### Question:   * Because these groups sum up to total biomass can I go a head and analyze these responses with four separate repeated measures ANOVAs, or will I be guilty of re-analyzing the same experiment?  It’s been suggested to me that I do a MANOVA followed by “protected ANOVA’s” but this doesn’t make sense to me because I am not interested in the multivariate response, rather I am interested in testing whether these specific variables have been affected by the treatments. I also don’t know how to do a repeated measures MANOVA correctly and it seems unduly complicated. Any general advice would be appreciated—I’ll try and design a more easily analyzable experiment next time.


title: Evaluate Quantile Regression results and t-stats at different levels

text: I have a few questions regarding quantile regression and how to interpret the results. I have several independent variables and I want to figure out which combination of variables have the highest predictive value. I have preformed quantile regression, in Eviews, on all possible combinations of the 10 independent variables I have. I would like to know which statistical measure provides the best information with regards to evaluate which combinations of independent variables I should use? The t-statistics seems to be more signification at the extreme quantiles (90, 80, 20, 10%), than the quantiles more centered in the middle (40, 50 ,60%). Why is it so?


title: Constructing a VECM with a mix of I(0) and I(1) variables

text: I've been using the Johansen Procedure to check and correct for cointegration in my model, by estimating a VECM instead of VAR. But now I want to estimate a new model, in which I expect the same cointegrated relationships, however, some of my variables are now stationary (because they are one of the former variables, split into growth and decay), is there a way to construct the Error Correction Term and insert it into my differenced data, and then estimate the VECM by ordinary VAR? And what are the implications of ignoring cointegrated relationships and just estimating the VAR? Is my model subject to bias or inconsistency? Thanks.


title: What factors do I have to consider when deciding if creating a composite score using four t-scores is defensible psychometrically?

text: Edit: I've conducted an item analysis on four variables - two quizzes and two psychological attributes, namely emotional intelligence and happiness. Its a within-subject design with students as participants. They all have Cronbach's alpha above .70 or have been revised to meet that criteria. I have a scenario where a selection panel chooses to select students based on their performance in all four composite variables. It was then decided that we should calculate an average standardized t-score across the four variables. Would creating a composite score out of four t-scores be psychometrically defensible? why or why not? I've also correlated all for variables and they were all statistically significant with the lowest being .30. Is that considered decent enough? the quiz scales and attribute scales are different but if using t-scores then would they be considered standardized? thanks for any help!=D


title: Find out the conditional probability

text: Consider I have the following probabilities: $$P(A|B) = 0.86 $$ $$ P(A|B^C) = 0.35 $$ $$ P(B) = 0.80 $$ $$ P(A) = 0.758$$ Is there necessary information given to calculate $P(B^C|A^C)$? If so please guide me how. Thanks in advance.


title: Meaning of p-value of logistic regression model variables

text: So I'm working with logistic regression models in R. Though I'm still new to statistics I feel like I got a bit of an understanding for regression models by now, but there's still something that bothers me: Looking at the linked picture, you see the summary R prints for an example model I created. The model is trying to predict, if an email in the dataset will be refound or not (binary variable `isRefound`) and the dataset contains two variables closely related to `isRefound` , namely `next24` and `next7days` \- these are also binary and tell if a mail will be clicked in the next 24hrs / next 7 days from the current point in the logs. The high p-value should indicate, that the impact this variable has on the model prediction is pretty random, isn't it? Based on this I don't understand why the precision of the models predictions drops below 10% when these two variables are left out of the calculation formula. If these variables show such a low significance, why does removing them from the model have such a big impact? Best regards and thanks in advance, Rickyfox ![enter image description here](http://i.stack.imgur.com/oiCrN.png) * * * ## EDIT: First I removed only next24, which should yield a low impact because it's coef is pretty small. As expected, little changed - not gonna upload a pic for that. Removing next7days tho had a big impact on the model: AIC 200k up, precision down to 16% and recall down to 73% ![enter image description here](http://i.stack.imgur.com/583nx.png)


title: Showing that the maximum likelihood estimator (MLE) exists but is not unique

text: I have a few questions with regards to a solution to the problem below: ![Question and Answer:](http://i.stack.imgur.com/aTeir.png)   1. How is it possible to have $max_{1\leq{i}\leq{n}}x_i-1<\theta<min_{1\leq{i}\leq{n}}x_i$? How can a value of $\theta$ be both greater than a larger value and less than a smaller value? Isn't that contradictory?   2. How is it that $max_{1\leq{i}\leq{n}}x_i-1 < min_{1\leq{i}\leq{n}}x_i $? For example, what if my $max_{1\leq{i}\leq{n}}x_i = 10$ and $min_{1\leq{i}\leq{n}}x_i = 1$? Would this then mean that $10<1$?   3. Finally, if $x_1,...,x_n= 1,0,1,0,1,0,...,1,0$, where the maximum is $1$ and the minimum is $0$, wouldn't it mean that $0<\theta<0$ and so implies that $\theta = 0$ and thus is unique? Thanks everyone!!!


title: spss: working with two binary/dummy variables

text: Am trying to set a few binary/dummy variables against each other, i.e. `propensity_to_dance` and `gender`. I assume that it' ok to be using a binary logistic regression in SPSS (or otherwise) to investigate the interaction on these two variables; I've previously been using that kind of regression to investigate the relationship between binary and linear variables. Is it correct to use this methodology for two binary variables as in the above?


title: Adding content information to matrix factorization-based recommender

text: I'm currently using a matrix factorization method to generate recommendations (for info on this, check: Matrix Factorization Techniques for Recommender Systems). At the moment, my rating estimate is given by the global mean, user and item biases, and the latent factors. I'm interested in adding content information to this system. In this case, both my items and users can be represented in a topic (word) space. I was thinking about adding the correlation between user and item (some value in [0,1]) as an additional input source to my rating estimation. My question is: how to best include this information? just add it to the sum? Is it not a problem if the correlation is a value between 0 and 1, and the global mean and biases are values that can be as big as 100 ? Thanks!


title: Introductory textbook on nonparametric Bayesian models?

text: I'd like to wrap my head around this topic but learning from white-papers and tutorials is hard because there are many gaps which are usually filled in textbooks. If it is important I have relatively strong mathematical background as I did my Ph.D. in applied mathematics (CFD to be more precise).


title: How to handle samples with many replicates for population estimation,when sample size is small?

text: I have a sample with very small sample size ($n<10$) and I would like to estimate the population mean. However, some of the samples are replicates (for instance, sampled at the same location). According to a previous discussion: Sampling with or without replacement? If my sample size $n$ is large, I could consider my case as sampling with replacement, which provides the same unbiased mean estimation but a larger variance estimation, as compared to sampling without replacement. However, in my case, the sample size is really small. Thus the estimated mean would be biased due to the replicates. Any suggestion about how to handle this case? Shall I exclude the replicates? Or is there any general rule about the proportion of replicates required to be a random sampling with replacement? Thanks


title: Maximum likelihood solution in classification problem

text: I have a couple of questions regarding the maximum likelihood solution in a classification problem (with only two classes $C_{1}$ and $C_{2}$) Basically, I have the following likelihood function: $$p({\bf{t}}|\pi, \mu_{1},\mu_{2},\Sigma)=\Pi_{n=1}^{N}[\pi N(x_{n}|\mu_{1},\Sigma)]^{t_{n}}[(1-\pi)N(x_{n}|\mu_{2},\Sigma)]^{1-t_{n}}$$ in which $\bf{t}$ is the target vector, $P(C_{1})=\pi$ and consequently, the probability $P(C_{2}) = 1-\pi$. $C_{1}$ and $C_{2}$ have mean $\mu_{1}$ and $\mu_{2}$, respectively and both share the same covariance matrix $\Sigma$. Finally, $t_{n}$ is equal to one if the input vector $x_{n}$ belongs to $C_{1}$ and $0$ otherwise. The first question is: Why the likelihood function seems to be the product of the joint distributions of $P(x_{n}, C_{1})$ and $P(x_{n}, C_{2})$: $$P(x_{n}, C_{1})=P(C_{1})P(x_{n}|C_{1}),\;\; P(x_{n}, C_{2})=P(C_{2})P(x_{n}|C_{2})$$ The second question deals with the calculation of the maximum likelihood for a covariance matrix $\Sigma$. In the book I'm following (Bishop, Pattern Recognition and Machine Learning), it describes the solution in the following way: ![enter image description here](http://i.stack.imgur.com/r1gJV.png) with $$S = \frac{N_{1}}{N}S_{1}+\frac{N_{2}}{N}S_{2}$$ $$S_{1}= \frac{1}{N_{1}}\Sigma_{n\in C_{1}}(x_{n}-\mu_{1})(x_{n}-\mu_{1})^{T}$$ $$S_{2}= \frac{1}{N_{2}}\Sigma_{n\in C_{2}}(x_{n}-\mu_{2})(x_{n}-\mu_{2})^{T}$$ How can I prove that the second and fourth term in the equation above can be expressed as a trace $\text{Tr}\\{\Sigma^{-1}S\\}$. I know that the first and third term give $-\frac{N}{2}\ln|\Sigma|$. Thanks a lot!


title: Metropolis Algorithm

text: I am currently doing my statistics thesis on modelling football data which requires quite a great knowledge of Bayesian Theory especially MCMC methods. However I have some minor problems regarding the Metropolis algorithm:   * Any distribution can be attributed to the proposal density function?   * Also, the initial values given to the parameters (i.e. the first step of the Metropolis algorithm) are random or there must be some thinking behind the choice of values? Thank you


title: Derivation of a the log-likelihood for a regression model where the outcome is a mixture between Poisson and a point mass at zero

text: Suppose $ \textbf{Y} = (Y_1, \dots, Y_n)'$ are independent and $$\eqalign{ Y_i = 0 & \text{with probability} \ p_i+(1-p_i)e^{-\lambda_i}\\\ Y_i = k & \text{with probability} \ (1-p_i)e^{-\lambda_i} \lambda_{i}^{k}/k! }$$ and $$\eqalign{ \log(\mathbf{\lambda}) &= \textbf{B} \beta \\\ \text{logit}(\textbf{p}) &= \log(\textbf{p}/(1-\textbf{p})) = \textbf{G} \mathbf{\gamma}. }$$ How do we derive the log-likelihood of $\textbf{Y}$? I know that it is $$L(\gamma,\mathbf{\beta} | \textbf{Y}) = \sum_{y_i=0} \log(e^{G_i \gamma}+\exp(-e^{\textbf{B}_i \mathbf{\beta}})) +\sum_{y_i >0} (y_i \textbf{B}_i \mathbf{\beta}-e^{\textbf{B}_i \mathbf{\beta}})-\sum_{i=1}^{n} \log(1+e^{G_{i} \gamma})-\sum_{y_i >0} \log(y_{i}!)$$ **Added.** I think that the likelihood function will be $$\prod_{Y_i=0} p_i+(1-p_i)e^{-\lambda_i} \times \prod_{Y_i >0} (1-p_i)e^{-\lambda_i} \lambda_{i}^{k}/k!$$ Then we know that $\lambda = \exp(\textbf{B} \beta)$ and $p = \frac{1}{1+\exp(-G \gamma)}$. So we just substitute these values and take logs?


title: Automated procedure for selecting subset of data points w/ strongest correlation?

text: Is there some standard procedure (such that one might cite it as a reference) for selecting the subset of data points from a larger pool with the strongest correlation (along just two dimensions)? For instance, say you have 100 data points. You want a subset of 40 points with the strongest correlation possible along the X and Y dimensions. I realize that writing code to do this would be relatively straightforward, but I'm wondering if there's any source to cite for it?


title: Sheppard's correction

text: Is there a good expository account of Sheppard's correction, written in a way that any ordinary mathematician can readily follow? http://mathworld.wolfram.com/SheppardsCorrection.html (I've thought of writing a Wikipedia article titled "Sheppard's correction", but I haven't done enough homework on it yet.) Sheppard's correction can be done for various functionals of a probability distribution. The Wolfram article treats cumulants. Curiously, each correction is $-1$ times the correction I'd expect to do for a uniform distribution. I'm not sure whether to believe this. Think of this: split the interval $(0,\theta)$ into bins of equal length. For a sample from the uniform distribution on the interval, suppose only the midpoint of the bin is reported. Then the variance of the reported data would be too small, since it would fail to report variability _within_ bins. The variability within bins would be just $1/12$ of the square of the bin length. But Sheppard's correction for normally distributed data tells you to add $-1/12$ of the square of the bin length. With the normal distribution, the error from binning is negatively correlated with the observation; with the uniform distribution it is uncorrelated with the observation. If we believe Wolfram, this multiplication by $-1$ seems to apply to all cumulants, not just the second cumulant. If so, that seems like a mystery to be investigated. I hesitate to believe Wolfram on this.


What are the standard statistical tests to see if data follows exponential or normal distributions?


title: Gamma function, unbiased estimation of standard deviation

text: According to the Wikipedia article, the following estimator of the standard deviation $$s=\sqrt{\frac{1}{n-1}\sum_{k=1}^n(x_i-\bar{x})^2}$$ for a normal variable, verifies $E[s]=C_4(n) \sigma$, where $$ C_4(n) = \sqrt{\frac{2}{n-1}}\cdot\frac{\Gamma(\frac{n}{2})}{\Gamma(\frac{n-1}{2})}=1 - \frac{1}{4n} - \frac{7}{32n^2} - \frac{19}{128n^3} + O(n^{-4}) $$ How could you prove this?


title: optimisation procedures before training SVM

text: I'm using the LIBSVM in Java for classification with 200 documents in inputs. I build/train the SVM using the same input training data. My response time for preprocessing of documents (tokenization, filtering and ngrams and then writing results out to files) and loading file into memory is 8 second. Could anyone recommend a good methods (maybe tools) for preprocessing and loading my result into memory ?? Thanks in advance for any comments.


title: Working through Lenth's maximum likelihood estimation in JS

text: I am attempting to write a program in Javascript to estimate where a radio- tracked animal is. Basically, it should take 3+ bearings as inputs, along with the locations they were taken, and output a good estimate of where the animal is. From what I've read, it seems like Lenth's 1981 method of Maximum Likelihood Estimation ('On Finding the Source of a Signal' PDF) is the best for this situation. I also need to find the area of a 95% confidence ellipse (and ideally, draw it on a map) so the folks collecting these bearings will know if the data they took was good enough. However, reading through Lenth's paper as a non-statistics person, I'm having a hard time working through the MLE. My goal is to use the data in Table 1 and calculate the estimate (7.23, 1.98), as Lenth said is what the unidirectional "clean data" technique yields. Then I'll feel confident with the math. As far as I understand, Fig. 2.6 in the PDF seems to be what I need to solve to calculate $x$ and $y$, though it's a little terrifying. Using the data in the first row of Table 1, I believe... $\theta_1= -2.5132741228718345$ (These values are in radians) $s_1 = -0.5877852522924732 $ $c_1 = -0.8090169943749473 $ $d_1 = [(x-9)^2-(y-3.2)^2]^{1/2} $ $s_1^* = (y-3.2)/d_1^3$ $c_1^* = (x-9)/d_1^3$ $z_1 = -0.5877852522924732*9 - (-0.8090169943749473*3.2) = -2.70121289$ So... $s_1s_1^* = (−0.5877852522924732)*(y-3.2)/([(x-9)^2-(y-3.2)^2]^{1/2})^3 $ And this is where I get discouraged because I have barely even started the summations and have no idea how to solve for the mess that $x$ and $y$ will be after the summations are completed and the matrices are multiplied. Is this even possible with simple math functions, or do you need a stats program like R to do this sort of math? I know there are some Javascript libraries (e.g. http://www.jstat.org/) with advanced statistics functions, but I have no idea if they'd be useful to me since I'm not really familiar with stats terminology. Any chance you folks might be able to point me into the right direction? _(My apologies if any of these questions are ridiculous. I'm very much out of my element with this sort of math.)_ \-- Edit: Apparently I neglected the bit about how to find a starting estimate and then using it to iterate over the formula until it converges to a single point ($x$, $y$). Thank you, naught101. I had a feeling I was doing this wrong.


title: How to identify a subset of effects as the drivers of significant differences?

text: Okay I will explain this as best as I can. I have an additive risk score for each of my observations (i.e. I'm putting 100 individual effects into a single variable). The means of this risk score variable are significantly different between sets of observations from two different populations. I now want to go a little deeper and ask the question of whether different sets of those 100 are driving the risk in each population. I don't want to test them individually for two reasons: I don't have the statistical power to test them individually, and when I keep them together I am also using information about their "clustering" in the same observation. I was thinking of writing a program to see which of the effects I can drop and still get a significant difference, by computing the additive risk score sequentially and seeing if it can still differentiate the two populations. I was thinking that I could do this iteratively, once I identify one effect that doesn't contribute as much, I drop it, and continue dropping until I have just those that are driving the difference. What do you think about this design of seeing which ones I can drop? Is this a reasonable thing to do? Also, if there is literature on this kind of analyses please let me know, I must be using my own made-up terms to search and they are not turning up any relevant results. I'm not bound by anything, so absolutely any other ideas that you may have as to how I can further explore this are more than welcome! I hope I have provided enough information to describe the problem, if not please let me know. Thank you for thinking with me :)


title: Why is this likelihood function equal to the noise PDF?

text: My professor has this slide up here: ![enter image description here](http://i.stack.imgur.com/ok8iy.png) Here, $y$ is an observed signal. $H$ is a deterministic transformation, which is assumed known. $f$ is the original signal (which we dont know), and $w$ is random gaussian noise. We are trying to recover $f$. I understand everything, except for, why $p(\mathbf{w})$ = $p(\mathbf{y}|\mathbf{f})$. That is, I understand that the multidimensional noise PDF is given by the above expression. But why is that expression, ALSO equal to the likelihood function, $\mathbf{y}$, given $\mathbf{f}$? I'm not seeing this...


title: Expected assortativity coefficient of random graphs

text: I am wondering whether there exist a closed-form expression for the expected assortativity coefficient (http://arxiv.org/pdf/cond-mat/0205405.pdf) of an Erdos Ranyi random graph model $G(n,m)$, where $n$ is the number of nodes and $m$ is the density (http://en.wikipedia.org/wiki/Erd%C5%91s%E2%80%93R%C3%A9nyi_model). A rough approximation might be provided by the correlation of a multinomial distribution, but it doesn't seem particularly accurate when the graph density is reasonably high. In fact, the degree of a node cannot be greater than $n-1$, whereas if we regard the degree vector of $G(n,m)$ to be distributed as $Multinomial(p_1, p_2, \ldots, p_n; m)$, where $p_i = m/n$, for $i = 1, \ldots, n$, the result is that a node might have a degree greater than $n-1$ with positive probability. This is clearly wrong! In any case, if we assume the degree vector of $G(n,m)$ to be distributed as $Multinomial(p_1, \ldots, p_n; m)$, the result is that the degree correlation (which represent the assortativity coefficient) is $\frac{1}{n-1}$. I think that many of you would probably claim that the best choice is to use a truncated multinomial, $P(d_1, \ldots, d_n | d_1 < n-1, \ldots, d_n<n-1)$, and to derive the correlation of this multivariate expression. Nonetheless, before undertaking such a tedious work I'd like to see whether someone already know a closed-form expression for the expected assortativity coefficient. Thank you very much for your help.


title: Acronyms to use for Bayesian posterior predictive distribution estimators

text: I am considering writing an article that discusses the Bayesian MMSE and MAP of the posterior predictive distribution. I was wondering if there are acronyms that have been used so that instead of always saying: "The MMSE of the posterior predictive distribution..." I might just say: "The MMSPE..." or "The MMSEP" or something. Any thoughts?


title: When is it a good idea to just use the average for imputation?

text: Suppose we have a data set **test** :               1     8     12     14     .     .     19      The ` . ` denotes missing values. When would it be better to use the average of the non-missing values to impute the missing values rather than assuming that the data comes from a normal distribution?


title: Eliminated items in Factor Analysis

text: I'd like to get your opinions on how to interpret items that had to be eliminated in factor analysis (FA). I've been researching consumer shopping motivations and ran a survey with 40 items, which included statements from previous research in the field (there is no developed go-to scale). I expected items would be reduced to 10-12 variables. I ran FA in SPSS (extraction method - principal components) and after all items with low factor loadings and communalities were eliminated, I was left with 25 items and 6 underlying variables (criteria for the Number of factors: eigenvalue>1). Now, what conclusion can I make about statements which dropped-out? Are these (eliminated) motivations not valid for the consumers I ran the survey with? Thanks a lot on your thoughts!


title: I don't understand the figure output from package lme4 in R using the effects library?

text: I am using **linear-mixed effect models** to analyse my data (the interaction between `on` and `arc`; both `arc` and `on` are continuous variables). I also used the `effects` library to plot the interaction of `on*arc`. However, I could not understand the figure: ![enter image description here](http://i.stack.imgur.com/ZQP0l.png) Could you please help in understanding what the 5 panels in the figure represent? Do they represent **quantiles** (e.g., the first 20% lowest arc)? Why there are **five panels** instead of say **four panels**? I understand that grey band represent the 95% **confidence intervals** , and that orange bars represent ranges of arc (high arc to low arc). I also wanted to check whether the location of the orange bars represent where the quantiles fall in arc.


title: World Cup: Is it worth to bet each match twice?

text: Let's assume we get **1 point for predicting the result** (win/lose/draw) and **3 points for predicting correct score for both teams** (e.g. 3:2). We've got two separate accounts where we can bet, my question is: (1) **does betting twice each match** (once per account, each bet with different score) allows you to get more points on one of the accounts than on account when you bet only once? (2) Is it possible to calculate **how probable** it is to get better result with two accounts? **Edit and clarification:** Betting doesn't cost money, it's just a company game. **My friend argued that if he could bet from two accounts he would be much more successful and I think it wouldn't increase his chances to win much since this is not entirely random but I'm lacking any skills to prove the gain would be insignificant.** What I mean in question 1 is that **if you bet on one account you choose the most probable outcome (according to your knowledge) so the second account would always get only less probable outcomes** which shouldn't increase chances since match results aren't totaly random (I was able to correctly guess outcome of almost half of matches and the account with my second guesses would do worse). So I've got two betting models: primary and secondary and I'd like to calculate how likely it is that using two models can give me better overall result on one of the accounts than using just primary model. I think it may be impossible but I'm a noob so maybe you guys can help me.


title: Unknown null-distribution, known alternative distribution

text: I'm rather new to statistics, so please bear with me. I would like to be able to say "to within 95% confidence, this data is exponentially-distributed." In this case, I know how the data _should_ be distributed and I want to be able to say it _is_ to some confidence. In other words: (Usually)   * $H_0$ = The data are exponentially-distributed.   * $H_a$ = The data are not exponentially-distributed.   * Choose a test statistic $T_\alpha$ such that Pr(T > $T_\alpha$) = $\alpha$ (under $H_0$).   * If T > $T_\alpha$, one can reject H_0 with confidence 1-$\alpha$. (What I'd like)   * $H_0$ = The data are not exponentially-distributed.   * $H_a$ = The data are exponentially-distributed.   * ???   * Reject H_0 with confidence 1-$\alpha$. I've also considered doing linear regression on the empirical CDF for the data, but I'm not sure how I would interpret the results of that in terms of confidence. I apologize if this has been answered before, but I've searched considerably and haven't found anything helpful. Thank you all for your help.


title: Normal Distribution with random mean and standard deviation

text: When trying to code this in R, I'm getting very confused about what to do. Apologies if my terminology is incorrect but I would be grateful for any advice. **The Problem:** I have been given two normal distributions for the mean and standard deviation of rat weights. So:   1. The mean rate weight is itself a normal distribution with a mean of 1.68 and a standard deviation, or confidence interval, of 1.81 (mean = 12.68, SD = 1.81)   2. The standard deviation is itself a normal distribution with a mean of 11.19 and a standard deviation of 3.2 (mean = 11.19, SD = 3.2) _Summary_ Mean Rate Weight: dnorm(Mean = 12.68, SD = 1.81) Standard Distribution of Rat Weight: dnorm(Mean = 11.19, SD = 3.2) **The Question:** In R, how to a code this to have a Monte Carlo run of 50,000 samples? Is the following example correct?                MC_Runs = 50000           Rat_Weight = rnorm(MC_Runs,                        mean = rnorm(MC_Runs,mean = 12.68, sd = 1.81),                        sd= rnorm(MC_Runs,mean = 11.19, sd = 3.2))


title: Repeating many ANOVA tests with a large number of factors

text: I have a large number of factors that categorize my numerical observations in different ways and no preconception which if any of these factors may be significantly explanatory, but I have a p-value threshold that I'm obliged to observe. I want to perform an analysis of variance for each factor and then select those factors that have significant p-values according to my threshold. Because I'm running so many ANOVA tests, I'm worried about the possibility of a false discovery. Is this concern legitimate? And if so, what techniques are there for adjusting my p-value threshold to take account of the fact that I'm performing so many tests?


title: Logistic Regression Model Validation

text: I am validating a logistic regression model. This is the first time i am validating a model. I am using split sampling method. I have split data randomly into two parts - 70% development and 30% validation data sets (70:30). Then i run logistic regression on development data set using SAS and rank their probabilities in descending order and split data into 10 groups (deciles). Check the percentage of responses in the upper deciles. My question - Is there any thumb rule to assess the model on the basis of percentage of responses in the upper 3-4 deciles? Someone said top 3 deciles should cover atleast 65% of responses. Is it correct? I have checked Hosmer and Lemeshow Goodness-Fit-Test. I have formulated an equation including intercept and coefficients derived from development dataset and run it on validation dataset. The code is shown below for reference - Data validation_Output; Set validation_Set; resp_xb1= -1.3844+(1.4708) _A_Flg+(2.9829)_ B_Flg+(.0317)* C + (-0.2372)*D +(-0.3359)*E; exppres1=exp(resp_xb1); p_resp1=exppres1/(1+exppres1); run; Then run PROC RANK and PROC SQL to calculate deciles on validation dataset. I have decile scores on development and validation sets. Should the significant variables in both the datasets be same? OR Concordant?


title: Prediction in support vector regression

text: Has anyone attempted prediction using support vector regression? I'm using LIBSVM, but I'm not sure how to use SVR in either univariate and multivariate time series. Say we have stock prices for $N$ days. For training inputs, $y$ are the stock prices for $N$ days, but what will we use for $x$?   1. Time series? For i.e. in one step ahead prediction $1,2,3...Z$ for $Z$ days?   2. (for one step ahead) sifting one day of $y$ values? To explain more:               matlab> model = svmtrain(training_label_vector,                               training_instance_matrix [, 'libsvm_options']);      For univariate: I use the stock prices for $N$ days in `training_label_vector` as a column vector and want to predict say next 30 days. I wonder which data I have to use in `training_instance_matrix`? For multivariate: say I have 22 more features (prices of other goodies), I use other features as column vectors in `training_instance_matrix`. But I'm not sure if I'm using the correct approach.


title: Why can log-linear models for binary data be trained using Poisson regression

text: Any arbitrary distribution of a multivariate binary variable can be represented as a log-linear model. That is, for $X = (X_1,\dots,X_d)$ a $d$ dimensional binary rv, the distribution can be written as $\log(\text{p}(x)) = \sum\limits_{A\subseteq V} \beta_A (\prod\limits_{i \in A} x_i)$ where $V={1,\dots,d}$. E.g. for $X = (X_1,X_2)$, $\log(\text{p}(x)) = \beta_\emptyset + \beta_1 x_1 + \beta_2 x_2 + \beta_{1,2} x_1 x_2$. Note, $\beta_\emptyset$ is a normalizing term. I've read that this model can be trained using Poisson regression. For example, in R               m_out = glm(count ~ x1 + x2 + x1*x2, family="poisson")      for the 2d case above. I'm unclear about why this is so: why does modeling the counts as being distributed as Poisson with mean = $\beta_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_{1,2} x_1 x_2$ work? In particular, let $\beta(x)=\log(\text{p}(x))-\beta_\emptyset$ if we write out the likelihood for our betas we have that: $L(\beta) = \prod\limits_{i=1}^m e^{\beta(x_i)+\beta_\emptyset}=\exp( \sum\limits_{i=1}^m (\beta(x_i)+\beta_\emptyset))$ so the log likelihood is : $l(\beta)=\sum\limits_{i=1}^m (\beta(x_i)+\beta_\emptyset) = \sum\limits_{A\subseteq V} y_A(\beta(x_A)+\beta_\emptyset) = m\beta_\emptyset+\sum\limits_{A\subseteq V} y_A\beta(x_A)$ where $x_A$ is the vector $x$ such that $x_i=1$ if $i \in A$ $x_i=0$ otherwise and $y_A$ is the empirical count of how many observations had the pattern $x_A$ and $\beta_\emptyset = -\log(\sum\limits_{A\subseteq V} e^{\beta(x_A)})$ is the normalizer. However, according to Wikipedia what we are doing with glm is modeling our counts $y_A \sim \text{Pois}(\beta(x_A)+\beta_0)$ where $\beta_0$ is an intercept term. Ignoring constants, this will lead to a log likelihood of: $\sum\limits_{A\subseteq V} (y_A(\beta(x_A)+\beta_0)-e^{\beta(x_A)+\beta_0})$. Is there a reason why maximizing both likelihoods will give the same coefficients (at least up to constants)?


title: Effect size in R (check significance level for bigger sample)

text: In my current analysis, I have not found statistical significance to some of my conditions. So I was wondering how to predict or calculate the effect size with more subjects and items. I have 20 subjects and each have taken 120 items. My dependent variable is Reaction Time (RT) and my independent variables are Related (2 levels -- control and experimental) and PrimeType (5 levels). In testing significance, I look at RT significance of Related for each Prime Type level. For example, the significance level for PrimeType 3 when compared in RT between Related yield a p-value = .075 So I wonder what packages or functions that I could use where I put in the current observation numbers and significance level and compare it or calculate it with a bigger sample. I want to check whether the the significance level would change or become statistically significant in a bigger sample size. Any ideas?


title: training period selection forecast (error analysis)

text: I have been lately testing the best training period length to perform a forecast. I have tested it for various days of training period length, among them 60 days and 30 days. My methodology is quite rigorous. I have the nrmse growth for the 6 hours of forecast for each training period length. The nrmse growth is telling me that the 60 days of training period length perform sligtly better at each step of the forecast than 30 days. On the other hand, my overall average error is saying that with 60 days of training period I have 10% of NRMSE and with 30 days 8% of NRMSE. So considering NRMSE growth 60 days is better, and considering the overall average 30 days is better. I do not know which one to choose. Can anyone enlighten me here? ***I am doing a 6 hours ahead forecast with a moving window every six hours, this for a period of 8 month.


title: Why is "stationarity" assumed in time series data?

text: A stochastic process is composed of a sequence of random variables ordered by time,   and a time series is just a realization of such a process.   The book that I'm reading says: "if we assume stationarity, then we can get expectation and variance from time series data". I don't understand this. Doesn't stationarity mean that every random variable in the stochastic process has same distribution? If not, every time series data comes from different distributions how they can calculate expectation and variance from them?


title: Best model for many independent variables

text: Let's assume I don't know of the existence of clocks, so I want to build a model to predict the current time of day based on a large amount of other things I can measure, for example pressure, temperature, noise, light, and so on. Some of these variables may be directly correlated with time, some may show some sort of peak shape and some may just be noise. There are about 10,000 variables in total. So as my model input I have a big matrix that gives me the measurements for each of the 10,000 variables at each of the 24 hour of the day. What would be the best way to go about building a model for this? Should I try to find a small set of variables for which I can build a linear model? Should I do a PCA and build a linear model on the first few principal components? Something else?


title: What to do when an independent variable is not significant, but it definitely should be!

text: I would appreciate your help with this problem. I am trying to replicate this study which is about finding the optimal level of bank capital. All this is done using Stata. The only difference is that I have quarterly data, and the authors have semi-annual data. We have a problem with the following regression:               Beta=alpha(i) + X(i,t-1)b + u(it)      X(i,t-1) represents a very much documented important predictor of beta. The process of estimation is the following:   1. Perform a unit-root test to make sure beta and X do not have a spurious link. _We performed the test and we reject the H0, therefore all good up to here._   2. Perform the regression using OLS, Fixed Effects and Random Effects. When we run the regressions, following every single step indicated in the study, we get that X is not significant. Not even at the 80% level! For any of the three methods. (Please see from the end of page 13 to page 15 for more details about the regression). Our panel data looks something like this: ![Panel data](http://i.stack.imgur.com/yrpMY.png) According to what is indicated in the paper, this is the regression we run:               . xtset bank quarters, quarterly          panel variable:  bank (strongly balanced)     time variable:  quarters, 1998q1 to 2013q1     delta:  1 quarter      \--               . xtreg beta ib(first).year LD.leveragefull, fe vce(cluster bank)      And this is what we get:                   Fixed-effects (within) regression               Number of obs      =       295     Group variable: bank                            Number of groups   =         5          R-sq:  within  = 0.2965                         Obs per group: min =        59            between = 0.0002                                        avg =      59.0            overall = 0.1387                                        max =        59                                                          F(4,4)             =         .     corr(u_i, Xb)  = 0.0000                         Prob > F           =         .                                             (Std. Err. adjusted for 5 clusters in bank)     ------------------------------------------------------------------------------                  |               Robust             beta |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]     -------------+----------------------------------------------------------------             year |            1999  |  -.2352546   .0679795    -3.46   0.026    -.4239961   -.0465131            2000  |  -.3410705   .1256874    -2.71   0.053    -.6900346    .0078936            2001  |  -.2187538   .1399625    -1.56   0.193    -.6073519    .1698443            2002  |  -.2628397    .094347    -2.79   0.050     -.524789   -.0008903            2003  |  -.1413272    .050201    -2.82   0.048    -.2807076   -.0019469            2004  |  -.2506764   .0697621    -3.59   0.023     -.444367   -.0569858            2005  |  -.3280419   .1220112    -2.69   0.055    -.6667993    .0107156            2006  |  -.2976486   .0996758    -2.99   0.040     -.574393   -.0209043            2007  |  -.1901107   .1236727    -1.54   0.199    -.5334811    .1532598            2008  |  -.1407684   .0840705    -1.67   0.169    -.3741854    .0926487            2009  |  -.0503666   .0505586    -1.00   0.376    -.1907398    .0900066            2010  |  -.0764699    .057999    -1.32   0.258    -.2375011    .0845613            2011  |  -.0035166   .0720702    -0.05   0.963    -.2036155    .1965823            2012  |   .0006204   .0492506     0.01   0.991     -.136121    .1373619            2013  |   .3996657   .0604211     6.61   0.003     .2319099    .5674215                  |     leveragefull |              LD. |   .0047472   .0068057     0.70   0.524    -.0141484    .0236428                  |            _cons |   .5461694    .069342     7.88   0.001     .3536452    .7386935     -------------+----------------------------------------------------------------          sigma_u |  .29632318          sigma_e |  .21633914              rho |   .6523096   (fraction of variance due to u_i)     ------------------------------------------------------------------------------      The authors of the study, on the other hand, get very significant values for leveragefull (X), and this makes absolute sense with the theory and practice. So the question is, what am I doing wrong? What can be the problem? How can I fix either the data or the regression in order to increase the significance of the results?


title: Variance stabilization "rule" for MCMC jumps...anyone?

text: I have an implementation of an MCMC algorithm (Metropolis-Hastings and Adaptive Metropolis-Hastings) that I want to modify to suit my needs (it's pyMC, if anyone is interested on the details). My problem is that the parameter space for my model is very high (10-20 parameters), so the burn-in period is crucial for each MCMC chain to converge. I've read that for a multivariate gaussian posterior an acceptance rate of 0.234 is ok (Gelman et. al, 2004), but in order to get those acceptance rates I need to stabilize the variances of the gaussian jumps between proposal values of each of my parameters in this burn-in period: can someone shine some light on rules (and proofs!) on how to stabilize them? So far I've only seen very few rules without proof. The standard rule for the adaptation of the variances in pyMC can be seen here, from lines 562 to 581, which is rather arbitrary. On the other hand, Ford (2006) on section 3.2 proposes to adapt the variances by some sort of rule-of-thumb which I never heard of before. I've also been thinking that by monitoring how the variance affects the acceptance rates one could detect the optimal variances for the jumps... ...help, anyone?


title: Are observations consistent with each other?

text: I have 60 apps. For each app, I look at 10 review of that app, and count how many of the reviews say the app is unreliable. Let $X_i$ denote the number of such reviews, for the $i$th app. I want to check whether the reviews seem to be pretty consistent for each app. For instance, if one review says the app is unreliable, does this make it more likely that the other reviews for that app also mention unreliability? So, one model for the situation is that the reviews don't reflect anything about the apps: there is a parameter $p$, and each of the 600 reviews is an iid Bernoulli($p$) random variable. In this model, $X_i \sim \text{Binomial}(10, p)$ for each $i$. The value of $p$ is not known a priori. An alternative model is that each app has some inherent degree of reliability, modeled by a parameter $p_i$ for the $i$th app, and each of the 10 reviews of the $i$th app is an iid Bernoulli($p_i$) random variable. In this model, $X_i \sim \text{Binomial}(10, p_i)$ for each $i$. The value of the $p_i$'s is not known a priori.   1. If I observe the values of $X_1,\dots,X_{60}$, is there any way I can determine which model seems to better reflect the data? Is there a good way to construct a hypothesis test to compute a $p$-value for the null hypothesis, that $X_i \sim \text{Binomial}(10, p)$ for all $i$?   2. Let's say that I reject the null hypothesis. Can you suggest a meaningful way to measure the degree of internal consistency across reviews of the same app? What would be a good measure of this sort of consistency/correlation? Ideally, something that is amenable to a relatively natural interpretation? (For instance, maybe I compute the sample standard deviation $\hat\sigma$ of $X_1,\dots,X_{60}$; I compute the expected standard deviation $\sigma_0$ under the null hypothesis, $\sigma_0 = \sqrt{10 \hat{p} (1 - \hat{p})}$, where $\hat{p} = (X_1+\dots+X_{60})/600$; and I look at the ratio $\hat{\sigma}/\sigma_0$?)


title: Best way to estimate a regression with reaction times

text: I have reactions times in a no-normal distribution, so I am using logarithms. Which can be the best way to estimate a regression? Later, I would like to plot the results according to the age of my subjects.