title: Conditional expectation for the reciprocal of a normal

text: We know the expected value of $1 \over X$, where $X$ is a normal random variable, does not exist. But suppose we condition on an interval not containing zero. For example, if $ \mu_X = 10$ and $\sigma_X = 1$, what is $E [{1 \over X} | ( 9<X<10 )] $ ? I know numerical integration or an approximation can be used, but can an analytic formula be found? If so, what technique is used to evaluate the integral?


title: How to use G Power 3 to calculate statistical power in mixed design ANOVA with unequal group sample sizes

text: In G power 3, ANOVA repeated measures within-between interaction: Only the total sample size is reported assuming equal sample size for the two groups. My questions are:   1. How would it work if the sample sizes are slightly different, for example: N1/ N2 = 1.16.    2. I have to input correlation between repeated measures. Is this the correlation between repeats after merging the data from the two groups?    3. Information about non-centrality parameter would be helpful.


title: Performance decay when testing on bootstrap data

text: I have a strategy that has a Sharpe ratio of 1.6 when back tested over the past 10 years. When I run this same strategy on re-sampled data, the performance of the strategy goes down to 1.32. Should I derive the expected return from the performance of the bootstrap results? Does the difference in performance tell me anything else?


title: SVM kernels choose without tests

text: I am having a big trouble to determine the kernel i should use in a non linear SVM without testing before, I want to know if there is any other ways to determine the best kernel without tests ? Is it linked to the data we are working on ?


title: What does it imply if accuracy and recall are the same?

text: I did a number of machine learning experiments to predict a binary classification. I measured precision, recall and accuracy. I noticed that my precision is generally quite high, and recall and accuracy are always the same numbers. I used the following definitions: $\text{Precision} = \frac{TP}{(TP + FP)}$ $\text{Recall} = \frac{TP}{(TP + FN)}$ $\text{Accuracy} = \frac{(TP + TN)}{(P + N)}$ I have some difficulties to interpret accuracy and recall. What does it mean if these two number are always the same in my case?


title: Reporting change when the baseline varies

text: I am discussing how growth of a business changes over a time period with a colleague. The issue is, that I want to extrapolate the impact of business investment beyond maintenance. The example is a hotel chain with n hotels in year one. In year two have $n + x - y$ hotels, where $x$ are new openings and $y$ are sold or closed hotels in the reference year. To compare the two years like for like I am tempted to just report on $n - y$, but this does not strike me as sensible as it will not consider effects such as new hotels making more revenue and cannibalizing some of the old hotels etc. Is there some common practice how these effects are usually mitigated?


title: How to create a score that gives same weight to N variables using different scales?

text: I'm analyzing a website, we have three variables Pageviews, Minutes spent on site and Entrances and want to produce a score that gives equal weight to each of them. At first I was simply going to multiply them but this would give different weight to each:               PageViews * Minutes * Entrances = Score     Score1 = 27 * 1.70 * 4 = 45.90      My second idea was to normalize the values dividing them by the maximum value found in each case:               ((pv / maxPV) + (ms / maxMs) + (e / mE)) / 3 = Score     Score1 = ((27 / 175) + (1.70 / 4.02) + (4 / 169)) / 3 =     Score1 = (0.15 + 0.42 + 0.02) / 3 = 0.19      Is this a good way to normalize the variables? If not how could I create a score that gives equal weight to each? Also how could I make one variable weight more or less than the others, or have a negative impact on the score?


title: Is it a truly random sample?

text: Any hints on how to solve this problem? > A sample of 400 items is taken from a normal population whose mean as well > as variance is 4. If the sample mean is 4.5, can the sample be regarded as a > truly random sample?


title: Multiple treatment nonparametric testing

text: I am a SAS programmer and am trying to help out a mate doing a project for work. He has tested 7 different treatments (growth cultures) across 31 days, I have a value from each treatment for each of the 31 days. The days go from Day 0 to Day 31, and the nature of the testing is that the values at day one are 0 and increase at day 31 to +1000. I have tried fitting distributions to the data but none are statistcally fitting (Minitab p values all <0.05). So we can't assume normality and use a simple ANOVA with Tukey's Studentized Range (HSD) Test. I am going to test these using non-parametric testing, using the Wilcoxon rank test on the differences between each treatment, so using day as variable to pair the data, use the differences and then test for a difference in median. The outcome of these tests is that all treatment medians are significantly different but also most of the distributions are significantly different using Kolmogorov-Smirnov Two-Sample Test (Asymptotic) testing. Does this mean that we can only use the significant values for median differences when the distributions are equal or does it impact the analysis at all? Can anyone suggest any other testing methods I could use here?


title: Understanding a paper on Gaussian Processes

text: I am reading and trying to understand the paper "Gaussian Processes for Signal Strength-Based Location Estimation" by Ferris et al. (A copy of the paper may be accessed here). At the moment I am stuck with Section II A. All notation used below is consistent with the paper.   1. The authors introduce two covariance functions in equations 3 and 4. Can anyone explain the difference between these two functions? (i.e. why they are useful.)   2. In equations 5 and 6, the authors use $\textbf{k}_*$, which is "the $n\times1$ vector of covariances between $x_*$ and the training inputs $\textbf{X}$," I don't understand how to compute it though. Any insights? (Is it found from equation 4?) I'm not very familiar with the theory of Gaussian Processes, and unfortunately don't have the luxury of time to spend understanding the theory. I would appreciate answers that focus on the "functional" aspects of the paper, as that is the reason I am reading it.


title: What is an adaptive copula?

