title: Eigenvalues of correlation matrices exhibit exponential decay

text: I have a data-set of $P$ samples of size $N$, and noticed that the eigenvalues of the correlation matrices $A^TA$, when presented in descending order, can in many cases be described as an exponential decaying function. That is, there is a good fit linear of from $i=1..N$ to $\log|\lambda_i|$. Moreover, for several data sets I found that the exponent of the decay is fairly constant. Is this a well-known fact or just one's tendency for finding patterns? Obviously from PCA / SVD it tells me something about the ability to approximate the data using low-dimensional matrix. Are there any solid mathematical results on the size of the exponent of this decay?


title: Bias and variance estimation with boostrap

text: The Wikipedia article about Jacknife estimation of the bias and variance of an estimator $\theta$ includes the following formulas: **Variance** of $\theta$: $ \operatorname {Var}(\theta )=\sigma ^{2}={\frac {n-1}{n}}\sum _{{i=1}}^{n}({\bar {\theta }}_{i}-{\bar {\theta }}_{{\mathrm {Jack}}})^{2}$ where ${\bar {\theta }}_{{Jack}}={\frac {1}{n}}\sum _{{i=1}}^{n}({\bar {\theta }}_{i})$ is the jacknife estimator. **Bias-correction** of $\theta$: $ {\bar {\theta }}_{{\mathrm {BiasCorrected}}}=N{\bar {\theta }}-(N-1){\bar {\theta }}_{{Jack}} $ * * * My question is: What are the corresponding formulas for bootstrap? Are they different?


title: Entropy of state distribution, complexity index and turbulence for sequences of varying length

text: I have a dataset of varying sequence lengths and want to calculate descriptive measured for it. The main cause of differing lengths is that some cases are right censored, and are dealt properly when defining the sequence through the `right="DEL"` argument in TraMineR's `seqdef`. My questions are:   1. Will the entropy of state distribution calculated by `seqstatd` with the option `with.missing=FALSE`, `weighted=FALSE` and `norm=TRUE` be affected by the variation in length?    2. Will the sequence turbulence calculated by `seqST` be affected by the variation in length?   3. Will the complexity index calculated with `seqici()` with the option `with.missing=FALSE` be affected by the varying sequence lengths?


title: What variables explain which PCA components?

text: Using this data:               head(USArrests)     nrow(USArrests)      I can do a PCA as thus:               plot(USArrests)     otherPCA <- princomp(USArrests)      I can get the new components in               otherPCA$scores      and the proportion of variance explained by components with               summary(otherPCA)      But what if I want to know which variables are most explained by what axis. i.e. is PCA1 or PCA2 mostly explained by `murder`, how can do do this? Can I say for instance PCA1 is 80% explained by `murder` or `assault`? I think the loadings help me here, but they show the directionality not the variance explained as i understand it, e.g.               otherPCA$loadings          Loadings:              Comp.1 Comp.2 Comp.3 Comp.4     Murder                         0.995     Assault  -0.995                          UrbanPop        -0.977 -0.201            Rape            -0.201  0.974


title: How to measure the variance from reading to reading

