title: Time Series Forecasting Method to use both Predicted and Predictor variables

text: I am learning Predictive modeling and building a Forecasting model to predict Insurance sales in US as a part of my academic project. I want to do Time Series forecasting. I have Y(t) as my response variable and x(t),X(t-1),.....X(0) as my exploratory input varaibles that are correlated to the the response predicted variable. I have Y(t-1),Y(t-2),..... for close to 100 observations. I want to build a Forecasting model that uses both the Y(t),Y(t-1)..so on and their corresponding X(t),X(t-1) to build the forecasting model to predict the Y(t=1).....so on. When I went through some of the documentation available online, I saw that usually we generate a Timeseries ID and give the resposne variable and the time series ID( which is a continuous variable starting from 1 to count of response variables) and this will not use any of the X(t) values. Is there any way where I can use both X(t),Y(t) to predict Y(t+1). Thanks in Advance, Sai


title: Sample/Null Value

text: I was running through question section on my stat book and the following question was posed: > If our sample value is close to the null/population value, we conclude that > ... I don't get if that is suppose to reject/accept or if there is an effect. Furthermore, can anyone explain what is meant by the sample value?


title: Skewness adjusted t-stat

text: I've some data that is divided into a series of groups and am testing whether the mean is different from 1. The data is highly skewed, with skewness ranging from 1.10 to 26! I did a one-sample t-test to see were the means different from one, in spite of the skewness. I also did a one-sample median test (Wilcoxon Signed Rank) and have produced p-values for both. Am I at risk of a TYPEII error when looking at the one sample t-test as the data is highly skewed? Would there be any benefit to using a skewness-adjusted t-statistic (as in Johnson (1978) bagusco.staff.ipb.ac.id/files/2014/01/modified-t-test.pdf)? If I'm showing a table of this, should I produce both the one sample mean and the one sample median statistics? Any other suggestions? EDIT: A few clarifications based on the answer below. The underlying data is a series of sequential bets. Each bet has a certain stake size and is a win or a loss. Some of the stake sizes are very small, some are quite large, so I standardised the first bet to 1. I'm interested in the path-dependent changes in stake size, so the ratio is ok, I think. The groups (which are nodes in the betting/result sequence) go from 'L' and 'W' for the first bet, to 'WWWW', 'WLLW' etc for the fourth bet. If stake sizes are independent of the result/path, the stake size ratio should still be 1 after 4 bets. If it's not, there's path-dependency. I want to test whether the stake sizes in each group are changing and whether they are statistically significantly different from 1. >   * Are you testing each group individually for mean=1? (my subsequent > discussion assumes so) > > > Yes. There's 16 groups. If the mean standardised stake is still 1, the stake > size changes are independent of the previous results. I expect there to be a > difference from one for those that had 3/4 winning bets and those that had > 3/4 losing bets. > >   * mean 1 sounds a little unusual. Are we dealing with ratios in any way? > > > Yes. As explained above. > >   * What sample sizes are these? > > > There's at least a few hundred observations in most groups, while some only > have 60/70 observations. > >   * Are the data discrete or continuous? > > > The data are continuous, but as mentioned above, it's stake/bet sizes so > highly skewed. Some agents bet small, some bet very big. > >   * Is it definitely the population mean you're interested in? > > > I'm not too sure. Maybe not. Could I use the Sign Test since it doesn't assume symmetry? If the median in each group is changing, that might answer the question of whether there's anything happening along the result path of these groups.


title: Multinomial logistic regression

text: I would like to get help from an expert or anyone who knows about this. As a beginner in SPSS, I've googled the steps to do multinomial regression. I've decide to use multinomial because from what I googled and understood, binomial is for situations that only have 2 categories like yes or no, which you could code 0 for yes and 1 for no. This is what I understand if you choose to use binomial. My situation is I have a set of questionnaires that use more than 2 answers. For example, `level of income: 1. $2000, 2. $3000, 3. $4000, 4. $5000`. That is the kind of question I have. So, obviously I choose to use multinomial regression. My IV's are 1) social capital, 2) training, 3) credit/loan; and my DV is effectiveness which can be measured by 1) income, 2) saving, 3) repayment rate. In my questionnaire, I only have two kind of question which "yes or no" and choose answer like above mention "1 or 2 or 3 or 4." So my questions are:   1. Did I choose the right regression to use, that is multinomial, or do I have to use binomial?    2. If it's multinomial, I need help with the steps like do I have to calculate all my DV questions and transform them into one, which will make new "item" in SPSS columns, then put it in the dependent column of multinomial regression? Because I have approximately 15 questions to measure my DV. Then in the factors columns what did I put? Every question of my IV's (like 30 questions separately) or just the total of the 30 questions? And what should I put on covariates?  http://sdrv.ms/JtOHu3 above is my questionaire sample link. income level, savings and repayment is to measure my dv ~ effectiveness of microfinance. other is my iv, that is credit/capital , training and social capital. feel free to improve my questionaire if you like.


title: Is there an R package for a multinomial repeated measures logistic regression with 3 levels of the dependent variable?

text: dependent variable: one of three forms of a word-ending independent variables: age, language, word stem with 3 levels random effects: speaker, sentence Each speaker produces many sentences, so it repeated measures is required.


title: Can we apply multiple regression on time series data?

text: It always create a doubt to me, whether we can apply linear or non linear multiple regression on time series data. If yes, should I consider year also an independent variable. Thanks,   Arushi


title: Using averages of two different size groups

text: I'm not really "into" statistics, so I could use a hand on this one. I have two small groups, one has 12 students and the another one has 6. Both groups answered a questionnaire and I'm trying to find out which of the groups did better. I'm not sure if I should be using the average of correct answers to do this comparison. Can't this be unfair with the small group? Should I be using just 6 students of the larger group?


title: estimating standard error for small sample sizes

text: I have generated a number of objects in 3D space (2D rectangles). I have sampled those objects along lines, calculating the number of intersections between the lines and objects. If I use increasing amounts of lines, the standard error in the number of intersections decreases, as expected. What I really want to know is: How many sampling lines do I need to use to get a good idea of the # of objects (or intersections). I have plotted a line graph of the mean number of intersections per line vs the number of sampling lines, along with the +/- 95% confidence interval for the mean number of intersections per line. At high numbers of lines, the graph looks good and stable. But, at low numbers of intersection lines, the graph is erratic and somewhat unstable due to the smaller sample size. How could I do a better job of calculating the mean and confidence interval for small sample sizes? Or maybe, can/should I analytically derive the number of intersections and confidence interval? Not sure how I would do that.![My Chart](http://i.stack.imgur.com/lca2h.png)


title: What is the difference between "priors" and "likelihood"?

text: I have been reading some papers regarding statistics and I seem to be confusing the terms priors and likelihood. Would it be possible to explain the difference between the two terms? I am interested in both a "down-to-earth" approach with examples and the mathematical and statistical aspects. Thanks.


title: How best to normalize count data to compare two distributions

text: Say I have a vector of length 1000. At each position (1 ... 1000) there is a count. I have two vectors with different range of counts such that in vector A the maximum number of counts at a position is 30, whereas in vector B the maximum is say 200 (i.e., there are more counts in B than A). So essentially I have two discrete distributions and when I plot the two distributions they have peaks and troughs (some regions along the vector have higher counts than others). I am not interested in the differences in raw counts between A and B but instead wish to compare the shape of the curves in order to test whether the same regions (based on the index) in A and B have the higher counts relative to other regions. My problem is that in order to compare the shapes I need to normalize the counts in each vector. I do not have have too much stats background and so any advice regarding an appropriate normalization procedure would be appreciated. The best, but probably not the best solution I can think of is to transform the vectors to have the same mean.


title: Choosing the threshold for computing the confusion matrix

text: Is there any way of selecting the best possible threshold to minimize the number of false positives and false negatives (to get the best possible minimized combination of false +ves and -ves)?


title: Calculating confidence intervals

text: Body mass index was compared for two groups, people with elevated triglycerides (above 1.7 mmol/L) and people with normal triglyceride levels. The 10 people in the group with elevated triglyceride had body mass index mean = 26.1 and standard deviation = 3.72 and the 15 people in the normal group had body mass index mean = 24.3 and S = 3.45. What are the confidence intervals for these means and without doing any further calculations, decide whether a test of the null hypothesis of no difference against a two-sided alternative at 5% would lead to rejection or not


title: Bootstrap or jack-knife for crossvalidation of predictive model?

text: Is a bootstrap or jack-knife method better for crossvalidation of a multivariate logistic regression based predictive model?


title: Small conditional probabilies textbook question

text: I have an example from a book I am working through (car starting problem, with fuel & dirty spark plugs), but need a little help please. We have probabilities: P(NS) = 0.016, P(NF) = 0.001, P(F) = 0.999, P(SPD) = 0.05 and P(SPC) = 0.95. Where S = car starts, NS = car doesn't start, F = there is fuel in car, NF = there is no fuel in car, SPD = spark plugs are dirty and SPC = spark plugs are clean. And probabilities:               ........   NS   S     NF   SPD   1    0     NF   SPC   1    0     F    SPD   0.1  0.9     F    SPC   0.01 0.99      (So for example, the first row: if the car has no fuel (NF) and spark plugs are dirtty (SPD) the car will not start with probability one) I want to calculate P(SPD|NS). My attempt:               P(SPD|NS) = P(SPD and NS) / P(NS) = P(NS|SPD)P(SPD) / P(NS)      So I need P(NS|SPD) – need for all fuel status (F or NF). I think this should be = 1*0.001 + 0.1*0.999. I am going in circles trying to write the above line as probabilities – what is the proper way to present this in terms of NS, SPD, F and NF (assuming it is right). Thanks.


title: Best distribution for insurance restitution

text: I wish to discover the best distribution to model **non-zero** restitution* values of insurance as a response variable in a regression. The problems I find here:   * Knowing that the data came from insurance restitution doesn't seem bring me any information about its distribution.   * Any restitution value is valid in theory, but in fact, all data have discrete values. So should I model it as discrete or as continuous? ![Histogram](http://i.stack.imgur.com/niE7I.jpg) So far, I thought of using zero-truncated negative binomial, gamma and log- normal. (*) Is restitution the word English speakers use for this?


title: Supervised learning approaches which can accommodate a supervisory signal composed of multiple dependent continuous variables?