text: My basic question is: What is an adaptive copula? I have slides from a presentation (unfortunately, I cannot ask the author of the slides) about adaptive copulae and I am not getting, what this means resp. what this is good for? Here are the slides: ![sl1](http://i.stack.imgur.com/0F76A.png) ![sl2](http://i.stack.imgur.com/F3H0r.png) Then the slides continue with a change-point Test. I am wondering what this is about and why I need this in connection to copulae? The slides end with an adaptively estimated parameter plot: ![sl3](http://i.stack.imgur.com/qJXPm.png) ![sl4](http://i.stack.imgur.com/jYIy9.png) This seems to show, that my estimates are lagged behind. Any other interpretations, comments would be great!


title: Good introductory biology/genetics readings for statisticians?

text: I'm looking for books on genetics or biology in general, at an introductory level, ideally where topics are tackled from a statistical point of view. I've been reading "Genes population and languages" by Luca Cavalli-Sforza these last few days and I'm enjoying it a lot, something similar would be perfect. To give you some background information, I'm a graduate student in biostatistics, recently switched to biostats from economics, so I don't know much about biology in general. I'm trying to get an idea of the field in order to decide wether to take courses in statistical genetics or spatial analysis. Thank you!


title: Difference estimate & confidence intervals for $\chi^2$ test between 2 proportions

text: I am using a chi-sq test between two proportions (http://statistic-on- air.blogspot.com/2009/07/comparison-of-two-proportions.html). I am using this over the z-test for proportions because I do not think my data is normal. My question is how do I get the difference (shift) estimate & confidence intervals for this test. For example,               Data1: 193/252=.77       Data2: 154/227=.68      P-value=.032 (I used the formula in the link above. I will also add my exact code below) Shift Estimate: I assume its the difference between proportions, ie .77-.68=.09. Is this correct? Confidence Interval: ?????????????????? Code for generating P-value (matlab):               % Pooled estimate of proportion      p0 = (n1+n2) / (N1+N2);     % Expected counts under H0 (null hypothesis)      n10 = N1 * p0;      n20 = N2 * p0;     % Chi-square test, by hand      observed = [n1 N1-n1 n2 N2-n2];      expected = [n10 N1-n10 n20 N2-n20];      chi2stat = sum((observed-expected).^2 ./ expected);      p = 1 - chi2cdf(chi2stat,1);      H=0; if(p<.05), H=1; end


title: How can I form ARIMA equation given MA and AR terms

text: ![enter image description here](http://i.stack.imgur.com/iRvX5.png) Above is output from SAS. What would be the corresponding ARIMAX equation? I would appreciate if someone could help me write the mathematical equation, preferably in the following form: $$ Y(t)= ay(t-1) + by(t-2) + \ldots + z $$ where $a, b, c$ are the coefficients, and $z$ is any error term.


title: Winning probability in a game with multiple players

text: I do not come from a mathematical background, and hope you can answer this (probably very basic?) question. I got a game group with friends where we play board games, normally 3-6 players each night. I have created an application which keeps track of game history, and keeps track of a Elo rating for each player. These Elo ratings can be used to determine the winning probability between two players. However, how can I use these winning probabilities to get the probability of a player winning a multi-way game? Let's say we have a three way game, two players with a Elo rating of 1200 and one with 1400. The two with 1200 have 50% chance of winning against each other and the 1400-player has a 76% chance of winning in a heads up match against a 1200-player. How can I find the probabilities that each player wins the three- way game?


title: In the output of lsmeans, is “estimate” is effect size of that contrast?

text: I am using lsmeans to run post-hoc analysis for contrasts of linear mixed model, and here is the output:               $contrasts      contrast            estimate       SE    df t.ratio p.value      aL2_like - english  450.7397 252.6889 75.99   1.784  0.1819      aL2_like - thai     224.6113 127.8683 59.66   1.757  0.1931      english - thai     -226.1284 274.6949 96.77  -0.823  0.6896      Can I ask if the estimate here refers to the 'effect size' of each contrast?


title: How can we evaluate the predicted values using Scikit-Learn

text: I am using AdaBoost Classifier to predict values I have. How can evaluate the accuracy of prediction model (I'd like to see how the accuracy of predicted values). You can check an example here: http://scikit- learn.org/stable/modules/ensemble.html#usage I found two options : using confusion matrix               from sklearn.metrics import confusion_matrix     cm = confusion_matrix(expected, y_1)      or using cross val score               scores = cross_val_score(clf_1, X_train, y_train)     print scores.mean()      There is also: AdaBoostClassifier.staged_score(X, y) AdaBoostClassifier.score(X, y) So, I am little bit confused. One last question: Should I use predict() or predict_proba().


title: Training a Hidden Markov Model, multiple training instances

text: I've implemented a discrete HMM according to this tutorial http://cs229.stanford.edu/section/cs229-hmm.pdf This tutorial and others always speak of training a HMM given an observation sequence. What happens when I have multiple training sequences? Should I just run them sequentially , training the model after the other? Another option is to concatenate the sequences to one and train on it, but then I will have state transitions from the end of one sequence to the start of the next one which are not real. Thanks.


title: Clustering groups that have replicated measures: hierarchical clustering on group-average VS regression tree

text: I measured 2 continous dependent variables (V1 and V2) on 10 occasions (10 replicates) for each of 4 groups. I aim to cluster my groups. i.e. I dont want to cluster replicates, since this could mix groups inside a cluster. I see two options:   1. hierarchical clustering on group-average (V1mean and V2mean)   2. multiple regression tree(MRT) with replicated measures as dependent variable and group as independent variable NB1: I am not interested in the hierarcical structure. Kmeans partitionning on group average would be also OK for me. I guess the MRT solution is more usefull since it give an idea on the significance of the resulting clusters. However, I never saw multiple regression tree used like that. Could you please confirm that I am not wrong. Here is a working exemple with R. Data simulation and visualisation:               gp=gl(4,10)     V1=rnorm(40,c(rep(1,10),rep(4,10),rep(8,10),rep(7,10)),rep(2,40))     V2=rnorm(40,c(rep(1,10),rep(4,10),rep(8,10),rep(7,10)),rep(2,40))     levels(gp)=paste("g",levels(gp),sep="")     df=cbind.data.frame(gp,V1,V2)     summary(df)     f.se=function(x)sd(x)/sqrt(length(x))     df.m=ddply(df,.(gp),summarize,V1m=mean(V1),V1se=f.se(V1),V2m=mean(V2),V2se=f.se(V2))          ggplot(df,aes(V1,V2,color=gp))+geom_point()+       geom_point(data=df.m,aes(V1m,V2m,color=gp),shape=4)+       geom_errorbarh(data=df.m,aes(V1m,V2m,xmin=V1m-V1se,xmax=V1m+V1se,color=gp),shape=4)+       geom_errorbar(data=df.m,aes(V1m,V2m,ymin=V2m-V2se,ymax=V2m+V2se,color=gp),shape=4)      Option 1: hierarchical clustering on group average               htree=hclust(dist(apply(df.m[,-1],2,scale),"euclidean"))     plot(htree)      Option 2: MRT               mrt=mvpart(data.matrix(apply(df[,-1],2,scale))~gp,data=df,xv="min")      We can see that MRT suggests a difference between gp 1 et 2, which is not obvious as regards the results of hierarchical clustering


title: Combinatorics: picking 4 numbers such that sum of 2 is equal to sum of other 2

text: Bus tickets in certain city contains four numbers u,v,w,x. Each of these numbers is equally likely to be any of the 10 digits 0,1,2....,9 and four numbers are chosen independently. A bus rider is said to be lucky if u+v=w+x. What proportion of riders are lucky?


title: Has anyone solved PTLOS exercise 4.1?

text: This an exercise given in **Probability Theory: The Logic of Science** by Edwin Jaynes, 2003. There is a partial solution here. I have worked out a more general partial solution, and was wondering if anyone else has solved it. I will wait a bit before posting my answer, to give others a go. Okay, so suppose we have $n$ mutually exclusive and exhaustive hypothesis, denoted by $H_i \;\;(i=1,\dots,n)$. Further suppose we have $m$ data sets, denoted by $D_j \;\;(j=1,\dots,m)$. The likelihood ratio for the ith hypothesis is given by: $$LR(H_{i})=\frac{P(D_{1}D_{2}\dots,D_{m}|H_{i})}{P(D_{1}D_{2}\dots,D_{m}|\overline{H}_{i})}$$ Note that these are conditional probabilities. Now suppose that given the ith hypothesis $H_{i}$ the $m$ data sets are independent, so we have: $$P(D_{1}D_{2}\dots,D_{m}|H_{i})=\prod_{j=1}^{m}P(D_{j}|H_{i}) \;\;\;\; (i=1,\dots,n)\;\;\;\text{Condition 1}$$ Now it would be quite convenient if the denominator also factored in this situation, so that we have: $$P(D_{1}D_{2}\dots,D_{m}|\overline{H}_{i})=\prod_{j=1}^{m}P(D_{j}|\overline{H}_{i}) \;\;\;\; (i=1,\dots,n)\;\;\;\text{Condition 2}$$ For in this case the likelihood ratio will split into a product of smaller factors for each data set, so that we have: $$LR(H_i)=\prod_{j=1}^{m}\frac{P(D_{j}|H_{i})}{P(D_{j}|\overline{H}_{i})}$$ So in this case, each data set will "vote for $H_i$" or "vote against $H_i$" independently of any other data set. The exercise is to prove that if $n>2$ (more than two hypothesis), there is no such non-trivial way in which this factoring can occur. That is, if you assume that condition 1 and condition 2 hold, then at most one of the factors: $$\frac{P(D_{1}|H_{i})}{P(D_{1}|\overline{H}_{i})}\frac{P(D_{2}|H_{i})}{P(D_{2}|\overline{H}_{i})}\dots\frac{P(D_{m}|H_{i})}{P(D_{m}|\overline{H}_{i})}$$ is different from 1, and thus only 1 data set will contribute to the likelihood ratio. I personally found this result quite fascinating, because it basically shows that multiple hypothesis testing is nothing but a series of binary hypothesis tests.


title: Nonparametric expected value estimation of sample from unknown distribution

text: I have a data sample (in this case an EEG data sample, but my question refers to any type of data samples of prior unknown distributions). I would like to do a nonparametric estimate of the expected value for my sample. I did some research, from what I understood I can do this using bootstrap sampling. I found a pdf here giving a formula for bootstrap expected value, hopefully it's correct. In case it's not, can someone please let me know how to do it once I have generated the samples by bootstrapping? Another possibility seems to be MCMC, but I would need to know the distribution from what I understood. I could do a kernel density estimation probably, but I think using bootstrapping might be less complex? I can use python, Matlab or R, in case you do this kind of thing often and have code at hand to share, I'd really appreciate it. Any other methods/suggestions are more than welcome.


title: How do I create clusters with a completely categorical data?

text: I am working on the project that requires data mining. I have been asked to use R. I have a dataset with all categorical variables and would like to form clusters on that. I am unable to figure out how to do so in R. Here is what I have done: I have converted all the variables to "factor" data type in R. But I am not able to see the underlying numbered levels. I also do not know how to use this with kmeans() to get the required result. My question is how do I form clusters on these factors. Here is what the data looks like: RowNum|EmpNum|EmpName|EmpOrganization|EmpTitle|EmpLeaderNumber|EmpDepartment|EmpAccesstoApplicaton|EmpAccessID The entire data is 14MB. The effort is to cluster people with similar access. So people with similar Title or in similar org might have similar access. I understand kmeans() isn't the best option, but that is what I would like to use for the first draft. I converted the EmpOrganization, EmpTitle etc to numeric data in excel using simple vlookup. It is easy to convert these to indicator variables using if statement in excel but I'm hoping that there is a more efficient way to do this in R itself.


title: How to compute this conditional probability in Bayesian Networks?

text: I met a problem related to conditional probability from the article "Bayesian Networks without Tears"(download) on page 3. According to the Figure 2, the author says $$P(fo=yes|lo=true, hb=false)=0.5$$ While I don't know how to calculate it and obtain the correct result. I attached the screenshot of Figure 2. from the article. I tried it as following $$P(fo|lo, \bar{hb})=\frac{P(lo, \bar{hb}|fo)\cdot P(fo)}{P(lo, \bar{hb})}=\frac{P(lo|fo)\cdot P(\bar{hb}|fo)\cdot P(fo)}{P(lo, \bar{hb})}$$ And $P(lo|fo)=0.6,\ P(fo)=0.15$, but how to know other 3 terms? Thank you! **P.S. According to some reference it has $$P(\mathbf{X})=\prod_{i=1}^{n}P(X_i|parents(X_i))$$ Thus $$P(fo|lo, \bar{hb})=\frac{P(lo, \bar{hb}, fo)}{P(lo, \bar{hb})}=\frac{P(fo)\cdot P(lo|fo)\cdot P(\bar{hb}|do)}{\sum_{{fo}^{'}}P(lo, \bar{hb}, {fo}^{'})}$$ $$=\frac{.15\times.6\times.3}{.15\times.6\times.3+.85\times.05\times.3}=0.6792453$$ While it is not equal to the author's result (i.e. 0.5)? What is wrong in my solution?** ![enter image description here](http://i.stack.imgur.com/i6wag.png)


title: Toeplitz variance matrix with nlme

text: I would like to specify a Toeplitz variance matrix for the random effects of my `nlme` model in `R`. Is it possible ? More precisely, my model looks like               lme(y ~ dose, data = dat, random = list(Lot = pdSymm(~ 0+dose )))      but I want a Toeplitz matrix instead of the unstructured matrix which is here specified by `pdSymm`. There's no available pdToeplitz-like argument in `nlme`, but maybe there is another way to do this ? I don't know if "Topelitz" is a standard terminology, this is the `SAS` terminology for a variance matrix with the only restriction that the diagonal entries are all equal. (EDIT) Sorry, my above definition of Toeplitz is erroneous. A Toeplitz matrix $\Sigma=(m_{ij})$ is the case when $m_{ij}$ depends on $i$ and $j$ only through $|i-j|$. This is not what I'm looking for, I really want a variance matrix with the only restriction that the diagonal entries are all equal.


title: What to conclude when you fail to find an association in an epidemiological study?

text: Normally when somebody finds an association in an epidemiological study people are quick to point out that it doesn't prove causality, that there are problems of missing co-founders, that it is at best hypothesis generating and at worst spurious. This leads to people not putting much weight on associations found in epidemiological studies. What if it goes the other way around? Say I already have a theory, maybe with some small earlier studies to back it up and even a good theoretical explanation for the effect. Then I do a big, well powered, epidemiological study and fail to find the association the theory predicts. How much weight can I put on the result now? Intuitively it seems to me that the result would be quite damning, despite it being only a epidemiological study. But I have long since learned not to trust my intuition when it comes to statistics. **Do all the weaknesses of associations found in epidemiological studies also apply to when you fail to find an association?**


title: Model for circular statistics

text: I am looking for advice on circular statistics. In particular, I'd like to know if any one had any advice/ references that deal with regression models for circular variables and whether it is possible to include random effects as well. At the moment I can fit very simple models in WinBUGS using wrapped cauchy distributions, but I don't know how to go on to the next step and add either fixed or random effects. Below is the WinBUGS code I have been using so far. I have tested it with simulated data and so far it has performed well, but trying to add in fixed/ random effects has, so far, not worked.               model{      for (t in 1:N) {           # likelihood for angles. We use the “ones” trick to sample from the      # Wrapped Cauchy distribution (see WinBUGS manual)           ones[t]<- 1      ones[t] ~ dbern(wc[t])      wc[t] <- (1/(2*Pi)*(1-rho[t]*rho[t])/(1+rho[t]*rho[t]-2*rho[t]*cos(theta[t]-      mu.t[t])))/ 300 # Density function for Wrapped Cauchy distribution           rho[t] <- lambda.t    # mean cosine for the circular distribution      mu.t[t]<- nu.t# mean direction for turns      }                     ###### priors for mean direction of angles      nu.t ~ dunif(-3.14159265359, 3.14159265359)     lambda.t ~ dunif(0,1) # prior for mean cosine of circular distribution      Pi <- 3.14159265359 # define Pi      }          ## Simulated Data ##     list(theta=c(1.57086666107637,0.624281203067249,4.83586153543422,5.52517105399153,0.250167755691792,5.24413183188724,0.175711907822086,0.503670499719972,0.00587906094477884,0.290131613934322,0.759047889069672,0.57973291007534,3.03128168541491,0.497790655905849,6.24730873150114,2.61159637947433,6.19811892339656,2.21476872674273,0.163464826891718,5.79300356573004,5.65352466175931,-0.0100726021401003,0.00574503925995024,0.260777171784755,5.8545805891331,6.09628602098184,6.07018161953988,5.90921466125829,0.0387070377090986,5.96019978900552,0.270388591408335,0.539775794451919,6.16303548945592,5.54317029065067,1.09867887761604,0.546155012914554,5.73154203573232,6.04837644493341,.242217723020124,0.201937287826239,6.19111529531002,0.602897213838987,5.53590129760264,0.304328180646957,6.12364810518025,0.0781317192586082,2.12148311222615,5.41742779164167,0.109722984863423,0.546244633029087,1.72483899231817,5.81142848191977,5.77431670621736,5.94852063016486,1.21880980868771,0.761391412364464,6.13385885651117,2.3278212791841,-0.00886837423371834,0.0509442654103693,0.919346146608449,0.22243489212092,0.0605109486858312,6.26215798187548,3.35930515203348,4.49262316826849,0.393662386151002,0.408276217352091,5.48604197934124,1.2319358669625,0.290890698266516,0.0356807866706245,5.01603150661483,2.13110569190685,5.58637984768018,0.705401496640296,0.474940761772081,5.58728776070886,6.12311166642116,0.00848809261322299,3.35074107197193,5.82089972193407,0.0531213061461832,5.97904289602246,4.31610462188531,5.61206825679503,0.184081838885041,0.288450927211418,0.594322121025956,1.07062485671203,0.400068367390392,5.08834932305335,4.35542895067301,6.08614182924595,6.14530696852739,5.25070254271081,5.91716602109256,1.78589020077607,6.23955405139402,6.09356179129423), N =100)


title: randomForest chooses regression instead of classification

text: I'm using the randomForest package in R and using the iris data, the random forest generated is a classification but when I use a dataset with around 700 features (the features are each pixel in a 28x28 pixel image) and the label column is named `label`, the `randomForest` generated is regression. I'm using the following line:               rf <- randomForest(label ~ ., data=train)      How come regression is used instead of classification? The data is read in through `read.csv()`.


title: Plotting a logistic GAM model in R - why is the scale not 0-1?

text: I have fit a simple binary logistic GAM model in R and have used the plot() function to plot the results of this model. The outputted graph shows a fitted line and a confidence interval, but the scale is clearly not 0-1. Does anyone know what is being plotted? Ideally I would like to get a graph of the predicted probability of the outcome versus the continuous predictor. Does anyone know how to get that out? Thanks


title: How to combine standard deviations?

text: I have a series of repeated measurements of light intensity at different points on a grid. I have made a histogram for each point showing the relative difference between modelling and measurement. My histograms are fairly close to normal with a mean of zero. I am using the standard deviation of this distribution to quantify the width of this distribution for each point on the grid. Ideally I would now like to combine all of these standard deviations into one figure to describe the entire grid. Is there a better way of doing this than simply finding the mean of all standard deviations? I don't want to assume that each point is ultimately a sample of the same population distribution (although in practice if there are significant advantages to doing so then pragmatically its probbaly not far from true). There is an equation at the bottom of the wiki page on standard deviations but I have never seen it before and no references are given. My question may be a repeat of Determining true mean from noisy observations, but I wasn't sure if that was assuming that each point was a sample from the same population.


title: Estimating Alpha parameter in a stable distribution

text: I am trying to estimate the alpha parameter of a supposed $\alpha$-stable distributed set of data. I have tried from the Hill estimator to more advanced fitting method, but they are or too approximated or too slow for my power of calculation. So after a lot of thinking i have found this way. I know that in a $\alpha$-stable distribution we have: $$ \lim_{x\rightarrow +\infty}f(x,\alpha,\beta)\sim -\alpha \gamma^\alpha \frac{\Gamma(\alpha)}{\pi}sin(\frac{\pi \alpha}{2})(1+\beta)x^{-(\alpha+1)} $$ and $$ \lim_{x\rightarrow +\infty}P(X>x_0)\sim \gamma^\alpha \frac{\Gamma(\alpha)}{\pi}sin(\frac{\pi \alpha}{2})(1+\beta)x^{-\alpha} $$ so plotting $P/f$ we must have straight line at x>>1 such that $$ \lim_{x\rightarrow+\infty}\frac{P(X>x_0)}{f(x,\alpha,\beta)}\sim -\frac{x}{\alpha} $$ and indeed i found a straight line at the tail of every data sample. Now i have a question: -because i found a straight line in the tail of every data sample, is this a general property of distributions?


title: Heteroskedasticity and standard deviation

text: ![enter image description here](http://i.stack.imgur.com/JOzkC.jpg)I am looking at state-wide data (entire population) of a school's grade as a function of the school's poverty index. The data appears to me to be an unconditional heteroskedastic distribution. I am drawing a line of best fit (LBF) through the data using linear regression (I changed this from 2-degree poly). What I am attempting to do is look at how each school is doing compared to how it was ‘predicted’ to do. In past years, a different test was used and the result was that the data was clearly linear, and moderately heteroskedastic. I then used SD as a measure so that I could look across multiple years. The question I have is how can I find SD across a LBF through heteroskedastic data? I just clustered the data based on Chernick's reply. You can see what it now looks like. Very interesting.


suppose $X$ has an F distribution with $\nu_1$ and $\nu_2$ degrees of freedom . it is well known that as $\nu_2$ approaches infinity. how can show limiting distribution of $Y=\nu_1X$ is chi-square with $\nu_1$ degrees of freedom


title: Variable creation in time series regression

text: I have a coworker who wants to create his own independent variable to add to a time series regression model because he believes that his variable will encapsulate more information. Is this advisable? His variable is essentially a weighted average, but I have concerns about the variable's interpretation (it is kind of odd), its behavior, its distribution etc. Any thoughts on this? The math major part of me felt slightly uncomfortable with this idea!


title: How can the formula for the expectation of a log-normal random variable be dimensionally sound?

text: If $X\sim\mathcal{LN}({\mu,\sigma^2})$, then $\mathrm{E}[X]=e^{\mu+\sigma^2/2}$. My question is: what right do we have to add a mean and variance together? If $X$ has physical dimensions, then the expression $\mu + \sigma^2/2$ is incoherent. So what gives? Thanks in advance!


title: Number of dependent variables

text: I have data from raters who each provided two ratings for a large set of categories. They were asked to rate how coherent each category was on a scale of 0 (members totally dissimilar) to 9 (members highly similar) and they did this twice: first, deciding if the category members had similar sizes, and, second, whether the category members has similar shapes (or something quite like this anyway - the two rating questions assessed different but related things). I averaged across the raters so now I have two ratings for each of my categories (size and shape) and there are various factors for the categories. I want to know whether category membership is best defined in terms of size, shape or a combination of both ... so I want to directly compare the size and shape ratings. This seems like a daft question but it's not clear to me how many dependent measures I have. I want to use the means of the raters in a repeated-measures ANOVA with rating as the dependent and instructions for rating instructions (size vs shape) as a repeated-measure factor. This allows me to directly compare the ratings. However, it could well be argued I have two dependent variables so I should be doing two separate ANOVAs (or a MANOVA). It seems like there could be a continuum from a situation which is definitely repeated measures using the same dependent variable (if I'd asked the raters to answer exactly the same question about shape similarity using exactly the same scale twice, 24 hours apart) to it definitely being two dependent variables (if I'd asked raters about the shape of the members and then about how expensive they were). My case seems to fall somewhere between. Correlation between the two types of ratings was low (-.16).


title: Struggling with non-normality in generalized linear model

text: Dear statistics experts, I am looking for correlations between certain measures of brain structural integrity (fractional anisotropy, given as ratio between two hemisphere ==> rational data range 0-1, normally distrubuted) and behavioral parameters from stroke patients which are not normal but in part extremely leftskewed. Having read numerous posts on assumptions of GLZM which I'd like to apply to assess the data, I still struggle with the assumption of non-normality of Y. Normality assumption in linear regression ; http://www.talkstats.com/showthread.php/55824-Adjusting-non-parametric- correlation-for-covariates ; Assumptions of generalised linear model ; Assumptions for robust linear regression, specifically ANCOVA type models ; ANCOVA and its disturbing assumptions For instance, it reads ".... In normal circumstances, the question isn't 'are my errors (or conditional distributions) normal?' - they won't be, we don't even need to check. A more relevant question is 'how badly does the degree of non-normality that's present impact my inferences?" I suggest a kernel density estimate or normal QQplot (plot of residuals vs normal scores). If the distribution looks reasonably normal, you have little to worry about. In fact, even when it's clearly non-normal it still may not matter very much, depending on what you want to do (normal prediction intervals really will rely on normality, for example, but many other things will tend to work at large sample sizes)...." I log transformed my dependent variable, can I use GLM normal distribution with LOG link function? Moreover, here it says: ".... Strictly speaking, neither, though the second is the thing you check. What is assumed normal is either the unobservable errors, or equivalently the conditional distribution of Y at each combination of predictors. The unconditional distribution of Y is not assumed to be normal. – Glen_b May 30 '13 at 13:38 " Normality of dependent variable = normality of residuals? or "....Since the residuals are just the y values minus the estimated mean (standardized residuals are also divided by an estimate of the standard error) then if the y values are normally distributed then the residuals are as well and the other way around. So when we talk about theory or assumptions it does not matter which we talk about because one implies the other." Normality of dependent variable = normality of residuals? So my question is: Is a GLZM applied to my data valid if the residuals look normally distributed but the Y is not? Would a non-significant Shapiro-Wilk test on the residuals together with a reasonable QQ plots allow to trust the confidence intervals estimated incl. the p-vals? I would be very grateful if somebody could help and clarify how to face this problem of non-normality in a generalized linear model. Having talked to two colleagues more familiar with modelling I got two different perspectives, A saying that normality of Y is crucial, B saying look at the residuals, GZLM should be robust to non-normal distributions. Thank you very much in advance. Kind regards, Robert Now I am confused with discrepant views on assumptions of GZLM. "


title: Moving average for continuous functions

text: In statistics, the moving average is usually defined for discrete data sets. Is there a moving average concept for continuous randomly fast-oscillating functions? I am seeking for the moving average determined in terms of integrals rather than sums.


title: Cross-Validation: Include IV which are significant in one model but insignificant in the other?

text: I ran regression analyses using SAS with two different data sets containing _different individuals_ but _exactly the same IV and DV_ : let's call them "low_deviance" and "high_deviance" and I would like to assess, whether there is a significant difference in the predictive power of the IV. Probably there is no such difference, but I just want to try out several methods checking this intuition... **What I have already tried:**   * _Beta Coefficients:_ 95 % confidence intervals show no significant difference between the beta coefficients of the two models/samples.   * _R²:_ Using the computer program “R2” (http://www.statpower.net/Software.html#R2) in a DOSBox I calculated 95 % confidence intervals for the determination coefficients of the two models/samples, too. These were overlapping again revealing no significant difference between the two models. **What I am trying to do now:**   * _Cross-Validation:_ Following Tabachnick & Fidell (2013: 140 ff.) I want to use the regression equation resulting from the data set “low_deviance” to calculate estimated values for the DV with the individuals from the data set “high_deviance”. Then I correlate the estimated values of the DV with the real values of the DV (1) for the data set “low_deviance” (estimated with the “correctly corresponding” equation gained from this data set) and (2) for the data set “high_deviance” (estimated with the regression equation gained from the “wrong” data set “low_deviance”). The higher the first correlation (data set used: “low_deviance”; regression equation used: “low_deviance”) the better is the predicive power of the model because the estimated values correspond well with the real values. The higher the second correlation (data set used: “high_deviance”; regression equation used: “low_deviance”) the better the regression equation of the data set “low_deviance” also predicts the real values of the data set “high_deviance”, i. e. the more the regression equation of the data set “low_deviance” is generalizable to the data set “high_deviance”. A low correlation would indicate that the regression equation from the data set “low_deviance” does not adequately predict the real values of the data set “high_deviance” which means that the groups of individuals (data sets, samples) differ in some relevant respect... And I could repeat the procedure vice versa. Now here is the question: **What do I do with IV which are significant in one model but insignificant in the other?** Do I always use all IV ignoring their significance in the different models/samples? Do I only use IV which yield significant beta coefficients in both models/samples? Alternative ideas? I will gladly clarify my question if necessary. Thanks a lot in advance!


title: How to solve the attached question?

text: Any help for solving the attached question (b), (c), (d) and (e) is appreciated.![enter image description here](http://i.stack.imgur.com/O2kCE.png)


title: What are some interesting examples of wrong or crazy inferences being drawn from Big Data?

text: I'm interested in well-known examples of Big Data misinterpreted, poorly analysed, or wrongly employed to unscientific and incorrect ends. Would appreciate any examples or observations. Thanks a lot!


title: Clustering procedure where each cluster has an equal number of points?

text: I have some points $X=\\{x_1,...,x_n\\}$ in $R^p$, and I want to cluster the points so that:   1. Each cluster contains an equal number of elements of $X$. (Assume that the number of clusters divides $n$.)   2. Each cluster is "spatially cohesive" in some sense, like the clusters from $k$-means. It's easy to think of a lot of clustering procedures that satisfy one or the other of these, but does anyone know of a way to get both at once? Thanks in advance for any pointers!


title: How to handle systematically missing values?

text: In my situation, one of two sources is not invoked if the confidence reported by the first source is higher than a threshold and hence it is missing in some examples. How can account for such missing data?


title: Gaussian MRF/Markov Network: the zero precision = no connection?

text: Gaussian MRF in Gaussian information form:   * edge potential: $exp(\frac{-1}{2} y_s\Lambda_{st} y_t)$    * node potential: $exp(\frac{-1}{2} y_t\Lambda_{t} y_t+\eta_ty_t)$ **Why:** > precision parameter $\Lambda_{st}=0$ <=> $x_s \bot x_t | rests$. according to books   * Koller's PGM: p.255    * Murphy's MLAPP :eq (19.28). * * * **My reasoning:** even $\Lambda_st=0$ ,the edge potential of st is > 0, that means there is connection in between, isn't it? But they say the connection is absent or $x_s \bot x_t | rests$ .


title: Fast/parallel alternative to GLS with nlme?

text: I am using the gls function from nlme to fit a fixed-effects model yet correct for spatial autocorrelation. My dataset has about 100,000 unique geographic observations, and running the following functions takes days on a R 64-bit edition installed on a quad-core Intel i7 3.4GHz CPU:               fit <- gls(frm, data, correlation = corSpher(form=~longitude+latitude,nugget=TRUE))       I tried the same on a Linux machine with 16 AMD Opteron processors, and it takes even longer. Is there any R package out there that will do the same job in a more efficient or parallel way?


title: Matching density with normal and uniform distributions

text: I will have data not too dissimilar from this: ![enter image description here](http://i.stack.imgur.com/E9sas.png) I'd like to match the density with a set of (truncated) normal distributions and one underlying uniform distribution. I know where the spikes are, these are my means for the normals, but I don't know what the variances are. Sort of like this decomposition of a similar dataset. ![enter image description here](http://i.stack.imgur.com/QMu04.png) So my only unknowns are variances of $k$ normals ($k$ is known) and the density of the uniform distribution. As for the truncation: Two of the means are 0 and 100 (boundaries), so I will only use half of the bell curve in each case. Edit: I've just realized that not all the normals will represent the same amount of data, so some truncation and rescaling will be in order. I suppose I'll have to do a grid search and compare the deviations from the density in each iteration -- but I don't really know where to start. Or if there's a package that would do that for me. I work in Stata and R, so any tips on on how to do it in either of the two is fine (preferably for the former). I have stumbled upon the `fmm` package (finite mixture models), but I'm not sure that fits (excuse the pun).


title: Hessian of Laplace distribution

text: The density of the Laplace distribution is given by: $$f(x;\mu,\sigma)=\frac{1}{2\sigma}\exp\left(-\frac{\vert x- \mu\vert}{\sigma}\right).$$ It is easy to see that this function is not differentiable at $\mu$. However, I am interested on some asymptotic normality results of the MLE $(\hat{\mu},\hat{\sigma})$, which require double-differentiability of the log likelihood function: $${\mathcal l}(\mu,\sigma) = \sum_{j=1}^n \log f(x_j;\mu,\sigma),$$ with respect to $(\mu,\sigma)$, for a random sample $(x_1,\dots,x_n)$. What I basically need is the existence of the Hessian matrix of the log-likelihood evaluated at the MLE, which means the existence of: $\frac{\partial^2}{\partial \mu^2}{\mathcal l}(\mu,\sigma) \Big\vert_{\mu=\hat{\mu}}$, $\frac{\partial^2}{\partial \sigma^2}{\mathcal l}(\mu,\sigma) \Big\vert_{\sigma=\hat{\sigma}}$, $\frac{\partial^2}{\partial \mu\partial\sigma}{\mathcal l}(\mu,\sigma) \Big\vert_{\mu=\hat{\mu},\sigma=\hat{\sigma}}$. Is there any reference to justify differentiability of the log-likelihood, or it is not differentiable at the MLE?


title: Good practices when doing time series forecasting

text: I've been working for months on short-term load forecasting and the use of climate/weather data to improve the accuracy. I have a computer science background and for this reason I'm trying to not make big mistakes and unfair comparisons working with statistics tools like ARIMA models. I'd like to know your opinion about a couple of things:   1. I'm using both (S)ARIMA and (S)ARIMAX models to investigate the effect of weather data on forecasting, do you think it would be necessary to use also Exponential Smoothing methods?    2. Having a time series of 300 daily samples I'm starting from the first two weeks and I perform a 5 days-ahead forecast using models built with auto.arima R function (forecast package). Then, I add another sample to my data set and I calibrate again the models and I perform another 5 days forecast and so on until the end of the available data. Do you think this way to operate is correct?  Thank you for your suggestions, although the target of our work is an engineering journal article, I would like to do a work as rigorous as possible from a statistical point of view.


title: Parameter selection in multiple kernel learning

text: I just got to apply multiple kernel learning to my data recently. I have data from three sources, so I want to learn three RBF kernels for each data source. But the MKL algorithms so far I know assume that the kernel parameters and the cost are fixed. When I used SVM with single RBF kernel before, I need to make a grid search for two variables (gamma and C) and do cross validation as well. I wonder how people in the application field usually deal with the kernel parameters. Take an n+1 dimension exhaustive search for n kernels? Use heuristic? Or make it into a convex problem to solve? Is there any method or tool related? Thanks.


title: Caret and coefficients (glmnet)

text: I am interested in utilizing caret for making inferences on a particular data set... Is it possible to do the following:   1. produce coefficients of a glmnet model I trained in caret. I would like to use glmnet because of the inherent feature selection as I do not believe glm has it?   2. other than the ROC metric, is there another metric that I can utilize to asses fit of the model? Such as adjusted R^2? The purpose of this analysis is to derive some inference on the effects of particular variables, rather than for prediction. I just like the caret package because it's been easy to work with thus far using matrices.


title: Transformation of skewed data shows strong linear relationship

text: The dataset comprises bandwidth usage for each customer. There is also a hybrid metric based on the distance covered by each traffic flow and aggregated to obtain 'Bit-Miles' for each customer (sum of traffic $\times$ miles for each flow). Clearly, there is some (causal) dependence among the above features but owing to skewness I had to resort to transformations (log, square root, z-scores etc.) After transforming both bandwidth and bit-miles, I see a strong linear relationship. Does this imply that bit-miles is a redundant feature (which seems counter- intuitive IMHO had there been no transformation)? Can I somehow prove/disprove that 'Bit-Miles' is redundant? Here are the plots before and after the transformation. ![before](http://i.stack.imgur.com/2UaAu.png) ![after](http://i.stack.imgur.com/jVbOi.png) Here is a small view of the data as well.               Name            Traffic(bps)            Bit-Miles     Customer1       729797243234.54         416983889869721.00     Customer2       411886504711.92         43841920479614.30     Customer3       253240650503.96         269534485841579.00     Customer4       251982742984.49         158900272002478.00


title: Confidence intervals and bootstrapping stochastic processes

text: I am currently using a stochastic method for prediction that only reports my parameter of interest $\widehat{T}$ and does not report confidence intervals, though I would like them. I understand that nonparametric bootstrapping involves   1. Resampling   2. Running a new estimate on the resampled data. However, intuitively I am also interested in the variance of this stochastic method if the data itself does not change. Would it be appropriate to rerun estimates on the same data ( _without_ resampling) and considering that my bootstrapped estimate CI for $\hat{t}$? If not, what is the significance of the method's variance without resampled data? I know for sure that I want to report it somehow, but I don't know what to call it.


title: How to choose the optimal number of latent factors in non-negative matrix factorization?

text: Given a matrix $V^{m \times n}$, Non-negative Matrix Factorization (NMF) finds two non-negative matrices $W^{m \times k}$ and $H^{k \times n}$ (i.e. with all elements $\ge 0$) to represent the decomposed matrix as: $$V \approx WH,$$ for example by requiring that non-negative $W$ and $H$ minimize the reconstruction error $$||V-WH||^2.$$ Are there common practices to estimate the number $k$ in NMF? How could, for example, cross validation be used for that purpose?


title: Queueing theory, and rate of processing as number of processors change

text: We are doing some testing of a major queue in our software system. During the entire testing, we have a "helper" that ensures that our queue stays full. Our intention is to figure out how quickly we are able to process the entries in the queue. We can vary the number of processors that are processing the entries in the queue. As we increase the number of processors, we find that we don't get much faster. ![How Number of Processors Affects Total Processing Rate](http://i.stack.imgur.com/CWPCu.png) Each blue point represents a 1-minute sample of number of items processed per second. For instance, with 9 processors going, we had a 1-minute period where we achieved a processing rate of about 22 items/second. In general, there seems to be a slight concave-down shape to this. In other words, we actually get faster, hit a peak, and then get slower as we add processors. That's confirmed in this next graph. ![How Number of Processors Affects Single-Server Processing Rate](http://i.stack.imgur.com/D83R6.png) This second graph shows that as we add processors, the number of items each processor handles per minute decreases. Now, I have a **strong** suspicion that there's a software/hardware problem at work here. But I want to rule one thing out: are there any results or theorems from queuing theory that would explain this phenomenon?


title: How could I generate an "explanation" for each prediction in a classification ?

text: I have a classification problem. My classes are 0 and 1. The dataset is a bit big, the training is done on 7 million lines and 100 + variables so I choose to use scikit learn and the logistic regression method (with weighted class). This performs rather well. My AUC is around 95 % with which I am satisfied for the moment. My problem that I need to be able to "explain" why each line predicted as being in class 1 is. So for example, a user can ask why is a line in class 1 and I should be able to give an answer like "It is because it has an important value on variable 1 and 2 and is very small along variable 3". The calculous could be done at query time or not, I do not really care. Has anyone already done something like this ? Is there any statistical / machine learning method to do that ? I guess the answer would be easier with a tree method. I am also interested by ideas for this kind of models. Thanks


title: How to visualize large amount of particles?

text: Often as a result of several simulations which are themselves intensive a huge amount of information i.e., points/line/plane depending to the subject of investigation are available. Although there are extensive list of multivariate analysis and data-mining techniques to summarize the results however I am often amazed by a simple visualization giving the most of the information in one shot.   Being familiar with some techniques in visualization e.g., application of OpenGL etc I would like to ask: 1- Is there any **software/technique/framework** /etc providing ability to visualize several million particles (points/lines/...)?   2- What is the most favorite application in the community of statistics? I'm not fan of **R**! What else? **Note:**   Let consider 3D data to be visualized **interactively** , so rotation and selection (slicing) being available **realtime**.


title: Distribution of independent binomial variables conditional upon the sum

text: Suppose that we have independent binomial variates with differing sizes and probabilities $X_i \sim Binomial(n_i,p_i)$, and $Z = \sum_iX_i$ is the sum. I understand that $Z$ is distributed poisson-binomial, and approximations to this distribution have been discussed here before. My question is what exact and approximate methods exist for estimating and generating variates from the distribution of $X$ given $Z$ (i.e. $P(X=x | Z=z)$ )? ## Update: So I've made some progress here. Regarding sampling from the distribution, Chen (2000) notes the connection between this model and sampling with probability proportional to size via maximum entropy. Here is a bit of R code:               library(sampling)     # generate via rejection     p <- c(.1,.05,.2)     n <- c(10,50,30)     z <- 10     a <- rbinom(1000000,size=n[1],p=p[1])     b <- rbinom(1000000,size=n[2],p=p[2])     c <- rbinom(1000000,size=n[3],p=p[3])     d <- cbind(a,b,c)     s <- a+b+c     d1 <- d[s==z,]     q <- p/(1-p)     pr <- q/sum(q)          #generate via PPS max-entropy     cn <-matrix(NA,nrow=10000,ncol=3)     for(ind in 1:10000){       q <- p/(1-p)       vals <- c(rep(q[1],n[1]),                 rep(q[2],n[2]),                 rep(q[3],n[3]))       i <- as.logical(UPMEsfromq(UPMEqfromw(vals,z)))       cn[ind,] <- c(sum(vals[i]==q[1]),                     sum(vals[i]==q[2]),                     sum(vals[i]==q[3]))          }     colMeans(cn)     colMeans(d1)          round(prop.table(table(cn[,2])),4)     round(prop.table(table(d1[,2])),4)      This seems to work, but I'm a little concerned that I am using $\pi_i = \frac{p_1}{1-p_i}$ for the inclusion probabilities, whereas in the paper we have $\pi_i = \frac{p_1 R(n-1,\\{i\\}^c)}{(1-p_i)R(n,S)}$. I have not tried to calculate the $R$ functions yet to see if there is a difference, as it is somewhat complicated from a computational perspective.


title: Comparison of transformed predictors

text: OK, I have been reading around on how to interpret transformed models, ... including this thread: Analysing log and square-root transformed variables I am interested in the OP's (of the link above) question in terms of using multiple transformations in the same model. Lets say I have proportion data on 3 land cover categories used as IVs in OLS regressions. The issue is that 2 of the 3 variables are skewed to some degree. One land cover variable is positively skewed; this is alleviated to some degree with the arcsin sqrt transformation, and the other is negatively skewed and this is alleviated by log10 transformation. Could someone explain to me 1) if this is OK, and 2) how to interpret the results as is done here: Interpretation of log transformed predictor I would like to be able to compare the relative influence of these predictors. An alternative would be to retain the log transform on one land cover variable, and not transform the other two. If this is done, how could I compare the relative effect of the land cover variables on a particular Y? Unfortunately, it seems what is good for the goose in this instance is not good for the gander. Thank you for any help. Patrick


title: Classic linear model - model selection

text: I have a classic linear model, with 5 possible regressors. They are uncorrelated with one another, and have quite low correlation with the response. I have arrived at a model where 3 of the regressors have significant coefficients for their t statistic (p<0.05). Adding either or both of the remaining 2 variables gives p values >0.05 for the t statistic, for the added variables. This leads me to believe the 3 variable model is "best". However, using the anova(a,b) command in R where a is the 3 variable model and b is the full model, the p value for the F statistic is < 0.05, which tells me to prefer the full model over the 3 variable model. How can I reconcile these apparent contradictions ? Thanks PS Edit: Some further background. This is homework so I won't post details, but we are not given details of what the regressors represent - they are just numbered 1 to 5. We are asked to "derive an appropriate model, giving justification".


title: Maximum Likelihood Curve/Model Fitting in Python

text: I have some 2d data that I believe is best fit by a sigmoid function. I can do the fitting with the following python code snippet.               from scipy.optimize import curve_fit     ydata = array([0.1,0.15,0.2,0.3,0.7,0.8,0.9, 0.9, 0.95])     xdata = array(range(0,len(ydata),1))          def sigmoid(x, x0, k):         y = 1 / (1+ np.exp(-k*(x-x0)))         return y          popt, pcov = curve_fit(sigmoid, xdata, ydata)      However, I'd like to use a maximum likelihood approach so I can report likelihoods. I think it's possible do to this using the statsmodels package, but I can't figure it out. Any help would be appreciated. ### Update: I think the approach may be to redifine the likelihood function as it is described here: http://statsmodels.sourceforge.net/devel/examples/generated/example_gmle.html The plot of the above code snippet looks like this: ![enter image description here](http://i.stack.imgur.com/21cbp.png) ### Update 2: Here's how to do it in R:               require(bbmle)          # this sigmoid function is used to make some fake data     rsigmoid <- function(y1,y2,xi,xmid,w){        y1+(y2-y1)/(1+exp((xmid-xi)/w))     }     counts <- round(rsigmoid(0, 1, 1:100+rnorm(100,0,3), 50, 10)*20,0)          # NOTE THAT THE SIGMOID FUNCTION IS REDEFINED AS AN R FORMULA     fit_sigmoid <- mle2(P1 ~ dbinom(prob=y1+(y2-y1)/(1+exp((xmid-xi)/w)), size=N),                   start = list(xmid=50, w=10),                   data=list(y1=0, y2=1, N=20, P1=counts, xi=1:100),                  method="L-BFGS-B", lower=c(xmid=1,w=1e-5), upper=c(xmid=100,w=100))


title: How to test significance of pre-post difference scores on 8 measures (1 group)

text: I have a set of pre and post treatment scores on 8 measures. I'd like to test which of these show significant improvement after treatment. If I use within- subjects t tests it will mean carrying out 8 t tests, and risking type I error. My textbooks all say that repeated subjects ANOVA is for more than two measures, but I only have 2 (pre & post). There is only one group - no control group or anything. What is the most appropriate test to use?


title: Polynomial Data Fitting for Two Unknown Equation

text: I have a kind of data and want to find the equation(poly coeff) of given data. For example equation for given sample data is simple `a^2*b+10`               a\b    5    10    15     ________________________     3|    55   100   145     4|    90   170   250     5|   135   260   385     6|   190   370   550      What is the correct algorithms?


title: Trouble using convert() in RTAQ package

text: I'm having a little trouble early on using the convert() function in RTAQ to convert .csv taq intraday trade data into an .RData format. I type this:               > from <- "2012-02-01"     > to <- "2012-02-29"          > ###convert data to .RData format     > convert(from, to, datasource = "/home/taylor/Desktop/trading",      datadestination = "/home/taylor/Desktop/trading", trades = T, quotes = F,      ticker = "AAPL", dir = F, extention = "csv", header = F)      and get this error:               [1] "no trades for stock AAPL"     There were 42 warnings (use warnings() to see them)     > warnings()     Warning messages:     1: In file(file, "rt") :     cannot open file '/home/taylor/Desktop/trading/2012-02-01      AAPL_trades.csv':          No such file or directory      or i type this:               > convert(from, to, datasource = "/home/taylor/Desktop/trading/aapl      2012-02.csv", datadestination = "/home/taylor/Desktop/trading", trades = T,     quotes = F, ticker = "AAPL", dir = F, extention = "csv", header = F)      and get this error               Error in setwd(datasource) : cannot change working directory      My hunch is that the raw data isn't in the anticipated format. In the documentation for RTAQ it mentions that the raw trade data has 9 columns (with no sample picture), however my data has 15 columns. An example row:               20120201,40000,793,AAPL,P,T...,100,457,,00,1,N,,AAPL,      The columns are Date, Time (HHMMSS), MS, Symbol, Exchange, Sale Condition, Trade Volume, Trade Price, Trade Stop Stock Indicator, Trade Correction Indicator, Trade Sequence Number, Source of Trade, Trade Reporting Facility, Symbol Root, and Symbol Suffix, respectively. If I have to change the format of the raw data, which columns do i drop/combine? Thanks in advance for any help.


title: How to identify quickly an aproximative number of clusters from a relatively small dataset

text: How is it possible to identify quickly (without doing many tests) an approximative number of clusters from a dataset which is not vary large, even if this value is not the correct number of clusters, I just want to identify a reasonable value representing the number of clusters from this small dataset. Note1: that I don't want to do many tests and/or cross-validate just to find an optimal number of clusters (time is important for me). Note2: I know that there is no way to automatically set the "right" K nor is there a definition of what "right" is. I just want a reasonably approximative value.


title: Priors for hierarchial Bayesian binomial model

text: I intend to use a hierarchical binomial model which would look like the following in BUGS language :               model{     for(i in 1:I){        y[i] ~ dbin(theta[k[i]],n[k[i]])     }     for(j in 1:J){        theta[j] ~ dprior(alpha, beta)     }     alpha ~ dhyperprior.alpha     beta ~ dhyperprior.beta     }      I wonder about the choice of the prior on $\theta$ and the hyperpriors. Do you know how to do this choice in order that: 1) There's no need to use MCMC (classical simulations suffice) and/or 2) This is a good "noninformative" choice, in the sense of the frequentist- matching property (the frequentist coverage of the posterior credible intervals is close to the credibility level)


title: Can you improve Adjusted R Squared in Multiple Regression with current data?

text: I am taking a theoretical approach to my multiple regression but the R squared and Adjusted R squared are both really low. I assume this means I can't improve my values with existing data and that the model is invalid? Can I do anything about it? I am trying to test the causes of the intensity of organised violence (DV) in Somali, Ethiopia between 1997-2007 against a set of independent variables of: vegetation cover; population increase; population density; presence of export ban; livelihood type; and land degradation. My unit of analysis is a spatial region. I have 200 observations. All but the vegetation cover are dichotomous dummy variables. I have created the dataset using satellite imagery, census data, and maps.


title: Are there guidelines for when to remove survey participants who blindly filled in the bubbles?

text: I just did a survey and I'm pretty sure some people just filled in the bubbles on the survey without reading the questions. A pretty obvious indicator are responses to items that are negatively coded (e.g. answering "strongly agree" to related items like "I'm a moron" and "I am not a moron"). Are there any methodological guidelines for finding and removing these responses? My results change quite a bit depending on my method for filtering people out. I've been toying with the idea of removing people whose standard deviation between all survey items is below a certain threshold. But, I'm unaware of validated guidelines for removing participants.


title: Iterated Conditional Mode approximation in E step of EM

text: I wanted to know what is the mathematical justification for using ICM as an approximation for the E step in an EM algorithm. As I understand in the E step the idea is to find a distribution that is equal to the posterior distribution of the latent variable, which guarantees that the likelihood increases or find the best possible distribution from some simpler family of distributions which guarantees that a lower bound of the likelihood functions increases. How does one mathematically justify the use of ICM in such an E-step? Any reference/derivations/notes would be very helpful. Thanks for your help.


title: Calculating joint probability and covariance given conditional probabilities

text: If we have the tables of conditionals $p(x|y)$ and $p(y|x)$, how can we calculate the joint probability $p(x, y)$? (Let us assume $x$ and $y$ are binary variables). Using these tables, how can we calculate the covariance of $x$ and $y$? Covariance of two variables are calculated using $E[xy]-E[x]E[y]$, but I am not sure how to evaluate this considering the probabilities. To clarify further, let me show the question as follows: ![enter image description here](http://i.stack.imgur.com/tz1xL.jpg)


title: Why does Okun's law regress change in GDP on change in unemployment rate?

text: Why does the Okun's law, the one that measures the relationship between unemployment rate and GDP, regress the percentage change in GDP on change in unemployment rate? Why cannot we just simply regress GDP on unemployment rate? What is the statistical difference between the two approaches?


title: Influential observations in Bayesian inference

text: In a (Bayesian in my case) inference, all observations do not contribute equally to the results. Some of them drive substantially the results while others have only negligible effects. When studying hierarchical models it is sometimes difficult to identify which data belong to the first class and to the second class. I guess that one possible empirical approach to identifying them could consist of inferring with a subset of the observations but I do not know how to do that properly. Moreover, I guess that some people have proposed less computationally costly approaches. Thanks for your expertise. Note: I am not comfortable in the choice of my title. Any suggestion to improve it?


title: normal distribution mean

text: Each week, John drives to his mother's house. The amount of time required for the trip varies and is normally distributed. On about 16% of trips, it takes him more than 54 minutes to reach his mother's house. On about 2.5% of trips, it takes him less than 33 minutes to reach his mother's house. Which is closest to the mean amount of time required to make the trip? Full disclosure: I'm studying for a test and this is a question from the practice test and I don't know where to start. I have the answer; I want to know how to get there.


title: Computing values of the Lilliefors c.d.f

text: How up-to-date is the assertion in this Wikipedia article that tables of the Lilliefors distribution have been computed only by Monte Carlo methods? (It would seem preferable to use deterministic numerical methods if possible.)


title: Normalizing data

text: i have a different sample size `n` for every month for example i have `13890` one month then `17756`, then `21425` the data every month for example `13890` is broken down into:               48 chairs, 12 tables, 2 couches etc...      next month we have similar metrics like               3 chairs, 23 tables, 4 couches etc..      i would like to know how i am doing with the furniture in relation to the total per month


title: What dependence is implied by a chi square test for independence?

text: I have a question on subject chi-square test for independence. I have, for example, two events A and B. If chi square test is not passed: is A dependent on B (A|B) or B on A (B|A)? Or does be valid both? (A|B and B|A). Thank you in advance.


title: Randomly distributed residuals or not?

text: Using Minitab I have checked for regression assumptions. Here I can see that the residuals (errors) of x4 are distributed normally. (y is response). However, are they distributed randomly? Looking at the "Versus Fits" graph they seem to be distributed normally, however I'm not very sure. What do you think? ![enter image description here](http://i.stack.imgur.com/godW0.png)               Regression Analysis: y versus x4           The regression equation is     y = 62,1 + 4,98 x4               Predictor    Coef  SE Coef     T      P    VIF     Constant   62,125    9,683  6,42  0,000     x4          4,979    2,634  1,89  0,063  1,000               S = 11,5115   R-Sq = 5,5%   R-Sq(adj) = 4,0%          PRESS = 8703,56   R-Sq(pred) = 0,00%               Analysis of Variance          Source          DF      SS     MS     F      P     Regression       1   473,5  473,5  3,57  0,063     Residual Error  61  8083,3  132,5       Lack of Fit    3   638,8  212,9  1,66  0,186       Pure Error    58  7444,6  128,4     Total           62  8556,9                1 rows with no replicates               Unusual Observations          Obs    x4      y    Fit  SE Fit  Residual  St Resid      35  3,00  49,00  77,06    2,21    -28,06     -2,48R      56  2,00  60,00  72,08    4,54    -12,08     -1,14 X          R denotes an observation with a large standardized residual.     X denotes an observation whose X value gives it large leverage.               Durbin-Watson statistic = 2,22949      Durbin-Watson is 2.22. Close to 2. So it means that the (standardized?) residuals are randomly and independently distributed?


title: How can I calculate the conditional probability of several events?

text: Could you inform me please, how can I calculate conditioned probability of several events? for example: P (A | B, C, D) - ? I know, that: P (A | B) = P (A intersection B) / P (B) But, unfortunately, I can't find any formula if an event A depends on several variables. Thanks in advance.


title: dispersion in summary.glm()

text: I conducted a glm.nb by               glm1<-glm.nb(x~factor(group))      with group being a categorial and x being a metrical variable. When I try to get the summary of the results, I get slightly different results, depending on if I use `summary()` or `summary.glm`. `summary(glm1)` gives me                   ...     Coefficients:                         Estimate Std. Error z value Pr(>|z|)           (Intercept)       0.1044     0.1519   0.687   0.4921           factor(gruppe)2   0.1580     0.2117   0.746   0.4555           factor(gruppe)3   0.3531     0.2085   1.693   0.0904 .         ---         Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1              (Dispersion parameter for Negative Binomial(0.7109) family taken to be 1)      whereas summary.glm(glm1) gives me                   ...     Coefficients:                         Estimate Std. Error t value Pr(>|t|)           (Intercept)       0.1044     0.1481   0.705   0.4817           factor(gruppe)2   0.1580     0.2065   0.765   0.4447           factor(gruppe)3   0.3531     0.2033   1.737   0.0835 .         ---         Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1              (Dispersion parameter for Negative Binomial(0.7109) family taken to be 0.9509067)      I understand the meaning of the dispersion parameter, but not of the line `(Dispersion parameter for Negative Binomial(0.7109) family taken to be 0.9509067)`. In the handbook it says, it would be the estimated dispersion, but it seems to be a bad estimate, as 0.95 is not close to 0.7109, or is the estimated dispersion something different than the estimated dispersion parameter? I guess, I have to set the dispersion in the `summary.nb(x, dispersion=)` to something, but I'm not sure, if I have to set the dispersion to 1 (which will yield the same result as `summary()` or if I should insert an estimate of the dispersion parameter, in this case leading to `summary.nb(glm1, dispersion=0.7109)` or something else? Or am I fine with just using the `summary(glm1)`? I would apprecciate your help!


title: How do I get the amplitude and phase for sine wave from lm() summary?

text: A simple sine curve could be written as $\text{amplitude}\cdot\sin(x+\text{phase})$. It can be also written in linear form as $a \cdot \sin(x) + b \cdot \cos(x)$. I run my analysis with R as:                fit.lm2 <- lm(temperature~sin(2*pi*Time/366) + cos(2*pi*Time/366))      summary(fit.lm2)          Coefficients:                            Estimate Std. Error t value Pr(>|t|)     (Intercept)             26.9188     0.1005  267.87  < 2e-16     sin(2 * pi * Time/366)   1.7468     0.1390   12.56  < 2e-16     cos(2 * pi * Time/366)   1.2077     0.1485    8.13 6.94e-11      The general form of the equation is $y = b_0 + b_1x_1 + b_2x_2$, thus, in my case, it can be written as $y = 26.9188x_0 + 1.7468x_1 + 1.2077x_2$. If I were to write it back to the simple sine form, $\text{amplitude}\cdot\sin(x+\text{phase})$, is it correct to say that: $\text{amplitude} = b_0 = 26.9188$ $\text{phase} = \arctan(\frac{b_1}{b_2})$ Is this the correct way how to do it? Thanks in advance.


title: Logistic regression prediction validation

text: I am try to evaluate a logistic regression with new data using R. I have created the logistic regression using the old data and evaluated that using chi^2 etc. i.e.               model <- glm(outcome ~ input + 0, data = training.data, family = binomial())      (note I do need the intercept = 0). And then the predicted probabilities:               predicted.probabilities <- predict(model, test.data, type = "response")      How can I now perform the same sort of tests with these new probabilities? and what are the best tests to do here?


title: Testing to see if a population is normal?

text: I'm an amateur statistician. I do it mostly for fun. I'm finishing up on my first semester stats class that covers up to hypothesis testing. I'm planning on beginning learning simple regression testing by myself so I can be prepared for the next level class when I have a free slot in my schedule again. Anyway, I have some data I would like to test to see if it's normal. Or, at least I can assume it to be normal with reasonable certainty. I have access to tons of samples but I cannot see the whole population. I understand that I can take a large enough sample from any type of population and treat the sampling distribution as normal by the central limit theorem but I'm actually trying to gain insight into the population itself. Reading around here I've realized that I need to do some hypothesis testing. From what I understand I need to set up the null hypothesis to be that the population is not normal and then try to disprove it. However, I don't know where to go from there. Can anyone suggest any reading I can do before I undertake this? Steps I should take? Thank you!


title: How to interpret glmnet?

text: I am trying to fit a multivariate linear regression model with approximately 60 predictor variables and 30 observations, so I am using the **glmnet** package for regularized regression because p>n. I have been going through documentation and other questions but I still can't interpret the results, here's a sample code (with 20 predictors and 10 observations to simplify): I create a matrix x with num rows = num observations and num cols = num predictors and a vector y which represents the response variable               > x=matrix(rnorm(10*20),10,20)     > y=rnorm(10)      I fit a glmnet model leaving alpha as default (= 1 for lasso penalty)               > fit1=glmnet(x,y)     > print(fit1)      I understand I get different predictions with decreasing values of lambda (i.e. penalty)               Call:  glmnet(x = x, y = y)                   Df    %Dev   Lambda       [1,]  0 0.00000 0.890700       [2,]  1 0.06159 0.850200       [3,]  1 0.11770 0.811500       [4,]  1 0.16880 0.774600        .        .        .       [96,] 10 0.99740 0.010730       [97,] 10 0.99760 0.010240       [98,] 10 0.99780 0.009775       [99,] 10 0.99800 0.009331      [100,] 10 0.99820 0.008907      Now I predict my Beta values choosing, for example, the smallest lambda value given from `glmnet`               > predict(fit1,type="coef", s = 0.008907)          21 x 1 sparse Matrix of class "dgCMatrix"                       1     (Intercept) -0.08872364     V1           0.23734885     V2          -0.35472137     V3          -0.08088463     V4           .              V5           .              V6           .              V7           0.31127123     V8           .              V9           .              V10          .              V11          0.10636867     V12          .              V13         -0.20328200     V14         -0.77717745     V15          .              V16         -0.25924281     V17          .              V18          .              V19         -0.57989929     V20         -0.22522859      If instead I choose lambda with               cv <- cv.glmnet(x,y)     model=glmnet(x,y,lambda=cv$lambda.min)      All of the variables would be (.). Doubts and questions:   1. I am not sure about how to choose lambda.   2. Should I use the non (.) variables to fit another model? In my case I would like to keep as much variables as possible.   3. How do I know the p-value, i.e. which variables significantly predict the response? I apologize for my poor statistical knowledge! And thank you for any help.


title: Meta analysis: comparing effects of two risk factors

text: I am conducting a systematic review to **compare the effect of intervention A vs. B on outcome Y (A, B, Y binary variables).** So far the collected studies are categorized as below:   1. 7 studies on both intervention A's effect on Y (A x Y) and B x Y    2. 60 studies on A x Y    3. 9 studies on B x Y  My questions are as follows:   1. In performing meta-analysis to compare A x Y and B x Y, is there any benefit (in terms of statistical strength of conclusion) to only pool the 7 studies that simultaneously study the effect of A & B on Y, in comparison to pooling 60 studies on A x Y and 9 studies on B x Y?    2. Of the 7 studies that simultaneously look at A x Y and B x Y, 4 studies found the effect of A x Y insignificant, but all found significant effect of B x Y. The effect sizes are similar. What are the risks of concluding that B has more significant effect on Y than A? I am aware of the risks associated with vote counting method, but in this case I am dealing with findings from identical set of studies.  I would very much appreciate any input or pointer to relevant literature/keywords. Thank you in advance.


title: How do you derive the Success-Run Theorem from the traditional form of Bayes Theorem?

text: In my industry it is common to test a sample of 20-30 and then use that data to draw conclusions about the reliability of the product with a certain confidence. We have tables for such things but it appears that for the case of 0 failures in the sample, the "Success Run Theorem" is used. In my references this appears as: $$R_c = (1-C)^{\frac{1}{(n+1)}}$$ where $C$ = confidence level, $R$ = reliability at confidence level $C$, and $n$ = sample size. However, I cannot find an explanation of how to get to the above equation from Bayes' theorem: ![enter image description here](http://i.stack.imgur.com/KH5f9.png) Every attempt to talk myself through Bayes theorem to arrive at the Success Run Theorem gets me confused. Even more confusing is when I try to extend my understanding to cases where some failures are observed in the sampling. Then I know to use this formula: ![enter image description here](http://i.stack.imgur.com/22Fvx.png) ![enter image description here](http://i.stack.imgur.com/9rogV.jpg) But again I don't understand where it comes from (binomial?) or how it relates to the above two other formulas, if at all. **My specific question would be how you go from Bayes Theorem (written as probabilities) to the Success Run Theorem (written as confidence, reliability, sample size)?** Thank you for helping a poor engineer lost in the world of stats.


title: R: Multivariate T distribution Sampling

text: In the question http://stackoverflow.com/questions/18153450/generating-random- variables-from-the-multivariate-t-distribution, I am confused by why the answer requires we modify the `Sigma` in such a way that we need to multiply is by `(D-2)/D`. Here sigma is the covariance matrix for me. The answer also mentions that the correlation matrix is defined when df > 2, shouldn't it be df>= 2? This is because the correlation coefficient can't be calculated when the data is continue itself one, there must be more than 1 series. Am I interpreting this correctly??


title: Machine-learning input data distribution

text: I'm trying to build a binary 1/0 ML classification algorithm, and was thinking about how to set up the input dataset. If the event I want to predict (the 1's) occur relatively less frequently in the total data than the 0's, does it makes sense to pare the dataset in such a way to get a more equal distribution of 1's and 0's? Would that be falsely representing the data to the algorithm? What are the disadvantages to doing so?


title: How do I test that two continuous variables are independent?

text: Suppose I have a sample $(X_n,Y_n), n=1..N$ from the joint distribution of $X$ and $Y$. How do I test the hypothesis that $X$ and $Y$ are **independent**? No assumption is made on the joint or marginal distribution laws of $X$ and $Y$ (least of all joint normality, since in that case independence is identical to correlation being $0$). No assumption is made on the nature of a possible relationship between $X$ and $Y$; it may be non-linear, so the variables are _uncorrelated_ ($r=0$) but _highly co-dependent_ ($I=H$). I can see two approaches:   1. Bin both variables and use Fisher's exact test or G-test.     * Pro: use well-established statistical tests     * Con: depends on binning   2. Estimate the _dependency_ of $X$ and $Y$: $\frac{I(X;Y)}{H(X,Y)}$ (this is $0$ for independent $X$ and $Y$ and $1$ when they completely determine each other).     * Pro: produces a number with a clear theoretical meaning     * Con: depends on the approximate entropy computation (i.e., binning again) Do these approaches make sense? What other methods people use?


title: Matching X'X with Wishart Samples in R

text: $X'X \sim Wishart(\Sigma,n)$, however I'm having a tough time producing this in R. Example:               data=cbind(rnorm(100,10,5),rnorm(100,5,2),rnorm(100,-4,3))          X=cbind(rnorm(1,10,5),rnorm(1,5,2),rnorm(1,-4,3))     t(X)%*%X     rWishart(10,99,cov(data))      The data generated from rWishart is not close to $X'X$. What am I doing wrong? The help documentation mentions $\Sigma$ should be a scaled matrix, however I'm unsure what this is.


title: Diagnosing collinearity in a Cox proportional hazards model

text: I am building a Cox Proportional Hazards Model to predict the survival outcome of seabird faced with predation pressure. I have 6 factor variables with two or three levels each that I have predicted to affect survival. Three of which are management relevant (they can be manipulated by wildlife managers to increase or decrease survival if significant). The ultimate goal of the model is prediction but I would like to include the management relevant variables as well even if not significant. How can I check for multicollinearity among my variables. I am using program R for the analysis.


title: Rolling quarters and four-quarters ended

text: I'm working on a financial project that in the past calculated its number quarterly. It would then also have a 4-quarters ended calculation as well. However, now we need to run the calculations on a rolling quarter basis, in other words every month we run the data for the previous 4 months. Does it make sense to do a 4-quarters ended on this data? For example, here is a set of current data ![enter image description here](http://i.stack.imgur.com/lXFNe.png) And here is a set of the new data and new type of reporting ![enter image description here](http://i.stack.imgur.com/En6Ot.png) My concern is that when we are reporting on distinct quarters, then it makes sense to do a 4-quarter ended calculation. But once we switch to rolling quarters, I don't think it makes sense to add up the data for the past 4 quarters since each has overlapping data.


title: How to show that method X is more reliable than method Y?

text: I compute a statistic s using two methods, X and Y. In order to measure the reliability of X and Y, I found split half reliability measures (alpha) using a bootstraping technique. So, for each of method X and method Y, I found the split half alpha value N times, where N is the number of bootstraps. Now, I have N alphas for X and N alphas for Y.               XAlphas = [xa1, xa2, xa3, ..., xan]     YAlphas = [ya1, ya2, ya3, ..., yan]      I'd like to show that X is more significantly more reliable than Y (meaning, the mean alpha for XAlphas should be higher than the mean alpha for YAlphas). What type of statistics test would I use for this? Is this a T Test?


title: Calculating prediction intervals when using cross validation

text: Are standard deviation estimates calculated via: $ s_N = \sqrt{\frac{1}{N} \sum_{i=1}^N (x_i - \overline{x})^2}. $ (http://en.wikipedia.org/wiki/Standard_deviation#Sample_standard_deviation) for prediction accuracies sampled from 10-fold cross validation? I'm concerned that the prediction accuracy calculated between each fold are dependent because of the substantial overlap between training sets (although the prediction sets are independent). Any resources that discuss this would be very helpful.


title: Combining two survey samples

text: I have two surveys of business owners. One is a sample (sample 1) of business owners who were _not_ members of the association, done using a random digit dialing approach. The other is a sample (sample 2) of business owners who were members of an association (used the list of members for the frame). Both samples were stratified by state. Sample 1 was weighted using state and gender. Sample 2 was weighted for non-response and weighted to reflect the association characteristics. As it turns out, there is a greater proportion of males in the association data (sample 2) than in the non-member data (sample 1 data). At first look, it appears that the association appeals to males more than females. We have 500 responses for each sample. I want to use a probit regression to find the factors (mostly demographic type of factors) influencing association membership. What is the best way to do this? I was hoping I could just combine the samples. However the samples were done independently. Clearly there is a huge difference in the size of the populations - sample 1 is huge (all owners in the country) and sample 2 is just the members. The survey instrument (questions) was almost identical - there were some extra questions for the members of the association (sample 2) related to the satisfaction with the association.


title: Part correlation and R squared

text: This is the result I've got after running linear regression analysis in SPSS: ![enter image description here](http://i.stack.imgur.com/GtViC.jpg) I am a bit confused why the sum of squared part correlations is not equal to (or less than) R squared, but rather exceeds the value of R squared (R squared = 77.7%, the sum of the squared part correlations = 88.6%). * * * **Edit** ![enter image description here](http://i.stack.imgur.com/DHcLx.jpg) The variance of each variable is represented by a circle of unit area. Total area of Y covered by the X1 and X2 areas represents the proportion of Y's variance accounted for by the two independent variables (areas B,C and D). B and D are portions of Y that overlap uniquely by X1 and X2. The area C is the overlap of both X1 and X2 with Y. The unique areas (B and D) are squared part correlation coefficient. The part correlation between Y and X1 from which the X2 have been removed - area B. R squared is the amount of variance in the dependent variable explained by all variables together. So, it sounds logical that the sum of the squared part correlation cannot exceed the value of R squared. Or I am wrong?


title: Concept of iterations in Adaboost

text: I can't seem to get my head around " **iterations** " in **Adaboost**. Are they analogous to **weak classifiers** that are used for **Boosting**? I've seen many examples of **Adaboost** where a programmers use a   1. **Single Classifier multiple Times** (initialized with different parameters resulting in different classifiers with different parameters).   2. **Multiple Weak Classifiers**. How does " **iterations** " fit in with above mentioned cases?


title: Strange linearity test for ANCOVA

text: There is an easy way to test the linearity hypothesis for a simple regression model $y \sim {\cal N}(\alpha+\beta x,\sigma^2)$: calling $H_0$ this model, perform a test against the ANOVA model $H_1$ obtained by considering the (numeric) covariate $x$ as a (qualitative) factor. Of course this requires to have several values of $y$ for each value of $x$. In R :               fit0 <- lm(y ~ x)     fit1 <- lm(y ~ factor(x))     anova(fit0, fit1)      This is quite sensible. But I have seen somewhere a linearity test for an ANCOVA model and I'm puzzling over what it means. The R code, denoting by $x$ the numeric covariate and $A$ the qualitative factor, is like:               fit0 <- lm(y ~ x + A + x:A)     fit1 <- lm(y ~ x + A + x:A + factor(x))     anova(fit0, fit1)      That means that $H_0$ is the classical ANCOVA model with interaction and $H_1$ is the bigger model obtained by adding the covariate $x$ considered as a factor. Does it make sense for you ? I cannot figure out how to interpret model $H_1$.


title: Wide or narrow confidence interval

text: I have a project in which I am looking at the correlation betwen fortnightly income and the amount saved each fortnight by low income workers. I am reporting the effect size and the confidence interval. One of my results is as follows: r (214) = 0.34, p < 001, 95% CI [0.22, 0.45] I know it is accepted that a narrow confidence interval is better than a wide confidence interval. However, I am unsure if this is a wide or a narrow confidence interval. Is there a quick way to determine what is 'wide' or 'narrow' in this context.


title: Is there a way in python/java/scala to convert/normalize log normal distribution into normal distribution?

text: We have a data set that looks like a lognormal distribution when we plot it. We would like to convert/normalize the distribution into normal distribution and see what feature weight got enhanced. It there a way to do it in python/java/scala?


title: One step ahead forecast with new data collected sequentially

text: Hi all I'm trying to do one step ahead forecast. Lets say I have 1000 data and fit an ARIMA model with it and then I do a forecast for one period ahead. When I get more data I would like to forecast another step using the new data without having to reestimate all coefficients and so on... This is my code but for some reason it's very slow for a bigger dataset and am not too sure that is doing what I want:               set.seed(1234)     y=ts(log(35+10*rnorm(1000)))     set.seed(4567)     new.data=ts(log(35+10*rnorm(10)))          library(forecast)     model = auto.arima(y)          onestep.for=forecast(model,h=1)     for (i in 1:10) {       data=c()       data=c(y,new.data[1:i])       newfit=Arima(data, model=model)       forec=forecast(newfit,h=1)       onestep.for=c(onestep.for,forec)     }


title: How to interpret results from non-parametric ANCOVA?

text: I 'm new to the CV and not very good at statistic:) I would much appreciate some help on a non parametric ANCOVA in R sm package. I do a pre post analysis on a set of pre/post variables of two groups (so `group` is a factor). The `pre`and `post` variables are numeric values (measures) or ratios of numeric measures. The `sm.ancova` output is simply a p-value and   * I do not really understand what the "models" are   * it does not say if there is an interaction factor:covariate. How should I interpret the output of the `sm.ancova` function in R? E.g. for `var1` to `var4`, `post` as response, `pre` as co-variate, `group` as factor:   * `var1` p-value: non significant for "equality model" and non significant for "parallel model",   * `var2` p-value: significant for "equality model" and non significant for "parallel model",   * `var3` p-value: non significant for "equality model" and significant for "parallel model",   * `var4` p-value: significant for "equality model" and significant for "parallel model". Edit: it seems that R sm package is about smoothing / "form of the relationship between y and a continuous covariate". Are there non-parametric alterniatives for ANCOVA when the distributional assumtions are not met? Some details: I compare 2 small gourps (group is a factor, categorical), **y** is post measure, and **x** is pre mesure (covariate). The question is if the effect of the treatement is the same on both groups.


title: How will two raters score videos?

text: I have 260 teacher conferences to score using a rubric. How do two raters need to score these for rater agreement? Do we need three raters?


title: What is the relationship between scale reliability measures (Cronbach's alpha etc.) and component/factor loadings?

text: Let's say I have a dataset with scores on a bunch of questionnaire items, which are theoretically comprised of a smaller number of scales, like in psychology research. I know a common approach here is to check the reliability of the scales using Cronbach's alpha or something similar, then aggregate the items in the scales to form scale scores and continue analysis from there. But there's also factor analysis, which can take all of your item scores as input and tell you which of them form consistent factors. You can get a sense of how strong these factors are by looking at loadings and communalities, and so on. To me this sounds like the same kind of thing, only much more in-depth. Even if all your scale reliabilities are good, an EFA might correct you on which items fit better into which scales, right? You're probably going to get cross loadings and it might make more sense to use derived factor scores than simple scale sums. If I want to use these scales for some later analysis (like regression or ANOVA), should I just aggregate the scales so long as their reliability holds up? Or is something like CFA (testing to see if the scales hold up as good factors, which seems to be measuring the same thing as 'reliability'). I've been taught about both approaches independently and so I really don't know how they relate, whether they can be used together or which one makes more sense for which context. Is there a decision tree for good research practice in this case? Something like:   * Run CFA according to predicted scale items      * If CFA shows good fit, calculate factor scores and use those for analysis     * If CFA shows poor fit, run EFA instead and take exploratory approach (or something) Are factor analysis and reliability testing indeed separate approaches to the same thing, or am I misunderstanding somewhere?


title: Testing equality of two or more contingency tables

text: Which options are there to test, whether two contingency tables of the same variables, say $X$ and $Y$, are equal across the levels of a third discrete variable, say $Z$? I can think of a log-linear model with terms $+X,+Y,+Z,+XY,+YZ,+XZ$. If it is significant against the saturated model, there is evidence of dependence of the contigency table on $Z$. Is this correct? What other tests are there available for this problem?


title: When Monte Carlo simulation can't be used to simulate a statistical system?

text: My question is simple. Which are the general conditions for which a Monte Carlo simulation can be used to represent a statistical system? Or conversely, which are the statistical system that cannot be represented, even in principle, by a Monte Carlo simulation?


title: Degrees of freedom for chi square test for particular model

text: Coming back to degrees of freedom for chi squared, I have found a wonderful answer here: How to understand degrees of freedom? However it is hard for me to apply it to my particular case. Let's say we observed counts of objects of 3 types:                A  B   C     35 121 344      The model, which has only one parameter, p, says that we should observe these objects with frequencies 0.01, 0.3, 0.69, which are some non-linear functions of p, let's say f(p), h(p), t(p). Given these frequencies, calculated using our model parameter p, we calculate expected counts and chi square statistic. Question: how many degrees of freedom for chi square distribution should we use in this case?


title: Ljung–Box test for a multivariate time series?

text: From Tsay's _Analysis of Financial Time Series_ , For a univariate weakly stationary time series $r_t$, its sample autocorrelation function $\hat{\rho}_l$ is defined as: ![enter image description here](http://i.stack.imgur.com/eU17a.png) and the Ljung-Box test is > ![enter image description here](http://i.stack.imgur.com/hkPZM.png) For a weakly stationary multivariate time series $r_t$, the sample lag-$l$ cross-covariance matrix $\hat{\Gamma}_l$ and sample cross-correlation matrix $\hat{\rho}_l$ of $r_t$ is defined as ![enter image description here](http://i.stack.imgur.com/fjSW3.png) ![enter image description here](http://i.stack.imgur.com/twLx3.png) The multivariate Ljung-Box test is: ![enter image description here](http://i.stack.imgur.com/hkEeJ.png) ![enter image description here](http://i.stack.imgur.com/GU9aF.png) I wonder how to see the test statistic for a multivariate time series is a generalization of the one for a univariate time series.   * Why do the coefficients $T(T+2)$ become $T^2$?   * Why does $\hat{\rho}_l^2$ become the trace of something?   * How is the representation in Kronecker product equivalent to (8.7)?


title: Neural network with and without cross validation

text: I have been running ANN (Neural network) on my data set, until last week that I figured out I will get more robust model using Cross validation, So That's why I have started using ANN with the aid of cross validation. For example using a 10-fold cross validation, all the dataset will be divided into 10 sunsets and each time one of the subsets is being used as test set while the rest is being used as training set. so 10 models in being built based on the 10-fold cross validation, and at the end we have all them combined in one model (?). Regarding all these process I was expecting a smaller R-sq (r-sq obtained from predicted values vs. actual values) using "ANN with cross validation" Vs. "ANN without Cross validation". but the strange thing is the Cross validation ANN R-sq are bigger than the simple ANN (without cross validation ANN) ones. I was thinking that maybe sth is going wrong?? I do not understand how, could you please help me with it?? Thanks guys