title: How do you use draws from a Markov chain (Monte Carlo) for (bayesian) inference?

text: Say $X_0, X_1, ...$ is a Markov chain with the state space is $S=\{0,1,2,3,...,N-1\}$ and the transition probability $P(X_{n+1}=a+1 \mod N | X_n=a)$ = $P(X_{n+1}=a \mod N | X_n=a)$ = $P(X_{n+1}=a-1 \mod N | X_n=a)$ = $1/3$. (Ie, the markov chain is a random walk on $S$ with 1/3 probability of staying put, 1/3 probability moving plus or minus 1, with wrap around at the ends.) The limiting distribution is $\pi=(1/N,1/N,...,1/N)$, the uniform dist on $S$. As I understand it, the point of MCMC is to create markov chain such that the stationary distribution is the distribution you want to sample from. So the above chain could be used to sample from the uniform dist on $S$ (silly example, but bare with me). From my understanding, you starting with $x_0$ in $S$, 'burn-in' for some number of iterations, and subsequent iterations are draws from $\pi$. Several questions: First, how do you determine the length of the burn-in period? Even for the trivial example above, the time to the stationary distribution seems highly- dependent on how "mobile" the chain is (for example if $P(X_{n+1}=a \mod N | X_n=a)$ = .9, then it would take longer to spread out the distribution). After burn-in, what role does the dependency between draws play? The literature leads me to believe that I can use draws after the burn-in period as draws from $\pi$, but clearly the individual draws are still dependent on each other according to the transition rule... if you use too few draws with 'stay-still' probability of .9 as above, that's not uniform at all. This seems to mean that any calculation I do with draws (like approximating expectation) must be done of sets of draws greater than the burn-in period, but by how much?


title: ANOVA or any other tests for differences in R

text: I have a timecourse dataset consisting of 3 groups of patients as follows: Group1: control Group2: Medication A to reduce pain after a medical procedure Group3: Medication B to reduce pain after a medical procedure I will like to find out if there are differences between the groups. I wonder if it's correct to carry out the following steps in R:   1. Check normality of data   2. ANOVA   3. pairwise.t.test I'm not really sure which ANOVA approach to use if data normally distributed/not normally distributed and how to find the differences if any. I will be happy for some help. Thanks


title: Non parametric anova or two non parametric ttest?

text: I have a statistical question about a project I recently got. We performed a flow cytometry assay on 3 different groups: controls (n=5), patients before treatment (n=5) and paired patients after treatment (n=5). My plan would be to compare the median expression of the same cell population in controls and patients before treatment using an unpaired non parametric t-test and to compare the median expression of the same cell population in patients before and after treatment using a paired non parametric t-test. Is this correct or should I use ANOVA? In case I should use ANOVA how do I handle the fact that some comparisons are paired and some are not? Any help would be greatly appreciated


title: How to understand if the points are constant around a line?

text: I would like to understand how can I check if some points around a line have constant distance between the point and the line. So image a normal chart(cartesian plane) with a line that is formed with points where its coordinates are: x > 0 and y > 0 so always numbers bigger then zero. Then, when I have this line i need to check if the distances of a list of points(x,y) are constant. With "constant" I mean that the distance between the points and the line is similar. Example: **LINE:** (x - y coordinates)               1 - 1     2 - 2     3 - 3     4 - 4     5 - 5     6 - 6      **POINTS:** (x - y coordinates)               1 - 1.1     2 - 1.9     3 - 2.9     4 - 4.2     5 - 5.1     6 - 5.8      How can I check if the difference between the point on the line and the point (of the list of points) are constant? Thanks


title: Understanding similarity sensitive hashing algorithm in AdaBoost

text: I'm a CS major and don't quite understand the mathematics behind a optimization problem coming from a machine learning algorithm. The algorithm is in Section 5 of the paper http://visl.technion.ac.il/bron/publications/BroBroOvsGuiTOG10.pdf **Problem Statement** given $P$ pairs of examples $(f_p, {f_p}')$ labeled by {1,-1}, where $f_p$ is a v-dimensional feature vector and +1 (reps., -1) indicates similar (reps., dissimilar) pairs, the goal is to find the $s$X$v$ matrix A and $v$X1 vector b such that $d_{A,b}$ reflects the desired similarity of the training examples. The distance $d_{A,b}$ is defined as: $d_{A,b}(x, x')=d_H(\mbox{sign}(Af+b), \mbox{sign}(Af'+b))$, where $d_H(y,y')=\frac{s}{2}-\frac{1}{2}\sum_{i=1}^s\mbox{sign}(y_i{y_i}')$ is the Hamming metric in the s-dimensional Hamming space of binary sequences of length s. Ideally, we would like to achieve $d_{A,b}(f,f')\leq d_0$ for similar pairs (P) and $d_{A,b} > d_0$ for dissimilar ones (N), where $d_0$ is some threshold. However, in practice ,there always exists false positives and false negatives. Thus, optimal A, b should miminize: $min\frac{1}{P}\sum_{(f,f')\in P}e^{\mbox{sign}(d_{A,b}(f,f')-d_0)}+\frac{1}{N}\sum_{(f,f')\in N}e^{\mbox{sign}(d_0-d_{A,b}(f,f')}$ **Algorithm** The learning of the optimal parameters A and b is posed as a boosted binary classification problem, where $\mbox{sign}(Af+b)$ acts as a strong binary classifier and each dimension of the linear projection $\mbox{sign}(A_ix+b_i)$ is a weak classifier. This way, the AdaBoost algorithm (http://www.cs.rochester.edu/users/faculty/stefanko/Teaching/09CS446/Boosting2.pdf) can be used to progressively construct A and b. Intuitively, this algorithm increases the weights of incorrectly classified examples so that the learner is forced to focus on the hard examples in the training set. Input: P pairs of examples $(f_p, {f_p}')$ labeled by $s_p$   1. Initialize ${w_p}^i=1/P$   2. For i = 1, ..., d do Set the ith row of A and b by solving the optimization problem $(A_i, b_i)=\min \sum_{p=1}^P {w_p}^is_p(2-\mbox{sign}(A_if_p+b_i))(2-\mbox{sign}(A_i{f_p}'+b_i))$ Update weights ${w_p}^{i+1}={w_p}^ie^{-s_p\mbox{sign}(A_if_p+b_i)\mbox{sign}(A_i{f_p}'+b_i)}$ **My Question** The optimization in the above algorithm is difficult, so the author tried to solve a simpler problem by setting $A_i=\arg \max \frac{{A_i}^TC_NA_i}{{A_i}^TC_PA_i}$ where $C_P$ and $C_N$ are covariance matrices of the positive and negative example pairs The author says that $A_i$ maximizing the above formula is the largest generalized eigenvector of ${C_P}^{-1/2}{C_N}^{1/2}$. Since this do not coincide exactly with the original optimization problem, the author select a subspace spanned by the largest 10 eigenvectors, out which direction and as well as the threshold parameter b minimizing the exponential loss are selected. I don't understand how the selection is done. Can anyone explain to me? Thanks a lot.


title: Sample size for a repeated measures design with missing values

text: I am planning a two arm study (a treatment group and a control). My endpoint is a continuous variable. Every patient will give between 1 to 3 samples, from different locations in the body. Most of the patients will give either 1 or 2 samples, only a minority will give 3. I know how to calculate a sample size for repeated measures, however, in this case not all patients will have the same number of measures, so the parameter in my formula, the number of measures, is not constant. I need a new method. Does anyone know how to approach this problem? The correlation by the way, between samples within a patient, is quite high, around 0.8.


title: What does empty output in Factor Analysis means in R?

text: factanal(x = ~V1 + V2 + V3 + V4 + V5, factors = 2, data = airfoil_self_noise)          Uniquenesses:        V1    V2    V3    V4    V5      0.898 0.005 0.005 0.995 0.394           Loadings:        Factor1 Factor2     V1 -0.319              V2  0.877  -0.474      V3          0.997      V4                     V5  0.754  -0.194       Here, I have no idea what empty space means in Factor1-V3, Factor1-V4, Factor2-V1, and Factor2-V4. Weren't they supposed to be filled with some numbers, like even with 0?


title: The disadvantage of using F-score in feature selection

text: F-score can be used to measure the discrimination of two sets of real-numbers and can be used for feature selection. However, I once read that > A disadvantage of F-score is that it does not reveal mutual information > among features. How to understand this statement, or why F-score has this kind of disadvantage.


title: Test a proportion with limited information

text: I was reading this article, and wanted to check whether 203/2000 was statistically different from 1/12. Is this a comparison between two means? How can I do this without knowing the standard deviations?


title: Endogeneity without correlation?

text: Is it, in general, possible that two variables are endogenous even if they are not correlated with each other at all?


title: Plotting more than three variables in matlab

