title: Inferring treatment start date for comparison group

text: I'm trying to work out the effect of treatment on a group using a counterfactual approach. However, I don't know how I can infer the treatment start date for the comparison group. For both groups I have a time point from which on they are eligible for treatment. I'm using kernel matching to match both groups on pre-treatment characteristics (i.e. from the time they became eligible). I don't want to use 1:1 matching because then I can't bootstrap the standard errors of the impact estimates. I want to measure outcomes relative to the treatment start. For that I will have to estimate a pseudo treatment start for the counterfactual. How can it be done? I thought of some sort of reverse matching (matching the treatment to the control as a weight on the treatment start date distribution), but not sure whether this could work.


title: A dataset with some very, very poor X variables - distribution and outliers

text: I have been handed a dataset with many poor variables and what appear to be outliers. Looking at histograms and box plots of the distributions, many cannot even form a box plot due to the distribution. Some of the distributions are just plain terrible with a chunk of single numbers and a second chunk of single numbers with an arbitrary gap. Nonreal example, variable x1, 500 entries =               50 entries with value 45,     200 entries with value 200,      250 variables with value 400.      Note that these may well be genuine. Others have dramatic outliers beyond an otherwise normal distribution. One problem is that I cannot legitimately exclude the outliers (if justified I would absolutely remove them) as I do not (and cannot) find their cause—but I feel I will have to remove most anyway. Also, I can only remove the outlier, not the whole row, otherwise my data set will have very, very few members—is this acceptable?


title: Where is density estimation useful?

text: After going through some slightly terse mathematics, I think I have a slight intuition of kernel density estimation. But I am also aware that estimating multivariate density for more than three variables might not be a good idea, in terms of the statistical properties of its estimators. So, in what sorts of situations should I want to estimate, say, bivariate density using non-parametric methods? Is it worth enough to start worrying about estimating it for more than two variables? If you can point to some useful links regarding application of estimation of multivariate density, that'd be great.


title: Hidden Markov Models methods for selecting optimal number of states

text: Package RHmm (R) I have a vector which I fit into a hmm model in an attempt to select an optimal number of states for a hidden markov model.               x<-c(-0.0961421466,-0.0375458485,0.0681121271,0.0259201028,0.0016780785,0.0311860542,           0.0067940299,0.0126520055,0.0357599812,0.0007679569,0.0409759326,0.0560839083,-0.0272581160,-0.0439501404,0.0321578353,0.0196158110,-0.0097262133,-0.0226182376,0.0119897380,-0.0099522863,-0.0359443106,-0.0039363349,-0.0476283592,-0.0383203835,-0.0518624079,0.0187455678,0.0950535435,0.0057115192,-0.0307805051,-0.0272725295,-0.0254645538,-0.0102565781,-0.0267986024,-0.0482906267,-0.0256826510,-0.0414746754,-0.0470666997,0.0284912760,0.1021992517,0.0875572274,0.0064152031,0.0200731787,-0.0091688456,-0.0575608699,-0.0442028942,-0.0277449185,-0.0115369429,0.0084710328,0.0745290085,0.0159369842,-0.0784550401,-0.0934970644,-0.0978390888,0.0160188869,0.0275268626,-0.0552651617,0.0033928140,0.0468507896,0.0374087653,0.0521167410,-0.0177752833,-0.0592673076,0.0514406681,0.0847486437,0.0738066194,-0.0098354049,-0.0572274292,0.0478305465,0.0096885221,-0.0445535022,-0.0153455265,-0.0105375508,0.0100704249,-0.0035215994,0.0243363762,0.0504443519,0.0570023276,0.0395103033,-0.0612817210,-0.0557737453,-0.0273657697,-0.0220077940,0.0083501817,0.0275081574,0.0323161331,0.0385741087,0.0175820844-0.0410599399,-0.0071019642,0.0431060115,-0.0107360128,-0.0007280372,0.0360799385,-0.0061620858  0.0164458899 -0.0050461344 -0.0578381588  0.0097198169  0.0027277926 -0.0127642317,     -0.0037062560, -0.0045482803,  0.0367596953, 0.0021176710,-0.0319243533,-0.0194663776,0.00 91915981,0.0061495737,-0.0090424506,0.0127655251,0.0161735008,0.0193814765,-0.0208605478,-0.0598025722,0.0022554035,0.0473633792,0.0247213549,-0.0063206694,-0.0201626938,0.0207952819,0.0379032576,0.0151612333,0.0038692090,0.0111271847,0.0497851603,0.0273431360,-0.0172488883,-0.0038909126,0.0264670631,-0.0065249612,-0.0467169856,-0.0255090099,0.0082489658, 0.0352569415,0.0272149172,0.0074228928,-0.0040191315,-0.0170611558,-0.0309531801,-0.0327952044,-0.0239372287,-0.0212792531,-0.0132712774,0.0086866983,-0.0007553260,0.0107026497,0.0065106253,-0.0321813990,-0.0081734233,0.0296845524,0.0268925281,-0.0025994962,-0.0038915206, -0.0126335449,0.0040244308,0.0227324065,0.0114903822,-0.0031516422,0.0031563335,0.0137143092,0.0026222849,0.0035802606,0.0111382363,-0.0008037881, -0.0282458124, 0.0056121633, 0.0254201390,0.0033781147,-0.0166139097,-0.0124559340,0.0088520417,0.0072600174, -0.0050320069,-0.0114740312,-0.0066160556, -0.0042080799, -0.0205501042,0.0027078715,  0.0122158472,-0.0206261771,-0.0267682015,-0.0107602258,0.0088477499,0.0165057256, 0.0106637013,0.0115216769,0.0278296526,0.0026376283,-0.0231543960,-0.0141964203)          #partitions test/train     nhs <- c(2,3,4) #number of possible states     S<-runif(length (x))<= .66     train<-print(S)          # mean conditional density of log probability of seeing the partial sequence of obs      for(i in 1:length(nhs)){     pred <- vector("list", length(x))         for(fold in 1:length(x)){             fit <- HMMFit(x [which(train==TRUE)],dis="NORMAL",nStates=nhs[i],             asymptCov=FALSE)             pred[[fold]] <-  forwardBackward(fit, x[which(train==FALSE)])        }     error[i] <- pred[[fold]]$LLH       }     nhs[which.max(error)]    # Optimal number of hidden states (method max log-likehood)      Every time I run the model trying to obtain the best number of states for the hidden markov model, I get a different number of states as I believe the model is trained over randmonly selected new values. This does not happen if I just fit the model.               #score proportional to probability that a sequence is generated by a given model     nhs <- c(2,3,4)     for(i in 1:length(nhs)){         fit <- HMMFit(x, dis="NORMAL", nStates= nhs[i], asymptCov=FALSE)         VitPath = viterbi(fit, x)        error[i] <- fit[[3]]     }     error<-c(error)      error[is.na(error)] <- 10000     nhs[which.min(error)]    # Optimal number of hidden states (method min AIC)      However the results of the two are very different. Which one is better? On one hand I have a model where I can test on new samples. On the other hand, the second provides best fit on seen samples however results are very different. In case of the model if I repeat the test given that the training/test set change (random) the resulting number of states also change. In this case what percentage sample / training should I use as to be certain that this choice will provide generalization in the number of states. What additional methods may I employ as to be able to select an optimal number of states Many thanks


title: Count data regression for data with upper bound

text: Is there someone who knows an econmetric model to describe count datas with an upper bound? Thank you!


title: Justification for using geometric weights in linear regression

text: In practical application, I have witnessed often the following practice. One observes a pair $(x_t, y_t)$ over time. Under the assumption that they are linearly related, we regress one against the other using geometric weights rather than uniform ones, i.e., the OLS minimizes $$\sum_{t=0}^\infty k^{t} (y_{T-t}- a x_{T-t}-b)^2$$ for some $k\in (0,1)$. This is very intuitive: we weight less observations far in the past. Compared to a "boxcar" weighting scheme, it has also the advantage of producing estimates that are changing smoothly over time, because observations do not fall abruptly off the observation window. However, I wonder if there is a _probabilistic_ model underlying the relationship between $x_t$ and $y_t$ that justifies this choice.


title: Using and/or operators in R

text: In R, what is the difference between:               if(x>2 & x<3) ...      and               if(x>2 && x<3) ...      Similarly:               if(x<2 | x>3) ...      and               if(x<2 || x>3) ...


title: How do I display detrended data from a linear regression?

