title: How to model and test a decision support system (e.g. a terrorist warning system)?

text: **[Edit]** In light of the feedback I have received, I have added some more information below, to help clarify what I am seeking to 'model'. Additionally, I have reformualted the problem from a "dog barking warning system" to something else which hopefully, has less ambiguity. Instead, I will repose the problem as follows: Lets assume that my neighbour is a "mad scientist", who claims to have invented a "terrorist event" forecasting machine. The machine has three colored bulbs - one of which will be illuminated, depending on the severity of a terrorist event forecasted by the machine. The light bulbs have the following interpretation:   * No light bulbs illumintaed means there is no percived threat   * Green bulb illuminated means there is a level 1 threat imminent   * Orange bulb illuminated there is a level 2 threat imminent   * Red bulb illuminated there is a level 3 threat imminent To avoid getting too pedantic, lets assume for the sake of argument that the following terms are defined and agreed upon:   1. Level 1 terrorist event   2. Level 2 terrorist event   3. Level 3 terrorist event   4. "imminent" terrorist event What I am trying to find out, is if there is a way I can design an experiment that can help me say with a degree of confidence, whether the scientists claims are statistically significant or not. if the claims are found to be statistically significant, then I would like to be able to add a CI (confidence interval) to the claim. So I can say something like - if the orange bulb is illuminated, then a level 2 terrorist event will occur withing a x% CI. Having said that - IIRC, CI are only meaningful for N~ RV.? As an aside, I was thinking that something akin to Fischer's tea experiment (or running a bernoulli test would be usefull, but my stats 'fu is not what it used to be ... The purpose of such a model (asssuming that the scientist machine really does work), is to act as a decision support system - i.e. decisions can be taken on the output of the machine - IF it can be depended upon.


title: Asymptotic power

text: I have found the term "asymptotic power of a statistical test" only related to the Kolmogorov-Smirnov test (to be precise: asyptotic power = 1). What does this term acctually mean? In my opinion it should be someting like this: "if the alternative hypothesis is true, than for every significance level alpha there exists a sample size n that the selected test would reject the null hypothesis". Is "my" definition correct? According to "my defintion" the majority of classical tests (t-tset, ...) should have the asymptotic power 1, not only KS test. Am I right? ;)


title: Does alpha need to be adjusted for running Mann-Whitney U test independently for two dependent variables?

text: I am comparing two group means on two dependent variables (aggression and stress levels) using Mann-Whitney U test by SPSS. My question is that do I need to adjust the alpha level for this?


title: Testing that a finite quantity is equally partitioned