text: This is a supervised learning problem. Ideally would like to work in R due to having an easy way to pre-process the input data, but could work around that as well. For each sample, input consists of tens of thousands of features. These are genomics data and will likely need to be reduced to a manageable amount, somehow, before being used to train the classifier. Supervisory signal consists of 4 dependent continuous values, representing relative composition of the sample. e.g. continuous between 0 and 1, all 4 summing to 1 for each sample:               Sub012  0.5940594   0.26732673  0.07920792  0.059405941     Sub013  0.5102041   0.34693878  0.08163265  0.061224490     Sub014  0.6521739   0.20652174  0.07608696  0.065217391      Wanted: a regression function capable of predicting the relative composition of a sample in terms of those same 4 dependent continuous values. The constraints on the supervisory signal are what is causing me pause: the dependence of the variables, being constrained between 0-1 and summing to 1. I was hoping someone might have attempted something similar and could point me in the right direction - packages or approaches which may work or definitely won't work - all thoughts welcomed. Thank you.


I have a large text and I need to select words randomly in proportion to their frequencies. A raw method is to put all the words in an array (N indexes) and to use a random (uniform) generator to generate numbers between 0 to N. For example, a word $w_1$ occurs 5000 times in the text and the word $w_3$ occurs 100 times. $w_1$ should be selected with a probability that is 50 times larger than the probability for $w_3$. This "raw" method does not seem to be effective (or feasible) if a I have a large text. Is there a better way to select the items randomly and in proportion to their frequencies?


title: Probit - multivariate

text: I have a data set of ~260 persons who answered 15 relevant items and were then passed on to get a specific treatment or not... I am trying to predict, whether they will recieve treatment or not. I use SPSS with all available add-ons. As far as I understand, I have to go for a Probit because the dependent is binomial. The measuring scale for my items is mixed; many are ordinal ("not at all", "little", "neutral", "quite", "positive"; which I coded 1-5); the rest of the variables are binominal except for two which are metric. However, in SPSS (and the "help"-program that comes with it) I only find a Probit asking for "response Frequency", "Total observed" and "Factor", primarily all referring to one item. Then, I can also put covariates. I Understand that I could create new variables counting the observed Frequencies and total cases observed. However, it would not make sense to do this with every variable!?! Is there a method to put many different variables without converting them - comparable to "optimal scaling" with variables of different scaling?


title: Expected value of the natural log of a ratio of variances

text: I am dealing with a one-way random effects model and am looking for the $E(\ln(\hat{\sigma}_\alpha^2/\hat{\sigma}^2))$ where $\hat{\sigma}_\alpha^2$ is the estimate of the between group variance and $\hat{\sigma}^2$ is the estimate of the error variance. I am basically just interested in what a good estimate of the bias would be. I know from Jensen's inequality that it is negatively biased.


title: What are the mean and variance for the Gamma distribution?

text: There are two forms for the Gamma distribution, each with different definitions for the shape and scale parameters. Rather than asking what the form is used for the gsl_ran_gamma implementation, it's probably easier to ask for the associated definitions for the mean and standard deviation in terms of the shape and scale parameters. Any pointers to definitions would be appreciated.


title: Sample size to achieve given confidence level

text: I am trying to estimate number of different unique visitors who visited a given website (online store). There are hundreds of millions of visits to the store and so this task is too difficult to handle by my database. Note that one visitor may have more than one visit so it is not enough to count number of all records, we need to find number of all unique visitors. Even when scaling the data into MapReduce cluster this task turns out to be quite difficult. I decided to sample the data. Each visitor has a Md5 hash (which is a uniform random hexadecimal number) and so I can select all users whose Md5 hash string finishes with for example '0'. By doing so I am sampling 1 in 16 users (1 hex number = 16 possibilities). In this way I am limiting number of records and the database now can select number of unique visitors (number of distinct visitor hashes). I can simply multiply this number by 16 to get the estimated total number of unique visitors. I did some tests on a small sample of a real data and it works pretty well. I am trying to find out what is the minimum sample size to get 95% confidence level that the number I arrived at is correct. Some pages are visited by millions of visitors while other have only few hundreds visits. If I will sample 1/16 users on a page which has only 20 users my estimate may be quite off. Sample data looks in the following way:               Visitor Hash | Website ID | Timestemp | Irrelevant data     009AB730 | 123 | 11111 |      009AB730 | 122 | 11112 |      009AB734 | 122 | 11112 |      0283AB22 | 122 | 11112 |      0283AB22 | 122 | 11112 |      .. repeated 1000 times     0283AB22 | 122 | 11112 |      0283AB20 | 122 | 11112 |       After sampling this data based on user hash ending with "0" we get:                009AB730 | 123 | 11111 |       009AB730 | 122 | 11112 |       0283AB20 | 122 | 11112 |       So it is irrelevant how many times a given user visited the website, what matters is whether his Hash/Id finishes with 0. This hash should be uniformly random and so we can assume that 1/16 of users IDs finishes with 0. Moreover it is easier to calculate number of unique IDs in the following smaller sample. That is why I was wondering what should be my minimum sample to achieve less than 5% error? I found this question but I am not sure if I can use the same calculations? There is no proportion in my problem. My first take on this problem would be to use the same technique and use the worst case scenario assumption by taking largest margin of error when the proportion is 0.5. Using exactly the same logic: the maximum margin of error at 95% confidence is $m = 0.98 / \sqrt{n}.$ Rearranging gives $n = (0.98/m)^2.$ And so if we want a margin of error of 5% = 0.05, $n = (0.98/0.05)^2 = 384$. Therefore we need 384 unique visitors for each webpage in which we want to estimate the number of unique visitors. Since this is the worst case scenario I hope this logic is not flawed? I am worried that I did not use the sampling 1/16 in this calculation. What if I would look for users whose ID finishes with '00' - I would expect that the minimum sample should be bigger in that case.


title: Sample size for actual vs backtested performance

text: I have 10 years of backtested simulated performance of some trading strategy (using historical prices), and N months of actual trading performance. I want to compare the two. How big does N have to be to get a statistically significant comparison? I'm comparing returns and Sharpe ratios.


title: Choosing the number of latin hypercube samples

text: I'm using latin hypercube sampling to draw parameter values for an epidemiological simulation model (about 30 parameters, following various distributions), and so far I haven't been able to find any good rule of thumb/heuristics/theory for choosing the number of samples (i.e. the number of equiprobable intervals).


title: Comparing means of three groups

text: When I want to compare the means between two groups I use a t-test. But how do I check whether the means of three groups are equal or statistically different? E.g. compare the average labour productivity of small, medium and large firms.


title: Egger's test in SPSS

text: How can I conduct an Egger's test using SPSS17? For each study included in the meta-analysis I know effect size and sample size of patients and controls groups.


title: Is winning a soccer match independent of previous wins\losses?

text: $\quad$ I have a friend of mine who is a bit of a gambler ask me this question. He is of poor mathematical background, but has a sense of logic and will probably accept a logical answer in the natural language. $\quad$ So my friend asked me if the probability of an outcome would change in consecutive trials of playing roulette. I said that it will not, and he needed no further verification. He only asked if it is because the roulette has no memory. I said he could put it that way if he liked. $\quad$ Then the person asked me if the same is true for soccer matches. I told him that, despite the influence on soccer by other factors, the same should hold. My friend insisted on the contrary, and I didn't know what to say to him. $\quad$ I would like to be able to explain to him why soccer matches are independent events. If this is not true, what is the dependence? Would considering a particular team change anything. Furthermore, what restrictions can we apply to ensure independence? Having the same opponent team? p.s. I wouldn't mind an answer including mathematics in addition to any other answers. p.p.s. I asked the exact same question on math.SE, where I am told it does not belong. I hope this is the appropriate place to ask.


title: Calculating hourly volatility and peak-to-average ratio in R

text: I have network traffic data in the following for for each hour of a ten day period as follows in a R dataset.                  Day   Hour         Volume          Category         0    00            100            P2P         0    00            50             email         0    00            200            gaming         0    00            200            video         0    00            150            web         0    00            120            P2P         0    00            180            web         0    00            80             email         ....         0    01            150            P2P         0    01            200            P2P         0    01             50            Web         ...         ...         10   23            100            web         10   23            200            email         10   23            300            gaming         10   23            300            gaming      As seen there are repetitions of Category within a single hour also. I need to calculate the volatility and the peak hour to average hour ratios of these different application categories. **Volatility** : Standard deviation of hourly volumes divided by hourly average. **Peak hour to avg. hour ratio** : Ratio of volume of the maximum hour to the volume of the average hour for that application. So how do I aggregate and calculate these two statistics for each category? I am new to R and don't have much knowledge of how to aggregate and get the averages as mentioned. So, the final result would look something like this               Category    Volatility      Peak to Avg. Ratio     Web            0.55            1.5     P2P            0.30            2.1     email          0.6             1.7     gaming         0.4             2.9


title: STFT statistical analysis

text: I am using the `evolfft` function in the `RSEIS` R package to do a STFT analysis of a signal. The signal is one hour long and it was acquired during 3 different conditions, in particular 0-20' control, 20'-40' stimulus, 40'-60' after stimulus. Visually, I see a change in the spectrogram during these 3 periods, with higher frequency and increased FFT power during the treatment, but I was wondering whether there was some kind of statistical analysis I could do to "put some numbers" on it. Any suggestion? EDIT: as suggested I will add an example of the data I'm dealing with ![STFT example](http://i.stack.imgur.com/5wPy0.png) The treatment is between 20' and 40', as you can see it produces an increase in the power of the FFT over a fairly wide range of frequencies. I have 50-60 of these STFT for each experiment (for 10 total experiments). I can average the spectra for each experiments and still get a similar type of pattern. Now, my problem is how to correctly quantify the data I have and possibly do some stats to compare before, during and after the treatment.


title: Negative binomial distribution fit

text: I am trying to fit a negative binomial distribution, in `R`, to my over dispersed data (out of 20 ,14 samples are 0, and rest are less than 5). The mean is $-0.8$ and the variance is $2.69$. The problem is if I use the `theta.md` function, I get theta around $0.5$, and if I use `glm.nb`, then I get theta as $.25$. Any idea why the difference is so large? My data are: 4 1 1 6 3 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 Thanks for the help.


title: Multiple regression with repeatedly measured independent variables?

text: **Design and hypothesis** : we measured `wellbeing` at Time-1 and Time-2, we want to see whether factor `A` (measured at Time-1 and supposed to be a stable factor over time) is a significant predictor of factor `B` (measured at Time-2). We also expect `wellbeing`, current or past, will contribute to `B`. **Question** : is it appropriate to do multiple regression with `wellbeing` measured at both time points (using the same instrument) as simultaneous predictor? - significant correlations amongst predictors are present, but multicollinearity diagnostics seemed fine... is there a better way to test the hypothesis that'd make good use of the _longitudinal_ design? **Many thanks!**


title: Time series and multiple variables

text: I created a time series in Excel (not ideal) using Holts-Winters to forecast daily loan values in a month and it works very well. I've been asked to build a similar model that integrates other variables but other than a multiple regression model (which has it's own issues) I'm not sure with my limited experience how I could achieve this and within Excel. Any suggestions would be most welcome Regards


