title: Differences in power in regression versus mean of a ratio

text: I've been futzing around with the following example, and feel that there must be a good explanation or reference that I am blindly missing. Supposer you have a linear relationship between two variables. Let's even assume that the intercept goes through 0. Why would there be more power in analysis by regression of y~x versus y/x ~ 1? And in particular, if x is biased, the decrease in power is even more dramatic. I can show this by simulation in the following R code, it strikes me that there is easy explanation here.               simBiasPower <- function(b=-1, ySD = 2, n=100, nsim=100){       sapply(1:nsim, function(i) {         unbiasedX <- runif(nsim, 0,10)         biasedX <- rexp(nsim,0.5)              unbiasedY <- rnorm(100, b*unbiasedX, ySD)         biasedY <- rnorm(100, b*biasedX, ySD)                   return(c(unbiasedSlopeP= coef(summary(lm(unbiasedY ~ unbiasedX)))[2,4],                  unbiasedRatioP = coef(summary(lm(unbiasedY/unbiasedX ~ 1)))[1,4],                  biasedSlopeP = coef(summary(lm(biasedY ~ biasedX)))[2,4],                  biasedRatioP = coef(summary(lm(biasedY/biasedX ~ 1)))[1,4]         ))         })          }          pMat <- as.data.frame(t(simBiasPower()))          library(plyr)     powerBySim <- colwise(function(acol) sum(acol<=0.05)/100)(pMat)          powerBySim      This produces the following estimates of power for the four analyses:                 unbiasedSlopeP unbiasedRatioP biasedSlopeP biasedRatioP     1              1           0.66            1         0.16      Note, I've tried this also forcing the intercept through 0 and thus not estimating an intercept for y~x, and the answer doesn't change. I've also changed the distributional form of the unbiased x, and it doesn't seem to matter very much. Similar drop in power. Of of the top of my head, when one divides by x, you reduce the leverage of points at the extreme end of the relationship - along the lines of reducing your variability in x in a regression context. But it feels like there is something more forehead-smackingly basic going on here.


title: K-L divergence is 0 for clearly different distributions. Why?

text: I have two datasets whose densities look clearly different. However, when I use the KLdiv function, it gives me a value of 0 and I don't understand why. Any suggestions on where I am going wrong? Datasets are at the end of this post. I did something like this:               library(flexmix)     KLdiv(cbind(density(A)$y, density(D)$y))      And it yields me:                    [,1] [,2]     [1,]    0    0     [2,]    0    0      JS divergence is zero too as a result of this:               0.5*(KLdiv(cbind(density(A)$y,0.5*(density(A)$y + density(D)$y)))) + 0.5*(KLdiv(cbind(density(D)$y,0.5*(density(A)$y + density(D)$y))))               [,1] [,2]     [1,]    0    0     [2,]    0    0      However, the pdfs of A and D are quite different: ![enter image description here](http://i.stack.imgur.com/K0ks8.png) Following are my two datasets:               A = c(53, 57, 54, 63, 66, 66, 75, 86, 77, 52494, 56722, 57097, 56754,      57399, 58145, 56582, 58429, 57868, 57349, 57569, 56507, 55925,      58177, 56455, 58162, 58253, 58340, 55172, 54767, 57750, 515,      510, 1686, 299, 571, 301, 300, 305, 300, 429, 975, 1364, 1377,      12122, 1255, 15845, 23767, 6883, 16844, 13867, 0, 0, 0, 8293,      6119, 7504, 6292, 6257, 6406, 1267601, 427055, 667541, 1209979,      760919, 1035907, 710933, 2217471, 725223, 888272, 1228833, 2153679,      1541844, 1310682, 1706112, 901743, 955793, 1900751, 54844, 2125783,      1109668, 2533008, 1476430, 3838432, 1594130, 889522, 884048,      1212327, 1407082, 2529584, 1410899, 1154139, 1216456, 1255570,      1550869, 2096869, 2360463, 2079124, 2239979, 5206460, 1350783,      2866784, 3487824, 610215, 2270840, 933909, 1049955, 901239, 1671539,      12292475, 2160519, 1730576, 1307651, 1659493, 1419815, 2333420,      2813600, 6260, 1557051, 2012477, 2383571, 1889906, 2497337, 3786140,      3233662, 4550030, 4765371, 4431086, 2143366, 1758742, 2520050,      2036522, 2945222, 4420937, 1452987, 4455005, 1838441, 1712123,      1213652, 1058007, 1156038, 11818454, 2913849, 5614698, 4188568,      4121445, 4067988, 5157146, 3858549, 3799965, 4465718, 5510926,      1426606, 4150290, 3995099, 5051489, 5518465, 6891987, 2599328,      2255332, 2921442, 2189249, 3552345, 4020428, 3752018, 4649134,      2982890, 3035131, 3337928, 7469273, 5306854, 3343986, 3690062,      4301378, 3280312, 3951362, 3074957, 3885109, 3965623, 3411739,      3275028, 2955563, 4281609, 3007306, 4932588, 1215474, 5038889,      3879138, 5686709, 99190451, 3857811, 21850538, 9272045, 23873337,      80095039, 99642585, 101780491, 60067705, 3546841, 2329745)          D = c(87983, 1130660, 721996, 881447, 0, 0, 56561, 148153, 82070,      88294, 81592, 303643, 1655103, 1660248, 1172794, 314841, 1585898,      339034, 464871, 1620502, 569710, 450522, 1152443, 1012775, 278482,      2019150, 8263736, 1913965, 1573202, 10456294, 13028222, 17432764,      19370911, 3683766, 12794272, 20638951, 14958438, 21847386, 33632736,      31652260, 23939674, 1432818, 37414993, 33164198, 34101979, 36387746,      41604477, 49173987, 39454209, 38942401, 44779660, 39630652, 37190054,      36059029, 45519990, 36309731, 44287159, 52460663, 45415958, 46605375,      70220155, 40907676, 0, 81266347, 71459813, 76837711, 112575185,      93436064, 109946521, 93868840, 122776950, 97212417, 108955821,      145587560, 216752084, 336859809, 185175516, 291129509, 349552593,      372988145, 292428463, 415094255, 185428625, 360299946, 387408971,      335578001, 393822458, 364598154, 416030537, 418917188, 448585165,      418281134, 380423194, 359740594, 510471342, 430134138, 463678777,      516374163, 512168586, 448045870, 491493753, 531045072, 499352572,      430560185, 507585264, 453550205, 436303557, 373408445, 169499536,      311790281, 468707162, 433377700, 335570055, 434688343, 358754832,      413128178, 358604964, 333417416, 103901186, 45295220, 107528512,      38698395, 68062269, 71930943, 63461445, 80691308, 84880325, 116910124,      78199074, 33932505, 98104155, 63883422, 2449889, 81687886, 112591823,      933357, 110882829, 167759933, 143331631, 128723052, 97515988,      89474519, 95248435, 77201398, 115732319, 150051651, 89815946,      119109327, 156889426, 104873286, 112453255, 103860299, 141042327,      161372922, 228847332, 279230413, 227022713, 245884032, 449327187,      399725637, 427264302, 385944417, 403338926, 384747124, 413700052,      373033651, 399978293, 455361833, 417924530, 473795493, 485446783,      551847677, 536782686, 551399627, 433065731, 434131355, 501258800,      468106521, 430447721, 408847961, 373866095, 360367708, 415444352,      439868563, 469183737, 572553168, 525927688, 417114578, 463202610,      335993664, 323369810, 319030248, 416164672, 393022155, 324094628,      324588328, 225682090, 218461312, 210500238, 204682022, 219859789,      244008146, 279700073, 231405316, 397901380, 351618973, 347947034,      475587874, 385887638, 441784135, 453616386, 425914531, 384140372,      394301503, 369180956, 359728515, 319432242, 340549134, 155250075,      390714221, 514657336, 357766885, 250607790, 313543688, 289122479,      289912982, 221289620, 212056585, 200354558, 172626254, 152561904,      141737182, 151121151, 135049494, 312876338, 155618873, 199563260,      243044037, 273332195, 302758331, 340552788, 340139019, 435174604,      467950978, 430091965, 345051079, 351816570, 263430776, 255793594,      164134822, 173371250, 254249143, 313999193, 148044484, 162452217,      130834412, 167723965, 200147738, 100376398, 103186023, 217542202,      302451572, 180647415, 262106015, 267346637, 186793712, 239398890,      312860449, 267037929, 159408789, 200059134, 132928248, 252577685,      255926249, 229046165, 113443872)


title: Non-parametric confidence interval for a small sample

text: I have a question about measuring a confidence interval for a binomial success rate using from a "small" universe. I should say that I don't think that the underlying distribution is binomial for reasons I will explain below. We rate success and failure of projects. Projects can either succeed or failure. Every year we review between 8 and 20 projects. We review 100% of the projects that are completed. We want to be able compare the success rate (for example 70% of projects were successful in year X) year to year and identify significant differences. I do not believe that the projects are independent-- if we carry out too many projects, it will lower the quality of each project. Likewise, I do not think that there is an underlying mean that drives the process; conditions vary from year to year. Since I can't use a binomial estimate, is there some sort of non-parametric method to identify the confidence interval? Is the sample too small? If we had a bigger sample (say 150), would it work better? Or perhaps with even with my earlier caveats, I should assume that the distribution is binomial? It also occurs to me that perhaps we could do a simple probit or logit on each observation (project) and use that to identify systematic differences from year to year. This question has been driving me crazy for the last two weeks. Any suggestions?


title: Troublesome residual plot from linear mixed model

text: I have fitted the following linear mixed model based on the results of an economic game: `lmer(TotalScore~perOOgivenP+Game+(1|Subject),REML=T,data=mdl1table)->m1` `TotalScore` is a integer.   `perOOgivenP` is a proportion between 0 and 1 (most of which are 0).   Game is numeric and tells us if it is the 1st or 2nd game played by the participant The QQ plot looks good so I am confident the residuals are normally distributed. The fitted vs. residuals plot does not look so good. To me it looks like the residuals are biased. I am not sure what Is causing this or what to do about it? Could it be caused by the amount of 0's in `perOOgivenP` (32/46 data points). `perOOgivenP` is the proportion of times a particular behaviour was made. Would anyone suggest making this binary as in 0 or 1 (actual value not 0). perOOgivenP is as follows:                [1] 0.500000 1.000000 0.333333 0.000000 0.000000 0.000000      NaN 0.166667 0.250000 0.800000     [11] 0.166667 0.000000 0.333333 0.000000 0.000000      NaN      NaN 0.000000 0.000000 0.000000     [21]      NaN 1.000000      NaN 0.000000      NaN 0.000000 0.000000      NaN      NaN      NaN     [31]      NaN 0.000000      NaN 0.000000 0.000000 0.000000      NaN      NaN      NaN 1.000000     [41]      NaN 0.000000 0.411765 0.000000 0.000000      NaN      NaN 0.000000 0.200000 0.000000     [51]      NaN 0.000000 0.333333 0.000000 0.250000      NaN 0.000000      NaN      NaN      NaN     [61] 0.000000      NaN 0.000000 0.000000      NaN 0.000000      NaN      NaN      NaN 0.000000     [71] 0.000000 0.000000 0.000000      NaN      NaN      NaN      NaN      NaN      NaN      NaN      ![Fitted vs. residuals plot](http://i48.tinypic.com/qohyc6.jpg)