text: Suppose I have a set of observations on a finite (continuously divisible) quantity $Q$ (say $Q = 100$ to be concrete). For each observation $i$, $Q$ is partitioned into three parts $(A_i, B_i, C_i)$, so my observations look like $(10,50,40)$   $(33,34,33)$   etc Now I want to test the hypothesis that $A_i$ and $B_i$ are equal in some central tendency in the sample $i=1..N$. First I thought two-sample t-test; but these are not normal (and let's suppose my sample is small). So then I thought two-sample Wilcoxon signed rank or rank sum on $A_i$ vs $B_i$; but $A_i$ is not independent of $B_i$. So my latest thought is one sample Wilcoxon signed rank on $p_i = (A_i/(A_i+B_i))$ vs $0.5$. But is there a standard test here?


title: Compare distributions of user ratings?

text: I have a bunch of user ratings, over time, for different products, and I'd like to figure out the best ways to analyse and compare them. Here's how the data looks:               Product, Week, 1 star, 2 star, 3 star, 4 star, 5 star     Product A, Week 1, 2, 0, 0, 10     Product A, Week 2, 3, 1, 1, 9     Product B, Week 4, 3, 3, 2, 1     ...      Except that I have about 25 products and 20 weeks of data. I've already plotted:   * Box and whisker plots for quick comparison of rating distributions   * Histograms for each product for quick comparison of rating distributions   * Line graph of mean ratings over time One stats question: **What calculation should I use to determine whether Product A has a statistically significantly lower mean rating than the rest of the group? (overall)** And one visualisation question: **Is there a standard graph format to use to compare the distribution of each product's reviews over time?** (as opposed to the overall distribution, for which I can use a box-and-whisker plot by product)


title: What makes Hoeffding's inequality an important statistical concept?

text: On his blog, Larry Wasserman has a post about what he planned to cover in his course last fall. He notes that he was abandoning some classical topics in favor of more modern issues. One topic that he mentions is Hoeffding’s inequality. What makes this result especially important for students and practitioners?


title: Determining correlation in certain subsets of a dataset in R

text: I have a large dataset that has many variables. I'm trying to determine which variables correlate strongly with one specific variable. When you look at the entire dataset as a whole, the correlation of different variables is pretty weak. I know, however, that within certain subsets of the data the correlation is strong. For example: When variable A is between X1 and Y1 and when variable B is between X2 and Y2, that resulting subset has unusually large instances of the variable I'm trying to optimize. How can I determine using R which subsets of data have unusually large instances of the optimization variable when there are hundreds of variables to test against?


title: hausman test, random effect or fixed effect?which must be accepted

text: What is the limit for Hausman test?                hausman onefe one                           ---- Coefficients ----                  |      (b)          (B)            (b-B)     sqrt(diag(V_b-V_B))                  |     onefe         one         Difference          S.E.     -------------+----------------------------------------------------------------     FDIsharein~P |   -.0061599     .0113644       -.0175243               .       TradeofGDP |   -.0026483     .0023264       -.0049748               .            liter |    .0001273    -.0025736        .0027009               .          laglife |   -4.424448     5.585048        -10.0095        1.810234        lagmortal |    3.693138     1.375929        2.317209        .4881797           lagpop |   -7.018966     -.212882       -6.806084        3.597954        lagpercap |   -.0096546    -.0085061       -.0011486               .     Corruption~x |    .0276846     1.140028       -1.112344               .     Unemployme~l |    .0138161     .1309029       -.1170867        .0175023     Healthexpe~f |    .3305049     .5088894       -.1783846               .     Urbanpopul~l |    .3873496     .0862762        .3010734        .1565548     Militaryex~P |   -.0498372    -.0845921         .034755               .          PHYSINT |    .0910831     .7238854       -.6328024               .     ------------------------------------------------------------------------------                                b = consistent under Ho and Ha; obtained from xtreg                 B = inconsistent under Ha, efficient under Ho; obtained from xtreg              Test:  Ho:  difference in coefficients not systematic                           chi2(13) = (b-B)'[(V_b-V_B)^(-1)](b-B)                               =       34.85                     Prob>chi2 =      0.0009                     (V_b-V_B is not positive definite)      That is my hausman test, should I accept random effect or fixed effect?


title: Estimating logistic regression coefficients in a case-control design when the outcome variable is not case/control status

text: Consider sampling data from a population of size $N$ in the following way: For $k=1, ..., N$   1. Observe individual $k$'s "disease" status   2. If they have the disease, include them in the sample with probability $p_{k1}$    3. If they do not have the disease, include them with probability $p_{k0}$.  Suppose you observed a binary outcome variable $Y_i$ and predictor vector ${\bf X}_i$, for $i=1, ..., n$ subjects sampled this way. The outcome variable is **not** the "disease" status. I want to estimate the parameters of the logistic regression model: $$ \log \left( \frac{ P(Y_i = 1 | {\bf X}_i) }{ P(Y_i = 0 | {\bf X}_i) } \right) = \alpha + {\bf X}_i {\boldsymbol \beta} $$ All I care about are the (log) odds ratios, ${\boldsymbol \beta}$. The intercept is irrelevant to me. My question is: **Can I get sensible estimates of ${\boldsymbol \beta}$ by ignoring the sampling probabilities $\{ p_{i1}, p_{i0} \}$, $i=1, ..., n$ and fitting the model as if it were an ordinary random sample?** * * * _I am pretty much certain the answer to this question is "yes"._ **_What I'm looking for is a reference that validates this._** There are two main reasons I'm confident about the answer:   1. I've done many simulation studies and none of them contradict this, and    2. It is straightforward to show that, if the population is governed by the model above, then the model governing the sampled data is  $$ \log \left( \frac{ P(Y_i = 1 | {\bf X}_i) }{ P(Y_i = 0 | {\bf X}_i) } \right) = \log(p_{i1}) - \log(p_{i0}) + \alpha + {\bf X}_i {\boldsymbol \beta} $$ If the sampling probabilities did not depend on $i$, then this would represent a simple shift to the intercept and the point estimate of ${\boldsymbol \beta}$ would clearly be unaffected. But, if the offsets are different for each person this logic does not quite apply since you will certainly get a different point estimate, although I suspect something similar does. **Related:** The classic paper by Prentice and Pyke (1979) says that logistic regression coefficients from a case-control (with disease status as the outcome) have the same distribution as those collected from a prospective study. I suspect this same result would apply here but I must confess I don't fully understand every bit of the paper. Thanks in advance for any comments/references.


title: What is the difference between two discrete uniform distribution with the same range but different number of categories?

text: Two random number generators with uniform distributions having min, max as (0,8) The first generates all integers between 0 and 8 uniformly. But the second generates only [0,2,4,6,8] uniformly. What is the difference between the second generator and the first? Smaller range, but scaled up? A discrete vs. continuous distribution? Higher variance? How would one model the second distribution?


title: Measuring relative variability for variables with different scales

text: I want to compare the relative variability of several sleep-related variables in the same group of subjects. For example, is there more variability in time spent in REM sleep compared to the number of eye movement in REM sleep? These variables have different scales, which does not allow me to simply compare their variances. I have computed several relative measures of variation such as the coefficient of variation (SD/M), coefficient of dispersion (Mean Absolute Deviation from Median/Median) as well as the coefficient of quartile variation (Q3-Q1/Q3+Q1). To provide some context, several of my variables are highly skewed. My major question is: Is one of these coefficients more appropriate than the others? I know that with skewed data I will probably want to avoid the traditional coefficient of variation. But what about dividing the 'mean absolute deviation from the median'by the median? Would it be better (e.g., more robust) to use the 'median absolute deviation from the median' instead (and divide this value by the median)? I have not seen much discussion of the coefficient of quartile variation. There is also the question of using of the two estimators proposed by Rousseeuw and Croux (1993)- Sn and Qn - could they be used here? These last two measures are not necessarily symmetrical around any measure of location, so I am not sure whether it would be appropriate to use them in this context. If anyone has any other helpful suggestions, that would be helpful as well. Thanks, Kevin


title: Question about reading an output for using ANOVA to compare two linear models

text: I tried to compare the following two models using "anova.lm()" in R:               Model 1: score ~ gpa + class      Model 2: score ~ gpa          Res.Df     RSS Df   Sum of Sq      F    Pr(>F)      1     90     213                                   2     91     201  1          96  14.07  <2.2e-16       Since there is only one DF, how should I report the result of this F test? Is " **F(df=1) = 14.07, p-value < 2.2e-16**" correct?


title: Probability distribution of functions of random variables?

text: I have a very stupid doubt: consider the real valued random variables $X$ and $Z$ both defined on the probability space $(\Omega, \mathcal{F},\mathbb{P})$. Let $Y:= g(X,Z)$, where $g(\cdot)$ is a real-valued function. Since $Y$ is a function of random variables it is a random variable. Let $x:=X(\omega)$ i.e. a realisation of $X$. Is $\mathbb{P}(Y|X=x)=\mathbb{P}(g(X,Z)|X=x)$ equal to $\mathbb{P}(g(x,Z))$?


title: How can multiple regression be performed as a sequence of univariate regressions?

text: Let's say $x$ is correlated to both $y_1$, and $y_2$. Why are the residuals of the nested regression of $x$ against $y_1$ and $y_2$, not equal to the residuals of the simultaneous (multiple) regression of $x$ against $y_1$ and $y_2$? To clarify: I take the residuals of the regression of $x$ against $y_1$ to get residuals $r_1$. I regress $r_1$ against $y_2$ to get residuals $r_2$. Why are these residuals $r_2$ not equal to the residuals of the multivariate regression of $x$ against both $y_1$ and $y_2$? Written in `R` code, we would say that,               lm(lm(x ~ y1)$residuals ~ y2)$residuals      is not equal to:               lm(x ~ y1 + y2)$residuals      I would like to understand this as would like progressively to extract the influence of explanatory variables from a dependent variable, so that I can "magnify" progressively the dependent variable's correlation to each subsequent factor. I am doing this in the context of PCA regression so specifically:   * `it30` = the 30 year point on the Italian yield curve   * `itpc1` = the first principal component of the Italian yield curve, calculated from maturity points 1y, 2y, 3y, ..., 30y.   * `itpc2` = the second principal component of the Italian yield curve I expect `it30` independently to have a relationship to `itpc1` (yield curve level) and itpc2 (yield curve slope). Another fact is that, due to the PCA, `itpc1` and `itpc2` are orthogonal, but I do not think that is important for this question. Indeed: ![enter image description here](http://i.stack.imgur.com/gG6ij.png) And: ![enter image description here](http://i.stack.imgur.com/QsgPR.png) ...so the 30y yield curve has a relationship to both `itpc1` and `itpc2`. Now if I take the residuals of the first regression and regress these against the second variable `itpc2`, I would expect there to be a relationship, and there does seem to be: ![enter image description here](http://i.stack.imgur.com/t7YEK.png) So it appears that my residuals from the first regression are linked to the second variable, as I would expect, that is, after accounting for the first correlation, that is, extracting `itpc1` from the data, there is still information related to the correlation with PC2. Interesting so far. Now I want to extract both `itpc1` and `itpc2` from `it30`, but I am wondering which approach to take, because of the following that I do not understand.... My question is, why is _this_ regression not a perfect straight line? ![enter image description here](http://i.stack.imgur.com/8bPEI.png) That is, if I progressively extract correlated variables from the dependent variable in a nested way, why are the residuals not equal to the regression which extracts them all in one shot? My objective is to understand to what extent each principal component affects a series. Yes I know I can do this using the eigenvector matrix but I am interested in the above behaviour of regressions. Not being a pure math genius, any intuitive explanation accompanying formulas would be appreciated.


title: Bayes' theorem in 1-d EM algorithm

text: I'm watching a video on the EM algorithm, ![enter image description here](http://i.stack.imgur.com/Z5mcu.png) When we use Bayes' Theorem to calculate $b_i$, how do I find $P(b)$ and $P(a)$ initially? It says we can estimate the priors $P(b)$ and $P(a)$ but that's once we have at least one value $b_i$. Should I simply let $P(b)$ and $P(a)$ be 1 to begin with? Also, why does it say we _could_ estimate the priors? Don't we have to in order to calculate $b_i$?


title: Are discrete single value prior distributions always lost in MAP estimation?

text: I’d like to illustrate my problem with a little (heavily abbreviated) excercise. I think it will help a lot to stress my point. Meet Mary, Tom and Jane. They all are programmers. Mary is a decent programmer. In writing five programs she usually makes about 3 mistakes. Tom on the other hand is quite bad. He makes about one mistake in every program. Jane is by far the best. In ten programs she writes you will only find one mistake. On one day the their boss checks five programs written by the same person and finds two mistakes. He wonders who might have written that program. He knows that Jane is the only one who has a full-time job and hence it is two times as likely that she has written these five programs. He thinks that he might use a Poisson distribution and some basic Bayesian point parameter estimation to find an answer. The likelihood-function for all three would be: $$f(X = 2 | \mu_i) = \frac{\mu_i^2}{2}·e^{-\mu_i}$$ with $$\mu_1 = 3, \mu_2 = 5, \mu_3 = 0.5$$ for Mary, Tom and Jane. Our prior distribution is based on information regarding part-time and full- time jobs: $$ P(\mu) = \left\{ \begin{array}{l l} 1/2 & \quad \text{if $\mu = 0.5$}\\ 1/4 & \quad \text{if $\mu = 3$ or $\mu = 5$}\\ 0 & \quad \text{else} \end{array} \right.$$ Now, let’s write the posterior distribution down: $$f(\mu_j | X = 2) = \frac{f(X = 2 | \mu_j) · P(\mu_j)}{\sum_{i=1}^3f(X = 2 | \mu_i) · P(\mu_i)}$$ Further, let us use MAP estimation (i.e. - for readers not familiar with this terminology - maximum a posteriori estimation whereby we use the mode of the posterior distribution.) When we want to find the mode of the posterior distribution for our little example we will find that Mary is our (un)lucky star. In determining this we weighted the likelihood $L(\mu) = f(X=2|\mu_j)$ by $P(\mu_j)$. Hence, this information (however insignificant) influences our maximum estimation. But now we want to use MAP to find our parameter estimate for $\mu$. In general it holds that when using MAP it actually doesn’t matter whether we include the denominator of the posterior distribution to find the maximum or not. Hence, we can only maximize the numerator $$\underset{\mu}{\arg \max}\;f(X = 2 | \mu_j) · P(\mu_j)$$ We would then form the weighted loglikelihood $$\log(f(X = 2 | \mu_j) · P(\mu_j))$$ and take the derivative of the loglikelihood. And here is where it get’s interesting for me. Whenever we have a discrete distribution that is only given in terms of a few values (in our example $\mu_i$) then $P(\mu_i)$ will become zero when we take the derivative. In our example the derivative will be: $$\frac{\partial\log(L(\mu))}{\partial \mu}=\frac{2}{\mu} - 1$$ which is equivalent to MLE (i.e. maximum likelihood estimation) although we used an informative prior. In short: What we can see is that although in our estimation of the mode of the posterior distribution we weighted the likelihood with the probability of the person having produced that program based on whether they have a full-time or part-time job this information is lost in estimating the parameter with MAP. Questions: (1) Do we always lose the a priori information for MAP parameter estimation when using only single discrete values in our formula or am I missing something here? (2) Is this situation remedied when I model the three values by a distribution/formula that is dependent on $\mu$ but which also only allows the three values for $\mu_i$? Thanks for any help!


title: Intuition behind the Calinski-Harabasz Index

text: Given $CH(k) = [B(k) / W(k) ] \times [(n-k)/(k-1)]$, where $n$ = # data points $k$ = # clusters $W(k)$ = within cluster variation $B(k)$ = between cluster variation. It is my understanding that the CH index can show the optimal number of clusters when doing k-means or hierarchical clustering; you would choose the number of clusters $k$ that maximize $CH(k)$. As $k$ increases, $B(k)$ increases, and $W(k)$ decreases. However, can someone explain to me the intuition behind the second part of the formula, namely $[(n-k) / (k-1)]$? Isn't it too punitive for cases where $n$ is very large, since increasing $k$ by 1 will drastically decrease the whole term?


title: Help on a probability theory question

text: Let $X \sim U(0,1)$ and let $0 < a < 1$, then find i) $P(X \le x|X > a)$ ii) $P(X \le x|X < a)$ I have the answers, but I'm wondering if someone can explain the process to me.


title: Multiclass transductive SVM software package?

text: Anyone know of one? I was looking at SVM-light but it seems like you can only use it in the binary classification setting. There is also SVM-multiclass but that doesn't have tSVM support. Any help appreciated. Thanks!


title: Sensible Transformations of Economic Indices like CFNAI and ADSBCI in Time Series Analysis

text: I am trying to fit an unobserved components model for revenue and transactions for a firm where I also use some exogenous variables that capture economic conditions. The UCM decomposes a time series into trend, seasonality, an idiosyncratic component, and up to 3 cycle(s) using a Kalman filter state space approach with the parameter estimated by MLE. The basic idea is \begin{equation} y_t = \tau_t + \gamma_t + \psi_t+\beta x_{t-k}+\varepsilon_t \end{equation} where   * $y$ is the dependent variable (revenue or transactions)   * $\tau$ is the trend component   * $\gamma$ is the seasonal component   * $\psi$ is the cycle   * $\beta$ is a vector of regression coefficients on explanatory variables $x$ For my explanatory variables, I have been using lags of the moving averages of _monthly_ Chicago Fed National Economic Activity Indices or the _daily_ Philly Fed Aruoba-Diebold-Scotti Business Conditions Index. Both of these indices measure the common component(s) of many different economic time series using a dynamic factor framework similar to PCA. Their average values are zero. When the indices are positive, the economic activity growth is above trend or average. When the negative, growth is below trend. The predictive performance of this approach is reasonable. The interpretation of the $\beta$s, however, is a bit awkward since I have essentially a growth rate on the RHS impacting a level on the LHS. This imposes a restriction that growth has the same effect on the level regardless of the level. This makes little sense to me. I tried using $\Delta_s \ln y = \ln y_t-\ln y_{t-12}$ as the outcome to make it a growth-on-growth model, but the fit and the forecast quality really suffered Are there any transformations of the indices or the outcomes that I can try to avoid that restriction?


title: Time Series estimation on specific lags on MA part

text: I need to estimate specific lag ARMA model. Here is an example.               ts1 <- arima.sim(list(order = c(0, 0, 5), ma = c(0, 0.7, 0, 0, 0.1)), n = 1000)      The above model is $X_t = Z_t + 0.7 Z_{t-2} +0.1 Z_{t-5}$, where $\{ Z_t\}$ is white noise. Is it possible to estimate this model in R? I tried `arima` function and I haven't found the option to specify the specific lag. If it is not possible in R, is there any specific reason?


title: Multivariate Gaussian with 3 Partitions

text: Given that we partition a Gaussian random vector $\textbf{x}$ into three groups, $\textbf{x}_a$, $\textbf{x}_b$, and $\textbf{x}_c$, with a corresponding partitioning of the mean vector $\mu$ and of the covariance matrix $\Sigma$ in the form $\mu = \begin{pmatrix} \mu_a\\ \mu_b\\ \mu_c\\ \end{pmatrix}, \enspace \Sigma=\begin{pmatrix} \Sigma_{aa} & \Sigma_{ab} & \Sigma_{ac} \\ \Sigma_{ba} & \Sigma_{bb} & \Sigma_{bc} \\ \Sigma_{ca} & \Sigma_{cb} & \Sigma_{cc} \\ \end{pmatrix} $ What is the expression for the conditional distribution $P(x_a|x_b)$ in which $x_c$ has been marginalized out? (Note the 3 groups are also multivariate)


title: Tail dependence and copulas

text: I have been given this formula for upper tail dependence and read that tail dependence depends on the copula and not the marginals: $$ \lambda_U = \lim_{a \to 1} \Pr[Y>F_Y^{-1}(a)\mid X>F_X^{-1}(a)] . $$ Would the inverse functions $F_X^{-1}(a)$ and $F_Y^{-1}(a)$ be copula functions ? I have an intuitive understanding of what couplas are, but my knowledge of stats is quite weak, I would appreciate if you provide simplistic answers where possible.


title: Appropriate Test? Is there a "t-test" for ratios when large number of data points and multiple runs?

text: N.B. I only have a very basic statistics background since I am in grade 9 in high school. Any help would be GREATLY appreciated. My Experiment: -Part 1: 2 variables (A and B) were each sampled in my homemade wind tunnel simultaneously (1000/s) for 10 seconds=10,000 data points in a run (Run 1) for Configuration 1 at Angle 1 -this was repeated 10 times=10 Runs -Next, a ratio (Mean A/Mean B) was calculated from the 2 variables for each run (Configuration 1 Angle 1 Ratios 1-10) -Part 2: the same procedure as Part 1 was performed for 7 different angles (i.e. Configuration 1, angle 1 for 10 runs, Configuration 1, angle 2 for 10 runs etc in the wind tunnel) Part 3: The same procedure as Part 1 and Part 2 was conducted for various configurations (Configuration 2, Angles 1-7 each with Runs 1-10, Configuration 3 etc.) -Now I want to compare whether these A/B Ratios between a) Configurations (i.e. Configuration 1 Angle 1 vs Configuration 2 Angle 1) and between b) Angles (i.e. Configuration 1 Angle 1 vs. Configuration 1 Angle 2) are "significant" (i.e. different from each other) using some kind of t-test.  I only know t-tests but not sure if I can use them here because of 1) ratios 2) large number of data points. I have SO MUCH data that I'm very confused. I hope this makes sense. It may seem like a simple problem to you, but if anyone can help I would REALLY appreciate it!


