title: How to properly consider uncertainty propagation in random variables

text: I have 34 input random variables and one output random variable, named `R`. From 33 input random variables and some dataset I determined the best predictive regression model. Based on this model I can predict `R` from any values of the input random variables. I also determined the best probabilistic models and their respective parameters for the 34 input variables. Based on this I generated random samples using LHS from each pdf. In order to account for the uncertainty both in selecting the best fit probabilistic model and their parameters values I generated the samples based on the following procedure:               Outer loop {          Varies the probabilistic models (best fit model and three other models that also fit well the data) and their parameters (based on 95% CI, assuming a Normal distribution)               Inner loop {              For each model and their parameters value, generate n samples      Now, I have several million samples of each input variable, and using the predict model also of the output variable R. Now, I want to compute:               g =  R - S (where S is the 34th input random variable)     beta = -inverse cdf of normal dist(g = 0)      I have two ways of doing this:               1) Compute beta for each inner loop. I'd get a pdf for beta and not an unique value.     2) Compute beta outside the two loops based on the wider R and S pdfs. I would get only one values, in general corresponding to a lower estimate of beta.      **What is the proper way to do it: option 1) or 2)?**


title: Choosing between different methods when the first one raises error message for linear regression

text: I have a linear regression problem $$Ax=b$$ My initial approach that helped to solve some of my questions was using **SVD** and obtaining the chi-square and some other values that I am interested but it is breaking down for some cases for instance if my regression problem is as following :               >>> coff=     array([[ 1.,  0.,  0.,  0.],            [ 0.,  0., -1.,  0.],            [ 1.,  0.,  0.,  0.],            [ 0.,  0.,  0., -1.],            [ 1.,  0., -1.,  0.],            [ 0.,  0.,  1., -1.],            [ 0.,  0., -1.,  0.],            [ 0.,  0.,  1., -1.]])      `coff` matrix is actually `A` matrix and `b` in regression problem is as following:               >>> b     array([-0.56168673,  0.8943901 , -0.56168673,  1.20952994,  0.33270337,             0.31513984,  0.8943901 ,  0.31513984])      Using singular value decomposition method :               print "==============================SVD calculation============================"     U, s, Vh = linalg.svd(coff, full_matrices=False)     print U.shape, Vh.shape, s.shape     print s     S = scipy.linalg.diagsvd(s, 4, 4)     print allclose(coff, dot(U, dot(S, Vh)))     Sh=scipy.linalg.inv(S)          for i in range(Sh.shape[0]):       if Sh[i,i]>1.0e+04 :          Sh[i,i]=0          Uh = scipy.transpose(U)     V= scipy.transpose(Vh)     aa=dot(V, dot(Sh, Uh))     aah= scipy.transpose(aa)     S_sq=dot(Sh,Sh)     V_sq=dot(V,Vh)     covar=dot(S_sq,V_sq)     #The least square problem results     res=dot(aa,b)          wt=zeros(s.shape[0],float)     for i in range(s.shape[0]):        wt[i]=0        if math.fabs(s[i])>1.0e-04:           wt[i]=1./(s[i]*s[i])          cvm=zeros((s.shape[0],s.shape[0]),float)     for i in range(s.shape[0]):         j=0         while j<=i:              cum=0.0            for k in range(s.shape[0]):              cum=cum+Vh[i,k]*Vh[j,k]*wt[k]                cvm[i,j]=cum           cvm[j,i]=cum           j+=1          print "SVD results for seventeen filters:\n",res     print "SVD's covariance matrix:\n",cvm               sig=zeros(cvm.shape[0],float)     for i in range(cvm.shape[0]):       for j in range(cvm.shape[1]):         if i==j:           sig[i]=math.sqrt(cvm[i,j])          print 'Variance:\n',sig               chi_square=0     v=dot(coff,res)          for i in range(b.shape[0]):         chi_square += (b[i]-v[i])**2               print "chi_square:\n",chi_square          reduce_chi=chi_square/(coff.shape[0]-coff.shape[1]-1)     print "Reduced-Chisquare:\n",reduce_chi      Well my approach is not the optimize approach but I need to see for instance what is the value for `Reduced-Chisquare` or `covariance` but it raises **singular matrix error** when I try to inverse `S` matrix but if I use the following procedure:               Least_squares,residuals,rank,Singular_values=np.linalg.lstsq(coff, b)      It doesn't give me any error and compute the regression problem. **My questions** : _First_ : why this problem occurs using `SVD` (what are the differences between SVD and least square in this case)? _Second_ (very programming question): How could I keep the first approach and use the second one just in case that one error raises in the first method?


title: An association rules algorithm that maintains the order of items

text: For example: if my dataset contains (A, B), but does not contain (B, A). Then the algorithm may generate the rule A -> B, but will not generate the rule B -> A. Is there an association rules algorithm, or variation of an algorithm, that respects the order of items in this manner?


title: Is there a standard notation for nested factors?

text: What is the standard notation for specifying nested factors, for example, if I am estimating a parameter $\large{\beta{{_\text{site}}_i}}$ for the effect of $i=1\ldots m$ sites, and $\large{\beta{{_\text{tr}}_j}}$ for $j=1\ldots n$ treatments, assuming that there are no treatments that occur at each site (except the control, in which case $\beta{_\text{tr}}=0$, how would I note that any given treatment $j$ is associated with one and only one site? Here is one representation of my model: $\large{\Theta_{k} = \beta{_0} + \beta_{\text{site}_{i}} + \beta_{\text{tr}_{ij}}+ \epsilon}$ But I am not sure that simply indexing $ij$ gets this point across, or does it? Some sources seem to use $\large{\beta_{\text{tr}_{ij}}}$ while others use $\large{\beta_{\text{tr}_{j(i)}}}$.


title: ANCOVA (GLM) versus Mixed Linear Model with covariate measured at both pre- and post-treatment (in SPSS)

text: I'm trying to analyze a pre/post study with 4 treatment groups. The Dependent Variable, **Speed** , measures how quickly a person's facial expression reaches peak intensity; this measure is normally distributed after log10 transformation. Another variable, **Intensity** of emotional expression, is also measured at both pre and post and might influence the **Speed**. **Intensity** is continuous and not normally distributed. **Question:** If Intensity does indeed significantly correlate with Speed, how should this be accounted for, statistically? IF Intensity was "fixed"/static, I would think that (using SPSS), I would choose the repeated-measures option of the General Linear Model, and enter Intensity a covariate. However, since Intensity is measured--and not identical--at both pre- and post-, can Intensity be included as a covariate? If so, how (Enter 2 covariates: Intensity at pre, and Intensity at post?; Create a change-score to use as the covariate?). In one of the treatment conditions, Intensity is significantly different at post-treatment. **Update:** I've been advised that mixed-effects modelling (mixed linear model) might be the way to go. After a quick look, the mixed approach seems appropriate, but please advise if you think that there are better/easier options. My stats and SPSS skills are fairly basic so, if explanations can be calibrated to novice-level, it'd be much appreciated.


title: Variance of exponential distribution when changing units

text: I was tutoring a Stat student taking a probability course and I came across a problem dealing with variance of exponential that has left me confused. $X \sim \mathrm{Exp}(\mathrm{rate}=\lambda)$ so $E(X)=1/\lambda$ and $\newcommand{\Var}{\mathrm{Var}}\Var(X)=1/\lambda^2$ Customers enter a store at a rate of 2 per hour, a customer has just arrived, what is the variability of the amount of time until the next customer arrives in minutes? $\lambda=2, \; E(X)=1/2, \; \Var(X)=1/2^2=1/4$ Here the expected value = 30 minutes (1/2 of an hour) and the variance = 15 minutes (1/4 of an hour) Now work on minutes scale: $\lambda=2/60=1/30, \; E(X)=30, \; \Var(X)=30^2=900$ Here the expected value = 30 minutes but the variance = 900 minutes I am confused why I am not getting the same variance when I switch from minutes to hour and vice versa.


title: Best Multiple Imputation Method for Multiple Surveys Mixed Together, Presented Randomly?

text: I am working with a dataset containing data from 15 different surveys. The surveys were presented all as one battery to participants, with questions from all surveys essentially placed into a pool and then presented to participants at one of 30+ 'random' assortments. Most of these surveys measure the DVs in which we're interested directly, but a good handful of the questions are demographic questions. The scales' scores cannot be computed in a manner that tolerates even one missing value. Thus, you can imagine there are a good deal of missing scale scores because of (almost always) one missing item score. In order to apply multiple imputation to this dataset, should we conduct MI to each individual survey (i.e., at the item-level) and then use the scale scores altogether, or conduct MI at the scale-level (i.e., with every different question from every different survey included in the model)? My impression is that item-level makes most sense for a number of reasons, not the least of which is the retention of a good deal of power. However, my colleagues argue that doing so ignores the entire covariance structure and therefore produces biased results. I understand their argument, but I don't think that sort of brute-force application to a number of unrelated surveys is the best approach, either. Does anyone have experience with this issue? Many thanks. (P.S. Solutions in SPSS and/or R are best understood, beyond mathematical models or plain-English descriptions.)


title: Is this normal to have big value of Chi Squared?

text: I am using chi squared for feature selection in text classification. However when I compute it I sometime have very big values. Like 100, 1000 or even 20000. Is this normal ? I wonder because I read that if the Chi Squared score equals 10.83 then the p-value is 0.01. So we know at 99.99% that the token and the class of a text document are dependant. What do this mean for 20000 ? Are we sure at 99.9999999 % or something like that ? * * * **---- EDIT 6/23/2014 ---** * * * Here is more information about what I am trying to do with chi-square. We have 64893 documents. For each document we tell if it belongs to a category (in this case we write a +1 next to it) or if it does not belong to it (in which case we write a -1 next to it). For each word in the corpus of document, we want to determine if it is a word which is helpful to put the document in the category. The NULL hypothesis is "document category has no influence over word frequency". The degree of freedom is 1. Here are a few example of the big results I receive: The first one have a Chi-Score of 32250 ![enter image description here](http://i.stack.imgur.com/stAsg.png) * * * The second one have a Chi-Score of 11410 ![enter image description here](http://i.stack.imgur.com/Lv6dk.png) **---- EDIT 6/24/2014 ---** Imagine I have the following documents:               I am very happy      I was happy yesterday     I am kind of sad     I am happy today     I am sad     they are happy now     yesterday, he was not happy     today he is happy     they will be happy tomorrow      For each document, we manually identify documents talking about someone happy **right now** by giving them the class +1               +1   I am very happy      -1   I was happy yesterday     -1   I am kind of sad     +1   I am happy today     -1   I am sad      +1   they are happy now     -1   yesterday, he was not happy     +1   today he is happy     -1   they will be happy tomorrow     +1   happy day today ! I am happy      For each document, we gather the statistics for each word like in the example below. _(Here I am not sure If I should put a 2 for the "happy" word because it appears twice or just a 1 because it appears)_               //  happy day today ! I am happy          category = "+1",           happy                  1     sad                    0     today                  1     kind                   0     of                     1     not                    0     will                   0     be                     1     tomorrow               0     very                   0     are                    0     now                    0     is                     0     day                    1     I                      1     was                    0     yesterday              0     they                   0     he                     0     am                     1        Now we want to determine for each word, how much they help to determine if a document belongs to the class. We compute the Chi-Square score for each word:               happy                               3.75     sad                                 3.75     today                   2.85714285714286     kind                    1.66666666666667     of                      1.66666666666667     not                     1.66666666666667     will                    1.66666666666667     be                      1.66666666666667     tomorrow                1.66666666666667     very                   0.740740740740741     are                    0.740740740740741     now                    0.740740740740741     is                     0.740740740740741     day                    0.740740740740741     I                      0.277777777777778     was                    0.104166666666667     yesterday              0.104166666666667     they                   0.104166666666667     he                     0.104166666666667     am                                     0      _If I put a "2" instead of a "1" for the happy word in the document "happy day today! I am happy" I obtain_               happy                   11.8518518518519     sad                                 3.75     today                   2.85714285714286     kind                    1.66666666666667     of                      1.66666666666667     .     .     .


title: What is the difference between a zero-inflated and a zero-truncated poisson?

text: I'm trying to make sense of a question which uses a zero-inflated poisson model given by: $$ f(x; \lambda,\omega) = \begin{cases} \omega + (1-\omega)e^{-\lambda} &\mbox{if } x = 0 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (1) \\ \frac{(1-\omega)e^{-\lambda}\lambda^x}{x!} & \mbox{if } x = 1,2,3,\dots \ \ \ \ \ \ \ (2)\end{cases} $$ We are given a table of data, x = 0, 1, 2, 3, 4 and the number of occurrences of each. In lectures we covered a very similar question, however it used a zero- truncated poisson, given by: $$ g(x; \theta) = \frac{e^{-\theta}\theta^x}{x!(1-e^{-\theta})} x = 1,2,3,\dots \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (3) $$ Obviously if I equate $(2)$ with $(3)$ I get $(1-\omega) = (1-e^{-\theta})^{-1}$. Now the question asks me to find the log-likelihood functions for $\lambda$ and $\omega$, then show that the MLE for $\hat{\lambda}$ is given by $$\frac{\hat{\lambda}}{1-e^{-\hat{\lambda}}} = \bar{x}$$ Where $\bar{x}$ is the expected mean for $x_i \neq 0$.) Now, this is the exact answer we got in lectures - and it's the answer I get here if I ignore (1) in my calculations. If I try to include both (1) and (2) in calculating the MLE, things get pretty messy. My question is, what happens to (1)? Is there actually a difference between these two distributions? Also, once I calculate $\lambda$ and $\omega$, and use (1) to get a value for the frequency of $x=0$, I get a value very close to zero, whereas the frequency in the original table is 97 (all the other frequencies calculated by the model tally closely, as expected). One other thing, if you are being kind enough to consider answering my question, is that I have no idea what a link-function is, and presume this question can be answered without resorting to such functions.


title: Taking into account Bayesian model uncertainty

text: I recently received a review of a paper from a Bayesian Statistics Journal. The Associate Editor wrote this mini-review (quoted below in full). > The paper is talking about Bayesian modeling of DNA sequences. However, the > model uncertainty is not taken into account in Bayesian inference; the > author used simulated annealing to search for a MAP model and then made > inference based on it. This is a big flaw of the paper and should be > corrected in the revision. The problem is that I have no idea what the reviewer has in mind exactly, or how one would go about doing this. I guess it is difficult for people here to give any precise advice about this issue without seeing the paper, but I would be happy to send anyone who is interested a draft. Regardless, the issue is that I'm using simulated annealing in the space of models (a very large but finite space) and when I have found a good candidate, I use that model, and only that model, in subsequent analyses. Obviously there is no very compelling reason for picking exactly this model; I could easily have picked a similar model. I think the reviewer wants me to integrate the uncertainty in the model into the analyses. I don't know exactly what this would involve, but I suspect this would pose substantial difficulty, both on the more theoretical side, as well as on the implementation side. The problem is already quite complicated computationally. My immediate reaction is that to do this, one would likely need a probability distributions on the space of models, and no such distribution is available, nor could I imagine how one could find such a distribution. I'd like a good reference for model uncertainty, preferably including an illustration of how one would take model uncertainty into account in Bayesian inference. Therefore, I'm writing to ask if someone could provide me with such a reference (or references).


title: Looking for impl of a Fuzzy Neural Network (FNN or SOFNN)

text: I'm looking for an implementation of FNN (or better yet, a SOFNN as described by Forecasting Time Series by SOFNN with Reinforcement Learning). Any language, though preference is Java, C#, C++ in that order.


title: What is a meaning of "p-value F" from Friedman test?

text: I'm a new R user and had just tried running friedman on non-normal and heteroscedastic data on seagrass. I am testing whether biomass is significantly different between sites across years. R (`friedman` function from `agricolae` package) returned result like this:               Friedman's Test     ===============     Adjusted for ties     Value: 0.01333333     Pvalue chisq : 0.9080726     F value : 0.01316678     Pvalue F: 0.9086942 0.9087002          Alpha     : 0.05     t-Student : 1.990847     LSD       : 17.34995          Means with the same letter are not significantly different.     GroupTreatment and Sum of the ranks     a   1   119      a   2   118       I know this means no significant difference given the p-value chisq. But what does p-value F mean?


title: If variable kernel widths are often good for kernel regression, why are they generally not good for kernel density estimation?

text: This question is prompted by discussion elsewhere. Variable kernels are often used in local regression. For example, loess is widely used and works well as a regression smoother, and is based on a kernel of variable width that adapts to data sparsity. On the other hand, variable kernels are usually thought to lead to poor estimators in kernel density estimation (see Terrell and Scott, 1992). Is there an intuitive reason why they would work well for regression but not for density estimation?


title: How to read large dataset in R

text: trying to read a large dataset in R the console displayed the follwing errors:               data<-read.csv("UserDailyStats.csv", sep=",", header=T, na.strings="-", stringsAsFactors=FALSE)     > data = data[complete.cases(data),]     > dataset<-data.frame(user_id=as.character(data[,1]),event_date= as.character(data[,2]),day_of_week=as.factor(data[,3]),distinct_events_a_count=as.numeric(as.character(data[,4])),total_events_a_count=as.numeric(as.character(data[,5])),events_a_duration=as.numeric(as.character(data[,6])),distinct_events_b_count=as.numeric(as.character(data[,7])),total_events_b=as.numeric(as.character(data[,8])),events_b_duration= as.numeric(as.character(data[,9])))     Error: cannot allocate vector of size 94.3 Mb     In addition: Warning messages:     1: In data.frame(user_msisdn = as.character(data[, 1]), calls_date = as.character(data[,  :       NAs introduced by coercion     2: In data.frame(user_msisdn = as.character(data[, 1]), calls_date = as.character(data[,  :       NAs introduced by coercion     3: In class(value) <- "data.frame" :       Reached total allocation of 3583Mb: see help(memory.size)     4: In class(value) <- "data.frame" :       Reached total allocation of 3583Mb: see help(memory.size)      Does anyone know how to read large datasets? The size of UserDailyStats.csv is approximately 2GB.


title: Random effect slopes in linear mixed models

text: In my data, the RT (gaze) of individuals (ID) is examined as a function of a visual conditions, the factor size (small, medium, large). Base model:               print(Base <- lmer(RT ~ Size + (1|ID), data=rt), cor=F)      Random effect:               print(NoCor <- lmer(RT ~ Size + (0+Size|ID) , data=rt))     print(WithCor <- lmer(RT ~ Size + (1+Size|ID), data=rt))      Addition of ID slopes improves the Base model. My question is, how can a significant random effect (Size/ID) be interpreted when there is no relationship between the random and fixed effect, i.e., when the correlation between the random factor and the fixed facor does not improve the model [the anova(NoCr, WithCor) does not show a significant improvement]?


title: How to reconstruct ensemble of trees from random forest?

text: I have a great prediction yet I am unsure how to uncover how the results were generated?


title: Weighting related attributes in hierarchical clustering

text: I have two questions about output of hierarchical clustering and improving the output. I'm trying to learn more about performing hierarchical clustering in R so I started looking at a simple dataset I created of sushi rolls at a local restaurant. I went though every roll on the menu and created a distinct list of the union of all ingredients. Then for each roll I put a 1 if it had that particular ingredient or a 0 for not having it. I then calculated distance based on Jaccard distance and created a dendogram using four different methods. So my first question is how to interpret these correctly. Since many rolls overlap on quite a bit of ingredients, single-linkage is producing clusters that aren't significantly different from each other? Complete-linkage and the two methods are better reflecting the diversity in the rolls where some are significantly different? **I'm not sure I really understand the subtleties between these four methods**. ![enter image description here](http://i.stack.imgur.com/f2dcn.png) My second question is **how would one handle related attributes**? Say for example there exists {Tuna, White Tuna, Spicy Tuna} as three different ingredients. I think we would agree that the distance between {Tuna,Avocado} & {Salmon, Avocado} should be larger than {Tuna, Avocado} & {Spicy Tuna, Avocado}. Are the attributes typically collapsed into {Tuna} or is there another way to reflect the relationship?


title: Logistic Regression - Two Dependent variables not playing well with each other

text: My name is Abhi and I trying to teach myself logistic regression by solving some of the problems available on the internet. I am using R and RStudio as the development environment **_Problem Statement_**   Given the age, sex, class(first,second,third), ticket_id for each passenger can you predict if he survived or died when the titanic sank. Assume passengers in first class are more likely to survive. Assume passengers in certain areas of the boat (those near to part which hit the iceberg) are more likely to perish. I ran logistic regression for both these variables independently. **Here are my results for passenger class**                 Deviance Residuals:          Min       1Q   Median       3Q      Max       -1.4390  -0.7569  -0.7569   0.9367   1.6673            Coefficients:                 Estimate Std. Error z value Pr(>|z|)         (Intercept)  1.44679    0.20743   6.975 3.06e-12 ***     Pclass      -0.85011    0.08715  -9.755  < 2e-16 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          (Dispersion parameter for binomial family taken to be 1)              Null deviance: 1186.7  on 890  degrees of freedom     Residual deviance: 1084.4  on 889  degrees of freedom     AIC: 1088.4      **Here are my results for Ticket Id**                 Deviance Residuals:          Min       1Q   Median       3Q      Max       -1.2624  -0.8909  -0.7040   1.0947   1.7415            Coefficients:                 Estimate Std. Error z value Pr(>|z|)         (Intercept)  0.19769    0.09974   1.982   0.0475 *       Ticket2     -0.91681    0.22155  -4.138  3.5e-05 ***     Ticket3     -1.46620    0.16052  -9.134  < 2e-16 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          (Dispersion parameter for binomial family taken to be 1)              Null deviance: 1186.7  on 890  degrees of freedom     Residual deviance: 1094.0  on 888  degrees of freedom     AIC: 1100          Number of Fisher Scoring iterations: 4      However when I add both passenger class & ticket id to the formula I dont get the results that I am expecting. It seems that both the features are cancelling each other out                 Deviance Residuals:          Min       1Q   Median       3Q      Max       -1.4431  -0.8649  -0.6772   0.9332   1.7806            Coefficients:                 Estimate Std. Error z value Pr(>|z|)         (Intercept)   1.1713     0.2171   5.396 6.80e-08 ***     Pclass       -0.5654     0.1110  -5.094 3.50e-07 ***     Ticket2      -0.4361     0.2421  -1.801   0.0716 .       Ticket3      -0.8311     0.2033  -4.089 4.33e-05 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          (Dispersion parameter for binomial family taken to be 1)              Null deviance: 1186.7  on 890  degrees of freedom     Residual deviance: 1067.8  on 887  degrees of freedom     AIC: 1075.8      I ran the vif test on both these features & they do not seem to be collinear. What am I missing here? If both these features give great results separately why dont they give good results when used together. Any help would be much appreciated


title: Standard error of the sampling distribution of the mean

text: I found an equation that says the standard error of the sampling distribution of the mean is: $$\sigma_{\bar{X}} = \sigma \cdot \sqrt{\frac{1}{n}-\frac{1}{N}}$$ And when the population size is very large, the factor $1/N$ is approximately equal to zero; and the standard deviation formula reduces to: $$\sigma_{\bar{X}} = \frac{\sigma}{\sqrt{n}}$$ However, I did not get how he found the first equation, may someone explain it? I found this equation here: http://stattrek.com/sampling/sampling- distribution.aspx I could explain the second equation in that way: \begin{equation} \begin{split} Var(\bar{X}) = Var\left(\frac{\sum_{i=1}^{n} X_i}{n}\right) \\ = \left(\frac{1}{n^2}\right)Var\left(\sum_{i=1}^{n} X_i\right) \\ = \left(\frac{1}{n^2}\right)\left(\sum_{i=1}^{n} Var(X_i)\right) \\ = \left(\frac{1}{n^2}\right)\left(\sum_{i=1}^{n}\sigma^2 \right) \\ = \left(\frac{1}{n^2}\right)\left(n\sigma^2\right) \\ = \frac{\sigma^2}{n} \end{split} \end{equation} However I don't know from where appeared the first equation.


title: Precedent for Bootstrap-like procedure with "invented" data?

text: Let's say I have a dataset with 1000 observations in 10 variables, "A" through "J." I have 1000 responses/measures for each of the first 8 variables, through "H," but only the first 500 observations for "I" are not missing, and only the last 500 observations for "J" are not missing -- there are no observations for which I have measures of both of the last two variables, I and J. Thus, if I calculate (pairwise) correlations, I have a full correlation matrix, with only the correlation between I and J missing. Let's say I want to run Principal Component Analysis, or some other such scaling procedure on this correlation matrix. What I think I would like to do is:   1. Randomly generate (perhaps from some distribution on [-1, 1], or perhaps via sampling from existing values in the rest of the correlation matrix) an "invented" correlation between I and J.   2. Put that in the correlation matrix.   3. Run PCA on the correlation matrix with this invented value.   4. Repeat steps 1 - 3 some large number of times.   5. Assess the collective results of this large number of PCAs, looking at the mean and variance of the loadings, scores, eigenvalues, etc., based on the "pseudo-bootstrapped" iterations. Questions:   1. Is there a better way to handle (a) missing value(s) in the correlation matrix?   2. Is there any precedent for replacing such (a) missing value(s) with random invented values? If so, what is it called?   3. Is this related to the bootstrap? Thanks a lot, in advance. Edit: Question 4. Is this a defensible approach to imputation?


title: Find range of possible values for probability of intersection given individual probabilities

text: If given probability of $A$ is $a$ and probability of $B$ is $b$, how do I find min/max probability of intersection? Max value of intersection would be $\min(a,b)$, how do I find the min?


title: Steps to perform time series analysis

text: I'm trying estimate an autoregressive model with an exogenous variable. It's about the impact of changes in oil prices on the economy. I'm planning on regressing gdp growth rate on its own lagged values and lagged values of oil price. I don't know from where to start. Here's what I think so far:   1. I should check if the gdp and oil price are stationary, right? (dickey-fuller test??)    2. Find the number of lags for gdp and oil price.( AIC criterion?.)    3. Estimate the model.    4. Test the signification of the estimators, heteroscadesticity and then serial correlation.  Should I test for the presence of cointegration? later I'm planning on using the same variables in a Markov switching model. I have the same question for Grangers causality test. I've been reading a few books and watching more than one video about times series, but I still can't find my way through it. Can you please guide me and see if there is something missing in those steps?


title: MATLAB - Maximum likelihood for periodic function (angles)

text: In the MATLAB stats tutorials there is a section called _"Fitting a More Complicated Distribution: A Mixture of Two Normals"_ http://www.mathworks.com/help/stats/examples/fitting-custom-univariate- distributions.html               pdf_normmixture = @(x,p,mu1,mu2,sigma1,sigma2) ...                          p*normpdf(x,mu1,sigma1) + (1-p)*normpdf(x,mu2,sigma2);     lb = [0 -Inf -Inf 0 0];     ub = [1 Inf Inf Inf Inf];     start = [pStart muStart sigmaStart sigmaStart];     paramEsts = mle(x, 'pdf',pdf_normmixture, 'start',start, 'lower',lb, 'upper',ub)      It uses MLE to estimate numerically the values " [...] = MLE(DATA,'pdf',PDF,'cdf',CDF,'start',START,...) returns MLEs for the parameters of the distribution defined by the probability density and cumulative distribution functions PDF and CDF. PDF and CDF are function handles created using @. They accept as inputs a vector of data and one or more individual distribution parameters, and return vectors of probability density values and cumulative probability values, respectively. If the 'censoring' name/value pair is not present, you may omit the 'cdf' name/value pair. MLE computes the estimates by numerically maximizing the distribution's log-likelihood, and START is a vector containing initial values for the parameters." I would like to apply the same methodology for fitting two or more normals to a univariate set of values that I have, but **within a periodic domain**. That is, angles that have values of 0° to 360° linked together as a circular range. I am not sure how to declare it in order to make MATLAB understand this kind of terminology. Would it be possible to change this implementation to add the circular range case? Regards, Ignacio


title: PCA prediction with centered vs. un-centered input data

text: I'm dealing spectral data. I want to classify 2 classes using principle components analysis. My PCA was constructed using `prcomp(data, center=TRUE)` in R. It works. My questions are:   1. When using unknown data (input) to predict a classification, should the data be centered before applying PCA?    2. PCA yielded two results which were a little bit different between centered input data and un-centered input data. If all of input data were un-centering (same condition) before predicting, is it fair to classify them?


title: How to determine if the data points are linearly separable from an SVM hyperplane

text: 1. How to know the data points are linearly separable from an SVM hyperplane?   2. How to get the optimal classifier during iteration process?   3. How to calculate the complexity of the SVM model?


title: How to generate the most common clickstream sequence

text: I have logs with the following information: date-time username view action action_data These logs are generated from a web-application which consists of several views where the users can perform a variety of actions which are logged (which elements they click, if they hover somewhere, switch views etc). My desired output would be one or several clickstreams with some diverging paths which describe the most common behaviour (this will in turn be used as input for which sections of the site to test). Since this will be used for testing, the generated clickstream doesn’t need to be definitive/straight but can include some diverging paths (40% of users click button A here and 60% press button B). I’ve been looking at Markov chains and associative analysis for solving this but can’t land in something definitive about which methods to use.


title: Truncate lognormal distribution with excel

text: I have a dataset where I extract the mean and the standard deviation, I want to generate a new synthetic dataset which is lognormal distributed based on the original dataset paramaters. Like a montecarlo simulation. My problem is that I need that synthetic dataset is truncated, with the maximum and minimum value of the orignal set, in other words I don't want any synthetic value greater than the maximum value and minimum value of the original datset. Thanks @whuber, following your guidelines it works perfectly with the positive side of my original dataset, long right tail distribution within the max and min limits. ![enter image description here](http://i.stack.imgur.com/63ddE.jpg)


title: Difference between cumulants and moments

text: In particular, is the $n$th cumulant equivalent to the $n$th central moment (i.e. about the mean)? There's little difference I can see between MGFs (moment generating) and CGFs (cumulant generating), apart from the former gives moments about the origin while the latter yields central moments.


title: How should I code these two variables?

text: I have two independent variables: "vaccination status" and "antibody titer assessment". Vaccination is coded as 0 and 1. The other variable is supposed to explain if the person has checked that whether the vaccination has increased his antibody titer sufficiently. The second variable has two states: yes or no. The problem is that people who are not vaccinated do not need to check their antibody titer (it is not logical). So I don't know if I should leave the cells in the column pertaining to tier check blank or not, when the persons' vaccination status is zero. Or should I enter zero even for the lack of antibody titer assessment, even if the person has not received any vaccination in the first place? The difference between the two is that if I fill those blank cells with zero, these two independent variables would become highly correlated (source of multicollinearity). Otherwise, if i leave those cells blank, I would loose data in listwise deletion unless I use other methods such as pairwise deletion or meansubstitution (in SPSS). What do you suggest? Do you have any other ideas that these two options? Thanks in advance. :)


title: Publications discussing "checking to see if outliers affect the result"

text: I have noticed it is somewhat common for people to collect data, perform a statistical test (e.g., t-test) and then also look for possible "outliers" in the data. Then they remove the outliers and also perform the test on the filtered data. If it is significant in both cases it is concluded that "the outliers don't affect the result". If it is significant only in the second case then usually this is noted and further analysis/discussion is done using only the "non- outliers". This seems to be considered acceptable as long as it is reported that you did it. I am wondering if anyone has thought deeply about this and discussed the logic of doing so. Also, has anyone investigated this process to determine the consequences for false-positive rate and power function. **Edit:** To clarify, when I am thinking of "outliers" these are data points that appear different from the rest but there is not necessarily a theoretical reason to treat them differently. The data itself _suggests_ they are different in some unknown way. In that way this question may differ from previous questions on the site. As an example of discussing the logic, Roger Ratcliff makes a well-reasoned argument in the case of reaction-time data: > The outliers' distribution most likely has a much greater standard deviation > than the distribution of the real process, giving a long and elevated tail, > because **outliers often come from spurious events such as loss of > attention, daydreaming, distraction, and so forth.** Determining whether any > particular group of extreme reaction times contains mostly real responses or > mostly outliers is extremely difficult. Therefore, finding methods that > minimize the effects of suspect observations is an important aim of > statistical methods in the analysis of reaction time. In fact, **the goal > for our models and empirical research should be to account for the middle > 85-95% of the observations in our reaction time distributions; these are the > data that are most likely to come from the real processes under > consideration and also most likely to be critical in testing hypotheses and > models.** Psychol Bull. 1993 Nov;114(3):510-32. Methods for dealing with reaction time outliers. Ratcliff R. (pg 2 of paper) I would agree with his assessment. But what if you do not have any a priori reason to suspect your data will contain outliers other than "individual variability"? **Edit 3:** Something I thought of regarding "the outliers don't affect the result". Clearly the inclusion/exclusion of outliers will affect any parameter estimation even if the p-value is still on the same side of the significance threshold. So in that sense the claim is false. Therefore it is not something that should be taken lightly (described in one sentence) because it will impede meta-analysis if different researchers use different definitions of outlier. Also, the exact values of the outliers would be of interest to future researchers since they may observe a similar phenomenon caused by the same process.


title: Writing the objective function and constraints for scipy.optimize.minimize from matrices

text: I'm trying to find the optimal threshold point from a ROC curve. I have two main constraints : tpr >=80 and fpr <=60. I tried three main minimization functions :               #where, tpr=True positive rate and fpr = false +ve rate           dis= math.sqrt(math.pow((1-arr[1]), 2) + math.pow((0-arr[0]), 2))     dis= math.pow((1-arr[1]), 2) - math.pow((0-arr[0]), 2)     dis= -1*math.sqrt(math.pow((1-arr[1]), 2) + math.pow((0-arr[0]), 2))      My first function assumes that the point with least distance from best point (1, 0) while meeting my constraints is the optimal point. Is that a correct way to think ? I find the point (0.84558824, 0.58104396) to be more optimal as it gives me better tpr with a higher fpr within my constraints but my function returns a value of x = 0.53708791, y = 0.80147059 as the optimal point. ![](http://i.imgur.com/TdqMPek.png) I then tried weighting the trp by a factor of 6. The reason for this is that I want to weight the "tpr" more than "fpr". This gave me (0.84558824, 0.58104396) which I think is the best point. dis = math.sqrt(((math.pow((0-arr[1]), 2))) + (6 * math.pow((1-arr[0]), 2))) The new graph : ![](http://imgur.com/7DwQkSL.png) There was a 3rd approach I tried, where I took, 1`*math.sqrt(math.pow((1-arr[1]), 2) + math.pow((0-arr[0]), 2))`, since I want to distance between `1 & tpr` to be low, but that between 0 and fpr to be small. ie, the best point lies left most & upward. **What is a more mathematical way to approach the problem ?** I also saw the minimize function with ability to specify constraints with bounds, but I am unable to formulate the problem . The three possibility for the mathematical formulation above, all look intuitive, but give different results wrt to the optimal point. To put things together, I used a pseudo heuristic as here               def min_fun(arr):          dis = math.sqrt(((math.pow((0-arr[0]), 2))) + (6 * math.pow((1-arr[1]), 2)))          if arr[1] > 0.80 and arr[0] < 0.60:               return dis           else:                #return max distance for ROC graph if conditions not met.                return math.sqrt(2)     #and then, just appending the list of minimum distances (candidate set), to get the best point later      a_list = zip(fpr,tpr,thresholds)     distances = []     for s in a_list:         distances.append(min_fun(s))     threshold_index = distances.index(min(distances))     threshold  =a_list[threshold_index][2]      I know grid search/DP/Scipy minimize are the right candiadtes here - to which, how do I frame this problem rigorously ? Related -   1. Python minimize function: passing additional arguments to constraint dictonary   2. Dynamically writing the objective function and constraints for scipy.optimize.minimize from matrices


title: Fishy significance test: Am I doing something wrong with ANOVAs?

text: Just when I thought I'd had a grip on how to do an analysis of variance this particular data set had me startled: it's a collection of response times (in ms) to a linguistic input. To be precise, it's part of a reading time experiment, and I'm trying to see if there's a significant effect of two factors. My experiment had a 2x2 factorial design, with 2 factors binary `Conflicting` and `ContextPresent`. The only value I'm interested in for the purpose of this question is one particular random variable, a response time. No transformations have been done on it, except removal of outliers via 3-sigma rule (I used mean + 3 * standard deviation to determine outliers.) So I run my anova in R:               > anova(lm(TextDisplay9.RT ~ Conflicting * ContextPresent, data=items.cropped))     Analysis of Variance Table          Response: TextDisplay9.RT                                 Df   Sum Sq Mean Sq F value  Pr(>F)     Conflicting                  1   111185  111185  7.0591 0.00808 **     ContextPresent               1    73591   73591  4.6723 0.03102 *     Conflicting:ContextPresent   1      352     352  0.0223 0.88128     Residuals                  651 10253667   15751      Jolly ho, I get pretty good results for a main effect on both my factors, and no interaction. That's fine. But here's the corresponding box-and-whiskers graph:               > ggplot(items.cropped,aes(Conflicting,TextDisplay9.RT)) + geom_boxplot() +       facet_grid(.~ ContextPresent)      ![Boxplot of suspiciously equal looking stuff](http://i.stack.imgur.com/R9d2V.png) So this plot actually makes it seem like there _shouldn't_ be a main effect of either variable! They're all too similar! Yes, the scale is rather squished, because of the outliers, but the means are really really close!               > with(items.cropped,mean(items.cropped[Conflicting=="semantic conflict" & ContextPresent == "mentioned in context",]$TextDisplay9.RT,na.rm=T))         [1] 431.8659         > with(items.cropped,mean(items.cropped[Conflicting=="semantic conflict" & ContextPresent == "not mentioned in context",]$TextDisplay9.RT,na.rm=T))     [1] 454.5305     > with(items.cropped,mean(items.cropped[Conflicting=="no semantic conflict" & ContextPresent == "mentioned in context",]$TextDisplay9.RT,na.rm=T))         [1] 407.485         > with(items.cropped,mean(items.cropped[Conflicting=="no semantic conflict" & ContextPresent == "not mentioned in context",]$TextDisplay9.RT,na.rm=T))     [1] 427.2188      Does it sound possible that there could be a main effect? Or did I somehow misuse ANOVAs? I could provide the data if needed! Thanks very much for any suggestions.


title: Standard deviation for a measurement that has two causes of uncertainty

text: I would like to know how to calculate the standard deviation for a measurement value $x$ when there are two separate sources of uncertainty, 1 and 2, each with known standard deviation, $\sigma_1$ and $\sigma_2$. I’ll give a couple of examples.   * I measure length with a ruler. Each measurement has an uncertainty ($\sigma_1$) introduced by simply using the tool. Additionally, I may know that the manufacturing process of such rulers isn’t perfect, so I could associate some standard deviation $\sigma_2$ introduced by the manufacturing process.   * I measure the the proportion of students in a class of size $N$ that answer a question correctly. There’s a standard deviation $\sigma_1$ associated with an assumed underlying binomial distribution. But I may also want to account for the variability in instructor effectiveness, which I choose model with a standard deviation $\sigma_2$. **Thus, my question:** Is there a way to combine $\sigma_1$ and $\sigma_2$ to get the standard deviation for such measurements? Or, what sort of reference book might be appropriate? This seems like a measurement idea. My gut tells me I can treat this like a propagation of uncertainty task where I am adding two numbers with known standard deviation and do something like $\sigma_\text{combined}=\left(\sigma_\text{binom}^2 + \sigma_\text{other}^2\right)^{1/2}$, but I can't seem to find a source of information on this (or I can't find the right terms to search).


title: Best way to estimate mean of RV from independent samples

text: Suppose we have N independently drawn samples from an unknown random variable $X$. What is the best way to estimate the expectation of $X$? For simplicity, we can assume that $X$ only returns values between 0 and 1. I know this might be a stupid question: the standard way to estimate the mean is just to just take the average of all the samples. But I got to thinking that this isn't optimal. For example, maybe it's better to take the median-of- averages? I.e. to split the samples into some number of equal-sized buckets, compute the average of each bucket, and then take the median of the averages? **Subquestion:** what is the best way to estimate the mean assuming $X$ is a normal random variable with standard deviation 1? Even for this goal I'm not sure which is better: median, average, or something else. **Update:** After thinking about @Nick's comment, I think there might be no "best" estimator for the mean. In particular, even for normal distributions, I think that taking the sample mean has advantages and disadvantages over taking the median-of-averages or just the median. Overall, I think median-based methods will give a more spread-out answer, but you'll have better certainty to be within some confidence interval than using mean, but I need to do the math to be sure of this. Maybe there are some good references on the topic?


title: Weibull distribution - is it “heavy tailed” with three parameters?

text: I am fitting my data to **Weibull** distribution (I would like to **prove** , that they follow a **power law** ). According to http://en.wikipedia.org/wiki/Heavy-tailed_distribution (and some articles), if "shape parameter" in Weibull is smaller then 1, I should have a heavy tailed distribution. If I use Weibull with two parameters to my data, I get "shape parameter" **greater** than **1** - so no heavy tail :-( But - with **three parameters-Weibull** I get "shape parameter" always **smaller** than **1** (heavy tail here). My question is: " **Do my data have a heavy tail, if I use Weibull with three parameters and shape parameter is smaller than 1**?"


title: Singularity issues in R when using mlogit package

text: I tried to use `mlogit` for discrete choice modeling, but I found that I cannot estimate a model for the alternative specific variables, which are defined only for one alternative--for example, in a dataset like below:               Num   Mode   Ftp     1      Bus    0     1      Walk   0     1      Car    1.0     2      Bus    0     2      Walk   0     2      Car    0.5      I have tried different shapes and different arrangements in the formula, but none of them works. Any thought?


title: R script for split plot augmented randomized complete block experiment design

text: I’m working on a data set in order to evaluate the impact of water stress on Rils population. The experimental design is split plot augmented design and it's as follows: The first factor (stressed/not stressed) is in Whole plot and there are no replicates. The second factor (genotypes) is in split plot and compound 152 lines, like way we have 8 checks repeated in each of 14 blocks and completed by 144 lines dispatched on 14 blocks so we have 8C + 10N in each block. The question is what is the model of Anova that explains this trial and how can I use R to calculate lsmeans.


title: Exploratory regression analysis for data with missing values

text: Recently I have performed an **exploratory regression analysis** , using `lavaan` R package and observed the following _output_ with some **warning messages** in it. I have the following **question** : What is the **potential negative impact** , if any, of these warnings on the _validity_ of the analysis and what can be done to **improve** the validity? Some people have already asked me why would I use essentially a SEM package for regression analysis and not simply `lm()` function. I'd like to clarify the logic behind my choice of using `lavaan` as follows (please feel free to express your opinion on my arguments/assumptions/facts here): A) While SEM is what I ultimately will use in a later phase of my study ( _PLS-SEM_ , not _CB-SEM_ , but this is a different topic), I simply use `lavaan` 1) to get familiar with it and 2) to test/use its FIML functionality, as my datasets contain **significant ratio** of _missing values_ (Determining probability distribution for datasets with missing values). I've already used _multiple imputation_ in another module, but wanted to use _FIML_ , as well. I do have _latent variables (LVs)_ , but I use SEM here in **exploratory fashion** - the above-mentioned model's components represent the LVs' _indicators_. B) Another reason to prefer `lavaan` to `lm()` is that `lm()` uses only _OLS_ estimation, whereas `lavaan` offers a **wide selection** of _estimators_. This is additionally attractive, since my research study's datasets do not conform to normal distribution.               Warning messages:     1: In lav_data_full(data = data, group = group, group.label = group.label,  :       lavaan WARNING: some observed variances are (at least) a factor 1000 times larger than others; use varTable(fit) to investigate     2: In lav_data_full(data = data, group = group, group.label = group.label,  :       lavaan WARNING: some cases are empty and will be removed:       1672 1673 1761 1762 1821 1866 1871 1892 2141 2167 2200 2205 2223 2228 2234 2276 2283 2322 2585 2588 2589 2590 2591 2592 2593 2602 2603 2606 2607 2608 2609 2918 2923 2924 3448 3451 3468 3469 3498 3499 4866 4867 4868 4869 4870 4881 4882 4889 4964 4965 4966 4967 4968 4969 4970 4971 4972 4973 4974 4975 4976 4977 5003 5004 5005 5135 5136 5137 5138 5139 5161 5162 5190 5194 5195 5205 5206 5207 5208 5209 5210 5211 5212 5213 5214 5215 5216 5223 5224 5225 5226 5227 5228 5229 5230 5231 5232 5233 5234 5235 5236 5237 5238 5239 5240 5241 5242 5243 5244 5245 5246 5247 5251 5252 5253 5254 5255 5256 5257 5258 5259 5260 5261 5262 5263 5264 5265 5266 5267 5268 5269 5270 5271 5272 5273 5274 5275 5276 5277 5278 5279 5280 5281 5282 5283 5284 5285 5295 5296 5300 5301 5302 5309 5310 5313 5314 5315 5316 5317 5318 5322 5326 5330 5331 5332 5333 5334 5335 5336 5337 5338 5339 5340 5341 5342 5343 5344 5345 5346 5347 5348 5349 5353 5354 5355 5356 5357 5358  [... truncated]     3: In lav_data_full(data = data, group = group, group.label = group.label,  :       lavaan WARNING: due to missing values, some pairwise combinations have less than 10% coverage     4: In lav_model_vcov(lavmodel = lavmodel, lavsamplestats = lavsamplestats,  :       lavaan WARNING: could not compute standard errors!          Model 4 - Analysis results:          lavaan (0.5-16) converged normally after 174 iterations                                                            Used       Total       Number of observations                         94146      100000            Number of missing patterns                        23            Estimator                                         ML       Minimum Function Test Statistic             5252.040       Degrees of freedom                                 5       P-value (Chi-square)                           0.000          Model test baseline model:            Minimum Function Test Statistic             8266.259       Degrees of freedom                                15       P-value                                        0.000          User model versus baseline model:            Comparative Fit Index (CFI)                    0.364       Tucker-Lewis Index (TLI)                      -0.908          Loglikelihood and Information Criteria:            Loglikelihood user model (H0)             -1449561.438       Loglikelihood unrestricted model (H1)     -1446935.418            Number of free parameters                         22       Akaike (AIC)                              2899166.876       Bayesian (BIC)                            2899374.833       Sample-size adjusted Bayesian (BIC)       2899304.916          Root Mean Square Error of Approximation:            RMSEA                                          0.106       90 Percent Confidence Interval          0.103  0.108       P-value RMSEA <= 0.05                          0.000          Standardized Root Mean Square Residual:            SRMR                                           0.154          Parameter estimates:            Information                                 Observed       Standard Errors                             Standard                             Estimate  Std.err  Z-value  P(>|z|)   Std.lv  Std.all     Covariances:       team.size ~~         prj.age          42.402                              42.402    0.148         license         -26.689                             -26.689   -0.386         restrict          2.529                               2.529    0.208         maturity          1.310                               1.310    0.138       prj.age ~~         license           3.123                               3.123    0.026         restrict         -8.233                              -8.233   -0.395         maturity          2.477                               2.477    0.152       license ~~         restrict          3.210                               3.210    0.639         maturity          0.360                               0.360    0.092       restrict ~~         maturity          0.089                               0.089    0.129          Intercepts:         commsize      146062.095                            146062.095    0.228         team.size         9.221                               9.221    0.712         prj.age         162.213                             162.213    7.311         license          13.953                              13.953    2.609         restrict          1.656                               1.656    1.764         maturity          1.758                               1.758    2.397          Variances:         team.size       167.488                             167.488    1.000         prj.age         492.298                             492.298    1.000         license          28.612                              28.612    1.000         restrict          0.882                               0.882    1.000         maturity          0.538                               0.538    1.000         commsize      410630472249.888                            410630472249.888    1.000                     tem.sz prj.ag licens rstrct matrty commsz     42487      1      0      0      0      1      1     14655      1      0      0      0      0      1     13632      1      0      0      0      0      0     6293       1      0      0      0      1      0     4901       0      0      0      0      1      0     4447       1      1      0      0      0      1     3786       0      0      0      0      1      1     1599       1      1      0      0      1      1     899        1      0      1      1      1      1     327        0      1      0      0      1      0     282        1      0      1      1      0      0     251        1      0      1      1      0      1     199        0      0      1      1      0      0     90         0      0      1      1      1      0     84         0      1      0      0      0      0     69         1      0      1      1      1      0     59         1      1      1      1      1      1     53         0      0      0      0      0      1     17         1      1      0      0      0      0     10         1      1      1      1      0      0     3          0      1      1      1      1      0     2          0      1      0      0      1      1     1          1      1      0      0      1      0                    tem.sz prj.ag licens rstrct matrty commsz     team.size 0.847     prj.age   0.061  0.065     license   0.016  0.001  0.019     restrict  0.016  0.001  0.019  0.019     maturity  0.514  0.020  0.011  0.011  0.605     commsize  0.644  0.061  0.012  0.012  0.488  0.682      The following is the R **source code** snippet, used to perform this analysis:               my.model.4 <- '       # means       commsize  ~ 1       team.size ~ 1       prj.age   ~ 1       license   ~ 1       restrict  ~ 1       maturity  ~ 1            # variances       team.size ~~ team.size       prj.age   ~~ prj.age       license   ~~ license       restrict  ~~ restrict       maturity  ~~ maturity            # covariances/correlations       team.size  ~~ prj.age  + license  + restrict + maturity       prj.age    ~~ license  + restrict + maturity       license    ~~ restrict + maturity       restrict   ~~ maturity     '          mcarFIML.fit <- sem(my.model.4, data=flossData, missing="fiml")     # alternatively, clear model's "means" section & use "meanstructure=TRUE"          message("\nModel 4 - Analysis results:\n")     summary(mcarFIML.fit, fit.measures=TRUE, standardize=TRUE)          # Missing data patterns and covariance coverage     inspect(mcarFIML.fit, 'patterns')      inspect(mcarFIML.fit, 'coverage')


title: Evaluating the $H_0$ of a contigency table

text: I have a 2 way contigency table with variables $A=\{a_1,a_2\}$ and $B=\{b_1,b_2\}$. I have the observed cell frequencies $O_1=a_1b_1$, $O_2=a_1b_2$, $O_3=a_2b_1$ and $O_4=a_2b_2$. My null hypothesis is that $A$ and $B$ are independent variables I am interested to evaluate what is the probability of observing each of the cells given that $H_0$ is true, $\text{Pr}(O_1=a_1b_1|H_0)$. How can I calculate this probability? Is there a bayesian approach or a simulation procedure that could help me obtain this probability?


title: Is it correct to apply Friedman+Nemenyi tests to several runs at the same dataset?

text: We have 7 classification models generated from 30 subsets sampled from the same data-set (via cross-validation). We would like to know if there are significant differences among the accuracy of the 7 models. Apparently, in a comment from here, there is no independence assumption for the Friedman test. Is that true?


title: using decompose function for high frequency data

text: I have a table as               Date       Time    Energy     1/1/2008   10:30   0.89     1/1/2008   11:30   0.76      and so on. The data is recorded for every half an hour. I wish to decompose this time series in R.


When the confidence interval of an odds ratio includes one, can the p-value be less than 0.05? For example,               odds ratio = 0.54  confidence interval ( 0.29 - 1.01 ) and      p value is 0.039.      In this situation can we consider the effect to be statistically significant or not? Or is something wrong?


title: How does one build/find a mathematical model for a convex relationship?

text: I think my various questions on this site about spline and polynomial regression boil down to this: I have many datasets to analyze where individuals grow over time and, depending on treatment group, either plateau out or decline. Reporting that the "shapes of the curves are significantly different" won't cut it, because to the casual observer that begs the question "in what way, and by how much", and I think the casual observer is right. What strategy should I follow in modifying that model to accommodate both peaking and plateauing groups? The ones that follow a plateau trajectory can probably be fit by a simple logistic model. For the growth-and-decline groups, is it a matter of fitting a weighted sum of an increasing and a decreasing logistic function? For example, in R or Splus `nls(y~a1*SSlogis(b1,b2,b3,b4)-a2*SSlogis(b5,b6,b7,b8))` (using that just an example, this is not intended as solely an R-specific question). More importantly, I don't want to reinvent the wheel. Curves with maxima and/or minima are such a common problem in biology, I'm sure there must be an established model for this case, I just don't know the correct name to search for. So, I'm hoping someone can point me in the right direction. PS: polynomials and splines seem to offer only a black-box approach; I'm looking for a model parametrized in a way that's interpretable and gives some insight into underlying mechanisms-- "e.g.: this coefficient is positive, therefore this treatment group's curve is shifted to the right of the control group's". That's supposed to be the forte of nonlinear models, but every example I can find is about non-decreasing or non-increasing trends.


title: How to get $P(X < x)$ and $P(X \geq x)$?

text: After doing some reading about the binomial distribution I found this about lower and upper bound probabilities We can get the **lower tail probability** of X as: $$P(X \leq x)=P(X=0)+...+P(X=x)$$   * **Stata:** as `binomial(n,k,p)`   * **R:** `pbinom(q, size, prob, lower.tail = TRUE)` And the **upper tail probability** can be obtained as the complement of lower tail probability to find the area above the cutoff x-value $$P(X>x)= 1-P(X \leq x)$$   * **Stata:** `1- binomial(n, k, p)`   * **R:** `pbinom(q, size, prob, lower.tail = FALSE)` ![enter image description here](http://i.imgur.com/glxH4.png?1) Am I right for the R and Stata commands? And also how can I get these probabilities? Are they complementary, as upper and lower bound probabilities? $$P(X < x)$$ $$P(X \geq x)$$


title: Meaning of return period on extreme events

text: If the distribution of the periods between an extreme event to another is a power law (as for example can be the return period of extreme earthquakes or flooding), the existence of the mean value is subordinated to the value of the exponent of the power law. In the cases where the exponent don't allow the existence of the mean or variance what is the meaning and the utility of the return period?


title: Intuitive explanation of a p value

text: I have a dilemma regarding p -value. The idea is that a p value results from a comparison between a test statistic and control group such that the null hypothesis is assumed and based on this comparative score, p, we can determine the strength of the assumption that the observed results were by chance, or determined by some factor. Here is my question, I was under the impression that the under the t-statistic this statistic is created by a region under the t-distribution. this region if it is large, determines the p value, if it is small it will also determine the p value. So if this is correct, this t distribution seems to imply that you just need the t distribution for a single event, and not two groups, but a single group that you are using the t distribution for. However, this is not correct because for a p value to exist you need two groups. can someone please help expain this ?


title: Interpretation and usage of Kolmogorov-Smirnov Test in Python

text: I have a sample of data, and I want to know weather it is Gaussian-distributed or not. The mean of my data is not zero. To check weather I'm using the K-S test correct, I generated some Gaussian- distributed data and added some bias:               data = stats.norm.rvs(size=10000) + 1      print stats.kstest(data, 'norm')      This gives a p-value of 0.0. If I subtract the bias, I get something like 0.7-0.8, depending on the seed, of course. Does the data need to have $\mu = 0$? If so, does $\sigma^2 = 1$? What if my distribution has a different, unknown $\sigma^2$?


title: Good resource to understand anova and ancova?

text: I am running experiments for a paper and am looking for an interesting book/website to understand properly how anova and ancova work. I have a good math background so I don't necessarily need a vulgarized explanation. I'd also like to know how to determine when to use anova versus ancova.


title: Can I use a likelihood ratio test when the error distributions differ?

text: This question might sound simple, but in fact I am absolutely unsure whether it is allowed to use an LRT for model comparison when the error distributions differ. For example, can I compare a model with a gaussian distribution with a gamma or a log-gamma model?


title: What does it mean to have two intercepts?

text: Suppose I am modeling the sales for multiple products for multiple stores. What is the difference between: \begin{align} \text{Sales} &= \beta_{0}(\text{Intercept × Store}) +\beta_{1}(\text{Intercept × Product}) \\ &\text{and} \\ \text{Sales} &= \beta_{0}(\text{Intercept × Store}\times\text{Product}) \end{align} What is the difference between these two regressions?


title: How does "thin" setting affect the number of samples in WinBUGS?

text: I do not understander how "thin" setting affects the number of samples in WinBUGS. Here is my case: Case 1: In `Model --> Update... -> Update Tool`, I set `updates` 5000, `refresh` 100, `thin` 1, click `update`. In `Inference -> Samples...-> sample` is 4500 in stats Case 2: In `Model --> Update... -> Update Tool`, I set `updates` 5000, `refresh` 100, `thin` 2, click `update`. In `Inference -> Samples...-> sample` is 4750 in stats Case 3: In `Model --> Update... -> Update Tool`, I set `updates` 5000, `refresh` 100, `thin` 5, click `update`. In `Inference -> Samples...-> sample` is 4900 in stats Case 4: In `Model --> Update... -> Update Tool`, I set `updates` 5000, `refresh` 100, `thin` 10, click `update`. In `Inference -> Samples...-> sample` is 4950 in stats Based on WinBUGS User Manual, thin is the samples from every kth iteration will be stored, where k is the value of thin. Setting k > 1 can help to reduce the autocorrelation in the sample, but there is no real advantage in thinning except to reduce storage requirements and the cost of handling the simulations when very long runs are being carried out. My questions are   1. In Case 1, why sample is 4500 not 5000?   2. How does the "thin" setting affect the sample size in Node Statistics? The following is my code.               model {     for(i in 1: N) {     CF01[i] ~ dnorm(0, 20)     CF02[i]  ~ dnorm(0, 1)     h[i] ~ dpois (lambda [i])     log(lambda [i]) <- beta0 + beta1*CF03[i] + beta2*CF02[i] + beta3*CF01[i] + beta4*IND[i]     }     beta0 ~ dnorm(0.0, 1.0E-6)     beta1 ~ dnorm(0.0, 1.0E-6)     beta2 ~ dnorm(0.0, 1.0E-6)     beta3 ~ dnorm(0.0, 1.0E-6)     beta4  <- log(p)     p ~ dunif(lower, upper)     }          INITS     list(beta0 = 0, beta1 = 0, beta2 = 0, beta3 = 0, p = 0.9)          DATA(LIST)     list(N = 15, lower = 0.80, upper = 0.95,          h = c(1,4,1,2,1,2,1,1,1,3,3,0,0,0,NA),          CF03 = c(-1.5,0.1,1.0,0.1,-0.7,0.6,0.2,     0.1,-0.3,1.9,-1.5,0.2,1.0,-0.3,0.8),          CF02 = c(NA,NA,0.3,0.1,-0.9,-0.1,-0.2,-0.7,     -0.9,2.3,1.4,1.2,1.2,-0.7,-1.5),          CF01 = c(NA,NA,NA,-0.1,-0.2,-0.3,-0.2,-0.2,-0.2,     -0.1,0.1,-0.2,-0.2,-0.1,0.1),          IND = c(1,1,0,0,0,0,0,0,0,0,0,0,0,0,0))


title: What is an "If-Then" rule?

text: Can someone give a concise, layman's explanation of an "If-Then" rule (as in rule-based systems). I am finding this term used frequently without anyone really defining it properly.


title: Testing assumptions of multiple regression

text: 1. When doing a multiple regression and testing for homoscedasticity some people look at raw observations and others the residuals. Which is correct?   2. Do you use raw data or residuals to test linearity?   3. Do you test the homoscedasticity for each IV against the DV or do you put all IVs in at the same time and then test for homoscedasticity?   4. When do you test the assumptions before running the analysis, after, or both?   5. What order do you do these things? Do you do any twice?     * test for linearity     * test for normal distribution     * test for equal variances     * run the multiple linear regression


title: How to compare two distributions, where the numbers in each distribution are proportions?

text: Could anyone suggest a statistical technique for comparing two distributions, where the numbers in each distribution are proportions? I have two distributions, each of which reports the proportion of entities in one of 12 classes. Unfortunately, I don't know much about the sample size for one of these distributions, so I'm comparing them on the basis of proportions: I can do so visually, but I'd also like to do so statistically. I'd use the chi-square test for equality of proportions, were the numbers in each distribution counts. In this case, the null hypothesis would be that the distribution of classes is essentially the same. However, I don't think this test is appropriate in this case. Could anyone suggest an alternative? Thanks in advance for any help.


title: Degrees of freedom for Gaussian Process

text: I am reading this paper on Generalised Wishart Process (GWP). It is about modelling covariance matrix of D - dimensional gaussian processes (GP) as GWP. I fail to understand interpretation of "degrees of freedom" $\nu$ for a GP. There are quiet a few mention of degree of freedom in the paper. In section 4, they define GWP as sum of outer product of $\nu$ D-dimensional gaussian processes. Based on this definition $\nu$ seemed to be no. of observations for the gaussian processes. But then in the same section, they also mention about real-valued degrees of freedom, which rules out previous explanation that $\nu$ = no. of observations/data points. Moving on, Section 5 Bayesian Inference says, degrees of freedom control how concentrated the prior is around the expected value $\Sigma(t)$ (This generic statement does make sense). Also, in section 5.1, they sample from posterior distribution over the gaussian processes with different degrees of freedom (which I fail to understand again). In section 5, they set $\nu$ to its minimum value D + 1 (so that the wishart distribution does not become degenerate). I'd appreciate if you can provide me with answers or any other resources which can help me understand this.


title: Predicting from the posterior of a random coefficients model

text: I've got a model that looks like this: $$Y_{ig} = \left(\beta_{everyone} + \beta_g\right)X_{ig} + Z_{ig}'\gamma + \epsilon_{ig}$$ in R, this is               library(lme4)     example = lmer(Y~ X + Z + (0 + X|g))      I want a conditional estimated distribution of $Y|G=g,X=x,Z=z$ For lm, gam, etc, this is just `predict(example)` as the mean of a normal distribution with the same variance as the residuals of the model divided by n-k. Is there anything different in a mixed model/random coefficients case? I don't think MCMC is appropriate, as I don't want to estimate the probability of being ins a certain group $g$ -- I want a distribution for a certain $g$. If indeed the same intuition is applicable, how would I code this? I'm relatively new to mixed effects models, so any advice is appreciated.


title: Forecasting of highly correlated time series

text: In time series forecasting using various models like AR,MA,ARMA, etc, we usually focus on the modeling of the data in the change of time. But when we have 2 time series that **Pearson correlation coefficient** show they are highly correlated, is it possible to model their dependency and forecast values of one from the other? For example, when one serie have a linear relationship with the other, it seems possible. But is there a **general method** for this kind of dependency analysis?


title: Confusion with Vowpal Wabbit's multiple-pass behavior when performing ridge-regression

text: I have encountered many peculiarities/misunderstandings of Vowpal Wabbit when trying to do online multiple-pass learning. Specifically, I need to solve a Ridge Linear regression problem, with `N=4e6` points and a total of around `K=2.38e5` features. Each point has sparsely populated features, usually between 10 and 100. Since `N` is too large to fit in my PC's limited memory, I decided to use Vowpal Wabbit's out-of-core SGD. Features are in text representation:               vw -d train.dat -c -f train.model --passes 10 --loss_function squared         --l2 0.00001 -l 0.05          Num weight bits = 18     learning rate = 0.05     initial_t = 0     power_t = 0.5     decay_learning_rate = 1     using l2 regularization = 1e-05     using cache_file = vw/Train_0.cache     ignoring text input in favor of cache input     num sources = 1     average    since         example     example  current  current  current     loss       last          counter      weight    label  predict features     0.028051   0.028051            3         3.0   0.0010   0.0070       24     0.089131   0.150212            6         6.0   0.0036   0.0256       20     0.131556   0.182465           11        11.0   0.0076   0.0297       47     0.102226   0.072896           22        22.0   0.0052   0.0368       16     0.095826   0.089426           44        44.0   0.0539   0.1463       64     0.094730   0.093608           87        87.0   0.0058   0.0605        8     0.084809   0.074887          174       174.0   0.0177   0.0914       13     0.073454   0.062099          348       348.0   0.2019   0.1518       18     0.065183   0.056912          696       696.0   0.0044   0.1866       20     0.061144   0.057106         1392      1392.0   0.8041   0.2107       12     0.060112   0.059079         2784      2784.0   0.0454   0.1215       15     0.054327   0.048543         5568      5568.0   0.0875   0.0016       58     0.052169   0.050009        11135     11135.0   0.0194   0.0785       41     0.048767   0.045365        22269     22269.0   0.6126   0.1835       42     0.046333   0.043900        44537     44537.0   0.0149   0.0000       55     0.045165   0.043997        89073     89073.0   0.4736   0.1643       15     0.043814   0.042463       178146    178146.0   0.0561   0.0189       51     0.042698   0.041581       356291    356291.0   0.0161   0.0412       54     0.042067   0.041436       712582    712582.0   0.0012   0.1215       25     0.041751   0.041435      1425163   1425163.0   0.0540   0.0677       33     0.042095   0.042439      2850326   2850326.0   0.1088   0.1141       41     0.000000   0.000000      5700651   5700651.0   0.0028  -0.0000       43 h     0.000000   0.000000     11401301  11401301.0   0.0006   0.1667       32 h     0.000000   0.000000     22802601  22802601.0   0.0534   0.1617       12 h          finished run     number of examples per pass = 3702833     passes used = 10     weighted example sum = 3.70283e+07     weighted label sum = 5.96004e+06     average loss = 0 h     best constant = 0.160959     total feature number = 1244815410      The `h` on the right end of the last three rows indicates that the loss is calculated on the holdout's validation-subset of data. Why am I getting an average loss of 0, when presumably the holdout predictions are incorrect, as indicated, again, by the difference between `current label` and `current predict` of the last 3 rows? This oddity is obscuring the fact of whether or not my algorithm has converged... I can bypass this obscurity with a work-around, by ignoring `vw`'s `\--passes` option and manually performing "passes" along with IO-intensive shuffling of data after each "pass":               for i in {0..10}     do         # Custom script that shuffles train.dat -> train.dat         ./shuffledata train.dat train.dat         vw -d train.dat -c -k --l2 0.00001 -l 0.05 -i train.model -f train.model     done      However, the above two codes (with the `\--passes` option and without the `\--passes` option + shuffling data script) performs almost exactly the same as the following snippet:               vw -d train.dat -c --l2 0.00001 -l 0.05 -f train.model      In other words, **it's as if passes are not needed at all for "convergence."** These experiences lead me to my question(s).   1. Why does multiple-pass (shuffled & non-shuffled) give the same result as single-pass in an online SGD setting?   2. Why is my holdout error 0?


title: How can the Kolmogorov-Smirnov test be used/interpreted?

text: The test of Kolmogorov-Smirnov (K-S) is a traditional test of normality, although Shapiro-Wilk test (S-W) is applied more frequently than K-S (Arango, 2012). I am not an expert in statistics, so my question concerns about the use of K-S test.   1. Is it possible to use K-S test to another purposes besides normality test?    2. Why some times we should use K-S test instead of S-W test? Is it related to sample size?


title: Macro and Micro Inter-Coder/Rater Agreement

text: I have a labeling task with 4 labels for which I'm using 2 coders/annotators/raters. This labeling effort leads to a confusion matrix like the one in https://docs.google.com/drawings/d/1T7pkLLkE7qvKanxqvZO_uQ_q3hz6S3fcTbOnWophx38/edit?usp=sharing I want to measure the agreement between the 2 annotators. How so? Kappa Coefficient works for binary labeling/classification but not for this case. I thought of having like a micro kappa for each 2 labels and then an averaged macro Kappa across all combinations. How dummy does that sound? Thanks


title: Pandas Statsmodels Time series seasonal forecasting

text: Using Stats models and Pandas (and requests for the data) I'm working on a forecast model.. my 1st step is just getting the Arma function working and understood. My data is available publically and is highly seasonal residential real estate unit sales data, I'm planning to see how a quarterly survey that we do can help with the forecast as a later step. Hence why I am changing the frequency to quarterly with dates that match the dates I have on the quarterly survey. So the code looks like:               #get the statewide actual data from Sheet 1 parse dates and select just unit sales     act = requests.get('https://docs.google.com/spreadsheet/ccc?key=0Ak_wF7ZGeMmHdFZtQjI1a1hhUWR2UExCa2E4MFhiWWc&output=csv&gid=1')     dataact = act.content     actdf = pd.read_csv(StringIO(dataact),index_col=0,parse_dates=['date'], thousands=',') #converts to numbers     actdf.rename(columns={'Unit Sales': 'Units'}, inplace=True)     actdf=actdf[['Units']]     actdfq=actdf.resample('Q',sum)     actdfq.index = actdfq.index + pd.DateOffset(days=15) #align the actual data dates to the   survey dates Eg the 15th of the quarter     actdfq=actdfq['2009':] # selcts the time periods for which we have surveys (actual results here) The survey would be shifted back by one     actdfqchg=actdfq['Units']     fig = plt.figure(figsize=(12,8))     ax1 = fig.add_subplot(211)     fig = sm.graphics.tsa.plot_acf(actdfqchg.values.squeeze(), lags=4, ax=ax1)     ax2 = fig.add_subplot(212)     fig = sm.graphics.tsa.plot_pacf(actdfqchg, lags=4, ax=ax2)      the data looks like:               2009-01-15     7867     2009-04-15     7483     2009-07-15    10109     2009-10-15    10648     2010-01-15     9678      The acf graphs look like: ![enter image description here](http://i.stack.imgur.com/PKuA1.png) So I don't really know what the ACF is telling me..The graph of auto correlation would tend to tell me that the 4q correlation is the strongest but still only .4? (correct?) and the 2q score of -.4 would indicate that the summer to winter correlation would be the weakest (which makes sense) and how to proceed with 1) a projection based on this actual data just using straight Arma Stats models capabilities. and 2) how to best incorporate survey data that attempt to predict the following quarter.. the survey is taken asking for predictions for the following quarter on a 5 point scale with a400+ participants, the straight correlation is not super strong but I think somehow I should be able to find the quarterly correlation to help inform the projection for the subsequent quarter..??


title: What is the word to mean "the actual data points", not the trends that it creates?

text: I want to talk about some data that is only relevant when looked at as a whole, to see the trends, but each individual piece of data is suspect, similar to web traffic data. What word should I use in this blank space? "This data is not supposed to be analyzed as **___** , it is only good for showing trends and getting a total statistical picture of the whole data set."


title: Forecasting a time series with weights

text: I'd like to forecast (or predict) a time series with weights. The following works using the regular `l`inear `m`odelling techniques of `lm` by applying a (sigmoidal) weight distribution to the input data, essentially weighing the latter data points more heavily than the former:               library("stats")     lm.weight.function <- function(x) {10 / (1 + exp(-x))} # Sigmoidal     lm.weights <- lapply(seq(-13, 14, length.out = 27), lm.weight.function)     lm.input <- as.data.frame(c(23957, 46771, 60767, 73284, 60296, 73122, 78304, 87154, 80459, 76885, 56479, 18809, 13453, 13951, 25140, 12035, 11920, 20683, 30357, 35019, 37732, 46150, 47856, 41931, 20985, 32526, 27283))     lm.input <- cbind(1:27, lm.input)     colnames(lm.input) <- c('x', 'y')     lm.model <- lm(formula = y ~ log(x), data = lm.input, weights = unlist(lm.weights))     predict.input <- as.data.frame(28:55)     colnames(predict.input) <- 'x'     predict.model <- predict(lm.model, predict.input)     plot(1:(27+28), c(lm.input$y, predict.model), type = 'l', xlab = 'x', ylab = 'y')      ![enter image description here](http://i.stack.imgur.com/5M4FP.png) Now I wish to do the same using the `forecast` package. However, I'm having difficulty specifying the `weights`:               library("forecast")     ts.weight.function <- function(x) {10 / (1 + exp(-x))} # Sigmoidal     ts.weights <- as.data.frame(lapply(seq(-13, 14, length.out = 27), ts.weight.function))     colnames(ts.weights) <- 'trend'     ts.input <- ts(c(23957, 46771, 60767, 73284, 60296, 73122, 78304, 87154, 80459, 76885, 56479, 18809, 13453, 13951, 25140, 12035, 11920, 20683, 30357, 35019, 37732, 46150, 47856, 41931, 20985, 32526, 27283), frequency = 1)     ts.model <- tslm(formula = ts.input ~ log(trend), weights = unlist(ts.weights))      The above prints an error:               Error in eval(expr, envir, enclos) :        ..1 used in an incorrect context, no ... to look in      How can I use `tslm` to forecast a time series with weights?


title: Showing $ λ_V(x) \leq min\{λ_1(x),\cdots,λ_n(x)\}$ Hazard function

text: > Suppose $X_1, \cdots, X_n$ are independent, nonnegative continuous > functions, each $X_i$ has hazard function $\lambda_i(x)$. > > If $V=\max\{X_1, \cdots, X_n\}$, I need to show that $\lambda_V(x)\leq > \min\{\lambda_1(x),\cdots, \lambda_n(x)\}$. I know that $\displaystyle\lambda_V(x)=\frac{f_V(x)}{1-F_V(x)}=\frac{\sum_j f_j(x)\prod_{i\neq j}F_i(x)}{1-\prod_{i=1}^n F_i(x)}$. However, I'm stuck here. Any help would be appreciated.


title: History: the role of statistics in astronomy

text: I recently boldly claimed in front of a group of fairly smart eighth grade students that astronomy contributed greatly to the foundations of statistics and many statistical concepts were invented for use in astronomy. However, looking to back that up, I was fairly disappointed. Errors, the mean and the median deviation from the mean may have been first observed in astronomy. However, even the concept of error propagation might stem more from classical mechanics than astronomy. Beyond these concepts, I was unable to find much more. Feigelson writes (http://arxiv.org/pdf/astro-ph/0401404.pdf): > Ptolemy estimated parameters of a non-linear cosmological model using a > minimax goodness-of-fit method. Al-Biruni discussed the dangers of > propagating errors from inaccurate instruments and inattentive observers. > While some Medieval scholars advised against the acquisition of repeated > measurements, fearing that errors would compound rather than compensate for > each other, the usefulnes of the mean to increase precision was demonstrated > with great success by Tycho Brahe. Can you suggest good references that have some more details on the historical links between astronomy and statistics? Thank you for the excellent answers!


title: Example of how the log-sum-exp trick works in Naive Bayes

text: I have read about the log-sum-exp trick in many places (e.g. here, and here) but have never seen an example of how it is applied specifically to the Naive Bayes classifier (e.g. with discrete features and two classes) How exactly would one avoid the problem of numerical underflow using this trick?


title: Why will a random forest not outperform a regression tree?

text: I have a training dataset with a binary response variable, 6 independent variables, and 21,000 observations. I've fit both an ordinary regression tree and a random forest (mtry = 2, ntree = 2000) and there is almost no difference between the two when each model is validated, using RMSE and predicted to actual ratio as goodness of fit metrics. Is this to be expected with a small number of independent variables, or am I not using the right metrics to measure goodness of fit?


title: Information theoretic central limit theorem

text: The simplest form of the information theoretic CLT is the following: Let $X_1, X_2,\dots$ be iid with mean $0$ and variance $1$. Let $f_n$ be the density of the normalized sum $\frac{\sum_{i=1}^n X_i}{\sqrt{n}}$ and $\phi$ be the standard Gaussian density. Then the information theoretic CLT states that, if $D(f_n\|\phi)=\int f_n \log(f_n/\phi) dx$ is finite for some $n$, then $D(f_n\|\phi)\to 0$ as $n\to \infty$. Certainly this convergence, in a sense, is "stronger" than the well establised convergences in the literature, convergence in distribution and convergence in $L_1$-metric, thanks to Pinsker's inequality $\left(\int |f_n-\phi|\right)^2\le 2\cdot \int f_n \log(f_n/\phi)$. That is, convergence in KL-divergence implies convergence in distribution and convergence in $L_1$ distance. I would like to know two things.   1. What is so great about the result $D(f_n\|\phi)\to 0$?   2. Is it just because of the reason stated in the third paragraph we say convergence in KL-divergence ( _i.e._ , $D(f_n\|\phi)\to 0$) is stronger? Thanks. NB: I asked this question sometime ago in math.stackexchange where I didn't get any answer.


title: Significance of $R^2$

text: I have set up a linear regression array in Excel and now want to test the significance of my $R^2$ value at a certain level of significance. I've only been able to find tables that give the critical $R^2$ value, but I want to test it at 99.73% level of confidence and none of them contain that specific value. Is there a way to calculate the critical $R^2$ value by hand?


title: which type of statistical analysis for survey of experts

text: I have a set of 33 questionnaires, each questionnaire has 4 responses, the purpose of analysis is to evlaute the existing flood management system in country, trend in response from 33 experts, and to assess the potential of the system to integrate strategic environmental assessment in the system at planning level, can any body suggest me which type of statistical analysis can help me in better interpretation of the data


title: Analyzing multiple choice survey questions

text: Consider the following survey question: * * * Q: Choose one or more of the following 5 items: A B C D E * * *   * How can one test which items are more frequently chosen for a sample with 100 individuals? Is it advisable to fit a distribution to the data?


title: What is this "maximum correlation coefficient"?

text: A typical image processing statistic is the use of Haralick texture features, which are 14. I am wondering about the 14th of these features: Given an adjacency map $P$ (which we can simply view an the empirical distribution of two integers $i,j < 256$), it is defined as: _the square root of the second eigenvalue of $Q$_ , where $Q$ is: $Q_{ij} = \sum_k \frac{ P(i,k) P(j,k)}{ [\sum_x P(x,i)] [\sum_y P(k, y)] }$ Even after much googling, I could not find any references for this statistic. What are its properties? What is it representing? (The value $P(i,j)$ above is the normalised number of times that a pixel of value $i$ is found next to a pixel of value $j$).


title: Accuracy changes with permutation of input data in Libsvm with precomputed kernel?

text: I'm doing quite simple SVM classification at the moment. I use a precomputed kernel in LibSVM with RBF and DTW. When I compute the similarity (kernel-) matrix, everything seems to work very fine ... until I permute my data, before I compute the kernel matrix. An SVM is of course invariant to permutations of input-data. In the below Matlab-code, the line marked with '<- !!!!!!!!!!' decides about the classification accuracy (not permuted: 100% -- permuted: 0% to 100%, dependant on the seed of rng). But why does permuting the file-string-array (named fileList) make any difference? What am I doing wrong? Have I misunderstood the concept of 'permutation invariance' or is it a problem with my Matlab-code? My csv-files are formatted as: LABEL, val1, val2, ..., valN and all the csv- files are stored in the folder dirName. So the string array contains the entries '10_0.csv 10_1.csv .... 11_7.csv, 11_8.csv' (not permuted) or some other order when permuted. I also tried to permute the vector of sample serial numbers, too, but that makes no difference.               function [SimilarityMatrixTrain, SimilarityMatrixTest, trainLabels, testLabels, PermSimilarityMatrixTrain, PermSimilarityMatrixTest, permTrainLabels, permTestLabels] = computeDistanceMatrix(dirName, verificationClass, trainFrac)     fileList = getAllFiles(dirName);     fileList = fileList(1:36);     trainLabels = [];     testLabels = [];     trainFiles = {};     testFiles = {};     permTrainLabels = [];     permTestLabels = [];     permTrainFiles = {};     permTestFiles = {};          n = 0;     sigma = 0.01;          trainFiles = fileList(1:2:end);     testFiles = fileList(2:2:end);          rng(3);     permTrain = randperm(length(trainFiles))     %rng(3); <- !!!!!!!!!!!     permTest = randperm(length(testFiles));          permTrainFiles = trainFiles(permTrain)     permTestFiles = testFiles(permTest);          noTrain = size(trainFiles);     noTest = size(testFiles);          SimilarityMatrixTrain = eye(noTrain);     PermSimilarityMatrixTrain = (noTrain);     SimilarityMatrixTest = eye(noTest);     PermSimilarityMatrixTest = eye(noTest);          % UNPERM     %Train     for i = 1 : noTrain         x = csvread(trainFiles{i});            label = x(1);         trainLabels = [trainLabels, label];         for j = 1 : noTrain             y = csvread(trainFiles{j});                         dtwDistance = dtwWrapper(x(2:end), y(2:end));             rbfValue = exp((dtwDistance.^2)./(-2*sigma));             SimilarityMatrixTrain(i, j) = rbfValue;             n=n+1         end     end          SimilarityMatrixTrain = [(1:size(SimilarityMatrixTrain, 1))', SimilarityMatrixTrain];          %Test     for i = 1 : noTest         x = csvread(testFiles{i});         label = x(1);         testLabels = [testLabels, label];         for j = 1 : noTest             y = csvread(testFiles{j});                         dtwDistance = dtwWrapper(x(2:end), y(2:end));             rbfValue = exp((dtwDistance.^2)./(-2*sigma));             SimilarityMatrixTest(i, j) = rbfValue;             n=n+1         end     end          SimilarityMatrixTest = [(1:size(SimilarityMatrixTest, 1))', SimilarityMatrixTest];          % PERM     %Train     for i = 1 : noTrain         x = csvread(permTrainFiles{i});                 label = x(1);         permTrainLabels = [permTrainLabels, label];         for j = 1 : noTrain             y = csvread(permTrainFiles{j});                         dtwDistance = dtwWrapper(x(2:end), y(2:end));             rbfValue = exp((dtwDistance.^2)./(-2*sigma));             PermSimilarityMatrixTrain(i, j) = rbfValue;             n=n+1         end     end          PermSimilarityMatrixTrain = [(1:size(PermSimilarityMatrixTrain, 1))', PermSimilarityMatrixTrain];          %Test     for i = 1 : noTest         x = csvread(permTestFiles{i});         label = x(1);         permTestLabels = [permTestLabels, label];         for j = 1 : noTest             y = csvread(permTestFiles{j});                         dtwDistance = dtwWrapper(x(2:end), y(2:end));             rbfValue = exp((dtwDistance.^2)./(-2*sigma));             PermSimilarityMatrixTest(i, j) = rbfValue;             n=n+1         end     end          PermSimilarityMatrixTest = [(1:size(PermSimilarityMatrixTest, 1))', PermSimilarityMatrixTest];          mdlU = svmtrain(trainLabels', SimilarityMatrixTrain, '-t 4 -c 0.5');     mdlP = svmtrain(permTrainLabels', PermSimilarityMatrixTrain, '-t 4 -c 0.5');          [pclassU, xU, yU] = svmpredict(testLabels', SimilarityMatrixTest, mdlU);     [pclassP, xP, yP] = svmpredict(permTestLabels', PermSimilarityMatrixTest, mdlP);          xU         xP          end      I'd be very thankful for any answer! Regards Benjamin


title: Determine whether a n-th finite moment of X exists

text: I have a question which asks: > Determine those values of the positive integer n for which a finite nth > moment of X about zero exists. How should I approach this question? Does it depend on the numbers of variables in X? I think that a first moment exists if the mean exists, and the second moment exists if the variance can be determined. Is that general concept correct?


title: Generating over-dispersed counts data with serial correlation

text: Could anyone provide some suggestions on how to generate over-dispersed counts data with serial correlations? I am using R software to conduct a simulation study. Any references on this subject will be much appreciated. Thanks for your help.


title: Cross-sectional panel data and a fixed effects model

text: Assume a dataset that has a time series of a dependent variable for a number of stocks, and then a time series of an independent variable over time. Would applying a fixed effects model on this mean that all unobserved characteristics of the stocks would be included in the intercepts, and all changes should therefore be due to changes in the independent variable(s)? With the assumption that the unobserved characteristics stay stable over time. What if an unobserved thing affects all the stocks and is thus not a characteristic of a specific stock? The model would look like this for example:               dependent(i,t) ~ a(i) + b1*independent(i,t)      where _i_ is the specific stock. Someone presented me this about fixed effects models and I'm trying to understand whether this really is true. For example, if I was to examine whether $X$ affects e.g. the liquidity of the market in general, would a fixed effects model rule out anything else than omitted things about the specific stocks? What if there is $Y$ that also affects all the stocks, I believe this would still be considered omitted variable bias? The story behind here is that this someone argues that all changes can be said to be result of $X$ thanks to using a fixed effects model. But if there is also $Y$ that has a similar effect on all stocks as $X$, would omitting $Y$ here overestimate the effect of $X$? Hopefully someone understands my questions, it is basically about the general characteristics of a fixed effects model, and what exactly can be omitted and what not.