text: I am here to seek opinion on how should i represent my data that i have collected. I am to create a presentation focusing on environment. I was told that a simple bar chart and a line graph is a bad visualization I have picked the following data set as a start   * Carbon dioxide emission   * forest area   * energy consumption   * population each of the data set has the following variable ( 55 countries and data ranging from 1990-2010) **Question** - Should i just focus on a single year to show the relation between the large population thus the need for higher consumption thus the reduce in forest area resulting in carbon dioxide emission ? Or should i include everything from 1990 - 2010 ? Also how should i visualize them ? what sort of plot of 3d graph should i be using . Any suggestion is welcome ! thanks in advance I did a couple of samples _This bar chart show the energy consumption of 55 from 1990-2010 with the color bar representing the years_ ![This bar chart show the energy consumption of 55 from 1990-2010 with the color bar representing the years](http://i.stack.imgur.com/PJE0m.jpg) _This graph is the displaying the same data as the first bar chart, just that this is in graph format_ ![This graph is the displaying the same data as the first bar chart, just that this is in graph format](http://i.stack.imgur.com/I3Ptz.jpg) _This is a snip shot of my data source_ ![This is a snip shot of my data source](http://i.stack.imgur.com/LS250.jpg)


title: Which repeated measure statistical test must I use in this case?

text: I have a set of general conditions A, B, C, D, E and F to solve a problem. Each of these conditions lead me to different number of experiments according to a parameter that i change in each general condition:   * condition A has 6 experiments   * condition B has 24 experiments   * condition C has 28 experiments   * condition D has 5 experiments   * condition E has 25 experiments   * condition F has 8 experiments Each experiment returns a 10-dimensional vector of accuracy used in a cross- validation machine learning approach. Now I want to compare the conditions in pairs: I want to know if its better to use condition A than B, condition C than D and Condition E than F. Which statistical test can I use in this case? My suspicion is a repeated- measure based statistical one.


title: Examples of hidden Markov models problems?

text: I read quite a bit of hidden Markov models and was able to code a pretty basic version of it myself. But there are two main ways I seem to learn. One is to read and implement it into code (which is done) and the second is to understand how it applies under different situations (so I can better understand how it relates to problems I might be working on). All the examples I have done so far have involved either some kind of DNA prediction or coin tossing. **I'm wondering if there are any resources to get other Markov problems** (language doesn't matter but hopefully with the answers as well so I can know if I'm right or wrong)?


title: Assumptions of Mann-Whitney U-test

text: I have 2 samples, which I got by conducting a 5-point Likert scale survey. I have 9 results for group A and 17 results for group B for 12 variables. I would like to test now if the groups differ significantly in one or more of these variables and chose to employ the t-test. Assumptions for the t-test are 1. Interval scaled 2. Independent samples 3. Equal variances of smaples 4. Normally distributed. My data is independent, and I have tested variance homogeneity with Levene’s test – nothing was significant – 0-hypothesis accepted – equal variances. I used the Shapiro-Wilk-test for the normality distribution and like expected when working with Likert the variables were – with a few exceptions – not normally distributed. I conducted a t-test nevertheless and came to, in my opinion, reasonable results. However, since there is no consensus whether to treat Likert as ordinal or metric and since the data is not normally distributed, I also conducted a Mann-Whitney-U-Test (MWU-Test). This one gave me almost the same results as the t-test (for the t-test 3 variables were significant and for MWU-test 2 of these 3 variables were significant). Now I read under assumptions for MWU-test that the population (the sample…) should have the same shape. I don’t understand that assumption. I have conducted boxplots (working with SPSS) but many of the plots do not have the same form. My questions:   1. How can I test the “shape” assumption and what are the implications if my samples do not fulfill this assumption? Is there a high risk of false results? I’m not sure since the results of the t-test and the MWU-test are so similar.   2. Which influence has the size of my samples on the tests (variances, t-test, MWU-test), because n<10 and n<20?   3. Would you go on working with the t-test results (3 significant variables) or the MWU-test (2 significant variables)?   4. After I have concluded that 2 or 3 of the variables are significant, how can I test which group is the “better” one, e.g. crossed mostly “increased” instead of “decreased”? Thank you so much already, I just started dealing with stats and SPSS, so every help is very welcome.


title: Can test a hypothesis with only nominal variables?

text: I have a hypothesis such as "People consult internet before making a purchase decision" for a thesis research. In order to test it, I created a game where people had the option to consult "internet option" or not. I have a similar questions, but I have now realised that all my variables for this particular point are nominal; (the answers look 0: did not consult internet, 1: consulted internet) However, I can't calculate a mean or do any of the tests like $t$-test, anova since it doesn't make sense. I was wondering academically, my results would be enough to confirm (not cannot reject) that people consult internet if I have 80% of the same that consulted internet? Several different questions show similar results.


title: Greater than 1 Naive Bayes Probabilities?

text: I am trying to train a Naive Bayes classifier. In addition to getting the most likely class as an output from the Naive Bayes classifier, I would also like to compute the probabilities associated with labels. I am making two assumptions: 1) conditional independence of features given the class label, and 2) independence of features. However, the math does not seem to be working out (I get greater than 1 probability for certain labels). Let's assume we are dealing with two features ($F_1$ and $F_2$). This is the probability I want to compute: $$P(C|F_1,F_2)$$ Where $C$ is the class. By Bayes rule: $$P(C|F_1,F_2) = \frac{P(F_1,F_2|C)P(C)}{P(F_1,F_2)}$$ Using the independence assumptions above: $$P(C|F_1,F_2) = \frac{P(F_1|C)P(F_2|C)P(C)}{P(F_1)P(F_2)}$$ Now, let's say we train the Naive Bayes classifier on the following data: ![enter image description here](http://i.stack.imgur.com/5BT1x.png) And we now want to classify a new observation $F_1=1$ and $F_2=1$. So let's 1st compute $P(C=A|F_1=1,F_2=1)$: $$P(C=A|F_1=1,F_2=1)=\frac{P(F_1=1|C=A)P(F_2=1|C=A)P(C=A)}{P(F_1=1)P(F_2=1)}=\frac{1*1*\frac{1}{2}}{\frac{1}{2}*\frac{1}{2}}=2$$ Clearly, I have gone wrong somewhere. However, I can't pinpoint it. Any insights would be highly appreciated!


title: ARIMA - SARIMAX modelling with R

text: I am really new to R and to time series. My field of studies is in the field of Networks and Telecommunication, but my summer internship is about trying to find a statistical model for some sets of data. The data consists of what is called "10-minutes-points", recorded over a year and which represent power consuption of a source substation. It means I have 6 * 24 * 365 = 52 560 points of data to process, one set for each source substation. It's been about a week I'm trying to found information about ARIMA models. This website and the report of my predecessor quite helped me getting in the subject, by I still encountered many problems. I found one might be due to the large size of the data set as explained here, the second one to the existence of exogenous data as mentioned there. My predecessor found the ARIMA model to be effective for short term predictions (up to 20-ish hours), and the SARIMAX for mid-term predictions (around a dozen days). I guess it is cause exogeneous data doesn't affect as much the core data on such short periods of time. I found this thread to be very interesting but I'm not sure I understand everything. In a first time I would like to know if my understanding of the general method to evaluate a model is correct :   1. first you plot your data and try to look for any trend/seasonality (the data I have showed to have a daily seasonality and a yearly one)   2. you use log in order to reduce the trend, and maybe differentiate to eliminate the seasonality (so I should use something like : `diff(data.ts, 144)` in my case to get rid of the daily seasonality (6*24 points a day) ?)   3. plot the acf/pcf of the differentiated time series and try to estimate a model from there   4. try to fit the model to my data with `fit <\- Arima(data, order=c(p,d,q), seasonal=c(P,D,Q))` but I don't where the seasonality (144) would appear in this function ?   5. study the residuals of fit to see if the model is correct (looking at the acf/pacf)   6. use fitted or forecast (I don't know which one is better) to predict future values Thing is, since the data set is huge, I always get significant spikes at many lags in the acf/pacf and I don't feel I can judge if a model is correct or not. Here is an example : `data = scan("auch.txt", skip=1) plot.ts(data)` ![Data](http://i.stack.imgur.com/weVCX.png) `data.ts = ts(log(data) data.diff = diff(data.ts, 144) plot.ts(data.diff)` ![Datadiff](http://i.stack.imgur.com/Ck7mu.png) Which seems somehow stationary to me. I then proceed to look at the acf/pacf, and had to differentiate once more because it wasn't stationary in fact : `tsdisplay(data.diff, lag.max=150) tsdisplay(diff(data.diff), lag.max=150)` ![Tsdisplay](http://i.stack.imgur.com/dgqG6.png) ![Tsdisplaydiff](http://i.stack.imgur.com/P7Kj7.png) And I really don't know how to handle these results, so I hoped I could find some help here, because I came across the website a lot during my researchs. Thanks in advance, and I apologies for any grammatical mistakes or vocabulary error ; English is not my native language. **Edit :** does anyone know why my pictures won't appear ? **Edit bis :** nvm in fact it might be me, because imgur is blocked on my work computer


title: Help with MLE regression

text: I have a data set containing two variables `x` and `y`. I want to estimate the parameters for a regression model. The regression model is: ![enter image description here](http://i.stack.imgur.com/hsL5M.png) The maximum likelihood function that estimates the parameters is: ![enter image description here](http://i.stack.imgur.com/eNOUu.png) Here is a link to the data set: http://www.filedropper.com/data_14 Is there a way to estimate the parameters of this function? I've tried and spent many hours with R but can't get it to work. Matlab will do just fine as well. Thanks in advance!


title: Grouping a "other" group with R and ggplot2

text: I'm trying to figure out, how to best deal with an "other" group or how to create one. I'm using R and ggplot2. # Data As an example: I'm analyzing the german Bundestagswahl (general election, every four years). There are at least 38 parties, that voters could vote for in 2009 (see bundeswahlleiter.de for numbers). Let's say, that I have an mysql database `bundestagswahl` with a column `party` and a column `votes`. Now I want to do some plots using ggplot. Normally I would write a sql query with votes per party. Next I would tell ggplot2 to draw a geom_bar. This would draw me at least 38 bars. ## Now my question How can I create an "other" group, e.g. all parties, which have less then 5% of all votes to not clutter up the graph? This way I would have only 5 bars (=less clutter). Of course I could write my sql-query to ignore parties, that have not at least 5%, but this way my counts wouldn't sum up to 100%. Over 6 % of all votes went to parties, which didn't made it into the parlament. How would you do it? Write another sql query? Or using reshape? Or some other way? ## Update To make it a bit clearer, what I want to do, here is a fake dataset of six parties (names A, B, C, D, E and F) and there number of votes:               > party <- c("A", "B", "C", "D", "E", "F")     > votes <- c(1000, 2000, 5, 1500, 15, 20)     > df <- data.frame(party, votes)     > df       party votes     1     A  1000     2     B  2000     3     C     5     4     D  1500     5     E    15     6     F    20      Now I want to do two things:   1. relabel all parties, which have less then 100 votes to "other".   2. add up those "other" parties to one entry. Wanted outcome:                 party votes     1     A  1000     2     B  2000     3 other    40     4     D  1500      How do I do this? Is maybe `combine_factor` of the packages `reshape` an option?


title: Is it correct to measure R coefficient for LOG LIN models

text: I have a model where the $y$ is very skewed and I convert it to log and run a log linear model. But, I have doubts about the way to measure the error, because in the original variables the error would be much bigger than in log variables. Are there any other ways to measure the error? I would prefer to use GEP instead of log linear. But the idea is to find simple coefficients to explain. Is there another way to do it?


title: Estimating Diversity of operon types using HMMs across metagenomes: Mann-Whitney? Kruskal-Wallis? or other?

text: Question that seems complex but probably isn't. I'm comparing 10 metagenomes (http://en.wikipedia.org/wiki/Metagenomics), that are supposedly unrelated (they are 10 geographically separate locations), for a specific group of operons (http://en.wikipedia.org/wiki/Operon). I'm using Hidden Markov Model protein profiles (HMMs) to find the specific genes in each of my metagenomes. The results are counts (hits) of the fragments of the operon to the HMM, which I assume to be representative of the protein. Some Proteins have more than 1 HMM. So, I have a list of HMM's and counts per Location. _**_** _*_ **** you can start reading from here if metagenomics is a mystery to you as Statistics is a mystery to me **_*_ ** _*_ **** To simplify : I'm comparing mug factories, there are Y different types of mugs, some mugs have 10 colours some have 5 some have 3 etc. I can detect the Chips of broken mugs, one chip at a time. each chip is of a different colour. I want to compare the diversity of mugs at each factory, to say one factory shows more diversity in mugs than the rest... or two of them show tendencies towards preferring 5 coloured mugs etc (A mug that has 5 colours). The snag is I'm not seeing the mugs as a whole, I'm only able to detect 1 chip at a time. So my counts are the chips. For example: for the three coloured mugs each chip is of a different colour, one 5 coloured mug has 5 different colours, another 5 coloured mug has 5 colours, 4 of which are the same as the first 5 coloured mug (so 6 possible colours for the 5 coloured mug) all different from the three coloured mug. The 10 coloured mug has 13 possible colours (also different from the last two mugs), etc (there are many types of mugs in any one of these factories- about 130 of them). So ... The way I've attempted to approach it is as follows:   * for each theoretical mug I grouped the colours took counts and averaged them according to the possible number of colours.    * For each factory I calculated the Shannon index H using the average (ignoring any whole mugs that gave me zeros)   * Compared the H's to each other...... (probably not correct to do so) After some reading, I came across the Mann-Whitney-U-test, In theory I can compare each 2 factories to each other, but instead of averaging them I use the colours as separate sample units (the colours become n s), and rank them? , So to compare the 10 factories I have to do the test for factory 1 to 2, 1 to 3, 1 to 4, etc... (Problem is some of my colour observations have less than 2 counts, many are 1 or 0) I could use the Kruskal-Wallis test (But I'm not sure about how there are different n's for each location), this will tell me that the factories are different, but not which is more different from the other(?right <---how do I find this out) What approach do you guys think I should take? Is there a better alternative I've not found? I'm not a statistician, and I'm doing all this in Excel, NOT familiar with R at ALL. Hope my question isn't too convoluted. Thanks.


title: What tool can we use to assess if there are enough data?

text: Imagine some paired data ${(x_i,y_i)}_{i=1}^n$ representing the results of two different measurement methods and the question is about the quantification of the bias between the two methods. Let us assume that the normality assumption is reasonable for $x_i-y_i$ hence we simply compute an estimate and a confidence interval of the mean difference. Of course the width of the confidence interval is an indicator of how precise is the quantification of the bias, but what other tool(s) can we use to assess if there enough data for the quantification to be reliable ? I have in mind a "cross-validation" (I quote because I am ignorant on this topic): for example we could assess whether the result changes when we drop a data value. Is there a standard way to perform such an assessment ? I am also interested in performing this assessment in the R software.


title: p-value vs. confidence interval obtained in Bootstrapping

text: I ran a simple randomized experiment with 1 control and 1 experimental condition (final N = 80). The dependent variable (frequency of a shown behavior) is clearly not normally distributed so I thought about bootstrapping my analysis (an independent t-test). The t-test _without_ bootstrapping resulted in a significant effect between the conditions (p < .05). Using SPSS, the p-value based on bootstrapping (5000 resamples) is only marginally significant (p < .10). However, the 95% confidence interval does not include zero. I just apply statistical methods and sometimes I don't know what is actually right (unfortunately!). But that's why I ask this question. When I tried to learn how Bootstrapping works, I thought that one has to look at the confidence intervals to detect whether an effect is not zero. In my example above, the 95% CI does not accompany with the bootstrapped p-value. So I don't know whether I should report the bootstrapped CI and/or the bootstrapped p-value and/or the typical unbootstrapped p-value. What would you say?


title: Random forest models and external validation

text: How can a random forest model be published so that it also can be externally validated? Is it possible to get an average of the node split values from the ensemble?


title: Intuition for higher moments in circular statistics

text: In circular statistics, the expectation value of a random variable $Z$ with values on the circle $S$ is defined as $$ m_1(Z)=\int_S z P^Z(\theta)\textrm{d}\theta $$ (see wikipedia). This is a very natural definition, as is the definition of the variance $$ \mathrm{Var}(Z)=1-|m_1(Z)|. $$ So we didn't need a second moment in order to define the variance! Nonetheless, we define the higher moments $$ m_n(Z)=\int_S z^n P^Z(\theta)\textrm{d}\theta. $$ I admit that this looks rather natural as well at first sight, and very similar to the definition in linear statistics. But still I feel a little bit uncomfortable, and have the following **Questions:** 1. **What is measured by the higher moments** defined above (intuitively)? Which properties of the distribution can be characterized by their moments? 2. In the computation of the higher moments we use multiplication of complex numbers, although we think of the values of our random variables merely as vectors in the plane or as angles. I know that complex multiplication is essentially addition of angles in this case, but still: **Why is complex multiplication a meaningful operation for circular data?**


title: References for how to plan a study

text: In an average (median?) conversation about statistics you will often find yourself discussing this or that method of analyzing this or that type of data. In my experience, careful study design with special thought with regards to the statistical analysis is often neglected (working in biology/ecology, this seems to be a prevailing occurrence). Statisticians often find themselves in a gridlock with insufficient (or outright wrong) collected data. To paraphrase Ronald Fisher, they are forced to do a post-mortem on the data, which often leads to weaker conclusions, if at all. I would like to know which references you use to construct a successful study design, preferably for a wide range of methods (e.g. _t_ -test, GLM, GAM, ordination techniques...) that helps you avoid pitfalls mentioned above.


title: How to get around Numerical Overflow in Stata?

text: I am working with a large data set (approximately 50K observations) and trying to running a Maximum likelihood estimation on 5 unknowns in Stata. I encountered an error message of "Numerical Overflow". How can I overcome this? I am trying to run a Stochastic Frontier analysis using the built in Stata command "frontier". The dependent variable is log of output and the independent variable is log of intermediate inputs, capital, labour, and utlities.


title: Is mutual information invariant to scaling, i.e. multiplying all elements by a nonzero constant?

text: I am trying to compare it to Euclidean distance and Pearson correlation


title: How could I predict the results of a simple card game?

text: I want to predict the results of a simple card game, to judge on average, how long a game will last. My 'simple' game is;   * Cards are dealt from a randomised deck to n players (typically 2-4)   * Each player gets five cards   * The top card from the deck is turned over   * Each player takes it in turns to either place a card of the same face value (i.e 1-10, J, Q, K, A), the same suit (i.e Hearts, Diamonds, Spades, Clubs) or any suit of magic card (a jack)    * If the player can place a card they do, otherwise they must take a card from the deck    * Play continues in turn until all but one player has no cards left I'm guessing that I could write code to play a mythical game and report the result, then run that code thousands of times. Has anyone done this ? Can they suggest code that does a similar job (my favoured language is R, but anything would do) ? Is there a better way ?


title: Numerical approximation of percentiles from arbitrary pdf

text: **Given an easily-computable probability density function $f(x)$, what algorithm can we use to numerically approximate percentiles?** For instance, we might be looking for $x$ such that given $X \sim f(x)$, $p(X \le x) = 95\%$. * * * Here is a use case that motivates this: We're estimating the shape $\alpha$ and inverse scale $\theta$ of a gamma distribution (from a sequence of samples). Via Wikipedia we find a conjugate prior of the following form, originally developed by Miller (1980): $$ f(\alpha, \theta) \propto \frac{\theta^{v \alpha - 1} p^{\alpha - 1} e^{-s\theta}}{\Gamma(\alpha)^n} $$ This somewhat unwieldy function represents our belief about likely values for $\alpha$ and $\theta$. We are only interested in $\theta$, so we find its pdf by integration: $$ g(\theta) \propto \int f(\alpha, \theta) d\alpha $$ Miller (1980) suggests that we need numerical integration to compute this function. Now we can easily plot this $g(\theta$) function by repeatedly approximating it for many $\theta$ values. From this density plot we can visually glean the most probably region for $\theta$. So far so good. But we also want to compute the $\theta$ values for the 5th and 95th percentiles, $G^{-1}(0.05)$ and $G^{-1}(0.95)$, where $G^{-1}(y)$ is the inverse cdf. So my question is, how do we approximate those given $g(\theta)$? There must surely be a way to numerically approximate this by harnessing some well-known algorithm, but I cannot find a way to express $G^{-1}$ in terms of $g$. Am I missing some intermediate step to make it work? The best I can come up with is getting a lot of samples for $g$ and summing them. References: Miller, Robert B. "Bayesian Analysis of the Two-Parameter Gamma Distribution" _Technometrics_ , 1980, 22(1), 65-69 I haven't done math in a while, so it's all rusty -- if there's mistakes or misphrasings in this question, please edit away!


title: comparing a variable (that is repeatedly measured) between groups in R

text: I have a rather simple question (I hope) but for some reason I can't find the right/ a working answer anywhere on the web. For my research I have measured mood on a scale 1-10, 6 times a day for 5 weeks. (thus a repeated measures design if i am not mistaken). The participants were assigned to three conditions based on anxiety scores (n= 15, 15, 5 respectively). What I want to do is comparing mood between the three groups to see if they differ. I know I can't use a normal anova because the mood measurements aren't independable. and because I have over 6000 entries even the smallest difference is significant. I also tried calculating the mean mood of each particpant and then conduct the anova using just the mean mood of the participants, but then I lose so much. What would be the correct way to test for group difference. By the way: I am using R and the data frame is in a long format (rows = individual mood entries, column = mood rating and group) I hope someone can help me out.


title: Kruskal-Wallis, rank data

text: I want to perform a Kruskal-Wallis test, but I can't figure out whether this test will automatically rank my data or whether I have to rank it myself before I run the test.


title: Determine sample size without knowing standard deviation

text: I want to calculate the required sample size for a call centre with a finite population. The survey will collect data on many variables. The objective is to gain data about a range of psychological aspects as well as durations, agent performance, call resolution etc. Therefore I cannot choose a sample size and calculate the standard deviation with a pilot, because I wouldn't know what to base the deviation on. I have read that an industry standard deviation can be used, however I have no such info. So if I cannot calculate it using a pilot, do not know the industry st.d should I guess the standard deviation in order to calculate the sample size? If so, how can I make an educated guess? Could I instead calculate the sample size using proportions? as described in (http://www.polarismr.com/research-lifeline/sample-size-calculator/) In this link (http://www.surveysystem.com/sample-size-formula.htm) there is a formula for calculating it ignoring the standard deviation. Is this recommendable? I am very grateful for replies. Please keep in mind that I am new to statistics and I might have overseen some important info above.


title: Parameter optimization of SVM

text: Currently I am using SVM to perform some classification task. I use libSVM with Matlab interface. From the practical guide of SVM (Link), we know that there are two parameters need to be tuned, namely L and C. L is the lambda parameter of the RBF kernel and C is the penalty parameter of the error term. Grid search is widely used and suggested by the practical guide. However, it is not clear to me when should I do the grid search to optimize the two parameters. For example, I want to apply a leave-one-out cross validation scheme to evaluate the SVM performance. I have got 50 data. And some pseudo-codes of the leave-one-out scheme shown as below:               load data;     nFold = 50;  %% leave-one-out for 50 data          %% Should I do the grid search here with part of the data?          for i=1:nFold            %% Partition the data into training (n=49) and testing (n=1)       dataPartition = cvpartition(data);             %% Or should I do the grid search here? Optimize the two parameters for each fold?            %% Build the SVM model from 49 training data       SVM_Model = svmTrain(dataPartition(training));            %% Apply the model to predict the labels of the testing data       SVM_Classified = svmClassify(dataPartition(testing), SVM_Model);          end      As shown in the pseudo-codes above, where should I apply the grid search? Options: A. Randomly choose a subset of my data for example 30 of 50. Separate my 30 data randomly into 80% training and 20% testing and do the grid search to optimize the parameters. B. Do the grid search within each cross-validation fold. Which way will have the less biased results? Or if you have a better idea to do this please tell me. Thanks a lot. A.


title: Addressing multicollinearity with key driver analysis

text: I am trying to determine the key drivers from a series of 30 Independent Variables (IVs) (attributes rated on 10 pt scale) on 3 Dependent Variables (DVs) (i.e. purchase intent). The 30 IVs are pretty highly correlated, but the Variance Inflation Factors (VIFs) are all under 10 (though most are over 5). I did a backward regression, but some of the beta coefficients had negative signs, which don't make sense logically. Is it better to just run separate simple regressions on each IV, or take out the IV's with particularly high VIFs and run a regression with all of the remaining in the model together? Any suggestions would be appreciated!


title: Finding stationary distribution of Markov chain with high number of states

text: I have transition matrix for the Markov chain of large size. Each entry of the matrix involves some complicated summations and products. I would like to find stationary distribution and analyze some other parameters the system. Is there any way to do the above tasks in simpler manner.


title: Removing seasonality from data

text: I have a dataset depicting weekly revenue over time for a computer company. The plot for the data looks like this: ![enter image description here](http://i.stack.imgur.com/DDoy8.png) I decomposed the data into its additive components using the `decompose` function in R and plotted the various components: ![enter image description here](http://i.stack.imgur.com/4HbER.png) Next I tried removing the seasonal component using the following code:               > RevenueDec <- decompose(Revenue)     > RevenueSeasonallyAdjusted <- Revenue  - RevenueDec$seasonal      However, I still get a seasonal component when I decompose 'RevenueSeasonallyAdjusted': ![enter image description here](http://i.stack.imgur.com/CC9V9.png) The y-axis has very small values but the seasonality exists nonetheless. Could you help me out here. EDIT: In the next step I tried using the `auto.arima` function on my seasonally adjusted data to get a forecast and got a plot like this: ![enter image description here](http://i.stack.imgur.com/1BPbL.png) Is this the correct approach to use, or should I try something different?


title: Clustering Two Variables With Disease Information

text: I was proposed a problem and I am not quite sure how to go about it. The problem is I want to find a relationship between two variables. For the simplified case there are only two variables, lets say Diabetes and Congestive Heart Failure. Each variables has a $0$ or a $1$, if you do not have the disease or if you do, respectively. This was presented to me as a clustering problem. We can see how the clusters form and see which disease is most dominant or just how they relate to each other. My assumption is there would be 4 clusters, cluster 1 is people who have neither, cluster two is people with Diabets, cluster 3 is people with Congestive Heart Failure and cluster 4 is people who have both. My questions are:   1. Is this a valid way to determine a relationship between the variables?   2. How would I create the affinity matrix?   3. How would I visualize the clusters? Has anyone has experience with this type of modeling? Any ideas outside of clustering would be appreciated! Any insight or help or references would be highly appreciated! Thank you in advance!


title: What is the distribution of a sum of a subset of probabilities, with each probability having the same distribution?

text: Suppose I have $k$ outcomes with probabilities, $p_i$, with $p_1+p_2+\dots+p_k=1$. Each probability has the same distribution. What would the distribution of a sum of probabilities be? For example, what would the distribution of $p_1+p_2$ be? The Irwin-Hall distribution is the distribution for the sum of uniformly distributed random variables. Each $p_i$ will be distributed between 0 and 1, but they won't be uniformly distributed. Does the Irwin-Hall distribution have anything to do with it? I don't require a fully derived answer, but a kick in the right direction would help.


title: Linear Model or Logistic Model: Can Someone Recommend a Book?

text: I have a huge data set that looks roughly like this:               x = [x1, x2, x3, ..., x800]     x is how much percent the real project cost is different from the planned project cost          y1 = [y1.1, y1.2, y1.3, ..., y1.800]     y1 is the estimated project cost in dollar (or another currency)          y2 = [y2.1, y2.2, y2.3, ..., y2.800]     y2 is the real projet duration          y3 = [y3.1, y3.2, y3.3, ..., y3.800]     y3 is the project manager (in coded form like "mk" and "op")          ....          y50 = [y50.1, y50.2, y50.3, ..., y50.800]     y50 is the estimated project duration      Now the aim is to predict the `x` (cost difference in percent) in dependency of other `y` variables. It is a priori not clear which variables influence it. And some of the `y` are not interdependent like for example the difference in percent `x`, the real costs and the estimated costs since `x` is calculated from the real and estimated costs. I am now wondering how is the _best practice_ way to analyse the data or how to build a model (linear or maybe logistic). The binary output for the logistic model could be **difference bigger than 10 %**. Since I have a lot of data sets (over 800) almost everything could get significant on a e. g. 5 % level. And It's not easy to handle/control the power if you do not have a simple t-test. What would you recommend? Is there a good book for "normal people" about that kind of statistics?


title: Multivariate normal distribution has peaks

text: I'm trying to calculate a bivariate normal distribution in matlab(with mvnpdf), but the pdf I obtain has a strange shape with several peaks. ![Results obtained](http://i.stack.imgur.com/vMLuj.png) This is the code I use:               s=[3,1]';     C_ee = [0.0473   -0.1446;            -0.1446    0.4440]          x1 = -5:.2:5; x2 = -5:.2:5;     [X1,X2] = meshgrid(x1,x2);     F = mvnpdf([X1(:) X2(:)],s',C_ee);     F = reshape(F,length(x2),length(x1));     figure(1);     surf(x1,x2,F);     caxis([min(F(:))-.5*range(F(:)),max(F(:))]);     axis([-5 5 -5 5 0 5])     xlabel('Re'); ylabel('Im'); zlabel('Probability Density');      I've noticed, that when I increase the value of C_ee (sigma matrix), for instance C_ee+0.05 the shape starts to look normal. I notice that this matrix is close to not being positive definite... but it still is. Could anyone explain this behaviour? Thank you for your time.


title: Forecasting an MA(1) process

text: Suppose $x_{t} = w_{t} + \theta w_{t-1}$ where $w_t$ is white noise with variance $\sigma_{w}^2$ Derive the minimum mean square error one-step forecast based on the infinite past and determine the mean square error of this forecast. Let $\tilde{x}_n^{n+1}$ be the truncated one step ahead forecast based on the n previous observations. Show that $E[(x_{n+1}-\tilde{x}_n^{n+1})^2=$ $\sigma_w^2(1+\theta^{(2+2n)})$ The solution for the first part states that: $x_{n+1} = \sum_{j=1}^{\infty}-\theta^{j}x_{n+1-j} + w_{n+1}$ and $\tilde{x}_{n+1}=\sum_{j=1}^{\infty}(-(\theta)^{j}x_{n+1-j})$ so the $MSE = E[x_{n+1}-\tilde{x}_{n+1}] = \sigma_{w}^2$. I am not sure how the expression for $x_{n+1}$ or $\tilde{x}_{n}^{n+1}$ was arrived at. Can someone explain how this was derived? I know $x_{n+1}$ is a linear process but I am not sure where the $-\theta$ plays a role. For the second part, the solution states $\tilde{x}_{n}^{n+1}= \sum_{j=1}^{n}-\theta^{j}x_{n+1-j}$ and $MSE = E(x_{n+1} - \tilde{x}_{n}^{n+1})^2 =E[\sum_{j=n+1}^{\infty}-\theta^{j}x_{n+1-j} + w_{n+1}]^2$. I am not sure how this was arrived at.


title: Estimate SVM a posteriori probabilities with platt's method does not always work

text: I have a problem.. I'm trying to create a multiclass SVM with probability output. The SVM is working so far, what means, that the accuracy is ok (see the last picture). But the probability estimation does not work so far. I use the one vs one multiclass encoding. Now as Platt (http://en.wikipedia.org/wiki/Platt_scaling) suggests, I'm fitting a sigmoid in each binary model by minimizing the negative log likelihood $$ \min_{A,B} -\sum_{i=1}^N t_i\log(p_i)+(1-t_i)\log(1-p_i),$$ where $ p_i = \frac{1}{1+\exp(Af_i+B)}$, $N$ the Number of samples and $f_i$ is the uncalibrated output (the score) of the SVM. These values are the leave one out crossvalidation values to prevent overfitting. But there is no much difference between trainings data, 10 fold cv and loo cv. The improved algorithm by  Lin et all works well, I can compare it to matlab's fminfunc. The problem: The fits don't fit sometimes, see the pictures below. There are 3 times 2 pictures. The first shows the featurespace of the binary SVM (for the one vs one coding) with the uncalibrated output as contour lines. The second picture shows a subplot of the probability density of the svm output for each class (estimated via ksdensity in matlab), the second subplot shall be the posteriori possiblity. They were calculated via Baye's theorem. There might be an error in the calculation, so I write it down: $$ p(y_1|f)=\frac{p(f|y_1)p(y_1)}{\sum_{j=1}^{\text{all classes}} p(x|y_j)p(y_j)} $$ where $y_i$ is the predicted ith class and $f$ is again the svm output. I can calculate the prior $p(y_j)$ simply by counting the positive and negative classes, right? In this case it's $0.5$ My problem is the likelihood. As far as I have understood in the bayesian case the likelihood is a probability density that shows us the distribution of the svm output by given class, correct? I thought this is exactly I'm showing in the first subplot, or am I wrong? My code follows:                       [pPlus,xi] = ksdensity(score(yBinary == 1),[-1.5:.0001:1.5]);             [pMinus,xi] = ksdensity(score(yBinary == -1),[-1.5:.0001:1.5]);                  evid = pPlus* Nplus/(Nplus+Nminus) + pMinus* Nminus/(Nplus+Nminus);             postPlus = pPlus *  Nplus/(Nplus+Nminus)./(evid); % likelihood * prior / evidence             postMinus = pMinus* Nminus/(Nplus+Nminus)./(evid);      Nplus is the number of positive classes, Nminus the number of negative. score is the SVM output, xi is simply the interval on the x-axis. pPlus is the positive likelihood. Is there an error in the formulars or do you need further code/data? I ask because there are some oscillations in the second picture. And in the 4th picture the sigmoid doesn't fit at all, whereas the sigmoid in the 6th picture is ok. Are the subplot pictures plausible, or might there be an error? Or is just the method not working? As I said above, the optimization works and gives different values for each sigmoid (though it might be impossible to see here), but fminfunc and the proposed algorithm support the same values. By the way the last picture is showing the predicted classes of the multiclass SVM, maybe you can understand the pictures of the binary svms better. If you have any ideas, you're welcome! Thank you very much. ![Binary SVM class 1 vs 2 with decision values](http://i.stack.imgur.com/oA69l.jpg) ![enter image description here](http://i.stack.imgur.com/Vtz4y.jpg) ![Binary SVM class 1 vs 3 with decision values](http://i.stack.imgur.com/YTcvm.jpg) ![enter image description here](http://i.stack.imgur.com/wrmbF.jpg) ![Binary SVM class 1 vs 4 with decision values](http://i.stack.imgur.com/cNTjy.jpg) ![enter image description here](http://i.stack.imgur.com/QwI52.jpg) ![featurespace result](http://i.stack.imgur.com/BN6re.jpg)


title: F statistic, F-critical value, and P-value

text: I'm very new to this area and am having difficulty understanding the concept of rejecting the null hypothesis based upon results from the ANOVA table. How do the calculated F and the critical value relate to the p-value? And if the calculated F is greater than 1, does that always indicate that the null hypothesis should be rejected, even if the p-value is less than the alpha? Sorry if these questions are signs of my ignorance, but I'm 57 and returning to school after a 35 year absence! Thanks for any help.


title: How close to zero should the sum of the random effects be in GLMM (with lme4)

text: I'm using the `lme4` package in R to do some logistic mixed-effects modeling.   My understanding was that sum of each random effects should be zero. When I make toy linear mixed-models using `lmer`, the random effects are usually < $10^{-10}$ confirming my belief that the `colSums(ranef(model)$groups) ~ 0` But in toy binomial models (and in models of my real binomial data) some of the random effect sum to ~0.9. Should I be concerned? How do I interpret this? Here is a linear toy example `               toylin<-function(n=30,gn=10,doplot=FALSE){      require(lme4)      x=runif(n,0,1000)      y1=matrix(0,gn,n)      y2=y1      for (gx in 1:gn)      {        y1[gx,]=2*x*(1+(gx-5.5)/10) + gx-5.5  + rnorm(n,sd=10)        y2[gx,]=3*x*(1+(gx-5.5)/10) * runif(1,1,10)  + rnorm(n,sd=20)      }      c1=y1*0;      c2=y2*0+1;      y=c(t(y1[c(1:gn),]),t(y2[c(1:gn),]))      g=rep(1:gn,each=n,times=2)      x=rep(x,times=gn*2)      c=c(c1,c2)      df=data.frame(list(x=x,y=y,c=factor(c),g=factor(g)))      (m=lmer(y~x*c + (x*c|g),data=df))      if (doplot==TRUE)       {require(lattice)        df$fit=fitted(m)        plot1=xyplot(fit ~ x|g,data=df,group=c,pch=19,cex=.1)        plot2=xyplot(y ~ x|g,data=df,group=c)        print(plot1+plot2)       }      print(colMeans(ranef(m)$g))      m     }      In this case the colMeans always come out $<10^{-6}$ Here is a binomial toy example (I would share my actual data, but it is being submitted for publication and I am not sure what the journal policy is on posting beforehand):               toybin<-function(n=100,gn=4,doplot=FALSE){       require(lme4)              x=runif(n,-16,16)       y1=matrix(0,gn,n)       y2=y1       for (gx in 1:gn)       { com=runif(1,1,5)         ucom=runif(1,1,5)         y1[gx,]=tanh(x/(com+ucom) + rnorm(1)) > runif(x,-1,1)         y2[gx,]=tanh(2*(x+2)/com + rnorm(1)) > runif(x,-1,1)       }       c1=y1*0;       c2=y2*0+1;       y=c(t(y1[c(1:gn),]),t(y2[c(1:gn),]))       g=rep(1:gn,each=n,times=2)       x=rep(x,times=gn*2)       c=c(c1,c2)       df=data.frame(list(x=x,y=y,c=factor(c),g=factor(g)))       (m=lmer(y~x*c + (x*c|g),data=df,family=binomial))       if (doplot==TRUE)        {require(lattice)         df$fit=fitted(m)         print(xyplot(fit ~ x|g,data=df,group=c,pch=19,cex=.1))        }       print(colMeans(ranef(m)$g))       m     }      ` Now the colMeans sometimes come out above 0.3, and definitely higher, on average than the linear example.


title: Downsampling excessively sampled curves

text: I'm looking for a quick (event if not so accurate) way to deal with over- sampled curves using R. Consider the following example in which `x` contains 1000 values in the range [0, 1) and `y` is some function of `x`.               x <- sort(runif(1000))     y <- sin(x)      I'm looking for a way to take only 100 points on this curve, preferably in equal intervals, without knowing the actual model that produces `y`'s from `x`'s. My naiive approach is to find indices of original `x` values that are closest to the new ones:               x.new <- seq(0, 1, 0.01)     y.new <- c()     for(x.val in x.new){       ix <- which.min((x.val - x) ** 2) #closest value in x       y.new <- c(y.new, y[ix])     }      Or, using a `sapply`               find.closest.y <- function(v){       return(y[which.min((v - x)**2)])     }          y.new <- sapply(x.new, find.closest.y)      I'm OK with the accuracy of this method, but I suspect that there is a way to acheive this using functions from standard R libraries. Am I right? Is there a simpler method to downsample curves?


title: Problem with estimating probability using the multivariate Gaussian

text: I'm working on a problem where I need to find the probability of a given data sample belonging to a give class using the Bayes Theorem for classification. From everything I've been able to find so far the way to do this would be to first calculate the conditional probabilities using multivariate Gaussians (shown below), then apply Bayes Theorem to get the desired posterior probability. $$P(x|k)=\frac{1}{\sqrt{(2\pi)^{p/2}|\boldsymbol\Sigma|^{1/2}}} \exp\left(-\frac{1}{2}({x}-{\mu_k})^T{\boldsymbol\Sigma}^{-1}({x}-{\mu_k}) \right)$$ The problem I'm having is that the Mahalanobis distance in the exponent term is consistently coming to extremely large values, making the conditionally probabilities far too small to work with programmatically. I've determined this is almost entirely due to the covariance matrix, which is populated with values that are incredibly small. Furthermore, when I just ignore the covariance matrix and solve using the Euclidean distance instead, I arrive at reasonable probability estimates with distances that appear to more accurately reflect what I would expect in classification. So I was wondering:   1. If there was a way to calculate these conditional probabilities or the overall posterior probability without direct exposure to the exponent, so they become feasible to work with programmatically.    2. If there is any statistical justification for using the Euclidean distances over the Mahalanobis distance in this case. * * * **[Edit]** Here are a few more details regarding the derivation of the covariance matrix and the dataset: The dataset I'm working with has been transformed into a lower dimensional space using the discriminant analysis (ULDA) technique described in this paper. In this form of discriminant analysis there is a shared pooled within- class covariance matrix corresponding to each class conditional probability (this being calculated over the reduced dimensional space). Interestingly, with this technique, if the total scatter matrix is used as the shared covariance matrix then the covariance matrix in the reduced space will be the identity matrix. This technique optimizes with respect to the total scatter matrix rather than the within-class scatter matrix, so maybe I'm misunderstanding what the covariance matrix should be in this case. Also, as mentioned in the paper, the ULDA technique produces features in the reduced space are uncorrelated to each other. I suppose another possibility would be that this particular dimensionality reduction reduces the usefulness of the within-class scatter matrix to the point that it gives worse results than just using the plain identity matrix as the covariance. * * * Any help or advice would be greatly appreciated, thanks in advance.


title: Clustering Variable

text: I have a set of 10 variables (v1..v10) which are continuous. I have another two control variables ctrl11 and ctrl12 which of course are categorical. ctrl11 can take any values from 1 to 50 and ctrl12 can take any value from 1 to 200. I am trying to shrink ctrl11 and ctrl12 to about 5-10 categories rather than the current 50 and 200 respectively. The end goal is to use the shrunk version of one of these two in a logit regression along with v1..v10, ie a total of 11 DVs. So my idea is to see how the variables v1..v10 cluster. Is this strategy sensible? How can I find which control variable they cluster along so that I shrink that variable and discard the other for use in the regression?


title: Algorithms that use multiscaler properties of data to cluster

text: I was thinking of devising a clustering algorithm (for fun and kicks) that would cluster data by looking at the distribution of the data at multiple scales. For example say my data was distributed on a 2d grid of 1000 x 1000. Are there algorithms out there that cluster this data by looking at the data by dividing this space up into say   1. 10x10   2. then 20x20   3. then 40x40   4. and so forth Would appreciate links to pseudo code or implementations in R, Java, Python, matlab -- preferably open source. Additional notes:   1. I am not looking for hierarchical clustering where I find clusters within clusters   2. I would be interested in other definitions of multiscalar


title: How to find almost identical records in dataset using an algorithm (for plotting purposes)?

text: There are 10 datasets, each contains around 2000 records with 21 features and a binary label.   There have been many attempts to clean the data and even showing that there were so many redundant, duplicate and almost identical records in the original version of dataset. In all versions of cleaned data, more than 50% of original records are removed which half of them are marked as duplicate or almost identical, I want to show these duplicated records and their connection in a plot, so I was wondering if there is any method to find the closest records? For example, using some distance relevancy method like Mahalanobis Method, Euclidean Distance, or even Nearest Neighbourhood can we find very close points (records) in the dataset? (by setting a cut-point value for the distance to mark values less than that as very/almost identical)


title: How to interpret output of CVlm() in R?

text: I am using 10 fold cross validation using the `CVlm()` function from the DAAG package. This is part of the result shown:               Predicted    2.47e-04 2.26e-04 -0.000359  0.000335     cvpred      -7.88e-05 3.72e-06 -0.000597  0.000322         y            1.47e-03 2.21e-03 -0.004676 -0.001969     CV residual  1.55e-03 2.21e-03 -0.004078 -0.002291          Sum of squares = 0    Mean square = 0    n = 48           Overall (Sum over all 48 folds)            ms      9.81e-06       What is the difference between “Predicted” and “cvpred”? If I change the seed in `CVlm()` cvpred changes but Predicted remains the same. Can someone tell me how Predicted is calculated?


title: Odds ratio for no. of hypertensives in different age groups

text: I have 6 age groups 18-20, 20-30, …, >60. I want to see the significant association between age and no. of hypertensives. I used Chi-square test for this and I found significant. Now can I compute odds ratio to each group so that to get in which age group risk is more? If I can use OR how to do it ?


title: What are the assumptions for multivariate outlier detection using Minimum Covariance Determinant (MCD)

text: Currently I am working on multivariate outlier detection. I would like to try some statistical methods, such as the Minimum Covariance Determinant (MCD) method, to identify the outliers. In my data set, some of the variables do not follow certain distributions, say the Gaussian distribution. Therefore, I am wondering whether it is valid to use the Minimum Covariance Determinant method in this case, to identify multivariate outliers. In another word, what are the basic assumptions for using MCD to detect multivariate outliers?


title: How can I sample from the conditional distribution?

text: I am learning Gibbs Sampling, in which there is a step named sampling from conditional distributions. I don't understand: 1. where is the conditional distribution from? From a general case, how can I get the conditional distribution? Any help will be greatly appreciated.


title: Should I use the use the confidence-t or confidence-norm function in Excel?

text: * n=6 (this relates to maximum values of measurements of force (in Newton) when conducting the same process with the same setting of parameters 6 times)   * There is `CONFIDENCE.NORM` in Excel   * and `CONFIDENCE.T`. When I refer to all these 6 values, I am using STDEVP (or whatever exactly it is called in the English version of Excel--I got various hits for the "whole" stdev variant).   1. What should I use for calculating a CI with `\alpha = 0.05` then, `CONFIDENCE.NORM` or `CONFIDENCE.T`?   2. If I refer to 4 of said 6 values and then use `STDEV.S`, would that be correct? Also, if I like to calculate the CI (alpha=0.05) for those 4 values, I have to switch the formula to using the stdev.s value, right?


title: How to develop a prediction model based on correlation in R?

text: I have two sets of data, say sales and profit, and I have calculated the correlation between these two data over different months using `R`. So currently I have somethings like this:               Month       | Jan | Feb | Mar | Apr | ....     sales       |  X  |  Y  |  Z  |  P  | ....     profit      |  A  |  B  |  C  |  D  | ....     correlation | .1  | .35 | .28 | -.47| ....      This correlation calculation has been done over period of years and now I want to develop a prediction model which can provide me prediction for coming months of coming years,say what will be the profit for month of Aug 2014? How can I develop such a prediction model and what else, apart from what I have as of now, will I need to develop the prediction model? I am new to statistics and predictive analysis and thus if anyone could provide me a guide of required steps and some info about how can I do those steps then it would be very helpful for me to get the understanding also. Thanks in advance. **NOTE:** If I could get help in respect to some non parametric bayesian method then it will be great.


title: Binary regression accuracy vs model fit in R

text: I ran two logistic regression models, one with a dataset including outliers and one without outliers, with multiple predictors. I checked each model's fit with the le Cessie – van Houwelingen – Copas – Hosmer unweighted sum of squares test for global goodness of fit from the rms package in R (following advice here).               model1 <- lrm(y ~ a + b + c + d, data1, method = "lrm.fit", model = TRUE, x = TRUE, y = TRUE, linear.predictors = TRUE, se.fit = FALSE)     residuals(model1, type = "gof")      For the model with outliers the p value was close to 0, indicating a lack of fit. For the model without outliers p value was 0.52, indicating that my model was not incorrect. I then ran 10-fold cross validation for both models with DAAG package and was surprised to get identical (poor) accuracy results for both = 0.56               cv10<-CVbinary(model1,nfolds=10)      I thought that the model created using the dataset without outliers, having a much better fit, will give me higher accuracy. Am I missing something here? I will be grateful for your help.


title: Bias correction of ridge regression

text: I'm using ridge regression for estimating the fair value of variable $Y$ against a vector of correlated predictor variables $X_i$, based on past observations of both $Y$ and $X_i$. Let's assume that both Y and Xi are standardized (mean $= 0$, std $= 1$). Problem is that ridge regression introduces bias into the estimate of $Y$, and the further the predictor vector is from zero, the greater the bias. In order to correct for the bias, I perform another regression of $Y$ onto one variable, which is the linear sum of predictors with coefficients calculated from the ridge regression. It seems that this procedure both benefits from the stable coefficients calculated using ridge and from the subsequent bias-correction, which looks like a double win. Is this procedure "valid"? Does anybody do the same? Is there any literature on such kind of bias correction?


title: Difference-in-Difference, number of Units

text: Is there a specific rule to the number of units in the treatment and control group. I am trying to carry out a difference-in-difference estimation and in the treatment group I have about 9000 individuals and in the control group I have about 800 individuals.


title: Multiple hypothesis testing with FDR in R - FDRtool and p.adjust

text: i am very new to this field, so please bear with me. I have 600 parameters, on which i perform simple pairwise x,y correlations.For adjustment of the p-values ive been stuck between two R-packages. FDR-tool gives me a set of qvalues,               FDR <- fdrtool(as.vector(a$P[!is.na(a)]), statistic=c("pvalue"), cutoff.method=c("fndr"))      whereas the p.adjust command in the basic packages gives me adjusted p-values.               FDR <- p.adjust(as.vector(a$P), method = "fdr", n = length(as.vector(a$P)))      "a" is an object created by the hcorr command of the package "Hmisc" The output is not the same, however. Can someone please explain me why that is? And how to proceed from here? I was told, the adjusted p-values (as obtained from p.adjust) can be treated as such, that i can just define a cutoff of 0.05 for false positives. I am not so sure about q-values, though, and if i am missing something fundamental on False Discovery Rates. Thank you very much!


title: High standard deviation: rule of thumb

text: I know this is a pretty vague question but I was wondering if there was a rule of thumb relation between the mean and the standard deviation of a given set of values by which you can say "Something's not right". I'm thinking something like: _If the standard deviation is half the value of the mean --> raise a flag_.


title: How should I evaluate the value of a correlation coefficient?

text: I am doing regression on a numeric dataset. This dataset is combined from 2 sources (one is from real-life data, one is from synthetic data). When I run cross-validation on this combined data, I get a correlation coefficient of 0.9. But when I train the model on the combined data, but test it only on real-life data, I get a correlation coefficient of 0.2. I know that correlation coefficient gives the quality of least squares fitting, so in that case why would my model performs well on combined data but performs bad on real-life data (even though I trained my model on it)? I read some book sections on correlation coefficient, but intuitively I can't comment on why I am getting this result.


title: Question about a marginal distribution

text: If I observe the following: $X \sim N(\mu_x,\sigma^2_x)$ $Y|X=x \sim N(x,\sigma^2_y)$ My objective is to calculate the marginal distribution of $Y$. (Since the variance term does not address some form of correlation $\rho$, the dependence between both random variables clearly needs to be addressed w.r.t. to the expected value of $Y$) By the usual formulas we will get $f(y)$ by: $f(y) = \int_{-\infty}^\infty f(x,y)dx$ where $f(x,y) = f(y|X=x) f(x) = \frac{1}{\sqrt{2\pi\sigma_x^2}}\frac{1}{\sqrt{2\pi\sigma_y^2}}exp\{-\frac{1}{2}[(\frac{y-x}{\sigma_y})^2 + (\frac{x-\mu_x}{\sigma_x})^2 ]\}$ I know the usual bivariate case and this looks like that case, where the correlation $\rho$ is zero and solving would be easy. But here in both terms the $x$ factor needs to be addressed and I'm not quite sure how to solve for $x$.


title: Estimating kappa of von Mises distribution

text: Is there a way to calculate an estimate of $\kappa$ from data for the von Mises distribution? It seems very easy to do in R, http://rgm2.lab.nig.ac.jp/RGM2/func.php?rd_id=CircStats:A1inv, but python doesn't have an A1inv function to calculate the ratio of the first and zeroth order Bessel functions of the first kind.


title: Detecting changepoint in ratio of rates of two Poisson processes

text: I'm interested in a changepoint detection problem of the following scenario: Consider two Poisson processes for which we have the event times. I'm interested in detecting a change in the _relative_ intensities of the processes, normalizing out fluctuations in the combined/joint intensity. Furthermore, one class is much rarer than the other, and as usual, it's the rare one which is the interesting business case (this is for a real application, not academic). For example in a web commerce application, the processes might represent clicks or page hits, and the interesting ones might be ones which converted to a sale, and could be but 1% or 0.1% or 0.01% of the total. You can imagine many other problems with a similar structure, such as detecting anomalies in internet communications (cyber security), industrial control (failure events) or financial transactions (risk). I'm interested in detecting a significant change (step or better yet ramp up) in 'hit rate' (if you count that as which fraction of clicks were converted to sales vs total rate /misses) but reasonably independent of both short term common dynamics or long-term trends of the overall volume. That is, a relative rate increase of both the rare and common process to the same degree is not interesting. I'd also like to avoid significant modeling or assumptions on the overall rate as that could vary for all sorts of unknown and uncontrollable reasons. I'd rather not impose arbitrary time boundaries for quantization/summation if not necessary but I'm open to any practical approaches. Early detection of a real change is very valuable in the business context. Computational efficiency should also be important as the data volume could be very large, as there would be potentially millions of these streams (each with a two-class pair of processes) and being alerted which one may be changing structure (and an estimate of when) is quite valuable. Does anybody know of any approaches, literature or software for this problem? It seems as though it should have been addressed somewhere already given the general nature of the problem.


title: Baseline predictors parametres

text: I've implemented baseline predictors model. It trains on data: "`user_number item_number rating_ui`" And then I need to predict raiting for "`user_number item_number`". I use following formulas for biases:               b_ui = average_rating + b_i + b_u;     b_i = (sum_(by all items)(rating_ui-average_rating))/(lamda_2+items_count);     b_u = (sum_(by all users)(rating_ui-average_rating-b_i))/(lambda_3+users_count);      I tried to use Netflix constants lambda_2 = 25 and lambda_3 = 10. But it gives not enough accuracy. How can I optimize lambda_2 and lambda_3 using training data? There are to much data to use gradient descent.


title: Guides to VARMA modelling in R

text: I'm looking at using a VARMA model to both determine the driver so value in some advertising campaigns and also to forecast future activity. I'm looking at the paper by Takada and Bass as a reference point and I'm looking through the documentation for the MTS package by Tsay (the accompanying book is on order). I've noticed that there does not seem to be much online for getting started with fitting VARMA moels, are there any open source docs that anyone knows of on how to get started? (While I wait for the book to arrive)


title: Stepwise regression modeling using multiply imputed data sets

text: After multiply imputing data, it is natural to estimate regression models on the data. When multiple predictors are available, sometimes stepwise regression is used for model building (forward inclusion or backward elimination of covariates). My question is how to use a similar procedure on multiply imputed data sets (e.g. using the function `mice` in R). At first glance, this appears to be difficult. An AIC criterion as in the `step`could be used for exaple on each multiply imputed data set. But it seems hard to pool data sets after stepwise regression in the end. Maybe there is some procedure available to do this.


title: survreg Ran out of iterations and did not converge

text: Survival data with censoring were generated by monte carlo method in R. Patients survival data was assumed to follow a weibull distribution. Parameters were extracted by fitting a parametric model to the simulated survival data with survreg. However, I got warning message: In survreg.fit(X, Y, weights, offset, init = init, controlvals = control, : Ran out of iterations and did not converge. The Scale and Intercept were NaNs extracted from survreg. My R code are as follows:               library(survival)     n<-100 # 100 or 500 the number of patients      foltime<-12 # calendar time when study ends     simulations<-1000     mean_time_simulation<-numeric()     mean_time_for_Hoyle<-numeric()     mean_time_simulation_Guyot<-numeric()     censored_proportion<-numeric()     mean_time_for_Guyot<-numeric()     mean_time_ls<-numeric()     mean_time_regression<-numeric()     mean_time_complete<-numeric()     tk_max<-numeric()     for(counter in 1:simulations){     shape_input<- 2     scale_input<-11.28      mean_time <- scale_input * gamma(1+1/shape_input)   # population mean time     mean_time_censored <- 3     te<-rweibull(n,shape_input,scale_input)     tc<-rexp(n, rate = 1/mean_time_censored )   # times when patients censored     t0<-runif(n,0,3)     t_end<-t0+te     t_calendar_censored<-t0+tc     it<-c(1*(t_end<t_calendar_censored)*(t_end<foltime))     time1<-(pmin(t_calendar_censored, rep(foltime,n))-t0)*    (t_end>pmin(t_calendar_censored,rep(foltime,n)))    +   te*(t_end<=pmin(t_calendar_censored,rep(foltime,n)))     fit<-survfit(Surv(time1,it)~1)     plot(fit,conf.int=FALSE)     censored_proportion[counter]<- (n - sum(it)) / n  #desired censoring rate is 90%      results<-summary(survfit(Surv(time1,it)~1),time=,censored=T,extend=T,rmean="individual")      xaxis<-c(seq(0,100,0.1))     yaxis_actual<-exp(-(xaxis/scale_input)^shape_input)       model_complete<-survreg(Surv(time1,it)~1,dist="weibull")     lambda_complete<-exp(summary(model_complete)$table[1])         gamma_complete<-1/exp(summary(model_complete)$table[2])     mean_time_complete[counter]<-lambda_complete*gamma(1+1/gamma_complete)     bias_mean_time_complete<- mean(mean_time_complete - mean_time)     MSE_complete<-(mean(mean_time)-mean(mean_time_complete))^2 +   (sd(mean_time_complete)/(simulations^0.5))^2     }


title: Probability of Achieving a Count Level in Time Series Data

text: I have some time-series data that displays a count value for every day: ![enter image description here](http://i.stack.imgur.com/zXl9P.png) These count values begin at 1 or -1 and will continue to count up (or down) if conditions in the time series are met. If the conditions are not met consecutively, then the count value can return to -1, 0, or 1 (to initiate the next count series). 0 values do not represent the beginning of a count series, but rather an instance when none of the directional count conditions are met; functionally, an NA value. My question is, with this data, is it possible to determine the probability of a future count value being reached once a minimum threshold is met? For example, what is the confidence that I can have that a count of -3 [10/16/2006 in the above data] will then make it down to a count of -6 [10/19/2006]? I've plotted the distribution of these counts for a much longer data series. The chart below shows that a count of up to -20 was reached at least once between 01-01-2000 and 12-31-2006! With regards to my probability question, I was thinking that maybe I could sum the area between -3 and -6 to determine the probability that a -3 will reach a -6, but that doesn't feel right. ![enter image description here](http://i.stack.imgur.com/oT1RG.png) I sincerely appreciate your help!


title: Dummy variables in logistic regression vs. svms

text: Suppose $y$ is a binary outcome variable and $x$ is a categorical predictor variable that takes three levels (1,2,3). In this case, you would create two dummy variables $x_2, x_3$. So $x_2=1$ if $x=2$ or $x_2=0$ if $x \neq 2$. Likewise $x_3=1$ if $x=3$ or $x_3=0$ if $x \neq 3$ A logistic regression model would be $$\text{logit}(p) = \beta_0+\beta_{1}x_{2} +\beta_{2}x_{3}$$ with $x=1$ as the reference level. In SVMs, do we create a dummy variable for each level (i.e. $x_1,x_2,x_3$)?


title: Arima with xreg, rebuilding the fitted values by hand

text: I'm using R to do some time series estimation. I'm trying to rebuild the fitted values from an Arima model by hand to use in an Excel spreadsheet using the estimated coefficients and the input data. I can use the fitted command, but I'm trying to understand more how it works. Ex:               library(MASS)     library(tseries)     library(forecast)          set.seed(1)     N = ts(mvrnorm(50, mu=c(0,0), Sigma=matrix(c(1,0.56,0.56,1), ncol=2),             empirical=TRUE), frequency=12)     head(N)          >            [,1]       [,2]     >[1,] -0.05270976  0.7239571     >[2,] -0.67232349 -0.6631604     >[3,] -0.20193415  0.8176053     >[4,] -0.54278281 -2.0458285     >[5,]  1.38279994  0.9405811     >[6,]  1.39979731  2.1717733          # Model: x(t) = a * x(t-1) + e(t)     fit = Arima(N[,1], order=c(1,0,0), include.constant=FALSE)          > fit       >Series: N[, 1]       >ARIMA(1,0,0) with zero mean               >     >Coefficients:       >         ar1       >       0.0293     >s.e.   0.1400       >     >sigma^2 estimated as 0.9791:  log likelihood=-70.42     >AIC=144.84   AICc=145.1   BIC=148.66          # Build the fitted values: x(t)=a * x(t-1)      pred  = fit$coef[1] * lag(fit$x, -1)      pred1 = fitted(fit)     head(cbind(pred, pred1))             >             pred         pred1     >[1,]           NA -2.255567e-05     >[2,] -0.001541849 -1.541849e-03     >[3,] -0.019666597 -1.966660e-02     >[4,] -0.005906915 -5.906915e-03     >[5,] -0.015877313 -1.587731e-02     >[6,]  0.040449232  4.044923e-02       In this case, `pred` and `pred1` match. However when I add in an `xreg`:               # Model: x(t) = a*x(t-1) + b*xreg + e(t)     fit1 = Arima(N[,1], order=c(1,0,0), xreg=N[,2], include.constant=FALSE)          >fit       >Series: N[, 1]       >ARIMA(1,0,0) with zero mean              >     >Coefficients:       >         ar1  N[, 5]       >       0.0860  0.5606       >s.e.   0.1401  0.1155       >     >sigma^2 estimated as 0.6675:  log likelihood=-60.85     >AIC=127.69   AICc=128.22   BIC=133.4          # Build the fitted values: x(t) = a*x(t-1) + b*xreg      pred2  = fit1$coef[1]*lag(fit1$x, -1) + fit1$coef[2]*fit1$xreg      pred21 = fitted(fit1)      head(cbind(pred2, pred21))          >              pred2     pred21     >[1,]         NA  0.4041670     >[2,]  0.4013329 -0.4112205     >[3,] -0.4296032  0.4325201     >[4,]  0.4410005 -1.2037229     >[5,] -1.1936161  0.5792684     >[6,]  0.6462336  1.2911169      In this case, `pred2` and `pred21` do not match, and the only thing changed was adding an `xreg`. The only time I cannot build out the fitted values by hand is when the AR part is included. I was able to do it when only MA parts were included with the `xreg`. I would really appreciate knowing how `Arima` treats `xreg` when generating the fitted values.


title: Misconceptions about random numbers in a range

text: I am currently working on a project the includes fitting. For the fitting I would like to try uniform random starting parameters. The possible fit parameters can lie in quite a large range, for example 10^-15 to 10^15 while others can lie in the range of -1 to 1. So quite diverse. When I try to create random numbers that reflect this range like               Number = Rnd(0, 1) * (UpperLimit - LowerLimit) + LowerLimit      This approach however is pretty useless for the large range. Here pretty much all resulting random numbers are very large. I have yet to see a number smaller than 10^12 for the above range. For the large range I can reduce the problem by using logarithms:               Number = 10^(Rnd(0, 1) * (Lg(UpperLimit) - Lg(LowerLimit)) + Lg(LowerLimit))      This now produces useful numbers over the whole logarithmic range yet it does not work for negative ranges or zero. I think I haven't really understood the basic concepts. If you plot the large numbers in a point plot they seem very uniform, yet they only reflect the large values. Does that mean that small numbers are less likely to turn up? The logarithmic numbers of course only look uniform in logarithmic scaling. I can't be the first person to run unto these problems. How are random numbers usually chosen? In the first step I don't care about quality of randomness. I just want to reflect larger and smaller number ranges including negative and 0 randomly. If I my explaination was to confusing please leave a comment and I will try to elaborate. As you may have noticed from the question I am note very well educated in statistics. So if you could keep it in Layman's terms I would be very grateful. Thanks in advance.


title: Three-parts correlation or equivalent?

text: I'm wondering what's the usual way for detecting "three-parts" association between variables. I'm not even sure about the name, so I apologize for any misunderstandings beforehand. Here is what I want. I suspect that buried behind some noise is a relation of this kind: $$s=\begin{cases} g(b) & \textrm{ if }f(a)\\ h(c) & \textrm{ otherwise } \end{cases}$$ where $a$, $b$ and $c$ are "input" signals, $c$ is an "output", $f$ has image $\{0,1\}$ and $g$ and $h$ are monotonic one-variable functions. In other words, $a$ acts as a switch so that in some cases $s$ and $b$ are correlated, and in others $s$ and $c$; but neither $g$, $h$ or $f$ are known beforehand.


title: Odds of drawing at least k identical values among m after n draws?

text: We draw $n$ values, each equiprobably among $m$ distinct values. What are the odds $p(n,m,k)$ that at least one of the values is drawn at least $k$ times? e.g. for $n=3000$, $m=300$, $k=20$. Note: I was passed a variant of this by a friend asking for "a statistical package usable for similar problems". My attempt: The number of times a particular value is reached follows a binomial law with $n$ events, probability $1/m$. This is enough to get odds $q$ that a particular value is reached at least $k$ times [Excel gives $q\approx 0.00340$ with `=1-BINOMDIST(20-1,3000,1/300,TRUE)`]. Given that $n\gg k$, we can ignore the fact that odds of a value being reached depends on the outcome for other values, and get an _approximation_ of $p$ as $1-(1-q)^m$ [Excel gives $p\approx 0.640$ with `=1-BINOMDIST(20-1,3000,1/300,TRUE)^300`]. _update: the exponent was wrong in the above, that's now fixed_ Is this correct? _(now solved, yes, but the approximation made leads to an error in the order of 1% with the example parameters)_ **What methods can work for arbitrary parameters $(n,m,k)$?** Is this function available in R or other package, or how could we construct it? _(now solved, both exactly for moderate parameters, and theoretically for huge parameters)_ I see how to do a simulation in C, what would be an example of a similar simulation in R? _(now solved, a corrected simulation in R and another in Python gives $p\approx 0.647$)_


title: What fields use modelling for analysis of a single time-series?

text: What fields of science (social or physical) use modelling for time-series analysis? In particular, time-series that cannot be replicated? Two examples that I can think of are climate modelling, and econometrics. Are there any others? I am asking because I am doing research on model independence in climate science, and am looking for references outside my field. I have had a look around, and couldn't find many, and was hoping to tap the knowledge resource here. I am reasonably sure that there is no "long list" to be requested, rather there are a few smallish research areas doing this kind of modelling. **Clarification** : By "cannot be replicated", I mean that only one time series can be relevant. eg. in climate science, we have only one earth, with which to validate the model. (in econometrics, there may be cases of both repeatable and not repeatable time-series)


title: LASSO/LARS vs general to specific (GETS) method

text: I have been wondering, why are LASSO and LARS model selection methods so popular even though they are basically just variations of step-wise forward selection (and thus suffer from path dependency)? Similarly, why are General to Specific (GETS) methods for model selection mostly ignored, even though they do better than LARS/LASSO because they don't suffer from the step-wise regression problem? (basic reference for GETS: http://www.federalreserve.gov/pubs/ifdp/2005/838/ifdp838.pdf - newst algorithm in this starts with a broad model and tree search that avoids path dependency, and has been shown to often do better than LASSO/LARS). It just seems strange, LARS/LASSO seem to get so much more exposure and citations than General to Specific (GETS), anyone have any thoughts? Not trying to start a heated debate, more looking for a rational explanation for why the literature does seem to focus on LASSO/LARS rather than GETS and few people actually point out shortcomings of LASSO/LARS.


title: Several bivariate correlation tests

text: I am looking at the relationship betwwen two interval level variables, being income and subjective well-being and I want to test the overall correlation as well as the correlation in three professional groups, namely doctors, high school teachers, and shop assistants. If I use bivariate correlation, I will have to run two tests, one for the overall relationship and the other for pofessional groups. I am using SPSS. Can this be done in one setup? Sorry, I am a newbie so regression may be an option but my focus is on correlation.


title: Is the following estimator for the normal distribution not asymptotically efficient?

text: Suppose that we have $X_1,...,X_n$ ~ i.i.d $N(\mu,\sigma^2)$, where $\mu$ is unknown and $\sigma^2 = 1. $ Is the following estimator: $\hat{\mu} = \frac{1}{n-1}\sum{x_i}$ consistent but not asymptotically efficient? There are different definitions of what it means to be asymptotically efficient and here, it means that if we have the form: $\sqrt{n}(\hat{\mu}-\mu)$ ~ $N(0,v(\mu))$, and $v(\mu) \geq \frac{1}{I_1(\mu)}$ where $I_1(\mu)$ is the Fisher information of one observation, then $\hat{\mu}$ is asymptotically efficient if $v(\hat \mu) = \frac{1}{I_1(\mu)}$. I tried to rearrange $\hat{\mu}$ but I cannot arrange it into the form above. Thank you!


title: Fitting a Poisson GLM in R - issues with rates vs. counts

text: I'm currently working on a project involving GLMs (and eventually GAMs) of some count data over time. Normally I'd do this in SAS, but I'm trying to move to R, and having...issues. When I fit a GLM to count data using the following: `cdi_model <\- glm(counts ~ exposure + covariate + month, data=test, family = poisson)` I get:               Deviance Residuals:          Min       1Q   Median       3Q      Max       -1.9825  -0.7903  -0.1187   0.5717   1.7649            Coefficients:                 Estimate Std. Error z value Pr(>|z|)         (Intercept)  1.97563    0.20117   9.821  < 2e-16 ***     exposure     0.94528    0.30808   3.068  0.00215 **      covariate   -0.01317    0.28044  -0.047  0.96254         months      -0.03203    0.01303  -2.458  0.01398 *       ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           (Dispersion parameter for poisson family taken to be 1)              Null deviance: 40.219  on 29  degrees of freedom     Residual deviance: 29.297  on 26  degrees of freedom     AIC: 137.7          Number of Fisher Scoring iterations: 5      Ignore for a moment the performance, or lack thereof of the model itself - mostly playing with syntax and the like at this point. However, when I try to fit rate data (counts/person-days) and use an offset like so: `cdi_model <\- glm(count_rate ~ exposure + covariate + months + offset(log(pd)), data=test, family = poisson)` I get 50+ warnings, all "1: In dpois(y, mu, log = TRUE) : non-integer x = 0.002082" etc. That is more than one for each observation (there's only 30 in the data set). Additionally, the model fit seems to go to pot. Output as follows:                Deviance Residuals:             Min          1Q      Median          3Q         Max       -0.0273656  -0.0122169   0.0002396   0.0072269   0.0258643            Coefficients:                  Estimate Std. Error z value Pr(>|z|)     (Intercept) -15.40110   15.12772  -1.018    0.309     exposure      0.84848   22.18012   0.038    0.969     covariate    -0.02751   21.31262  -0.001    0.999     months       -0.01889    0.95977  -0.020    0.984          (Dispersion parameter for poisson family taken to be 1)              Null deviance: 0.0068690  on 29  degrees of freedom     Residual deviance: 0.0054338  on 26  degrees of freedom     AIC: Inf          Number of Fisher Scoring iterations: 9      Despite this, if I plot the predicted rate against the actual data, the fit doesn't look that much worse, and the actual effect estimate doesn't seem to change all that much. Anyone have an idea what's going on - or if everything's going right and I'm missing something due to inexperience?


title: High-dimensional dependent binary explanatory variable

text: I am dealing with a data set containing roughly $n=4000$ **binary** observations $Y_1, \ldots, Y_n$ with $p=1000$ **binary** explanatory variables. I suspect that a lot of these explanatory variables are not relevant to the prediction of the observations. Moreover, it is clear that there are certain groups of explanatory variables that are very highly correlated. What are the usual approaches to deal with these kinds of situation. Indeed, quick simulations show that classical logistic regression and SVM do not work too well. Is there any method that can tackle the identification of these groups of highly correlated binary variables at the same time that the model is fitted to the data. I am looking something in the spirit of the LASSO that can do variable selection and model fitting at the same time.


title: How to uniformly project a hash to a fixed number of buckets

text: Hi Fellow Statisticians, I have a source generating hashes (e.g. computing a string with a timestamp and other information and hashing with md5) and I want to project it into a fixed number of buckets (say 100). sample hash: 0fb916f0b174c66fd35ef078d861a367 What I thought at first was to use only the first character of the hash to choose a bucket, but this leads to a wildly non-uniform projection (i.e. some letters apppear very rarely and other very frequently) Then, I tried to convert this hexa string into an integer using the sum of the char values, then take the modulo to choose a bucket:               import sys          for line in sys.stdin:         i = 0         for c in line:             i += ord(c)         print i%100      It seems to work in practice, but I don't know if there are any common sense or theoretical results that could explain why and to which extent this is true ? [Edit] After some thought I came to the following conclusion: In theory you can convert the hash into a (very big) integer by interpreting it as a number : i = h[0] + 16*h[1]+16*16*h[2] ... + 16^31*h[31] (each letter represents an hexadecimal number). Then you could modulo this big number to project it to the bucket space. [/Edit] Thanks !


title: questionnaire stacked bar plot with R

text: I have a file containing questionnaire data (categorical data, 4 possible values). I have imported this to R to get the corresponding data frame. Now I've done some messy things with the table function to obtain a new data frame which contains instead of the whole data listing only the frequency counts in right order of every variable, i.e. question in the questionnaire (maybe someone could tell me a good way to do this too). What I would like to obtain a graphical visualisation of the frequencies as a stacked bar plot. That's what it should look like: for each question, for stacked bars, the size of each bar corresponding to the frequency count (the frequency data in the transformed data frame that I constructed is already ordered). Is this possible? Thanks in advance!


title: Which distribution to use?

text: Here is a sequence of sample data where user is using some functionality:   * user_1 -> func_1   * user_2 -> func_2   * user_3 -> func_3   * user_1 -> func_1   * user_1 -> func_1   * user_1 -> func_1 Problem: How to calculate distribution of most popular functionality. So that I would get back result something like:   * func_1 -> 70%   * func_2 -> 20%   * func_3 -> 10% While at first glance it looks trivial problem to solve, the result can be very easily "hacked". Imagine we would want to track most popular video on youtube. If we just count how many times each video is viewed we can get wrong results, because you can have one crazy person who has watched the same video hundreds of times. This one person is kind of outlier. I am thinking towards some kind of weighted function, where each user importance is considered.


title: The distribution of Gaussian density at a fixed point given random mean

text: Let's say I have random variable $\mu \sim N(\mu ; 0, \xi^2)$. Is there anything useful known about the distribution of random variable $d = N(x ; \mu, \sigma^2)$, that is, the distribution of density at a fixed point $x$ given that the mean is a normally distributed random variable?


title: Metropolis sampling of constraint surface

text: I have a four variable $(x_1, x_2, v_1, v_2)$ state space for a system of ODEs. I would like to build a random sample of the initial conditions for these ODEs --- in other words, I wish to build a sample of vectors $\{(x_1, x_2, v_1, v_2),... \}$ from some PDF. However, I require that these vectors satisfy a constraint function $C(x_1, x_2, v_1, v_2)=E$, leaving a three dimensional surface to sample; and the PDF is defined only on this surface. This constraint surface is analogous to saying the vectors all have the same energy. It is therefore finite and bounds the origin. Suppose for example that I want a "uniform" sample of vectors satisfying this constraint (Note: by "uniform" I mean something akin to equally spaced along the surface, like electrons on the surface of a conductor). Without the constraint this would be a simple application of the Metropolis algorithm. With the constraint I have tried the following: I used a proposal distribution that is a Gaussian in the 4D space $(x_1, x_2, v_1, v_2)$ and took a random step in this direction (off the surface). I then attempted to project this new point onto the constraint surface by translating it back down along the normal $n$. This was, however, a nightmare computationally and only works for very well-behaved surfaces. Does anyone have a better approach?


title: Using AIC to distinguish between models using multiple datasets

text: I want to use AIC to compare three candidate models (labeled by m), each having K_m parameters. However, I have M datasets over which I can make the comparison. My ultimate goal is to report the "relative goodness" of each of the three models for a single fit. How to I make use of the multiple datasets? One idea is to find the Akaike weights (See Anderson and Burnham, 2002) for each dataset and average the weights over all the M datasets (perhaps, weighting by the number of points in each dataset?). Another approach would be to say that, for model m, I have a single M*K_m parameter model that I fit over the M datasets. In this case the AIC value, ACIC_net, is the sum of the AIC values for a fit to each of the M datasets. In this case, I would use AIC_net to compare the three models. How should I proceed?


title: How to get Standard Deviation from Mean Squared Error

text: So I've been working on the same question in my homework for days. I've done it over and over and I know I'm doing something wrong. I can't get a straight answer from anybody about it. I just want to understand. Right now, I'm pretty sure I've calculated my common standard deviation of the 5 varieties incorrectly. My margin of error is just much too high (which I'm guessing is due to an error in calculating my common standard deviation) compared to the point estimates that I found with my linear contrasts. The question that asks to calculate the information looks like so: What is the best estimate for the common standard deviation of the 5 varieties? (Hint: This value can be found (with minor calculation) from a value in the completed ANOVA table. You will need this value for parts (c), (d) and so on.) My math looks like so: MSE = 63.838 = s^2 s = √63.838 = 7.990 How do I calculate standard deviation from an ANOVA table? Is this not the way? EDIT * PARTS (c) and (d) c) Using the Bonferroni method for simultaneous inference with family-wise confidence of 95%, summarize all significant differences. Begin by identifying what would be the resulting individual confidence level for each interval. Note: You don’t have to construct all 10 confidence intervals. You can simply find E, the margin of error (it will be the same for all pairs since there are equal sample sizes), and then find all pairs of sample means that differ by at least E. Make sure to state clearly what E is and how it was calculated. Then, summarize the differences. d) Suppose varieties A and B were planted in Texas, varieties C and D in Florida, and variety E in California. Construct contrasts to make comparisons among the yield means that will address the following questions: i. Is the combined mean yield of the Texas varieties different from the combined mean yield of the Florida varieties? ii. Is the combined mean yield of Texas varieties different from the mean yield of variety grown in California? iii. Is the combined mean yield of Florida varieties different from the mean yield of variety grown in California? iv. Is the combined mean yield of the Texas and Florida varieties different from that of the variety grown in California? Use the Bonferroni method to construct simultaneous 96% confidence intervals for the above contrasts to assist in answering the following questions: Which of the contrasts are significantly different from zero? Which region has the most effective yield?