title: How do I interpret a significant interaction effect that has a low effect size?

text: I am comparing measures of a DV x with two group factors (2-way ANOVA) - y (with 2 levels) and k (with 3 levels). Both present significant main effects (p<.001) and high partial eta squared (above 0.6). The interaction effect is also significant (p<0.001) but the partial eta squared is 0.029. Post-hoc testing for the interaction of y*k yields significant differences on the DV for all comparisons. In this case, shouldn't the effect size be greater? How do I interpret this? Thank you for your help! P.S. - Note: (I have a big sample; N>600)


title: Combining results of two samples to find combined standard deviation

text: I was making an attempt to solve the following problem sets, however, my answer **15.913** seemed quite different from that in the answer key **16.84**. I have checked my answer a few times and re-doing the question as well, but I was consistently pulled back to the same outcome of **15.913**. Appreciate some guidance and advice please. **Question** > A recent study of home technologies reported the number of hours (to the > nearest hour) of personal computer usage per week for a sample of 10 > persons. Excluded from the study were people who worked out of their home > and used the computer as a part of their work. > > 8, 5, 6, 8, 6, 1, 5, 6, 9, 4 > > In a study with sample size of 30 people who worked out of their home and > used the computer as a part of their work, the sample mean and standard > deviation of the number of hours personal computer usage per week are 42.0 > and 6.4, respectively. If we wants to combine the results of these two > samples, find the standard deviation of the number of hours personal > computer usage per week of in the combined sample. **My Attempt** > **New Sample** > > n=10, sample mean = 5.8, Sum(xi) = 58, Sum(xi^2) = 384 > > **Old Sample** > > n=30, sample mean = 42, s= 6.4, Sum(xi) = (30)(40) = 1260, Sum(xi^2) = > (30)(40^2) = 52920 > > **Combined Sample** > > n=40 > > sample mean = (1260+58)/40 = 32.95 > > Sum(xi) = 1260 + 58 = 1318 > > Sum(xi^2) = 384 + 52920 = 53304 > > s = SqRt{[53304 - (40)(32.95^2)]/[40-1]} = 15.91314 The answer offered in the answer key **16.84** seemed quite different from my **15.91314**. Appreciate any advice if there might be some steps which I may be doing wrongly.


title: One Bayes classifier to rule them all?

text: I listen for tweets matching keywords. When I receive a tweet, it may or may not be part of the set of things I'm interested about. Think about the Gemini Awards, a show similar to the Oscars, in Canada. #Gemini also happens to be an astrological sign, and many, many, many people talk astrology on Twitter. Now, imagine I have thousands of those TV show keywords I'm watching for on Twitter / Facebook / et al. I thought of using Naive Bayes Classifier(s) to determine if the interaction's part of my result set or not. I want to know if I should create / maintain one classifier per TV show, or one global classifier. On the one hand, having one classifier per TV show would mean:   * I can parallelize **training** per TV show (membership tests can be parallelized for both approaches);   * Training is very fine grained;   * The classifiers answer one of two categories: in or out;   * A message that mentions two or more TV shows would be classified "in" by multiple classifiers;   * When something is really outside of what we want, we have to train all N classifiers on this;   * New TV shows can be trained very quickly into their resulting domain. On the other hand, one classifiers means:   * Only one show at a time can be trained;   * Only one kind of answer can be returned per message: show X or no shows;   * Adding a new TV show would imply a longer training period, to rebalance the common words into the different categories;   * Training for irrelevancy is very easy, since there is only one classifier to talk about. I'm wondering about general performance and storage characteristics for one or more Bayes Classifier. Maybe I'm also taking a completely wrong approach and should be using a different algorithm to determine membership?


title: Variant of discriminant analysis for known multiple independent classifications?

text: I have a large data set: over 100,000 data points, each with 60 dimensions. I want to display the data in 2D to visibly maximize the separation between classes, which I know for each point. I asked a similar question before, and Linear Discriminant Analysis (LDA) was suggested -- a solid answer. Now, besides knowing that each data point falls into a certain class (say 1/2/3), each of these data points also has a second type of classification (say A/B/C). I would like to visualize how well these points cluster by one class, and then see how similar the substructures of the clusters are when each cluster is broken up based on the other class. For example, if we choose the numerical class for the "outer" clustering and the alphabetical class for the "inner" clustering, I would want to find the 2D representation of the data which makes it look maximally like this: ![A sample classification diagram](http://i.stack.imgur.com/T0Z9v.png) Bonus points if we can find a way to guarantee that the A/B/C subclasses are all oriented the same way with respect to each other, for each main class. Note, however, that I would also expect to be able to perform the clustering in the opposite order (with the alphabetical class forming the "outer" clusters and the numerical class forming the "inner" clusters) because the classifications are independent. The problem with LDA in this context is that, when considering only one class at a time, it maximizes the between-class separation at the expense of potentially blurring the boundaries of the other class in the resulting visualization, because it doesn't know about that the other class exists. So, when I visualize LDA on one classification, the distinctions for the second classification are not as clear as they could be. However, if I do a 9-way LDA (with each class being a combination of the two classifications, like 1A), I see 9 clusters, but they are not grouped together in any easily visualizable way; this is expected, however, because the classifications are independent (not to mention the fact that there are now 8 LD dimensions). I think about it this way. Say I start at center of the figure above. If I know which instance of the numerical class I want, I know which direction to move in: down for 1, up-left for 2, up-right for 3. The goal is to show that, no matter which numerical cluster I went to, I have the same directional options if I want to then move from that cluster's center to an alphabetical cluster. In each case, it's up-left for A, down for B, and up-right for C. And of course I would expect the same to be true if the classifications were performed in reverse (although this would require a different figure). Are there any standard techniques for doing this kind of hierarchical cluster visualization? If not, any advice on how to get started? **UPDATE:** The best solution I've found so far is to perform two separate LDAs: one for the alphabetical class and one for the numerical class. Then I can visualize the first dimension of each in the same graph, like this: ![Best-so-far classification diagram](http://i.imgur.com/QO7uY.png) While this gives me something like the picture I was looking for, I'm not sure that it proves what I set out to prove, which is that the data represent these two types of classifications independently. While this graph makes it look like the data bear this out, to make the argument robust I'll have to show that the alphabetical LD and the numerical LD are uncorrelated.


title: Meta analysis: How to assess the quality of each of the individual studies in a systematic way?

text: I have a meta analysis at hand. I have reported on the quality of the included studies in a narrative form. But after I was requested about "quality assessment". I guess it might be something else, a statistical model for example, something like weighting. I googled it but to my surprise, no good links were found in the first half an hour of searching. So I am asking here. Plus from the vague comments of the reviewers, it appeared to me that there is a "quality assessment" of the meta-analysis itself too. I mean something analytical and other than PRISMA. Do you have any advice or link regarding which tests or analyses should I do for this quality assessment purpose?


title: Is it permissible to run post-hoc tests after a non-significant ANOVA?

text: We ran a 3-way ANOVA getting a non-significant result, but in doing so we saw that one condition looked like it had a strong effect. We tried a Bonferroni- corrected t-test of just that group, and it turned out significant. Is it valid to consider this? Is there some other correction needed to deal with the fact that we tested that group only after we saw that it "looked" like an effect? (And does the answer to this sort of question depend on between vs. within factors? We have a mix. I ask this because SPSS doesn't seem to give the option for post hoc tests under the GLM procedure except for between-factor effects. Not the best way to learn, but it's a clue, nonetheless.)


title: How do I calculate uplift in a marketing campaign?

text: I have a campaign that had two interventions A and B that were run in 3 markets each. I, also, had a control in 15 markets where no interventions were run. For the markets where intervention A was run we had 5 units sold in the time period before intervention and 5 units during intervention. For the markets where intervention B was run we had 10 units sold in the time period before intervention and 15 units during intervention. Now, the control markets had 100 units prior to intervention and 200 during. So...                           Prior         During     Control      100            200     Group A        5              5     Group B       10             15      The marked increase in the Control was due to seasonality. Now, I want to calculate uplift (or downlift with this example data) by saying Group A should have expected to have 10 unit sales during intervention thus                (5-10)/5  = -1      While Group B has an uplift of               (15-20)/10 = -0.5       One of my colleagues wants to divide by the control values to account for the seasonality while I think we should subtract by the expected unit sales to account for the seasonality. Sorry, I cannot give more details.


title: Repeated measures ANOVA: Mauchly's test undefined

text: I'm doing a two-way between-within ANOVA in SPSS. I have two groups with 9 subjects each (so total = 18), and 24 levels of one repeated measure. I understand why Mauchly's test of Sphericity has no meaning when there are are only 2 levels of a repeated measures factor, but I notice (using General Linear Model.....repeated measures in SPSS) that Mauchly's test of Sphericity also appears to be undefined (or gives the useless output of Mauchly's W = '.0' , p = '.') when the number of levels of a repeated measure is equal to or greater than the number of cases (subjects). In these instances, even though Mauchly's statistic is not calculated, Greenhouse-Geisser, Huynh-Feldt, and Lower-Bound Epsilon values are calculated. I would be really happy if someone could provide some insight on why Mauchly's statistic is not calculated in these cases and what should be done to assess sphericity in the absence of Mauchly's statistic.