text: I have some values being read by physical sensors. Two of the values from nice lines which are very clear, but the third is all over the place. The problem is calculating a metric to show that the data is almost unusable. I can't use the variance of the data set as the variance is roughly the same as the other readings, it just jumps around a lot more. For a better idea of what I'm trying to say look here: ![enter image description here](http://i.stack.imgur.com/I6XWR.png)


title: Why are RBMs symmetric?

text: I'm doing some experiments with RBMs and note that they use symmetric weights between the input and hidden layer. Why is this? I'm particularly interested in the _intuition_ behind this design decision - why would asymmetric weights not work for example?


title: Generalized log likelihood ratio test for non-nested models

text: I understand that if I have two models A and B and A is nested in B then, given some data, I can fit the parameters of A and B using MLE and apply the generalized log likelihood ratio test. In particular, the distribution of the test should be $\chi^2$ with $n$ degrees of freedom where $n$ is the difference in the number of parameters that $A$ and $B$ have. However, what happens if $A$ and $B$ have the same number of parameters but the models are not nested? That is they are simply different models. Is there any way to apply a likelihood ratio test or can one do something else?


title: Licenses renewals prediction

text: I have come recently to the following real world problem concerning licence renewals of a software product. I have just rudimentary knowledge of the basics in this field and I mostly interested in solving this problem, but I would love to learn something along the way. The input is vector of renewed licenses $(r_i)_{i=1}^n$ indexed by days and a vector of expiring licenses for each day $(e_i)_{i=1}^{n+k}$. That is $r_i$ is number of licenses renewed on day $i$ and $e_i$ is the number of expiring licenses on day $i$. The renewals depend very much at the weekday. E.g. there are lot of renewals on Mondays and Sundays and not very much on Thursdays. The problem is to predict or extrapolate the vector of renewals $(r_i)$ for $i=n+1, \ldots, k$. What I find interesting about this problem is the fact that $r_i$ does not depend only on $e_i$ but also on neighbouring values since people often decide to renew the license before it expires or after it expires. In fact, you know that on average 95% of people renews the license in 60 days centred around the day their license expires and even more - you have a histogram of these renewals $(h_j)_{j=-30}^{30}$. If I understood the data correctly, the value $h_j$ represents the part of renewals that happen on average $j$ days before (or after depending on sign of $j$) the license expires for $j=-30,\ldots, 30$ and the values for $j=-30,30$ cover the renewals that happened more than one month after or before the expiry. In particular $\sum_{j=-31}^{31} h_j = 1$. **Q1:** What is the best approach for this type of problem? After some thought I came up with the following. Let $c_{ij}$ be the proportion of licenses that expire on day $j$ and are renewed on day $i$ to the total number of licenses that expire on day $j$. Then $$ r_i = \sum_{j=1}^n c_{ij}e_j $$ or in matrix form $r = Ce$. Moreover we have $$ h_j \approx \frac{1}{n} \sum_{i=1}^n \frac{c_{i+j,i}}{t_i}, $$ where $t_i = \sum_{j=1}^n c_{ji}$ is the quotient of renewed licenses to all expiring licenses at day $j$. This gives me $n+60$ equations for $n^2$ variables. **Q2:** Is hidden Markov model the right approach here? I know nothing about them so if this approach is the right one I would appreciate some links to materials from which I can learn the theory as well as links to implementations. We can cut down the number of variables by setting $c_{ij} = 0$ for $|i-j| > 30$ but that still leaves us with very underdetermined system. Right now I am about to try out replacing $c_{ij}$ by $h_{i-j}t_j$. This will give me $n$ equations for $n$ unknowns. Which I can hopefully solve and then I will just extrapolate the values $p_i$ for $i=n+1, \ldots, k$ and estimate the renewals by $r_i = \sum_{j=-30}^{30} h_{i-j}t_je_j$. _Update:_ The solutions oscillate towards $i\to n$. Disclaimer: Crossposted at math.stackexchange


title: What are the differences between survival analysis and Poisson regression?

text: I'm working on a classical churn prediction problem using the number of visits of a given user to a site and I thought that Poisson Regression was the right tool for modelling the future engagement of that user. When then I came across a book about survival analysis and Hazard Modelling and I don't know which technique is best. I don't want to be researching both topics at the same time, so what is best for modelling user engagement using past data and demographics?


title: Calculating Average Waiting Time

text: I am trying to calculate average wait time for one of the projects that I am working on. The way the project works is that there are pathways that students follow and each pathway has several steps. At each step students get a challenge, or challenge sets or they earn a badge. Students can earn multiple badges (at different steps) within the same pathway. My task is to calculate the average time the students have to wait in a pathway to earn a badge, and to find out whether it has some effect on the success of the pathway or not. Below is the example that will help understand my question.               PathwayID    TotalSteps    BadgeEarned@Step         1            4                2                                       4         2            7                3                                       5                                       7         3            5                2                                       4         4            9                4                                       6                                       8      They way I am calculating the average wait time right now, is by dividing the average of BadgeEarned@Step by total number of steps in that pathway. So according to these calculations, Average Wait time for the pathways will be below.               PathwayID    AvgWaitTime         1            0.75         2            0.71         3            0.60         4            0.67      Wait Time for Pathway 4 appears to be smaller than wait time for Pathway 1, which is not actually the case, as in pathway 4, students have to wait longer until they even earn their first badge. Right now the way I am calculating the average wait time, there is some bias toward the pathways with more number of steps. I will appreciate if someone can give me an idea to calculate average wait time, that is not biased by the length of the pathway.


title: SAS error: variance covariance matrix singular

text: I am now working with a classmate on an estimate of an Almost Ideal Demand System for tobacco consumption. We have chosen to make an AIDS regression using goods tobacco and alchol, and the CPI price index. Once we modeled the programme on SAS, we got an error that says **ERROR: The iterative procedure stopped because of the singularity of the variance- covariance estimate.** How could I solve this problem? What could have we have ommitted or done wrong? Should we search for more data (thing we want to avoid because the project is urgent?)? Any help will be truly appreciated!


title: Can you use proportions as a covariate in a Cox proportional hazards model?

text: In R's `survival::coxph` function, can I mix a covariate representing proportions (in the range 0.0-0.5) with an integer covariate (in the range 1-15), or should I transform the first one also to integers (0-50)?


title: How to Implement an Empirical Bayes Analysis in BUGS/JAGS/Stan

text: My data is a set of $N$ observations $y_i$. I would like to estimate $\mu$ and $\sigma$ in the following model: $y_i \sim \mathrm{Normal}(\theta, \sigma)$ $\theta \sim \mathrm{Normal}(\mu, \frac{\sigma}{N})$ $\mu \sim \mathrm{Normal}(??, ??)$ $\sigma \sim ???(??,??)$ I would like to perform an Empirical Bayes analysis. Would someone be able to help me and explain how I would go about parametrizing the priors for empirical bayes?


title: When to use robust standard errors in Poisson regression?

text: I am using a Poisson regression model for count data and am wondering whether there are reasons _not_ to use the robust standard error for the parameter estimates? I am particularly concerned as some of my estimates without robust are not significant (e.g., p=0.13) but with robust are significant (p<0.01). In SAS this is available by using the repeated statement in `proc genmod` (e.g., `repeated subject=patid;`). I've been using http://www.ats.ucla.edu/stat/sas/dae/poissonreg.htm as an example which cites a paper by Cameron and Trivedi (2009) in support of using robust standard errors.


title: Difference between sample and sequence in Law of Large Numbers

text: I'm comparing the Chebychev (Weak) Law of Large Numbers (LLN) to the Kolmogorov (Strong) LLN in an econometric textbook and both definitions start off differently. The Chebychev LLN begins with If $x_i, i = 1, . . . , n$ is a sample of of observations such that . . . while the Kolmogorov LLN begins with If $x_i, i = 1, . . . , n$ is a sequence of independently distributed random variable such that . . . I know they refer to convergence in probability and convergence almost surely respectively but what is the significance of using a sample in one and a sequence in the other? Is it harder to prove convergence in a sequence as opposed to a sample?


title: What are some examples of "real world" problems that can be modeled using vector autoregressive models?

text: I would like to stay away from finance and economic examples since they are too abstract for me to understand. Are there any "real world" examples for example healthcare, exam marks, environmental science and so on? From a managerial perspective I would like a basic understanding of how VAR could potentially be used to model real world phenomenon.


title: An analytical framework for considering the geometric mean

text: Is there an analytical method of looking at the geometric mean that will allow one to break it down to its various components? The focus of the question is more for financially related returns, but I am open to consider other fields as well to see if it is applicable. So the aim of the question is to get some sort of results along the lines of               GEOMETRIC MEAN = f(ARITHMETIC MEAN,X ,Y,Z)      where $X,Y,Z$ are all factors/variables that help determine the value of the geometric mean. I know there is a difference in the way its calculated etc, but I guess what I am looking for is what specific factor drives the difference between arithmetic and geometric mean. Most of the time, in financially related contexts arithmetic mean tends to be larger in magnitude than geometric mean due to the variability of return values, i.e. where $x_i$ does not equal $x_j$ where the arithmetic mean is the $\sum \limits_{i=1}^n\frac{x_i}{n}$, and the geometric mean is $\prod \limits_{i}^n {(x_i+1)}^{\frac{1}{n}} - 1$, but just saying that all the returns are not the same is what causes the difference to not quite feel exact enough. Is there something that is better that will show me exactly what causes the gap between geometric and arithmetic means? Approximations using Taylor series would be ok too...


title: Are there any good popular science book about statistics or machine learning?

text: There a bunch of really good popular science books around, that deal with real science, as well as the history and reasons behind current theories, while remaining extremely enjoyable to read. For example, "Chaos" by James Gleick (chaos, fractals, nonlinearity), "A brief history of time" by Stephen Hawking (physics, origin of the universe, time, blackholes), or "The Selfish Gene" by Richard Dawkins (evolution and natural selection). Some of these books present arguments (Dawkins), and some don't (Gleick). But they all serve to make it easy for those of us without in-depth scientific educations to understand otherwise difficult concepts. Are there any such books that focus mainly on Statistics, or machine learning? **Please include a summary of what each book covers.**


title: What is the relationship between substitution cost matrix and indel costs?

text: In TraMineR I might compute dissimilarities as such:               properties <- matrix(c(# left, married, child ,divorced       0, 0, 0, 0,  # parent       1, 0, 0, 0,  # left       0, 1, .5, 0, # marr       1, 1, 0, 0,  # left+marr       0, 0, 1, 0,  # child       1, 0, 1, 0,  # left+child       1, 1, 1, 0,  # left+marr+child       .5, 1, .5, 1 # divorced     ), 8, 4, byrow=TRUE)          sm <- as.matrix(dist(properties))     indel <- .5*max(sm)     dOM <- seqdist(biofam.seq, method="OM", indel=indel, sm=sm)      Why do I provide both indel costs (indel) and substitution cost matrix (sm)? I thought they were one and the same? Please help me to understand the relationship between the substitution cost matrix and the indel costs.


title: Getting positive definite matrix for chol update

text: I am making a square-root UKF implementation. I also use `cholupdate` function in Matlab. However `cholupdate` needs a positive definite matrix. The positive definiteness is tested using `[R,p] = chol(A)` where `p` produces `0`. Then after running `cholupdate`, generates an error saying `A` is not positive definite. My problem is.   1. Does `[R,p] = chol(A)` return `p=0` means `A` is positive-definite always?   2. What is the best method for making a positive definite matrix. Note that here `A` is a square-root matrix where $AA^T$ gives a covariance matrix. Any help is greatly appreciated.


title: What exactly is the name of the type of regression analysis where you try to see if the model is significant over *multiple* start/end values?

text: An example is here: http://www.reddit.com/r/askscience/comments/ine4x/regarding_the_recent_lapse_of_global_warming_in/c2554al I'm sure it's related to robust statistics. But I'm sure that there's a label more specific than that. == Okay, so I tried to run regression analyses on global warming datasets over the last 10 years. What I wanted to show was that there was no significant warming trend over the last 10 years. However, because there are certain anomalous years, I wanted the argument to be more convincing, so I tried _different_ start-year values. My point was to show that the premise behind the article (http://www.physorg.com/news/2011-07-global-linked-sulfur- china.html ) was basically correct. I'm not a climate change denier by any means as I work with climate scientists myself (I hate climate change deniers just as much as any other scientist) - however - the point is that China has released so much sulfur into the atmosphere that climate change had practically been stalled over the last 10 years (that being said, I have no doubt it will resume once again once China reduces the amount of sulfur- containing coal it uses) I'll use the dataset from http://www.ncdc.noaa.gov/cmb-faq/anomalies.php \- and particularly - the ftp://ftp.ncdc.noaa.gov/pub/data/anomalies/annual.land_ocean.90S.90N.df_1901-2000mean.dat dataset (the land+ocean anomalies), which is more scientifically rigorous == Anyways, by performing regression analyses on this data... For the 1997-2010 period, for the coefficient of x, I get a value of 0.00722 degrees per year. This is not statistically significant, as the p-value is .13 1998-2010: 0.00515 degrees per year/p-value of 0.19. Not statistically significant. For the 1999-2010 period, for the coefficient of x, I get a value of 0.01182 degrees per year. This barely meets statistical significance, as the p-value is 0.0457 2000-2010: 0.0089 degrees per year/p-value of 0.17. Not statistically significant 2001-2010: 0.0016 degrees per year/p-value of 0.759. Not statistically significant 2002-2010: -0.00167 degrees per year/p-value of 0.784. But here we have a negative slope 2003-2010: -0.00141 degrees per year/p-value of 0.86. Again, negative slope


title: Kriging on log transformed rainfall data

text: I am beginner in R. I had found in the literature that prior to performing kriging on the data, the distribution has to be investigated to check if it is Gaussian. So, in order to check if the data follows Gaussian I plotted the quantile- quantile plot of rainfall for all days using data from 50 stations. So, inorder to check if the data follows gaussian i plotted the quantile- quantile plot of rainfall for all days using **N=50 stations** data. Shown below is the q-q plot for **good rainy day**. ![enter image description here](http://i.stack.imgur.com/M1flG.png) The q-q plot of **log transformed using transformation function : log(x+1)** ![enter image description here](http://i.stack.imgur.com/aazRE.png) Q-Q plot for **bad(low) rainy day** and its transformed plot are provided at to save space.. http://s20.postimage.org/a3kocwfgd/image.png http://s20.postimage.org/4qvtyrrjx/image.png I have few basic questions:   1. Though the fit of log-transformed data is good on good rainy days,it is not that on days where rainfall is scanty.   2. How to back transform the data after interpolation?   3. How to choose the block size in block kriging? any guidelines. I would like to compare the kriged value against the gridded rainfall from TRMM(25 km by 25 km)?


title: Selecting GAM with/without random effects - residual plots vs. AIC

text: My question is about fitting GAMs with a random effect in mgcv, using s(x, bs="re"). I understand that determining the random effects structure should occur before determining the fixed effects structure in model selection, so I begin with a GAM with all possible smoothers and parametric terms:               SBall.gam2<-gam(count~offset(vol_offset)                 +s(logdepth, by=StageF,bs="cs", fx=FALSE,k=-1)                 +s(logdepth,by=StnF,bs="cs", fx=FALSE,k=-1)                 +StageF*StnF,                 family=Tweedie(p=1.1),                 data=SBall,method="ML")      This model shows heterogeneity in the residuals when they are plotted against StnF. In this study I do not have any main hypotheses about StnF, so I choose to use it as a random effect instead:               SBall.gam2a<-gam(count~offset(vol_offset)                 +s(logdepth, by=StageF,bs="cs", fx=FALSE,k=-1)                 +s(logdepth,by=StnF,bs="cs", fx=FALSE,k=-1)                 +StageF+s(StnF, bs="re"),                 family=Tweedie(p=1.1),                 data=SBall,method="ML")      This model (SBall.gam2a) significantly improves the homogeneity of residuals when plotted against StnF. However, when I compare the two models using AIC (should I even do this? And is ML ok to use vs. REML?), the first model is better:               > AIC(SBall.gam2, SBall.gam2a)                   df      AIC     SBall.gam2  83.08053 2401.317     SBall.gam2a 47.59405 2799.378      In my opinion, it would be better to use the second model which includes the random effect, as it better meets the assumption for homogeneity of residuals, and then proceed with further model selection of smoother and parametric terms. Is this correct? A portion (one station of 9) of the dataset provided below. Many thanks, Raeanne               out <-     structure(list(SBall.StnF = structure(c(9L, 9L, 9L, 9L, 9L, 9L,      9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L,      9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L), .Label = c("1", "2", "3", "4",      "5", "6", "7", "8", "9"), class = "factor"), SBall.StageF = structure(c(1L,      2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L,      3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L), .Label = c("2",      "3", "4", "5", "6", "7"), class = "factor"), SBall.count = c(0,      0, 19, 0, 0, 0, 0, 20, 10, 0, 0, 0, 30, 0, 5, 10, 50, 10, 0,      0, 100, 150, 225, 175, 0, 50, 150, 50, 100, 0), SBall.vol_offset = c(1.59591002173839,      1.59591002173839, 1.59591002173839, 1.59591002173839, 1.59591002173839,      1.5945614038671, 1.5945614038671, 1.5945614038671, 1.5945614038671,      1.5945614038671, 1.61955021036791, 1.61955021036791, 1.61955021036791,      1.61955021036791, 1.61955021036791, 1.61526331179942, 1.61526331179942,      1.61526331179942, 1.61526331179942, 1.61526331179942, 1.64905470098532,      1.64905470098532, 1.64905470098532, 1.64905470098532, 1.64905470098532,      1.61691429450892, 1.61691429450892, 1.61691429450892, 1.61691429450892,      1.61691429450892), SBall.logdepth = c(1.90308998699194, 1.90308998699194,      1.90308998699194, 1.90308998699194, 1.90308998699194, 1.77815125038364,      1.77815125038364, 1.77815125038364, 1.77815125038364, 1.77815125038364,      1.60205999132796, 1.60205999132796, 1.60205999132796, 1.60205999132796,      1.60205999132796, 1.30102999566398, 1.30102999566398, 1.30102999566398,      1.30102999566398, 1.30102999566398, 1, 1, 1, 1, 1, 0.698970004336019,      0.698970004336019, 0.698970004336019, 0.698970004336019, 0.698970004336019     )), .Names = c("SBall.StnF", "SBall.StageF", "SBall.count", "SBall.vol_offset",      "SBall.logdepth"), row.names = c(NA, 30L), class = "data.frame")


title: Use of cosine similarity on DGGE data

text: I am using the software package GelCompar to analyze some DGGE images of processed milk samples from a logitudinal study of 30 sheep. The DGGE technique works by separating bacterial DNA in each milk samples as it migrates through a polyacylamide gel in an electric field. This produces a set of bands in each sample as different bacterial genera migrate to different regions of the gel according to differences in their DNA sequence. This then provides a snapshot of community diversity. Normally, you might expect to use a binary variable i.e. band present/absent to analyse this type of data. However, this becomes very subjective when you have complex banding patterns as I do with my DGGE images. I am therefore basing my analysis on a curve-based method. The standard protocol from the literature appears to be to use a Pearson correlation followed by UPGMA for clustering for this type of data. However, I also have the option of using the Cosine coefficient. I was considering whether the cosine coefficient may be more suitable than Pearson for my data, as it would do some sort of log transformation before analysing it I am guessing from what I have read. Could anyone tell me anymore about cosine and its differences to Pearson as I am trying to see what is the most appropriate for my data at the moment! Thanks


title: Understanding the bivariate Poisson distribution

text: Researching bivariate Poisson over the web is no easy task unless you can make sense of the Greek symbols. I am familiar with Poisson and have a deep understanding of it. So could someone explain bivariate Poisson with its formula so that I can fathom it. My whole model it seems should be based on bivariate Poisson and not Poisson. I have found this paper to be helpful but am still struggling: "Bayesian and Non-Bayesian Analysis of Soccer Data using Bivariate Poisson Regression Models" by D. Karlis and J. Ntzoufras. I understand that the univariate model treats each teams goal scoring as a stand-alone independent process. Thus, the probability of the home team scoring ($\text{HG}$) 1 goal is the same whether the away team scores ($\text{AG}$) 0,1,2 or 10 goals, i.e.: $$p(\text{HG}=\alpha, \text{AG}=x) = p(\text{HG}=\alpha, \text{AG}=y) \quad \text{for all $\alpha$, $x$ and $y$}$$ So the expression above is independent of $x$: $$p(\text{HG}=\alpha,\text{AG}=x)=f(\alpha) \quad \text{for all $x$}$$ With the bivariate distribution, the probability function is a function of both the home and away scores. This allows a dependency between the home and away scores to be included in the model. But I would appreciate it if someone gave a complete example with a working out as it will help my get a thorough understanding of bivariate Poisson.


title: What would a robust Bayesian model for estimating the scale of a roughly normal distribution be?

text: There exists a number of robust estimators of scale. A notable example is the median absolute deviation which relates to the standard deviation as $\sigma = \mathrm{MAD}\cdot1.4826$. In a Bayesian framework there exist a number of ways to robustly estimate the _location_ of a roughly normal distribution (say a Normal contaminated by outliers), for example, one could assume the data is distributed as a t distribution or Laplace distribution. Now my question: **What would a Bayesian model for measuring the _scale_ of a roughly normal distribution in a robust way be, robust in the same sense as the MAD or similar robust estimators?** As is the case with MAD, it would be neat if the Bayesian model could approach the SD of a normal distribution in the case when the distribution of the data actually is normally distributed. _edit 1:_ A typical example of a model that is robust against contamination/outliers when assuming the data $y_i$ is _roughly_ normal is using a t distribution like: $$y_i \sim \mathrm{t}(m, s,\nu)$$ Where $m$ is the mean, $s$ is the scale, and $\nu$ is the degree-of-freedom. With suitable priors on $m, s$ and $\nu$, $m$ will be an estimate of the mean of $y_i$ that will be robust against outliers. However, $s$ will not be a consistent estimate of the SD of $y_i$ as $s$ depends on $\nu$. For example, if $\nu$ would be fixed to 4.0 and the model above would be fitted to a huge number of samples from a $\mathrm{Norm}(\mu=0,\sigma=1)$ distribution then $s$ would be around 0.82. What I'm looking for is a model which is robust, like the t model, but for the SD instead of (or in addition to) the mean. _edit 2:_ Here follows a coded example in R and JAGS of how the t-model mentioned above is more robust with respect to the mean.               # generating some contaminated data     y <- c( rnorm(100, mean=10, sd=10),              rnorm(10, mean=100, sd= 100))          #### A "standard" normal model ####     model_string <- "model{       for(i in 1:length(y)) {         y[i] ~ dnorm(mu, inv_sigma2)       }            mu ~ dnorm(0, 0.00001)       inv_sigma2 ~ dgamma(0.0001, 0.0001)       sigma <- 1 / sqrt(inv_sigma2)     }"          model <- jags.model(textConnection(model_string), list(y = y))     mcmc_samples <- coda.samples(model, "mu", n.iter=10000)     summary(mcmc_samples)          ### The quantiles of the posterior of mu     ##  2.5%   25%   50%   75% 97.5%      ##   9.8  14.3  16.8  19.2  24.1           #### A (more) robust t-model ####     library(rjags)     model_string <- "model{       for(i in 1:length(y)) {         y[i] ~ dt(mu, inv_s2, nu)       }            mu ~ dnorm(0, 0.00001)       inv_s2 ~ dgamma(0.0001,0.0001)       s <- 1 / sqrt(inv_s2)       nu ~ dexp(1/30)      }"          model <- jags.model(textConnection(model_string), list(y = y))     mcmc_samples <- coda.samples(model, "mu", n.iter=1000)     summary(mcmc_samples)          ### The quantiles of the posterior of mu     ## 2.5%   25%   50%   75% 97.5%      ##8.03  9.35  9.99 10.71 12.14


title: Determining if experiment data can be pooled from multiple experiments

text: We have a dataset with Time and Treatment as independent variables, and growth volume as a dependent variable. The experiment was repeated on three occasions, with multiple time points and multiple observations per combination of factors for each experiment instance (so 8 observations at time = 12, treatment = 1, experiment1, etc). There's a subset of common factors for time points and treatment between all the experiments that I am focusing on. Note that the samples were not the same physical samples between experiments (so it's not strictly a repeated measures design?). **We'd like to figure out if we can pool the data from the individual experiments, or if the variation between batches is too significant to merge them.** Originally this seemed like a design for MANOVA, but when I put all of the data in the format I figured I need for R, it looks like there's actually only one dependent variable (growth volume) if we treat experiment as an independent variable. I also read in some other CrossValidated answers that there are regression methods that are preferable to MANOVA nowadays anyway, but I wasn't sure where to start with that.


title: How to deal with this design?

text: I have an experiment with 3 biological factors (`dev`, dev_stage; `reg`, brain_region; `gen`, gen_type) plus a batch factor (`bat`, batch). Unfortunately I have only rather superficial understanding of ANOVA or other statistics, so I'm not sure what is the right formula for my data (I'd consider all 3 biological factors plus their higher interactions with batch effect as an additive factor.). So I hope some luck with the right ANOVA formula here :) As you can see below,   1. `dev` could not be separated from `bat`;    2. At `dev` P28, `gen` and `bat` are unbalanced. So based on many discussions among others seen on the web (mainly here), I guess that bat should be treated as a random effect, and the formula is               aov( x ~ dev * reg * gen, random = ~1|bat)       However, it seems that `random = ~1|bat` was not working since R returned a warning               Warning message:     In lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) :       extra arguments random are just disregarded.       Futhermore, I am not sure how to deal with the unbalance between `gen` and `bat` at `dev` P28... Also, I am not sure if the `gen` and `reg` factors could/should be treated as repeated measures from each batch. Repeated measures are meant to be from the same subject, but to me subject and batch are similar to each other in the sense that they both "shift" the measures in a systematic way, although hidden factors are meant to be controlled within subject and are biasing among batches. But I'm not sure if it's reasonable to compare different levels of `gen` within `bat`. Anyway, I tried this formula:               aov( x ~ dev*reg*gen + Error( bat / (reg*gen)))       Then `dev` became much less significant compared to that of `aov( x ~ dev * reg * gen, random = ~1|bat)` or `aov( x ~ dev * reg * gen + bat)`, and affected less genes than did `reg`. But according to PCA plot, `dev` is the strongest factor, much stronger than `bat`, and `reg` is the second. Also I don't know how to do multiple comparisons in this case, since Tuckey's HSD could not handle the returned values of `aov` when `Error` is used.               ..........brain_region.gen_type.dev_stage....batch     .anteromedial_thalamus........1.....E15.5.11/17/06     .anteromedial_thalamus.......22.....E15.5.11/17/06     ........frontal_cortex........1.....E15.5.11/17/06     ........frontal_cortex.......22.....E15.5.11/17/06     ..somatosensory_cortex........1.....E15.5.11/17/06     ..somatosensory_cortex.......22.....E15.5.11/17/06     ..ventrobasal_thalamus........1.....E15.5.11/17/06     ..ventrobasal_thalamus.......22.....E15.5.11/17/06     .anteromedial_thalamus........1.....E15.5..5/10/07     .anteromedial_thalamus.......22.....E15.5..5/10/07     ........frontal_cortex........1.....E15.5..5/10/07     ........frontal_cortex.......22.....E15.5..5/10/07     ..somatosensory_cortex........1.....E15.5..5/10/07     ..somatosensory_cortex.......22.....E15.5..5/10/07     ..ventrobasal_thalamus........1.....E15.5..5/10/07     ..ventrobasal_thalamus.......22.....E15.5..5/10/07     .anteromedial_thalamus........1.....E15.5...6/1/07     .anteromedial_thalamus.......22.....E15.5...6/1/07     ........frontal_cortex........1.....E15.5...6/1/07     ........frontal_cortex.......22.....E15.5...6/1/07     ..somatosensory_cortex........1.....E15.5...6/1/07     ..somatosensory_cortex.......22.....E15.5...6/1/07     ..ventrobasal_thalamus........1.....E15.5...6/1/07     ..ventrobasal_thalamus.......22.....E15.5...6/1/07     .anteromedial_thalamus........1.......P28.10/19/07     .anteromedial_thalamus.......22.......P28.10/19/07     .anteromedial_thalamus.......22.......P28.10/19/07     ........frontal_cortex........1.......P28.10/19/07     ........frontal_cortex.......22.......P28.10/19/07     ........frontal_cortex.......22.......P28.10/19/07     ..somatosensory_cortex........1.......P28.10/19/07     ..somatosensory_cortex.......22.......P28.10/19/07     ..somatosensory_cortex.......22.......P28.10/19/07     ..ventrobasal_thalamus........1.......P28.10/19/07     ..ventrobasal_thalamus.......22.......P28.10/19/07     ..ventrobasal_thalamus.......22.......P28.10/19/07     .anteromedial_thalamus........1.......P28.11/30/07     .anteromedial_thalamus........1.......P28.11/30/07     .anteromedial_thalamus.......22.......P28.11/30/07     ........frontal_cortex........1.......P28.11/30/07     ........frontal_cortex........1.......P28.11/30/07     ........frontal_cortex.......22.......P28.11/30/07     ..somatosensory_cortex........1.......P28.11/30/07     ..somatosensory_cortex........1.......P28.11/30/07     ..somatosensory_cortex.......22.......P28.11/30/07     ..ventrobasal_thalamus........1.......P28.11/30/07     ..ventrobasal_thalamus........1.......P28.11/30/07     ..ventrobasal_thalamus.......22.......P28.11/30/07.      I'm looking forward to your help!