title: Estimate global mean and stdev from samples taken at regular intervals

text: I have data from a load test of an automated system with several thousand data points that cover roughly 1 week of operation. I need to compare several algorithms for this system to see which is faster, if any (the metric used is the response time of the algorithm in milliseconds). Since it's technically impossible to store the response time of each operation of the algorithm, the simulator computes the mean, stdev and sample count of all the operations performed during the last minute, in order to reduce the amount of data generated. **Question 1:** Is it possible (or wise) to estimate the global mean/stdev based on the mean/stdev of each sample? Since each sample has a different sample count, I assume that some kind of weighting should be used for this estimation, but I'm not sure how to do this robustly, and I don't want to inadvertently alter the analysis. **Question 2:** A coworker suggested that computing the global quartiles (Q1, Q2, Q3) from the samples might be a more sensible approach to convey the general behavior of each algorithm. Which method should be preferred? I'd like to ensure that the analysis is statistically robust and that I'm not seeing "things that aren't there". **UPDATE 1:** I've found out that the global mean can be estimated by computing a weighted mean of the means of the samples and the global stdev can be obtained using the pooled variance estimate. Is this correct? Nevertheless, I still don't know if it is possible to do the same with the median as I haven't found much literature regarding weighted medians. I still don't know which measure is more appropriate for my data. **UPDATE 2:** I've learnt that Quartiles should be preferred to mean and standard deviation when data does not follow a normal distribution, although I don't know how to determine if my data follows it or not. I've found out some methods to compute the weighted median that I assume can be applied to any quartile but ¿does it make any sense to compute the weighted median from the collected means?


title: Checking normality of residuals in 2x2 mixed methods ANOVA

text: I am running a 2x2 mixed method ANOVA (with one between-group factor and one within-group factor - time 1 and time 2). I want to check whether the assumption of normality of residuals is met. When I run the ANOVA in SPSS I get two sets of residuals if I save the unstandardised residuals - one for time 1 and one for time 2. I was wondering why there are two sets of residuals, and which set of residuals should be normally distributed (or is it both?).


title: How to setup and interpret ANOVA contrasts with the car package in R?

text: Let's say I have a simple 2x2 factorial experiment that I want to do ANOVA on. Like this, for example:               d   <- data.frame(a=factor(sample(c('a1','a2'), 100, rep=T)),                       b=factor(sample(c('b1','b2'), 100, rep=T)));     d$y <- as.numeric(d$a)*rnorm(100, mean=.75, sd=1) +            as.numeric(d$b)*rnorm(100, mean=1.2, sd=1) +            as.numeric(d$a)*as.numeric(d$b)*rnorm(100, mean=.5, sd=1) +            rnorm(100);        1. In the absence of a significant interaction, by default (i.e. `contr.treatment`) the output of `Anova()` is the overall significance of `a` over all levels of `b` and of `b` over all levels of `a`, is that right?   2. How should I specify a contrast that would allow me to test the significance of effect `a` with `b` being held constant at level b1, of effect `a` with `b` being held constant at level b2, and of the interaction `a:b`?


title: Are there any ways to glean useful analysis information from a comments field?

text: I have a dataset of data and two of the fields are user comments - (I did not design this form) I would like to know if there is anything I can do to get any information from a comments field - are there any best practices out there to make use of this type of data in a report / analysis ? Otherwise is this data useless?


title: How to get equivalent ANOVA in R for a model with and without a random effect?

text: My question is: Is there a way to either force `Anova()` to somehow analyze `gls` objects (which internally are almost identical to `lme` objects), or force `Anova()` to honor the `test.statistics='F'` argument, or at least do a valid type-II sum of squares by hand on an `lme` and a `gls` object? Why: I'm trying to get Anova output in the same format for an `lm` or `aov` or `gls` object and an `lme` object that uses the same fixed effects formula but in addition has random effects. If I use `Anova()` from the `car` package, I get F-statistics for `aov` and `lm` objects but Chi-square statistics for `lme` objects, and it doesn't work at all for `gls` objects [1]. If I use `anova.gls()` and `anova.lme()` then they do both return F-statistics, but they use type-III or type-I sum of squares and I'm trying to use type-II. Thank you kindly. [1]: Gives error `Error in eval(expr, envir, enclos) : object 'y' not found` where y is the response variable... this can be traced to the attribute for `model.matrix()` for `gls` objects failing to have an `assign` attribute.


title: Significance test of whether participants have more friends of the same sex

text: I've run into a simple problem - but no idea how to access it correctly.   I've 85 people asked concerning their social network; first example: how many friends do you have. Second, what gender is each of this friend. Then I do a table which displays the average number of female friends and the average number of male frieds of each interviewed person: separated for the male and female interviews and for all together. I got the following table:   $$ \begin{array} {rrrrr} &&\text{avg no} &&\text{answers}& & \text{Sdev in}& \\ &&\text{of friends}& & & &\text{no of friends}& \\ & & m &f& m& f& m& f \\ \text{Sex of Interviewed}& M &5.57 &4.61& 54& 54& 3.543& 2.609 \\ & F &4.84 &6.42& 31& 31& 2.734& 3.264 \\ & All& 5.31& 5.27& 85& 85 &3.273& 2.978 \end{array} $$ and I see, that for the male interviewed ("M") the avg number of male friends is higher than the avg number of female friends, and for the female interviewed ("F") it is oppositely. Now what test for significane of the difference were correct here?   I could t-test for the means-difference in the M and F-group separately, but this seems a loss of infomation to me. Or the table suggests something like a chi-square; but how should I apply it here? [update]   Another approach is to determine the percent of male friends for each respondents, and then determine the average of this percentages for male and female respondents separately: $$ \begin{array} {cccc} &&\text{avg "% of }&N&\text{sdev of} & \\ &&\text{friends are male"}& &\text{"%..."} & \\ \text{sex of respondent}&M&52.746&54&21.087& \\ &W&40.868&31&17.235& \\ &all&48.414&85&20.487& \\ &&&&& \\ &f=7.101&&&& \\ &sig=0.009&&&& \\ \end{array} $$ The comparision of that means is significant due to the f-test - but is this a better sensical approach? [update 2]   This another idea to use the chisquare-rationale. I (re-)expand the averages to the sums: "sum of male/female contacts per respondent" and compute the chisqare based on the indifference-table. $ \qquad \small \begin{array} {c | cc |c} \text{Sum} &\text{m} &\text{f} &\text{all} & \\ \hline \text{M} &301&249&550& \\ \text{F} &150&199&349& \\ \hline \text{all} &451&448&899& \\ \\ \\ \text{Indifference} &\text{m} &\text{f} & \\ \text{M} &275.92&274.08& \\ \text{F} &175.08&173.92& \\ \\ \\ \text{Residual} &\text{m} &\text{f} & \\ \text{M} &25.08&-25.08& \\ \text{F} &-25.08&25.08& \\ \\ \\ \text{Chisq} &\text{m} &\text{f} & \\ \text{M} &2.28&2.3& \\ \text{F} &3.59&3.62& \\ \end{array} $ $ \qquad \chi^2 =11.79$ On the other hand, here -I feel- is the chisquare "inflated" because we have such a big N (which is actually an overall sum). Then the significance should be considered critically. Then gaian - what is the most sensical one? [update 3]   Here I show a table using an "homophily"-index: 0 means completely heterophil, 1 means completely homophil (in terms of same sex between respondent and his reported friends - requires at least one response/one friend per respondent) $ \qquad \begin{array} {c|cc|c} &&\text{avg of hom} &\text{sdev} &\text{semean} &\text{N} & \\ \hline \text{sex of respondent} &\text{M} &0.53&0.21&0.03&54& \\ &\text{F} &0.58&0.16&0.03&30& \\ \hline &\text{All} &0.55&0.19&0.02&84& \\ &&&&&& \\ &\text{f=1.297} &&&&& \\ &\text{sig=0.258} &&&&& \\ \end{array} $   I've got another test-value f and another significance level; well here I ask, whether male and female respondents are differently homo/heterophil which is another question than before. However, it is more precisely focused to an interesting indicator. The semean shows, that (only) female respondents seem to deviate significantly (5%-level) from indifference (which means hom=0.5) It might, anyway, have a little drawback in that the index for each respondent is based on another number of responses and thus has a more or less reliable value for each of that respondents. But this seems to be a too sophisticated problem here, so I think, I'll stay with that type of measuring.   Thanks so far to all respondents here!


title: Vector error correction model

text: I am working on finding the relationship between equity, gold, crude oil, and currency value. I used the Johanson cointegration test. Trace test indicates 4 cointegrating equations at the 0.05 level. For using VECM the rank for cointegration in Eviews has a maximum of 3. What is the problem in this model? How do you do the analysis when both number of variables and cointegration are same?


title: R fitting Poisson distribution with weighting

text: I am wanting to fit a Poisson distribution but I am wanting to weight some of the values as more important. Is there a way to do this in R with `fitdistr()` or some equivalent function? Here is what I have at the moment:               randoms <- rpois(15,10)     weighting <- seq(1, 100, by=1)     fit <- fitdistr(randoms, "poisson")      I would like to use the 'weighting' data as a way to emphasize the observation at the end of the data set.


title: Inverse probability of treatment weighted (IPTW) estimator for a binary outcome

text: Recently, there are several estimators have been proposed to estimate the average treatment effect (ATE) in observation studies, such as IPTW, doubly- robust estimator, etc. It fully makes senses to me for using these estimators to estimate ATE when outcome variable is continuous. However, it does not make any sense to me at all to use these estimators for a binary outcome. There are many researchers have published their papers in statistical journals such as "statistics in medicine" to demonstrate how to use these estimators to estimate ATE (risk difference, OR, RR) when outcome variable is "dichotomous /binary" variable. The IPTW estimator is $\hat{\Delta} = \frac{1}{N}\{\sum_{i=1}^{N}\frac{Z_iY_i}{\hat{e}_i} - \sum_{i=1}^{N}\frac{(1-Z_i)Y_i}{1-\hat{e}_i}\}$, where $N$ is sample size, $\hat{e}_i$ is an estimated weight, $Z_i$ is the treatment assignment (control/trt or 0/1) for $i^{th}$ subject and $Y_i$ is the outcome (0/1 or No/Yes) of $i^{th}$ subject. My question is that $Z_i$ and $Y_i$ are categorical variables and $\hat{e}_i$ is a continuous variable, which is between 0 and 1. How can categorical variables be divided by a continuous variable? many thanks in advance


title: How can I validate a logistic regression model using averaged parameter estimates?

text: Let me say thanks in advance. I'm working with a set of data that contains reported coyote sightings. I use 2/3 of the data for model calibration along with an equal number of pseudo absences. I developed all possible models and ranked them according to their AIC weights. I chose the top models, who's weights summed to 75%, and created averaged model estimates for each of my parameter. Now I would like to test the accuracy of my model using the 1/3 of the data I held out. I assume I will need another equal amount of pseudo absences to include in the validation. My problem is that I have averaged parameter estimates. All the software packages that I have come across have settings where you train your model using a specific data set and test it against a specific data set. However, I want to validate my model using averaged parameters, rather than the parameters calculated from running a simple regression. I'm assuming I need to transform my outputs into percentages and choose an error threshold. Then I can produce a simple omission/commission report. But as this is for my Master's thesis, I want to be sure I'm using a respected and widely used method. Can anyone point me in the right direction? My advisers are used to traditional validation methods so I'm doing a bit of outsourcing. Again, thanks for any and all advice.