title: taking the log of (Only) some variables in a regression

text: My question is as follows: Ratio=Ratios + log (numbers) +dummy variables + volatility I have this type of regression in a paper published by the Federal Reserve Bank.Can someone tell me why we took the log of variables that represent numbers (such as size,and fees) whereas for the fractions such as (market value of equity/book value of equity) we didnd't take log. Also we didn't take it for the dummies .Finally the volatility which is a number was without log . Best, Wesso


title: Difference between Factorization machines and Matrix Factorization?

text: I came across the term Factorization Machines in recommender systems. I know what Matrix Factorization is for recommender systems but never heard of Factorization Machines. So what's the difference?


title: Can you help me interpret GARCH results?

text: Can you please explain what the next GARCH results mean? where in these results can i know how well does the model forecast? ![garch](http://i.stack.imgur.com/lGFwf.jpg)


title: Bayesian Inference for Poisson process

text: I have a historical dataset of about 1500 events representing structure fires in the Bay Area over the last few years. I'm building a small dashboard that displays these events by month, and would like to build in a prediction that shows the number of events expected to occur for the rest of the month for comparison's sake (does this look like a busier or slower month). The data is somewhat cyclical by season: higher in December/January and May/June, and lower in the Spring and Fall. Naturally this seems like a Poisson process. Currently I'm doing a simple estimation by determining the average number of events in the historical data during the current month, and then multiplying it by the fraction of the month remaining (so if historical average is 30 events in December and it's the 15th, predict $30\times\frac{16}{31}$ as the number of events for the rest of the month. This is probably good enough, but academic curiosity has taken hold and I'd like to go further. Ideally, I'd like to take in to account recent trends, so that if the recent history is noticeably busier than historically, the projection takes that into account. It seems like using Bayes' rule would be useful here, where the prior is based on the historical average, evidence is the recently observed data (say 30-60 days worth), and then solve or analyze for the more likely λ. Is this the right, or a good, way to go? If so, can you point me towards setting up the equation, or what other information do I need to figure out first?


title: Parameters of energy function for TSP

text: I am reading this paper by Hopfield et al. On page six, the authors defined the energy function of the Traveling-Salesman-Problem (TSP) mapped onto an artificial neural network as follows: $$E = \frac{A}{2} \sum_X \sum_i \sum_{j \ne i} V_{Xi} V_{Xj} + \frac{B}{2} \sum_i \sum_X \sum_{X\ne Y} V_{Xi} V_{Yi} + \frac{C}{2} (\sum_X \sum_i V_{Xi} - n)^2 + \frac{1}{2D} \sum_X \sum_{Y\ne X} \sum_i d_{XY} V_{Xi} (V_{Y, i+1} + V_{Y, i-1})$$ Here the first three terms maintain correctness of the function while the final term calculate actual distance of the tour. Later on page seven, the authors picked the values of A, B, C and D for simulation. I don't understand on what basis they choose those values. Any help?


title: Deep belief networks or Deep Boltzmann Machines?

text: I'm kinda confused. Is there a difference between Deep belief networks and Deep Boltzmann Machines? Are they different stuff or same thing?! If so, what's the difference?


title: What is best for comparing variability for environmental time series data: repeated-measures ANOVA or linear mixed models?

text: I've collected environmental data on the dust level using direct-reading instruments. The outputs from the instruments are concentrations recorded at each minute over 1 hour. The data collected at many locations and locations are grouped into Groups and Department. The data looks like this in wide format               Depart1  Group1  Location1:  3, 3, 4, 3.5, 4, 6, etc...     Depart1  Group1  Location2:  2, 7,5, 6, 3, 5, 5, etc..     Depart2  Group2  Location3:  8, 4, 3, 6, 6, 3, 6, etc..      I want to find the within and between variability among the Groups. But for some groups, I have data at 3 locations per group and sometimes I only have data at 1 location per group. Thus, I have uneven datasets for Groups. I am wondering what is the best statistical model to answer the variability question. I've looked up ANOVA for repeated measurements but how do you deal with the uneven grouping? Would a linear mixed model with time being considered as a fixed effect and Group and Location being random effects work here? At the end, I am only interested in within and between variability for the Group not the Location.


title: How does one show that there is no unbiased estimator of $\lambda^{-1}$ for a Poisson distribution with mean $\lambda$?

text: Suppose that $ X_{0},X_{1},\ldots,X_{n} $ are i.i.d. random variables that follow the Poisson distribution with mean $ \lambda $. How can I prove that there is no unbiased estimator of the quantity $ \dfrac{1}{\lambda} $?


title: Merging several classifiers - does it make sense?

text: all! I have a bit theoretical question - imaging, you trained several classifiers of the same type on the several datasets, that are in some sense close to one another (for example, they have intersection, or the ranges of parameters, that are being estimated, are close to each other) Does it make any sense to merge the parameters of these classifiers and in this way to build a new classifier? Does something like this (a method to derive a new classifier from aset of already existing classifiers) exist already?


title: Least squared regression where the coefficients switch signs upon the addition of new variable

text: I am conducting a least square regression using the python library numpy. When I run OLS (ordinary least square regression) over just the variable IE6 I get this output (with the key takeaway being that the coefficient is negative):               =============================================================================     variable     coefficient     std. Error      t-statistic     prob.     ==============================================================================     const          -0.632626      0.184070     -3.436882      0.001209     IE6            -11.845141      8.695507     -1.362214      0.179360     ==============================================================================     Models stats                         Residual stats     ==============================================================================     R-squared             0.036488         Durbin-Watson stat   1.635514     Adjusted R-squared    0.016825         Omnibus stat         12.433490     F-statistic           1.855627         Prob(Omnibus stat)   0.001996     Prob (F-statistic)    0.179360          JB stat              43.784881     Log likelihood       -47.743508         Prob(JB)             0.000000     AIC criterion         1.950726         Skew                 0.052388     BIC criterion         2.026484         Kurtosis             7.538025     ==============================================================================      Yet when I add two additional variables percent and percent_m, the sign of the coefficient switches to positive. Apparently, I don't understand the mechanics of OLS well enough, because this switch perplexes me. Here is the output of the second regression:               ==============================================================================     variable     coefficient     std. Error      t-statistic     prob.     ==============================================================================     const          -1.511714      0.417273     -3.622840      0.000725     percent           1.852263      0.710854      2.605689      0.012313     percent_m          -2.657285      1.233789     -2.153760      0.036533     IE6             9.625211      6.499429      1.480932      0.145442     ==============================================================================     Models stats                         Residual stats     ==============================================================================     R-squared             0.282134         Durbin-Watson stat   1.896240     Adjusted R-squared    0.219711         Omnibus stat         8.682105     F-statistic           4.519698         Prob(Omnibus stat)   0.013023     Prob (F-statistic)    0.003660          JB stat              14.649562     Log likelihood       -40.238818         Prob(JB)             0.000659     AIC criterion         1.774071         Skew                 0.366255     BIC criterion         1.963466         Kurtosis             5.521377     ==============================================================================


title: Longitudinal models in R and WINBUGS or JAGS

text: I've tried to use R to fit some longitudinal models, mostly via `lmer` and `nlme` packages. However, it seems that many standard models are lacking, such as antedependence models or factor analytic models for covariance matrices. These models are readily available in SAS. Would anyone recommend other packages for the job in R? I don't really care if I use to work in a frequentist or bayesian world as long as I have more modeling flexibility. I would also be interested in doing that in WINBUGS/JAGS.


title: Determining statistical significance of a repeated measures t-test

text: ### Scenario: An industrial/organizational psychologist is interested in determining whether adding 15-minute breaks increases worker productivity. She selects a sample $n$ and measures productivity (on a continuous scale) before and after introducing the intervention. The researcher runs a repeated measures t-test. ### Question   * How can I work out whether the intervention is effective?


title: How to smooth time-series NDVI data using polynomial regression

text: I have a time-series NDVI image. The image has 26 bands. 26 bands mean that images taken in 8-day time interval and counted in Julian days (97 to 297). For example; first band of the image is NDVI image at Julian date of 97, second is 105 etc until Julian day of 297. The image size is 12618, 4144. I have chosen 30 point samples from each land cover (30 x 5 LC=150) from all band of time series NDVI image using stratified random sampling method. The image below is samples chosen from desert area of NDVI images.Totally I have chosen 50 samples. First column is Julian date (day number when the image taken). Now I'd like to do smooth NDVI image by polynomial regression. NDVI image response to Julian date to represent the seasonal changes in NDVI as a function of Julian day. I expect the output image to have different values in each land cover type and the same dimension with original one. ![enter image description here](http://i.stack.imgur.com/EFQh8.png). When I tried to do it in Sigmaplot, response variable should be one (50 samples), expanatory variable (Julian days) can be multiple variables. But in my case, this is reverse. So I would like to ask: Can I use all the samples of NDVI at once to get regression coefficient to predict NDVI value of desert? If so, how should I prepare the data for Minitab. Or do I need to run the regression for each sample point (150 samples) one by one get their regression coefficients? Is there any method to simplify these processes to save time? Thanks a lot


title: Comparing different deep learning models?

text: Does anyone know a paper that describes the differences and compares the different deep learning architectures? like Stacked autoencoders, deep believe networks, maxout networks ... etc.


title: Which component to choose when variable load heavily on more than one component

text: I am doing factor analysis with Principal component extraction method. Below is the Structure matrix Extraction method: PCA Rotation method: Promax with kaiser normalization ....................................................................COMPONENT............................................................. ............................................................. 1 ............ 2 ......... 3 .............4 ..........................................   Effeciency class ............................ .893.......... .236........ .139......... .876 Pre sales support......................... .839 ......... -.224 ........ .811.......... .326 I have two variable Eff_class and Pre_sales support and four component..Now variable seems to load alomost equally on two components. How to decide which component to choose.


title: Agreement statistics for confusion matrix with nominal categories

text: I have a decision support system that recommends Treatments (`T1` to `T8`) to patients, where `T`'s are nominal (i.e. no ordering or rating between them). I compare the system recommendations with the recorded treatments in a confusion matrix, where the columns list computer recommendations and the rows list recorded treatments:               cmg <- matrix(c(1639, 116, 49, 35, 138, 0, 0, 236,                      150, 274, 27, 21,  28, 0, 0,  73,                       22,  24, 58,  9,  94, 0, 0,  30,                       33,  27, 31, 21, 146, 0, 0,  49,                       14,   9,  5,  1,  49, 0, 0,  22,                        1,   0,  1,  1,   7, 0, 0,   6,                       11,   0,  0,  1,  14, 0, 0,  21,                      201,  11,  8,  5,  49, 0, 0, 253),                    ncol=8,dimnames = rep(list(c("T1","T2","T3","T4","T5", "T6","T7", "T8")),2))      Based on `cmg`, I need to provide some **agreement statistics** and report how much the system agrees with recorded treatments. I checked this link, which suggests finding overall and specific agreements for nominal data. Here, specific agreements are calculated by collapsing the matrix to a binary confusion matrix [2*2] for each `Ti`. Also, they suggest calculating Cohen's Kappa both for `cmg` and for each collapsed matrix separately. I do these as below:               require(psych)     # Overall agreement     overall_agg <- sum(diag(cmg))/sum(cmg)          # Overall Cohen's Kappa for cmg     unweighted_kappa <- cohen.kappa( cmg, n.obs=sum(cmg) )          # initialise containers     spec_agr_guideline <- list()             collapsed_mat_guideline <- list()      unweighted_kappa_psych <- list()          # loop through all treatments         for (i in seq(1,nrow(cmg)) ) {       # Specific agreements       spec_agr_guideline [i] <- 2*cmg[i,i] / (sum(cmg[i,]) + sum (cmg[,i]) )       # Collapsed positive agreement confusion matrices per treatment       collapsed_mat_guideline[[i]] <- matrix(c(cmg[i,i],             sum(cmg[i,])-cmg[i,i],                                              sum(cmg[,i])-cmg[i,i],  sum(cmg)-sum(cmg[i,])-sum(cmg[,i])+cmg[i,i]),                                             ncol=2)       # Calculate unweighted Cohen's Kappa per collapsed (binary) confusion amtrix       unweighted_kappa_psych[[i]] <- cohen.kappa( collapsed_mat_guideline[[i]], n.obs=sum(collapsed_mat_guideline[[i]]) )     }      **Questions:** **1)** Do I report the magnitudes of the Cohen's Kappas, or their p-values or something else? How can I interpret the Kappa results?E.g. for `cmg`, I get the output below:               unweighted_kappa     Cohen Kappa and Weighted Kappa correlation coefficients and confidence boundaries                       lower estimate upper     unweighted kappa  0.34     0.36  0.38      **2)** It is said that the usefulness of the Kappa values are open to debate. Are there any other agreement stats that I need to report in my (nominal data) case? For instance, in the same link, they suggest Stuart-Maxwell, McNemar tests, etc..


title: Effect size for indirect effect with a binary outcome

text: Is it possible to obtain the effect size of a mediation with a binary outcome? I am looking to report an effect size specifically for the indirect effect. I have already reported the effect sizes of the individual effects. For example: x = continuous IV - How connected to nature M = continuous mediator - How motivated Y = binary DV - Recycled yes or no I am using the PROCESS macro from Hayes to calculate the bootstrapped indirect effect. I can test the significance of the indirect effect using the PROCESS macro (by looking at the bootstrapped bias-corrected confidence intervals), but the program does not calculate an effect size for the indirect effect when the outcome is binary (it does when the outcome is not binary). Is there a formula for calculating the indirect effect effect size for logistic regression mediation? References would be a great help.


title: Proof of consistency of Maximum Likelihood Estimator(MLE)

text: I would appreciate some help comprehending a logical step in the proof below about the consistency of MLE. It comes directly from Introduction to Mathematical Statistics by Hogg and Craig and it is slightly different than the standard intuitive one that makes use of the Weak Law of Large Numbers. So here goes: Assume that $\hat{\theta_n}$ solves the estimating equation $\frac{\partial l(\theta)}{\partial \theta}=0$. We also assume the usual regularity conditions. Denote $\theta_0$ the true parameter which by assumption is an interior point of some set $\Omega$, Then $\hat{\theta_n} \xrightarrow{P} \theta_0$ **Proof** Let $\mathbf{X}=(x_1,x_2, \ldots, {x_n})$, the vector of observations. Since $\theta_0$ is an interior point in $\Omega$ , $(\theta_0 -a, \theta_0 +a) \subset \Omega $ for some $a >0$. Define $S_n$ to be the event $$S_n= \{ \mathbf {X} : l(\theta_0 ; \mathbf{X}) > l(\theta_0 -a ; \mathbf{X}) \} \cap \{ \mathbf{X}: l(\theta_0; \mathbf{X}) > l( \theta_0 +a ;\mathbf{X}) \} $$ **(The log- likelihood is maximised on the true parameter, all good so far)** But on $S_n$, $l \left(\theta \right)$ has a local maximum, $\hat{\theta_n}$ such that $\theta_0-a<\hat{\theta_n}<\theta_0+a$ and $l^{\prime} \left(\hat{\theta_n} \right)=0$ That is $$ S_n \subset \{ \mathbf{X}: | \hat{ \theta_{n} } \left( \mathbf{X} \right) -\theta_{0} | < a \} \cap \{ \mathbf{X}: l^{ \prime} \left( \hat{\theta_n} \left( \mathbf{X} \right) \right) =0 \} $$ It is precisely at this point that I find their proof a little obscure. How come they consider $S_n$ to be a subset of that other set? Their explanation is unclear. Of course the proof is not complete at this point but if I have this clarified, I can take it from there. Thank you in advance.


title: Cronbach's Alpha low- significant results

text: One of my variables has a low Cronbach's Alpha. I continued to run an planned contrast to test my hypothesis and the results were significant. I understand that the results should be interpreted lightly. My question is however, what does it mean if the Cronbach's Alpha is low yet significant results have been found? I am writing an honours thesis so I am still learning. Thank you.


title: Maximum likelihood estimation of dlmModReg

text: I'm studying R package dlm. So far it seems very powerful and flexible package, with nice programming interfaces and good documentation. I've been able to successfully use dlmMLE and dlmModARMA to estimate the parameters of AR(1) process:               u <- arima.sim(list(ar = 0.3), 100)     fit <- dlmMLE(u, parm = c(0.5, sd(u)),                   build = function(x)                     dlmModARMA(ar = x[1], sigma2 = x[2]^2))     fit$par      Now I'm trying to use similar code to estimate the parameters of simple linear regression model:               r <- rnorm(100)     u <- -1*r + 0.5*rnorm(100)     fit <- dlmMLE(u, parm = c(0, 1),                   build = function(x)                     dlmModReg(x[1]*r, FALSE, dV = x[2]^2))     fit$par      I expect fit$par to be close to c(-1, 0.5), but I keep getting something like               [1] -0.0002118851  0.4884367070      The coefficient -1 is not estimated correctly. However, the strange thing is that the variance of the noise is returned correctly. I understand that max-likelihood estimation might fail given bad initial values, but I observed that the likelihood function returned by dlmLL is very flat in the first coordinate. So I wonder: can such model be estimated at all using dlm? I believe the model is "non-singular", however I'm not sure how the likelihood function is calculated inside the dlm. Any hint greatly appreciated.


title: Coupled Hidden Markov Models

text: Does anyone know where can I find a pseudo code or Matlab Code for CHMMs please? I am using categorical observations. Thanks.


title: Solving PCA with correlation matrix of a dataset and its singular value decomposition

text: Suppose I have a $d \times n$ matrix $\mathbf X$ (each entry point has $d$ dimensions) and after some manipulation of data (i.e. summarizing the data $\mathbf X$) I get its $d \times d$ symmetric, quadratic correlation matrix $\rho$ (defined by Pearson). Then by **SVD definition** we know that matrix $\mathbf X$ can be decomposed in three matrices: $\mathbf X = \mathbf U \mathbf \Sigma \mathbf V^\top$ Hence here is the question: I want to know if its correct to suppose that $\mathbf X \mathbf X^\top = \rho$ then doing some algebra we can get: $\rho = (\mathbf U \mathbf \Sigma \mathbf V^\top)(\mathbf V \mathbf \Sigma \mathbf U^\top)$ so $\mathbf X \mathbf X^\top = \rho = \mathbf U \mathbf \Sigma^2 \mathbf U^\top$ I want to know if this equivalence is correct, or in which cases is correct? Can anyone give an explanation about it?. Note I am using correlation Matrix instead of covariance matrix, (we know we can get eigenvalues and eigenvectors from a correlation matrix, and that would solve PCA)


title: Intepretation of crossvalidation result - cv.glm()

text: My logistic model has been suspicious due to enormous coefficients, so I tried to do a crossvalidation, and also do a crossvalidation of simplified model, to confirm the fact that the original model is overspecified, as James suggested. However, I don't know how to interpret the result (this is the model from the linked question):               > summary(m5)          Call:     glm(formula = cbind(ml, ad) ~ rok + obdobi + kraj + resid_usili2 +          rok:obdobi + rok:kraj + obdobi:kraj + kraj:resid_usili2 +          rok:obdobi:kraj, family = "quasibinomial")     [... see http://stats.stackexchange.com/q/48739/5509 for complete summary output ]          > cv.glm(na.omit(data.frame(orel, resid_usili2)), m5, K = 10)     $call     cv.glm(data = na.omit(data.frame(orel, resid_usili2)), glmfit = m5,          K = 10)          $K     [1] 10          $delta     [1] 0.2415355 0.2151626          $seed       [1]         403         271  1234892862 -1124595763  -489713400  1566924080   147612843       [8]  1879282918  -694084381  1171051622  2063023839 -1307030905  -477709428  1248673977      [15]  -746898494   420363755  -890078828   460552896  -758793089  -913500073  -882355605     [....]     Warning message:     glm.fit: algorithm did not converge      I guess the delta is the mean fitting error, but how to interpret it? Is it a good or bad fit? BTW, the algorithm did not converge, maybe due to the enormous coefficients (?) I tried a simplified model:               > summary(m)          Call:     glm(formula = cbind(ml, ad) ~ rok + obdobi + kraj, family = "quasibinomial")          Deviance Residuals:          Min       1Q   Median       3Q      Max       -2.7335  -1.2324  -0.1666   1.0866   3.1788            Coefficients:                   Estimate Std. Error t value Pr(>|t|)         (Intercept) -107.60761   48.06535  -2.239 0.025335 *       rok            0.05381    0.02393   2.249 0.024683 *       obdobinehn    -0.26962    0.10372  -2.599 0.009441 **      krajJHC        0.68869    0.27617   2.494 0.012761 *       krajJHM       -0.26607    0.28647  -0.929 0.353169         krajLBK       -1.11305    0.55165  -2.018 0.043828 *       krajMSK       -0.61390    0.37252  -1.648 0.099593 .       krajOLK       -0.49704    0.32935  -1.509 0.131501         krajPAK       -1.18444    0.35090  -3.375 0.000758 ***     krajPLK       -1.28668    0.44238  -2.909 0.003691 **      krajSTC        0.01872    0.27806   0.067 0.946322         krajULKV      -0.41950    0.61647  -0.680 0.496315         krajVYS       -1.17290    0.39733  -2.952 0.003213 **      krajZLK       -0.38170    0.36487  -1.046 0.295698         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           (Dispersion parameter for quasibinomial family taken to be 1.304775)              Null deviance: 2396.8  on 1343  degrees of freedom     Residual deviance: 2198.6  on 1330  degrees of freedom     AIC: NA          Number of Fisher Scoring iterations: 4      and it's crossvalidation:               > cv.glm(orel, m, K = 10)     $call     cv.glm(data = orel, glmfit = m, K = 10)          $K     [1] 10          $delta     [1] 0.2156313 0.2154078          $seed       [1]         403         526   300751243  -244464717  1066448079  1971573706 -1154513152       [8]   634841816 -1521293072 -1040655077   505710009  -323431793 -1218609191  1060964279      [15]  1349082996   -32847357 -1387496845   821178952  -971482876  1295018851  1380491861      Now it converged. But the delta seems more or less the same, despite of the fact that this model looks much more sane! I'm confused by the crossvalidation now... please give me a hint on how interpret it.


title: What does it mean to develop a proper implementation of the grammar of graphics?

text: I want to develop an implementation of the grammar of graphics in python, and I need to wrap my head around a few things. I am actually having some difficulties in understanding some of the book's broadest concepts. For what I understand, a correct implementation of the grammar means that a graphic is defined with its data characteristics instead of being composed (as it is done in the usual plotting libraries). When a plot command is issued, a graphic almost "renders itself". Is this interpretation correct and sufficient to build a proper implementation of the grammar or am I missing some deeper mathematical concepts? Edit: The first approach I have in mind is to define a class for a graph, and then have a way to compose into instances of this class: scales, coordinates, aesthetics etc. Also, should I enforce lazy data loading and calculations until the moment a plot is actually required?


title: Correlation between task type and measures

text: I'm pretty new to statistics in general; maybe someone can edit the title to provide a better idea of my question. The data I'm working with represents refined data extracted from browser log files gathered in a lab study. The probands solved three different types of tasks (let's call them a, b and c) and there's a seperate dataset for every task type. The datasets contain the same kind of data, e.g. the duration of a mousemove. How can I work out if there's a correlation between the task type and a variable in the data? Thanks in advance.


title: One sided confidence interval for the mean

text: I want to test the following hypothesis: $$ H_0: \mu \geq \mu_0 \\ H_a: \mu < \mu_0.$$ I use one sample t-test: $$t_0 = \frac{\bar{X} - \mu_0}{S/\sqrt{n}}.$$ Since under the alternative we expect that $\bar{X}$ is small, we should reject for small values of $t_0$. Thus, given that the size of the test is equal to $\alpha$, we reject if $$t_0 \leq -t_{1-\alpha}(n-1).$$ My aim is to construct a confidence interval for $\mu$. Given my results, can I construct a confidence interval for $\mu$ based on the fact that we don't reject if $$t_0 > -t_{1-\alpha}(n-1),$$ thus the confidence interval should be constructed by solving the following inequality $$ \frac{\bar{X} - \mu}{S/\sqrt{n}} > -t_{1-\alpha}(n-1)?$$ I am not sure that my way of reasoning is consistent since I use t-test (with $\mu_0$) in order to construct the rejection region and I use $\frac{\bar{X} - \mu}{S/\sqrt{n}}$ (with $\mu$ instead of $\mu_0$) in order to determine the confidence interval.


title: Difference in means in multiple-choice poll

text: Let's say that N randomly chosen persons where asked a question where the answer could be in either of X categories. For example, 500 persons where asked which of the top 5 political parties they support the most. Each person can give only one answer. How do I determine if the leading party which, for example, got 33 % is truly larger than the 2nd largest party that got 29 % (on a certain confidence level). Calculations like these solve my question when there are only two choices, but how about my case when there are several? http://www.stat.wmich.edu/s216/book/node85.html


title: How to derive Poisson distribution from gamma distribution?

text: Let $T_1, T_2, \dots$ be iid sequence of exponential random variables with parameter $\lambda$. The sum $S_n = T_1 + T_2 + \dots + T_n$ is a Gamma distribution. Now as I understand the Poisson distribution is defined by $N_t$ as follows: $$N_t = \max\{k: S_k \le t\}$$ How do I formally show that $N_t$ is a Poisson random variable? Any suggestions appreciated. I tried to work out a number of proofs but cannot get to the final equation. **References** http://en.wikipedia.org/wiki/Exponential_distribution http://en.wikipedia.org/wiki/Gamma_distribution http://en.wikipedia.org/wiki/Poisson_distribution


title: How to interpret mean of Silhouette plot?

text: Im trying to use silhouette plot to determine the number of cluster in my dataset. Given the dataset **Train** , i used the following matlab code               Train_data = full(Train);       Result = [];       for num_of_cluster = 1:20           centroid = kmeans(Train_data,num_of_cluster,'distance','sqeuclid');           s = silhouette(Train_data,centroid,'sqeuclid');           Result = [ Result; num_of_cluster mean(s)];       end       plot( Result(:,1),Result(:,2),'r*-.');`      The resultant plot is given below with xaxis as **number of cluster** and yaxis **mean of silhouette value**. How do i interpret this graph? How do i determine the number of cluster from this? ![enter image description here](http://i.stack.imgur.com/KSqC5.jpg)


title: How to interpret decreasing AIC but higher standard errors in model selection?

text: I've got a problem choosing the right model. I have a model with various variables (covariables and dummy variables). I was trying to find the best size for this model, so I first started by comparing different models with AIC. From this it followed, that the minimum AIC was reached when allowing all variables to stay in the model (with the whole bunch to interact with all dummies). When I compute the summary of the model, all effects are absolutely not significant and the standard errors are very high. I was a bit confused, when comparing the "best" (on AIC) model with a smaller model with any interaction. The smaller model had small standard errors and nice p-values... But the AIC is higher compared to the big model. What might be the problem? Overspecification? I really need help in this, because I have absolutely no idea how to handle this! Thanks a lot


title: How to test correlation between variables at two points in time?

text: I am currently analyzing some data for my Master Thesis. I measured the dependent and independent variable in 2012 and also 2013. Now I want to see whether X (independent) has an impact on y (dependent). I thought about doing a correlation. However, I am not sure about this because it is measured at two points in time. Could anybody help me out? Also, I use SPSS and would be very grateful if you could provide me with steps to be taken! Regards, Laura


title: Which algorithm to compute p-value of logrank test with three or more groups is best?

text: There seem to be two different algorithms for comparing three or more survival curves using the logrank test. Algorithm A. Found in books by Altman and Machin. Computed by GraphPad Prism. This method uses an algorithm that is easy to understand, computing chi-square from the discrepancy between observed and expected numbers of deaths (with df=# of groups minus 1). Basically:               Chi2 = SumForAllCurves[(Oi - Ei)^2 / Ei]      Algorithm B. Computed by SAS and SPSS and NCSS. Also computes chi-square (also with df=# groups minus 1), but using a more complicated equation. Basically:                   Chi2 = U'*(V^-1)*U, where         U, V - vector and covariance matrix defined at p.6 of the same pdf      Here are sample data as an Excel file. Method A computes chi-square = 4.094; df = 2; P = 0.1291   Method B computes chi-square = 4.844; df = 2; P = 0.0888 Why the difference? Is one a simplification? Do they make different assumptions? Is one outmoded? (Note that both are variations of the logrank method that gives equal weight to deaths at any time. You can get different results by weighting deaths differently using Gehan-Wilcoxon, or Tarone-Ware or Peto-Peto..., but those choices don't explain the distinction I've seen.)


title: When to use time dummies in multiple regression

text: When is it appropriate to use time dummies in multiple regression analysis? I am using pooled ols, random effects model, and fixed effects model. I have a period of 3 years. I don´t know whether it is corrected but I think that time dummies should be included when the period under investigation is large, at least larger than 3 years. Then by using time dummies you can correct for time trends that are not observed by the model? Am I right or? And is it always included in fixed effects model?


title: Fitting Noisy Data

text: I've written the following code to fit $$f=a*sin(b*x+c)$$ for `x=linspace(0,4*pi,100)`. I know the form of the required function so I can use `fminsearch` to find the best coefficients. To be able to use `fminsearch` it's need to write a matlab function that for given values of `a`, `b` and `c`, calculates the square of the disparity. Its first argument is a vector of coefficients, its last argument is a vector of the given data values. The other arguments are vectors of free variables.                function fittingg      x=linspace(0,4*pi,100);      y=sin(x);      bestcoeffs=fminsearch(@fun,[1 1 1],[],x,y);      yfit=bestcoeffs(1)*sin(bestcoeffs(2)*x +  bestcoeffs(3));      y_noise=awgn(y,0.5);      figure      plot(x,y,x,yfit,'r','LineWidth',2);      hold on      plot(x,y_noise,'o','MarkerEdgeColor','k',...                               'MarkerFaceColor','b',...                         'MarkerSize',8);           function outt=fun(coeff,X,Y)      a = coeff(1);      b = coeff(2);      c = coeff(3);      Y_fun = a.*sin(b .* X +c);      DIFF = Y_fun - Y;       SQ_DIFF = DIFF.^2;      outt = sum(SQ_DIFF);      ![output](http://i.stack.imgur.com/VRilu.png) I have been also trying to fit and plot that `y_noise=awgn(y,0.5)` but I didn't reach a conclusion. Could you please help me with this?


title: Why is a Mixed model a non-linear statistical model

text: I was wondering if anyone could elaborate more on this statement which I came across whilst reading a book on non-linear mixed effects model. I know that you can have linear, generalised linear and non-linear mixed effect models. The book says "mixed models are non-linear statistical models, due mainly to the presence of variance parameters". How do I explain this to someone who understands very little statistics? For instance, how can a linear mixed effect model be non-linear? Thanks


title: Non-linear modelling with several variables including a categorical variable

text: I am trying to model some data regarding a predator prey interaction experiment (n=26). Predation rate is my response variable and I have 4 explanatory variables: predator density (1,2,3,4 5), predator size, prey density (5,10,15,20,25,30) and prey type (3 categories). I started with several linear models (GLM) and found (as expected) that prey and predator density were non-linearly related to predation rates. If I use a log transformation on these variables I get really nice curves and an adjusted $R^{2}$ of 0.82, but it is not really the right approach for modelling non- linear relationships.               model <-glm(rates ~ log(pred) + log (prey) + type)      Therefore I switched to non-linear least square regression (`nls`). I have several predator-prey models based on existing ecological literature e.g.:                ### Holling's type II functional response          model1 <- nls(rates ~ (a * prey)/(1 + b * prey),     start = list(a = 0.27,b = 0.13), trace = TRUE)          ### Beddington-DeAngelis functional response          model2 <- nls(rates ~ (a*prey)/(1+ (b * prey) + c * (pred -1 )),     start = list(a=0.22451, b=-0.18938, c=1.06941), trace=TRUE, subset=I1)       These models work perfectly, but now I want to add prey type as well. In the linear models prey type was the most important variable so I don't want to leave it out. I understand that you can't add categorical variables in nls, so I thought I try a generalized additive model (GAM). The problem with the gam models is that the smoothers (both spline and loess) don't work on both variables because there are only a very restricted number of values for prey density and predator density. I can manage to get a model with a single variable smoothed using loess. But for two variables it is simply not working. The spline function does not work at all because I have so few values (5) for my variables (see model 4).               model3 <- gam(rates~ lo(pred, span=0.9)+prey)     ## this one is actually working but does not include a smoother for prey.          model4 <- gam(rates~ s(pred)+prey)     ## this one gives problems:      A term has fewer unique covariate combinations than specified maximum degrees of freedom      My question is: are there any other possibilities to model data with 2 non- linear related variables in which I can also include a categorical variable. I would prefer to use `nls` (`model2`) with for example different intercepts for each category but I'm not sure how to get this sorted, if it is possible at all. The dataset is too small to split it up into the three categories, moreover, one of the categories only contains 5 data points. Any help would be really appreciated.


title: Random Variable

text: Three components are randomly sampled, one at a time, from a large lot. As each component is selected, it is tested. If it passes the test, a success (S) occurs; if it fails the test, a failure (F) occurs. Assume that 80% of the components in the lot will succeed in passing the test. Let X represent the number of successes among the three sampled components. What are the possible values for X? And There Probabilities ?


title: Does Neyman-Pearson Lemma consider the case when the likelihood ratio equals the critical value?

text: Here are three different versions of Neyman-Pearson lemma. They differ in that the first two (books) ignore the case when the likelihood ratio equals the critical value, while the last one (Wikipedia) doesn't. I feel not ignoring is natural, because ignoring will leave the set where the likelihood ratio equals the critical value to be undecided. So I was wondering why some ignore while some don't? Thanks and regards!   1. In Casella and Berger's Statistical Inference ![enter image description here](http://i.stack.imgur.com/t0oEy.png)   2. In Bickle and Doksum's Mathematical Statistics Vol 1 ![enter image description here](http://i.stack.imgur.com/Ep7dQ.png) ![enter image description here](http://i.stack.imgur.com/YEyk4.png) Especially in (c), the event excludes the case when the likelihood ratio $L(x, \theta_0, \theta_1)$ equals the critical value $k$. So how can we say the $\phi_k$ ahd $\phi$ are equivalent up to a set of probability zero?   3. In Wikipedia: > when performing a hypothesis test between two point hypotheses $H_0: θ = > θ_0$ and $H_1: θ = θ_1$, then the likelihood-ratio test which rejects $H_0$ > in favour of $H_1$ when $$ \Lambda(x)=\frac{ L( \theta _0 \mid x)}{ L > (\theta _1 \mid x)} \leq \eta $$ where $$ P(\Lambda(X)\leq \eta\mid > H_0)=\alpha $$ is the most powerful test of size $α$ for a threshold $η$.


title: What's the difference between likelihood and confidence in claim being true

text: I'm reading the IPCC report on climate change from 2007. In their uncertainity guide they make a distinction between _likelihood_ and _levels of confidence_. What's the difference between the terms? http://www.ipcc.ch/pdf/assessment- report/ar4/wg1/ar4-uncertaintyguidancenote.pdf


title: Getting p-values for "multinorm" in R (nnet package)

text: How do I get p-values using the `multinorm` function of `nnet` package in `R`? I have a dataset which consists of “Pathology scores” (Absent, Mild, Severe) as outcome variable, and two main effects: Age (two factors: twenty / thirty days) and Treatment Group (four factors: infected without ATB; infected + ATB1; infected + ATB2; infected + ATB3). First I tried to fit an ordinal regression model, which seems more appropriate given the characteristics of my dependent variable (ordinal). However, the assumption of odds proportionality was severely violated (graphically), which prompted me to use a multinomial model instead, using the `nnet` package. First I chose the outcome level that I need to use as baseline category:               Data$Path <- relevel(Data$Path, ref = "Absent")      Then, I needed to set baseline categories for the independent variables:               Data$Age <- relevel(Data$Age, ref = "Twenty")     Data$Treat <- relevel(Data$Treat, ref="infected without ATB")       The model:               test <- multinom(Path ~ Treat + Age, data = Data)      # weights:  18 (10 variable)      initial value 128.537638      iter 10 value 80.623608      final  value 80.619911      converged      The output:               Coefficients:              (Intercept)   infected+ATB1   infected+ATB2   infected+ATB3    AgeThirty     Moderate   -2.238106   -1.1738540      -1.709608       -1.599301        2.684677     Severe     -1.544361   -0.8696531      -2.991307       -1.506709        1.810771          Std. Errors:              (Intercept)    infected+ATB1   infected+ATB2   infected+ATB3    AgeThirty     Moderate   0.7880046    0.8430368       0.7731359       0.7718480        0.8150993     Severe     0.6110903    0.7574311       1.1486203       0.7504781        0.6607360          Residual Deviance: 161.2398     AIC: 181.2398      For a while, I could not find a way to get the $p$-values for the model and estimates when using `nnet:multinom`. Yesterday I came across a post where the author put forward a similar issue regarding estimation of $p$-values for coefficients (How to set up and estimate a multinomial logit model in R?). There, one blogger suggested that getting $p$-values from the `summary` result of `multinom` is pretty easy, by first getting the $t$values as follows:               pt(abs(summary1$coefficients / summary1$standard.errors), df=nrow(Data)-10, lower=FALSE)                    (Intercept)   infected+ATB1   infected+ATB2   infected+ATB3    AgeThirty     Moderate 0.002670340   0.08325396      0.014506395     0.02025858       0.0006587898     Severe   0.006433581   0.12665278      0.005216581     0.02352202       0.0035612114      According to Peter Dalgard, "There's at least a factor of 2 missing for a two- tailed $p$-value. It is usually a mistake to use the $t$-distribution for what is really a $z$-statistic; for aggregated data, it can be a very bad mistake." According to Brian Ripley, "it is also a mistake to use Wald tests for `multinom` fits, since they suffer from the same (potentially severe) problems as binomial fits. Use profile-likelihood confidence intervals (for which the package does provide software), or if you must test, likelihood-ratio tests (ditto)." I just need to be able to derive reliable $p$-values.


title: Is this a nested model?

text: I have a dataset with explanatory variable $X$ and response variable $Y$. Besides, there is a two-level factor $S_{j}, j=1,2$. The random error variance might be different between these two strata. I am applying three models to a data set. Model 1, simple linear regression model: $Y_{i}=\alpha+\beta\times X_{i}+\epsilon_{i}$ where $\epsilon_{i} ~N(0,\sigma^{2})$ Model 2, linear regression model with varying variance per sratum: $Y_{i}=\alpha+\beta\times X_{i}+\epsilon_{i}$ where $\epsilon_{i} ~N(0,\sigma^{2}_{j})$. (The coefficients of this model can be estimated using generalized least squares) Model 3: stratified linear model: $Y_{i}=\alpha+\beta_{1}\times X_{i,j=1}+\beta_{2}\times X_{i,j=2}+\epsilon_{i}$ where $\epsilon_{i} ~N(0,\sigma^{2})$ I have read that model 1 and 2 are nested, and the null hypothesis for the likelihood ratio test is H0: $\sigma_{1}=\sigma_{2}$. Similarly, model 1 and 3 are nested, with H0: $\beta_{1}=\beta_{2}$. According to my understanding, a nested model means that we are able to assign some of the parameters to zero. In these examples, I find it quite difficult to assign some parameters to zero. Can someone help me to clarify this? Can I find some formal definition of what a nested model is? I could think of re-wiring $\sigma_{2}=\sigma_{1}+\delta$ or $\beta_{2}=\beta_{1}+\delta$. Thus, both null hypothesis become $\delta=0$. However, my concerns is that the $\sigma$ or $\beta$ in model 1 are estimated from the entire dataset, while $\sigma_{i}$ or $\beta_{i}$ are estimated from strata. Is this a problem?


title: Proper regression model for ratio data

text: I've been thinking about a regression model with dependent variable a ratio in the range [0-600+]. About 50%+ are 0 values though, 40%+ values [0-5], and a very small percentage of the values are higher than 5, so the picture is quite skewed. I've been wondering what kind of transformation I need to do in order to linearize the model. Is the natural log transformation appropriate?


title: How to read files (text and excell) into R?

text: Heloo, pleas help me How to read files (text and excell) into R ? pleas write full instruction. Thank for help


title: Probabilities of events related to drawing cards

text: from a deck of 52 cards, 9 cards are drawn at random without replacement. If X and Y are the number of hearts and diamonds respectively. What is the joint probability of X and Y and also what is the probability of X being greater than Y?


title: Estimate sample size if two properties are measured and their probabilities are known

text: I'm sorry if this is some kind of standard exercise, but I just don't know how to search for it. I have infinite amount of samples and I can measure two properties $A$ and $B$ for each sample. Chance to measure $A$ is $P(A)$ and for $B$ it's $P(B)$, both properties are independent. I take N samples and on $a$ of those I measure property $A$, but property $B$ is encountered $b$ times. If $P(A)$, $P(B)$, $a$ and $b$ are known, how can I get $N$ and what's the error/deviation for this quantity? I can also formulate the question as a problem with some numbers: I have infinite amount of boxes, in each box there is a candy wrapped in paper. 34% of the candies are made of chocolate and 10% of the papers are blue. I've opened some boxes and seen 16 blue papers and have eaten 45 chocolates. How many boxes have I opened? What changes if $P(A)$ and $P(B)$ are not precise but experimentally measured from $N_A$ and $N_B$ measurements instead? That means - before my _measurements_ my brother opened some $N_A=522399$ boxes and that's where the $P(A)=34%$ came from. He's colourblind tho, but my sister saw him opening some $N_B=92587$ boxes and noticed that 10% of the papers were blue. Any changes if there are more than two properties?


title: PCA Using prcomp in R

text: I'm trying to do principal component analysis (PCA) in R using the prcomp function. My input is a large matrix of 1,188 observations (rows) and 15,462 features (cols). I input this to the function with option scale = TRUE. The PCA completed normally, and the loadings ($rotation of the PCA data frame) is a matrix of dim 15,462 x 1,188 with 1,188 PC's. This is bizarre as there are 15,462 features so there should be 15,462 PC's computed. Am I right? Another weird observation is that the projected data points onto the PC's, which I suppose can be extracted via the $x attribute of the PCA data frame is a square matrix of dim 1,188 x 1,188. It should be of the same dimension as the input matrix, shouldn't it? At this point, it looks like prcomp has swapped what should be returned by "rotation" and "x" with each other? I have experimented prcomp using other smaller datasets such as from the example with the USArrests dataset, where the dimensions returned are correct and as expected. Has anybody encountered this problem before? Thanks


title: boot package, cross-validate a linear model, cost argument

text: I would like to cross-validate a linear model using 10-fold CV and calculate the average error rate using the RMSE. I am wondering if the cost argument is the right way to do it. At least for me being not a statistician the help file is not clear. It states "A function of two vector arguments specifying the cost function for the cross-validation. The first argument to cost should correspond to the observed responses and the second argument should correspond to the predicted or fitted responses from the generalized linear model. cost must return a non-negative scalar value. The default is the average squared error function."               library(boot)     data(mammals, package="MASS")     mammals.glm <- glm(log(brain) ~ log(body), family= gaussian, data = mammals)     ( cv.err <- cv.glm(mammals, mammals.glm, K=10, cost = function(y, yhat) sqrt(mean((y-yhat)^2))) )      PS: this related posting deals with a binary model


title: Finding multiple topics in short texts

text: I am building a natural language understanding component for a dialog system. The input is a natural language sentence (usually a short one), and the output should be a set of zero or more classes from a pre-specified set. For example:               input: "I offer a salary of 20000 per month for work as a programmer with a company car"     output: ["OFFER(Salary=20000)", "OFFER(Job=programmer)", "OFFER(Car=yes)"]      I collected some examples from real dialogs, in order to train a classifier. My current approach is:   * Features: I take the words of each sentence as features (I also tried word bi-grams and letter n-grams but the performance was worse).   * Training: I train a separate binary classifier for each class. As positive examples, I take all sentences that have this class, and as negative examples, I take all sentences that don't have this class.   * Running: I run each of the binary classifiers on the new sample, and return the set of classes that correspond to the classifiers that answered 'yes'. This approach lead to too many false positives. Since my dataset contains many sentences with many classes (about 6 classes per sentence), the classifier gets confused between words that belong to different classes. For example, it takes the words "work as a programmer" as positive signal to the class "OFFER(Salary=20000)", because they appeared together in many training samples. So, I edited my training data and separated each sample sentence to several sub-sentences, such that each sub-sentence matches just one class:               input: "I offer a salary of 20000 per month"     output: ["OFFER(Salary=20000)"]     input: "for work as a programmer"     output: ["OFFER(Job=programmer)"]     input: "with a company car"     output: ["OFFER(Car=yes)"]      The new approach lead to too many false negatives: since every sentence is a negative example of all classes except its own class, the classifier takes the words "a salary of 20000" as negative signal to the class "OFFER(Job=programmer)". Therefore, for compound sentences such as "I offer a salary of 20000 per month for work as a programmer with a company car", the classifier finds a single class and misses all other classes. I also created a mixed dataset - half of the one-class samples and half of the multi-class samples (the other half I used for testing), and got much better result. However, this approach seems like a hack, and I would like to know if there is a more professional approach. Any suggestions for solution will be welcome... NOTE: I use the winnow classification algorithm, which currently has the best score (I also tried Naive Bayes, SVM, neural network and perceptron). -- Some numbers:   * Performance of classifier trained on multi-class sentences and tested on single-class sentences: Precision=69% Recall=97% F1=81%   * Performance of classifier trained on single-class sentences and tested on multi-class sentences: Precision=99% Recall=32% F1=48%   * Performance of classifier trained on a mixed dataset (1/2 single-class and 1/2 multi-class), and tested on the other half: Precision=88% Recall=89% F1=89%


title: Density of a scaled log transformed random variable?

text: If $X$ and $Y$ are two independent random numbers from the interval $(0,1)$, then what is the PDF of $U = -2 \ln X$ and $V = -2\ln Y$?


title: Machine Learning using Python

text: I am considering using Python libraries for doing my Machine Learning experiments. Thus far, I had been relying on WEKA but have been pretty dissatisfied on the whole. This is primarily because I have found WEKA to be not so well supported (very few examples, documentation is sparse and community support is less than desirable in my experience), and have found myself in sticky situations with no help forthcoming. Another reason I am contemplating this move is because I am really liking Python (I am new to Python), and don't want to go back to coding in Java. So my question is, what are the more (1) comprehensive (2) **scalable** (100k features, 10k examples) and (3) well supported libraries for doing ML in Python out there? I am particularly interested in doing text classification, and so would like to use a library that has a good collection of classifiers, feature selection methods (Information Gain, Chi-Sqaured etc.), and text pre- processing capabilities (stemming, stopword removal, tf-idf etc.). Based on the past e-mail threads here and elsewhere, I have been looking at PyML, scikits-learn and Orange so far. How have people's experiences been with respect to the above 3 metrics that I mention? Any other suggestions? Thanks in advance. Andy


title: Credit Risk and Concentration

text: I am working with a UK credit-union and we are looking to build a model to assess our credit risk and changes to this over time. We have a number of loans to borrowers who each have a credit rating (lets say these are banded A, AA, AAA, B, BBB, etc...). The way I have proceeded so far is to get the history of defaults in each rating band and build a conjugate model (Gamma/Poisson) to measure the risk profile of each band. I can then use the expected default and variance to build up some measure of expected loss, risk profiles, confidence levels, etc.... The loans are to individual members and as the credit union is not tied to a specific organisation/industry I have treated the individual loans as independent and identically distributed with no strong, identifiable correlation. One factor that worries me is that some members might have abnormally higher loan balances than others. Say, for example I have 300 loans in category B that average £500, but I have 15 loans in category AAA that average £10,000. The exposure in each category is the same, but as the average loan per borrower is higher in one category, should I recognize in my model that a few defaults in this category will have a larger financial effect than a higher number of defaults in the other category? If so, how can I model this? Is a simple Poisson model acceptable for this? I guess I'm asking how I should treat concentration factors when I'm currently assessing the risk as if its homogeneous with regard the borrowers and the amounts the borrow. Any help would be very welcome!


title: How to derive p-value from standard error between two group means?

text: Data was collected on binary variable and continuous variable. The binary variable represents two groups.   * Group 1 has a sample size $n_1$ and mean $\mu_1$ on the continuous variable   * Group 2 has a sample size $n_2$ and mean $\mu_2$ on the continuous variable The significance test is as follows:   * H0: μ1 = μ2 versus Ha: μ1 > μ2, using an $\alpha$ significance level.  The first sample mean was $z$ pooled standard errors above the second sample mean.   * How can the p-value be derived from the standard error $z$?   * What if $n_1=15$, $n_2=30$, $\alpha=.10$ and $z=1.6$, what would the p-value be?


title: Jeffreys prior for geometric distribution?

text: What is the Jeffreys prior for the geometric distribution?


title: Repeated measurement model for Dichotomous variable in spss?

text: I have a repeated data in 7 diffrent time and 4 group. my variable is Dichotomous (yes/no). How I can do it with spss? what non-parametric test could I use? I appreciate your answers Leili


title: Is spectral clustering ever similar to k-means?

text: Does there exist a flavor of spectral clustering (graph construction, type of laplacian, et cetera) that results in a clustering (on arbitrary data sets) that is similar to the k-means clustering "in spirit", i.e. finds a partition of the data into isotropic high density regions? (Of course, I don't expect exact equivalence.) Whenever I've seen spectral clustering, it has seemed to go more after manifolds (i.e. non-isotropic clusters), but that may be because commonly spectral clustering is presented using k-NN type graphs. Of course, a k-NN (or epsilon-neighborhood) graph contains no large-scale information on the relative positions of points, hence it is in no way "encouraged" towards isotropy. Thanks.


title: How to find number of samples needed for a 95% confidence interval

text: I wish to perform a number of tests, what are the minimum number of tests I need to perform to be 95% sure the mean of the test is significant to within ±10% when I repeat the test ? As one example, I wish to sample the mean time to set up a phone call and assume the set-up time is normally distributed. How many phone calls do I need to perform so the ±10% of the mean set-up time falls in a 95% confidence interval ?


title: What's a component in gaussian mixture model?

text: What is the relation between a dimension and a component in a Gaussian Mixture Model? And what are the meanings of dimension and component? Thank you. Please correct me if Im wrong: my understanding is the observed data have many dimensions. Each dimension represents a feature/aspect of the collected data and has its own Gaussian distribution. I don't know where "component" fits into this picture and what it means.