title: When should I use $E(X)$ and when to use $\mu$  for the mean?

text: I see that $E(X)$ and $\mu_x$ are used to refer to the mean of some distribution or variate, but I am not clear on what circumstances each is appropriate. $E(X)$ seems to always be more often in the context of a random variable $\mu_x$ seems to appear when referring to the mean of some distribution. for example my books lists the formulas for all the standard distributions similar to for example Binomial mean : $\mu_x=np$ Some times the distinction seems to be whether its the subject of interest, or just being used as a valued variable. eg; $$E(T_z) = \mu_T(1+\beta t)$$ in that case both $E(T_z)$ and $\mu_T$ both represent the expectation of the variable. I guess in that sense one is a function that takes an argument, and the other is a variable to be assigned to, and read from.


title: How to make representative sample in dataset with two variables?

text: I have two variables and 1000 cases. How can I statistically find representative cases from total of 1000, based on statistical properties of both variables and correlation between them. Perhaps something based on T-test and 95% (or 99%) interval but for both variables? I would like to know which statistical method can find cases that have both values (simultaneously) statistically the most significant. I know that this deals with sample distribution and estimating the proportions.


title: How to calculate the hat matrix for logistic regression in R?

text: I want to calculate the hat matrix directly in R for a logit model. According to Long (1997) the hat matrix for logit models is defined as: $$H = VX(X'VX)^{-1} X'V$$ X is the vector of independent variables, and V is a diagonal matrix with $\sqrt{\pi(1-\pi)}$ on the diagonal. I use the `optim` function to maximize the likelihood and derive the hessian. So I guess my question is: how do i calculuate $V$ in R? _Note:_ My likelihood function looks like this:               loglik <-  function(theta,x,y){     y <- y     x <- as.matrix(x)     beta <- theta[1:ncol(x)]     loglik <- sum(-y*log(1 + exp(-(x%*%beta))) - (1-y)*log(1 + exp(x%*%beta)))     return(-loglik)     }      And i feed this to the optim function as follows:               logit <- optim(c(1,1),loglik, y = y, x = x, hessian = T)      Where x is a matrix of independent variables, and y is a vector with the dependent variable. _Note:_ I know that there are canned procedures for doing this, but I need to do it from scratch


title: How to detect patterns between fields of a distribution in SPSS?

text: Suppose I have the following simplified distribution:               time | value     1    | 2     2    | 4     3    | 8     4    | 16          1    | 1     2    | 3     3    | 9     4    | 27          1    | 40     2    | 20     3    | 10     4    | 5          1    | 12     2    | 1     3    | 99     4    | 23423      These are all part of the same dataset (so one x can have multiple values here, eg. time = 1 corresponds with value = 2,1,40,12). I separated them because the first 3 have an obvious pattern within their slope (2, 3 and 0.5) and the last does not have a pattern in its slope. Time is in days, and the values represent a quantity. Now, how do I use SPSS to find the elements of a distribution that together form a pattern (any kind of pattern) in their slope? And how can I make it also include elements that almost follow this pattern, but have slight variations (a variation that I can set)? Any feedback is appreciated. Thanks. **Update:** It looks like SPSS is not the best tool for the task. I am interested in the R language however. If anyone could recommend any books regarding this field of pattern recognition in R, that would be great.


title: How to build an incremental classification strategy?

text: I want to build a tool that allows users to label data instances by themselves. So the user selects some instances and then train a classifier. Afterwards the user might choose to select more instances. However I don't want the classifier to be build allover again because it's a time consuming process. Is there a good strategy for building such an incremental classification tool? One idea I have in mind is to build many classifiers and then to combine them into one. So in the first iteration the user selects some instances and the classifier is built. Then in another iteration he selects another instances and _a new_ classifier is built. Now suppose a new instance is needed to be classified, I let the instance go through all the classifiers and take the majority vote. I guess this is called bagging. However how valid is this technique given that in each iteration different numbers of instances are selected? Rather than the exact number of instances as in bagging.


title: How to make a Neural network understand that multiple inputs are related (to the same entity)?

text: I am not sure if this is the right place to ask this but here goes: Sometimes times two or more inputs of a neural networks can often be related to a single "real world" entity. E.g : `Height` and `weight` of a person to predict the probability of disease in population or `price` and `volume` of a stock to predict the market. When a single training set contains data about a number of these entities, how can we make a neural network understand that two inputs (or more) are related to the same entity? Amongst all the people I have asked, the general consensus seems to be:   * Neural Networks do not work this way   * It is not possible   * Such a grouping of data is not required   * Neural Networks are supposed to find the relationship amongst inputs, you are not supposed to feed it the relationships   * The training data set need to be reworked / reconfigured   * I have never heard of such a thing So, obviously this is not in the mainstream. Has anyone heard of any research in this direction? **P.S**. If you agree with the above opinion (it can't / shouldn't be done) please provide a reason why.


title: Odds ratio traps 1 when p value significant

text: I have a set of imputations on which I'm using PROC GLIMMIX to build a regression model for a binary outcome. I'm using PROC MIANALYZE to roll up the output, part of which includes upper and lower confidence limits which can be exponentiated to obtain a confidence interval for the odds ratio. GLIMMIX also yields Type III F tests of the model effects, but it's generally conceded that chi-square tests are more appropriate for these models. A good approximation can be obtained by multiplying the F value by the numerator degrees of freedom. That's tenable in my sample seeing that it's very large. I found it difficult to roll up the F tests in MIANALYZE, so I used a data step to compute the chi-square and average over the 10 imputations, and then simply tested as a chi-square at the numerator degrees of freedom. Here's the rub: A couple of my model effects are significant when I roll up the chi-square, but the odds ratios trap one. I don't know how to resolve the paradox. Is there some variance due to imputation not accounted for in the roll up of the chi- square tests? How do I get around this problem?


title: What is the magnitude of bias in censored regression when OLS is applied?

text: If my dataset comprises few censored variables (<1%) and I fit the OLS regression using a heteroscedastic resistant estimator (the residuals are not terribly heteroscedastic to begin with)- are the results valid?


title: Adjusting for Confounding Variables

text: I have data from human participants in a study. There are more females in the study (60%) and males are older. I have a binary categorical variable $O$. If those who are $True$ for $O$ are older, do I need to correct for sex and/or age? Maybe those $True$ for $O$ contain more men. What concepts or methods should I use to determine this? I'm using R.


title: Why are there more bayesians than frequentists in astrophysics?

text: I heard somewhere that in astrophysics, there are more bayesians than frequentists (although in general, I think that bayesian community is smaller). Is there something special about "astro" data what makes it more suitable for bayesian methods? **EDIT** : The motivation for my question was also this kaggle competition, where the bayesian approach outperformed others. Therefore it seems that there could be something special about astro data.


title: What if residuals are normally distributed, but y is not?

text: I've got a weird question. Assume that you have a small sample where the dependent variable that you're going to analyze with a simple linear model is highly left skewed. Thus you assume that $u$ is not normally distributed, because this would result in normally distributed $y$. But when you compute the QQ-Normal plot there is evidence, that the residuals are normally distributed. Thus anyone can assume that the error term is normally distributed, although $y$ is not. So what does it mean, when the error term seems to be normally distributed, but $y$ does not?


title: Non-normal observations in regression modelling

text: I read an article that says the dependent variables in a regression model must be normally distributed. The way i understand it, is that the observations for the regression model must then be normally distributed. Or in other words if i choose sample data from a population then the sample must be normally distributed. But then i read answers to a similar question in another forum and apparently the the observations or dependent variables does not need to be normally distributed, just the error terms after the model is implemented must be normally distributed. So now i am not sure, does the observations (dependent variables) has to be normally distributed as well as the error terms or is it just the error terms.


title: Non-normal data and non-parametric tests

text: I have two non-normal variables (one DV, one IV) and a few 7-point Likert scale IVs (normally distributed). The non-normal variables are centrality scores from network analysis - DV is from the collaboration matrix; IV is from the friendship matrix. Transforming the data did not work out - unfortunately. Here is what I would like to test: DV: centrality score - collaboration IV: centrality score - friendship annual budget number of employees years of operation I think I have to go with a non-parametric test. The problem - at least for me - is that the 2 non-normal variables are continuous and have different ranges. I thought about assigning a rank instead of the centrality scores, but I am not sure which test to use. I have the feeling I am overlooking something quite simple.


title: Regression on a non-normal dependent variable

text: I need to do a regression with a non-normal DV for which no proper non-linear transformation (that I know of) exists: ![enter image description here](http://i.stack.imgur.com/7fPjE.png) It is a score ranging from 10 to 50, with a high peak at 10, a drop at 11 and a regular decline from 11 to 50. The distribution of residuals is not normal. ![enter image description here](http://i.stack.imgur.com/pSjZ3.png) There are more or less 500 subjects in the study, and the covariates are all dichotomous. I thought of using n-tiles and performing an ordinal regression on the resulting variable, but then I end up with a high proportion of empty cells -- not because of a low sample size, but rather because of the several covariates that need to be in the model. The best I could do for now is to remove covariates and use tertiles (instead of, say, quintiles) to minimize the number of empty cells. I am not sure, though, to what extent avoiding empty cells is THAT important. Compared to the original analysis on a raw score, regressing on tertiles with less covariates seems like a lot of sacrifices. Based on the decent sample size, and given the fact that the distribution of residuals is not THAT far from normal, I am wondering if the results of the regression would be reliable as they are.


> **Possible Duplicate:**   >  Regression on a non-normal dependent variable If the dependent variable in a multiple regression is not normally distributed, what can be done to correct it?


title: Multiple comparison in the unequal variance case

text: I have found procedures such as Tamhanes T2, Dunnets T3 and the Games & Howell procedure that deal with unequal variances in the one-way model. However, I have a Randomized complete block design, which is basically a two-way model. And variances differ quite a lot in the treatment group, from 2 to 6. In SAS for mixed models 2nd ed. (p. 369) Ramon C. Litell et. al., uses a unequal variance model, proc mixed data=TV ic; class age sex; model time=sex|age/DDFM=KR OUTP=R; repeated / group=age; lsmeans age sex / diff adjust=Tukey; run; However, I am not sure this is correct since the multiple comparison test (Tukey), uses a pooled estimate for the variance and which would give deflated p-values when the variances are unequal. Is this the correct way of performing multiple comparison under unequal variance?


title: Let F be a distribution function. Prove that X is a RV

text: Let $F$ be a distribution function on $(\Omega,\, \mathfrak{F},\, P)=\left((0,1),\ \mathfrak{B}(0,1),\ \lambda\right)$ where $\lambda$ denotes Lebesgue measure. Define $X:\Omega \to \mathbb{R}$ by $X(\omega) = \sup(y \in \mathbb{R}: F(y) < \omega)$. 1 Show that $\forall x \in \mathbb{R}, (\omega: X(\omega) \leq x) = (\omega: \omega \leq F(x))$ 2 Show that $X$ is a RV in $(\Omega, \mathfrak{F}, P)$ and that $F_X = F$. My attempt: LHS = $(\omega: X(\omega) \leq x)$ = $(\omega: X(\omega) \in (-\infty, x])$ = $(\omega: \sup(y \in \mathbb{R}: F(y) < \omega) \in (\infty, x])$ RHS = $(\omega: \omega \leq F(x))$ = $(\omega: \omega \in (-\infty, F(x)])$ = $(\omega: \omega \in (-\infty, P(X^{-1}( \ (-\infty,x] \ ))))$ = $(\omega: \omega \in (-\infty, P(LHS)))$ I'm stuck. Help please? :( Cross posted : http://math.stackexchange.com/questions/888054/let-f-be-a- distribution-function-prove-that-x-is-a-rv


title: What to do when some time points have heavily skewed responses and some do not in a repeated measures study?

text: Typically, when one encounters continuous but skewed outcome measures in a longitudinal design (say, with one between-subjects effect) the common approach is to transform the outcome to normality. If the situation is extreme, such as with truncated observations, one might get fancy and use a Tobit growth curve model, or some such. But I am at a loss when I see outcomes that are normally distributed at certain time points and then heavily skewed at others; transformation may plug one leak but spring another. What might you suggest in such a case? Are there "non-parametric" versions of mixed effects models that I am not aware of? Note: an applied example would be knowledge test scores pre/post a series of educational interventions. Scores begin normal but then cluster at the high end of the scale later on.


title: Predicting time series with OpenBUGS

text: I have a number of fairly short time series (about 4–100 observations) which I need to forecast into the future. I decided to use Bayesian inference, because there is external information about each time series (the time series are mutually similar), and prediction probabilities would be very desired. I chose to create a Bayesian model resembling the classical AR model: $$ y_{i+1} = b_0 + b_1 y_i $$ The prior for $y_i$'s in the future should be a normal distribution derived from all time series. However, I have a problem combining these two sources of information (from the autoregressive model and prior) together in an OpenBUGS model (via _BRugs_ ). This is a fairly simple model implementing only the autoregressive model:               model {         b0 ~ dunif(-20,20)         b1 ~ dunif(-20,20)         sd ~ dunif(0, 10)         tau <- 1/(sd*sd)              for (i in 2:N) {             mu[i] <- b0 + b1*y[i-1]             y[i] ~ dnorm(mu[i], tau)             z[i] <- y[i]         }              for (i in (N+1):(N+K)) {             mu[i] <- b0 + b1*z[i-1]             z[i] ~ dnorm(mu[i], tau)         }     }      where `N` is the number of observations, `K` number of predicted values, and `z[(N+K):(N+K)]` are the predicted values. It works pretty well with simple time series such as (1,2,4,8,16), but fails with real time series, as it tends to grow exponentially (it is not stable). I belive that restricting y's to a normal prior distribution would alleviate the problem, but I do not know how to extend to model to do that. It should be something it the sense of:               z[i] ~ dnorm(mu0, tau0)      where `mu0` and `tau0` are defined externally, but that would not be possible here, as `z[i]` is already defined. How to construct such a model?


title: Cramer-Rao lower bound questions

text: I've been reviewing questions from a statistics exam of the last year. There is a question with the probability density function below $$\displaystyle f(x,\theta) = \frac 1{2\theta^3}x^2e^{-\frac x\theta}$$ where $\displaystyle0<x<\infty$ and $\displaystyle0<\theta<\infty$. Question is as follows: a) Find the maximum likelihood estimator for $\displaystyle\theta$. b) Find the minimum variance estimator $\widehat{\theta}$ for the given pdf. c) Is the maximum likelihood estimator obtained in a) efficient and consistent? I have found the answer to a) as $\widehat\theta=\frac {1} {3n}\sum x_n$. However, I couldn't find the Cramer-Rao lower bound to the end. I've shown the equation below $$\displaystyle \frac {d^2\log(f(x,\theta))} {d\theta^2} = \frac 3 {\theta^2}-2\frac x {\theta^3}$$ I should take the expectation of this, but I'm not sure what should I replace $E[x]$ with. Should I just make it equal to $\displaystyle \theta$ and go further with the solution? By the way, could you also hint me how can I find the variance of the maximum likelihood estimator above for c)? I'd appreciate it Thanks in advance


title: Best value to store ratio data and compare it to time period average

text: I'm looking for a simple way to store ratios. For a time component, I must store the average ratio between two behavior. For example the number of people that turn left compared to the number of people that turn right. I have to detect unusual behavior (people that turn right abnormally). How should I mathematically compare the average ratio against the analyzed ratio, and how should I display the difference on a graph ? Thanks a lot in advance.


title: Computing by hand the optimal threshold value for a biomarker using the Youden Index

text: I have an empirical estimate of a ROC curve, that is, a plot of the sensitivity versus 1-specificity over all possible cut-off values of the marker. Based on an empirical ROC curve, I would like to determine the optimal cut-off point that represents a better trade-off between sensitivity and specificity. I have read that the Youden Index can be used in that purpose. Here is an example:               oneMinusSpecificity <- c(1.00000000,                               0.636363636,                              0.436363636,                              0.315151515,                              0.163636364,                              0.096969697,                              0.072727273,                              0.006060606,                              0.000000000)     sensitivity <- c(1.00000000,                      0.91566265,                      0.77108434,                      0.66265060,                      0.39759036,                      0.33734940,                      0.28915663,                      0.07228916,                      0.00000000)      which results in the following ROC curve. The vertical line represents the Youden index = largest distance between "1 - specificity" and "sensitivity". ![enter image description here](http://i.stack.imgur.com/FpnOk.png) I find               > youden <- max(sensitivity - oneMinusSpecificity)     > youden     [1] 0.3474991      **How can I calculate the optimal cut-off using this information?**


title: Good games for learning statistical thinking?

text: Are there any games that get the player "think like a statistician"? For example, lightbot gets you to "think like a programmer" (in a very basic way). Are there any games - designed for entertainment or teaching - that can help get one comfortable with basic concepts like correlation, p-values, least squares, variance, different kinds of probability distributions, regression to the mean ... One example would be this correlation-guessing game. (I'm asking because I'm thinking about developing such an application, and am trying to get a broad view of what prior work exists)


title: Where can I find the articles discussing about the suggested training-data partitioning ratio?

text: Right now, I'm doing my thesis which use multilayer perceptron learning method to train a model. What I learned from my class is the purpose to partition the data which can be separated into 3 group   * Training dataset - This set is to train the model.   * Validation dataset - This set is to find the best model parameters.   * Testing dataset - This set is to evaluate the performance of the model. But the problem is I cannot find the reference to support why I have to partition the samples in my thesis with certain ratio. I try to google it but can't find those articles. It's probably because I googled it with the incorrect technical terms. Back to the question, Can anyone here suggest the articles discussing the best ratio with supported evidence?


title: Distance between the product of marginal distributions and the joint distribution

text: Given a joint distribution $P(A,B,C)$, we can compute various marginal distributions. Now suppose: \begin{align} P1(A,B,C) &= P(A) P(B) P(C) \\\ P2(A,B,C) &= P(A,B) P(C) \\\ P3(A,B,C) &= P(A,B,C) \end{align} Is it true that $d(P1,P3) \geq d(P2,P3)$ where d is the total variation distance? In other words, is it provable that $P(A,B) P(C)$ is a better approximation of $P(A,B,C)$ than $P(A) P(B) P(C)$ in terms of the total variation distance? Intuitively I think it's true but could not find out a proof.


title: Non-conjugate prior

text: Can someone explain why the integral in the posterior density may not be "analytically tractable" if the prior we choose is non-conjugate?


title: How to bootstrap for the standard error of a two sample test

text: I want to calculate the standard error of the area under the curve computed with the Mann-Whitney U-Test (MWUT). Hanley and McNeil 1982 give a formula for it (they call the AUC calculated with the MWUT "W"). However I'm finding this hard to implement so I'd like to bootstrap it to check/replace my work. My question is how do I sample from the two distributions to perform the bootstrap? Shall I resample from them independently or resample from the whole set together and keep track of the case labels? Another concern I have is how the resampling will affect the MWUT, given that there will be tied ranks.


title: Why does converting data into a Time Series change the values and how do you change them back?

text: I imported data and converted into a simple time series.                 claims <- read.csv("C:/claimsrawdatabyweek.csv")             # time series with 52 data points per year         em<-ts(claims[,4],start=c(2009,44),frequency=52)      It changed the values from               3496     2453     2311     2559     3149     5103     4951     4165     4340     4868     3862     2811     3583     3469     3171     2660     2772     2223     2115      to               138     105     97     109     130     155     154     146     148     153     142     119     139     137     132     112     116     93     82      How do I calculate what these new values indicate?   I am using the tbats forecast function `fitem<-tbats(em)` and it is giving me forecast values that are base on these numbers. Even when I don't convert the data and use tbats, it automatically converts to these numbers. Stated simply, I must convert the new ts values into the actuals. Why does 3496 equal 138? Any help is appreciated.


title: ARIMAX model's exogenous components?

text: Does anyone know, considering an ARIMAX model that fitting a stationary process Y, then do the exogenous components for the model need to be (weakly) stationary? I think exogenous components can be any process, even non-deterministic ones, am I right?


title: How to find whether the k nearest neighbors of an instance belongs to same class or not

text: I am working on imbalanced datasets which have two classes: majority and minority class. Here I want to find whether the minority class examples are majority class. How can I do that? I have implemented a code which finds k nearest neighbours for the minority class. This is my code:               A=[1 2 1;3 4 1;5 6 1;];      B=[1 2 2;3 14 2;15 16 2;17 18 2;1 2 1;3 4 1;5 6 1;];     %A=A1;B=B1;     [row,col]=size(A);     [row1,col1]=size(B);          dist=zeros(row,row1);     neararray = zeros(row,row1);     neararray1 = zeros(row,k);     k=5;          for i=1       for j=1:row1             dists(i,j)=sqrt(sum((A(1,1)-B(j,:)).^2));%for each element       end            [y,index]=sort(dist(1,:));       neararray(i,:)=index';       for k=1:5       neararray1(i,k)=neararray(i,k);       end     end                     Input:              A is:         feature1,feature2,class              1     2     1              3     4     1              5     6     1         B is:         feature1,feature2,class              1     2     2              3    14     2             15    16     2             17    18     2              1     2     1              3     4     1              5     6     1      My doubt is this approach correct? For example `A(1,1)` I will get nearest neighbours as: `1,5,6,7,2 But How can I say that it `A(1,1)` belongs to which class?


title: Correlation vs. independence: statistical power, available tests

text: Correlation is frequently applied, but in fact often the interest of the experimentalist is not so much to find a significant correlation as a significant _dependence_ of the variables, and clearly independence is not the same as lack of correlation. To give a particular example, consider two genes and their expression across a number of individuals. If genes are co-regulated (and therefore possibly functionally related), their expression values will not be independent. Whether they will be correlated, however, is another matter. Now, there are tests out there that can be seen as an alternative to a test for a significant correlation -- test of variable independence for continuous variables, such as the Longest Increasing Subsequence Independence Test (R package LIStest) or a non-parametric test described here. I have very little background and no experience with these tests. My question is: are these tests comparable in power to correlation tests? Do you have any practical experience? What would you recommend?


title: How to model change over four time points on an ordinal variable where some participants have missing data at some time points?

text: ### Context I am accumulating a data set that is starting to look like this:   1. I have about fifty individual animals (n=50) that I am studying over, say, a year   2. Every three months or so, I would measure each one and assign a score about their health. So at the end of the year I would theoretically have made four measurements/samples.   3. The scores are 0,1,2,3,4, with 0 being the best and 4 being the worst. I.e., I think I have ranked, discrete data?   4. However, everytime I only get to measure about 30 out of the 50 individuals. So at the end of the year, some animals may have been measured everytime, while others may have only been measured one or two times. ### Question Based on the above, what would be a good statistical test to see if the entire group of 50 individuals have changed over the four measurements? I am guessing some sort of repeated measurement, ranked non-parametric test might be suitable?


title: Choosing a measure of similarity to quantify similarity between individuals on a set of personality scales

text: I have a bunch of users. Each user has a number of personality attributes, such as "fitness level" or "eco-consciousness", rated on a scale from 1 to 5. I want to calculate how similar two users are, so I can show each user a sorted list of "most similar users". This seems to be a classic IR problem, and I've seen three different metrics used, but no discussion of why to choose one over the other:   * **Simple Arithmetic.** The scores are already normalized to the same scale, so I can just add each user's scores up, and compare the sums to see who is most similar.   * **Cosine Similarity.** Treat each user as an n-dimensional vector, where each scale is one dimension. Calculate the cosine of the angle between two users' vectors; cosines closer to 1 (smaller angles) are more similar.   * **Euclidean Distance.** Each user is an n-dimensional vector again, but this time, calculate the distance between endpoints. Users that are close together are similar. What are the advantages and disadvantages of each method? How does that change if the scores are not normalized to the same scale (i.e. if I add an "age" attribute)?


title: Is it feasible to predict a numerical variable with only 100 measures in 70 000 observations?

text: A quick explainer: I have a data set of 5 numeric variables, and a numeric target variable. The total observations / rows are about 70 000, but I only have about 90 measurements of the target variable. Cluster analysis of the 70 000 returns from 5 to 8 clusters. In reality I expect there to be about that many groups or classes as well, and that the target variable would naturally group within them. I wish to predict values of the target variable for those observations where it is missing, and have some There are no missing values in the 5 variables for the 70 000 data set. A slightly larger total of observations (100 000) would exist, with a greater number of variables (nearly 30), if I allowed missing values. There would be an increased number of target variable measurements as well, to 114. I have been attempting to use self-organising maps to model and then predict values, hoping to take advantage of its simplicity, plus I was introduced to them via a proprietary software package! Is it feasible to predict values for the target variable from this data set?


title: How do I reconcile how EViews estimates an AR(1) model with OLS?

text: I discuss estimation of the following $ARX(1)$ model: $Y_t=α+βX_t+u_t$ where $u_t=ρu_{t-1}+ϵ_t$ Substituting the the value of $u_t$ in the first equation into the second, we have that: $Y_t-α- βX_t =ρ(Y_{t-1}-α-βX_{t-1})+ϵ_t$ Rearranging, we have: $Y_t=α(1-ρ)+ρY_{t-1}+ βX_t-ρβX_{t-1}+ϵ_t$ Can we obtain estimates for the parameters in equation 2 using OLS? Well, yes. However, notice that if we did this, we’d be ignoring the fact that the parameters in the model are restricted, in particular the parameter of $X_{t-1}$, $- ρβ$, is the product of the coefficients on $Y_{t-1}$ and $X_{t}$. We should utilise this information. The OLS estimator will be unbiased, consistent, and BLUE (best linear unbiased estimator), but by utilising the information regarding the restriction of parameters, we can find a better estimator than OLS. In other words, the non-linear estimator produced by the Marquardt algorithm will be superior to OLS. Unsurprisingly, EViews estimates all ARMAX models using the Marquardt algorithm. Consider the following 3 situations: $1.$ If $ρ=0$, then our model becomes: $Y_t=α+ βX_t+ϵ_t$ The best thing to do in this case is to ignore $Y_{t-1}$ and $X_{t-1}$, since we know, in theory, they are irrelevant in explaining $Y_{t}$. The Marquardt algorithm is not necessary. $2.$ If $β=0$ and $ρ ≠1$, then our model becom$Y_t=α(1-ρ)+ρY_{t-1}+ϵ_t$ The best thing to do in this case is to ignore $X_t$ and $X_{t-1}$ since we know they are irrelevant variables. $3.$ If, in the population, $X_{t-1}$ is uncorrelated with the other regressors, in particular $X_{t}$, then we can just regress $Y_t$ on a constant, $Y_{t-1}$, $X_{t}$. Again, this is extremely unlikely. Especially so because we are dealing with time series variables that are in all but the most rare circumstances autocorrelated. If, in the sample, $X_{t-1}$ is uncorrelated with the other regressors, a very unlikely scenario, then the OLS estimates of the $α(1-ρ)$, $ρ$ and $β$ when we ignore $X_{t-1}$ are exactly the same as from the Marquardt algorithm, by the Frisch-Waugh-Lovell theorem, as long as the estimate of $ρ$ is not precisely zero (which is, in any case, unlikely). We can obtain the estimate of $α$ from the Marquardt algorithm by simply multiplying the OLS estimate of the constant by 1 minus the OLS estimate of $ρ$. This is not so much a question. I would like your expertise regarding whether you think my assertions are correct. It is very difficult to find answers to this online. It is also quite an interesting problem. Your help is greatly appreciated. Thank you. Christian


title: Agreement of clustered data

text: I have the following situation: I have analyzed several data curves from a group of patients (16 curves per patient) with different analysis methods and want to test for the agreement of the methods. So far, I have neglected the potential correlation within the patients and was thus able to compute ICC (of agreement) values, which yielded very reasonable results. Unfortunately, I have reason to believe that the data are correlated within the patients. I am now looking into mixed linear models and generalized estimating equations, both of which deal with this situation of clustered data. My question is: Is there any way to calculate something similar to the ICC (or CCC) which tests agreement or what would you guys use as a measurement of concordance between the methods? I guess I can get beta values, but something stronger would be nice, that you can actually attach some statistical significance to.


title: Does measurement error include effects of moderator variables?

text: I am not very sure how to interpret measurement error correctly: as a constant, as a bias, or as a moderator factor? Is it presumed that measurement error includes moderator effects? Let it be in the context of meta-analysis of effect sizes.


title: How to test periodical patterns of events in a timeline

text: _**Imagine the situation:_** there is a lecture of a famous scientist in auditorium which begins in - say 9:00. If students do not arrive in time - the door to auditorium will be closed and nobody allowed in. There are many styles of arrival ( _random, wave-like, cummulative..._ ). _With my skills I am only able to discriminate between cummulative and random pattern by measuring the time gaps between subsequent arrivals. In cummulative pattern the gaps will be smaller as the time of the lecture approach. In random pattern there is no correlation between the length of time gaps and remaining time._ **How do I test various (more complicated) patterns? (f.e. how do I prove that pattern is wave-like and not random statistically). The approach of measuring time between subsequent arrivals will not work, because it expands and contracts multiple times.** **Are there any statistical tests which can work with such data?** ![Each line is the arrival of one student](http://i.stack.imgur.com/3C53k.jpg)               # Example data     # wave-like pattern     wave.pat<-c(0,5,7,8.5,9.5,10.5,11,11.25,11.5,11.625,11.75,11.875,12,12.25,     12.5,13,14,15,16.5,18.5,23.5)     wave.like.data<-c(wave.pat, c(wave.pat+28.5), c(c(wave.pat+28.5)+28.5),      c(c(c(wave.pat+28.5)+28.5))+28.5)     # cumulative pattern     cum.pat <- numeric(84)     cum.pat[1] <- 5     for (i in 2:7)     cum.pat[i] <- cum.pat[i - 1]+2     for (i in 8:17)     cum.pat[i] <- cum.pat[i - 1]+1.5     for (i in 18:30)     cum.pat[i] <- cum.pat[i - 1]+1     for (i in 31:55)     cum.pat[i] <- cum.pat[i - 1]+0.5     for (i in 55:84)     cum.pat[i] <- cum.pat[i - 1]+0.25     cum.pat[83]<-64.25     cum.pat[84] <- 64.50     cum.pat.data<-c(cum.pat + 44.5)     # random pattern     my.ch<-seq(0,109, 0.125)     random.pat.data<-sample(my.ch, 84, replace=FALSE)          # SOURCE DATA FRAME     my.data<-data.frame(timing=c(wave.like.data,cum.pat.data,random.pat.data),                         type=rep(c("wave.like", "cumm.like", "random"), each=84))          # plot the data     library(lattice)     stripchart(timing ~ type, data=my.data, pch="|", ylim=c(0,4))


title: Did I correctly specify my model with nested and interacting random factors in glmer?

text: I have the following data: I planted seed bags in the vicinity of 10 mother plants into three directions and four distances (10, 31, 56, or 100 cm). The germination success (0/1) was assessed by extracting one seed bag per microsite for the next three years. I am interested in the effect of distance on germination success and possible differences in germination between years. The data have a lot of 0 values. I originally used the following syntax:               glmer(germination~distance+(1|plant/direction/year),family=binomial,data=seed)      but I was told that distance, year and direction can interact and not always in the way as specified by this type of nesting. Would this be the correct one?               glmer(germination~distance+(1|plant/direction)+(1|year),family=binomial,data=prot)      Or would it be useful to also test interactions between fixed and random effects? From a biological point of view, I do not expect them. I am afraid of making the model more complex than necessary.


title: Is there a bias correction for effect size in a data mining context?

text: Given $K$ possible 'treatments' of some kind, and independent observations of some response under those treatments, say $X_{i,k}$ for $i=1,\ldots,n_k$ and $k=1,\ldots,K$, I am faced with the classical data-mining dilemma. The task is to simultaneously:   1. Find the treatment, say $k^*$, which maximizes effect, and   2. Estimate the effect size of treatment $k^*$. Here you can think of effect size as a measure of location, say mean. The naive approach to this problem is as follows:   1. Estimate the effect size of each treatment. For example, if the effect size is the mean, let $\bar{X_k} = (\sum_i X_{i,k}) / n_k$.   2. Pick the treatment which maximizes the estimated effect, _i.e._ let $k^* = \arg\max_k \bar{X_k}$.   3. Estimate the effect size of $k^*$ by $\bar{X_{k^*}}$. The problem, of course, is that the estimate of the optimal effect size is positively biased because the choice of optimal treatment is not independent of the estimated effect sizes, which are unconditionally unbiased. Additionally, the goals are often somewhat broader: instead of picking the treatment which maximizes effect, we are tasked with sorting the treatments by estimated effect size and estimating the effect sizes. Because of this, I am looking for some kind of function $f$ which 'unbiases' the estimated effect sizes, but preserves monotonicity. That is, $f$ takes a $K$-vector to a $K$-vector, and   1. $f$ preserves monotonicity: if $\bar{X_{i_1}} \le \bar{X_{i_2}} \le \ldots \le \bar{X_{i_K}}$ and $(Y_1,Y_2,\ldots,Y_K) = f(\bar{X_1},\bar{X_2},\ldots,\bar{X_K})$, then $Y_{i_1} \le Y_{i_2} \le \ldots \le Y_{i_K}$.   2. $f$ (approximately, possibly under weird assumptions) unbiases the effect size estimates: if $(Y_1,Y_2,\ldots,Y_K) = f(\bar{X_1},\bar{X_2},\ldots,\bar{X_K})$, and $Y_{i_1} \le Y_{i_2} \le \ldots \le Y_{i_K}$, then $Y_{i_j}$ is a (nearly) unbiased estimate of the effect size of the $j$th smallest effect treatment. Probably such a function would have to be aware of some additional information, like the number of observations $n_k$ and some measure of the spread of the response under each treatment. Probably the monotonicity requirement is questionable if the standard errors of the effect estimates vary wildly. Sweep that under the rug if possible. I am also interested in approaches to this problem when the response observations are not independent, rather they are paired..


title: Self study help with question on Hypothesis Test for 2 related samples

text: Good day, I am currently attempting the following question on Hypothesis Testing for 2 related samples. My eventual answer turned out to be **t=2.5584** , while the rejection region is **2.2060**. However, when I consulted the answer key, I seem to be wrong, the answer provided was |t| = 4.536; Reject H0 if |t| > 2.3060. Appreciate some advice on where I might be going wrong. **_Question_** > Nine experts rated two brands of Colombian coffee in a taste-testing > experiment. A rating on a 7-point scale (1 = extremely unpleasing, 7 = > extremely pleasing) is given for each of four characteristics: taste, aroma, > richness, and acidity. The following data display the total ratings – > accumulated over all four characteristics. > > ![enter image description here](http://i.stack.imgur.com/OBFpF.jpg) > > A test at the 0.05 level significance is conducted to test whether there is > a difference in the mean of total ratings between the two brands. Assume > that the populations of both brands are normally distributed. Find the test > statistic for this test and the rejection region of this test. **_My Attempt_** > ![enter image description here](http://i.stack.imgur.com/3ZEl3.jpg)


title: R-square from rpart model

text: How can I extract the R-square from a fit rpart model? `rsq.rpart(fit)` plots the two graphs, but I simply want to extract the R-square value for the full tree. I assume this is fairly obvious, but numerous searches didn't really lend anything useful. Any help you can provide will be greatly appreciated.


title: Alternatives to using Coefficient of Variation to summarize a set of parameter distributions?

text: ## Background I have a model with 17 parameters, and I currently use the coefficient of variation ($\text{CV}=\sigma/\mu$) to summarize the prior and posterior distributions of each parameter. All of the parameters are > 0. I would also like to summarize these pdfs on a normalized scale (in this case standard deviation normalized by the mean) so that they can be compared to each other, and with other statistics presented in similar adjacent plots (sensitivity, explained variance). I will include density plots for each parameter separately, but I would like to summarize them here. However, the sensitivity of the CV to $\mu$ causes the following confusion that, although easily explained in text, would be preferable to avoid.   1. the posterior CV of one parameter is greater than the prior because the mean has decreased more than the variance (parameter `O` in figure).   2. one of the parameters (`N`) is in units of temperature. It has a 95% prior CI of (8,12 Celsius $\simeq$ 281-285K); when I present the data in units of Kelvin which is only defined for positive values, the CV is <1%, if presented as C, the CV is closer to 40%. To me, it seems that neither of these CVs provides an intuitive representation of the CI. ## Question Are there better ways to present this information, either as a CV or as another statistic? ## Figure As an example, this is the type of plot that I am planning to present, with posterior CV in black and prior CV in grey. For scale, the CV of parameter `O` is 1.6. ![alt text](http://i.stack.imgur.com/CVeO1.png)


title: Checking if improvement is significant. Test design

text: I have designed two simple tests to check my algorithm. For a set of 40 pictures, 20 people will judge if:   Test 1: Object A or B is seen in the picture (i know the correct answer)   Test 2: Object A, B or C is seen in the picture (i know the correct answer) Afterwards i will change the pictures with my algorithm and rerun the test with the same group of people. To check if there is an improvement i want to compare the number of correct answers (from test 1 + 2) before and after the use of the algorithm. Can i mix both test results and use McNemar's test to check if there is a significant improvement (that improvement is not just accidentally)? I am not sure if this is the right way. Many thanks.


title: logistic regression model for process with memory

text: Say there is a stream of binary variables x(t) where t is a discrete index. Say the model is that p=p(x(t)=1) is logistic regression on past realization: log(p/1-p) = x(t-1)b1 + x(t-2)b2 + ... x(t-L)bL + b0. Now to fit the model, one can prepare the design matrix X which contains consecutive lags of order L so that X is (N-L)L in size: X = [x(1) ... x(L); x(2) ... x(L+1); ... x(N-L-1) ... x(N-1)] and the dependent variable Y is just [x(L+1), x(L+2), ... x(N)] in order to maximize likelihood, standard function such as glmfit in MATLAB assumes that the variables X and Y are independent, however the model is such that they are not. How can this be solved? Thanks.


title: Predicting likelihood of attending next event

text: I work in a sports stadium and I would like to run a predictive model on all people who have previously bought tickets to events, to predict the likelihood of them returning to the next event, or indeed any event in the future. We have hundreds and thousands of customers, and approximately 40 historical events in which we can use as binary data on whether or not they attended each event. Does anyone know of a method I can use to run such a predictive model? I am not really sure where to start, although I believe this can be done. Thanks very much in advance.


title: What are the data needed to conduct Cost-Effectiveness Analysis?

text: I'm aiming to conduct a cost-effectiveness analysis (CEA) of the introduction of a new vaccine for HPV (Human Papillomavirus), in a developing country ( I think this will be done from the societal perspective ) I have a dataset from a national prevalence survey (sociodemographic, habits, eaten foods, a variable indicating if the women has HPV or not... ) I have determined the risk factors and the prevalence of the disease given this dataset, but the dataset has nothing related with the vaccine. From what I have read, I think that I must have effectiveness of the vaccine (QALY), the total cost... I'm really new to this field of analysis. My questions are about the first steps of the CEA: What are all the data I need? How can I estimate them given my dataset? If I will not use my dataset, from where can I get needed information?


title: Self study question on Confidence Interval

text: I am attempting a self study question on confidence interval. My answer worked out to be $20.4 \pm 1.411$. However, the model answer appeared to be $20.4 \pm 1.486$. I am pretty confidence that I am correct. I'd appreciate if anyone could give me some advice and validation please? **_Question_** > The owner of a big egg farm wants to estimate the mean number of eggs laid > per chicken. A sample of 25 chickens shows they laid an average of 20.4 eggs > per month with a standard deviation of 3.6 eggs per month. Assume the number > of eggs laid per chicken per month has a normal distribution. Construct a > 95% confidence interval estimate for the population mean number of eggs laid > per chicken per month. **_Attempt_** ![enter image description here](http://i.stack.imgur.com/ESJ5x.jpg)


title: Selecting a feature modeling approach for text classification

text: I am new to text processing. Currently I am trying to determine which type of feature vector I need for a classification problem. I am mainly deciding between binary feature modeling and statistics-based approaches, such as term frequency/inverse document frequency (tf-idf) or chi square. In terms of classification, I have 100 document related to computer science, 100 related to biology, and 100 that are associated with other disciplines. I would like to build a system that, given a new document, can determine which of those groups it belongs in.


title: Forecasting asset returns using index models in R

text: How do you forecast returns and the associated risk in R using index models? How do you represent risk in multi index models as a single value in R?


title: R - QQPlot: how to see whether data are normally distributed

text: I have plotted this after I did a Shapiro-Wilk normality test. The test showed that it is likely that the population is normally distributed. However, how to see this "behaviour" on this plot? ![enter image description here](http://i.stack.imgur.com/NpI0O.png) **UPDATE** A simple histogram of the data: ![enter image description here](http://i.stack.imgur.com/3nOTu.png) **UPDATE** The Shapiro-Wilk test says: ![enter image description here](http://i.stack.imgur.com/39ABw.png)


title: Matrix factorization and gradient descent for recommender systems; user bias?

text: I've been reading about using Matrix Factorization techniques to do collaborative filtering. A popular thing to do seems to be to add user and item biases into the ratings prediction. What I don't understand is why, in some papers (like Factorization meets the Neighbourhood, Koren, Y., 2008), these parameters are learnt. Shouldn't the biases be static whilst training the model? They show tendencies in the actual ratings, right? Thanks!


title: MMSE estimation with random system matrix

text: Suppose our system is given by $y=Ax+n$ where $y$ is observed vector, $x$ is signal vector, $n$ is white noise vector and $A$ is a matrix of i.i.d zero mean random variables. Usually in linear minimum mean squared error estimation $A$ is assumed to be known. I want to estimate assuming $A$ is unknown except that $A$ is independent of $x$ and all elements of $A$ have zero mean. Then LMMSE is given by $$\hat x=Cov[yx]Cov^{-1}[yy]\cdot y$$ Therefore we need to calculate $Cov[yx]=E[y^Tx]$ (assuming zero mean of noise). My question is this $E[yx]=E[x^TxA^T]=E[x^Tx]E[A^T]=0$. So how can I proceed? Thank you.


title: Estimate of censored poisson process

text: I have a set of processes, each of which has number of events and the total length of time. I'm trying to model them as independent Poisson processes with there own rates. The rate of the ith process may be estimate thus $$ \lambda_{i} = \frac{n_{i}}{T} $$ however many of the processes do not generate any events and so have a $\lambda$ of zero We can use the estimates of lambda to generate simulated processes, with the goal of comparing the empirical recruitment curve (when the first event was emitted for each $i$) to simulated ones under the model. From here. Draw $X_{i} \sim Pois(\lambda_{i}T ) $ and use to generate $X_{i}$ uniform random numbers on $(0\space T)$. Take the smallest value to be the recruitment for each $i$. This seems to give a biased estimate, since _some_ processes with a non-zero observed rate are simulated with a zero rate, but _no_ processes with a zero- observed rate are simulated with a non-zero rate. It seems like this is a similar question to this one. Although, unlike lifetime data, in this case some of the processes may truly have a zero rate. EDIT I'm adding a simulation in matlab that seems to suggest that simulating the recruitment as $n_{i}$ uniformly distributed events rather than $Pois(n_{i})$ uniformly distributed events is the correct thing to do               N = 500;     L = nbinrnd(0.5,0.2,N,1);    %real distribution of lambdas (could be any distribution though)     nEvt = poissrnd(L);          %actual number of events        rEvt = [];                   %actual recruitment           for i = 1:numel(nEvt)         if nEvt(i)>0             rEvt(i) = min(rand(nEvt(i),1));   %simulate the true recruitment         else             rEvt(i) = inf;         end     end          nS = 2000;                         %shuffles     xsim = nan(N+1,nS);                %save the xpostion of the ecdf     xsim2 = nan(N+1,nS);               %save the xpostion of the ecdf     for s = 1:nS         nEvtSim = poissrnd(nEvt);         %number of events simulated as a poisson draw         nEvtSim2 = (nEvt);                %number of events simulated SAME as number of actual events         rEvtSim = nan(numel(nEvtSim),1);           rEvtSim2 = nan(numel(nEvtSim),1);                for i = 1:numel(nEvtSim)             if nEvtSim(i)>0                 rEvtSim(i) = min(rand(nEvtSim(i),1));             else                 rEvtSim(i) = inf;       %never had an event             end                  if nEvtSim2(i)>0                 rEvtSim2(i) = min(rand(nEvtSim2(i),1));             else                 rEvtSim2(i) = inf;       %never had an event             end          end              [f,xsim(:,s)]= ecdf(rEvtSim);   %work out the ecdf         [f,xsim2(:,s)]= ecdf(rEvtSim2);   %work out the ecdf     end     xsim = sort(xsim,2);     h1 = stairs(xsim(:,nS*[0.025 0.975]),[f f],'color',[0.9 0.9 0.9]); %95% confidence interval     hold on          xsim2 = sort(xsim2,2);     h2 = stairs(xsim2(:,nS*[0.025 0.975]),[f f],'color',[1 0.9 0.9]); %95% confidence interval          [f,x]= ecdf(rEvt);     h3 = stairs(x,f,'color','k');  %real data recruitment     hold off          legend([h3 h1(1) h2(2)],{'actual data','poiss(N) events','N events'})


title: Fisher Information

text: Suppose $X$ is a one dimensional Normal variable. Find the $2 \times 2$ fisher information matrix $I(\theta)$ for $\theta = (E(X), E(X^2))$. Basically, what I did was to write down the pdf in terms of $E(X)$ and $E(X^2)$...and so I got Let $d = E(X^2)$ (so $\sigma^2 = d - \mu^2$) $(\frac{1}{2(d-\mu^2)})^{1/2} e^{[-(X-\mu)^2]/[2(d-\mu^2)]}$. And then I tried to compute the matrix from that in the usual way....but it is literally taking me forever to compute this matrix (especially when it comes to computing the first and second derivative of $\log f(x,\mu, d)$, with respect to $\mu$). So I was just wondering if there was an easier approach for this problem. I do know that the fisher matrix for the normal distribution with $\theta = (\mu ,\sigma)$ is $$\begin{bmatrix} 1/\sigma^2 & 0 \\\ 0 & 2/\sigma^2 \end{bmatrix}$$ So I was wondering if there was any way to manipulate this to get the fisher matrix when the parameter is $\theta = (E(x), E(x^2))$. Thanks in advance


title: Difference between pooled t-stat and normal t-stat

text: Suppose I have data like this -               Val Bin     -1  Y     5   N     -2  N     4  Y      so forth - where Bin is a binary value. I want to construct a 95% confidence interval for the average difference between the values that have a Y associated with them and the values that have a N associated with them. Note - I have more values associated with Y than with N(which is why I'm not clear about using a normal T-stat). I'm unclear as to whether I should use a pooled t-stat or a normal t-stat. What should I be using?


title: Modeling noise in a feed-forward neural network

text: I am very new to neural networks and have looked at some to the online and youtube tutorials. I would like to model an XOR output for two inputs, one of which is noise. In other words, my output should be an XOR of noise and data. I have a couple of questions: i) since this is a simple well-known model, there must be a standard neural network configuration that is used. Wikipedia shows a diagram of a 2-3-1 feed- forward network (2 inputs, 3 hidden nodes in the hidden layer, one output): http://en.wikipedia.org/wiki/Feedforward_neural_network. Are there other configurations and which is best? ii) is noise usually modeled as one of the inputs? iii) the data is in fact a set of N observations, x1 ... xn, x is either 0 or 1. The input is sum(xi). Should the xi be grouped and modeled as a single input or individually. In other words, should this be treated as a 2 input problem or a N +1 (noise) input problem. The activation function is logistic(b0 + b1x1 + b2x2 ...). Any help to get me started would be appreciated. I plan to model this using the R nnet library.


title: How to correct for gender in Mann-Whitney U test?

text: My statistical knowledge is very limited. I have two groups (patients and controls) and a measured concentration value for each indivdual. Not normally distributed: Mann-Whitney   How do I correct for gender in Mann-Whitney? (I want to know if the difference between patients and controls is due to difference in gender distribution, since males have higher values than females).


title: Interpretation of regression coefficients obtained from applying left inverse of regressor matrix in an underdetermined system?

text: If $X^\dagger$ is the pseudo-inverse of $X$, $\beta = X^\dagger y$ is the least squares solution for $\beta$ when $y=X\beta$. In the overdetermined case, applying $X^{\dagger,L} = (X^TX)^{-1}X^T$ corresponds to the minimum norm residual (superscript $L$ used to indicate left inverse). In the underdetermined case, applying $X^{\dagger,R} = X^T(XX^T)^{-1}$ corresponds to the minimum norm $\beta$ (superscript $R$ used to indicate right inverse). If I solve for $\beta$ in an underdetermined system using $X^{\dagger,L}$ rather than $X^{\dagger,R}$, what does the solution correspond to, conceptually?


title: How to code separate data for a regression parameter in nonlinear regression?

text: Say I have a bunch of measurements over a range of my independent variable (example data at the end) and I want to fit them to an equation such as a classic pharmacological inhibition curve: $Y=B+\frac{M-B}{1+\exp{\left( \ln{X} - \ln{IC_{50}}\right)}}$ I cannot natively use data with $X=0$ in determining the parameters with NLR as $\ln{0}$ is not defined, but theoretically as $X\rightarrow 0$, $Y\rightarrow M$. Experimentally, it's often easier to combine a negative control and a "baseline" into one treatment by having an untreated group. How can I use this untreated group in determining $M$ in a single NLR analysis? I specifically use SPSS, but any advice with a more theoretical explanation is fine. Note 1: I know I can just substitute an arbitrarily small value for $X=0$, such as $X=10^{-20}$, but this seems a little... cheaty. Also, I feel like my definition of "arbitrarily small" could affect the regression outcome, especially if weighted regression is involved. Is there a less arbitrary way? Note 2: If it makes it easier, this doesn't need to be about NLR but instead can be about any asymptotic independent link function in GLM. EDIT: Maybe I chose a bad example. Here's another one: $B+\frac{M-B}{1+\exp{\left(H\left( \ln{X} - \ln{IC_{50}} \right) \right)}}$ which is not readily simplified. The equation doesn't matter- what I want to know is this: Can I specify some data to be used in determining one parameter while other data is used to find other parameters, all while including the uncertainty contained in all measurements in my final model fit? Edit 2: so far everyone is fixated on the example equations, which really aren't important for my question. For the pedants out there, consider this equation that I am making up: $A+\frac{B}{1+C\left( \ln{\left( X \right)} - D \right)}$ * * * Here's an example dataset:               X   Y     0   0.9732462     0   0.9857081     0   1.00479     0.001   0.9811251     0.001   0.982346     0.001   1.028089     0.003   0.9731305     0.003   0.9784695     0.003   0.9815493     0.009   0.9254836     0.009   0.9342329     0.009   0.9448794     0.027   0.8100319     0.027   0.7948713     0.027   0.7979656     0.081   0.5945458     0.081   0.5938009     0.081   0.6131012     0.243   0.3620556     0.243   0.3847672     0.243   0.3633385     0.729   0.2013999     0.729   0.1919975     0.729   0.2240313     2.187   0.1351639     2.187   0.1432508     2.187   0.1508913


title: Testing for equality of multivariate distributions

text: I have two distributions $F$ and $G$ that I conjecture to be identical, mathematically. Essentially, before investing time in mathematically proving that $F$ and $G$ are equal, I'd like to sanity check it by doing a simulation experiment. So, given that I can easily simulate $X_i \sim F$ and $Y_i \sim G$, I'd like some test that checks $F = G$. The $X_i$ and $Y_i$ live in $\mathbb R^p$, and I'm in the process of checking that the marginals are equal in distribution by simulating many data-sets and checking that the Kolmogov-Smirnov test statistics are approximately uniformly distributed, but I'd like some assurance that I'm not being deceived by just looking at the marginals. In addition to the marginals, I'm checking the distribution of the principal components of the generated data/some nonlinear transformations of the data, and so far everything points to equality. If someone could point me to a method - hopefully with an existing implementation in R, or otherwise something easy to implement - that would be fantastic. **EDIT** : To be clear, I'm asking for a multivariate test. I know the marginals aren't sufficient, I'm only checking them because if the marginals don't match then I don't need to check the joint.


title: SPSS-independent-sample-t-test

text: I want to test mean differences between 4(0,1,2,3) groups for 6 different cell-types. In other words, I would have 6 comparisons to do for each cell- type (0-1,0-2,0-3,1-2,1-3,2-3). I have run a repeated-measures ANOVA before but it only gives me the main effects and interactions. To correct my t-tests for multiple comparisons, I would use Bonferroni, does SPSS let me do that, I just cannot seem to find it except as part of the means-testing using ANOVA? Or do I have to do it manually and how would i do that? alpha-level here is 0.05.Thanks


title: Fast multivariate unimodal density estimator

text: I have a sample $\boldsymbol{x}_i$ for $i$ in $1,\dots, n$, from a $d$ dimensional density $f(\boldsymbol{x})$ and I would like to estimate this unknown density. In addition I know that $f(\boldsymbol{x})$ is unimodal, but it can be skewed or fat tailed. Given that $d$ is around 20, I think that using a standard kernel density estimator (KDE) is not an option (I can use a sample size $n$ of around $10^4$). In addition given that I know that $f(\boldsymbol{x})$ is unimodal, I don't think that all the flexibility of KDEs is needed in my case. Hence I was looking for a parametric density estimator, that can be fitted reasonably fast. I've read a bit about multivariate skewed normal or student-t distributions but I would like to know if there are other options out there (maybe constrained KDEs or semi-parametric density estimators?). Thanks!