title: Modeling vacancy rate

text: I have 100 geographical regions in a country. For each region the total number of houses and the number of vacant houses have been collected yearly over 20 years. I have also some other economic indicators at the country level (GDP, interest rate etc.). Now, given the forecasts for these indicators for next year, I want to forecast next year vacancy rate. I have first used an auto-regressive mixed-effect model in `R` (package lme4) where the vacancy rate (computed as the ratio of vacant houses over the total number of houses) in a region depends on the last year's vacancy rate, the mean vacancy rate of neighboring regions, GDP and interest rate. The problem with this model is that the vacancy rate can go outside the range [0,1], which obviously does not make sense. I need to restrict the range of the vacancy rate: a simple fix is restricting _ex post_. Does anybody have experience with such models? I think that I can use some mixed multinomial logit model probably. I would appreciate if you could provide guidance along with some `R` code. Regards


title: Correlation between two Decks of cards?

text: I've written a program to simulate an overhand card shuffle. Each card is numbered, with suit going from `CLUBS, DIAMONDS, HEARTS, SPADES` and the rank from Two up to Ten then Jack, Queen, King and Ace. Thus the Two of Clubs has a Number of 1, the Three of Clubs a 2 ....Ace of Clubs is 13... Ace of Spades is 52. One of the methods to determine how shuffled the cards are is to compare it to an unshuffled card and see if the order of cards is correlated. That is, I might have these of cards, with the unshuffled card for comparison:               Unshuffled          Shuffled            Unshuffled number   Shuffled number     Two of Clubs        Three of Clubs      1                   2     Three of Clubs      Two of Clubs        2                   1     Four of Clubs       Five of Clubs       3                   4     Five of Clubs       Four of Clubs       4                   3      The correlation by the Pearson method would be: 0.6 With a large set of cards (all 52) you might see patterns emerge. My Hypothesis is that after more shuffling you'll get less correlation. However, there are lots of ways to measure correlation. I've tried my hand at Pearson's correlation but I'm not sure if this is the right correlation to use in this situation. Is this a suitable correlation measure? Is there a more suited measure? **Bonus Points** I sometimes see this sort of data in my results: ![Sample Card Correlation](http://i.stack.imgur.com/KcJh1.png) Clearly there is some correlation but I don't know how you measure the separate 'trendlines'?


title: Making feature vector from Gabor filters for classification

text: My aim is to classify types of cars (Sedans,SUV,Hatchbacks) and earlier I was using corner features for classification but it didn't work out very well so now I am trying Gabor features.   code from here   Now the features are extracted and suppose when I give an image as input then for 5 scales and 8 orientations I get 2 [1x40] matrices. **1. 40 columns of squared Energy.** **2. 40 colums of mean Amplitude.** Problem is I want to use these two matrices for classification and I have about 230 images of 3 classes (SUV,sedan,hatchback). I do not know how to create a [N x 230] matrix which can be taken as vInputs by the neural netowrk in matlab.(where N be the total features of one image). My question:   1. How to create a one dimensional image vector from the 2 [1x40] matrices for one image.(should I append the mean Amplitude to square energy matrix to get a [1x80] matrix or something else?)   2. Should I be using these gabor features for my purpose of classification in first place? if not then what? Thanks in advance


title: Find and fill values for the missing points using close/similar vectors

