title: Categorical variables in Predictive models - Data structure question

text: I'd like to know how to best structure categorical variables. Based on some reading I did, most researchers suggest making the categorical variable into a factor. Let's say my categorical variables are something like               Variable ( no ordering here, 9999 is just missing)     1     2     3     -1     9999          Anothervariable  ( ordering here 5 is worse compared to 4 and so on )     1     2     3     4     5      Question   How do I structure the two variables above during data reduction?   Let's say I am using Variable clustering to reduce data. Is it better to make it a factor for variable clustering, or should I leave it alone as a numeric variable?


title: How to compute p-value for quadratic fit of (25) discrete data points?

text: **Preamble** : I've had a couple of statistics courses @ uni, but can't really say I got "hooked" (I read NNT's Black Swan _before_ starting the degree :), hence the depth of my statistical knowledge remained (so far) rather shallow and mechanical, lacking deeper insight and a "feeling" for the matter; maybe I'm just lacking real-world practice and application... I've got a set of data (only 25 data points); in `gnumeric` I chose a "quadratic fit" (see for instance here) which returns an $R^2$ value of $0.39$ and a quadratic fit expression of $-0.13 x^2 + 4.65 x - 0.57$ (rounded). Regarding _how to report results_ , from here I gather that three things generally ought to be reported (for simple regressions and the like):   * the **R-squared** value (goodness of fit _viz._ coefficient of determination)   * the equation (though I don't think this is always informative/interesting)   * the **p-value** (regarding the statistical _significance_ of the fit), where it'd be enough to write one of two things in a report: "p < .05" or "n.s." (not significant). What I'm missing from `gnumeric`'s output for _quadratic_ regressions is the **p-value** (this can easily be obtained for linear cases, though). Thus, I was wondering whether a $\chi^2$ goodness-of-fit test would help (after classifying the continuous fit values according to the discrete, original 25 data points). Barking up the wrong alley?


title: Homework: Bayesian Data Analysis: Priors on both binomial parameters

text: The following is a problem from _Bayesian Data Analysis 2nd ed_ , p. 97. Andrew Gelman has not included its solution in the guide on his website and it has been driving me crazy all day. Literally all day. For some data $y$, modeled as a binomial distribution with population $N$ and probability $\theta$ parameters, both of which are unknown. The problem sets up the question with this information: (1) Setting a prior on $N$ is difficult, as it only takes on positive natural numbers, so it is treated as $\Pr(N|\mu) = Poisson(\mu)$, where $\mu$ is unknown. (2) To define the prior on $(N, \theta)$, we have $\lambda=\mu\theta$. (The logic here is that it may be easier to formulate a prior considering the unconditional expectation of the observations, rather than the mean of the unobserved $N$.) (3) A potential noninformative prior is $p(\lambda, \theta) \varpropto 1/\lambda$. The part of the problem that I am hung up on is how to transform the variables and determine $p(N, \theta)$. The approach that I have attempted is to write $p(N,\theta|\lambda)p(\lambda, \theta)$, and eliminate the unwanted $\lambda$ via integration, that is $p(N,\theta)=\int_0^\infty C\mu^N/(exp(\mu)\lambda N!)d\lambda$, and substituting out $\mu$ with the relation $\mu=\lambda/\theta$. This approach reduces to $p(N,\theta)=C/(N+1)$, where $C$ is the constant of proportionality introduced from (3). This result concerns me, because it implies that the joint probability of some values of $\theta$ and $N$ only depends on $N$, and not on $\theta$. Furthermore, some vague bells are ringing from my quite decrepit multivariable calculus, attempting to remind me about Jacobians and coordinate transformations, but I am uncertain that this integration approach is even appropriate. I appreciate your help and insight.


title: How to analyse Self-Assessment Manikin data