text: I'm just having a look through Hadley's very excellent book about his ggplot2 R package. He has some code to remove a linear trend in the diamonds dataset, like so:               d <- subset(diamonds, carat < 2.5 & rbinom(nrow(diamonds), 1, 0.2) == 1)     d$lcarat <- log10(d$carat)     d$lprice <- log10(d$price)          detrend <- lm(lprice ~ lcarat, data = d)     d$lprice2 <- resid(detrend)          qplot(lcarat, lprice, data = d)     qplot(lcarat, lprice2, data = d)      Produces these graphs Unadjusted... ![enter image description here](http://i.stack.imgur.com/OV174.png) Detrended... ![enter image description here](http://i.stack.imgur.com/JWxiC.png) I'd like to see what the actual values of lprice would be without the effect of lcarat. Plotting residuals vs lcarat shows the right shape, but the points are shifted toward y = 0 (look at the range of the y-axis units). To get what I want, does it make sense to simply plot residuals + mean(lprice)? i.e. shift the previous graph up by mean(lprice).               qplot(lcarat, lprice2 + mean(lprice), data = d)      ![enter image description here](http://i.stack.imgur.com/hrrDr.png) **Does it make sense to do this? Is there a name for what I'm trying to do?**


title: Goodness-of-fit for contingency tables

text: I have a question on testing mortality table. Suppose I am given a simple mortality table:               age  | prob of dying  | prob of surviving     ---------------------------------------     20   |  0.01          |   0.99     21   |  0.02          |   0.98     22   |  0.03          |   0.97     23   |  0.04          |   0.96     ...      I want to test whether the table fit my observed data, i.e., whether the Observed values               age  | actual dead    | actual survivors     ---------------------------------------     20   | 0              |  397     21   | 1              |  189     22   | 0              |  136     23   | 2              |  100     ...      fit the expected values               age  | actaul dead       | actual survivors     ---------------------------------------     20   | 3.97 (397*0,01)   |  393.03     21   | 3.8   (190*0.02)  |  186.2     22   | 4.08  (136*0.03)  |  131.92     23   | 4.08 (102*0.04)   |  97.92     ...      How do I do that? What method should I use? Should I use Chi-Square, although my expected data is quite small? or there is another method?


title: Test whether data set approximates normal distribution using mean and median

text: I am of the knowledge that to test whether a data set approximates a normal distribution, the median and the mean should be approximately equal. So my question is to what degree should the difference between the median and the mean be accepted?


title: What to do about Seasonality Patterns in ACF, Time Series Data

text: I'm dealing with a time series data and I'm trying to construct a time series model for this particular dataset. I'm new to R and tried using the the `auto.arima` function under the forecast package:               fit <- auto.arima(tsdata, xreg=cbind(CSS2$Month,CSS2$DayID,CSS2$Year),                        stepwise=FALSE, approximation=FALSE)     summary(fit)     resid(fit)     acf(resid(fit))      However, I noticed some problems in the ACF plot and the PACF plot from the resulting model. Both of those appear to show some trends or seasonality. The data that I'm dealing with do have some strong seasonal patterns (on a weekly basis and on a monthly basis), and I thought that this would be captured using `auto.arima`. I do have the most recent version of the forecast package. I should mention that I have daily data for 3 years, so there are a total of 1095 observations. Any suggestions would be much appreciated, thank you! ACF Plot PACF Plot


title: Formula for computing the Pearson $\chi^2$, comparison with R

text: I suspect this question is really about basic definition, but I could not find the ressource I need to solve my problem. I want to understand why the pearson $\chi^2$ test statistic, and corresponding residuals, are computed the way they are in R. First, some tests:               >d<-data.frame(x=1:10000,y=sample(c(rep(1,100),0),10000,replace=T))     >M<-glm(y~x,family="binomial",data=d)     >d$p<-predict(M,type="response")     >chisq.test(table(d$p,d$y))          Pearson's Chi-squared test          data:  table(d$p, d$y)     X-squared = 10000, df = 9999, p-value = 0.4953          Warning message:       In chisq.test(table(d$p, d$y)) :       l'approximation du Chi-2 est peut-être incorrecte      Ok, now an alternative that gives consistent results computed as in this answer               >sum(residuals(M,type="pearson")^2)     [1] 10000.75      Considering the formula $$\chi^2 = \sum_{i \in \text{observations}} (O_i - P_i)^2/P_i$$ where the $O_i$ are the observed values and $P_i$ are the probabilities given by the model, I would have, perhaps naïvely, calculated               >sum((d$y-d$p)^2/d$p)     [1] 14      which provides another result. This is because the résiduals are different:               >head(residuals(M,type="pearson")^2)               1           2           3           4           5           6      0.001286902 0.001286924 0.001286946 0.001286968 0.001286989 0.001287011      >head((d$y-d$p)^2/d$p)     [1] 1.653989e-06 1.654045e-06 1.654101e-06 1.654157e-06 1.654213e-06 1.654269e-06      Going further, one discovers that the formula for the residuals is actually (see below for code) $$\chi^2 = \sum_{i \in \text{observations}} (O_i - P_i)^2/(P_i(1-P_i))$$ which is much nicer (to my opinion) as it results in a constant $\chi^2$ as a function of $P_i$. But where does that come from? Thanks! ----- how I found the last formula:               > getAnywhere(residuals.glm)     A single object matching ‘residuals.glm’ was found     It was found in the following places       package:stats       registered S3 method for residuals from namespace stats       namespace:stats     with value          function (object, type = c("deviance", "pearson", "working",         "response", "partial"), ...)     {         type <- match.arg(type)         y <- object$y             r <- object$residuals         mu <- object$fitted.values             wts <- object$prior.weights         switch(type, deviance = , pearson = , response = if (is.null(y)) {             mu.eta <- object$family$mu.eta             eta <- object$linear.predictors                 y <- mu + r * mu.eta(eta)             })             res <- switch(type, deviance = if (object$df.residual > 0) {             d.res <- sqrt(pmax((object$family$dev.resids)(y, mu,                 wts), 0))             ifelse(y > mu, d.res, -d.res)         } else rep.int(0, length(mu)), pearson = (y - mu) * sqrt(wts)/sqrt(object$family$variance(mu)),             working = r, response = y - mu, partial = r)         if (!is.null(object$na.action))                 res <- naresid(object$na.action, res)         if (type == "partial")             res <- res + predict(object, type = "terms")         res     }     <bytecode: 0x000000000a27cfc0>     <environment: namespace:stats>     > M$family$variance     function (mu)     mu * (1 - mu)     <bytecode: 0x000000000a191810>     <environment: 0x0000000008b31710>


title: Is it normal distributed?

text: My Error's plot from a linear regression ![enter image description here](http://i.stack.imgur.com/bprYB.png). Is it Normal distributted? If not, why? What that negative big bar means? The code used is:                 library(MASS)     sresid <- studres(reg3)      hist(sresid, freq=FALSE,           main="Distribution of Studentized Residuals")     xfit<-seq(min(sresid),max(sresid),length=40)      yfit<-dnorm(xfit)      lines(xfit, yfit)      Edited: ![enter image description here](http://i.stack.imgur.com/mWffB.png)


title: cox.zph says the standard deviation is zero

text: I need to compare my observed point vs random points that have been constrained to be relevant to that observation. In my case, I am looking at habitat selection along transects of animal movements. I have several types of variables categorical and continuous. Here is an example of 1 of my many models:               coxph(Surv(faketime,use)~bogB+CCB+estFB+mYRLWM+strata(LocNum)+cluster(mID),data=vector, robust = T)      bogB, CCB, estFB = binary vegetation values mYRLWM = continuous data strata(LocNum) = the "step" along my transects so it is an indicator that links the 1 observed with the corresponding random. Cluster(mID) = the animal id as all animals were followed for a while use = observed data faketime = should be the random points (i think, im confused on this one) I run the model               Call:     coxph(formula = Surv(faketime, use) ~ bogB + CCB + estFB + mYRLWM +          strata(LocNum) + cluster(mID), data = vector, robust = T)                         coef exp(coef) se(coef) robust se      z      p     bogB   -0.0142     0.986   0.0965    0.0909 -0.156 0.8800     CCB    -0.6356     0.530   0.2717    0.2002 -3.175 0.0015     estFB   0.1777     1.194   0.0771    0.0702  2.532 0.0110     mYRLWM  0.0967     1.101   0.2437    0.0489  1.977 0.0480          Likelihood ratio test=12.1  on 4 df, p=0.0165  n= 9966, number of events= 906           Call:     coxph(formula = Surv(faketime, use) ~ bogB + CCB + estFB + mYRLWM +          strata(LocNum) + cluster(mID), data = vector, robust = T)            n= 9966, number of events= 906                      coef exp(coef) se(coef) robust se      z Pr(>|z|)        bogB   -0.01417   0.98593  0.09652   0.09093 -0.156   0.8762        CCB    -0.63559   0.52963  0.27165   0.20016 -3.175   0.0015 **     estFB   0.17766   1.19442  0.07707   0.07018  2.532   0.0114 *      mYRLWM  0.09666   1.10148  0.24372   0.04890  1.977   0.0481 *      ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1                 exp(coef) exp(-coef) lower .95 upper .95     bogB      0.9859     1.0143    0.8250    1.1783     CCB       0.5296     1.8881    0.3578    0.7841     estFB     1.1944     0.8372    1.0409    1.3705     mYRLWM    1.1015     0.9079    1.0008    1.2123          Concordance= 0.531  (se = 0.062 )     Rsquare= 0.001   (max possible= 0.626 )     Likelihood ratio test= 12.11  on 4 df,   p=0.01654     Wald test            = 23.76  on 4 df,   p=8.934e-05     Score (logrank) test = 11.07  on 4 df,   p=0.02577,   Robust = 12.02  p=0.01719            (Note: the likelihood ratio and score tests assume independence of          observations within a cluster, the Wald and robust score tests do not).      But now I wanting to evaluate the model. Any better suggestions would be appreciated, but In reading, people have said use the cox.zph to look at the residuals vs predicted and when I do, I get this:               > cox.zph(mooseSSFB[[2]])            rho chisq   p     bogB    NA   NaN NaN     CCB     NA   NaN NaN     estFB   NA   NaN NaN     mYRLWM  NA   NaN NaN     GLOBAL  NA   NaN NaN     Warning message:     In cor(xx, r2) : the standard deviation is zero      I don't understand where to even start to figure out why this is happening. Before trying the cox.zph I simply extracted the predicted and residual values and plotted those and the graph looks awful and Ive no way to understand or interpret that, so I didn't know if I was even doing that right. How in the heck do you interpret that???               resB<-residuals(mooseSSFB[[2]])     predB<-predict(mooseSSFB[[2]])     plot(resB, predB)      ![residuals vs predicted of model 2](http://i.stack.imgur.com/f815y.png)


title: How do I perform a Mixed model analysis on my data in SPSS?

text: In my thesis I'm trying to discover which factors influence the CSR (corporate social responsibility, `GSE_RAW`) behavior of companies. Two groups of possible factors / variables have been identified: company-specific and country-specific. First, company-specific variables are (among others)   * `MKT_AVG_LN`: the marketvalue of the company   * `SIGN`: the number of CSR treaties the company has signed   * `INCID`: the number of reported CSR incidents the company has been involved in Second, each of the 4,000 companies in the dataset is headquartered in one of 35 countries. For each country, I have gathered some country-specific data, among others:   * `LAW_FAM`: the legal family the countries' legal system stems from (either French, English, Scandinavian, or German)   * `LAW_SR`: relative protection the countries' company law gives to shareholders (for instance, in case of company default)   * `LAW_LE`: the relative effectiveness of the countries' legal system (higher value means more effective, thus for instance less corrupted)   * `COM_CLA`: a measurement for the intensity of internal market competition   * `GCI_505`: mesurement for the quality of primary education   * `GCI_701`: measurement for the quality of secondary education   * `HOF_PDI`: power distance (higher value means more hierarchical society)   * `HOF_LTO`: country time orientation (higher means more long-term orientation)   * `DEP_AVG`: the countries' GDP per capita   * `CON_AVG`: the countries' average inflation over the 2008-2010 timeframe In order to make an analysis on this data, I "raised" the country-level data to the company-level. For instance, if Belgium has a `COM_CLA` value of 23, then all Belgian companies in the dataset have their `COM_CLA` value set to 23. The variable `LAW_FAM` is split up into 4 dummy variables (`LAW_FRA`, `LAW_SCA`, `LAW_ENG`, `LAW_GER`), giving each company a 1 for one of these dummies. This all results in a dataset like this:               COMPANY MKT_AVG_LN ... INCID ... LAW_FRA LAW_SCA ... LAW_SR LAW_LE COM_CLA ... etc     ----------------------------------------------------------------------------------        1      1.54          55          0       1          34     65     53        2      1.44          16          0       1          34     65     53        3      0.11           2          0       1          34     65     53        4      0.38          12          1       0          18     40     27        5      1.98         114          1       0          18     40     27        .       .             .          .       .           .      .      .        .       .             .          .       .           .      .      .      4,000    0.87           9          0       1           5     14     18      Here, companies 1 to 3 are from the same country A, and 4 and 5 from country B. First, I tried analyzing using OLS, but the model seemed very "unstable", as is shown below. The first model has a r-squared of .516: ![OLS 1](http://i.stack.imgur.com/9l9iB.png) Adding only two variables changes many of the beta's and significance levels, as well as the r-squared (.591). Of course the r-squared increases when variables are added, but this is quite an increase from .516: ![OLS 2](http://i.stack.imgur.com/3jy4O.png) Eventually, it was suggested in another post that I should not use OLS here but mixed models, because of the categorical countly-level data. However, I am confused as to how perform this in SPSS. The examples I found online are not comparable to mine, so I don't know what to fill in, amongst others, in the below mixed model dialogue: ![SPSS Mixed Model SPSS dialogue](http://i.stack.imgur.com/dnjjI.png) Could somebody using SPSS please help me explain how to perform this analysis so that I may come to a regression model (CSR = b1*MKT_AVG_LN + b2*SIGN + ... + b13*CON_AVG) so that I can conclude wheter CSR is determined by company- features or country-features (or by neither or both)? Any help is greatly appreciated!


title: Statistical analysis of multiple algorithms over a single dataset

text: I have a dataset $X = \{x_1, ..., x_n\}$. I also have three algorithms, $A_1, A_2, A_3$, that each take a single data point as input and produce some measure of how well they performed. If I apply each of the algorithms to each of the data, I get a matrix of performance measures, $M$, where $M_{a,i}$ is the value of $A_a(x_i)$. This matrix is therefore $3 \times n$. What statistical methods are most appropriate to analyse the differences in performance of these algorithms and the significance of any differences? If it matters, the "better" algorithm is that which gives performance measures closer to zero. I'm interested in determining whether any algorithm is _significantly_ better than another.


title: Not sure whether to treat my data as nominal or scale

text: The data I have is based on percentages. I have several variables on percent of children that have been immunized. So one column that is one kind of immunization, and is ordered by state. So for example in my data it shows that in Massachusetts, 53.3% of children have gotten the HPV4 vaccine. I want to correlate them with another column that has percent of children that have asthma by state, like in Massachusetts, 9.8% of kids have asthma. I'm not sure whether to treat the data as nominal or scale


title: Bayesian uninformative priors vs. frequentist null hypotheses: what's the relationship?

text: I came across this image in a blog post here. ![enter image description here](http://i.stack.imgur.com/n2EuU.jpg) I was disappointed that reading the statement did not illicit the same facial expression for me as it did for this guy. So, what is meant by the statement that the null hypothesis is how frequentists express an uninformative prior? Is it really true? * * * **Edit:** I'm hoping someone can offer a charitable interpretation that makes the statement true, even in some loose sense.


title: How to transform an accuracy distribution for a violin plot

text: I am trying to find the best way to visualize different distributions of accuracy. Accuracy here is a value in the interval [0,1], 0 meaning not accurate, 1 maximum accuracy. I have different methods to compare, so I decided to use violin plots. ![accuracy violin plots](http://i.stack.imgur.com/AaZaE.png) Distributions are clustered near 1 but they have also a long tail (the first one is cut at 0.45). How can I transform the data (e.g. in log-scale), in order to visualize better the differences between these plots? I want to focus in the interval [0.8, 1], but I want also to retain the long tails. I do not want to use boxplots because in this case I saw that they do not show correctly the distributions (also because the lower quartile is already 1). I add also the corresponding boxplot. ![boxplot](http://i.stack.imgur.com/cgkz8.png) Thanks


title: Cross Correlation (xcorr in MATLAB) vs Pearson Correlation (corr in MATLAB)

text: Are the values of `xcorr(x,y)` in MATLAB correlation values or not? I'm asking this because in MATLAB `xcorr(x,y,'coeff')` normalizes values. Is it normalizing covariance values to get correlations? I'm confused if cross correlation values are necessarily between -1 and 1 like Pearson correlation values. Also, I see from running an example that `xcorr(x,y,0,'coeff')` != `corr(x,y)`. Could someone explain this?


title: How to improve classification performance based on multiple known classification results

text: I am working on a classification problem, which may contains a unknown number of data classes, typically 5-50 classes in each sample. I had several classification algorithms, each gives me a classification output based on a given sample. However, these classification outputs are almost unlikely to agree with each other completely. When I looked into these classification results, I noticed that there always exists a better classification output in a subset of a given sample. This means that I shall be able to get a better classification output if I can use these existing classification outputs properly. However, the problem is I do not know how to do this job. So far, what I have is 500 manually annotated sample data, and corresponding classification results for these training data. Is there any simple way that I can train a system to automatically form a new classification output for me? Another thing that bothered me is that I am completely lost in formulating my task in mathematics. I know my following description is problematic, please help me to correct it. Suppose we have a training sample composed of $m$ 2D data points $\Omega =\{d_1,d_2,\cdots,d_m\}$ with a corresponding known label set $C_{\rm gt} = \left\{l(d_i|{\rm gt})\big|i\in\Omega\right\}$, where $l(d_i|{\rm{gt}})=k$ means the label of data point $d_i$ is of the $k$th class in ground truth (gt). Assume I have $n$ classification algorithms $A_1(\cdot),\cdots,A_n(\cdot)$ and each takes the entire dataset as input and generates a classification result $C_j = \left\{l(d_i|A_j)\big|i\in\Omega\right\}$, which is the label collection for all data points using algorithm $A_j$. And my objective is to find a function $A_{new}=f(A_1,\cdots,A_j)$ such that $$A_{new} = \arg\min_{f(\cdot)}\sum_{i\in\Omega}\left\|l(d_i|A_{new})-l(d_i|\rm{gt})\right\|.$$ This is some formulation that I learnt in signal estimation, but in this problem a label is quite different from a random variable: 1) it is discrete rather than continuous 2) even if one algorithm assigns a label $k'$ to a data point which is marked as class $k$ in the groundtruth, it does not necessarily mean it is wrong, because what we really care about is whether the $k$th class in the groundtruth has been marked as one single class in our output, and we do not care about which label it used in the output. Finally, I have no idea how to construct a function $f(\cdot)$.


title: How to classify different types of fingerprints using svm two class problem in matlab?

text: Here I am working on fingerprint recognition. I am using SVM for fingerprint feature verification and classification. I want to classify fingerprint images using two class problem of SVM that means in one class I want right loop, left loop and tented arch type images and in another class I want arch type and whorl type images and then I should get the result in which class testing fingerprint image is there, but I am not getting exact way for implementing it in MATLAB so please tell me how to use SVM for fingerprint classification in MATLAB.


title: Correlation between series of series

text: I want to study some data describing $k$ different systems. For a given system, these data can take two forms:   * _Global_ metric: a single value describing the whole system.   * _Local_ metric: a series of values describing each element in the system. This is _not_ a temporal series. Systems can have different sizes (contain different numbers of elements). There are several global and local metrics. My goal is to identify relationships between metrics. For _global_ ones, I simply use the correlation coefficient: each metric leads to a series of $k$ values, so it's pretty straightforward. However, which tool can I use to perform an equivalent operation on _local_ metrics? For a given system of size $n_i$, a local metric leads to $n_i$ values. For the whole dataset, I therefore have a series of $k$ vectors, each one with a (possibly) different size $n_i$. Hence my question. If some tool exists for this purpose, an additional question would then be to know how to compare a global metric with a local one. * * * Here's a small example:   * System 1: $n_1$=5      * Global metric 1: 4     * Global metric 2: 10     * Local metric 1: 4, 45, 12, 62, 47     * Local metric 2: 1, 6, 7, 9, 3   * System 2: $n_2$=3      * Global metric 1: 12     * Global metric 2: 2     * Local metric 1: 26, 23, 21     * Local metric 2: 9, 6, 7   * System 3: ... To compare global metrics 1 and 2, I simply process the correlation between series (4, 12, ...) and and (10, 2, ...). My questions are:   * How to compare local metrics 1 and 2: ((4,45,12,62,47), (26,23,21), ...) and ((1,6,7,9,3), (9,6,7), ...)   * How to compare global metric 1 and local metric 1: (4, 12, ...) and ((4,45,12,62,47), (26,23,21), ...) * * * The only method I could think of is the following:   1. Process the normalized histogram of each local metric, for each system.   2. For each local metric, consider each pair of systems in the dataset, and process the distance between their histograms (using Earth Mover's distance). Thus, for a given metric, I obtain a series of $k(k-1)/2$ values.   3. To compare two local metrics, process the correlation between their series of distances. The problem with this approach is its results are more difficult to interpret. If I find a strong correlation for two local metrics, it doesn't mean they are directly related. Instead, it means that if one metric is similarly distributed over the elements of both systems, then this also holds for the other metric (and vice-versa). Any feedback and advice welcome! (including those related to the title of this question)


title: Multiple bivariate normal distributions and conditional expectation

text: The problem I have is one that ideally seems to be suitable by considering it as a multivariate normal pdf problem however this way seems impossible. Consider three Normal distributions $X\sim N(\mu_{x},\sigma_{x})$, $ Y\sim N(\mu_{y},\sigma_{y})$, $Z\sim N(\mu_{z},\sigma_{z})$. Now consider the fact that we have the covariance (calculated from sample results) of $X$ and $Z$ and the covariance of $Y$ and $Z$ only. i.e. we cannot calculate any correlation coefficient between $X$ and $Y$. So I know the correlation $\rho_{xz}$ and $\rho_{yz}$ only. Now I have the condition $X = a$ and $Y = b$. How can I calculate the expected value of $Z$? For $X,Z$ we have $E(Z|X=a) = \rho_{xz} \frac{\sigma_{z}}{\sigma_{x}}(a -\mu_{x})+\mu_{z}$. For $Y,Z$ we have $E(Z|Y=b) = \rho_{xz} \frac{\sigma_{z}}{\sigma_{y}}(b -\mu_{y})+\mu_{z}$. Call the first expectation $e_{z|x}$, and the second $e_{z|y}$. If these two values are different (which seems probable) what can I do to estimate $e_{z|x,y}$ (i.e. expected value of Z given values for X and Y)? My initial thought is simply to weight the two expectations by the correlation coefficients. So $e_{z|x,y} = \frac{|\rho_{xz}|e_{z|x} +|\rho_{yz}|e_{z|y}}{|\rho_{xz}|+|\rho_{yz}|}$, to get some sort of a mean estimate from the two above estimates. This assumes the correlations are not zero and their sum is not zero. I am not in any way knowledgeable in this sort of thing so doubt I am right. Please let me know if there is a better, or even a correct, way to do this, remembering that We do not know the covariance of $X$ and $Y$, not even an estimate from samples. So I don't think we can consider this case a multivariate problem of 3 dimensions. Edit: I think Glen_b's answer seems pretty conclusive, but something is niggling me. If we can guess the correlation sign (not the value, just the sign) of X and Y (e.g. we guess $\rho_{xy}$ is positive),it seems the above method would have some merit. Is this just wishful thinking? Secondly, I'm starting to doubt my own doubts now, because thinking some more, if $\rho_{yz}$ is close to one and $\rho_{xz}$ is close to zero, or very small, then taking the weighted average seems inappropriate since the YZ correlation is far superior to the XZ, so weighting the average is in effect diluting our accuracy and our certainty of getting the right value. Is this thinking valid?


title: Predicting total number of bugs based on number of bugs revealed by each tester

text: Assuming n testers were independently testing the same application for a given period. Each tester found a given set of bugs (Some of the bugs were detected by more than one tester). For example: Tester 1 found bugs {1,2,3,4,5} Tester 2 found bugs {3,5,6,7} Tester 3 found bugs {1,3,5,8,9,10} (The serial number is assigned to each bug once it is detected by some tester, according to the detection order, nevertheless, the testers are unaware of each other, and of the assigned serial numbers.) Assuming all bugs have equal probability to be detected. (but without assuming that all testers are equally talented) Can I calculate the probability of having x undetected bugs? How?


title: Comparing prices between different stores (non-parametric models?)

text: I have thousands of product info in a database with category, title, description, make, model, price and store name. I have grouped those products by store names and it seem to span to 4 stores. So thousands of products in those 4 stores. I want to compare every store's product with rest of stores' products to see if a product in store A is present in Store B or Store C or Store D. After reading about Ranking, I think it's non-parametric models which I am going to need. But by reading Non-parametric Statistics I found a lot of methods there and I have never used any of them before. I have used very basic stats like mean, median, variance, regression etc. I am looking for your opinion about what non-parametric method will best estimate product comparisons.


title: Normality assumption and sample size

text: I know this is a very debated topic, even on this site, but I still couldn't find an answer to my problem. Recently I am working with large samples (300, 400 and more). For now, I am trying to use simple techniques, such as correlation, T-tests, and ANOVA, all of which require the normality assumption (from what I have read so far in the textbooks, online etc.) I also read that, if the sample size is large enough, the normality assumption is not so much of a problem and these techniques are robust to the violation of normality. Is normality a problem given my sample size? Should the data be at least bell-shaped, even if the tests fail to accept normality? Or could I get away even with extremely skewed or lumped data when using these techniques? Should I apply parametric techniques or should I just stick to the non- parametric ones, which, from what I know, have lower power? LATER EDIT for those who want to find more about my data: I have a variable which represents the number of days a user has been employed in the program (mean=176, median=167, stdev=87, IQR=113, Skewness=0.61, kurtosis=-1.64) sample=340users. The histogram does show that maybe the variable has the potential to split into 'fairly normal' groups (however I have not found this factor). This is the variable that for now i am trying to explain, in terms of 'The Number of Weeks The User has worked in the first 8 weeks', which takes values from 1 to 8, so I assume it is ordinal, and has a negative skew(most of them have worked 8 weeks out of 8). So, the main question would be is there a relationship between how much a user will stay employed and the amount of work he does in the first 8 weeks? Also, later I would compare the length of employment with other possible factors, for which currently I don't have the data (education, gender, age, ..)and try to build a more 'elaborate' model, but right now I try to analyse what I have with 'simple' statistics. For now, I have done a Spearman correlation between length of employment and weeks worked out of first 8, had a coefficient of 0.450 (which means low correlation), so I am trying to see if those who have worked an amount of no of weeks in the first 8 differentiate themselves in terms of length of employment from those who have worked fewer/more weeks. So I studied the distributions of length of employment for each of the groups (a group means a certain no of weeks they have worked in the first 8). In each group I have around 25-40 cases (except weeks worked=0 where I have 9cases and weeks worked=9 where there are 156). The normality test (Shapiro-Wilk, which I knew were suitable with such sample size) showed that only 2 groups out of 9 are normally distributed. So my thought was to drop ANOVA and t-tests and head for Mann-Whitney U Test and Kruskal-Wallis Test. However now I am re-considering maybe looking at ANOVA and t-tests because the groups do not look _that_ non- normal. Thanks.


title: One-way repeated measures anova

text: I tried finding an answer to this question on this and other sites but to no avail - if I am missing something please excuse my inability to locate the answer! Basically I have several dependent variables (`dv`) for two land management treatments (`treat`) measured quarterly (`season`). I know that I need to run a RM-ANOVA with season as a random effect, and so I am using the `lme` test in R. What I cant work out is whether I should be including season as a random variable like so:               lme <- lme(dv ~ treat, random = ~1|season, data=data)      or whether I should be nesting season within treatment like so:               lme <- lme(dv ~ treat, random = ~1|treat/season, data=data)      Many thanks for your assistance!


title: Multiple regression analysis with spatial data as independent variable

text: In my PhD thesis I am working on spatial modeling of different chemical parameters in groundwater, and for spatial modeling I am also using the multiple statistical approach. I have a question about multiple regression analysis. (Or it is better to use polynomial regression?) The equation for spatial regression modeling is: $Y = α + β_1x_1 + β_2x_2 +.... + β_ix_i + ε$ For my dependent variable, I have concentrations of calcium in groundwater, which were measured from different sampling points in the entire research area. For the independent variable, I choose the spatial data that influence the distribution of calcium in groundwater. I have lithology, vegetation, slope, climatic conditions (temperature, precipitation), depth of soil, ... The problem is that lithology and vegetation are categorical data (lithology = 3 categories from `1` to `3`, where `1` means clastic rocks, `2`= carbonate rocks and `3`= metamorphic and igneous rocks; and vegetation = 4 categories (`1`= bare rocks, `2`= agriculture land, `3`= grassland, `4`= forests); all others variables are numerical and continuous. Do you have any idea how to solve the problem with categorical data in multiple regression analysis? Might it be better to use some other method? Best regards and thank you very much for your help.


title: How do you calculate the statistical significance of a correlation between $Y$ and $X$?

text: Can one just fit a linear model between $Y$ and $X$, and then take the p-values of the $t$-test done on the slope (or beta) of the linear model? Since presumably correlation would be in the same direction as the slope/beta of the linear regression model?


title: How we can draw an ROC curve for decision trees?

text: Normally we cannot draw an ROC curve for the discrete classifiers like decision trees. Am I right? Is there any way to draw an ROC curve for Dtrees?


title: Which non-parametric test should I run on ordinal data from a small sample?

text: I have one group of six raters who scored recordings of two separate six- person groups of participants, each group under separate conditions. While the participants were divided into two independent groups, the ordinal data came from the same group of raters. Due to the small sample size (i.e. 12 ordinal scores to compare both participants), which test could I run. I've read that I'd need at least 20 scores to run a Wilcoxon... Thank you for your time and consideration in answering this question!


title: Influence of two factors on non-normal outcome. Which test?

text: I measured two variables (outcomes) $Y_a$ and $Y_b$ for several subjects. $Y_a$ and $Y_b$ are continuous variables defined on $[0,100]$ and have a bimodal distribution. I tried using a transformation but came up with nothing worthy. I have two independent variables $X_1$ (factor with 2 levels) and $X_2$ (factor with 3 levels) characterizing each subject. I would like to investigate if $X_1$ and $X_2$ influence $Y_a$ or $Y_b$. In normally distributed outcomes I would have used `lm(Ya ~ X1 + X2)`. What is the non- parametric equivalent test? Is it the Kruskal-Wallis? Do I have to recode $X_a$ and $X_b$ in one variable?


title: Coupling and Total variational distance

text: Suppose we have two distributions: $\mu$ and $\upsilon$ on $\{1,2,3\}$. $\mu(1) = 1/2, \mu(2) = 1/3, \mu(3) = 1/6,\upsilon(1) = 1/3, \upsilon(2) = 1/6, \upsilon(3) = 1/2$. Could anyone explain me (as simple as possible since I'm at a beginner level in stats) how I can compute the total variation distance of $\mu$ and $\upsilon$? Also, how one should go about constructing a coupling $(X,Y)$ of $\mu$ and $\upsilon$ such that $P(X\not=Y)$ is equal to the total variational distance between $\mu$ and $\upsilon$? Much thanks for helping me understand the concept and problem.


title: Testing Classification on Oversampled Imbalance Data

text: I am working on severely imbalanced data. In literature, several methods are used to re-balance the data using re-sampling (over- or under-sampling). Two good approaches are:   * SMOTE: Synthetic Minority Over-sampling TEchnique (SMOTE)   * ADASYN: Adaptive Synthetic Sampling Approach for Imbalanced Learning (ADASYN) I have implemented ADASYN because its adaptive nature and ease to extension to multi-class problems. My question is how to test the oversampled data produced by ADASYN (or any other oversampling methods). It is not clear in the mentioned two paper how they performed their experiments. There are two scenarios: 1- Oversample the whole dataset, then split it to training and testing sets (or cross validation). 2- After splitting the original dataset, perform oversampling on the training set only and test on the original data test set (could be performed with cross validation). In the first case the results are much better than without oversampling, but I am concerned if there is overfitting. While in the second case the results are slightly better than without oversampling and much worse than the first case. But the concern with the second case is if all minority class samples goes to the testing set, then no benefit will be achieved with oversampling. I am not sure if there are any other settings to test such data. Waiting for your inputs.


title: What are ad-hoc classifiers?

text: I'm reading a paper that mentioned the term: ad-hoc classifiers. I look into the web but couldn't find an answer. I wish someone could help.


title: Data visualisation- summarise 190 means and response rates

text: I have a very large data set which I would like to summarise in as small a space as possible, preferably one side of A4. The data are from a customer satisfaction survey and are Likert-type scales, 5 scales for each work area, with 190 work areas in total. I would also like to represent the response rate on the visualisation somehow, because response rates are very variable and I want management to look at these as well as the actual scores. If necessary I don't mind somehow reducing the 5 scales down to one (using factor analysis or some such thing). One or two sides of A4 to go to the senior management team who are of course very busy with lots of other things and decidedly non-technical. Use of colour is no problem, in fact would probably be seen as a boon. It's just occurred to me that representing the order of the work areas, rather than their absolute value, would be OK, but again I don't want to lose the response rate information. Hope this question isn't too vague, any ideas gratefully received. I am using R and anticipate that this work will involve my learning ggplot2, which I have not as yet got around to.


title: Variance of an ARIMA forecast

text: Suppose I have an $ARIMA(p,d,q)$ model: $ \left(1 - \sum_{i=1}^p \alpha_i B^i \right) \nabla^d X_n = \left(1 + \sum_{i=1}^q \beta_i B^i\right)\epsilon_n$ where $B$ is the shift operator (or lag operator), $\nabla$ is the backward difference operator and $\{\epsilon_n\}_{n=0}^\infty$ is a white noise process. Further, suppose I have observed values $x_0, x_1 , \ldots, x_m$, with $m > \max(p+d, q)$. I want to estimate a future value $\hat x_j$, where $j > m$, using Box- Jenkins. Is there, in general, a closed-form expression for the variance of the estimator $\hat X_j$? _(Not homework. I'm asking this becauseanother question might benefit from this one being answered.)_


title: Exercise 4.2 of Jaynes Probability theory

text: Assume we have three hypotheses: $$A\equiv\text{we have a box with 1/3 defective}\\B\equiv\text{we have a box with 1/6 defective}\\C\equiv\text{we have a box with 99/100 defective}$$ given the prior information $X$, we take: $$P(A|X)=\frac{1}{11}(1-10^{-6})\\P(B|X)=\frac{10}{11}(1-10^{-6})\\P(C|X)=10^{-6}$$ The threshold fraction ($f_t$) is defined like this: as the number of tests ($m$) goes to infinity ($m\rightarrow\infty$) with $f=m_b/m\rightarrow\text{constant}$, $e(C|DX)$ tends to $+\infty$ if $f>f_t$ and $-\infty$ if $f<f_t$. $m_b$ is the number of defective draws, $D$ is the data and $e(A|B)$ is the evidence for $A$ given $B$ defined as: $$e(A|B)=10\log_{10}\frac{P(A|B)}{P(\overline{A}|B)}$$ The exercise is this: given the above what is the threshold fraction? In the book he gives the plausibility flow diagrams in Fig.4.2 using the definition of evidence for two hypotheses at a time: $$e(C|DX)=e(C|X)+10\log_{10}\frac{P(D|CX)}{P(D|\overline{C}X)}$$ where in a two hypothesis case between $C$ and $A$: $\overline{C}=A$. So as shown in Fig.4.2 we easily see that the last term becomes $b=4.73$ if we get a bad draw and $g=-18.24$ if we get a good one. So assuming we get $n_b$ bad draws and $n_g$ good ones the total evidence added to the prior evidence ($e(C|X)$) is: $bm_b+g(m-m_b)=(b-g)m_b+gm=m((b-g)f+g)$ for this to be zero $f$ has to be: $f_t=-g/(b-g)$ placing the values above we get: $0.794079$. This is very similar to the given solution: $0.793951$ but not identical. It also only resulted from considering $C$ and $A$. When I include $B$ through equation 4.39 in the book: $$P(D|\overline{C}X)=\frac{P(D|AX)P(A|X)+P(D|BX)P(B|X)}{P(A|X)+P(B|X)}$$I get $b=7.359974$ and $g=-16.619425$. Placing these in the above equation to find $f_t$ I find: $f_t=0.693070$ which is very much different with the desired answer. I would be very grateful if anyone can help me in understanding how Jaynes derived $f_t=0.793951$?


title: Filtering random guessers from daily surveys

text: Assume we're polling or surveying people about their predictions regarding some future event. This survey is done every day, and the event they are asked to predict is the same every time. For example, each morning they are asked "Do you think it will be sunny this afternoon?" Assume for the purposes of this question that the respondents are being compensated to give an answer, but there is no incentive to be accurate and there is no downside/penalty for being wrong. Some of these people will just guess to get the polling process over with each day and thus earn their compensation. Others will make a legitimate attempt. Some people will vary between these two. Assume that you can compare predictions to the actual outcome after the event has happened. That is, that you can measure each respondent's accuracy. How could you probabilistically tell those that were just guessing from those that were giving their true opinion, keeping in mind that we don't want penalize poor performance. We wish to maintain as much diversity of opinion as possible. But we do wish to eliminate random guessers.


title: Perceived circularity in the Metropolis-Hastings algorithm - Where is my error in reasoning?

text: If I understood it correctly, the Metropolis-Hastings algorithm allows one to sample from a distribution without an analytical representation, which comes in handy, for instance in the Bayesian statistics literature. Quite generally, the algorithm seems to work as follows (this is similar to Koop (2004, p. 93): Bayesian Econometrics):   1. Choose a starting value, $\theta^{(0)}$.   2. Take a candidate draw, $\theta^{*}$, from the candidate generating density, $q(\theta^{(s-1)}; \theta)$.   3. Calculate an acceptance probability based on the candidate draw $\theta^{*}$ and the previous draw $\theta^{(s-1)}$: $\alpha(\theta^{(s-1)}, \theta^*)$.   4. Set $\theta^{(s)}=\theta^*$ with probability $\alpha$ and set $\theta^{(s)}=\theta^{(s-1)}$ with probability $1-\alpha$.   5. Repeat steps 2 to 4 often enough to reach convergence. Later on, they argue that the acceptance probability is a function of the posterior density evaluated at point $\theta^*$ or point $\theta^{(s-1)}$. This, I don't understand. It seems that you need to know the very same posterior density function to compute the acceptance probability that the algorithm is meant to sample from because it is unknown. This seems circular to me. Also, I looked at the Metropolis-Hastings chapter from Robert/Casella (2010): Introducing Monte Carlo Methods with R. Here, I found for instance this example to draw from a beta distribution:               a=2.7; b=6.3; c=2.669 # initial values     Nsim=500000     X=rep(runif(1),Nsim) # initialize the chain     for (i in 2:Nsim){       Y=runif(1)       rho=dbeta(Y,a,b)/dbeta(X[i-1],a,b)       X[i]=X[i-1] + (Y-X[i-1])*(runif(1)<rho)     }          #Here, the algorithm is just an alternative to sampling directly from rbeta     summary(X)       Min. 1st Qu.  Median    Mean 3rd Qu.    Max.      0.00255 0.18970 0.28480 0.30010 0.39520 0.92870      summary(rbeta(Nsim, a, b))         Min.  1st Qu.   Median     Mean  3rd Qu.     Max.      0.001079 0.189500 0.284500 0.299900 0.395000 0.916600      Note that they have to call `dbeta` to compute the acceptance probability. They even mention in their book that this example is just for educational purposes because you can directly sample from `rbeta` in this case. But what are real cases in which I would know a density function (here `dbeta`) but couldn't sample from it (here with the function `rbeta`)?


title: Robust ANCOVA literature

text: I am looking for some source of information, some articles about Robust ANCOVA. Can you suggest me something? And do you know any source that directly gives information about how to perform various robust ANCOVA with R? I have seen a citation about Maxwell et al. 1993, an article called _Analysis of Covariance_ but I do not find it! I know Quinn and Keough, _Experimental Design and Data Analysis for Biologists_ but they do not give much information about Robust ANCOVA.


title: The Central Limit Theorem in Quantile Estimation

text: I don't have a very strong background in statistics so I have a few conceptual questions and there is a strong possibility I'm missing something obvious. Suppose I'm interested in estimating the 99th percentile of body weight in the United States and I have data from every town in every state. I could simply aggregate all the data and find the 99th percentile, but I'm not sure what statistical power this number would have. The data aggregated all together isn't really normally distributed. It has a very large kurtosis excess, so I don't believe any kind of confidence interval estimate would hold. But perhaps I'm missing a basic statistical concept. The second option would be to find the 99th percentile in each town, and apply the central limit theorem, treating each town as an independent random variable. I know the CLT applies to other quantiles, but it doesn't hold for extreme quantiles. I have used MATLAB simulations to prove to myself that it will hold for the 0.99 quantile. However, body weight in each town is not identically distributed. You can imagine that low income, urban, areas will have heavier individuals. So I must apply the Lyapunov or Lindeberg-Feller CLT. Is this a valid thing to do? It seems like these Central Limit Theorems are making statements about the distribution of a random variable divided by the sum of variances, rather than a mean or quantile. How would confidence interval estimates change for these theorems? Any references or insight is greatly appreciated. Thanks.


title: Data visualization: comparing two segments of a bipartite time series

text: Hopefully, this question will make sense to someone. I have a multi-valued continuous time series: $(x, y_1, y_2, y_3, y_4...y_n)$, where $x$ is the time axis. This time series can be divided into 2 contiguous segments of equal length: $(x, y_1, y_2, y_3, y_4...y_n)$ for $0 < x < t_1$   $(x, y_1, y_2, y_3, y_4...y_n)$ for $t1 < x < t_2$ where $[0, t_1]$ is a period of same length as $[t_1, t_2]$ My question is: What would be the best way to **visualize** the differences between the 2 segments? To make things less cluttered and more feasible in some kind of 2D chart or graph, it is fine to assume the series are $(x, y_1, y_2, y_3)$ series.


title: Why divide by square root of K in calculating Stouffer Z?

text: If Stouffer's method of combining k probabilities is based on the variance of the sum being equal to the sum of the variances, and the variance of z is 1, why does the Stouffer formula divide by the square root of k rather than by k?


title: Handling outliers when comparing two means in a repeated measures design

text: I am doing a simple study that involved taking a measure at time point 1 and time point 2 (12 weeks later). While the sample was a class, not all members were present at both time points, so I have 20 date points at time 1 and 21 date points at time 2. The measure has a score, and I am taking the means and doing a simple t-test to determine if the intervention caused any increase in the measure at time point 2. Questions: Do I need to throw out outliers if they are more than 2 standard deviations higher than the mean? When I do the t-test, do I need to look at one- or two-tailed distributions? My hypothesis is that the intervention will increase the mean at time 2, so I think I should consider a one-tailed distribution. Lastly, I am assuming that I have to do a paired t-test since it is a repeated measures design.


title: Impact of lagged values on identity variable

text: Let's say I'm working with the following simplified macroeconomic accounting identity                Y_t = C_t + I_t + G_t,      meaning that GNP in time `t` equals consumption in `t`, `C_t`, plus investment in `t`,`I_t`, plus government spending in `t`, `G_t`. I'd like to study the impact of past government spending on current consumption, that's                C_t = f( G_{t-1}, G_{t-2}, ... )      so I can make a statement like **"past government spending accounts for x% of current consumption"**. What techniques should I use to run such an analysis while still meeting the equation identity above? Thank you.


title: Bayesian rating system with multiple categories for each rating

text: I'm implementing a rating system to be used on my website, and I think the Bayesian average is the best way to go about it. Every item will be rated in six different categories by the users. I don't want items with only one high rating to shoot to the top though, which is why I want to implement a Bayesian system. Here is the formula:               Bayesian Rating = ( (avg_num_votes * avg_rating) + (this_num_votes * this_rating) ) / (avg_num_votes + this_num_votes)      Because the items will be rated in 6 different categories, should I use the average of the sums of those categories as "this_rating" for the Bayesian system? For instance, take one item with two ratings (scale of 0-5):               Rating 1:       Category A: 3       Category B: 1       Category C: 2       Category D: 4       Category E: 5       Category F: 3       Sum: 18          Rating 2:       Category A: 2       Category B: 3       Category C: 3       Category D: 5       Category E: 0       Category F: 1       Sum: 14      Should "this_rating" be simply the average of the sums listed above? Is my thinking correct, or should a Bayesian system be implemented for each category as well (or is that overthinking it)?


title: How to derive errors in neural network with the backpropagation algorithm?

text: From this video by Andrew Ng around 5:00 ![enter image description here](http://i.stack.imgur.com/2TnT1.png) How are $\delta_3$ and $\delta_2$ derived? In fact, what does $\delta_3$ even mean? $\delta_4$ is got by comparing to y, no such comparison is possible for the output of a hidden layer, right?


title: 'mean' of logit-normal distribution

text: $ X \sim N(\mu,\sigma^2) $ $ Y = \frac{\exp(X)}{1+\exp(X)} $ Y has a logit-normal distribution. When I have a correct estimate of the mean, say $\bar{Y}$, of the logit-normal distribution, how can use it to get a correct estimate of $\mu$ ? I don't have full sample from logit-normal distribution, but only the correct estimate of the mean. I know that $ \log(\dfrac{\bar{Y}}{1-\bar{Y}}) $ is not the right answer due to Jensen's inequality.


title: Inference from simulated quantile function

text: I generated a quantile function $\hat X$ using Monte Carlo simulation. The random variable I simulate is the mean value of 5 draws from an i.i.d. range statistic $Y$. I.e., I have $Y(\sigma) \sim \sigma F()$, and I simulated the value of $X(\sigma=1) \sim \sum_{1}^{5} y(1)_i / 5$. Is it valid for me to use this quantile function for statistical inference on real-world samples from $X$? In particular, given a sample $x_i$:   1. Is it correct to say that $x_i$ is an estimator of the mean of $\bar X$, and therefore our best guess is that $\hat\sigma = x_i/\overline{\hat X(1)}$?   2. Is it correct to say that our confidence intervals for $x_i$ are given by the simulated quantile function? E.g., the 90% confidence range on $x_i$ is $[\hat X_{.05}(\hat \sigma), \hat X_{.95}(\hat \sigma)]$?


title: Hypothesis testing that one time-series of a measure of entropy doesn't belong to a population

text: Disclaimer: I'm a software engineer, not a statistician, so please forgive any blunt error :-) I have a set of time-series "curves", each measuring the entropy of a given artifact. Now, I'm standing over the following premises (please criticize them as you see fit):   1. In order to approximate the upper bounds of the Kolmogorov complexity $K(s)$, of a string $s$, one can simply compress the string $s$ with some method, implement the corresponding decompressor in the chosen language, concatenate the decompressor to the compressed string, and measure the resulting string's length.   2. For this purpose, I've used the _bzip2_ application, setting its compression level to the supported maximum (-9).   3. If one is only interested in a time-series analysis of a set of evolving strings, calculating the compressed _deltas_ is enough to present a relative measure of entropy between any two strings (at least that's my interpretation after reading Cilibrasi05).   4. For that, I used the _diff_ unix tool, with the (--minimal) parameter, again followed by a _bzip2_ compression, with the aforementioned settings. I'm doing this to analyze the evolution of the entropy in a software artifact (code, model, whatever). I'm _not_ worried with the absolute values, but with the _relative_ increase (or decrease) in entropy. Now here comes the problem:   1. I've done this for a set of 6 artifacts, which ought to belong to the same population, but I don't know how to provide statistical evidence of that (the corresponding of doing a two-tailed t-test of two samples).   2. One of the artifacts evolution _should_ be different from all the others. We're talking something like an exponential v.s. sub-linear growth. How do I provide statistical evidence of that? Again, the disclaimer of being a software engineer. Although I would appreciate every academic reference (papers, books, etc.) you could handle, I'm looking for something pragmatic that I can use in the next few days, like a script in R, or something in SPSS. P.S. I'm sorry for asking for a _recipe_ , instead of a theoretical explanation.


title: What criteria are used to compare feature-based classification techniques?

text: When comparing feature-based classification techniques, what characteristics about the different processes should be considered? I'm comparing different classification techniques to try to figure out what should be considered when selecting a classification technique. In particular, I'm comparing the following:   * Linear classifiers   * K nearest neighbors   * Gaussian mixture models   * Bayesian networks   * Rule induction / decision trees * * * **Edit** When I am considering which classifier I should use, I think about the following:   * Training time   * Classification time   * Classification accuracy   * Generates a probability vs. giving a class What other things should be considered?


title: What are the pros and cons of using the logrank vs. the Mantel-Haenszel method for computing the Hazard Ratio in survival analysis?

text: One way to summarize the comparison of two survival curves is to compute the hazard ratio (HR). There are (at least) two methods to compute this value.   * Logrank method. As part of the Kaplan-Meier calculations, compute the number of observed events (deaths, usually) in each group ($Oa$, and $Ob$), and the number of expected events assuming a null hypothesis of no difference in survival ($Ea$ and $Eb$). The hazard ratio then is: $$ HR= \frac{(Oa/Ea)}{(Ob/Eb)} $$   * Mantel-Haenszel method. First compute V, which is the sum of the hypergeometric variances at each time point. Then compute the hazard ratio as: $$ HR= \exp\left(\frac{(Oa-Ea)}{V}\right) $$ I got both these equations from chapter 3 of Machin, Cheung and Parmar, Survival Analysis. That book states that the two methods usually give very similar methods, and indeed that is the case with the example in the book.  Someone sent me an example where the two methods differ by a factor of three. In this particular example, it is obvious that the logrank estimate is sensible, and the Mantel-Haenszel estimate is far off. My question is if anyone has any general advice for when it is best to choose the logrank estimate of the hazard ratio, and when it is best to choose the Mantel- Haenszel estimate? Does it have to do with sample size? Number of ties? Ratio of sample sizes?


title: Why use a z test rather than a t test with proportional data?

text: When testing the difference between 2 proportions why do we use a z-test rather than a t-test? Further, is there a simple way to conduct an omnibus test for significant differences between more than 2 proportions (in the form of percentages). Is there an equivalent to a one-way ANOVA for this? I imagine that you could use logistic regression (assuming you have the original data for the proportions in the form of 0s and 1s) Are there other options?


title: Computing Out of Bag error in Random Forest

text: I am implementing a Random Forest classifier as a side-project, and I am a bit unclear on what the correct approach is to compute, say, the OOB estimate for the classifier error rate. My understanding is that typically, for each tree in the forest, one creates a training sample from the original sample by taking Examples with repetition, and what is left out can be used to compute out-of-bag estimates. The part I am unclear about is how to aggregate the errors across the different out-of- bag samples. The naive approach would be for each tree to count how many OOB examples are mis-classified, and compute the average mis-classification rate over all of them (total mis-classified / total Examples out-of-bag). However, it seems to me that in essence, this would be computing the average classification error of each of the individual trees, and missing the fact that the forest is taking a majority vote over the verdict of each tree, compensating for "weaker" trees. A more complicated way would be to take each OOB Example, look up for each tree if it was included or not in the training, and take a majority vote over all trees that didn't use that Example for training. This is computationally more painful, but seems like it would properly take into account the fact that a forest is more than the sum of its parts. Can anybody enlighten me on which of the two approaches is the correct way to approach this, or, if neither is correct, what I should be doing instead?


title: Is there any proper way to fix a sample to adjust for known demographic overrepresentation?

text: My spouse frequently works with (expensive, hard to obtain) data samples; for example route information for commuting bicyclists collected using a smartphone app. More often than not, these samples suffer from some kind of known demographic over-representation that they'd like correct for various applications. Mindful of Karl Roves' and friends "corrections" to "obvious" democrat oversampling for the Nov. 2012 election polls which led to rather embarrassingly incorrect predictions, is there any theoretically appropriate way of doing this? I'm not even sure what to call what I'm looking for -- is this what in some places is called reject inference?


title: In decision tree construction, can a good splitter have low information gain?

text: I have a data set with a candidate splitter variable that is a natural choice from the business perspective. It has two values, and the distributions of the target when conditioned on the two values of this variable are quite distinct (visually). Yet the information gain of this variable is very small (almost 80% smaller than some other variables whose conditional distributions are far less differentiated.) Is it possible that a variable with low information gain is nevertheless a good splitter in some other sense?


title: perturbing discrete probability distribution

text: What is the best way to slightly perturb a given discrete probability distribution ? Adding a zero mean Gaussian noise to the probability distribution and re- normalizing it such that it sums to 1 is one way. But it does not ensure that all the probabilities are positive. Is there an other way to achieve this ?


title: Can feature selection be considered a way to observe relationship between variables like correlation?

text: In correlation we can observe relationship between a pair of variables, let me call it X1 and Y. Now, considering I have the predicting variables X1, X2, ..., Xn and the variable Y. Does the following assumption holds: If the variables Xi,...,Xj are observed to better classify Y, where 1 <= i,j <= n, can I claim that the variables Xi,...,Xj have more association with Y rather the remaining variables Xk, k not belonging to the i,j interval.


title: Modifying the time granularity of a state sequence

text: I created a sequence object from my SPELL-formatted data set. The sequence length of the sequence object is 1440 (i.e., 1-min intervals for a day). Is there any easy way for TraMineR to convert the sequence length from 1440 to 288 (i.e., 5-min intervals), using the status in the first minute of every time interval to represent an individual's status during that 5-min period?


title: How to connect nodes in full conditional model of GRN?

text: I have problem in intuitive understanding of full conditional model in building gene regulatory network. it's basically asking that "can correlation between two genes be explained by all other genes in the model?" . and say that, Node i and j are connected if and only if Xi ⊥ Yj | Xrest . my question is that, if Xi is independent of Yj given Xrest, why we should draw edge between them ?


title: Estimators, sufficiency, consistency, and bias

text: > A random variable is said to have the Pareto distribution with parameters > $\alpha$ and $\beta$, $P(\alpha, \beta)$, if its cumulative distribution > function is given by > > $$F(x)= 1 - (\frac{\beta}{x})^{\alpha},$$ if $x \geq \beta$, and $0$ > otherwise. > > a) If $X_1, ... ,X_n$ are iid $P(\alpha, \beta)$ find the probability > density function of min$(X_1, ... , X_n)$. Let's rewrite the random variables as $Y_1, ... ,Y_n$, where $Y_1 < Y_2 < ... < Y_n$. Then $F_{Y_1}(y) = P(Y_1 \leq y, ... , Y_n \leq y)$   $$= 1 - P(Y_1 \geq y, ... , Y_n \geq y)$$ $$= 1 - P(Y_1 \geq y) ... P( Y_n \geq y)$$ $$= 1 - [P(Y_1 \geq y)]^n$$ So it suffices to find $P(Y_1 \geq y)$. But $P(Y_1 \geq y) = 1 - P(Y_1 \leq y) = 1 - (1 - (\frac{\beta}{y})^{\alpha}) = (\frac{\beta}{y})^{\alpha}$ So, $f_{Y_1}(y) = (\alpha n)\beta^{\alpha n}y^{-\alpha n -1}$ > b) If $X_1, ... ,X_n$ are iid $P(\alpha = 3, \beta)$, show that $\beta = min > (X_1, ... ,X_n)$ is: > > i. a biased estimator for $\beta$. Compute the bias of $\hat{\beta}$ $E(Y_1) = \int^{\infty}_{\beta} 3n \beta^{3n} y^{-3n} dy$ $$=lim_{u \rightarrow \infty} [\int^u_{\beta} 3n \beta^{3n} y^{-3n} dy]$$ $$=lim_{u \rightarrow \infty} [\frac{3n}{-3n+1} \beta^{3n} y^{-3n+1} ]^u_{\beta}$$ $$= \frac{-3n}{-3n+1 \beta}$$ Now we compute the bias $E(Y_1 - \beta)$ $$= E(Y_1) - E(\beta)$$ $$= [\frac{-3n}{-3n+1 \beta}] - [\int^{\infty}_{\beta} 3n\beta^{3n+1}y^{-3n-1} dy] $$ $$= [\frac{-3n}{-3n+1 \beta}] - [\frac{3n \beta^{3n+1} y^{-3n}}{-3n}]^{\infty}_{\beta}= [\frac{-3n}{-3n+1 \beta}] +\beta$$ > ii. a consistent estimator for $\beta$ This one seems confusing to me. According to the textbook, an estimator is consistent if it approaches the estimated value in probability. However, we have $\lim_{n \rightarrow \infty} P[|Y_1 - \beta| < \epsilon]$...which doesn't seem right to me, because there is no $n$ in the expression $|Y_1 - \beta|$. So I'm guessing if I'm on the wrong track... > iii. a sufficient statistic for $\beta$. It is sufficient by Neyman's theorem, because we can write the pdf as $f_{Y_1}(y) = [3n] [\beta^{3n}y^{-3n-1}]$. Are my answers correct? If not, can you give me a hint? Also, can you help me with part ii? Thanks in advance


title: Intuitive understanding of Local Probability Distribution

text: I'm learning Bayesian network. I have problem in intuitive understanding of Local Probability Distribution. Can anybody explain to me what it is?


title: Which test to use for a set of data

text: I have data from groups of professionals with in the same medical field. They are grouped according to years in practice. Each group going up in 5 year increments.ie 0-5 yrs experience 6-11 etc. Individually all were asked a set of questions from which I was only interested in one response based on my hypothesis. I'm wanting to see, if the more experienced in practice an individual is, the better they are at spotting a disease. What is the best statistical test to use when comparing the nine age groupings against each other in relation to the frequency each group has for the 1 question I zoomed in on?


title: How to compare two classifications' classes' size?

text: I want to know if, between two classifications, one provides more homogeneous classes size than the other. I devised to use the chi square test (classes' size are far more than the rule of thumb of 5 or 10 frequencies), to compare each classification to the most homogeneous classification. In example, if I have two classifications such as:                     1     2     3     4     5     A    31    33    45    25    56     B    29    30    41    30    60      Over 190 items classified, I compared them (using chi square test) with                     1     2     3     4     5     H    38    38    38    38    38      Now I have two pvalues. May I assume that, if A's pvalue is greater than B's pvalue, A's classes are more homogeneous than B's, with respect to size? Thanks


title: covariance of RVs under a nonlinear transformation

text: I have a multivariately distributed random 3-vector $$\mathbf{x}=\begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix}\sim\mathbf{N}\left(\begin{pmatrix}2\\-1\\3\end{pmatrix},\begin{pmatrix}4&1&0\\1&2&1\\0&1&30\end{pmatrix}\right)$$ and I'd like to compute $\text{cov}(x_1x_2,x_1x_3)$. I don't think the fact that $\text{cov}(x_1,x_3)=0$, i.e. the fact that $x_1$ and $x_3$ are independent, implies the covariance I want also vanishes, but I'm stuck about how to proceed because of the non-linearity of the transformation. Suggestions? I assume the problem simplifies immensely because of the normality assumption and that the general question (without the distribution assumption) would have no closed form solution, but I still don't know how to proceed.


title: Appropriate analysis on non normal data with different variance

text: Is the Kruskal-Wallis test valid when applied to a subset of data? Is there a more appropriate test for comparing the effects of two variables on non- parametric, incomplete data with different variances? Super simple version I analyse my data using method A. A has 6 values. I have two variables L and P L has 3 values: L1, L2, L3. P has 4 values, P1-P4. I have data for 9 combinations L1P1, L1P2, L1P3, L1P4, L2P1-4 and L3P1. I have 10 samples for each combination. I am interested in the effect of P on L for a given value of L. (is L1P1 different to L1P2-4? and is L2P1 different to L2P2-4?) I also want to compare L1-3. The distribution of values is non-normal and the variance differs between combinations. I then compare the results of A1-A6 (This is likely to be a qualitative assessment so I don't think I need a statistical test here.) My hypothesis is that A correlates quantitatively with the experimental data and predicts the correct value of P, furthermore the value of P will be constant for L1 and L2. More detailed version. Including where I have got to with my understanding of stats. I have 4 proteins. (Not really, it is the same protein in different conformations but for the sake of clarity and to avoid jargon I have treated them as different in the explanation below) I have 3 ligands/drugs. I have 9 initial structures. 2 ligands combined with each protein + final ligand in one protein. I have conducted 90 simulations. 10 of each initial structure. Each simulation is stochastic/determinist/reproducible to binary accuracy (not random) consisting of 476 data points. But the initial conditions were randomised (initial velocity of each atom using a seed value to randomise), therefore each run is reproducible within itself but different compared to other runs of the same initial structure. I have biological data (inhibition values) for my compounds. I want to know how well my simulations predict that data. There are several parameters that I have measured in each simulation. I want to compare how well those parameters predict the experimental data. One flaw in my experimental data is that it cannot distinguish between the proteins for the first two drugs. My simulations were trying to predict which protein is preferred. My parameters are of two types: single points and averages. Some of the parameters are single measurements taken during the simulation, lowest energy is one example. Other parameters are averages of a value throughout the simulation, e.g. average energy. Ideally I believe I would use a 2-way anova to conduct the analysis. I have three problems with this approach. Firstly, the variance of the ten simulations of each initial structure differs. This differing variance is seen for all parameters I have measured (Bartlett's test for equal variances). Secondly, I have no reason to believe these should be normally distributed measurements, this is borne out by various normality tests (D'Agostino & Pearson omnibus normality test) I have used on the data; (the average value parameters are closer to a normal distribution than the single point parameters) Finally the software I am using (Graphpad) doesn’t like the fact I am missing data. It needs/wants the other three proteins for the final ligand. For the first two reasons I cannot use a 1-way anova. So I have been using a Kruskal-Wallis test, however I have concerns with how I have been applying it. I have tabulated my data in 9 columns x10 rows and have conducted the test on everything. I want to test each of the drugs alone. (each drug against its four proteins to determine the preferred protein) Then I want to compare drugs. (each drug in its preferred protein against the other drugs in their preferred protein) My concern is that if I do a pairwise comparison (Dunn's Multiple Comparison Test) for the first step, within the test spanning all initial combinations then the variance of the other drug and its proteins are affecting the results. If I perform the Kruskal-Wallis comparison on each drug separately then compare the 3 best I think I might be doing the equivalent of a multiple T-test and hence type 1 errors. My co-supervisor has suggested that Bayesian might be able to answer these questions, but I have no experience at all with the technique and based on my reading I have no idea how to form the priori.


title: Mann-Whitney for non-normal distributions with n>20?

text: I have 12 groups of datasets (mostly non-Normal, quite right-skewed), and I would like to see if there is a significant difference within each pair. (But the data itself is not a paired dataset, as in dataset #1 and dataset #2 that I am comparing are from independent populations.) My data also has outliers (following the 1.5*IQR rule), so I understand that the t-test is not valid because (x bar) and s are not robust against outliers. The main alternative I see is to use the rank test Mann-Whitney U test. (Let me know if I am not on the right track.) In Zar's _Biostatistical Analysis_ (4th ed, pg. Ap100): "For the Mann-Whitney test involving n1 and/or n2 larger than those in this table, the normal approximation may be used." The table goes to n1=20, n2=40. In Moore, McCabe, Craig's _Introduction to the Practice of Statistics_ (6th ed., pg. 432): For sample sizes 15≤ n ≤39, " _t_ procedures can be used except in the presence of outliers or strong skewness." For samples sizes ≥40, " _t_ procedures can be used even for clearly skewed distributions." Does this mean if n for one of my populations ≥20 (or ≥40?), I will have more accurate/powerful results using the 2-sample t-test (assuming unequal variance) instead? What if n1 ≤20 but n2 ≥20? If I need to use the Mann- Whitney U test for some, should I be using Mann-Whitney U test for all in order to compare the results? Are there other hypotheses tests I should be using instead? For reference, the populations that I am comparing have counts:               n1 = 15, n2 = 20      n1 = 7,  n2 = 27      n1 = 12, n2 = 11      n1 = 13, n2 = 4      n1 = 22, n2 = 47      n1 = 25, n2 = 15      n1 = 20, n2 = 21      n1 = 12, n2 = 27      n1 = 22, n2 = 22      n1 = 26, n2 = 14      n1 = 32, n2 = 48      n1 = 48, n2 = 36       Your thoughts are appreciated. Thank you so much! FYI I have been using Excel's 2-sample t-test assuming unequal variance and Vassarstat's online Mann-Whitney test: http://vassarstats.net/index.html (under the heading _Ordinal Data_ ).


title: Is there a standard name for probabilistic graphical models like this?

text: I have some data I need to analyze, and some prior knowledge I'd like to apply. I have discrete time points, and noisy, continuous outputs $x(t)$ that I observe at those time points, and would like to infer the discrete state $y(t)$ that is most likely to correspond to $x(t)$. If it matters, I know the distributions $x(y_1),\cdots,x(y_n)$, and that there are only six possible hidden states, and have lots of annotated data to train the model on. Were it that simple, I think a hidden Markov model and the Viterbi algorithm would do the trick. But I know that the output at time $t_i$ is affected not only by the hidden state at time $t_i$, but by the three previous hidden states, like so: ![link to picture of model](http://i.stack.imgur.com/VbGXt.png) I left the diagram undirected, since the correlations between successive hidden states are not causal, but I could live with rightward arrows all around if I have to. Anyway, does anyone know what such a model might be called, and/or if it's even plausible to try to learn the parameters for it? It seems simple enough that I feel like someone has got to have done it before, but I haven't been able to locate anything. I think it's basically a linear-chain conditional random field, but those extra edges make me feel lost and confused. **EDIT:** If it's any easier to parse/think about, here's an example where the current observed state depends on the past two, rather than four, hidden states: ![Simplified model](http://i.stack.imgur.com/xMQ3L.png) This is sort of the minimal extension of a linear chain CRF, but I still can't find an example of something like it. Any help?


title: How do I compute a cutoff based on sensitivity/specificity when the characteristics of my sample is different from the population?

text: I have a dataset containing the performance of a novel instrument to screen for disease A. The novel instrument uses a scoring system to score the subject to determine if they have disease A. I then proceed to use logistic regression to use the novel instrument and demographic variables such as age and gender to classify them into either having or not having disease A and comparing it with the gold standard result. Since the prevalence of disease A is different in my sampling and in the actual population, how should I handle this problem? I see a couple of possible solutions:   1. Conduct the logistic regression and use the predicted values for each subject to compute for the sensitivity and specificity of the regression.   2. Conduct the logistic regression in the sample. Sample from my existing sample to create a new sample that corresponds to the population proportion of people with and without disease A. Perform the prediction to compute for the sensitivity and specificity of the regression. My questions are:   * Should I adopt either strategy (1) or (2)?   * During my search for answers, I came across: Case weighted logistic regression. Following the thread on the mailing list, I see a way to make replicate weights for each subject (via `survey` package). Why should the proportion (in my case) of people with and without disease A affect the results of the logistic regression? Is the replicate weights method a valid way of overcoming this issue or would strategy (2) work as well?


title: Create a matrix of tf-idf values from documents

text: I have a set of `documents` like:               D1 = "The sky is blue."     D2 = "The sun is bright."     D3 = "The sun in the sky is bright."      and a set of `words` like:               "sky","land","sea","water","sun","moon"      I want to create a matrix like this:                  x        D1           D2         D3     sky         tf-idf       0          tf-idf     land        0            0          0     sea         0            0          0     water       0            0          0     sun         0            tf-idf     tf-idf     moon        0            0          0      Something like the example table given here: http://www.cs.duke.edu/courses/spring14/compsci290/assignments/lab02.html. In the given link, it uses the same words from the document but I need to use the set of `words` that I have mentioned. If the particular word is present in the document then I put the `tf-idf` values, else I put a `0` in the matrix. Any idea how I might build some sort of matrix like this? Python will be best but R also appreciated. I am using the following code but am not sure whether I am doing the right thing or not. My code is:               from sklearn.feature_extraction.text import CountVectorizer     from sklearn.feature_extraction.text import TfidfTransformer     from nltk.corpus import stopwords               train_set = "The sky is blue.", "The sun is bright.", "The sun in the sky is bright." #Documents     test_set = ["sky","land","sea","water","sun","moon"] #Query     stopWords = stopwords.words('english')          vectorizer = CountVectorizer(stop_words = stopWords)     #print vectorizer     transformer = TfidfTransformer()     #print transformer          trainVectorizerArray = vectorizer.fit_transform(train_set).toarray()     testVectorizerArray = vectorizer.transform(test_set).toarray()     #print 'Fit Vectorizer to train set', trainVectorizerArray     #print 'Transform Vectorizer to test set', testVectorizerArray          transformer.fit(trainVectorizerArray)     #print     #print transformer.transform(trainVectorizerArray).toarray()          transformer.fit(testVectorizerArray)     #print      tfidf = transformer.transform(testVectorizerArray)     print tfidf.todense()      I am getting very absurd results like this (values are only `0` and `1` while I am expecting values between 0 and 1).               [[ 0.  0.  1.  0.]      [ 0.  0.  0.  0.]      [ 0.  0.  0.  0.]      [ 0.  0.  0.  0.]      [ 0.  0.  0.  1.]      [ 0.  0.  0.  0.]      [ 1.  0.  0.  0.]]         I am also open to other libraries for calculating `tf-idf`. I just want a correct matrix which I mentioned above.


title: Statistical analysis on categories before text classification

text: I want to classify text by different topics. However, one of the current problems is that there are several topics/categories that are quite intuitively independent and statistically standalone, but there are also several other topics/categories that are quite intertwined with each other. In order to make a good decision on the final categories, I realize I'd better carry out some statistical analyses on my raw text collection and explore the statistical nature of the collection. I was thinking using clustering to have a clearer view of my text collection. But I want to know, what are options that I get in this step to make a better decision on my categories? Any help and suggestion would be appreciated.


title: Comparing standard deviations between variables with very different ranges

text: I have a number of variables for which I want to compare the standard deviations. Each of these variables is actually the data for a particular wavelength of light, and I want to produce a plot showing the standard deviation for each wavelength. However, the values that I am calculating my SD from are much higher for some wavelengths than others - for example, at 500nm I get values of around 50, and at 1000nm my values are more like 1. In this case, comparing standard deviations seems rather silly as a SD of 5 at 500nm is actually less variable than a SD of 5 at 1000nm. In this situation would it make sense to divide the SD by the mean to give a sort of percentage difference standard deviation measure? Or is there some other better way (or 'proper') way of doing this? If I were to do this, would it make statistical sense?


title: Difference between bias-variance dilemma and overfitting

text: I'm wondering what difference it makes whether we talk about bias-variance dilemma where fitting a regression line to the given dataset reduces bias and increases variance or whether we talk about overfitting where the model is fitted to the dataset and does not generalize to new datasets. Essentially overfitting implies the same tradeoff between biasedness and variance doesn't it?


title: Probability of 2 uniformly randomly generated, n-digit numbers (n>=4) having the same last 4 digits?

text: I apologize if question is too simple! For example, is the probability of generating 2, 8-digit numbers that have the last 4 digits match the same as the probability of the last 4 matching with a 12 digit number? **For an 8 digit number:** 10^8 total possible numbers and 10*10*10*10*1*1*1*1 = 10^4 / 10^8 = 1/10,000 probability of generating 2 with the same last 4. **For 12 digit:** 10^12 total possible. 10^8 / 10^12 = 1/10,000 probability again of generating 2 with the same last 4? **Is this correct?**


title: Conditional Expectation has the minimum mean prediction error

text: Consider $\{X_t\}$ as a general time series including random variables $X_t$. Assume that we have observed $X$'s until time t. The goal is to come up with a function of the observed $X$'s to predict a future X, let's say unobserved $X_{t+h}$. Condsider $h(\vec{X}) = h((..., X_{t-1}, X_t))$ as our proposed function. We define mean prediction error as: $$E\{(X_{t+h} - h(\vec{x}))^2 | \vec{x}\}$$ **Goal: I want to show that among all possible functions h, the function $E(X_{t+h}|\vec{x})$ has the minimum mean prediction error**. Thank you for your help


title: Are categorical variables standardized differently in penalized regression?

text: In penalized/regularized regression (lasso, ridge, etc.) the predictors are typically standardized to be centered at 0 and often to have variance 1. Are categorical predictors treated differently. If so, why? What are the consequences of using the same standardization? Is a reference available?


title: Can there be endogeneity?

text: I am performing a regression in which I suspect there could be endogeneity between the main explanatory variable and the dependent. As a first step I perform an OLS regression and then, to control for the endogeneity, I use a 2SLS regression with an instrumental variable. The coefficient for this variable in my OLS regression is not statistically significant. Does it still make sense to perform a 2SLS? Can it be that there is some 'counterbalancing' effect for which the relationship between the dependent and the explanatory variable is actually there but non-observable and not detected with the OLS? I apologize if this is a stupid question but I am not an expert in this field and need some hints to carry on with my analysis :)


title: Statistical grounding for difference between median as a proportion of overall visible spread

text: I'm helping someone with prep for their stats assessment and they are being instructed to use the following to make an informal inference about the population: ![Figure](http://i.stack.imgur.com/ok0Fd.png) I was wondering if there is any statistical basis to these figures or how it was decided this is a good indicator.


title: Continuous time markov chain backward/forward equations

text: ![enter image description here](http://i.stack.imgur.com/DOMHz.jpg) Using Kolmogorov's forward and backward equations, show that $p_{11}(t) + p_{21}(t) + p_{31}(t) = 1$ and $p_{21}(t) = p_{31}(t)$ where $p_{ij}(t) = P(X(t) = j | X(0) = i)$. **My attempt:** I can show the first part but not the second, that is, for all $t \ge 0$, the forward equations for $p_{11}$, $p_{21}$, and $p_{31}$ are as follows: \begin{gather} p_{11}'(t) = \lambda p_{12}(t) + \lambda p_{13}(t) - 2\lambda p_{11}(t) \cdots (1)\\ p_{21}'(t) = \lambda p_{22}(t) + \lambda p_{23}(t) - 2\lambda p_{21}(t) \\ p_{31}'(t) = \lambda p_{32}(t) + \lambda p_{33}(t) - 2\lambda p_{31}(t) \end{gather} The backward equations for $p_{11}$, $p_{21}$, and $p_{31}$ are as follows: \begin{gather} p_{11}'(t) = \lambda p_{21}(t) + \lambda p_{31}(t) - 2\lambda p_{11}(t) \cdots (2) \\ p_{21}'(t) = \lambda p_{11}(t) + \lambda p_{31}(t) - 2\lambda p_{21}(t) \\ p_{31}'(t) = \lambda p_{11}(t) + \lambda p_{21}(t) - 2\lambda p_{31}(t) \end{gather} Subtracting $(1)$ from $(2)$ yields, \begin{gather} 0 = p_{21}(t) - p_{12}(t) + p_{31}(t) - p_{13}(t) \end{gather} Noting that $1-p_{11}(t) = p_{12}(t) + p_{13}(t)$, we have: \begin{gather*} 0 = p_{21}(t) + p_{31}(t) - (1-p_{11}(t)) \\ \implies p_{11}(t) + p_{21}(t) + p_{31}(t) = 1 \end{gather*} However, I'm not exactly sure how to show $p_{21}(t) = p_{31}(t)$ from algebraic manipulations... (I can show it by actually solving the differential equation but the question wants me to show it through manipulation)


title: Cross validation - how many samples?

text: I'm not a statistics guy but want to investigate a problem like this: I have a number (hundreds) of trends that I wish to use to come up with an average trend. My method was to be by using the average value of the individual trend at each time step, to come up with the average, then fit this new average trend with a smoothed line. From this original set of hundreds, I then wish to take a subset of the original and an average of those the same way, fit a smoothed trend to this average and compare the resulting trend to the original and repeat that procedure. I think this method is "Repeated random sub-sampling validation" (Wikipedia). My questions are:   * is there a rule (of thumb?) as to how many sub-samples I should take from the original (i.e. how many times should I repeat the process)   * how large this sub-sample group should be in relation to the original group (how may trends from the original to include in the smaller subsets)? I hope this makes sense.


title: Naive Bayes non-Dictionary Term in Test Document

text: Using Laplacian Smoothing we can get rid of 0 probabilities if a term occur in spam and does not occur in ham class or vice versa. My question is about what if a term in test document does not occur in training dataset(i.e. in dictionary). For example if we extend example in page 44 in -> http://www.stanford.edu/class/cs124/lec/naivebayes.pdf as follows: (Consider only change is test document part) Test 5 Chinese Chinese Chinese Tokyo Japan Greek ? P(c|d5)'s calculation requires P(Greek|c). However since Greek doesn't exist in dictionary we didn't calculate it before. What should P(Greek|c)'s and P(Greek|j)'s value?


title: How to best save intermediate results computed from an R data.frame?

text: I have a set of data created by users answering a questionnaire. I imported their answers from a .csv file and got them as a data frame with one user per row and one question per column. However, the questions were not homogenous. I have to first evaluate some questions, which gives me an ordered list of the users' preferences for 10 predefined categories. Then I want to evaluate the remaining questions, and for this I have to use some information of this ordered list (for example, which is the category the user ranked highest?). I calculated the score for each category. It is currently kept in a matrix, which looks like that:                     cat1.score cat2.score    ...    cat10.score      user1   2.50       2.25        ...     3.20     user2   3.85       2.05        ...     2.40      and I plan to create lists and sort them, so I'll get for the first user a list like `preferences <\- list("cat10", "cat1", "cat2", ...)` (assuming that the scores not shown are lower than 2.25). But I am not sure how to structure the information. My plan was to create a new data frame, which will have the same data as the matrix, and its eleventh column will hold the list of the categories ranking for the user. I tried lots of ways to construct such a data frame, but couldn't do it. Now I am very new to R, so I guess that I haven't thought of all ways yet, and I could try a lot more. But as I saw how hard it is to do, I guessed that maybe I am trying to do something which makes little sense - if it was a good practice, R would have probably had a convenient mechanism to do it, or the tutorial books would have had an example. So, my question is, is this a good way to structure my intermediate results? And if not, what is a better way? I get one such list per user, and I really need it as ordered data (for each user, I will later have to access it as `preferences[2]` and get the category the user liked second most, or similar). To make it clear, I know which data structures in R can contain a list and which can't. My question is not what the language will let me do, but what is the sensible thing to do here.