text: I have a matrix which looks like the attached image below, and I would like to know what's the statistically best educated guess for the missing steps (filled in orange) for the vector $\mathbf N$. Please note that:   * The values are _binary_ : they are always `0` or `1`.   * Underneath the `0`-columns, we're supposed to have _a shape_.   * Vector N+-1 is more similar to N than N+-2. ![http://f.cl.ly/items/2d3A0s3r0s442E3K312i/data.png](http://i.stack.imgur.com/BLMLK.png)


title: Making sense of the output of svd() in R

text: I created a matrix of dimensions 6X5 and applied the svd() function on it.               >matrix_A       V W X Y Z     A 4 1 5 0 0     B 0 1 3 0 1     C 0 1 0 5 4     D 0 1 0 4 0     E 0 1 2 0 4     F 5 0 4 1 0          >svd_A <- svd(matrix_A)      When I see the 'U' matrix, it has dimensions 6X5.               >svd_A$u                [,1]        [,2]        [,3]       [,4]       [,5]     [1,] -0.6337425  0.27376062 -0.01574796 -0.2086500 -0.5249861     [2,] -0.2648734 -0.01058103  0.33757481 -0.6686728  0.6057323     [3,] -0.2315021 -0.80402040 -0.06816000  0.2435956  0.1600463     [4,] -0.1117109 -0.40493588 -0.51738041 -0.4942192 -0.2799135     [5,] -0.2618584 -0.26685759  0.71943981  0.1595892 -0.3116502     [6,] -0.6273517  0.20810652 -0.30963419  0.4245959  0.3954564      Now, I have read that the 'U' matrix must by a square matrix of dimension 6X6 (here). Further, the singular value matrix ('D') must be of 6X5 and the 'V' matrix must be of 5X5 dimension. Now, what is the R's svd() doing? Why is my U matrix of dimension 6X5? Is there another way of factorizing the matrix other than the method I know? (there most probably is) Update: The manual says the dimension of U would be c(n,nu) where _nu_ is the number of left single value vectors to be computed between 0 and _n_ and the default _nu_ is *nrow(matrix_A)*. So, by default, in my above example, _n_ = _nu_ = 6. Hence it should have been a 6X6 matrix. I got this doubt because I saw an example Here and the V dimensions weren't matching with the output of svd(). R's output for a:                    [,1] [,2] [,3]     [1,]    3    1    1     [2,]   -1    3    1      was:               svd(a)     $d     [1] 3.464102 3.162278          $u                [,1]       [,2]     [1,] -0.7071068 -0.7071068     [2,] -0.7071068  0.7071068          $v                [,1]          [,2]     [1,] -0.4082483 -8.944272e-01     [2,] -0.8164966  4.472136e-01     [3,] -0.4082483  5.231953e-16      Here V has dimension 3X2 while in the example in the link, V has dimensions 3X3.


title: Why do PCA and Factor Analysis return different results in this example?

text: The following question is about an example from "The Elements of Statistical Learning" by Hastie, Friedman and Tibshirani $X_1 = Z_1 $ $X_2 = X_1 + 0.001 * Z_2 $ $X_3 = 10 * Z_3 $ Where $ Z_1, Z_2 $ and $ Z_3 $ are independent standard normal variables. According to Ex 14.15 of "The Elements of Statistical Learning": With PCA, "the leading principal component aligns itself in the maximal variance direction $X_3$" With Factor Analysis: "the leading factor essentially ignores the uncorrelated component $X_3$, and picks up the correlated component $X_2 + X_1$" Why? I thought that they are both "powered by" the same matrix decomposition? What have I missed?


title: What is the intuition behind beta distribution?

text: Disclaimer: I'm not a statistician but a software engineer. Most of my knowledge in statistics comes from self-education, thus I still have many gaps in understanding concepts that may seem trivial for other people here. So I would be very thankful if answers included less specific terms and more explanation. Imagine that you are talking to your grandma :) I'm trying to grasp the **nature** of **beta distribution** - what it should be used for and how to interpret it in each case. If we were talking about, say, normal distribution, one could describe it as arrival time of a train: most frequently it arrives just in time, a bit less frequently it is 1 minute earlier or 1 minute late and very rarely it arrives with difference of 20 minutes from the mean. Uniform distribution describes, in particular, chance of each ticket in lottery. Binomial distribution may be described with coin flips and so on. But is there such **intuitive explanation** of **beta distribution**? Let's say, $\alpha=.99$ and $\beta=.5$. Beta distribution $B(\alpha, \beta)$ in this case looks like this (generated in R): ![enter image description here](http://i.stack.imgur.com/7YgUM.png) But what does it actually mean? Y-axis is obviously a probability density, but what is on the X-axis? I would highly appreciate any explanation, either with this example or any other.


Is it feasible to use t-SNE to reduce a dataset to one dimension? Suppose that I have a matrix, $X$, can I reduce it to a column vector, $Y$ with t-SNE? Suppose that $X$ has 100 columns, how much information can I expect to lose by reducing it to just a column vector with t-SNE (if possible)? References: L.J.P. van der Maaten and G.E. Hinton. Visualizing High- Dimensional Data Using t-SNE. Journal of Machine Learning Research 9(Nov):2579-2605, 2008. • L.J.P. van der Maaten. Learning a Parametric Embedding by Preserving Local Structure. In Proceedings of the Twelfth International Conference on Artificial Intelligence and Statistics (AI-STATS), JMLR W&CP 5:384-391, 2009. • L.J.P. van der Maaten. Barnes-Hut-SNE. In Proceedings of the International Conference on Learning Representations, 2013. All text available at: http://homepage.tudelft.nl/19j49/t-SNE.html


title: lm regression in R - inconsistent R² for intercept-free model?

text: I am carrying a linear regression on some data. One of my variables is a factor (categorical). Using regression with an intercept leads to difficult interpretation, since one of the factor levels is taken as the intercept, and the remaining levels are given relative to that. Removing the intercept give me an effect of each level of the factor, which I what I want. As far as I know, both models are precisely equivalent. They produce identical predictions (on the training set) - to within ~3e-15. However, their R² scores vary wildly.               # MWE     library(car)     int <- lm(fscore ~ 1 + partner.status + conformity + fcategory,               data = Moore)  #with intercept     nint <- lm(fscore ~ 0 + partner.status + conformity + fcategory,                data = Moore) #w/o intercept     summary(int)$r.squared     summary(nint)$r.squared  #R² values are not remotely the same     max(predict(int)-predict(nint)) #Predictions are essentially identical      Why are the models not identical? Is it because R² is a comparision of the model to "no model", and that "no model" corresponds to "y=0" and "y=mean(fscore)", for nint and int, respectively?


title: Jags error with dgamma

text: I am trying to sample from a Gamma distribution in JAGS               gd[i] ~ dgamma(k,r)      where **k** and **r** are have priors such that they are positive. However I get the following error when I compile the model:                 Error in node gd[9]     Unobserved node inconsistent with unobserved parents at initialization      Any ideas? **EDIT: BACKGROUND** What I am trying to do is to estimate the parameters of the CIR process with Bayesian techniques. In order to assess the estimator, I have simulated some data with given parameters on a quarterly basis (timeStep=90):               jdata <- list(timeSerie=c(35000, 47123.3999346341, 55373.904244638, 71755.7487795017, 91246.3370619453, 76969.2358694695, 83390.030308266, 80657.1222156083, 76784.466470129, 83962.4642428193, 85335.5796326587, 87018.8279086056, 74059.9912730579, 69683.5404626939, 90356.0122717258, 87240.7793078735, 74276.2931903815, 39899.5948334073, 37911.738710607, 73811.7536157058, 69334.5390010963, 51326.297331755, 42711.7358571988, 32409.9939465308, 16663.7353817099, 22421.1721803643, 38604.0551895535, 50181.3269231991, 44346.3517857709, 52127.2156791002, 54740.1181745155, 41869.456839522, 50665.5626313423, 51121.4339441224, 39059.432945341, 44350.0205383802, 68629.3291332604, 49933.2126408239, 64358.5131164748, 102975.634845413, 84658.5787221384, 52243.061086781, 56988.7184089767, 43548.6558054764, 38536.145535317, 43115.4155553003, 50328.6208374212, 47191.4284861612, 42603.6038389972, 43749.2977769786, 43583.8022320709, 36672.3782472362, 26377.9033658583, 41746.6472558272, 26086.7868276075, 47228.5390765096, 75765.4757503756, 70939.107521347, 63949.8857556708, 55153.6701237129),timeStep=90)      Then I wrote a model with JAGS. The problem is that the CIR process has scaled noncentral chi square distribution and this is not directly supported in JAGS. So inspired by `R` implementation, I came up with               model {         for (i in 2:length(timeSerie)){             ncp[i] <- 2*c*et*timeSerie[i-1]             pd[i] ~ dpois(ncp[i]/2)             cd[i] ~ dgamma(pd[i]+1E-10,c)             gd[i] ~ dgamma(df/2,4*c)             timeSerie[i] ~ dsum(cd[i],gd[i])         }              #parameters         et <- exp(-kappa*timeStep)         c <- 2*kappa/(sigma^2*(1-et))         df <- 4*kappa*mu/sigma^2              #CIR parameter priors          lkappa ~ dunif(-10,-2)         kappa <- exp(lkappa)         mu ~ dunif(30000,100000)         lsigma ~ dunif(0,5)          sigma <- exp(lsigma)     }      Trying to run the model (`jags.model` in `R`), I get the following error                 Error in node gd[18]     Unobserved node inconsistent with unobserved parents at initialization      (Now it is 18 instead of 8, because I use other simulated data). I tried to debug JAGS (gdb and valgrind) to no avail: I don't have enough knowledge to do it...


title: t-test that means are the same?

text: I think this is simple, and I am wrong-thinking it.... t-tests are used to test the null hypothesis that the means are equal; $H_0 : \beta_1 = \beta_2$, $H_1 : \beta_1 \neq \beta_2$ But I want to test the other way around; $H_0 : \beta_1 \neq \beta_2$, $H_1 : \beta_1 = \beta_2$ Do I just swap the hypotheses and 1-p the test result? or is there something subtle that I am missing?


title: How to analyse pre-post quasi-experimental study with an additional baseline measure?

text: I am having a hard time with what seems to be a very basic research design issue. The study has 3 variables:   * `Self-control` and `adjustment` are the main ones, since they were measured in pre and post tests, before and after a year long social and emotional program in two different high schools. (I guess this part of the design would be quasi-experimental)   * The third variable is `values in interpersonal relationships`. This was measured only one time at the beginning of the program since it tends to be more permanent. Also, the program wasn't designed to have an effect on this, except values can also affect the adjustment and self control of the students. The results show the groups tend to the individualist values vs. collectivistic ones (descriptive design).  I wanted to add hypothesis about this last variable to my paper, but a teacher told me it would be better not to include those, since those are more descriptive hypotheses.   * What is the best and correct way to include all this information on the paper, so it´s statistically sound?


title: How can I speed up calculation of the fixed effects in a GLMM?

text: I'm doing a simulation study which requires bootstrapping estimates obtained from a generalized linear mixed model (actually, the product of two estimates for fixed effects, one from a GLMM and one from an LMM). To do the study well would require about 1000 simulations with 1000 or 1500 bootstrap replications each time. This takes a significant amount of time on my computer (many days). `How can I speed up the computation of these fixed effects?` To be more specific, I have subjects who are measured repeatedly in three ways, giving rise to variables X, M, and Y, where X and M are continuous and Y is binary. We have two regression equations $$M=\alpha_0+\alpha_1X+\epsilon_1$$ $$Y^*=\beta_0+\beta_1X+\beta_2M+\epsilon_2$$ where Y$^*$ is the underlying latent continuous variable for $Y$ and the errors are not iid.   The statistic we want to bootstrap is $\alpha_1\beta_2$. Thus, each bootstrap replication requires fitting an LMM and a GLMM. My R code is (using lme4)                   stat=function(dat){             a=fixef(lmer(M~X+(X|person),data=dat))["X"]             b=fixef(glmer(Y~X+M+(X+M|person),data=dat,family="binomial"))["M"]             return(a*b)         } I realize that I get the same estimate for $\alpha_1$ if I just fit it as a linear model, so that saves some time, but the same trick doesn't work for $\beta_2$. Do I just need to buy a faster computer? :)


title: Probability of drawing some number when randomly selecting subsets of size L

text: I am trying to find a way to solve the following problem : > I have some set _A_ , which contains a set of _N_ numbers, _X_ , where _Xi_ > is a number within _X_. The numbers in _X_ are unique and they are in sorted > order. > > I will draw _K_ random subsets of size _L_ from _X_. There are no > restrictions on how many subsets I can draw, the only restriction is that > the numbers in the subset are contiguous. Subsets are drawn with > replacement. > > For example, if my set was (1,2,3,4,5,6,7,8,9,10) then valid subsets of size > 3 is 1,2,3 or 8,9,10 or 3,4,5 and so on > > What is the probability of selecting _Xi_ in anyone of these _K_ subsets > given some subset size _L_? If the subset size were 1, then the problem is solved by the binomial distribution, or a poisson distribution when N is very large. Also, when the size of the subset =N, then everything has a 1/N probability of being selected. But I am not sure how to solve the problem when L is neither 1 nor N. This is because the positions are not independent because the subset represents contiguous numbers. Are there any suggestions for how to tackle this problem? Is it similar to other problems?


title: Does Cox Regression have an underlying Poisson distribution?

text: Our small team was having a discussion and got stuck. Does anyone know whether Cox regression has an underlying Poisson distribution. We had a debate that maybe Cox regression with constant time at risk will have similarities with Poisson regression with a robust variance. Any ideas?


According to this document: _The estimated logarithm of the hazard ratio is approximately normally distributed with variance (1/d1) + (1/d2), where d1 and d2 are the numbers of events in the two treatment groups._ Do you have a reference for this statement? Or at least can you say me which estimator is used? Thank you in advance, Marco


title: Online calculator for power analysis: what value to give?

text: Suppose, I have a correlation coefficient of 0.6234 between Variables A and B, on a sample size of 20. My alpha threshold is 0.05. I want to do a power analysis on this result. I found a calculator but it is asking for "Correlation p H1". Does this refer to the correlation coefficient "r" (0.6234 as above) or the p value of r (which I have calculated as 0.00330, two tailed).


title: How do I "vectorize" calls to dpois?

text: I am not sure if this is an instance of vectorizing the operations in R, but this is where I am stuck: I want to get:               dpois(1, 0.1)     dpois(2, 0.2)     dpois(3, 0.3)      and I tried:               dpois(1:3, 0.1:0.3)      and               do.call(dpois, list(x = 1:3, lambda = 0.1:0.3))      both do not work. It there a R-ish way of doing this?


title: Proper wording for lottery skewed by contribution size

text: I have an idea for a lottery that is skewed in a fairly simple way: The chance of winning is proportional to the square of each player's contribution, divided by the sum of the square of all players' contributions. It might be expressed this way:               Chance of winning = contribution^2 /                          (contribution1^2 + contribution2^2 + ... + contributionN^2)      I'm looking for a concise and accurate way of describing this, and any help is appreciated.


title: Interater reliability

text: What inter-rate reliability test is best for continuous data? I am doing a study with one variable with continuous data, now the measurement involves measurements done by two people. I would wish to do inter-rater reliability test for the data, so far I have collected a few samples and a sample data I have given below, what test would I use? Rater A 23.1 22.0 21.8 24.1 20.2 22.1 23.8 Rater B 23.0 21.5 22.0 23.9 19.8 21.9 22.9


title: Ratio scale explanatory variables in ordinal regression

text: I have a dependent variable which is ordinal, so based on logistic distribution. I need to run ordinal regression, but I have ratio scale independent variables as well as categorical: so how can I run this technique? Is it suitable? If yes, then, in which software?


title: Plotting Multiple Proportions With Standard Error

text: The Data: The observed probability (proportions) of three mutually exclusive events for five species. What is the best way to plot these data in R along with their standard errors? I'd like to avoid a "beside" bar plot with error bars (3 bars for each species). I was hoping to use a stacked bar plot, but I'm unsure how to incorporate standard error! Thanks!


title: How to choose a good operation point from precision recall curves?

text: Is there any standard method to determine an "optimal" operation point on a precision recall curve? (i.e., determining the point on the curve that offers a good trade-off between precision and recall) Thanks


title: Probability that one chi-square distributed variable is greater than another based on ratios

text: Supposing $\ X $ ~ $\chi_1^2 $ and $\ Y $ ~ $\chi_{10}^2$ and $\ X $ and $\ Y$ are independent. How would I calculate the probability that $\ X $ is bigger than $\ Y $? I know that the ratio of $\ X/1$ and $\ Y/10$ would be $\ F$ distributed, but I don't know if that's helpful or not. Can anyone help? Thanks


title: Regression when the value of one independent variable is dependent on another independent variable

text: I need to predict the number of retweeets a tweet receives as a function of (1) whether there is a hyperlink within the tweet's text and (2) the position of the hyperlink within the tweet. If a tweet doesn't contain a hyperlink then the position of the hyperlink within the tweet is a null value. Example data for this situation would be:               numberOfRetweets     hyperlinkInTweet?     positionOfHyperlink     (dependent var)    (binary, 0=no, 1=yes)   (integer, 1 to 135)               5                      0                    null          10                      1                     50          35                      1                     12           0                      0                    null           0                      1                      5           :                      :                      :      If a tweet doesn't contain a hyperlink (hyperlinkInTweet = 0) then it cannot have a position within the tweet text where the hyperlink begins (positions range from 1 to 135 inclusive). I've toyed with the idea of using zero in place of null, but I have reservations about doing this. Exploratory analysis of the data indicates that about half of the tweets contain a hyperlink. I'm wondering if this is a nested regression problem, but would be grateful for a second opinion on the best course of action for this problem. Many thanks.


title: How do I enter a continuous variable as a random effect in a linear mixed effects model?

text: I collected data on the growth of juvenile fish from 4 different types of crosses using multiple distinct family blocks and I am trying to see if cross type has an effect on growth using linear mixed effects models. I have one fixed factor (Length), 3 categorical random factors (sire, dam and sire/dam interaction) and a continuous random factor (density per tank). As I am relatively new to R and new to mixed effects models, I was wondering if I would have to code the continuous random factor in R differently from the categorical random factors using lme4? So far I have               model5=lmer(Length~(1|Dam)+(1|Sire)+(1|Sire:Dam)+(1|Density))     model6=lmer(Length~Cross+(1|Dam)+(1|Sire)+(1|Sire:Dam)+(1|Density))      is density coded in properly, or would I have to alter it due to the fact that it is a continuous variable?


title: How to choose window width for smoothing a periodogram?

text: I would like to firstly smooth a periodogram with Daniell kernel, and then test if it is zero at some frequency of interest, following Shumway and Stoffer's book Time Series with R: ![enter image description here](http://i.stack.imgur.com/9042w.png) I wonder how to choose the window width for Daniell kernel used to smooth the periodogram, for my time series? (e.g., in the quote, the half width is 9, and thus specify 2*9=18 the width, to calculate the critical value C for the test.) Thanks!


title: Basic knowledge about cognitive robotics

text: I am interested to learn about cognitive robotics and its application in situation assessment using natural language processing esp monitoring of human activities/disaster management.The processing would involve data visualization techniques. Can anyone suggest what is the starting point,which books to follow,the programming language/platform and some links to reference/example codes to begin with.Do I need to use ROS(robotic Operating system)?I came across GOLOG and situational calculus.Can anyone suggest where to download GOLOG & how to use it?Any other help would be highly solicited.


title: How to interpret model diagnostics graphics after R linear regression?

text: I am interested in understanding the graph plots we get after running `lm()` command (for linear regression) in R like, for example               lm.mod1 = lm(y ~ x1 + x2)      I then get the do the summary by:               summary(lm.mod1)      I get the result as:                Residuals:       Min      1Q  Median      3Q     Max      -750.32 -160.54  -49.83  115.83 2923.74           Coefficients:                                Estimate Std. Error     t value Pr(>|t|)         (Intercept)               -345.1552      37.0393   -9.319   <2e-16 ***              x1                52.9091       2.4929    21.224   <2e-16 ***              x2                8.9669        0.5395    16.620   <2e-16 ***          Residual standard error: 274.4 on 1985 degrees of freedom     Multiple R-squared: 0.2059, Adjusted R-squared: 0.2051      F-statistic: 257.3 on 2 and 1985 DF,  p-value: < 2.2e-16      I then do the plotting by               par(mfrow = c(2,2))     plot(lm.mod1)      I get 4 graphs (I can't post the graphs since I am a new user and my experience level is below 10. :/) My questions are :   1. How do they calculate F-statistics and t-value?   2. Could someone explain me the what do we interpret with the last two graphs i.e. $\text{Scale-Location vs. (Standardized residuals)}^{1/2}$ and $\text{Residuals vs. Leverage}$. What do you mean by Leverage?   3. What do you mean by Cook's Distance? I saw it on wikipedia but I didnt get it.    4. How could we suggest if our model is a good model or not?


title: t-test on two variance calculations

text: I have survey data that was recorded in a many different towns. I want to see for which variables is the data significantly more or less variable _within towns_ than _across all towns_. Currently, I'm thinking of conducting a Welch's t-test (for samples with unequal variance) on the difference between $\bar{X}_{1} =$ the mean of town- level variances and $\bar{X}_{2} =$ the variance of all observations. The t-test formula is: $$ t = \frac{\bar{X}_{1} - \bar{X}_{2} }{ \sqrt{\frac{s_{1}^{2}}{n_{1}} +\frac{s_{2}^{2}}{n_{2}} } } $$ Is this the right approach? I'm not sure if the t-test is right, and even if the t-test is right am I calculating the variance of variances correctly? I also wonder if there is a better way to approach this using ANOVA. In particular, I'm confused by the fact that in this test $\bar{X}_{i} = s_{i}^{2}$, which seems strange to me.


title: How to interpret the result of Fisher's unit root test

text: Following are the results from Fisher-type unit-root test for RDI (dependent variable). **How do you interpret it?**               Based on Phillips-Perron tests:          Ho: All panels contain unit roots           Number of panels       =    100     Ha: At least one panel is stationary        Avg. number of periods =  10.74          AR parameter:    Panel-specific             Asymptotics: T -> Infinity     Panel means:     Included     Time trend:      Included     Newey-West lags: 1 lag                                               Statistic               p-value          Inverse chi-squared(196)     P       207.1519                 0.2788     Inverse normal               Z        2.0005                  0.9773     Inverse logit t(389)         L*        0.5211                0.6987     Modified inv. chi-squared    Pm        0.5633                 0.2866          P statistic requires number of panels to be finite.     Other statistics are suitable for finite or infinite number of panels.


title: Churn predictive modelling

text: I'd like to build a predictive model for predicting churn for a website. Here is the information I have for each customer : What they did :   * visit the website    * Buy something   * Do not read thenewsletter   * Read thenewsletter and when :   * the date of the action For instance for a customer X :   * 90 days ago : visit the website    * 50 days ago : do not read the newsetter    * 20 days ago : do not read the newsletter The goal is to predict if a customer is gonna be inactive based on this rule : If the customer did not visit, buy or read our email on the last 12 months, he's lost. For now, I'm gonna use this data structure (example with the customer X) :               User_id | Last_action_1 | Last_action_2 |Last_action_3 |       X       Do not read       Do not read       Visit      The problem is that I lost the time information, I just retain the order. What is it recommanded to keep it ? I think about these ideas, but **I think there is better practices** :               User_id | Last_action_1 |Last_action_1_date |  Last_action_2 |Last_action_2_date| Last_action_3 |Last_action_3date          X       Do not read      20 days ago           Do not read       50 days ago      Visit         90 days ago


title: Longitudinal data using functional principal components mixed effect models

text: I have collected speech perception performance in children during one-year followed-up (four time points). A battery of speech perception measures (i.e., five measures) were used at each test interval. I want to examine the effects of three variables on speech perception performance (i.e., gender, chronological age, maternal educational level and cognitive abilities). To reduce the number of speech perception tests we are dealing with，the original speech perception scores will be reduced to a standardized speech perception score using the functional Principal Component Analysis. Then, the mixed- effect models will be used to examine the effects of three variables on speech perception performance during one-year followed-up. Is it feasible? If not, is there any other methods?


title: Band pass filter giving ‘wrong’ turning point

text: I am trying to run a Christiano FitzGerald band pass filter to estimate a long-run trend (with period in excess of 70 years). My data are the demeaned natural log of a commodity price index. My problem is that the asymmetric CF filter output does not look right around the start and end points. In particular my long run trend has a peak a few years before the end point, which I think shouldn’t be there (by reference to previous work and also just from looking at the data, which seem to be still trending up). ![cf filter output, looks wrong](http://i.stack.imgur.com/hqaiv.jpg)![hp filter, looks good but not standard practice](http://i.stack.imgur.com/yN8Qf.jpg) I have included R code below to replicate the issue in case that is useful. Any help would be very much appreciated! This is my first attempt at using this filter (this filter and the minimum period appear to be best practice for this kind of data).               library(mFilter)     library(ggplot2)     library(reshape2)     d<-c(1.17, 1.11, 0.68, 0.83, 0.98, 1.01, 0.76, 0.52, 0.34, 0.36,      0.18, 0.12, 0.07, -0.17, -0.47, -0.43, -0.01, -0.16, -0.51, -0.79,      -0.69, -0.76, -0.86, -0.85, -0.78, -0.8, -0.89, -0.95, -0.83,      -0.84, -0.96, -1.22, -1.41, -1.25, -0.94, -0.43, -0.58, -0.97,      -0.82, -0.47, -0.44, -0.56, -0.55, -0.32, -0.31, -0.46, -0.34,      -0.3, -0.25, -0.31, -0.57, -0.6, -0.42, -0.06, -0.3, -0.5, 0.12,      0.3, 0.33, 0.24, 0.49, -0.18, -0.06, -0.15, -0.11, 0.02, 0.14,      -0.3, -0.31, -0.28, -0.39, -0.94, -0.51, -0.74, -0.35, -0.37,      -0.28, -0.22, -0.24, -0.36, -0.37, -0.32, -0.39, -0.45, -0.46,      -0.48, -0.39, -0.21, 0.01, 0.02, 0.01, -0.06, -0.09, -0.04, -0.01,      0, -0.02, 0.03, 0, -0.03, -0.05, -0.06, -0.07, -0.08, -0.11,      -0.13, -0.15, -0.15, -0.18, -0.2, -0.24, -0.23, -0.25, -0.04,      0.86, 0.8, 0.8, 0.81, 0.73, 1.03, 1.36, 1.33, 1.18, 1.06, 0.98,      0.9, 0.22, 0.46, 0.25, 0.4, 0.59, 0.39, 0.33, 0.19, 0.11, 0.14,      0.27, 0.19, -0.16, 0.11, 0.5, 0.32, 0.31, 0.46, 0.72, 1, 1.12,      1.19, 1.47, 0.99, 1.23, 1.38, 1.35)     d<-data.frame(orig.series=d)     d$year<-as.numeric(row.names(d))         filter.output<-cffilter(d$orig.series,pl=70,pu=1500)     d$long.run.trend<-filter.output$cycle     d<-melt(d,id='year')     gp<-ggplot(d) + geom_line(aes(x=year,y=value,colour=variable))     gp


title: Does the regression line tell anything about mean

text: I am trying to find the association of BMI on the age at onset of hypertension in diabetes vs non-diabetes patients. Or in other words, does diabetes (if the patient is already diabetic before the hypertension diagnosis) modify the association of BMI and age at onset of hypertension. To find this, I have tried a multiple linear regression model with onset age of hypertension as dependent variable and BMI as independent variables adjusted for other covariates. I am getting a slope of -0.33 * BMI for non-diabetic hypertensives and -0.65 * BMI for diabetic hypertensives. i.e. slope is almost two fold for the diabetic hypertensives compared to non-diabetics. So I was concluding that the average rate of change in onset age of hypertension with respect to increase in each unit of BMI is more in diabetic patients than in non-diabetic patients. Now, counter intuitive to this result the mean onset age of diabetic hypertensive patients (55 years) is higher than non-diabetic hypertensives (46 years). Shouldn't the diabetic hypertensive patients have lower mean onset age than the non-diabetic hypertensives since the slope is steeper for diabetic hypertensives? Can anyone tell me if there is any flaw in this result and conclusion. Thanks much Here is the plot ![enter image description here](http://i.stack.imgur.com/A87f7.png)


title: Pooled regression vs non-pooled regression

text: I hope this is not a stupid question to ask. Let's say I have 5 years' data of 10 banks for 7 ratios; 1 ratio is dependent the other 6 is independent. If I ignore the time and bank effects in this cross-sectional data it is called "pooled data." Is this pooled data different from the not-pooled data in terms of results? Since I ignore the fixed and random effects and pool everything together, the result should be same with a regression result of simple linear regression. However, Stata has different menu options for pooled and nonpooled OLS regression, and the estimates are also not the same. I am confused now.


title: Variable Importance for different methods

text: Is it possible to compare variable importance, obtained for different models, built on different datasets with simular variables, using caret package? The property are different, but close, and the variables used for modelling are the same. For example, one can calculate importance for pls and random forest model, then scale them to obtain values from 0 to 1. Would it be valid to compare such values for two different Y? Is there any way to compare two different statistical models, built on two different datasets in terms of variable influence?


title: How to draw neat polygons around scatterplot regions in ggplot2

text: **How do I add a neat polygon around a group of points on a scatterplot?** I am using ggplot2 but am disappointed with the results of `geom_polygon`. The dataset is over there, as a tab-delimited text file. The graph below shows two measures of attitudes towards health and unemployment in a bunch of countries: ![scatterplot with density2d](http://i.stack.imgur.com/1lzPT.png) I would like to switch from `geom_density2d` to the less fancy but empirically more correct `geom_polygon`. The result on unsorted data is unhelpful: ![enter image description here](http://i.stack.imgur.com/WX148.png) How do I draw 'neat' polygons that behave as contour paths around the min-max y-x values? I tried sorting the data to no avail. **Code:**               print(fig2 <- ggplot(d, aes(man, eff, colour=issue, fill=issue)) +      geom_point() + geom_density2d(alpha=.5) + labs(x = "Efficiency", y = "Mandate"))      The `d` object is obtained with this CSV file. **Solution:** Thanks to **Wayne** , **Andy W** and others for their pointers! The data, code and graphs have been posted to GitHub. The result looks like this: ![result](http://i.stack.imgur.com/hvh2j.png)


title: Mean service time of a $M/E_2/1$ queueing system?

text: Consider a $M/E_2/1$ queueing system, where the customer arrival rate is $\lambda$ and the service time distribution has a gamma distribution with parameters $2$ and $\mu$, i.e. with p.d.f. $\mu^2te^{-\mu t}$ , $t ≥ 0$ (1) How can I determine the mean of the service time distribution? (2) What is the traffic intensity $\rho$ of the system in terms of the parameters $\lambda$ and $\mu$? * * * My reasoning thus far: wouldn't I simply take the mean of a Gamma(2,$\mu$) distribution and thus just say the mean service time is $\frac{2}{\mu}$ but I'm guessing there has to be something more to it than that? As for the triffic intensity $\rho$ I do not know what it would be for an erlang-2 queueing system? My thinking is (using $c=1$ for the number of servers): $\rho = \frac{"mean.service.time"}{c*"mean.customer.interarrival.time"} = \frac{"mean.service.time"}{c*\frac{1}{"arrival.rate"}} =\frac{\frac{2}{\mu}}{1*\frac{1}{\lambda}}=\frac{2\lambda}{\mu}$ Note: Follow up question: Mean length of time spent queueing in $M/E_2/1$ system? _I am looking at question 5 for one of the 2012 RSS exams:_ http://www.rss.org.uk/uploadedfiles/userfiles/files/GD3_2012_final%20(web%20version).pdf


Often in research articles you read the researchers have controlled for certain variables. This can be done by methods such as matching, blocking, etc. But I always thought controlling for variables was something done statistically by measuring several variables that could be of influence and performing some statistical analysis on those, which could be done in both true and quasi experiments. So, for instance you would have a survey or other test in which you´d measure the independent variable and some possibly confounding variables and do some analysis.   * Is is possible to control for variables in quasi experiments?   * What is the link between methods such as matching and statistically controlling for variables?


title: Online algorithm for computing variance of complex numbers

text: I was wondering if there is an online algorithm like Welford 1962 algorithm exists for computing the variance of complex numbers. Check these links for more information: Accurately computing running variance Algorithm for calculating variance UPDATE: If $X$ is a complex-valued random variable, with values in complex space, then its variance is $E[(X-\mu)(X-\mu)^t]$, where $X^t$ is the conjugate transpose of $X$. This variance is also a positive semi-definite square matrix. Check Variance.


title: Reference for implementing generalized likelihood ratio test to determine online whether time-series mean has shifted

text: What is a reference that describes the "generalized likelihood ratio" test to determine **online** (i.e., meaning that we add an observation, then check, then add an observation, then check) whether the mean of a time-series has shifted? I need descriptions of algorithms to be implemented into computer code or reference implementations in python, C++, Java or anything really. ### References   * Basseville, M., & Nikiforov, I. V. (1993). Detection of abrupt changes: theory and application (Vol. 104). Englewood Cliffs, NJ: Prentice Hall.


title: Confusion related to variance calculation

text: I was referring to this wikipedia article for calculating the variance of a variable X which is the linear combination of other variables. As given there $$ Var(\sum_{i=1}^{i=N}a_iX_i) = \sum_{i=1}^{i=N}a_iVar(X_i) + 2\sum_{i<j}a_ia_jCov(X_i,X_j) $$ Now suppose I have different parameters suppose something I have values from only half of the N variables and the others are missing at that time I will use something like $$ Var(\sum_{i=1}^{i=N/2}b_iX_i) = \sum_{i=1}^{i=N/2}b_iVar(X_i) + 2\sum_{i<j}b_ib_jCov(X_i,X_j) $$ Similarly at times I might have values from only 1/3 of the predictors. At that time I will have something like $$ Var(\sum_{i=1}^{i=N/3}c_iX_i) = \sum_{i=1}^{i=N/3}c_iVar(X_i) + 2\sum_{i<j}c_ic_jCov(X_i,X_j) $$ So based upon availability of data, I will have different parameters for the weighted average Now the constraints are $$ sum(a_i) = 1 sum(b_i) = 1 sum(c_i) = 1 $$ How can I impose these constraints and learn all the parameters at once. Rather than taking a chunk of data learning first suppose parameters $a_i$. Then take another chunk and learn parameters suppose $b_i$ and so on. What will my objective function look like?


Most measurements have measurement errors.   1. Should measurement errors and uncertainty be accounted for in every statistical test?   2. Are duplicates better than single measurements if the sample size is large enough?


title: Post hoc ANOVA following a significant ANOVA?

text: I am analyzing a data set using a one-way, within-subjects anova with 12 levels. Each of the levels represents a unique trial type in a reaction time experiment. However, the trials can be reasonably grouped together into 4 equal groups, let us call them A through D. My primary interest is in showing that A differs from the controls (B-D). Following a significant Anova, linear contrasts show that this is the case. However, it seems that I should also demonstrate that B-D are not different from each other. Parsimony suggests I could run another Anova (rather than pairwise tests) to demonstrate this. However, I have never seen this done. Is this simply convention or is there a more rigorous reason that I am missing? Thanks!


title: Does a Neural Network actually need an activation function or is that just for Back Propagation?

text: I have a feed forward neural network (1 hidden layer with 10 neurons, 1 output layer with 1 neuron) with no activation function (only transfer by weight + bias) that can learn a really wonky sin wave (using a 2in1out window) with production usable accuracy trained via stochastic climbing in a couple seconds:               for (int d = 0; d < 10000; d++)         data.Add((float)(Math.Sin((float)d * (1 / (1 + ((float)d / 300)))) + 1) / 2);      I'm probably just drunk, but if you don't use an activation function do you lose that universal function approximator status? Or is it just for gradient descent / back propagation etc. to act as a differentiable function? Alternatively, have I probably just overlooked a bug and am actually secretly activating without knowing it? source in C# (draws on a form)


title: Simulating unfair die with java/R (programming done)

text: Task : Unfair die(6 sides) is being rolled n times. Probability of 1 is p1, probability of 2 is p2 and so on. Write a computer program, that for given n (n<100), the probability of set (p1,p2,p3,p4,p5,p6) and $x \in [n,600n]$ would find the probability of sum of dice values is less than x. Program cannot work more than 5 minutes. This is an extra question that will give me extra points, but so far nobody has done it. I quess beginnes computer scientist like me can learn from this code also, since i found 0 help with bias dice in the web and came up with roulette like solution. I kind of wanted to show the world my way also. I have 2 solutions - using geometrical and statistica probability. My question is: 1)Is it correct when i do it like this or am i going wrong somewhere ? 2) Which one you think gives me better answer geometric or statistical probability ? My intuition says it is geometric, because it is more reliable. i think it is correct answer that my code is giving me - more than 0.99..... usually. I prefer Java more since it is much faster than R with loops, but i gave R code also for statistical , they are very similar i hope it is not a problam. Java code :               import java.util.ArrayList;               public class Statistical_prob_lisayl2_geometrical {              public static double mean(ArrayList<Double> a) {             double sum=0;             int len = a.size();             for (int i = 0; i < len; i++) {                 sum = sum + a.get(i);             }             return (sum/len);         }              public static double geom_prob(double p1,double p2,double p3,double p4,double p5,double p6){             ArrayList<Double> prob_values = new ArrayList<Double>();                  int repeatcount = 1000000;             int[] options = {1,2,3,4,5,6};             int n = 50;             double[] probabilities = {p1,p2,p3,p4,p5,p6};             for (int i = 0 ; i < repeatcount ; i++ ) { // a lot of repeats for better statistical probability                 int sum = 0; //for each repeat, the sum is being computed                 for (int j = 0; j < n ; j++ ) { // for each repeat there is n cast of dies and we compute them here                     double probability_value=0; // the value we start looking from with probability                     double instant_probability = Math.random(); // we generate random probability for dice value                         for (int k = 0; k < 6; k++ ) { // because we have 6 sides, we start looking at each probability like a roulette table                             probability_value = probability_value + probabilities[k]; // we sum the probabilities for checking in which section the probability belongs to                             if (probability_value>instant_probability) {                                 sum = sum + options[k]; // if probability belongs to certain area , lets say p3 to p4, then p3 is added to sum                                     break; // we break the loop , because it would give us false values otherwise                             }                         }                     }                 double length1 = (600*n)-n-(sum-n); //length of possible x values minus length of sum                 double length2 = 600*n-n;                 prob_values.add( (length1/length2) ); // geometric probability l1/l2                      }             return mean(prob_values); // we give the mean value of a arraylist, with 1000000 numbers in it         }         public static double stat_prob(double p1,double p2,double p3,double p4,double p5,double p6){             ArrayList<Double> prob_values = new ArrayList<Double>();                  int repeatcount = 1000000;             int[] options = {1,2,3,4,5,6};             int n = 50;             double[] probabilities = {p1,p2,p3,p4,p5,p6};             int count = 0;             for (int i = 0 ; i < repeatcount ; i++ ) {                 int sum = 0;                 for (int j = 0; j < n ; j++ ) {                     double probability_value=0;                     double instant_probability = Math.random();                         for (int k = 0; k < 6; k++ ) {                             probability_value = probability_value + probabilities[k];                             if (probability_value>instant_probability) {                                 sum = sum + options[k];                                     break;                             }                         }                     }                 int x = (int)Math.round(Math.random()*(600*n-n)+n);                 if( x>sum ) {                     count = count + 1;                   }             }             double probability = (double)count/(double)repeatcount;             return probability;         }              public static void main(String[] args) {             System.out.println(stat_prob(0.1,0.1,0.1,0.1,0.3,0.3));             System.out.println(geom_prob(0.1,0.1,0.1,0.1,0.3,0.3));          }          }      R code:               repeatcount = 100000     options = c(1,2,3,4,5,6)     n = 50     probabilities = c(1/10,1/10,1/10,1/10,3/10,3/10)          count = 0     for (i in 1:repeatcount) {     sum = 0     for (i in 1:n) {         probability_value=0         instant_probability = runif(1,0,1)         for (k in 1:6){             probability_value = probability_value + probabilities[k]             if (probability_value>instant_probability) {                 sum = sum + options[k]                  break             }             }         }         x = runif(1,n,600*n)         x         sum         if ( x> sum ) {             count = count + 1            }     }     count     probability = count/repeatcount     probability


title: Clustering large movie dataset using k-medoids?

text: I have to cluster a movie dataset of 10000 movies. A movie has attributes like Genres, Actors, Directors, Year. Earlier I thought that we can use a simple clustering algorithm like k-medoids and the distance can be pre-computed between two movies by subtracting genres & actors.               Initialise d(movie1, movie2) = 0          d(movie1, movie2) -= number of common genres     d(movie1, movie2) -= number of common actors     d(movie1, movie2) -= 1 (if they have a common director)     d(movie1, movie2) -= 1 (if they belong to same decade)      Is this approach correct? Is k-medoids fast enough to cluster this dataset (I doubt it isn't)? If it isn't fast enough any better clustering algorithm and strategy to cluster this dataset?


title: Literature on IV quantile regression

text: In the last months I have read intensively about quantile regression in preparation for my master thesis this summer. Specifically I have read most of Roger Koenker's 2005 book on the topic. Now I want to expand this existing knowledge to quantile regression techniques that allow for instrumental variables (IV). This seems to be an active field of research which is growing at a fast pace. Can maybe someone suggest to me:   * papers or other literature on IV quantile regression   * a short overview of these different statistical techniques   * the pros and cons of different techniques I'm looking mainly for literature to get me started and to have a good overview of what is out there. So the first point is the important one. The second and third would be nice to have! Also my interest lies mainly on the cross-sectional methods but also panel-methods are welcome. Thanks in advance.


title: Identifiability in linear regression

text: If we have a generative model: $X_2=X_1a_1+\varepsilon$, where $\varepsilon \sim \mathcal{N}(0,\sigma_2^2)$, do we have $X_1=X_2a_2+\varepsilon '$, where $\varepsilon \sim \mathcal{N}(0,\sigma_1^2)$ meaning unidentifiability? The context is the search of a structure (more complex than 1-1) between correlated covariates. I need to know if there is a chance to find the true generative model or if I have to consider that I will only find a group of equivalent models. Practically speaking, when I compare with the BIC if found different values... but it is not a proof... In the true model $\varepsilon$ is independent of $X_1$ but $X_2$ depends on $\varepsilon$. So $X_2|X_1$ is Gaussian, $X_2$ is Gaussian (or Gaussian mixture according to $X_1$) but the dependencies seems to disallow a simple permutation. It is not clear that $X_1|X_2$ is Gaussian (under the hypothesis of the true model (the first)) The context can be compared to SEM. The objective is to find the structure of the generative SEM knowing only $X$ the set of covariates. So the question is : "are there several (linear) structures leading to the same joint distribution for $X$ ?" This question stands for $X$ all Gaussian and then (if necessary) $X_1$ following Gaussian Mixture models. $P(X)=P(X_2|X_1)P(X_1)=P(X_1|X_2)P(X_2)$ the question is to define if $P(X_1|X_2)$ is a linear regression $i.e.$ $P(X_1|X_2)$ is Gaussian.


title: How do I fit a constrained regression in R?

text: I have a formal model from which I'm deriving some parameters that I would like to estimate. I haven't done this kind of thing before, and I'd like to have some help to solve this issue in R. I observe $z_i$ (money) and four characteristics of population. I want to regress these characteristics on $z_i$. My regression model should be $$z_i = \frac{\beta N}{\gamma} \ \textbf{1}_{\{\text{type 1}\}} + (\tilde{\alpha} + \tilde{\beta}) \frac{\beta N}{\gamma} \ \textbf{1}_{\{\text{type 2}\}} + (1 - \tilde{\beta}) \frac{\beta N}{\gamma} \ \textbf{1}_{\{\text{type 3}\}} + (\tilde{\alpha}) \frac{\beta N}{\gamma} \ \textbf{1}_{\{\text{type 4}\}} + \epsilon_i$$ where $\epsilon \sim N(\mu, \sigma^2)$ I want the weighting parameters $\tilde{\alpha}$ and $\tilde{\beta}$ to be obtained from a OLS regression, and I want all my coefficients to sum to 1. Parameters in the equation above are not identified, so I propose to solve this by re-parametrization:$$D_1 = \frac{\beta N}{\gamma} \quad D_2 = (\tilde{\alpha} + \tilde{\beta}) \frac{\beta N}{\gamma} \quad D_3 = (1 - \tilde{\beta}) \frac{\beta N}{\gamma} \quad D_4 = (\tilde{\alpha}) \frac{\beta N}{\gamma}$$ From this I can find $$\frac{D_2}{D_1} = \tilde{\alpha} + \tilde{\beta} \qquad \frac{D_3}{D_1} = 1 - \tilde{\beta} \qquad \frac{D_4}{D_1} = \tilde{\alpha}$$ I want to estimate the equation above in R to get $\tilde{\alpha}$ and $\tilde{\beta}$. How can I do that? Thanks


title: Correlated Random Effects Probit vs. GEE Population-Averaged Probit

text: My question relates to recent work on correlated random effects probit models (see these slides from Wooldridge) and comparing them to GEE population averaged probit models: Is one approach better as accounting for within- subject correlations than the other? CRE (more mixed effect?) models account for within-subject variation by including a within-subject variable ($\bar{x}_{i}$) capturing the mean of each time-variant explanatory variable ($x_{it}$). The statistical significance of this parameter is then a test of traditional RE assumption that the individual effects are independent of the explanatory variables. It also serves as a correction for violations of this assumption, but requires (I think) assuming that the individual effect is normally distributed around the mean of the individual's values. In GEE models, from what I understand, within subject correlations are accounted for via the correlation structure of the error term. What I don't quite understand is if/why the betas are not affected by within-subject correlations and the extent to which we need to worry about that. I understand that GEEs can't recover the individual effect, but I don't follow why much of the action is in the errors and there is little discussion of how within-subject variation might bias the population estimate of the betas.


title: How to interpret a log transformed (X+1)?

text: Assume that I have a log transformed model as follows: Model 1: Y = a + b*ln(x) -> Interpretation: A 1% increase in X is associated with an average b/100 units increase in Y. If I add 1 to X to avoid having 0 values and get: Model 2: Y = c + d*ln(x+1) Should I interpret the model as "a 1% increase in (X+1) is associated with an average d/100 units increase in Y?" Or there are some better ways to interpret the model? Thanks.


title: which test to use here to test the significance of overlap with respect to a control?

text: i'm not sure what statistical test to use to compare a scenario like this: I have a list of candidate regions A (560k regions) and I want to know if my candidate regions enrich for a particular set of attributes X (14k items). So I did an overlap with A and X and the overlap is 3300 of X overlapping A. I also generated a randomly distributed, size match control regions B (560k regions) and did the overlap between B and X and found 1800 of X overlapping B. So I want to test if the 3300 overlap is statistically significant. Which test to use here? I kinda feel I can use Fisher's exact test but not sure if it's valid. Thanks!


title: What does the | in |Z| mean in mathematical expressions for distribution statistics

text: I am doing some self study on statistics and noticed that in the notes that I was using the $|Z|$ expression as attached in the photos below. I am confused with the "$|$" that is being used. The only logical conclusion that I was able to draw was that $P(|Z|\le 2)$ represents $P(-2 < Z < 2)$. I would appreciate some advice on what the "$|$" symbol means. ![enter image description here](http://i.stack.imgur.com/zNT9s.jpg) ![enter image description here](http://i.stack.imgur.com/Nk6x9.jpg)


title: Unbiased estimator of weighted sum of two poisson variables

text: Suppose that $X_1$ and $X_2$ are two random variables sampled from a Poisson distribution with parameter $\mu$. Let $T_1=\bar{X}$ be the sample mean and let $T_2=(1/3)X_1 +(2/3)X_2$. Are T1 and T2 both unbiased estimators?


title: Why does policy iteration algorithm converge to optimal value? (reinforcement learning)

text: I was reading Andrew Ng's lecture notes on Reinforcement learning and was trying to understand why policy iteration converged to the optimal value function $V^*$ and optimum policy $\pi^*$. Recall policy iteration is: $ \text{Initialize $\pi$ randomly} \\ \text{Repeat},\{\\ \quad Step1: Let \ V := V^{\pi} \text{ \\for the current policy, solve bellman's eqn's and set that to the current V}\\ \quad Step2: Let \ \pi(s) := argmax_{a \in A} \sum_{s'}P_{sa}(s') V(s')\\ \} $ Why is it that a greedy-algorithm leads to a the optimal policy and the optimal value function? (I know greedy algorithms don't always guarantee that, or might get stuck in local optima's, so I just wanted to see a proof for its optimality of the algorithm). Also, it seems to me that policy iteration is something analogous to clustering or gradient descent. To clustering, because with the current setting of the parameters, we optimize. Similar to gradient descent because it just chooses some value that seems to increase some function. These two methods don't always converge to optimal maximums and I was trying to understand how this algorithm was different from the previous ones I mentioned. * * * These are my thoughts so far: Say that we start with some policy $\pi_1$, then after the first step, for that fixed policy we have that: $ V^{\pi_1}(s) = R(s) + \gamma \sum_{s'}P_{s\pi_1(s)}(s')V^{\pi_1}(s')$ $V^{(1)} := V^{\pi_1}(s)$ Where V^{(1)} is the value function for the first iteration. Then after the second step we choose some new policy $\pi_2$ to increase the value of $V^{\pi_1}(s)$. Now, with the new policy $\pi_2$, if we do do the second step of the algorithm the following inequality holds true: $R(s) + \gamma \sum_{s'}P_{s\pi_1(s)}(s')V^{\pi_1}(s') \leq R(s) + \gamma \sum_{s'}P_{s\pi_2(s)}(s')V^{\pi_1}(s')$ Because we choose $\pi_2$ in the second step to increase the value function in the previous step (i.e. to improve $V^{(1)}$. So far, its clear that choosing $\pi_2$ can only increase V^{(1)}, because thats how we choose $\pi_2$. However, my confusion comes in the repeat step because once we repeat and go back to step 1, we actually change things completely because we re-calculate $V^{2}$ for the new policy $\pi_2$. Which gives: $V^{\pi_2}(s) = R(s) + \gamma \sum_{s'}P_{s\pi_2(s)}(s')V^{\pi_2}(s')$ but its is NOT: $V^{\pi_1}(s) = R(s) + \gamma \sum_{s'}P_{s\pi_2(s)}(s')V^{\pi_1}(s')$ Which seems to be a problem because $\pi_2$ was chosen to improve $V^{(1)}$, and not this new $V^{\pi_2}$. Basically the problem is that $pi_2$ guarantees to improve $R(s) + \gamma \sum_{s'}P_{s\pi_1(s)}(s')V^{\pi_1}(s')$ by doing $\pi_2$ instead of $pi_1$ when the value function is $V^{\pi_1}$. But in the repeat step we change $V^{\pi_1}$ to $V^{\pi_2}$, but I don't see how that guarantees that the value function improves monotonically at each repeat because $\pi_2$ was calculated to improve the value function when the value functions stay at $V^{\pi_1}$, but step 1 changes $V^{\pi_1}$ to $V^{\pi_2}$ (which is bad because I $\pi_2$ only improved the previous value function we had).


title: Can a machine learning algorithm be evaluated based on a random sample?

text: I am trying to evaluate how well (or bad) a semi-supervised algorithm is performing on a given dataset. The algorithms assigns one of 10 labels to each data point. The dataset is huge, and it's not possible to obtain gold standard labels to each data point. Up until now, I have only come across papers where the evaluation is done against the entire dataset (e.g. using SVM to assign 'positive' or 'negative' labels in a sentiment analysis task, where each data point is a document). Since I can't obtain gold-standard labels for the whole dataset, I was wondering whether it is scientifically valid to evaluate the performance on a randomly selected subset of the data. If yes, then how do I ensure that the subset is indeed a true representation of the entire dataset? I want to provide a valid argument that good performance on such a subset _implies_ good performance on the entire dataset. I don't have a statistics background, so any reference for this will be extremely helpful.


title: Confidence interval for probability equal to 0

text: I meet this real life problem. The simplified version is this: Assuming I have 10000 apples, either good or bad. What is the sample size if I want to be 95% sure that all the apples are good?


title: cumulative distribution function , cdf problem

text: I cannot understand how step 2 transformed to step 3, anybody help me please ??? ![enter image description here](http://i.stack.imgur.com/bsghu.png)


title: Looking for simple examples of how to calculate type III and type IV SS

text: I have data collected on five species of fish at half a dozen locations in a lake over four years. The categories are not (at all) fully crossed, and I have a lot of empty cells due to logistical challenges with the sampling. I'm interested in comparing the effect of location and year on the variance. Despite the empty cells, it seems to me (after extensive reading) that it should be possible to do that comparison by using a corrected SS (Type III or Type IV), and calculating intraclass correlations. But I don't really understand how a Type III or a Type IV SS are calculated, and as a result, I can't assess if those are the right tools to answer my questions. So what I'm looking for is **simple examples of how one calculates Type III and Type IV SS, preferably by hand** (because I learn best by doing arithmetic and drawing plots), **or in R** (but I can muddle through SAS if necessary). I've looked a lot but haven't found anything really helpful. Any suggestions (tutorials, examples, book chapters) would be welcome. (I'd also welcome any other suggestions for strategies to tackle this question given the (albeit cursory) description of the data above!)


title: High value of denorminated degree of freedom mixed model in R?

text: There are three breeds of total 80 animals, timed with 200 repeated observations/animal; response variable is milk yield (MY), predictor is day from birth (DFC) . I performed a mixed model on this data set, considering animal as random effect, also I included the interaction term between DFC and breed to see if the slopes are different for each breed, and got the following results: model <- lme(MY~ DFC + DFC*Breed, random=~1|Animal, data= my_data) Linear mixed-effects model fit by REML Data: na.omit(centred_phuong) AIC BIC logLik 95866.09 95929.37 -47925.05 Random effects: Formula: ~1 | Cow_code (Intercept) Residual StdDev: 1.274623 2.591664 Fixed effects: Red_MY ~ Centred_DFC + Centred_DFC * Breed Value Std.Error DF t-value p-value (Intercept) 2.1921929 0.2475080 20049 8.85706 0.0000 DFC 0.0479475 0.0003760 20049 127.51737 0.0000 Breed2 0.4672407 0.3501083 78 1.33456 0.1859 Breed3 -0.3924527 0.3500558 78 -1.12111 0.2657 DFC:Breed2 0.0086141 0.0005292 20049 16.27611 0.0000 DFC:Breed3 -0.0086937 0.0005311 20049 -16.37037 0.0000 I am concerning why the DF for slopes (DFC, DFC: Breed2, and DFC:Breed3, and the intercept) ) is very high. As this may explain for always significant P values even the variation/diferrence is very small. Is there any way to adjust/ or to solve this problem of high DF? Many thanks!


If the Kullback-Leibler distance between f and g is zero then f = g almost surely. I am really having trouble solving this. It should be simple!


title: Need help to choose the right statistical test for difference among 4 types

text: Which statistical test can I use to test if DNA damage is different between 4 cell types with respect to each other (comparing 2 cell lines at one time). The data does not seem to be normally distributed, and I have 3 experiments for each condition.


title: Reverse coded items factor analysis: all load onto one factor

text: I have 21 items that should load onto four factors. Several of these items had been reverse-coded. However, in a factor analysis, all the reverse coded items load onto one factor (when they actually should load onto four separate ones). I had reverse coded those items prior to running the CFA, but a previous question did suggest I am correct to do this, so what might be the issue here?


title: Mutual information versus correlation

text: Why and when we should use Mutual Information over statistical correlation measurements such as "Pearson", "spearman", or "Kendall's tau" ?


title: Error to report with median and graphical representations?

text: I've used a wide array of tests for my thesis data, from parametric ANOVAs and t-tests to non-parametric Kruskal-Wallis tests and Mann-Whitneys, as well as rank-transformed 2-way ANOVAs, and GzLMs with binary, poisson and proportional data. Now I need to report everything as I write all of this up in my results. I've already asked here how to report asymmetrical confidence intervals for proportion data. I know that standard deviation, standard error or confidence intervals are appropriate for means, which is what I'd report if all my tests were nicely parametric. However, for my non-parametric tests, should I be reporting medians and not means? If so, what error would I report with it? Associated with this is how best to present non-parametric test results graphically. Since I largely have continuous or interval data within categories, I'm generally using bar graphs, with the top of the bar being the mean and error bars showing 95% CI. For NP tests, can I still use bar graphs, but have the top of the bar represent the median? Thanks for your suggestions!


title: Longitudinal data: time series, repeated measures, or something else?

text: **In plain English:** I have a multiple regression or ANOVA model but the response variable for each individual is a curvilinear function of time.   * How can I tell which of the right-hand-side variables are responsible for significant differences in the shapes or vertical offsets of the curves?    * Is this a time-series problem, a repeated-measures problem, or something else entirely?    * What are the best-practices for analyzing such data (preferably in `R`, but I'm open to using other software)?  **In more precise terms:** Let's say I have a model $y_{ijk} = \beta_0 + \beta_1 x_i + \beta_2 x_j + \beta_3 x_i x_j + \epsilon_k$ but $y_{ijk}$ is actually a series of data-points collected from the same individual $k$ at many time-points $t$, which were recorded as a numeric variable. Plotting the data shows that for each individual $y_{ijkt}$ is a quadratic or cyclical function of time whose vertical offset, shape, or frequency (in the cyclical case) might significantly depend on the covariates. The covariates do not change over time-- i.e., an individual has a constant body weight or treatment group for the duration of the data collection period. So far I have tried the following `R` approaches:   1. **Manova**              Anova(lm(YT~A*B,mydata),idata=data.frame(TIME=factor(c(1:10))),idesign=~TIME);       ...where `YT` is a matrix whose columns are the time points, 10 of them in this example, but far more in the real data. **Problem:** this treats time as a factor, but the time-points don't exactly match for each individual. Furthermore, there are many of them relative to the sample size so the model gets saturated. It seems like the shape of the response variable over time is ignored.   2. **Mixed-model** (as in Pinheiro and Bates, _Mixed Effect Models in S and S-Plus_ )              lme(fixed=Y~ A*B*TIME + sin(2*pi*TIME) + cos(2*pi*TIME), data=mydata,          random=~(TIME + sin(2*pi*TIME) + cos(2*pi*TIME))|ID), method='ML')      ...where `ID` is a factor that groups data by individual. In this example the response is cyclical over time, but there could instead be quadratic terms or other functions of time. **Problem:** I'm not certain whether each time term is necessary (especially for quadratic terms) and which ones are affected by which covariates.     * Is `stepAIC()` a good method for selecting them?      * If it does remove a time-dependent term, will it also remove it from the `random` argument?      * What if I also use an autocorrelation function (such as `corEXP()`) that takes a formula in the `correlation` argument-- should I make that formula for `corEXP()` the same as the one in `random` or just `~1|ID`?      * The `nlme` package is rarely mentioned in the context of time series outside Pinheiro and Bates... is it not considered well suited to this problem?   3. **Fitting a quadratic or trigonometric model to each individual, and then using each coefficient as a response variable for multiple regression or ANOVA.** **Problem:** Multiple comparison correction necessary. Can't think of any other problems which makes me suspicious that I'm overlooking something.   4. As previously suggested on this site (What is the term for a time series regression having more than one predictor?), there are **ARIMAX and transfer function / dynamic regression models**. **Problem:** ARMA-based models assume discrete times, don't they? As for dynamic regression, I heard about it for the first time today, but before I delve into yet another new method that might not pan out after all, I thought it would be prudent to ask people who have done this before for advice.


title: Do we have a problem of "pity upvotes"?

text: I know, this may sound like it is off-topic, but hear me out. At Stack Overflow and here we get votes on posts, this is all stored in a tabular form. E.g.:               post id     voter id    vote type     datetime     -------     --------    ---------     --------     10          1           2             2000-1-1 10:00:01      11          3           3             2000-1-1 10:00:01      10          5           2             2000-1-1 10:00:01       ... and so on. Vote type 2 is an upvote, vote type 3 is a downvote. You can query an anonymized version of this data at http://data.stackexchange.com There is a perception that if a post reaches the score of -1 or lower it is more likely to be upvoted. This may be simply confirmation bias or it may be rooted in fact. How would we analyze this data to confirm or deny this hypothesis? How would we measure the effect of this bias?


title: Confidence interval for sample mean (possible error in text)

text: I am working through a problem on CI for a sample mean and I cannot get the answer listed in the text I have. I am wondering if I am missing something or whether the text may be incorrect. We are given a sample of 10 scores:               45,38,52,48,25,39,51,46,55,46      I get a mean of 44.5 and SD of 8.68 which is the same as the solutions listed. I also get a standard error of the mean of 2.75 which is the same as the solutions. However I cannot get the final answer. We are asked to calculate a 95% CI for the mean. When I calculate the margin of error (standard error of the mean * tcritical) I get the answer `2.75*2.262=6.2`. (I choose df= 9). Thus the CI I get is `38.3 - 50.7`. However, the text says it is `28.35-60.65`. Am I missing something? I would be really grateful for any comments. Thanks.


title: Density-based spatial clustering of applications with noise (DBSCAN) clustering in R

text: this question started as "Clustering spatial data in R" and now has moved to DBSCAN question. As the responses to the first question suggested I searched information about DBSCAN and read some docs about. New questions have arisen. DBSCAN requires some parameters, one of them is "distance". As my data are three dimensional, longitude, latitude and temperature, which "distance" should I use? which dimension is related to that distance? I suposse it should be temperature. How do I find such minimum distance with R? Another parameter is the minimum number of points neded to form a cluster. Is there any method to find that number? Unfortunately I haven't found. Searching thorugh Google I could not find an R example for using dbscan in a dataset similar to mine, do you know any website with such kind of examples? So I can read and try to adapt to my case. The last question is that my first R attempt with DBSCAN (without a proper answer to the prior questions) resulted in a memory problem. R says it can not allocate vector. I start with a 4 km spaced grid with 779191 points that ends in approximately 300000 rows x 3 columns (latitude, longitude and temperature) when removing not valid SST points. Any hint to address this memory problem. Does it depend on my computer or in DBSCAN itself? Thanks for the patience to read a long and probably boring message and for your help.


title: How to compute marginal effects in a Tobit model using R?

text: is anybody familiar with the `tobit()` command using the package AER? I'm searching for a command to compute the marginal effects for y (not for the latent variable y*). It seems to be $\phi(x\beta/\sigma)\beta$, where $\phi$ is the std.normal cumulative distribution function. But how can I compute those effects with R?