text: I am trying to analyse data collected through the Self-Assessment Manikin (SAM), which uses three pictorial scales to assess the three dimensions of _pleasure_ , _arousal_ and _dominance_. Participants are explicitly instructed that they can tick either of the five pictures or _any of the four spaces in between_ , making it essentially a 9pt Likert item. ![The SAM Arousal scale](http://i.stack.imgur.com/NPKKL.png) Unsurprisingly, however, people mostly tick the boxes themselves. Out of 96 responses, only 5 ticked any of the spaces in between. I am wondering if this is a problem that can influence my results. So far, I've tried collapsing the data to a 5pt scale by aggregating `1` and `2`, `3` and `4`, `6` and `7`, and `8` and `9`. Is this a sensible step to take before applying any of the usual Likert statistics and visualizations or is it better to leave the data as is? Does collapsing the scale to 5 values give me a better justification to treat the data as interval?               function(x) ifelse(x > 5, floor(x/2)+1, ceiling(x/2))


title: What type of study is this?

text: Suppose someone is trying to estimate the diagnostic prevalence of ADHD in a cohort of individuals born in 1985. What type of study is this? This is not a prospective cohort study because we are not comparing two groups (i.e an exposed and unexposed group). All individuals in this study are exposed (born in 1985).


title: Compare increase rate over time 2 groups

text: If you have two groups lets say Y any Z that change over "time", how do you compare if the increase rate of x and y over time is same or not. How would i set this up in R


title: Possible values for sensitivity, specificity, precision and accuracy

text: I have some results where the tester claims the following values:   Sensitivity: 0.525, Specificity: 0.925, Precision: 0.516, Accuracy: 0.907 Where   Sensitivity=TP/(TP+FN),   Specificity=TN/(TN+FP),   Precision=TP/(TP+FP),   Accuracy=(TP+TN)/(TP+TN+FP+FN) I'm having trouble seeing how these reported values are possible; are they?


title: How to determine the accuracy of logistic regression in R?

text: I'm curious about how to understand the accuracy of my model which I computed with `glm( family = binomial(logit) )`. In some articles it is mentioned that we should perform chisq test with residual deviance with it's DoF. When I call summary() of my glm module. **"Residual deviance: 9109.9 on 99993 degrees of freedom"** Therefore when I perform pchisq test with these inputs: **1-pchisq(9110, 99993)** it returns 1. Hence it is much more greater than our significance level. So we are curious about why does it return 1, is it a perfect model ? In addition to these, here's the output of my Logistic Regression Model               Logistic Regression Model          lrm(formula = bool.revenue.all.time ~ level + building.count +          gold.spent + npc + friends + post.count, data = sn, x = TRUE,          y = TRUE)                                Model Likelihood     Discrimination    Rank Discrim.                                  Ratio Test            Indexes          Indexes            Obs         1e+05    LR chi2    1488.63    R2       0.147    C       0.774          0          99065    d.f.             6    g        1.141    Dxy     0.547          1            935    Pr(> chi2) <0.0001    gr       3.130    gamma   0.586         max |deriv| 8e-09                          gp       0.011    tau-a   0.010                                                    Brier    0.009                                              Coef    S.E.   Wald Z Pr(>|Z|)     Intercept      -6.7910 0.0938 -72.36 <0.0001      level           0.0756 0.0193   3.92 <0.0001      building.count  0.0698 0.0091   7.64 <0.0001      gold.spent      0.0020 0.0002  11.05 <0.0001      npc             0.0172 0.0057   3.03 0.0024       friends         0.0304 0.0045   6.82 <0.0001      post.count     -0.0132 0.0042  -3.17 0.0015       This is validation with bootstrap's output                 index.orig training   test optimism index.corrected    n     Dxy           0.5511   0.5500 0.5506  -0.0006          0.5518 1000     R2            0.1469   0.1469 0.1465   0.0005          0.1465 1000     Intercept     0.0000   0.0000 0.0002  -0.0002          0.0002 1000     Slope         1.0000   1.0000 0.9997   0.0003          0.9997 1000     Emax          0.0000   0.0000 0.0001   0.0001          0.0001 1000     D             0.0149   0.0149 0.0148   0.0000          0.0148 1000     U             0.0000   0.0000 0.0000   0.0000          0.0000 1000     Q             0.0149   0.0149 0.0148   0.0001          0.0148 1000     B             0.0086   0.0086 0.0086   0.0000          0.0086 1000     g             1.1410   1.1381 1.1365   0.0016          1.1394 1000     gp            0.0111   0.0111 0.0111   0.0000          0.0111 1000      And this is the output of my calibration curve:               n=100000   Mean absolute error=0.002   Mean squared error=5e-05     0.9 Quantile of absolute error=0.002      ![Calibration Curve](http://i.stack.imgur.com/2AUEX.png) Thanks.


title: Maximum Likelihood question

text: I need to obtain the functional (i.e.: statistic obtained from the empirical distribution function) associated to the most powerful test to contrast $H_0: f_o$ vs. $H_1: f_1$, using a sample of size $n$. In particular, my problem is that I don’t know how to construct this functional when I am not given a particular density function. Using the Neyman Pearson lemma, I know that the maximum power test can be performed by computing the likelihood ratio test as: $$ \Lambda(x_1,...,x_n)=\dfrac{{\sup}f_0(x_1,...,x_n)} {{\sup}f_1(x_1,...,x_n)} $$ However since I don’t have a specific density function, I don’t know if this expression can be further operated or simplified in any way. **EDIT** : Since I got no answers or comments, my guess is that there may be something wrong with my question. **EDIT 2** : Per comments I understand now that I should use non parametric statistics and instead of conventional maximum likelihood, the appropriate ratio to use should be based on empirical likelihood. In the introduction of this paper: Empirical Likelihood Based Hypothesis Testing, there is a derivation of an statistic based on Empirical likelihood ratio to test $H_0: F = F_0$. However, since my data do not have the prerequisite of $F$ being a continuous distribution I am not sure if such reasoning can be applied to my original question. Any hint will be greatly appreciated!


title: Can we make probabilistic statements with prediction intervals?

text: I've read through the many excellent discussions on the site regarding interpretation of confidence intervals and prediction intervals, but one concept is still a bit puzzling: Consider the OLS framework and we've obtained the fitted model $\hat y = X\hat\beta$. We're given a $x^*$ and asked to predict its response. We compute $x^{*T}\hat\beta$ and, as a bonus, we also provide a 95% prediction interval around our prediction, a la Obtaining a formula for prediction limits in a linear model. Let's call this prediction interval PI. Now, which of the following (or neither) is the correct interpretation of PI?   1. For $x^*$ in particular, $y(x^*)$ lies within PI with 95% probability.   2. If we're given a large number of $x$s, this procedure to compute PIs will cover the true responses 95% of the time. From @gung's wording in Linear regression prediction interval, it seems like the former is true (although I could very well be misinterpreting.) Interpretation 1 seems counterintuitive to me (in the sense that we're drawing Bayesian conclusions from frequentist analysis), but if it's correct, is it because we're _predicting_ the realization of a _random variable_ vs. _estimating_ a _parameter_? (Edit) Bonus question: Suppose we knew what the true $\beta$ is, i.e. the process generating the data, then would we be able talk about probabilities regarding any particular prediction, since we're just looking at $\epsilon$? My latest attempt at this: we can "conceptually decompose" (using the word very loosely) a prediction interval into two parts: (A) a confidence interval around the predicted mean response, and (B) a collection of intervals which are just quantile ranges of the error term. (B) we can make probabilistic statements on, conditional on knowing the true predicted mean, but as a whole, we can only treat prediction intervals as frequentist CIs around predicted values. Is this somewhat correct?


title: A random variable that is invariant under null & alternative hypotheses

text: Suppose a random variable $X$ has probability density function $f_0$ and $f_1$ under null and alternative hypothesis, respectively. Is it possible to find another random variable, $g(X)$, which has the same density under null & alternative hypotheses? Thanks in advance! P.S. Not considering trivial cases such as $g(X)$= uniform [0,1].


title: How to find yearly return level form hourly data?

text: I have hourly values from which I extract every data over a threshold (Peak- over-threshold method). I then filter those extreme values so that i only get one value in a given 48 hours period. This leads to a Generalized Pareto (GP) distribution. I estimated the GP parameters with the standard Matlab built-in function. Now, I wish to estimate yearly (1-yr, to 100-yr) return level from the estimated GP parameters. I understand that I need to use the inverse GP function. My problem (in terms of understanding) comes from the quantiles I need to find for yearly return level considering I have hourly values. My thouht was to use (inside Matlab but the question is not related to Matlab):               m=1:1:100     gpinv(1-1/m,k,sigma,threshold)      My aim was to estimate the return level for return period between 1 yr and 100 yr. But I do not get the results I expected (the return level seems to low for any given return period).


title: What are the tests to measure agreement between two methods of measurement when the unit of measurement is different between the two methods?

text: Bland-Altman plot measure the agreement between two different methods which measures the same variable. As far as I understand Bland-Altman can measure only if the two methods have the same unit of measurement. I want to measure the agreement between two methods which measure the same variable (a continuous variable) but in different units, how can I proceed?


title: Oddly large R squared values in meta regression (metafor)

text: I am using the metafor package in R. I have fit a random effects model with a continuous predictor as follows               SIZE=rma(yi=Ds,sei=SE,data=VPPOOLed,mods=~SIZE)      Which yields the output:               R^2 (amount of heterogeneity accounted for):            63.62%     Test of Moderators (coefficient(s) 2):      QM(df = 1) = 9.3255, p-val = 0.0023          Model Results:                           se    zval    pval   ci.lb   ci.ub         intrcpt  0.3266  0.1030  3.1721  0.0015  0.1248  0.5285  **     SIZE     0.0481  0.0157  3.0538  0.0023  0.0172  0.0790  **      Below I have plotted the regression.The effect sizes are plotted proportionally to the inverse of the standard error. I realize that this is a subjective statement, but the R2 (63% variance explained) value seems a lot larger than is reflected by the modest relationship shown in the plot (even taking weights into account). ![enter image description here](http://i.stack.imgur.com/3JNmM.jpg) To show you what I mean, If I then do the same regression with the lm function (specifying study weights in the same way):               lmod=lm(Ds~SIZE,weights=1/SE,data=VPPOOLed)      Then the R2 drops to 28% variance explained. This seems closer to the way things are (or at least, my impression of what kind of R2 should correspond to the plot). I realize, after having read this article (including the meta-regression section): (http://www.metafor-project.org/doku.php/tips:rma_vs_lm_and_lme), that differences in the way the lm and rma functions apply weights can influence the model coefficients. However, it is still unclear to me why the R2 values are so much larger in the case of meta-regression. Why does a model that looks to have a modest fit account for over half the heterogeneity in effects? Is the larger R2 value because the variance is partitioned differently in the meta analytic case? (sampling variability v other sources) Specifically, does the R2 reflect the percent of heterogeneity accounted for _within the portion that cant be attributed to sampling variability_?. Perhaps there is a difference between "variance" in a non-meta-analytic regression and "heterogeneity" in a meta-analytic regression that I am not appreciating. I'm afraid subjective statements like "It doesn't seem right" are all I have to go on here. Any help with interpreting R2 in the meta-regression case would be much appreciated.


title: Non-central scaled Student's t cumulative density function required (alternatively the pdf)

text: I need to cite the pdf(density) or cdf(distribution function) of a non-central scaled Student's t distribution. There is an article about the non-central Student's t distribution http://en.wikipedia.org/wiki/Noncentral_t-distribution. But I have found no article that states the pdf or cdf and which also considers a scaling parameter. Note that the distribution should be scaled independently of the shape parameter $\nu$. In particular a parametrization as in the `gamlss.dist` package, function `TF2`, would be nice - i.e. mean $\mu$, degrees of freedom $\nu$ and a scale $\sigma$ that corresponds to the standard deviation. Has anyone seen something alike in a paper or book? Many thanks!!! Jo


title: Intercept term in logistic regression

text: Suppose we have the following logistic regression model: $$\text{logit}(p) = \beta_0+\beta_{1}x_{1} + \beta_{2}x_{2}$$ Is $\beta_0$ the odds of the event when $x_1 = 0$ and $x_2=0$? In other words, it is the odds of the event when $x_1$ and $x_2$ are at the lowest levels (even if this is not 0)? For example, if $x_1$ and $x_2$ take only the values $2$ and $3$ then we cannot set them to 0.


title: relative risk/ contingency tables

text: I have this question :the number of gun-related deaths per 100,000 people in 1994 was 14.24 in the United States, 4.31 in Canada, 2.65 in Australia, 1.24 in Germany, and 0.41 in England and Wales. Use the relative risk to compare the United States with the other countries. Interpret. I am not sure if I need to compare usa with the other countries at the same time. do something likt that: pi1/(pi2+pi3++pi5) or just compare usa to each country at a time like pi1/pi2 ; Pi1/pi3; etc...


title: Are there any good publicly accessible wisdom of crowd data sets?

text: I am looking for publicly accessible data sets to explore some wisdom of the crowd ideas. An ideal data set would contain predictions from a number of individuals (who were not collaborating with each other, or at least only minimally) on events with known outcomes. Predictions in terms of probability of events occurring (e.g., 80% chance of depression in the US in 20xx), binary-valued predictions (e.g., there will be a depression in the US in 20xx), or real values of quantities (e.g., average number of goals per match in World Cup 2014 is 2.75) are all okay. These predictions need to be paired with their corresponding objective ground truth, so that I can evaluate how different methods of aggregating data can produce estimates that are the closest to the ground truth. Any domain are okay. This includes but not limited to politics, economy, business, military, sports. Missing data (some predictors did not make predictions for all questions) are okay. I'd like the data to be readily downloadable (e.g., I don't want to write code to scrap a web site).


title: Interpreting online sample size stat tool in order to (gain) confidence intervals

text: Newbie, so please bear with me.. I'm trying to get my head around how I can use the awesome tool I found and what it does not give me any info about. The case: Using this site: http://hedwig.mgh.harvard.edu/sample_size/js/js_parallel_quant.html I have made a pilot study with 7 measurements and stated that I want a power of 80%. It gives me this output: > The provided parameters were: significance level (adjusted for sidedness) = > 0.025, standard deviation = undefined, number of patients = 7, power = 0.8, > difference in means = undefined, location of mean in one group as a > percentile of the other group = undefined. > > The variable calculated was the minimum detectable difference. > > A total of 7 patients will enter this two-treatment parallel-design study. > The probability is 80 percent that the study will detect a treatment > difference at a two-sided 0.05 significance level, if the true difference > between treatments is 2.661 times the standard deviation. My questions:   * is 2.661 the z-score; **___ __**?   * I consider 2.661 as a min. resolution that I'm able to differentiate between two different means. Is that a wrong assumption; **___ ____**? Would plotting the mean with top and bottom error bars with length of 2.661 x STD be the 95% confidence interval; **___ ____**? * * * A concrete example were I wish to use this knowledge: I wish to know something about the packing of spheres in a cylinder by dividing the total volume of spheres with the cylinder volume. I have filled up a die with sphere particles and weighted the content afterwards (n=7). Counting the particles and comparing it to the mass, gives me the knowledge that a similar experiment must have had an average of 188,7 particles with std = 2.2 particles. 85% of the particles, or above, have a diameter in the range of 1000-1400 µm. Plotting the average diameter of the particles against the packing gives me this: ![Plot](http://i.stack.imgur.com/WTRjW.png) Would it be fair to say that another sample with average in region 1 and 4, and a qual or smaller STD, would be detected by my experiment_ ___?


title: The best measure of reliability for interval data between 0 and 1

text: I have 6 sets of interval data each of which between 0 and 1. Each set, calculated by a computer program, is related to the degree of similarity between some sounds (pairwise). What do you think in the best inter-rater reliability measure I can use to see how close the 6 judges are? If I want to explain the data in each set, it can be: 0.98, 0.01, 0.5, ... which shows 'sound1' and 'sound2' are very similar (0.98), 'sound1' and 'sound3' are much different (0.01) and so on. Thank you so much.


title: How do I clean up inconsistent survey data?

text: I have some survey response data that are inconsistent and I'm not sure what is the best way to handle them. This is the nature of the inconsistency:   1. Have you ever had a food poisoning episode that you think might have been caused by the food in a particular restaurant ? a) Yes b) No   2. Did any of these episodes prevent you from eating at that restaurant again? a) Does not apply, never have had a food poisoning episode b) Yes c) No   3. Following these episodes, did you come back to your normal condition after 24 hours? a) Does not apply, never have had a food poisoning episode b) Yes c) No I believe that the intent of the designers of the survey was that anyone answering "No" to 1) should answer "Does not apply" to 2) and 3), but that's not how it worked. A large group of respondents picked "No" for 1) and "No" for at least one of 2) or 3). What is the best way to deal with this inconsistency? Where are some good references (books, papers, websites, etc.) to look at discussions and solutions for issues of this kind, and data cleaning in general?


title: Modeling null transitions in the Hidden Markov Model for use with the Viterbi algorithm

text: I've implemented the classic HMM model from Rabiner's tutorial for gesture recognition and it has worked well. Now, I'm trying to implement the HMM Threshold Model which calls for an HMM with **null transitions**. Null transitions do not provide any output symbols and the threshold model has 2 states which are connected to all the others through null transitions. How do I model this so that I can still use the Viterbi algorithm to evaluate the likelihood of an observed sequence?


title: Can a categorical variable be divided by a numerical variable?

text: I am reading P.R. Rosenbaum, Model-based direct adjustment, Journal of the American Statistical Association, 1987 (82), pp. 387-394. However, there is one equation, which I can't understand. $$ d = \frac{1}{N}\{\sum_{s=1}^{S}\sum_{i=1}^{N_s}\frac{z_{si}r_{si}}{\hat{e}_s} - \sum_{s=1}^{S}\sum_{i=1}^{N_s}\frac{(1-z_{si})r_{si}}{1-\hat{e}_s}\}, $$ where $d$ represents the difference, $\hat{e}_s$ is the probability, $z_{si}$ is the indicator function and $r_{si}$ is a binary response such as yes/no. My question is that How can a binary categorical variable $r_{si}$ be divided by a continuous variable $\hat{e}_s$? For example, how can a categorical variable (sex) male = 1 & female = 0 be divided by 0.23? It does not make any sense to me at all. many thanks in advance.


title: Calculating Total Probabilities from a Poisson Random Variable

text: First time post and I'm not sure if I'm posting in the right place but here goes: I enjoy Sports and was reading a book over Christmas (Mathletics by Wayne Winston) and to cut a long story short have started ranking sports teams based on least Squares Regression (based on the number of points scored or conceded) to generate Power Rankings. I've also done a little bit of reading around the subject such as on "The Probability of Winning a Football Game (Stern, The American Statistician, Aug 1991). I've noted that this works out fairly well in predicting results with some leagues (based on margins of error) and has led me to look at occurrences to do with probabilities of scoring goals. Below is a sheet that I have pulled from one of my Excel tables for a given league. The means are based on my least squares data from previous results, used to predict the likely number of goals based on a teams power rankings (offensively and defensively). The probability of scoring a goal I have treated as a Poisson Random Variable and have used the Poisson function in Excel to generate the numbers below:               mean    3.07554771 3.389358782       goals   team 1          team 2          team 1 wins      0   0.04616434      0.03373030       1   0.14198062  0.11432408  0.00478905      2   0.21833409  0.19374267  0.03232532      3   0.22383230  0.21888780  0.07650522      4   0.17210173  0.18547233  0.09649483      5   0.10586142  0.12572645  0.07898926      6   0.05426364  0.07102201  0.04731158      7   0.02384149  0.03438844  0.02248027      8   0.00916570  0.01456934  0.00895759      9   0.00313217  0.00548675  0.00310669      10  0.00096331  0.00185966  0.00096076      CUMULATIVE  1.00    1.00         tie 0.159139717          team 1 wins 0.371920564          team 2 wins 0.468939719           (sorry for the layout - I haven't got the hang of the formatting on here yet when pasting from Excel) My query is could I use this to predict the probability of an over under (e.g. over/under 5.5 total goals for instance)? I've got this far and had a total brain fade (high school statistics was a long while ago and the topics we learnt were never that interesting or relevant at the time). For what it's worth I've been able to interpret the means fairly well against totals (64% for this league) and I am just looking at ways of refining this (e.g. putting a probability out there as opposed to in this case 6.45, is over 5.5, I'd go over). I find it is always fun to have an opinion, but it's better when you can back it up with a statistics.


title: categorical variables in factorial design

text: One dependent variable, One independent continuous variable and One independent categorical variable (two levels only). Factorial design is easy except I don't know how to add center points. I am planning to have -1 and +1 , representing low and high, but I can't add a zero 0 for the center point because the categorical variable has two levels only. How should I do it?


title: calculating the expected value and variance of a log AR(1) process

text: I have an AR(1) process that looks like this: $$ \ln(g_t) = (1 - \rho_g)(\ln(\mu_g) - c) + \rho_g\ln(g_{t-1}) + \epsilon^g_t $$ where $|\rho_g| < 1$, $\epsilon^g_t \sim N(0, \sigma^2_g)$, and $c = \cfrac{1}{2} \left( \cfrac{\sigma^2_g}{1 - \rho^2_g} \right)$ and I want to find $E(g_t)$ and $Var(g_t)$. * * * I thought about starting out with a simple substitution $y = \ln(g_t)$, to make it a standard AR(1) process that I _can_ find the variance of: $$ y_t = B + \rho_g y_{t-1} + \epsilon^g_t $$ where $B = (1 - \rho_g)(\ln(\mu_g) - c)$. I can find the expected value: \begin{align} E(y_t) &= E(B + \rho_g y_{t-1} + \epsilon^g_t) \\ &= B + \rho_g E(y_{t-1}) \\ \mu_y &= B + \rho_g \mu_y \\ \mu_y &= \cfrac{B}{1 - \rho_g} \\ &= \ln(\mu_g) - c \end{align} and variance of this (using stationarity): \begin{align} Var(y_t) &= Var(B + \rho_g y_{t-1} + \epsilon^g_t) \\ &= Var(\rho_g y_{t-1} + \epsilon^g_t) \\ \sigma^2_y &= \rho^2_g \sigma^2_y + \sigma^2_g \\ &= \cfrac{\sigma^2_g}{1 - \rho^2_g} \end{align} but I can't see how this transformation helps me find $\sigma^2_g$. Apart from that I don't really know where to begin.


title: Estimate event probability from yes/no data

text: I need to predict (estimate) probabilities of (rare) events when the training data only contains the yes/no indicator. I.e., my target (dependent) variable is binary (logical). What I need is not just to predict yes/no, but estimate the probabilities of yes/no for each observation. If I use logistic regression, then the model output is, indeed, an estimate of the probability. What if I am using a different model, e.g., `vw`? (because, e.g., it is faster and outperforms logistic regression as a binary classifier). **So, I have a model which produces a score for each observation and I want to convert the score to probability.** It is natural to use total variation distance to evaluate the probability prediction, which motivated my previous question. The accepted answer there suggests Liblinear with L1 loss, but that produces a binary classifier, not a probability estimator. So, how do I calibrate model scores so that they actually estimate the event probabilities? I now train a single-independent-variable logistic regression to map the scores to probabilities. Can I do better?


title: Trouble training Neural Network

text: I'm trying to use Encog to define an artificial neural network in order to process this dataset (6 inputs, 2 yes/no outputs), but I can't get any lower than ~65% error. The NN is feedforward with backprop and sigmoid for activation. My steps were:   1. Normalize the dataset: temperatures from 35-42 to [-1,1], no to -1, yes to 1. Pick 80/120 random entries to use for training (the rest was supposed to be used for cross-validation - pun intended - later).   2. Use a NN with 6 input neurons, 2 outputs, bias, and backpropagation.   3. Start playing with learning rate, momentum and hidden layers/neurons count. The "best" result was achieved with 6 hidden neurons, 0.1 learning rate and 0.7 momentum. Tried all the common heuristics (half the input neurons, somewhere between the input and output count, using a large number of hidden units to facilitate the finding of a minimum, even if it raises generalization error, etc). Also, tried combinations of high/low learning-rate/momentum, but with no success. Funny enough, I did find a paper where the author seems to have solved this problem rather well, but is splitting the dataset by the 2 outputs a valid choice? I mean, the author used Levenberg Marquardt, which can only be used for 1 output neuron, so how is it supposed to encode the result in this way? The only way I can think of is training the NN to recognize output disease 1, and then do the same for another NN for output disease 2. I hope this was not too confusing, but I'm just getting started in this field of study.


title: How to show that there is no relationship

text: I have some data which I am using to show that there is no relationship between two variables. (Or only a weak one.) In a previous writeup, I included the scatterplot showing no visible relationship, as well as the Pearson's and Spearman's correlation coefficients, which were both low. One of the reviewers commented that the statistics in the paper were weak and we 'only showed some correlations'. Is there any better way to show that there is little/no relationship between these variables? I have transformed the data using every reasonable transformation (log, square root, exponentiation), and even averaged the values for each datapoint, though it seems to me that that is falling prey to the ecological fallacy. Nothing. Not even a hint of a pattern anywhere. Obviously, I can't show everything I tried. I want to show that there is no (or at least little) relationship, because common wisdom is that there should be a strong one. The fact that there is not a strong relationship is a surprising result. I know that you can't prove the null hypothesis, but I would like to show as much as possible that any other options are unlikely. **How do I convincingly show that there is no relationship between the variables?** (Other than hypothesis-testing my correlations, which I am planning to do, but that doesn't show that there is no relationship.) Note that both variables have heavy-tailed distributions - does that make a difference to the answer? Below is a scatterplot of my data. ![](http://i.stack.imgur.com/QCQeJ.png)


title: What if your random sample is clearly not representative?

text: What if you take a random sample and you can see it is clearly not representative, as in a recent question. For example, what if the population distribution is supposed to be symmetric around 0 and the sample you draw randomly has unbalanced positive and negative observations, and the unbalance is statistically significant, where does that leave you? What reasonable statements can you make about the population based on a biased sample? What is a reasonable course of action in such a situation? Does it matter when in our research we notice this imbalance?


title: How to interpret p value of regression coefficient which is nearly 0?

text: When regression coefficient is nearly 0 (in fact in the real model it's exactly 0), what's the meaning of p value (<0.05) of the coefficient? For example, I did a multiple variable regression with simulated data in R with lm(). Generate simulation data with the equation $$ y=2x_1^2+3x_2^2+3x_1+5 $$ The terms $x_1x_2$ and $x_2$ coefficients are zero. Using the data to do regression.               xmesh=mesh(seq(-4,4,0.1),seq(-4,4,0.1))     x1=as.vector(xmesh$x)     x2=as.vector(xmesh$y)     y=2*x1^2+3*x2^2+3*x1+5     model=lm(y~x1+x2+I(x1^2)+I(x2^2)+I(x1*x2))      summary(model)      The result is :               Call:     lm(formula = y ~ x1 + x2 + I(x1^2) + I(x2^2) + I(x1 * x2))          Residuals:            Min         1Q     Median         3Q        Max      -8.871e-12 -4.500e-15 -7.000e-16  5.700e-15  4.194e-12           Coefficients:                   Estimate Std. Error    t value Pr(>|t|)         (Intercept)  5.000e+00  3.301e-15  1.515e+15  < 2e-16 ***     x1           3.000e+00  7.545e-16  3.976e+15  < 2e-16 ***     x2          -3.348e-15  7.545e-16 -4.438e+00 9.22e-06 ***     I(x1^2)      2.000e+00  3.609e-16  5.542e+15  < 2e-16 ***     I(x2^2)      3.000e+00  3.609e-16  8.314e+15  < 2e-16 ***     I(x1 * x2)  -9.377e-16  3.227e-16 -2.906e+00  0.00367 **      ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 1.429e-13 on 6555 degrees of freedom     Multiple R-squared:     1,  Adjusted R-squared:     1      F-statistic: 2.313e+31 on 5 and 6555 DF,  p-value: < 2.2e-16       We can see that the coefficients of term $x_2$ and $x_1x_2$ are nearly 0, and p-value<0.01. I think that lm() did the significance test of coefficient based on t-test with NULL hypothesis $\beta=0$. So p-value<0.05 should mean that the coefficient is significantly different to 0. However the coefficient should be 0 in my model. I am confused. How to interpret these two coefficients' significance? **Add a new test** $y=2x_1^2+3x_1+0.001x_2+5$               > y2=2*x1^2+3*(x1)+5+0.001*x2     > model3=lm(y2~x1+x2+I(x1^2)+I(x2^2)+I(x1*x2))      > summary(model3)          Call:     lm(formula = y2 ~ x1 + x2 + I(x1^2) + I(x2^2) + I(x1 * x2))          Residuals:            Min         1Q     Median         3Q        Max      -9.237e-12 -1.700e-15 -1.000e-16  2.200e-15  2.757e-12           Coefficients:                   Estimate Std. Error    t value Pr(>|t|)         (Intercept)  5.000e+00  2.840e-15  1.761e+15   <2e-16 ***     x1           3.000e+00  6.492e-16  4.621e+15   <2e-16 ***     x2           1.000e-03  6.492e-16  1.540e+12   <2e-16 ***     I(x1^2)      2.000e+00  3.105e-16  6.441e+15   <2e-16 ***     I(x2^2)     -2.722e-16  3.105e-16 -8.770e-01    0.381         I(x1 * x2)  -3.226e-16  2.776e-16 -1.162e+00    0.245         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 1.229e-13 on 6555 degrees of freedom     Multiple R-squared:     1,  Adjusted R-squared:     1      F-statistic: 1.257e+31 on 5 and 6555 DF,  p-value: < 2.2e-16       You can see that in the new test, the coefficients and std Error of term $x_1x_2$ and $x_2^2$ are essentially zero. Their p-value are large enough to accept the null hypothesis that $\beta=0$, it's a good result. How to interpret the p-value of essentially zero coefficients in the two tests?


title: Standard deviation of errors

text: I have a multiple regression model (information derived from gretl): $$Y=‐3,859921 \ln(P) + 1,707514\, (A) + 3,578656$$ $$\quad\quad\space(1,216387)\quad\quad\quad(1,259650)\quad\quad(0,323867)$$ the standard errors of the parameter estimates are underneath. How to calculate the standard deviation of errors? I really appreciate your answers!!!


title: How can I represent Conditional Probability Table based Bayesian Networks in Winbugs?

text: Winbugs seem to support either stochastic or deterministic relationship between variables. However, many Bayesian Networks represent relationships between variables using conditional probability tables. The "visit to Asia", "burglar alarm", "smoking & cancer" examples are classic introductory material. However, conditional probability tables used in this approach does not correspond to stochastic or deterministic relationships in Winbugs. Google searches bring examples of Gibbs sampling on Bayesian Networks, but these are mostly algorithms that I'd have to implement in either R or another programming language (for example see: http://www- users.cselabs.umn.edu/classes/Spring-2010/csci5512/notes/gibbs.pdf) Is there a way of using Winbugs for complex Bayesian Network inference? I need to express causal relationships between different variables (a continuous variable and a discrete one, a categorical one etc), so I need to perform MCMC based inference on hybrid Bayesian networks. I am not so sure if I can use Winbugs for this Finally, would it make sense to try to transform the conditional probability table based, or causal relationships into a stochastic relationships? If I turn CPT based relationships between variables into regression for example, would that be an acceptable way of performing inference on Bayesian Networks using Winbugs?


title: How to test if two populations are different with the x axis a time variable

text: I have two populations: Class `A` and Class `B`. Both populations have two variables: `percent_correct` (i.e. number right answers / total answers) and `week`. I have graph comparing the two populations with the `percent_correct` on the y-axis measured over `time` (x-axis measured in weeks). The two populations have similar trends. However, I would like to do a statistical test to see if they are statistically similar. What test do I use? I'm confused as I'm looking at the similarities over time and not at one fixed point. Should I use the $r$ correlation? Or what do you advise? Also the variable `percent_correct` (y-axis) is independent of the variable `week` (x-axis). Sample data looks like this:               Class percent_correct date;       A        10%          04/01/2013       B        50%          04/01/2013       A        5%           04/10/2013       B        3%           04/10/2013       A        11%          04/20/2013       B        14%          04/22/2013       A        0%           04/28/2013       B        2%           04/27/2013


title: Predicting high frequency finance time series with HMM

text: I have a the following time series                 Price      BrokerID 632 Behaviour  BrokerID 680 Behaviour ...BrokerID XYZ Behaviour            5.6          IP                       SP                          5.7          BP                       IP       5.8          SP                       BP       5.83         IP                       SP      where `IP` is idle position, `BP` is buying position, and `SP` is selling position. I want to use Broker behaviour as the known variable and price as the hidden variable and predict it using HMM. But my question is how to find the emission matrix between a character vector (broker behaviour) and price numeric vector?


title: Confidence interval for classification error: binomial assumption vs. bootsrap resampling

text: I am developing a classifier using a set of N patterns, where N~1000. I am using K-fold cross-validation (with K=5) and computing the probability of classification error p (typical value is p=0.03). I am also computing estimates of the 95% confidence interval; I do this by assuming the classifier output is binomial distributed (see Brown, Cai, & DasGupta, Anirban, “Interval Estimation for a Binomial Proportion”, _Statistical Science_ , 16, 2, 2001, pp. 101-133). Now to my question. I have also heard mention of using bootstrapping to measure confidence intervals. I don't have a feel for why you would need to use a resampling method; can we not assume the number of errors to be binomial distributed? Does the cross-validation screw this up? Is resampling only necessary if N and p are such that directly estimating the confidence interval becomes difficult/impossible? This seems like a question that would come up any time the performance of a classifier is being analyzed.


title: What PDF should be fit to a rank histogram?

text: A Rank Histogram (or Talagrand Diagram) is a neat way of measuring whether your numerical model is giving appropriate variance. It's used for weather and climate forcasting, where you only have one observational series, and many model series (an ensemble). It's described pretty clearly here. Basically, you take a handful of runs of your model, then for each timestep/gridpoint, whatever, you calculate how the observations rank relative to the model ensemble. Then you plot a histogram of those ranks. If the histogram is u-shaped, then variance is too low (obs are rank high or low too often), if the histogram looks kind of gaussian, then the variance is too high (obs rarely ranks high or low), and if the histogram is flat, then your variance is spot-on (obs have a similar variance as the ensemble). Examples from http://www.eumetcal.org/resources/ukmeteocal/temp/msgcal/www/english/msg/ver_prob_forec/uos4b/uos4b_ko1.htm: ![Example rank histograms](http://i.stack.imgur.com/4nSdP.jpg) So the question is, **which distribution** should I fit to this kind of data, and **why**? The latter part of the question is more important, because I _think_ that the correct distribution is the Beta-binomial distribution, but I'm unfamiliar with this area.


title: How to model a network analysis problem

text: I have a weighted graph in which the nodes represent users and weighted undirected edges represent the tie between a pairs of users. For a piece of content $c$, and a node $A$ in the graph, given that $A$ has viewed the content and a few of $A$'s neighbors have viewed it as well, the problem is to infer whether $A$ liked the piece of content or not. I am currently thinking of using a Naive Bayes approach. Let's say $P(A_c = 1 | U_{1_c} =1 , U_{2_c} = 1,U_{3_c} = 0, U_{4_c} = 1,...,U_{nc})$ where $A_c = 1$ indicates that $A$ likes content $c$ and $U_{1_c}$ to $U_{n_c}$ indicate whether all of $A$'s neighbors liked the content $c$. In the above expression we are assuming that $U_{1_c}$ ..... $U_{n_c}$ are independent and have no effect on each other but the fact is that at some previous moment in time user $U_3$ might have adopted the content because he saw user $U_2$ and $U_4$ (who let us assume are also neighbors of $U_3$) adopt the content. My question is whether this assumption will affect the accuracy of the model? And if so could you suggest a better model for this problem. I have barely any statistics background so please forgive me for any blatant fallacies in the above equations or model.


title: neural network for rainfall forecasting

text: I'm new to machine learning, and I have been trying to figure out how to apply neural network to rainfall forecasting. I have found resource related to my query, but I seem to still be a bit lost. I think a basic explanation without too much detail would help. Let's say I have 10 climate features values for each month over a few years, and I want to predict rainfall for new month. How I can apply neural network to this problem is what I am trying to find out. how can i give input in this because we have 10(features) * 12 months * 10 years.   please help me .


title: How to treat rare / new items in the validation of a recommender system?

text: I'm working on movie recommendation algorithm. The data set consists of about 40 million ratings (user, film, rating). I want to separate the ratings into two groups - training set and probe set. The training set will be used for training, probe set for checking that the algorithm doesn't overfit the data. My first plan was to try to follow this approach: http://www.netflixprize.com/community/viewtopic.php?id=332. They create the probe set by taking the last 3 ratings from every user (and every user is guaranteed to have 20+ ratings). But the problem is, what if there is some really new film whose all ratings are in the probe set? The training algorithm wouldn't have any info about this film.


title: How to distance and to MDS-plot objects according their complex shape

text: Suppose I have **four basal forms** of signal (blue, purple, red, green). I also have created **_transition_ forms** between each other. If you carefully look on the picture below, you can see that for example blue signal ( _A1_ ) slowly transforms into purple ( _A5_ ) - horizontally, or into red ( _E1_ ) - vertically. Taken in general, **the closer is the difference between any two signals the more similar they are**. I am looking for some **method/algorithm/technique** which is able to:   1. **Extract as much information as possible about complexity of signals**   2. **Map signals into 2D (or 3D) according to their similarity** ![enter image description here](http://i.stack.imgur.com/EdnIJ.png) _Link to source data:_ each signal is coded into JPG image (25 jpg files, each 100x180 pixels) https://www.dropbox.com/sh/ynidcsjdymrh85f/AAAHVtYSG0GUvX3CEdDWWG42a * * * **_I have spent some time trying to solve this issue, so here I add my approach, which does no yield desired result:_** At first, I've set working directory, **load all 25 jpg-files into a R list**...               setwd("add directory where you downloaded jpg-files here")          # required libraries     library(jpeg)     library(raster)     library(asbio)          sgnl.vctr<-c("A1.jpg","A2.jpg","A3.jpg","A4.jpg","A5.jpg",                  "B1.jpg","B2.jpg","B3.jpg","B4.jpg","B5.jpg",                  "C1.jpg","C2.jpg","C3.jpg","C4.jpg","C5.jpg",                  "D1.jpg","D2.jpg","D3.jpg","D4.jpg","D5.jpg",                  "E1.jpg","E2.jpg","E3.jpg","E4.jpg","E5.jpg")           sgnl.list <- list()      # sgnl.list contains all 25 signals     for (i in 1:length(sgnl.vctr)){     sgnl.list[[i]] <- readJPEG(sgnl.vctr[i])     }                                                             I had some problems with pixel values (range from 0 to 1), therefore I **recclassify them into binary**. (either 0 or 1).               # reclassification of values     for (i in 1:25) {     sgnl.list[[i]][1:100,1:180,1][sgnl.list[[i]][1:100,1:180,1] > 0.5]<- 2     sgnl.list[[i]][1:100,1:180,1][sgnl.list[[i]][1:100,1:180,1] <= 0.5]<- 1     sgnl.list[[i]][1:100,1:180,1][sgnl.list[[i]][1:100,1:180,1] == 2]<- 0     }      Then, I've **extracted** binary (0-1) vectors from each jpeg-file. _If someone know how to shorten procedure bellow, please edit the R-code._               # A row     A1<-as.vector(sgnl.list[[1]][1:100,1:180,1])     A2<-as.vector(sgnl.list[[2]][1:100,1:180,1])     A3<-as.vector(sgnl.list[[3]][1:100,1:180,1])     A4<-as.vector(sgnl.list[[4]][1:100,1:180,1])     A5<-as.vector(sgnl.list[[5]][1:100,1:180,1])     # B row     B1<-as.vector(sgnl.list[[6]][1:100,1:180,1])     B2<-as.vector(sgnl.list[[7]][1:100,1:180,1])     B3<-as.vector(sgnl.list[[8]][1:100,1:180,1])     B4<-as.vector(sgnl.list[[9]][1:100,1:180,1])     B5<-as.vector(sgnl.list[[10]][1:100,1:180,1])     # C row     C1<-as.vector(sgnl.list[[11]][1:100,1:180,1])     C2<-as.vector(sgnl.list[[12]][1:100,1:180,1])     C3<-as.vector(sgnl.list[[13]][1:100,1:180,1])     C4<-as.vector(sgnl.list[[14]][1:100,1:180,1])     C5<-as.vector(sgnl.list[[15]][1:100,1:180,1])     # D row     D1<-as.vector(sgnl.list[[16]][1:100,1:180,1])     D2<-as.vector(sgnl.list[[17]][1:100,1:180,1])     D3<-as.vector(sgnl.list[[18]][1:100,1:180,1])     D4<-as.vector(sgnl.list[[19]][1:100,1:180,1])     D5<-as.vector(sgnl.list[[20]][1:100,1:180,1])     # E row     E1<-as.vector(sgnl.list[[21]][1:100,1:180,1])     E2<-as.vector(sgnl.list[[22]][1:100,1:180,1])     E3<-as.vector(sgnl.list[[23]][1:100,1:180,1])     E4<-as.vector(sgnl.list[[24]][1:100,1:180,1])     E5<-as.vector(sgnl.list[[25]][1:100,1:180,1])      **Vectors were compared** by `Kappa` function to obtain **total agreement value**. See this link: http://stackoverflow.com/questions/24534192/how-to- compare-all-possible-combinations-of-objects-in-r-by- loop/24534794#comment37992299_24534794 _(Many thanks to @digEmAll)_               # looping of:     #            Kappa(sgnl.list[[1]][1:100,1:180,1],     #                  sgnl.list[[2]][1:100,1:180,1])$ttl_agreement     M <- rbind(A1,A2,A3,A4,A5,                B1,B2,B3,B4,B5,                C1,C2,C3,C4,C5,                D1,D2,D3,D4,D5,                E1,E2,E3,E4,E5)          res <- outer(1:nrow(M),                  1:nrow(M),                  FUN=function(i,j){                    # i and j are 2 vectors of same length containing                     # the combinations of the row indexes.                     # e.g. (i[1] = 1, j[1] = 1) (i[2] = 1, j[2] = 2)) etc...                    sapply(1:length(i),                           FUN=function(x) Kappa(M[i[x],],M[j[x],])$ttl_agreement )                  })          row.names(res) <- c("A1","A2","A3","A4","A5",                         "B1","B2","B3","B4","B5",                         "C1","C2","C3","C4","C5",                         "D1","D2","D3","D4","D5",                         "E1","E2","E3","E4","E5")     colnames(res) <- c("A1","A2","A3","A4","A5",                        "B1","B2","B3","B4","B5",                        "C1","C2","C3","C4","C5",                        "D1","D2","D3","D4","D5",                        "E1","E2","E3","E4","E5")      Finally, similarity matrix (`res object`) was used for **multidimensional scaling**...               # mds based on ttl_agreement matrix     d <- as.dist(res)     mds.coor <- cmdscale(d)     plot(mds.coor[,1], mds.coor[,2], type="n", xlab="", ylab="")     text(jitter(mds.coor[,1]), jitter(mds.coor[,2]),          rownames(mds.coor), cex=0.8)     abline(h=0,v=0,col="gray75")      However, as you can see (left plot), **four basal signals were not separated as I expected**. Does anybody know better solution which would lead to desired outcome (right plot)? ![enter image description here](http://i.stack.imgur.com/rPlJe.png)


title: Why did I get significant difference result before using post-hoc test but I got an insignificant result after using post-hoc test?

text: I have done an one-way anova test for to compare among 4 groups. I did the homogeneity test for variance and the result showed the variance was not equal but there was a significant difference between group (p=0.013). Thus, I did the Dunnett's T3 test for post-hoc test. However, the p values between all the groups were more than 0.05. How would this happened? Does it mean there are no significant difference between the 4 groups? Thank you for answering my question.


title: Strange results from Chow test - surprisingly significant p-value

text: I have a set of data, sampled once per year from 1990 to 2008. It looks more or less linear. However, doing a Chow test results in a significant p-value (0.003 and below), wherever I test for a structural break. What is going on? (There was a similar question once before, but it didn't really receive a satisfactory answer.) Here is a screenshot showing what I am doing. The p-value is very small wherever I test to see if there is a structural break. An answer would be much appreciated!


title: Which clustering technique to use for a temporal dataset?

text: I have seen a similar question but thought I'd ask my own to hopefully garner some usefull feedback. Basically, I have a large temporal dataset, consisting of domestic smart energy meter use collected at 30 minute intervals for three months. I have data for about 120,000 households and I'm looking to apply an appropriate clustering technique to create customer profiles. I thought about using K means but I understand that it relies on each reading being independent of the previous which in temporal (and energy use) data is not the case. Am I right in this? Any suggestions are greatly appreciated and if you need any extra information, please don't hesitate to ask!


title: Discrete Time Survival Analysis - Correct Way to Write Survival Function

text: I have seen several ways to write (and calculate and interpret) a survivor function in discrete time survival analysis and I wonder which is correct or if they both are, but the interpretation and/or setup of the problem is different and I am missing it. Here is an example.   * Customers open an account at a bank.   * Record the month of opening and month of close and the discrete random variable T is in the set {0,1..,23}. After 23, there is right censoring. T denotes the number of month boundaries crossed between open and close. For example:   * A customer that opens in January 2012 and closes in January 2012 is labeled as 0.   * A customer that opens in January 2012 and closes in Feb 2012 is labeled as 1 An interval is of the form [a,a+1) where a=0,..,A. The following hazards, h(t), are calculated. For example, a customer opens and closes in the same month with probability 0.019194. **Question:** What is the proper way to construct the survivor function?   1. I am seeing that some sources set S(t) as the probability that the event will occur AFTER the period t: S(t) = **Pr(T>t)**. This is the yellow column below. It is calculated as S(t) = $\prod_{t=0}^{t}(1-h(t))$   2. Others will set the first period S(0) = 1 and then continue. This seems to be saying S(t) = **Pr(T>=t** ). It is calculated as S(t) = $\prod_{t=0}^{t}(1-h(t-1))$ In the continuous case, I guess it doesn't matter between Pr(T>t) and Pr(T>=t), but in the discrete case it does. ![enter image description here](http://i.stack.imgur.com/8KDxw.jpg)


title: Why are gaussian processes called smoothers

text: I have this confusion, why are gaussian processes called smoothers. I mean I know they are also used for regression. But why are they called smoothers. Any guidance will be much appreciated


title: Understanding F-test to compare variances (R's var.test)

text: I gave run the following codes in R console:               #for built-in matrix ToothGrowth     tooth.1mg <- subset(ToothGrowth, dose==1)     var.test(len~supp,data=tooth.1mg,alternative="two.sided")      I have found the result:               F test to compare two variances     data:  len by supp           F = 2.4176, num df = 9, denom df = 9, p-value = 0.2046          alternative hypothesis: true ratio of variances is not equal to 1           95 percent confidence interval:     0.6004952 9.7332038           sample estimates:     ratio of variances      2.41759       But I do not understand what the result expresses. Is that we fail to reject `null hypotheses` since `p-value>sig.level`? If so, then I found in most cases by implying the command `var.test` that we fail to reject `null hypotheses` since `p-value > sig.level`. Or, if I focus on confidence interval, it doesn't include hypothesized value which concludes that we reject `null hypotheses` in favor of `alternative hypotheses` at 0.05 significance level.


Does a random variable of t-square converge almost surely to another random variable? Here is the random variable $Tn$: $$ Tn=\frac{Z_{0}}{\sqrt{\frac{1}{n}\left ( Z_{1}^{2}+Z_{2}^{2}+...+Z_{n}^{2}\right )}} $$ Will there exist an RV $X$ such that $Tn^2$ converge almost surely to $X$ as $n$ tends to infinity? Why? Thank you very much. $Tn^2$ converges to a distribution $$ n^{\lim_{}\propto }P(Tn^2\leq x)=P(\sqrt{x}> Z> -\sqrt{x}) $$ According to the definition of converge almost surely, there is an event $A$ with $P(A)=1$ s.t. for all $w$ in $A$, $Xn(w)$ converges to $X(w)$ as $n$ tends to infinity. Can I say that there is no $P(A)=1$ exists in $Tn^2$ unless it converges to infinity?


title: Detecting a (statistically) significant change in time series trend

text: Suppose I have two time series (Series1 and Series2) which are identical from timePeriod -200 to timePeriod 0. Say that in timePeriod 0 they are both equal to 100. Series1is equal to 200 at timePeriod 50. Series2 is equal to 250 at timePeriod 50. They both go up and down over time, but Series2 is always equal to Series1+timePeriod number. (Ex. at timePeriod 20, it equals Series1+20.) Obviously there is a significant difference between the _levels_ of these 2 series at timePeriod 50. However, is there a way I can test to see whether there is a (statistically) significant difference in the _trends_ of these two series over time? (Preferably, the method could be used with different numbers and time periods)


title: Comparing AUC vs accuracy

text: I understand this question has been asked many times however, i am unable to understand the answers well enough and apply to my situation. I have attached 2 screenshots of my model. There are 5 class labels Firstly this is a model with a lower accuracy but with a higher ROC area. ![enter image description here](http://i.stack.imgur.com/vZPNO.png) Then i have a model with a higher accuracy but a lower ROC area. ![enter image description here](http://i.stack.imgur.com/krWAq.png) How do i decide which model is a better one if the cost for misclassification for all classes is the same


title: Dispersion parameter in GLM output

text: I have ran a glm in R, and near the bottom of the `summary()` output, it states               (Dispersion parameter for gaussian family taken to be 28.35031)      I've done some rummaging on Google and learnt that the dispersion parameter is used to fit the standard errors. I'm hoping somebody could provide more details on what the dispersion parameter is and how it should be interpreted?


title: Calculating $R^2$ two different ways

text: I am calculating my $R^2$ for a model two different ways. I used `pred = predict(model, newdata = mydata, se.fit = T)` and then calculated `1 - (mean(pred$se.fit))/var(mydata)` and this gave me a value of `0.9903`. However, then I calculated `rmse.model = sqrt(mean((mydata - pred$fit)^2))` and then found `1 - rmse.model^2/var(mydata)` and this gave me a value of `0.0372`....so I'm not sure what is happening. Thanks for any help.


title: How many hypotheses? And secondary outcomes?

text: Just wanted your opinion on whether you think its appropriate to have a 2nd hypothesis when you have included in your study a secondary outcome? And also do you think its appropriate to have more than one secondary outcome? I believe it should be fine but wanted to double check with you guys.


title: Use of standard error of the mean when the components are not identically-distributed

text: Let {$X_1$, ..., $X_n$} be a random sample of size n in which the different elements are measurements drawn from $m$ different populations with different distributions. From the Central Limit Theorem, as long as $n \gg m$ and $n \to \infty$, the values of $\bar{X}$ obtained by repeating this experiment will follow a normal distribution irrespective of the properties of the different $m$ starting distributions. Does this mean that it is valid to use the standard error of the mean as representative of the uncertainty on the $\bar{X}$? The different $\bar{X}$ are not independent measurements of the same variable, but the fact they they follow a normal distribution would make me believe it might be legitimate to do so.


title: Conjugated priors (Pareto and Beta), name of the unconditional distribution?

text: $$F_X(x)=\begin{cases} \quad\dfrac{\alpha}{\alpha+\theta}\left(\dfrac x\omega \right)^\theta &\text{ if } x<\omega \\ \\ 1-\dfrac{\theta}{\alpha+\theta}\left(\dfrac\omega x\right)^{\alpha} &\text{ if } x>\omega \end{cases}\quad\text{ where } 0\le x< +\infty.$$ I derived it by having $F_{X|Y}(x|y)=\left(\dfrac{x}{y}\right)^\theta$, which is $Y$ with a multiplicative "shock/noise" give by a $\mathrm{Beta}(\theta,1)$ (max of $\theta$ uniforms) and $Y$ follows a Pareto distribution $F_Y(y)=1-\left(\dfrac{\omega}{y}\right)^\alpha$. I was not able to find a classification for it. Does it have a name/family?


title: How to verify linearity assumption in linear regression with categorical predictors?

text: I have used simple linear regression, and I'm now checking that the model meets the assumption of linearity. The model used a continuous response variable and categorical explanatory variables. How can I asses linearity when using categorical explanatory variables?


title: perfect collinearity among multiple continuous variables

text: When there is a perfect collinearity among more than two continuous variables, how do you deal with it and how are the regression results interpreted? I have three independent variables which represent the percentage of different races within different cities, Say $x_1$ is the percentage of Hispanics, $x_2$ is the percentage of Blacks and $x_3$ is the percentage of Whites. Logically, $x_1 + x_2 + x_3 = 1$. The dependent variable is collage attendance among the population in different cities. I have eliminated one of the independent variables ($x_1$) and estimated the following equation: $\hat{y}=\beta_0 + \beta_1 X_2 + \beta_2 X_3 $   1. Is this the best way to deal with this problem?   2. How shall I interpret the coefficients of $\beta_1$ and $\beta_2$. If $\beta_1$ is equal to 3 for example, shall it be interpreted as: with 1 unit of increase in percentage of Blacks in the population of the city, the collage attendance increases by 3 units? Or shall I consider the omitted variable (X1: Hispanics percentage) as the base group and say that a unit of increase in Blacks percentage increases the collage attendance three times more than a unit of increase in Hispanics?


title: Seeking a continuous, parametric, bimodal sampling distribution for proportions

text: I am seeking a parametric probability model whose _pdf_ has the following characteristics: (1) it is supported on a variate axis that is bounded between 0 and 1; (2) it is continuous; and (3) it is capable of having two (or more) modes that are not necessarily located at 0 or 1 (e.g., in contrast to the beta distribution whose _a_ and _b_ are both <1). I am thinking of some sort expansion on the beta or Kumaraswamy distribution analogous to the bi-Weibull expansion on the Weibull distribution. My intention is to use such a model as a sampling distribution for sample proportions that vary continuously. My simulated data tend toward the low side of the variate range, typically ranging between 0 and 0.01, with a center of mass that in the neighborhood of 0.001. See the illustration below for examples. Has anyone run across such an obscure distribution that meets this challenging combination of criteria? ![Results of a Monte Carlo simulation exploring the temporal distribution of 30 random samples.  Four timeslices are shown.  The x axis describes the sample proportion for each single simulation run, the histogram, rug plot, and kernel density estimate describe the distribution of sample proportions for 2000 runs](http://i.stack.imgur.com/nkucw.png)


title: Use of survival analysis for trigger mining

text: I am tasked with analyzing data to find "triggers" to an event. Specifically, this is transaction data from a bank (e.g., checking account daily balances, daily over draft fees, daily number of checks cleared etc) and the event of interest is the checking account being closed by the customer. It sounds like I am really needing to do a feature selection (important features might be something like "having three over drafts in 6 months"). I was thinking about cox regression (possibly with time varying covariates) and some variable selection (the paper here was an inspiration). In this way significant variables could be considered triggers. The business will be using these triggers in monitoring software (if a trigger event is detected), flag the customer as an attrition risk, and then do something to keep them from leaving. Couple questions:   1. Does this seem like the correct approach or are there others?   2. It seems like it will be important to construct variables to determine the best triggers. For example, is the best trigger: a. Number of overdrafts in 3 months?   b. Number of overdrafts in last 7 days?   c. 3 or more overdrafts in the past 2 months   d. (...) Any suggestions on how to determine these triggers?


title: Typo in Wikipedia for Edgeworth expansion?

text: In Wikipedia for Edgeworth expansion, > Let $F_n$ denote **the cumulative distribution functions** of the variables > $Y_n$. ... > > Again, after inverse Fourier transform, **the density function** $F_n$ > follows as $$ F_n(x) = \Phi(x) + \sum_{j=1}^\infty \frac{P_j(-D)}{n^{j/2}} > \Phi(x)\,. $$ Is the same notation $F_n$ used for two different purposes? In the second,   * whose density function is it? I.e. which random variable's density function?   * Why does the cdf $\Phi$ of a standard normal distribution appear in the density function of some distribution?


title: I'm sampling a population, looking for a small proportion of individuals

text: I want to find some particular individuals in a population, ones with a special genetic makeup for instance. I estimate the population to be 10000 and I'm hoping to find one of 45 individuals. I sample 2000 but still haven't found any. Assuming my sampling is random, how do I put a figure on the likelihood of that happening?


title: How can you use HMMs and ANNs for on-line handwriting recognition?

text: > I've asked this question on cs.stackexchange before. It has a 20-hours > remaining bounty there. On-line handwriting recognition is the task of converting a series of $(x(t),y(t))$ coordinates to symbols and words. In contrast to off-line handwriting recognition, where you only have a bitmap of the written symbols, you have the way it was written in on-line recognition. The paper Npen++: An On-Line Handwriting Recognition System describes a system that makes use of Hidden Markov Models (HMMs) and Artificial Neural Networks (ANNs) to recognize symbols. > [...] each character is modelled with three states representing the first, > middle and last part of the character. Hence, the MS-TDNN can be regarded as > a hybrid recognizer that combines features of neural networks and hidden > Markov models (HMMs). Source: Npen++: An On-Line Handwriting Recognition System, page 6. One of my problems is that I don't know what a TDNN (Time delay neural network) is. At the moment, I think I've understood what feed-forward neural networks are and I know that they are trained with backpropagation. But all examples I have seen so far have a fixed input vector (e.g. colors of a fixed- size image; positions on a fixed-size board). How exactly are the lists of coordinates mapped to the input neurons? Is the following "mind model" I have of those HMMs correct? ![enter image description here](http://i.stack.imgur.com/q4pSA.png) When I understand it correct, then you have one HMM for each symbol. How do you get the probabilities $p_1, p_2, p_3, p_{1,2}, p_{2,3}$? How do you get the probability which character to recognize given the trained HMMs and a list of coordinates? When you can get a probability for every symbol with HMMs only, how does the ANN help?


title: Do infrequent examples screw up classifiers? If so, when is it okay to remove the infrequent examples from the data?

text: It's hard to think of a more eloquent way of phrasing this question - I'm basically wondering if a classifier trained on data where examples of some of the classes are infrequent/rare would be a bad model? I'm mainly interested in decision trees (C4.5). I think the answer is no, but that you will get a high error, because you will usually classify members of the infrequent classes as instances of the more frequent classes. This has been my experience so far. I'm also wondering when it's okay to remove these examples and when it's considered bad practice (i.e. doing it just to lower the error). I'm guessing that it's okay to remove these if there's a good reason to do so, and you explain that reasoning when you report your results. I'm not really interested in building the best classifier, I'm more interested in understanding relationships between the variables and the structure of the data. But all my variables are categorical and it's non-linear data, so decision trees have so far been the best tool I've found to do this. (SVMs and ensemble methods are more accurate, but you can't really see the internal model structure, which you get with decision trees.) thanks.


title: Confidence interval for the intercept in logistic regression

text: Some major commercial statistical packages (e.g., SPSS) do not report a CI for the intercept term in logistic regression. [Based on answer below R does provide CI for intercept] Why might confidence intervals for the intercept term not be included by default? UPDATE: Based on feedback, Confidence Intervals for odds associated with intercept term are reported in some stat. packages. And, obviously, they can be easily computed manually, knowing the standard error. Therefore, reformulated question: What is the interpretation of intercept's CIs (for odds) in logistic regression?


title: Test whether two multinomial samples come from the same distribution

text: This should be a very elementary question, yet I cannot figure out where I am going wrong. The matrix below contains data on the colour distributions of balls in two urns. I am looking for a formal method that can tell me whether the contents of the two come from the same population distribution.               freqs = c(25,94,85,47,13,1685)     data = matrix(freqs, nrow=2)     dimnames(data) = list("treatment"=c("Urn1","Urn2"), "outcome"=c("Blue","Green","Red"))      Plotting the (frequency-based) MLE's per urn, I can qualitatively observe that the colour distributions of Urn1 and Urn2 look pretty dissimilar.               toplot<- as.matrix(rbind(data[1,],data[2,] ))     barplot(toplot, beside = TRUE, col = c("green", "gray"), las=2);       ![enter image description here](http://i.stack.imgur.com/aSWTG.png) I have seen the $\chi^2$ independence test used to check 'association' between two sample sets like mine. When I run the test (below) I get the p_value < 2.2e-16 (below), which accepts (?) the null hypothesis that the colour distribution of sample set Urn1 is independent of the colour distribution of Urn2. I had expected to see a test result that indicates the two sample sets come from independent / different population distributions. I think I am mixing concepts here. Am I trying to use $\chi^2$ test for something that it is not meant for ? If so, which method should I use for my simple comparison?               result <- chisq.test(data)          #   Pearson's Chi-squared test     #     #data:  data     #X-squared = 884.9506, df = 2, p-value < 2.2e-16


title: How to use Levene test function in R?

text: I'm a newbie to statistics and R and I have a trouble with using Levene function (I would like to check the equality of variance of two samples). The documentation says that I should run: > levene.test(y, group) But I have no idea what I should put as y and group? I have two different samples which of I would like to check the equality of variance. Should I put one of the sample's values as y and the second as group parameter? Any hints?


title: How should I design a database for multi-year data analysis with JMP?

text: I'm building a database of economic and political indicators provided by the World Bank and other similar sources. I would like to run uni- and multivariate analysis of these indicators over the past 20 years of data. I'll be using JMP to analyze the data. Here's an excerpt of the database structure:               ID          | Pop. | GDP | (20ish other columns)     ------------------------------------------------     Afghanistan | xx   | xx  | xx     Albania     | xx   | xx  | xx      Algeria     | xx   | xx  | xx     Andorra     | xx   | xx  | xx     etc.        | xx   | xx  | xx      What would be the best way to input this data so that JMP can deal with all 20 years of data? Should I have 20 separate tables, one per year? Or should I have a column for each year all in one table (`GDP_1990`, `GDP_1991`, `GDP_1992`, etc.)


title: How can I calculate the hazard ratio (HR) from one-sided $p$ values on stratified log-rank?

text: Recently, I have been preparing to perform a meta-analysis about the treatment of pancreatic cancer. When I extracted time to event (survival-type) data from the publications, I encountered some obstacles. Most articles included showed the two-sided $p$ value for the log-rank test and the number of observed events in two arms and these data could be used to calculate the hazard ratio. But in one article, the author reported a one-sided $p$ value on a stratified log-rank test. My question is how do I get the HR from this? The title of the article is "Gemcitabine Alone Versus Gemcitabine Plus Radiotherapy in Patients With Locally Advanced Pancreatic Cancer: An Eastern Cooperative Oncology Group Trial"


title: How to compare group means for two samples with t-test using R?

text: ### Context I need to compare two groups of students. Students of these groups did some work, which later was evaluated. Now I have the values of the accuracy of every student. Here are the results: > 0.2 0.065 0.123 0.075 0.181 0.054 0.185 0.106 0.142 and > 0.269 0.357 0.2 0.221 0.275 0.277 0.253 0.127 0.246 I would like to check if the first group was better in terms of accuracy than the second group. As the samples are small I thought about using t test to compare the samples. H0 hypothesis was that the accuracy of both groups are equal H1 hypothesis was that the accuracy of the 1st group was greater than the accuracy of the 2nd group. I used the following R function:               t.test(precisionAdHoc$adhoc, precisionH4U$H4U, "g", 0, FALSE, TRUE, 0.95)      where `precisionAdHoc$adhoc` holds the data from the first group and `precisionH4U$4HU` holds the data from the second group. I got the following result:               data:  precisionAdHoc$adhoc and precisionH4U$H4U      t = -4.3687, df = 16, p-value = 0.9998     alternative hypothesis: true difference in means is greater than 0      95 percent confidence interval:      -0.1701334        Inf      sample estimates:     mean of x mean of y      0.1256667 0.2472222       ### Questions   1. Did I use the `t.test` function correctly? If not, what did I get wrong? and how should I correct it?   2. If the function was correct, how should I interpret the results?


What are the major machine learning theories that maybe used by Twitter for suggesting followers?


title: How can I implement multiclass k-NN?

text: I want to implement k-NN to use in a multi-class dataset. I found "A k-Nearest Neighbor Based Algorithm for Multi-label Classification" but didn't get the algorithm. Do you know any clear explanation of it? What I have in mind it works as follows:   1. Calculate posterior probabilities for each class (simply by dividing number of samples who are labelled as class_i to the number of total samples).   2. Then, for each test sample find it's nearest neighbours (e.g. assume k = 5, so we find top 5 nearest neighbour among all samples in the training set).    3. Then take majority class among these 5 neighbours as the class of the test sample (so if 3 of the 5 nearest neighbour has the Classx, then we'll classify the test sample as Classx). Would that be a correct implementation?


title: Exploratory analysis of spatio-temporal forecast errors

text: **The data:** I have worked recently on analysing the stochastic properties of a spatio-temporal field of wind power production forecast errors. Formally, it can be said to be a process $$ \left (\epsilon^p_{t+h|t} \right )_{t=1\dots,T;\; h=1,\dots,H,\;p=p_1,\dots,p_n}$$ indexed twice in time (with $t$ and $h$) and once in space ($p$) with $H$ being the number of look ahead times (equals something around $24$, regularly sampled) , $T$ being the number of "forecast times" (i.e. times at which the forecast is issued, around 30000 in my case, regularly sampled), and $n$ being a number of spatial positions (not gridded, around 300 in my case). Since this is a weather related process, I also have plenty of weather forecast, analysis, meteorological measurments that can be used. **Question:** Can you describe me the exploratory analysis that you would perform on this type of data to understand the nature of the interdependence structure (that might not be linear) of the process in order to propose a fine modelling of it.


title: Extremely new to R: Support.CEs package - help with design matrix please!

text: I've only started using R one day ago, still so much to get my head around. I need to create a choice experiment, and I have been following the example of H.Aizaki I think I have created a successful choice experiment design using the input below.               des1 <- rotation.design(attribute.names = list( Pepper = c("Black", "Mountain Hot", "Mountain Mild"), Packaging = c("Refill", "Grinder", "Shaker"), Label = c("Australia", "Tasmania", "Import"), Price = c("0.50", "1.25", "3.00")),nalternatives = 2, nblocks = 1, row.renames = FALSE,randomize = TRUE, seed = 987)     des1      and have converted it into a questionnaire format too. I need to make a design matrix now. I have used:               desmat1 <- make.design.matrix(choice.experiment.design = des1, optout = TRUE, categorical.attributes = c("Pepper", "Packaging", "Label"), continuous.attributes = c("Price"), unlabeled = TRUE)     desmat1[1:100, ]      But it only lists the first 27 lines (9 questions per survey, 3 options per question, so only 1 survey has been listed) before it says 'NA' e.g.               BLOCK QES ALT ASC Mountain.Hot Mountain.Mild Grinder Shaker Tasmania          ..............          25        1   9   1   1            0             1       0      0        1          26        1   9   2   1            0             1       0      1        0          27        1   9   3   0            0             0       0      0        0          NA       NA  NA  NA  NA           NA            NA      NA     NA       NA          NA.1     NA  NA  NA  NA           NA            NA      NA     NA       NA          NA.2     NA  NA  NA  NA           NA            NA      NA     NA       NA      Ideally I need 100 surveys (so 2700 lines I think???). I would love some help! Thank you in advance!


title: Estimating $n$ and $p$ for Binomial distribution, repeated counting of partly hidden population

text: A brief motivation: $n$ critters live in an aquarium, where sadly they often hide in, under or behind things. When the aquarium is observed, each critter is only seen with probability $p$ (independently), so the number of critters spotted is $X\sim \operatorname{Bin}(n,p)$. The observer takes repeated counts $X_1, X_2, ... , X_k$. Assume the time interval between counts is long enough that the critters have had time to rearrange themselves, so that $X_i$ and $X_{i+1}$ may be regarded as independent, but not so long that critters are born or die, so that $n$ is fixed. A first issue - I'll call it case (i) - is how $n$ and $p$ might be estimated solely from the vector of observations $(x_1, ..., x_k)$. That's basically this question here. There are method-of-moments and maximum likelihood estimators available. Neither seems good! How could a Bayesian estimation proceed without additional information, is there any sensible prior for $n$? A more interesting, and practical, scenario provides a sort of "anchor" for $n$. Suppose in case (ii) that a while before the observations were made, it was known there were $n_0$ critters, though there may have been births or deaths since. In the absence of any observations, in some sense $n_0$ remains the "best guess" for $n$. Once observations $x_i$ are available, this guess can be updated to reflect new information. This sounds more amenable to a Bayesian approach. Is there a good choice of prior? On the other hand I can't see how knowledge of $n_0$ would affect classical estimates of $n$ - not least because I haven't specified the relationship between $n_0$ and $n$. (It's unrealistic for the observer to know much more than that $n$ is likely near $n_0$.) Now consider case (iii), where critters may die but can't give birth in captivity. The fact there were $n_0$ critters in the past now provides an upper bound on $n$. An obvious non-informative prior for $n$ has probability $\frac{1}{n_0}$ for $n = 1, 2, ..., n_0$ and zero otherwise; perhaps another prior is more useful, but the constraint $n \leq n_0$ means it should look different to case (ii). For classical estimation, the likelihood of $n$ above $n_0$ is zero, so clearly MLE calculations proceed differently. Maybe other estimators can give absurd results with $\hat{n} > n_0$ - could their performance be improved by capping them at $n_0$? I hope the three cases are similar enough to bring together as a single question. In all three, a rough-and-ready $\hat{n}$ would be the maximum observation $x_i$, an approach that would be reasonably successful, albeit downwardly biased (it can't be an overestimate), if the sample size $k$ were large, and particularly if the probability $p$ of a critter being seen was near 1. But case (iii) is distinctive in that, if $x_i = n_0$ for any $i$, we can conclude with certainty that $n$ is $n_0$. I'd be interested in any practical approach, classical or Bayesian, to estimating $n$ and $p$, how such an approach is adapted between the cases, and how well such approaches can be expected to perform. Intuitively they will all struggle when $p$ is low, for instance, but some quantification of that would be brilliant. I'm sure the scenario described here is a well-studied practical problem - features like the partially-observed population of unknown size remind me of the German tank problem - but I can't find a name or discussion of it.


title: Multiclass classification with SVM a question about the feature vectors

text: I was told I had to direct my machine learning questions to this site. So here it goes. I'm trying to do Multiclass classification with SVM. I have 7 classes. Now I was wondering if the following is possible. I'm thinking of creating 7 SVMs for 1 vs all approach. Am i allowed to create 1 kind of feature vector per class? So e.g. Class 1 vs rest ==> Use type feature vector 1 (designed for class 1) Class 2 vs rest ==> Use type feature vector 2 (designed for class 2) Class 3 vs rest ==> Use type feature vector 3 (designed for class 3) And then assign the class-label with the highest confidence (probability), to the datapoint. Is this cheating ? Or is this allowed ? Is this common practice ?


title: Can someone enlighten me on what is "Neglected Nonlinearity"?

text: I ask this question out of curiosity earlier today when i was trying to test for heteroscedasticity in `R`, i accidentally mistook `white.test` of `tseries` package for `white.test` of `bstat` package. i found out later that the former tests for **Neglected Non Linearity** while the latter tests for **Heteroscedasticity** now this is something new, i hadn't heard about the **_"neglected"_** part before, can someone please enlighten me about the **_"neglected"_**?


title: Feature selection based on cost function

text: Suppose that we are searching for best `features` using an optimization algorithm for a classification model (MLP,SNM,Regression,etc...). We should set a cost function for this propose. This is one of the well-known cost functions:               NN = Selected_f     Z = MSE(1+beta*NN)      We should set `beta` in this equations. What `beta` value (between 0 and 1) do you think is suitable? I know this depends on problem and data set that I'm using, but as you know likes optimization algorithm parameters, maybe we have some proposed value for this. ( we don't know anything about how many features is good for this problem - we have 21 features) Can we only use `MSE` for cost function? Thanks.


title: Unimodal or bimodal data (MATLAB)?

text: I am trying to figure out what I did wrong or what I could do to get accurate results. I have _n_ vectors of data, and I am trying to decide whether each dataset is unimodal or bimodal. I assumed that it could be a mixture of Gaussians, so in MATLAB I attempted:               fit1 = gmdistribution(data(:,x),1)  %x is the index for each vector     fit2 = gmdistribution(data(:,x),2)       I used the minimum AIC of the two to decide whether it is unimodal or bimodal. The results return all to be bimodal, but they all aren't. Also by looking at the histograms of the data, they all don't seem bimodal. My question is why am I getting inaccurate or false positive results? Is there a better method? Did I do something wrong? SOME EDITS: using the gmdistribution function and the mean and variance from fit1 and fit2, I created random samples using:               Y = random(obj,length(original_data) %obj contains unimodal dist.     Y2 = random(obj2, length(original_data) %obj2 contains bimodal     hist(Y, length(original_data)     hist(Y2, length(original_data)      I compared each to the original vectors of data and the bimodal seems to provide samples that most resembles the original data. Is this a proper way of testing the AIC or dip test? I also used the hartigans dip test in matlab and the p values were close to 0 (less than 0.05), so I assumed this also means bimodal. I checked the kurtosis also, but they didn't return all bimodal. Any inputs or direction?


title: Can a single instrument reduce OVB and measurement error in your model?

text: My endogenous variable is very likely to be measured with error, but also confounded by unobservables. I have an instrument that rises with the value of the endogenous variable, at the same time it is exogenous to the the error in the outcome equation.


title: How does support vector regression work intuitively?

text: All the examples of SVMs are related to classification. I don't understand how an SVM for regression (support vector regressor) could be used in regression. From my understanding, A SVM maximizes the margin between two classes to finds the optimal hyperplane. How would this possibly work in a regression problem? Thank you in advance!


title: Asymptotic distribution of Kernel density estimator

text: For my research I am looking for proof of the asymptotic distribution of the univariate Kernel density estimator as proposed by Rosenblatt 1956 and Parzen 1962. A proof is for example given here and is also included in Nonparametric Econometrics by Pagan and Ullah. But since I am not writing about an economics topic I would like to cite a more "original" source. Does somebody know the paper in which the proof occurred first? Or an more basic math book would also be helpful.


title: Finding optimal hyperplane

text: I have a set of vectors $\{V_i\}$ in $n$-dimensional space. There is a number corresponded to each vector $\alpha_i = f(V_i)$ ($\alpha_i$ can be negative). I want to find a hyperplane which would maximize the difference between sums of $\alpha_i$ on the different sides of the space, divided by plane. What is the best way to do this?


title: Specifying a glmm for panel data

text: I'm trying to predict counts based on variables sampled on a monthly basis as well as a few that are not related to time. In several places I've read that the MCMCglmm package in R would be appropriate for this purpose, but I'm having a little difficulty understanding exactly how the panel aspect should be incorporated in a glmm model. Is it as simple as entering the data in wide format with each column representing the variable sampled during that month? Or do I add the data in long format with the id as a fixed effect, the columns as interaction terms between month:variable, month as a random effect, and y as the count for the next month? Or is there a better way to go about this?


title: $\chi^2$ distance and multivariate gaussian distribution

text: I want to know how to write the $\chi^{2}$ distance between two multivariate Gaussian distributions $f$ and $g$ in terms of their parameters only. The parameters of $f$ is the vector $\mu_{1}$ and a covariance matrix $\Sigma_{1}$. The parameters of $g$ is the vector $\mu_{2}$ and a covariance matrix $\Sigma_{2}$.


title: Truncating a posterior predictive distribution in JAGS

text: I have run into an error associated with truncating a distribution in JAGS. In my minimum reproducible example, I have data for 9 observations and would like to find a posterior predictive distribution for the 10th observation. To do this, I include the 10th observation as an NA and estimate its posterior predictive distribution as the variable `pi10`.               jagsdata <- data.frame(Y = c(47, 126, 68, 43, 67, 80, 61, 9, 26, NA))          model.string <- "     model{      for (k in 1:10){         Y[k]  ~ dlnorm(Z[k], tau.sp[k])         tau.sp[k] ~ dgamma(0.01,0.01)         Z[k] <- beta.o + beta.sp[k]      }       for (g in 1:10) {         beta.sp[g] ~ dnorm(0, 0.0001)       }       beta.o    ~ dgamma (2, 0.04)       pi10   <- Y[10]     }     "     writeLines(model.string, con = 'jagstest.bug')          library(rjags)     j.model  <- jags.model(file = "jagstest.bug",                             data = jagsdata,                             n.adapt = 500,                             n.chains = 4)     mcmc.object <- coda.samples(model = j.model,                                 variable.names = c('pi10'),                                 n.iter = 5000)      This works, but I would like to truncated the distribution of Y, for example by using the `T(1,200)`. However replacing line 4 above with                   Y[k]  ~ dlnorm(Z[k], tau.sp[k])T(1,200)      gives the error:               Unobserved node inconsistent with unobserved parents at initialization      Although Y with a normal distribution does not give an error.                   Y[k]  ~ dnorm(Z[k], tau.sp[k])T(1,200)      I have read through the JAGS manual section 7 and some examples online, but it is not clear to me how to implement this or why I am getting this error. Suggestions appreciated.


title: Mean survival time for a log-normal survival function

text: I've found plenty of formulas showing how to find the mean survival time for an exponential or Weibull distribution, but I'm having considerably less luck for log-normal survival functions. Given the following survival function: $$S(t) = 1 - \phi \left[ {{{\ln (t) - \mu } \over \sigma }} \right]$$ How does one find the mean survival time. As I understand it, $\sigma$ is the estimated scale parameter, and that exp($\beta$) from a parametric survival model is $\mu$. While I _think_ I can manipulate it symbolically to get t all by itself after setting S(t) = 0.5, what's especially stumping me is how to handle $\phi$ in something like R when it actually comes down to inputting all the estimates and obtaining a mean time. Thus far, I've been generating the survival function (and associated curves), like so:               beta0 <- 2.00     beta1 <- 0.80     scale <- 1.10          exposure <- c(0, 1)     t <- seq(0, 180)     linmod <- beta0 + (beta1 * exposure)     names(linmod) <- c("unexposed", "exposed")          ## Generate s(t) from lognormal AFT model          s0.lnorm <- 1 - pnorm((log(t) - linmod["unexposed"]) / scale)     s1.lnorm <- 1 - pnorm((log(t) - linmod["exposed"]) / scale)          ## Plot survival     plot(t,s0.lnorm,type="l",lwd=2,ylim=c(0,1),xlab="Time",ylab="Proportion Surviving")     lines(t,s1.lnorm,col="blue",lwd=2)      Which yields the following: ![enter image description here](http://i.stack.imgur.com/kUMEL.png)


title: Does logistic regression use crosstabs?

text: I am performing a binomial logistic regression with my outcome being a count with an upper bound which is 45. My data is of the following form: ![enter image description here](http://i.stack.imgur.com/oE2qX.png) Now I am wondering whether the actual analysis happens on a crosstab or not? I am confused since I have 45 categories and also continuous variables so that would make to many cells, increasing the chance of empty cells? I want to know, because I want to know how many degrees of freedom I have, and whether this depends on my sample size, or the cells of a crosstab.


title: Problem with the formulation of a gaussian copula likelihood function

text: I recently got to hear about copulas which to me sounded like a nice tool to model relationships between variables. I decided to try to implement the likelihood function for a bivariate Gaussian copula with normally distributed marginals in R for use in MLE estimation or Bayesian estimation. In this case I would expect that this likelihood function would have the same likelihood as a bivariate normal distribution. This doesn't seem to be the case (as demonstrated below) and my question is: **Where do I go wrong implementing the likelihood function for the Gaussian copula with normal marginals? Am I maybe misunderstanding something fundamental?** Here is the code for my likelihood function in R:               # x1 and x2 is the bivariate data, rho is the correlation     loglike_fun <- function(x1, x2, mu1, mu2, sigma1, sigma2, rho) {       # transforming the data according to a Gaussian copula       qx1 <- qnorm(pnorm(x1, mu1, sigma1, log.p=T), 0, 1, log.p=T)       qx2 <- qnorm(pnorm(x2, mu2, sigma2, log.p=T), 0, 1, log.p=T)       qx <- cbind(qx1, qx2)            loglike <- sum(dmvnorm(qx, c(0,0), matrix(c(1, rho, rho, 1), ncol=2), log=T))       loglike <- loglike + sum(dnorm(x1, mu1, sigma1, log=T))       loglike <- loglike + sum(dnorm(x2, mu2, sigma2, log=T))       return(loglike)      }      Now if I generate bivariate normal randomly distributed variables with `mu = 0`, `sigma = 1` and `rho = 0.75` as               x <- rmvnorm(100000, c(mu, mu), matrix(c(sigma, rho, rho, sigma),ncol=2))      I would expect the following to have the maximum likelihood:               loglike_fun(x[,1], x[,2], mu1 = 0, mu2 = 0, sigma1 = 1, sigma2 = 1, rho = 0.75)     ## -526812.6      But this is not the case as, for example, the following parameters consistently have a larger likelihood:               loglike_fun(x[,1], x[,2],mu1 = 0, mu2 = 0, sigma1 = 1.7, sigma2 = 1.7, rho = 0.9)     ## -484621.1      Any pointers to where I go wrong are much appreciated!


title: Algorithm to detect if a particular event (date) has an impact on a time series

text: I have many time series and a set of dates (key events eg public holidays or big sales days). There are too many for me to look as manually and determine if either an increase or decrease occured as a result of the event. Is there some algorithm that I could use? I am currently using R.


title: ARIMAX and xreg variables

text: Just wondering how to setup xreg variables for ARIMAX models? I am particularly interested in whether I should be grouping together events for dummy variables. For example should I create 1 dummy variable and store all public holidays or should I create a dummy variable for each public holiday (this could lead to many dummy variables). BTW - I am using R and the Forecast package.


title: Unique subsequences of the human genome

text: I am working on a project involving genomes. I am going to divide a genome into subsequences, and then do something with those subsequences. I want to know if the subsequences that I produce are going to be unique, or the probability that they will be unique. So assuming we are dealing with the base pairs of DNA, then we have the "characters" A,C,T, and G. Lets say I have subsequences of length 140. How would I figure out the probability that every subsequence is unique? Ideally I would know this versus length of the given genome (human, cat, ant...), so I can decide based on the probability of repetitive subsequences which genome to use in this project.


title: How to define fit robustness?

text: Suppose I want to fit a function to measured data minimizing the chi^2. "It is known" that if it has too many correlated parameters the fit will not only have very large expected statistical errors on the fitted parameters but many times it will fail (e.g. stuck in a minimum more than 10 standard errors away from the "true" minimum). Here's my question... How can the "robustness" of a function to be fitted be defined and/or studied? Example. Suppose my data are described by the function y=a*x+b*sin(x). x is in the range [-0.1,0.1] and a,b ~1. I take 6 measurements in this range. I can compute the expected correlation matrix and it is like 0.9xxx between a and b. However if the errors on the measured y are 1e-7 I will fit this function without any problems. If the errors on the y are about 1e-2 there is no possibilty to converge, and the actual chi^2 minimum will be tipically much farther then the expected error. How can I study a-priori whether a function to be fitted is "robust" or not?