title: Modeling a time series using various approaches (including auto.arima() in package forecast)

text: I have tried to model a time series (https://dl.dropboxusercontent.com/u/13029929/ts.txt) with various approaches, but haven't been satisfied. The original time series and its ACF and PACF look like ![enter image description here](http://i.stack.imgur.com/sZbkV.png)   1. It doesn't look weakly stationary to me, because the mean and variance don't seem constant. I want to test its stationarity by the augmented Dickey-Fuller test, but the test gives a very small p-value, rejecting the existence of a unit root. (I wonder if I misuse the adfTest() function?) > adfTest(tmp,lags=20,type=c("c"))              Title:  Augmented Dickey-Fuller Test          Test Results:   PARAMETER:     Lag Order: 20   STATISTIC:     Dickey-Fuller: -12.4317   P VALUE:     0.01           Description:  Thu May  1 00:47:45 2014 by user:           Warning message: In adfTest(tmp, lags = 20, type = c("c")) :      p-value smaller than printed p-value      Its ACF and PACF seem to suggest a AR(20) model, so I fit it with AR(20), but the residual series didn't pass Ljung-Box test (p-value = 1.11e-15):              out = arima(tmp, order=c(20,0,0))          Box.test(out$residuals, lag = max(log(length(tmp)), 20)+10, fitdf = 20, type = "Ljung-Box")        2. Because the series doesn't look stationary to me, I then difference the series once, and the differenced series seems to have constant mean, although still nonconstant variance. ![enter image description here](http://i.stack.imgur.com/zf2Ss.png) It looks like a MA(6) process, but after fitting MA(6), the p-value of Ljung- Box test on the residual series is still too small (p-value < 2.2e-16)   3. I then think of taking logarithm of the original series, to reduce variability: ![enter image description here](http://i.stack.imgur.com/3Wm8x.png) and then difference the logarithm of the original series (see below): ![enter image description here](http://i.stack.imgur.com/CSXJb.png) I am not sure if it is worth to take the logarithm, and what models to try with the differenced logarithm series, as I tried several ARMA models on it, the residual series doesn't pass the Ljung-Box test.   4. I tried `auto.arima()` in package `forecast` written by Prof Hyndman, on my time series. It suggests ARIMA(3,1,3), while, suprisingly, my `adfTest` doesn't say the series has unit root (see part 1). But the residual series still has spikes, and doesn't pass Ljung-Box test > out1 = auto.arima(tmp) > > out1              Series: tmp      ARIMA(3,1,1)                              Coefficients:          ar1      ar2     ar3      ma1       0.9747  -0.2607  0.1460  -0.9811     s.e.  0.0063   0.0081  0.0061   0.0020          sigma^2 estimated as 0.1172:  log likelihood=-9881.62          AIC=19773.24   AICc=19773.25   BIC=19814.53      > Box.test(out1\$residuals, lag = max(log(length(out1\$residuals)), 4)+10, > fitdf = 4, type = "Ljung-Box")              Box-Ljung test          data:  out1$residuals     X-squared = 61.8217, df = 16.257, p-value = 3.101e-07      The residual series and its ACF look like: ![enter image description here](http://i.stack.imgur.com/S8BfW.png)   5. It is similar to apply `auto.arima()` to logarithm of the original time series, which suggests ARIMA(3,1,3) I really appreciate if you could suggest some models (and transforms on the original series) to try. Thanks!


title: Aggregated dispersion measure

text: I am using the function `betadisper {vegan}` to test for differences in, in my case, raw genetic dispersion between sites within species. I would like to derive a _single measure of dispersal_ for each species, thus pooling over sites (factor levels). This _dispersion measure_ would be the `betadisper's` _Average distance to median_ for all sites within a species. To achieve this, I have tested to run `betadisper` with a distance matrix and a grouping factor with a single level. When plotting this, Im seeing something I can't understand; the centroid and where the spiders converge is not the same point (see 2nd plot B). However, doing the exactly same thing, but with the `varespec data {vegan}` (see below), i.e. changing the number of levels from two ("grazed", "ungrazed") to only one (lets say "A") the plot looks as you would expect (see 1st plot). The only difference here is the distances used; genetic p-distance (raw diff between basepairs) and bray-curtis. Any pointers on why my 'one-factor (1B) plot' looks strange in regards to spider convergence and centroid position compared to plot 1A, or if it's even legit to derive a single measure of dispersion as I do here, would very much appreciated, thanks!               #Dummy code          require vegan     data(varespec)          ## Bray-Curtis distances between samples     dis <- vegdist(varespec)          ## First 16 sites grazed, remaining 8 sites ungrazed     groups <- factor(c(rep(1,16), rep(2,8)), labels = c("grazed","ungrazed"))          ## Pool over sites, i.e. factor with a single level, e.g. A     group <- 1:length(groups)     group[1:24] <- 'A'     group <- as.factor(group)          ## Calculate multivariate dispersions     mod <- betadisper(dis, groups) # plot 3A          ## Calculate multivariate dispersions using pooled data     mod2 <- betadisper(dis, group) # plot 3B          ## Average distance to median for two-level factor:     tapply(mod$distances, mod$group, mean)          ## grazed  ungrazed      ## 0.3925879 0.2706457          ## Deriving a single measure of dispersion, i.e. Average distance to median for one-level factor     tapply(mod2$distances, mod2$group, mean)          ## A      ## 0.4254586       This is vegans varespec data; plots as you would expect. **A** is one-level factor and **B** is two-level factor ![vegans varespec data](http://i.stack.imgur.com/dOUP5.jpg) This is for one of my species. I do the same procedure as above. **A** is two-level factor and **B** is one-level factor ![enter image description here](http://i.stack.imgur.com/rD2El.jpg)


title: R - how to let glmnet choose lambda range when using caret?

text: To fit a lasso model using `glmnet`, you can simply do the following and `glmnet` will automatically calculate a reasonable range of `lambda` values appropriate for the data set:               glmnet(x, y, alpha = 1)      I know I can also do cross validation natively using `glmnet`. However, I would like to use the `caret` package so I can train and compare multiple models in a unified fashion. The problem is I don't know how to fix only the `alpha` parameter when calling `glmnet` from `caret` because `caret` wants to tune over both `alpha` and `lambda`.               train(y ~ ., data = train, method = 'glmnet', trControl = ctrl, tuneGrid = data.frame(alpha = 1))     Error in train.default(x, y, weights = w, ...) :        The tuning parameter grid should have columns alpha, lambda      Is there any way in general to specify only one parameter and allow the underlying algorithms to take care of the unspecified parameters according to the default methods of the algorithms? EDIT: The reason I would like to let `glmnet` to choose the `lambda` is that the appropriate range for `lambda` can vary a lot and `glmnet` does a good job of choosing a good range:               library(ElemStatLearn)     library(glmnet)          dat <- prostate          train   <- subset(dat,  train, select = -train)     train.x <- as.matrix(subset(train, select = -lpsa))     train.y <- train$lpsa          print(glmnet(train.x, train.y, alpha = 0)$lambda)     print(glmnet(train.x, train.y, alpha = 1)$lambda)      You can see that the ranges of `lambda` are very different. To use custom model in `train`, I thought about running `glmnet` first and then assign the calculated `lambda` to the grid:               glmnetGrid <- function(x, y, len = NULL) {       library(glmnet)       lam <- glmnet(x, y, alpha)$lambda       expand.grid(lambda = lam)     }      But then I need to pass `alpha` to the grid function and I don't think the framework allows that. Also it seems inefficient to run `glmnet` one more time just to get the `lambda`. Is there a way to simply let `glmnet` run and assign its `lambda` to the grid so `caret` can use the `lambda` values in subsequent analysis?


title: How are user opinions compared?

text: Consider a situation where people are asked to rate a product on a scale of 1 to 5. Product A has a rating of 3.95 based on 15,000 ratings. Product B has a rating of 4.2 based on 1000 ratings. From the data given, can we determine which product is the preferred one? Is there a direct formula to solve problems of this type? Can anyone explain in layman's terms? PS: Please retag as appropriate. I couldn't find a relevant tag for this question.


title: The relationship between trading profitibility and forecasting accuracy

text: is forecasting accuracy and trading profitability related? How to explain the existence of the relationship? Thank you.


title: What is the distribution of the difference of two-t-distributions

text: ... and why ? Assuming $X_1$,$X_2$ are independent random-variables with mean $\mu_1,\mu_2$ and variance $\sigma^2_1,\sigma^2_2$ respectively. My basic statistics book tells me that the distribution of the $X_1-X_2$ has the following properties:   * $E(X_1-X_2)=\mu_1-\mu_2$   * $Var(X_1-X_2)=\sigma^2_1 +\sigma^2_2$ Now let's say $X_1$, $X_2$ are t-distributions with $n_1-1$, $n_2-2$ degrees of freedom. What is the distribution of $X_1-X_2$ ? **This question has been edited:** The original question was _"What are the degrees of freedom of the difference of two t-distributions ?"_. mpiktas has already pointed out that this makes no sense since $X_1-X_2$ is not t-distributed, no matter how approximately normal $X_1,X_2$ (i.e. high df) may be.


title: Number of times to run a lengthy experiment

text: I have this long running experiment. Each time I run it I get a new goodness value, since the algorithm has random variables in it. So I need to report the mean and the std of some n runs. What should n be? I need to be able to defend n based on some statistical ideas. Some kind of scientific reference (a book, a paper) would be wonderful, too. I provide more details as you say, thanks for the answers: In computer vision, an important challenge is to recognize objects from images. Different computer algorithms are developed for this purpose. To see how good a new algorithm is, one sometimes constructs a test and a training set of images, say 1000 images for each, train the algorithm with the training images, and produce a success rate using the test set. If out of the 1000 objects in the test images, 800 is recognized by the algorithm, the success rate is said to be 80 percent. Now, my algorithm analyses, say, 1000 RANDOM points in the image, and using that analysis, tries to recognize the objects in the image. Each time I run the algorithm, I get a different success rate, since 1000 points are produced RANDOMLY. So I think its best to report some kind of summary statistics (e.g. the main and std deviation) of the success rate. Also, one sometimes needs to say, "well in addition to my algorithm, I tried these, say, 10 algorithms on the same dataset, and this table shows that mine is the best in this and this way..." Some of these algorithms may need to run more than once, too. So one can really have a long experiment. So, as I said before, at least how many times should I run the long running experiment? Thx.


title: Fitting normal distribution to given quantiles

text: Hey everyone, thanks for taking the time to look at this question. This is pretty simple, and I understand statistics quite well, but I think I'm not wrapping my head around the words, it seems to me as if there's something missing. Here's the question. > The manager of the aerospace division of General Aeronautics has estimated > the price it can charge for providing satellite launching services to > commercial firms. Her most optimistic estimate (a price not expected to be > exceeded more than 10 percent of the time) is 2 million. Her most > pessimistic estimate (a lower price than this one is not expected more than > 10 percent of the time) is 1 million. The expected value estimate is 1.5 > million. The price distribution is believed to be approximately normal. > > What is the expected price?   >  What is the std dev of the launch price?   >  What is the probability of receiving a price less than 1.2 million? So what is the actual probability of a 2 million dollar price? It just says that they don't expect to exceed 2 million more than 10% of the time. Does that mean 90% of the time they expect the price to be less than 2 million? Thanks for your suggestions.


title: Half annual covariance $\mathrm{Cov}(X,Y)_{\text{6 months}}$ from annual covariance?

text: I have annual returns and standard deviations for two funds, $r_{a}$, $r_{b}$, $\mathrm{SD}_{a}$ and $\mathrm{SD}_{b}$ but I do not have individual data, just the annual data. The annual correlation between the prices of the funds is 0.7. If I had the individual data I could use $\mathrm{Cov}(A,B) =\sum_i (\bar{r}_{a} - r_{i a}) (\bar{r}_{b} - r_{i b}) \>,$ but now I am lost without the data. Are there some approximations to get the half annual data from the annual data or some formula to do it directly?


title: Algorithm to create Bayesian priors from measurements

text: I am trying to devise a simple algorithm to create a Bayesian prior from measurements obtained from time series data. Firstly, I presuppose that the data can take on one of five possible "shapes" or "patterns" and then measure the Euclidean distance of each time series data point from its equivalent point in time in each of the five possible patterns, summed over a specified period of time for each pattern. This will give five distinct numbers with the lowest one indicating which pattern the data is most likely to resemble due to the fact that this distance metric is the smallest "error measurement." However, this is not a sufficient measure in and of itself hence my desire to use it as a prior in a naive Bayesian classifier. My problem is how to use these numbers to form a prior? I am thinking of using some form of ratio i.e. the numbers are 1 2 3 4 5 for a total of 15, therefore 15/1=15, 15/2=7.5, 15/3=5, 15/4=3.75, 15/5=3 and so the prior for 1 would be 15/(15+7.5+5+3.75+3), prior for 2 would be 7.5/(15+7.5+5+3.75+3) etc. Is there any justification for using an approach such as this, or is there perhaps a better approach?


title: Comparing numbers of p-values from many linear models

text: My current dataset has three conditions, and we've measured the activity levels of 10,000 genes in each condition. Replicated 8 times. Using 10,000 linear models, we determine for each pair of conditions (ie for each of three contrasts) the number of genes with significantly different activity levels. This is standard procedure for this kind of microarray data. We find:   * 2000 genes have significantly different activity levels between A and B   * 1500 genes have significantly different activity levels between A and C   * 100 genes have significantly different activity levels between B and C This suggests that conditions B and C are more similar to each other than to C. PCA suggests the same result. Is there any way for us to quantify the extent to which "conditions B and C are more similar to each other than to C (ie to put a p-value on it?) Thanks for your help, and apologies if this question is trivial. Kind regards,


title: Fill colors in Gnuplot

text: I am doing some histogram rowstacked plots using the gnuplot command               plot 'file.dat' using (100.*$2/$6) title column(2), for [i=3:5] '' using (100.*column(i)/column(6)) title column(i)      I want to change all the colors used for fill the bars in the histogram. For change the first box color, I simply add _lt rgb "color"_ before the comma, but I don't know how to modify colors in the loop to make each box have a different color. Modifying the line's style does nothing.


title: Computing average value ignoring outliers

text: This is more of a general statistics question, though if it matters I'm writing PHP code. Let's say I'm trying to compute the average value of a toy that is commonly bought and sold on the secondary market, and I have a set of price values culled both from auctions and from user-entered "price paid" data. The data points that represent auctions are pretty reliable, but I also get the occasional "garage sale" type of data point, where someone may have paid a buck to buy something from Aunt Polly at a garage sale. The problem is that the `$1` type of data points aren't really valuable to me, as they don't really indicate value--Aunt Polly didn't know any better, and didn't care. Similarly, I may occasionally get a data point coming from a jokester entering `$9000` for a toy that is really only worth `$9`. So, when computing value, what's the best way to factor these types of anomalies out of otherwise useful data? I've read about outliers, and something about generally ignoring anything that is more than 2.5 standard deviations outside the rest of the data, but I'm looking for the full recipe, here. Thanks so much!


title: In GPML, why does scaler in Squared Exponent covariance function gets very small?

text: I'm running regression using GPML with a covariance function being a sum of a Gaussian noise and a Squared Exponent (SE). Input is in R4 and both the input and the output are normalized. I run optimization of hyper-parameters by minimizing the negative log marginal likelihood. covSEiso is k(x_i,x_j)=sf^2*exp-(x_i-x_j)^2/(2*ell) After optimization, ell is reasonable (1.28), but sf^2 is extremely small (5.6445e-005). The noise variance is also reasonable (0.3715). Therefore, the covariance matrix is practically 0 and thus the regression result is just the prior mean. What in the data or elsewhere could be causing the scaler (sf) to become so small? Note that it's the WHOLE covariance matrix that is almost zero, not JUST the covariance entries which would indicate lack of "relationship" between different entries; the variance is also practically zeroed out. Thx, David


title: Reporting regression statistics after logarithmic transformation

text: I'm a bit troubled about how to report linear regression statistics after log transformation of the dependent variable. I suppose I should report the transformed coefficient, but would they be easily interpreted? In the regression plot should I use the transformed values or the original ones? while the seconds would be clearer than it would make no sense to plot the regression line produced after the transformation. Any direction?


title: What transfer function should I use to model a gradual permanent step positive intervention with increasing rate of growth?

text: ![enter image description here](http://i.stack.imgur.com/xl4OI.png) Hello, I want to ask what transfer function I should use to model a gradual permanent step intervention with increasing rate of growth. The relevant graph is presented at the attached image. I want to model a step intervention that leads to an increase in the level of the series. I consider three cases, an abrupt (sudden) permanent step intervention, a gradual permanent positive intervention with decreasing rate of growth and a gradual permanent positive intervention with increasing rate of growth (see graph). In all the cases I create a binary dummy variable to model the intervention. 1) In the simple case of the abrupt (sudden) permanent step intervention I run the regression Yt = intercept + [omega] * Dummy. Example data are as follows: … 20 20 33.5 33.5 33.5 33.5 33.5 33.5 33.5… 2) In the case of a gradual permanent step positive intervention with decreasing rate of growth I run the regression Yt = intercept + [omega/ (1-delta*B)] * Dummy. Example data are as follows: … 20 25 28 30 32 33 33.5 33.5 33.5 … 3) In the case of a gradual permanent step positive intervention with increasing rate of growth what transfer function should I use to model it? Yt=intercept+[Transfer Function] * Dummy Example data are as follows: … 20 20 20.25 21.5 23 25 29 33.5 33.5 … Thanks in advance, Andreas


title: Time Series Theorem

text: Consider the following theorem: > **Theorem.** For any sequence of autocovariances $\{ \lambda_k \}$ generated > by a stationary process, there exists a function $F$ such that $$\lambda_k = > \int_{(-\pi, \pi]} e^{ik \lambda} dF(\lambda)$$ where $F$ is the unique > function on $[-\pi, \pi]$ satisfying > >   1. $F(-\pi) = 0$ > >   2. $F$ is non-decreasing and right continuous > >   3. For any $0 \leq a \leq b \leq \pi$ we have $F(b)-F(a) = F(-a)-F(-b)$. > > > > Then $F$ is called the spectral distribution function. So basically, given a set of autocovariances (e.g. $\{E(X_{t}X_{t+1}), E(X_{t}X_{t+2}), \dots, E(X_{t}X_{t+10}) \}$), we can express each of these autocovariances in terms of the spectral distribution function? **Added.** What is the purpose of the this theorem?


title: Simultaneous confidence interval for multiple proportions in R

text: I want to perform in R the analysis "simultaneous confidence interval for multiple proportions", as illustrated in the article of Agresti et al. (2008) "Simultaneous confidence intervals for comparing binomial parameter", Biometrics 64, 1270-1275. Unfortunately I did not succed to find any R example, and I have no idea how to implement this process in R. As a case study, I need to apply such analysis on the following simple prolbem: I did an experiment in which 12 participants had to choose between 3 conditions when provided with 3 stimuli.               Stimulus  Condition1  Condition2 Condition 3     A            9          1          2     B           10          2          0     C            8          2          2      I want to prove that it is not by chance that Condition 1 is preferred rather than the other two conditions. Can you please give me an example in R on how to perform such analysis? In addition I would need to know how would you report such analysis in a paper. Thanks in advance


title: Do statisticians assume one can't over-water a plant, or am I just using the wrong search terms for curvilinear regression?

text: Almost everything I read about linear regression and GLM boils down to this: $y = f(x,\beta)$ where $f(x,\beta)$ is a non-increasing or non-decreasing function of $x$ and $\beta$ is the parameter you estimate and test hypotheses about. There are dozens of link functions and transformations of $y$ and $x$ to make $y$ a linear function of $f(x,\beta)$. Now, if you remove the non-increasing/non-decreasing requirement for $f(x,\beta)$, I know of only two choices for fitting a parametric linearized model: trig functions and polynomials. Both create artificial dependence between each predicted $y$ and the entire set of $X$, making them a very non- robust fit unless there are prior reasons to believe that your data actually are generated by a cyclical or polynomial process. This is not some kind of esoteric edge case. It's the actual, common-sense relationship between water and crop yields (once the plots are deep enough under water, crop yields will start diminishing), or between calories consumed at breakfast and performance on a math quiz, or number of workers in a factory and the number of widgets they produce... in short, almost any real life case for which linear models are used but with the data covering a wide enough range that you go past diminishing returns into negative returns. I tried looking for the terms 'concave', 'convex', 'curvilinear', 'non- monotonic', 'bathtub', and I forget how many others. Few relevant questions and even fewer usable answers. So, in practical terms, if you had the following data (R code, y is a function of continuous variable x and discrete variable group):               updown<-data.frame(y=c(46.98,38.39,44.21,46.28,41.67,41.8,44.8,45.22,43.89,45.71,46.09,45.46,40.54,44.94,42.3,43.01,45.17,44.94,36.27,43.07,41.85,40.5,41.14,43.45,33.52,30.39,27.92,19.67,43.64,43.39,42.07,41.66,43.25,42.79,44.11,40.27,40.35,44.34,40.31,49.88,46.49,43.93,50.87,45.2,43.04,42.18,44.97,44.69,44.58,33.72,44.76,41.55,34.46,32.89,20.24,22,17.34,20.14,20.36,24.39,22.05,24.21,26.11,28.48,29.09,31.98,32.97,31.32,40.44,33.82,34.46,42.7,43.03,41.07,41.02,42.85,44.5,44.15,52.58,47.72,44.1,21.49,19.39,26.59,29.38,25.64,28.06,29.23,31.15,34.81,34.25,36,42.91,38.58,42.65,45.33,47.34,50.48,49.2,55.67,54.65,58.04,59.54,65.81,61.43,67.48,69.5,69.72,67.95,67.25,66.56,70.69,70.15,71.08,67.6,71.07,72.73,72.73,81.24,73.37,72.67,74.96,76.34,73.65,76.44,72.09,67.62,70.24,69.85,63.68,64.14,52.91,57.11,48.54,56.29,47.54,19.53,20.92,22.76,29.34,21.34,26.77,29.72,34.36,34.8,33.63,37.56,42.01,40.77,44.74,40.72,46.43,46.26,46.42,51.55,49.78,52.12,60.3,58.17,57,65.81,72.92,72.94,71.56,66.63,68.3,72.44,75.09,73.97,68.34,73.07,74.25,74.12,75.6,73.66,72.63,73.86,76.26,74.59,74.42,74.2,65,64.72,66.98,64.27,59.77,56.36,57.24,48.72,53.09,46.53),                        x=c(216.37,226.13,237.03,255.17,270.86,287.45,300.52,314.44,325.61,341.12,354.88,365.68,379.77,393.5,410.02,420.88,436.31,450.84,466.95,477,491.89,509.27,521.86,531.53,548.11,563.43,575.43,590.34,213.33,228.99,240.07,250.4,269.75,283.33,294.67,310.44,325.36,340.48,355.66,370.43,377.58,394.32,413.22,428.23,436.41,455.58,465.63,475.51,493.44,505.4,521.42,536.82,550.57,563.17,575.2,592.27,86.15,91.09,97.83,103.39,107.37,114.78,119.9,124.39,131.63,134.49,142.83,147.26,152.2,160.9,163.75,172.29,173.62,179.3,184.82,191.46,197.53,201.89,204.71,214.12,215.06,88.34,109.18,122.12,133.19,148.02,158.72,172.93,189.23,204.04,219.36,229.58,247.49,258.23,273.3,292.69,300.47,314.36,325.65,345.21,356.19,367.29,389.87,397.74,411.46,423.04,444.23,452.41,465.43,484.51,497.33,507.98,522.96,537.37,553.79,566.08,581.91,595.84,610.7,624.04,637.53,649.98,663.43,681.67,698.1,709.79,718.33,734.81,751.93,761.37,775.12,790.15,803.39,818.64,833.71,847.81,88.09,105.72,123.35,132.19,151.87,161.5,177.34,186.92,201.35,216.09,230.12,245.47,255.85,273.45,285.91,303.99,315.98,325.48,343.01,360.05,373.17,381.7,398.41,412.66,423.66,443.67,450.39,468.86,483.93,499.91,511.59,529.34,541.35,550.28,568.31,584.7,592.33,615.74,622.45,639.1,651.41,668.08,679.75,692.94,708.83,720.98,734.42,747.83,762.27,778.74,790.97,806.99,820.03,831.55,844.23),                        group=factor(rep(c('A','B'),c(81,110))));          plot(y~x,updown,subset=x<500,col=group);      ![Scatterplot](http://i.stack.imgur.com/yeMsy.png) You might first try a Box-Cox transformation and see if it made mechanistic sense, and failing that, you might fit a nonlinear least squares model with a logistic or asymptotic link function. So, why should you give up parametric models completely and fall back on a black-box method like splines when you find out that the full dataset looks like this...               plot(y~x,updown,col=group);      My questions are:   * What terms should I search for in order to either find link functions that represent this class of functional relationships? or   * What should I read and/or search for in order to teach myself how to design link functions to this class of functional relationships or extend existing ones that currently are only for monotonic responses? or   * Heck, even what StackExchange tag is most appropriate for this type of question!


title: Best practices when treating range data as continuous

text: I am looking at whether abundance is related to size. Size is (of course) continuous, however, abundance is recorded on a scale such that               A = 0-10     B = 11-25     C = 26-50     D = 51-100     E = 101-250     F = 251-500     G = 501-1000     H = 1001-2500     I = 2501-5000     J = 5001-10,000     etc...       A through Q... 17 levels. I was thinking one possible approach would be to assign each letter a number: either the minimum, maximum, or median (ie A=5, B=18, C=38, D=75.5...). What are the potential pitfalls -- and as such, would it better to treat this data as categorical? I have read through this question which provides some thoughts -- but one of the keys of this data set is that the categories are not even -- so treating it as categorical would assume the difference between A and B is the same as the difference between B and C... (which can be rectified by using logarithm - thanks Anonymouse) Ultimately, I would like to see whether size can be used as a predictor for abundance after taking other environmental factors into consideration. The prediction will also be in a range: Given size X and factors A, B, and C we predict that Abundance Y will fall between Min and Max (which I suppose could span one or more scale points: More than Min D and less than Max F... though the more precise the better). Thanks.


title: Building a ML ordered logit regression model

text: I am building a ML ordered logistic regression. First of all, I really don't know if this is the best way to fit a model to my data, as I am not too confident in ML ordered logit regressions, compared to ML linear regressions. The relevant variables are as follows: $Y_{ij}$ = support for the European Union, where 0 = no support, 1 = some support and 2 = much support. $X_1$ = constitutional patriotism as a sumscore of three variables. Possible outcomes are $0,1,...,9$ $X_2$ = cultural patriotism as a sumscore of two variables. Possible outcomes are: $0,1,...,6$ The ML model consists of individuals nested with seven countries. Apart from the variables above, I have added a lot of control variables. First off, I want to create an interaction term including constitutional patriotism (CON) and cultural patriotism (CUL), where both variables are treated as continuous. The problem though, is that the two scales are not comparable. Hence, I thought about scaling the variables like this: $X_{1 scaled} = X_1/9$ and $X_{2 scaled} = X_2/6$ By doing it this way, I secure that both variables ranges from 0 to 1. Hence, a one unit increase in $X_1$ is comparable to a one unit increase in $X_2$ as they both have minimum = 0 and maximum = 1.   * Do you think it would be more appropriate to standardize the variables instead, so they have mean 0 and variance 1? Next I will have to choose between ML linear regression and ML ordered logistic regression. As said before, I am not too confident in the practical use of ordered logit. I wonder how much of a difference it makes, taking into account that ML ordered logit adds substantially more complexity to the model. I will, of course, choose a design that both makes sense, and delivers the best fit to the data. If I run a ML linear regression with my ordinal DV, the stadardized residuals are distributed as follows: ![Standardized residuals](http://i.stack.imgur.com/ZKgsf.png) This seems to be a problem, as the residuals are clearly not normal distributed.   * Does this call for a ordered logit, or could I work around the problem in some other way? Finally:   * Do you know of any methods to check how well my model fits the data in a ordered logit regression? I have thought about calculating the predicted distribution of individuals choosing $Y_{ij}=0$, $Y_{ij}=1$ and $Y_{ij}=2$ compared to the actual distribution found in the data. Also, I want to investigate if the ordered logit violates the proportional odds assumption. However, I am not sure how to do this in a ML framework. Any suggestions are much appreciated!


title: I have survey data where people respond to multiple items. I want to find the avg and SE on each item, controlling for within subject variation.

text: I have survey data from 650 respondents. Each participant rated 11 items on the same scale. I would like to know, at the population level, how the average of each of these items compare. For example, I would like to know that item 1 had an average of 2.4, item 2 had an average of 3.2, item 3 an average of 1.3, etc. I would also like to have confidence intervals around each of those items. I am imagining some subjects may just rate all items higher or have greater dispersion than other subjects. I would like to eliminate this within subject noise. I thought of a couple options. First, I could standardize each response within subjects, so each subject would have a mean of 0 and a standard deviation of 1. I was also thinking that maybe I could include subject level fixed effects, although that seems like it would only deal with the means being different, not the standard deviation. I also thought maybe I could cluster based on respondent, but that does not seem quite right and to be honest I am not very familiar with the technique. How would you suggest I deal with this? Your assistance is very much appreciated!


title: Visualising successive proportions

text: I am trying to visualise some consumer data, which has 4 categories. Users are free to switch between different categories. I would like to visualise the last three or four switches for each individual. So we would start with a plot with a column with 4 stacked proportions. After that we would have 16 as each category breaks down into what people did on the previous occasion, then 64, and so on, until the bins become too small to be useful. I am thinking somewhere between a marimekko chart and a stacked barchart or a dendro gram should work, but I don't even know what that would be called! If anyone can help with the type of plot I should be using, and, if you want to be extra nice, how to implement it in R then I would be very grateful.


title: How does Scheffe test deal with data snooping?

text: In Wikipedia > hypotheses suggested by the data, if tested using the data set that > suggested them, are likely to be accepted even when they are not true. This > is because circular reasoning would be involved: something seems true in the > limited data set, therefore we hypothesize that it is true in general, > therefore we (wrongly) test it on the same limited data set, which seems to > confirm that it is true. Generating hypotheses based on data already > observed, in the absence of testing them on new data, is referred to as post > hoc theorizing. > > The correct procedure is to test any hypothesis on a data set that was not > used to generate the hypothesis. For Post Hoc analysis of ANOVA, > **Henry Scheffé's simultaneous test** of all contrasts in multiple > comparison problems is the most[citation needed] well-known remedy in the > case of analysis of variance.1 It is a method designed for testing > hypotheses suggested by the data while avoiding the fallacy described above. So I was wondering how Scheffé's test avoids the fallacy of data snooping? Thanks!


title: How to combine time-series based features with different frequencies

text: I have 3 features which I want to use in my classifier. They are all time- series data-based. However, they are all at different frequencies and there have different matrix dimensions. I was wondering if someone could give me some pointers on how to combine these three features? Appreciate your help.


title: Combining models for prediction based on residual performance

text: I have never read or seen someone do this before, so I wanted to pose the question here. Suppose I fit a basic linear model, $\text{price of house} = \beta_0 + \beta_1*\text{taxes} + \beta_2*\text{bedrooms} + \beta_3*\text{number of garage doors}$. I fit this using OLS and I get a model that has combined $R^2$ of say .68... and I look at the residuals from the actuals vs predicted plot, and I see that the model does pretty well for houses that have less than 3 bedrooms, and maybe more than 2 garage doors.... I then build a robust linear regression model, or maybe a neural net, and I get a model that now performs worse as a whole, with a combined $R^2$ of say .53 BUT when I look at the residuals from the actuals vs predicted plot, I see that this model, model #2 was actually better a predicting the sales price for houses that model #1 did a poor job of predicting, say the outliers.... this scenario goes on and on, and lets say I have identified 5 different models, that for different parts of the population, do a better job of predicting those parts of the populations and these parts are not mutually exclusive. Is there are known approaches or good documentation on say using diagnostics to identify where your model may be poorly performing (on what types of data), and then some standardized approach to say joining predictions from multiple models to get the highest $R^2$? I am not suggesting averaging model predictions, or what is generally seen in ensembles, at least in how I have seen them defined. I am saying a tactic to apply different types of models to different parts of a dataset, (through semi-automatic identifcation of which model is best applied where), and then joining up all of the best predictions. I do not think segmenting the population first, and then building seperate models for each segment of the population is the best answer here, but is certainly one way to pseudo tackle it. Hopefully this concept makes sense, as I think it would be very powerful. Thanks to everyone in advance for providing your thoughts and or supporting documentation / literature.


title: Problem using rsFit to calcualte Hurst Exponent

text: I am trying to calculate the hurst exponent using the Rescaled range Analysis. I'm using the rsFit function of the fArma package. To test it I create a bronwian motion "r" and the hurst should be .5 but it is close to 1. Here is the R code:               library(fArma)          set.seed(100)     r<- cumsum(rnorm(500))     rs<-rsFit(r, levels = 50, minnpts = 3, cut.off = 10^c(0.7, 2.5),       doplot = TRUE, trace = FALSE, title = NULL, description = NULL)     rs      Any idea why the hurst exponent is close to 1 and not .5 for a pure random walk. Also - when i set r <- rnorm(1000) i.e. just white noise the hurst exp should be 0 but it is not. Thank you!


title: HMC: How many dimensions is too many?

text: From what I have read Hamiltonian Monte Carlo is the "goto" `MCMC` method when your problem is high dimensional. Practically speaking, how many dimensions 10's, 100's, 1,000's, 10,000's, 100,000's, ..., is too many? Computational cost will no doubt become an issue and I suppose the model being used is important to consider but those things aside, is there a practical limit to the number of dimensions when it comes to getting good samples from the desired distribution using `HMC`? Also, how can we monitor the convergence (or lack there of I guess) for problems where the number of parameters is far too many to check trace plots, running means, autocorrelations, etc. for individual parameters? **Update:** Found this post which mentions some non-visual diagnostics


title: Bhattacharyya distance for three histograms

text: There is a paper “Auto White Balance Based on the Similarity of Chromaticity Histograms” mention about automatic white balance. One of the key point of this algorithm is how to measure the similarity between three histograms (r, g, b channels). They choose Bhattacharrya distance to measure it. Equation of Bhattacharrya distance for two histograms (1) $$d(H_1, H_2) = \sqrt{1 - p(H_1, H_2)}$$ $$p(H_1, H_2) = H_1(0) * H_2(0) + H_1(1) * H_2(1) + … + H_1(255) * H_2(255)$$ However, no matter how hard I try, I can't find a way to apply Bhattacharrya distance on three histograms. Could I just alter the second equation to $$p(H1, H2, H3) = H_1(0) * H_2(0) * H_3(0) + H_1(1) * H_2(1) * H_3(1) + … + H_1(255) * H_2(255) * H_3(255)$$


title: SPSS linear regression - variables that are percentages

text: I have some data about music styles played at festivals. It looks like this: ![data](http://i.stack.imgur.com/SGXAa.png) Based on this (and other) data, I'd like to make a model to explain the price of a ticket for each festival. However I'm not sure how to handle this data in SPSS. Should I treat each style as a variable? And if I did, the styles aren't normal distributed. I can't take the natural logaritm to transform the data since there are zeros in every style. How should I proceed?


title: Explain the difference between multiple regression and multivariate regression, with minimal use of symbols/math

text: Are multiple and multivariate regression really different? What _is_ a variate anyways? :)


title: Normalizing data before applying MDS with strain criterion

text: The features of my dataset are like below:   `• BI-RADS assessment:` 1 to 5 (ordinal)   `• Age:` patient's age in years (integer) ranges from 18 to 96   `• Shape:` mass shape: round=1 oval=2 lobular=3 irregular=4 (nominal)   `• Margin:` mass margin: circumscribed=1 microlobulated=2 obscured=3 ill- defined=4 spiculated=5 (nominal)   `• Density:` mass density high=1 iso=2 low=3 fat-containing=4 (ordinal)      When I run `MDS` with `"strain"` criterion on such a dataset without normalizing it first, I get a result as follows: ![enter image description here](http://i.stack.imgur.com/sLg0B.png) However, if I normalize the data the result is as follows: ![enter image description here](http://i.stack.imgur.com/KmBzz.png)   The second results is pretty similar to results that I have got for other criteria and also for the `PCA` even I didn't normalize the data for them also. So, my question is: Why does normalizing data make difference for "strain" criterion? Thanks in advance...


title: Computation of distribution parameters for the maximum of two random variables

text: Let $X$ be a beta distributed variable with parameters $a$, $b$. Let $Y$ be a beta distributed variable with parameters $c$, $d$. Let $Z = \max(X, Y)$. Does anybody know of a fast way to compute the mean and variance of $Z$ given $(a, b, c, d)$? $X$, $Y$ don't necessarily need to have beta distributions, but they should be similar -- in $[0, 1]$. What I really want to find is a way to quickly approximate the distribution of the sample max of a bunch of random variables. Let $Z = \max(\{X_1, X_2,\dots , X_n\})$. If I know parameters (not necessarily beta-distribution) for $X_1, X_2,\dots, X_n$, can I quickly compute a set of parameters (same distribution as $X_i$) to form a distribution that closely approximates $Z$? * * * **EDIT:**   To clarify, the purpose of this has to do with programming AI for board games. If we are evaluating a position in a chess game and we determine one move leads to a 60% chance of winning while all others have a 20% chance of winning, then we value the position as a 60% win. This is the minimax algorithm... However, what I'm curious about is whether this can be improved upon. In the simple case, we approximate: $$ \mu_c = \max(\{\mu_1, \mu_2,\dots, \mu_n\}) $$ where $\mu_c$ is the expectation from the current position and $\mu_i$ is the expectation from each of the sub-trees. It seems like we should be able to do better than this. Can we find a function f that replaces max and closely approximates parameters for the distribution of the sample maximum? $$ (\mu_c, \sigma_c) = f\left(\{(\mu_1,\sigma_1), (\mu_2,\sigma_2), \dots, (\mu_n, \sigma_n)\}\right) $$ In this case, we still want $\mu_i$ to be the expectation from each of the nodes. However, $\sigma_i$ does not necessarily need to be variance/standard deviation; just a parameter that represents uncertainty. Just as: $$ \max(\{a, b, c\}) = \max\left(\{\max(\{a, b\}), c\}\right) $$ We should expect: $$ f\left(\{a, b, c\}\right) \approx f\left(\{f(\{a, b\}), c\}\right) $$ So $f$ does not necessarily need to take parameters for more than two nodes. One thing, however, is that $f$ should be relatively fast to compute. If the integrals must be calculated, I think doing double exponential integration might work well, but I'm not sure --   * http://www.johndcook.com/double_exponential_integration.html   * http://www.codeproject.com/Articles/31550/Fast-Numerical-Integration


title: When are Shao's results on leave-one-out cross-validation applicable?

text: In his paper Linear Model Selection by Cross-Validation, Jun Shao shows that for the problem of variable selection in multivariate linear regression, the method of leave-one-out cross validation (LOOCV) is 'asymptotically inconsistent'. In plain English, it tends to select models with too many variables. In a simulation study, Shao shows that even for as few as 40 observations, LOOCV can underperform other cross-validation techniques. This paper is somewhat controversial, and somewhat ignored (10 years after its publication, my chemometrics colleagues had never heard of it and were happily using LOOCV for variable selection...). There is also a belief (I am guilty of this), that it's results extend somewhat beyond the original limited scope. The question, then: how far do these results extend? Are they applicable to the following problems?   1. Variable selection for logistic regression/GLM?   2. Variable selection for Fisher LDA classification?   3. Variable selection using SVM with finite (or infinite) kernel space?   4. Comparison of models in classification, say SVM using different kernels?   5. Comparison of models in linear regression, say comparing MLR to Ridge Regression?   6. etc.


title: Why does the random forest OOB estimate of error improve when the number of features selected are decreased?

text: I am applying a random forest algorithm as a classifier on a microarray dataset which are split into two known groups with 1000s of features. After the initial run I look at the importance of the features and run the tree algorithm again with the 5, 10 and 20 most important features. **I find that for all features, top 10 and 20 that the OOB estimate of error rate is 1.19% where as for the top 5 features it is 0%.** This seems counter-intuitive to me, so I was wondering whether you could explain whether I am missing something or I am using the wrong metric. I an using the randomForest package in R with ntree=1000, nodesize=1 and mtry=sqrt(n)


title: Different use of neighbors in kNN classification algorithm

text: Suppose for each sample S in training set TS I already know the k nearest neighbors. The standard algorithms usually give S the same label as the majority of S's neighbors has. What other use of nearest neighbors could I make? For instance, if I have not 3, not 5, but more neighbors, would it make sense to create a temporary model based on them and only then classify S? If so, what model would you recommend?


title: Are there any examples where Bayesian credible intervals are obviously inferior to frequentist confidence intervals

text: A recent question on the difference between confidence and credible intervals led me to start re-reading Edwin Jaynes' article on that topic: Jaynes, E. T., 1976. `Confidence Intervals vs Bayesian Intervals,' in Foundations of Probability Theory, Statistical Inference, and Statistical Theories of Science, W. L. Harper and C. A. Hooker (eds.), D. Reidel, Dordrecht, p. 175; (pdf) In the abstract, Jaynes writes: "...we exhibit the Bayesian and orthodox solutions to six common statistical problems involving confidence intervals (including significance tests based on the same reasoning). In every case, we find the situation is exactly the opposite, i.e. the Bayesian method is easier to apply and yields the same or better results. Indeed, the orthodox results are satisfactory only when they agree closely (or exactly) with the Bayesian results. NO CONTRARY EXAMPLE HAS YET BEEN PRODUCED". (the SHOUTING being mine not Jaynes') The paper was published in 1976, so perhaps things have moved on, so my question is, are there examples where the frequentist confidence interval is clearly superior to the Bayesian credible interval (as per the challenge implicitly made by Jaynes). Examples based on incorrect prior assumptions are not acceptable as they say nothing about the internal consistency of the different approaches.


title: SVD in linear regression

text: I was reading the book Elements of Statistical Learning and came across the section that tried to interpret ridge regression using singular value decomposition (SVD) of the design matrix, $X$. Specifically, I found the following: $X=UDV^{T}$, where matrix $U$ is $N\times p$, $V$ is a $p\times p$ orthogonal matrix, and $D$ is a $p\times p$ diagonal matrix. I am confused because from Wikipedia, the orthogonal matrix has to be a square matrix. In this case matrix $U$ does not qualify. Later I tend to believe that $U$ contains orthogonal columns only, and that results in $U^{T}U=I$, but $UU^{T}\ne I$. This seems to make sense because I found in the book $X \hat{\beta}=X(X^{T}X)^{-1}X^{T}Y=UU^{T}Y$, and $UU^{T}Y$ should not be equal to $Y$ So my question becomes: are there two versions of SVD I can do? One results in both $U$ and $V$ being orthogonal and square matrix, and the other like this? Or is there anything wrong with my argument? Any guidance is appreciated. Update after receiving initial answer: After reading @BabakP 's answer, I thought testing the algorithm using software is a good idea. So I tried svd function in Matlab. The result shows a square U matrix in dimension NxN, a diagonal matrix D in dimension Nxp, and a square V matrix in dimension pxp. Example below:               A=[ones(10,1) randn(10,1)];     [U,S,V]=svd(A);     >> size(U)          ans =          10    10          >> size(S)          ans =          10     2          >> size(V)          ans =           2     2       So does this mean R and Matlab give two different versions?


title: Confidence Interval in Monte Carlo integration

text: I want to integrate $\int_{\mathbb{R}_+}\mathbb{1}_A(x) d\mathbb{P}(x)$, in other words I am interested in $\mathbb{P}(A)$. I did this numerically with two Monte Carlo steps. First, I drew, say a sample of $10^3$ ($\mathbb{P}$-distributed). From this sample I calculated the mean. Repeating this $N$ times, i got a new mean (better) and calculated the 95% confidence interval using the variance of these $N$ estimates for $\mathbb{P}(A)$ with the Student-t distribution. I see now, that some values do not fall into the confidence intervals (too many to be random). So I fear that my approach is flawed and underestimates the variance. Of course, there is always the possibility of having a wrong implementation. But then I increased the number of repetitions dramatically and it looked better. I am not sure what to do. What is the standard procedure to compute the confidence interval here? Maybe using a $\sigma$-error/Chebyshev-inequality, i.e. just one step, not two? I am grateful for any hints.


title: What can we conclude about the distribution of the sum of two random variables?

text: If we know, for independent random variables $X$ and $Y$, $P(X>x)\leq0.05$, and $P(Y>y)\leq0.05$, can we say anything about $P(X+Y>x+y)$? Can we be certain that it is less than $0.05$? Under what conditions could we say that?


title: Can I use z-scores to compare results from different measures?

text: Can I compare outcome data from different literacy measures with different scales (delivered to different subjects) by transforming the data into z scores? Ideally I'd like to use regression and ANOVA techniques to determine if there are significant differences between groups that got different types of literacy interventions. Is this remotely plausible? If not, is there another legitimate way to compare the scores from different measures?


title: Comparing confidence limits

text: I have a data set that is broken into the four categories, each with four data points. I am asked to compare the confidence limits of the four categories. How do I go about doing this? I know to create confidence limits you need t tests, however they only compare two means, not four.


title: GLM vs square root data transformation

text: I am currently analysing some pretty awful/awkward data on the abundance of fish under three different "Hydro-Regimes" (5 abundance measurements for each regime - Short/Medium/Long). The current analysis plan had been a one way ANOVA. Looking at the residuals vs. fitted and the normal Q-Q after plugging it into a linear models indicates right skew. To correct this, I tried to Sq. Rt transform the outcome variable, which gives a reasonable answer, with "short" being significantly different to "long" and "medium" after a post-hoc Tukey. However, plotting the data shows the error bars of "long" and "medium" not to overlap. I also looked at using a Poisson GLM, which I have read is good for count and skewed data, but I am not sure if this is the right way to go. Any suggestions?


title: Central limit theorem and the law of large numbers

text: I have a very beginner's question regarding the Central Limit Theorem (CLT): I am aware that the CLT states that a mean of i.i.d. random variables is approximately normal distributed (for $n \to \infty$, where $n$ is the index of the summands) or the standardized random variable would have a standard normal distribution. Now the Law of Large Number states roughly speaking that the mean of i.i.d random variables converges (in probability or almost surely) to their expected value. What I don't understand is: If, as the CLT states, the mean is approximately normally distributed, how then can it also converge to the expected value at the same time? Convergence would imply for me that with time the probability that the mean takes a value which is not the expected value is almost zero, hence the distribution would not really be a normal but almost zero everywhere except at the expected value. Any explanation is welcome.


title: Determining the "variability" of a benchmark

text: I have a software benchmark which is quite noisy. I am trying to for the bugs which are causing the noise, and I need to be able to measure it somehow. The benchmark is comprised of a number of subbenchmarks, for example:               "3d-cube": 31.56884765625,     "3d-morph": 21.89599609375,     "3d-raytrace": 51.802978515625,     "access-binary-trees": 15.09521484375,     "access-fannkuch": 45.578857421875,     "access-nbody": 8.651123046875,      The times are in milliseconds. The times typically vary between runs. For example, on my machine, the "3d-cube" benchmark tends to take around 35ms, but I've seen it go as high as 44ms, and 31ms (above) is uncharacteristically low. My aim is to change the benchmark so that minor improvements to the run-time can be visible in a benchmark result. What I need is a number that tells me whether I have reduced the "variability" of the benchmark. ### My own solution I run it the benchmark 1000 times, the took the sum of the differences between each subbenchmark's mean and its actual run-times. In pseudo-code:               v = 0     for s in subbenchmarks:       x = mean of all iterations of s       for i in iteration         v += absolute_value(results[s][i] - x)      I'm sure this isn't statistically valid (having asked someone), but what is a "correct" way of measuring this "variability" so that I can reduce it.


title: Extract Information from Traded volumes

text: Assume that we want to extract information from a constantly updating table like stock trading volume, where a banker can buy or sell at specific price. The table has the following 4 columns: price, available, traded, buy/sell. The fist column gives the price that the amount given by the 2nd column is available to be traded or has already been matched as given by the 3rd column. Finally the forth column indicates of the available amount is available to be bought or be sold. > price | available | traded | buy/sell This table is constantly updating so we have a time series of the various changes (3rd column only increases of course), is there any specific method that I could refer to, to extract information or make assumptions about the movements from the second column (available) to the third(traded), taking under consideration if its buy or sell?


title: Integrating Gaussian white noise over a Gaussian density

text: I have encountered the following integral: $$ \int_{- \infty}^{+ \infty} X(\theta) \frac{\exp \big \{-\frac{\theta^2}{2\sigma^2} \big \}}{\sigma\sqrt{2 \pi}}d\theta $$ where, for fixed $\theta = \theta_0$, $X(\theta_0) \sim N(0, 1)$. $X(\theta_1)$ and $X(\theta_2)$ are independent if $\theta_1 \neq \theta_2$. Hence $X(\theta)$ is Gaussian white noise indexed by $\theta$. I am not sure about how to interpret this expression. In particular, by thinking about the integral as the limit of a sum: $$ \frac{1}{N}\sum_{i=1}^N X(\theta_i), \;\;\;\;\;\;\; X(\theta_i) \sim N(0, 1) \;\;\;\;\text{and}\;\;\;\;\; \theta_i \sim N(0, \sigma) $$ for $i = 1, \dots, N$, I would say that it should converge to zero in probability, but I'm not completely sure.


title: Moderation with latent variables in Stata 13

text: I am wondering if anyone has insight regarding the creation of interaction terms (moderators) in a model with latent predictors using Stata 13 - SEM. I have what I want to be the moderator set up currently as a mediator and need to test competing models. I am not sure if I need to run a multiple group model or if there is a cleaner approach. Any thoughts would be helpful.


title: Direct Sampling from posterior distribution

text: Why is direct sampling from the posterior distribution intractable?


title: Invariance in neural networks

text: I have tried to read about tangent propagation in neural networks (although I guess it could be applicable to other methods) which is a procedure to create models that are invariant to certain transformation. So basically the argument goes like this: Start with an input vector $x_{n}$. Then use a linear transformation (maybe a rotation or a translation) $s(x_{n},\xi)$ such that when $\xi=0$, it returns the same input: $s(x_{n},0)=x_{n}$. So when we apply a given transformation $s(x_{n},\xi)$ we move $x_{n}$ in some way (according to the value of $\xi$.) If we plot every new value for $x_{n}$ depending on our choice of $\xi$ then we can draw a curve $M$. Then we ask for the tangent vector of this curve in the original point $\xi = 0$, which is given by: $$\tau_{n}=\displaystyle \frac{\partial s(x_{n},\xi)}{\partial \xi}\bigg|_{\xi=0}$$ At this point, I think I understand what is going on. However, I can't justify the following calculation: Let's suppose I want to know the derivative of the output vector $y_{k}$ with respect to $\xi$ at the original point $\xi=0$. That's a fair question to ask since the output vector depends on the input vectors. However, Bishop's book (page 264) says that $$\frac{\partial y_{k}}{\partial \xi}\bigg|_{\xi=0}=\sum_{i=1}^{D}\frac{\partial y_{k}}{\partial x_{i}}\frac{\partial x_{i}}{\partial \xi}\bigg|_{\xi=0}=\sum_{i=1}^{D}J_{ki}\tau_{i}$$ where $J_{ki}$ is the $k,i$ element of the Jacobian matrix. As I understand it, we are replacing $x_{n}$ by its transformation $s(x_{n},\xi)$, so I think $y_{k}$ depends on $s(x_{n},\xi)$ and $s(x_{n},\xi)$ (as it indicates) depends on $x_{n}$ and $\xi$. Actually, I don't understand how the expression above can be correct, since $\tau_{i}$ is defined in the first equation based on $s(x_{n},\xi)$, instead of $x_{i}$ What is the reason for this description of the change of $y_{k}$ when we change a bit $\xi$? _Note: I found out the originalpaper by Simard in which it's explained the general idea but unfortunately he doesn't carry out the calculation of $\frac{\partial y_{k}}{\partial \xi}\bigg|_{\xi=0}$, so no help from there._


title: Aggregation of correlated variables

text: I've been trying to aggregate _correlated_ time series, by using Alexander's proposal that you can see here: http://bit.ly/1hIPwiI. Her proposal to find a random variable $Y=\sum_{i=1}^N X_i$, where each $X_i$ represents a stochastic process, is based on solving an integral of dimension $N-1$. The problem that I see is that the integral must be solved by iterative methods (I use Monte Carlo) when the dimension is high. Even so, when $N>5$ my results are really disappointing; with large errors in some point estimates. The question that I ask here is, Is there any other method to solve the problem of portfolio aggregation (summation of random variables) when the variables are correlated? In other words, if I have the marginal distributions of $X_1, \ldots, X_N$, how can I find the distribution of$Y$?


title: Minimum number of observations for logistic regression?

text: I'm running a binary logistic regressions with 3 numerical variables. I'm suppressing the intercept in my models as the probability should be zero if all input variables are zero. What's minimal number of observations I should use?


title: A better bayesian way of modelling autoregressive mixtures

text: I have a JAGS hierarchical model which includes a temporal sub-model for the primary vote share between four party groups (LNP, Labor, Green, and Other). For each day in the temporal model, the vote shares sum to one. The temporal model is an auto-regressive random walk: the vote shares for the four party groups today are much the same as they were yesterday.               #### -- temporal model (a daily walk where today is much like yesterday)     for(day in 2:PERIOD) { # rows         for (party in 1:PARTIES) { # columns             tmp[day, party] ~ dnorm(walk[day-1, party], walkPrecision[party])         }     }          ## -- impose a sum-to-one constraint ... total of all parties sums to one every day     for(day in 1:PERIOD) { # rows         walk[day, 1:PARTIES] <- tmp[day, 1:PARTIES] / sum(tmp[day, 1:PARTIES ])     }          ## -- constrained priors for the day-to-day variance of the temporal model     for(party in 1:PARTIES) { # for each party         sigmaWalk[party] ~ dunif(0, 0.005)  ## uniform prior on std. dev.           walkPrecision[party] <- pow(sigmaWalk[party], -2)        }          ## -- uninformative priors for first day in the temporal model     for (party in 1:PARTIES) { # for each party         tmp[1, party] ~ dunif(0.0001, 0.9999) # fairly uninformative     }      My question - is there a better way of modelling this than I have above.


title: How to combine two independent repeated experiments with different success probabilities?

text: Repeating an experiment (about which I asked before) with $n$ possible outcomes $t$ times independently, where all but one outcomes have probability $\frac{1}{n+1}$ and the other outcome has the double probability $\frac{2}{n+1}$, I also obtain an independent information where the outcome with double probability shows up 50% more likely than the other outcomes (i.e., with probability $\frac{3}{2n+1}$ and the others' probability is $\frac{2}{2n+1}$). _How do I combine the two results?_ **Alternative formulation:** Given is the probability space $(\Omega\times\Omega, \mathcal{P}(\Omega\times\Omega), \mathrm{p})$ with $$\mathrm{p}(\omega, \omega') = \left\{\begin{array}{cc} \frac{6}{(n+1)(2n+1)} & \mbox{ if } \omega = \omega_0 \mbox{ and } \omega' = \omega_0\\ \frac{4}{(n+1)(2n+1)} & \mbox{ if } \omega = \omega_0 \mbox{ and } \omega' \ne \omega_0\\ \frac{3}{(n+1)(2n+1)} & \mbox{ if } \omega \ne \omega_0 \mbox{ and } \omega' = \omega_0\\ \frac{2}{(n+1)(2n+1)} & \mbox{ otherwise } \end{array}\right.$$ where $\omega_0\in\Omega$ is unknown (and $|\Omega| = n$). The goal is to find $\omega_0$ given $t$ samples from $\Omega\times\Omega$. Currently I'm counting how often each value shows up in the first coordinate and add the number of times it shows up in the second coordinate multiplied by a weighting factor of $\log_2(\frac{3}{2})$. The value with the highest number should be $\omega_0$ (if $t$ is big enough). _Is this the correct weighting factor rsp. the correct way to find $\omega_0$?_ PS: I'm also thankful for anyone finding better tags for my question.


title: Selecting trustworthy rows based on minima and maxima

text: I'm quite new to this, but I have geographical data and I want to calculate latitudinal shifts. I have in each row a minimum, maximum, median, average and standard deviation of latitude BEFORE 1980. In the same row, I have the same variables but from AFTER 1980. Goal is to select the rows that have trustworthy values. I want to compare the data before and after 1980 that don't have values that are too different in terms of minimum and maximum latitude. Tricky thing is that latitude can be positive OR negative.. Another thing is that the goal is to look for shifts in latitude, so differences in minimum and maximum should be allowed. I'm looking for a method/formula to delete rows that have for example a minimum latitude before 1980 of -23, and after 1980 of +10 (which would mean that there is an error or bias of some kind).


title: Will there be Volume 2 of Mathematical Statistics by Bickel and Doksum and Robust Asymptotic Statistics by Rieder?

text: 1. It was mentioned that Volume 2 of Mathematical Statistics by Bickel and Doksum was scheduled to be published in 2003. I was wondering if I miss something, but why it is not yet available? Is there some change of plan?   2. Volume 1 of Robust Asymptotic Statistics by Helmut Rieder was published in 90's. But it seems no volume 2 yet. Am I correct? Thanks?


title: How to analyze pre-post attitudinal data (ordinal) collected from a group of 60 students?

text: I have a survey with 25 items that were administered pre and post an intervention. Each item had four possible responses:   * Not Appealing at all   * Not Appealing   * Somewhat Appealing   * Very Appealing. The survey was given to a class of 60 students. I have entered the data into SPSS and created a case for each student, with Pre and Post responses. As the data are ordinal in nature it seems like I should use a non-parametric test. Is the Wilcoxon-Mann-Whitney test appropriate? Thanks.


title: Maximization of Log-likelihood Numerically and Local Maximum Problem

text: I am maximizing the log-likelihood function for a generalized linear model without using per-written functions like `glm` due to the form of the model I have. So to maximize the log-likelihood, I am using the `optim` function. That works fine. But how can I double check if the estimated parameters are actually occurring at the global maximum and not at a local maximum? Is there any method that guarantees the global maximum?


I'm trying to fit a Dirichlet process Gaussian mixture model to 2-D data using a stick-breaking prior. I have a normal-inverse Wishart prior on my cluster means and covariances, where $$\{\mathbf{\mu_i}, \mathbf{\Sigma_i}\} \sim \mathcal{N}\mbox{IW}(\mathbf{\theta_0}, \lambda_0, \mathbf{\Sigma_0}, \nu_0) \leftrightarrow \mathbf{\mu_i} \sim \mathcal{N}(\mathbf{\theta_0}, \lambda_0^{-1}\mathbf{\Sigma_i}), \mathbf{\Sigma_i}\sim \mbox{IW}(\mathbf{\Sigma_0}, \nu_0)$$ To generate data, I set $\mathbf{\theta_0} = \mathbf{0}$, $\lambda_0 = 0.01$, $\mathbf{\Sigma_0}=\mathbf{I}$, and $\nu_0=10+\dim{\mathbf{\theta_0}}$. I'm just playing with a toy example so $\dim{\mathbf{\theta_0}}=2$. I'm learning the means and covariances through Gibbs sampling and am setting the model hyperparameters to the true hyperparameters. I've noticed that whenever my cluster means are far from $\mathbf{0}$, my model learns the correct number of clusters and learns the correct means and covariances. On the other hand, when I have one or more clusters that are close to $\mathbf{0}$, the learning algorithm uses more means than necessary to model the cluster that's close to $\mathbf{0}$ and cannot sample out. Moreover, the means that are not used to represent clusters congregate very tightly around $\mathbf{0}$ (regardless of whether a cluster is close to $\mathbf{0}$). Could someone please explain why this happens and how to get around the problem of having more means than necessary to represent one cluster in this simple toy example?