title: How is interpolation related to the concept of regression?

text: Explain briefly What is meant by interpolation.How is it related to the concept of regression? interpolation is art of reading between the lines of a table and in elementary mathematics the term usually denotes the process of computing the intermediate values of a function from a set of given or tabular values of that function. I can't give the answer of second question. Please help


title: Using the EM Algorithm for unimodal distributions?

text: I've really only seen EM used for mixtures where one can point out multiple modes visually - e.g, the classic mixture of gaussians example. I would like to use EM for a mixture of an empirically defined, sharply peaked distribution and something that is more uniform - does anyone have an idea as to how much confidence I should put in the resulting estimates, or prior experience with similar applications?


title: Saturation in ARIMA (et al) models?

text: I've been learning X12-ARIMA by looking at data from a friend's service company, and wondering how to model the capacity of the company. That is, if the company is limited by a particular resource to only be able to handle 1,000 customers a week, how do I keep my ARIMA model from happily predicting 1,200 customers next summer? (This isn't an issue with time series like GDP or stock prices, which don't have a hard cap.) It doesn't seem that you can do anything in the optimization phase (which is simply choosing parameters), nor with exogenous variables (which drive the process, not react to it). Maybe changing the ARIMA model to a State Space representation would help? (Any recommendations on an R package to do this? I've looked at several and DLM's many matrices confuse me at this point.)


What do "marginal" and "conditional" mean in "marginal models" and "conditional models"? Are they related to marginal distributions and conditional distributions? Thanks!


title: componentwise boosting based on fisher scoring

text: componentwise boosting dates back at least to Bühlman and Yu (2003), where in each boosting iteration a set of base-learners (e.g. simple linear models) depending on a subset of the covariates are used. This implies some variable selection. Let us consider the case, in which we use base-learners with only one covariate. If we have $P$ covariates, we then have $P$ base learners in each boosting iteration. At the end of each boosting iteration, we select the best base-learner among this $P$ base-learners for the update. Since only one covariate is in cluded in each base-learner, we update only the parameters for one covariate. In the literature one can find componentwise boosting based on functional gradiend descent on the one hand and componentwise boosting based on one-step Fisher scoring (e.g. Tutz and Binder 2006; Tutz and Binder 2007 ...) on the other hand. My question refers to the latter, in which a penalized maximum likelihood estimation (ridge penalty) is used to update the parameters by one- step Fisher scoring. Here comes my question: If we have only _metric covariates_ and use componentwise boosting with _only one single predictor in each step_ , do we then need **penalized** maximum likelihood estimation? My guess is no, because in my view penalization is only needed if we use componentwise boosting, where we want to update a _set of predictors_ in each step (and not only one single predictor). Another reason where I see the need for penalization is the case, where we have a single multicategorial covariate (with $K$ categories) that has to be updated. Each multicategorial covariate has $K-1$ associated parameters and thus updating this multicategorial covariate means updating all $K-1$ associated parameters. Here, one can do penalization to update only a subset of the $K-1$ associated parameters, namely the 'important' ones, while the 'unimportant' ones are shrinked towards zero due to penalization. Since we want to update **one single metric predictor** in each boosting step, we don't need penalization, right?


title: Frequency weighted least square estimation

text: I recently realised an RLS Algorithm which fit a stepresponse of a system with an underdamped complex pole pair. Now i realised that there is a bias for some kind of responses i get. I found out that using an estimation based on weighting the frequency range where the input spectrum has the most power, gives me a much more accure and stable parameter set. (tryed on Matlab Toolbox arx command , with focus: simulation) Now the question, is there a descrpition how to implement this method? Which things of the LS / RLS algorithm do i have to alter and how?Are there Papers, Sources available from Mathworks? Bye and thank you!


title: How to simulate data for those distributions which are not defined in MATLAB

text: Suppose there is any distribution whose parameter values are known. So how to simulate data in this case?


title: Multiple Comparisons Correction on Peason's Chi-Squared Tests

text: I have a dataset consisting of 1000 samples (rows). Each sample consists of 300 categorical data points (columns). The process of generating these data points has a known probability of generating an erroneous result, $pErr$. Let $pErr = 0.1$ for this example. As each of the 300 columns consists of 1000 categorical data points, the number of data points, which can be considered as errors in each column, is $1000 * pErr = 100$. Consider this value $Expected$. I need to know the lowest number of errors per column, which is considered statistically significant. I run $\chi^2$ tests, starting with $Observed = Expected + 1$, and increasing $Observed$ by $1$, until the test statistic is greater than the $\chi^2$ value, for one degree of freedom. In this example, an $Observed$ value of 119 is statistically significant: test statistic = 4.01 > $\chi^2$ = 3.84. I then examine each of the 300 columns for categorical data points occurring less than 119 times, and consider them to be errors. **Question:** Do I need to apply a multiple comparisons correction to the analysis? If so, where? **Edit:** Each data point (each column for each row) is independent. I am not interested in any overall significance or testing. I am only interested in each column as an entity. The number of errors for each column are considered separately from any other column.


title: Lower tail of power law distribution

text: I have noticed that in some datasets that can be modeled using a power law distribution, there appears to be a kind of second tail trailing off at the lower end which appears as an almost horizontal line in log-log (for instance, see the figure below; source). ![enter image description here](http://i.stack.imgur.com/eOo0z.jpg) Is there any theory as to why this might happen?


title: Using linear regression for bias correction

text: I have a predictor and the ground truth. I have found that I can use linear regression for bias correction, so instead of using the predictor's estimate directly, we can use the one obtained from the LR between the ground truth and the predictor. I am actually finding it confusing what bias correction actually means. Any pointers?


title: Mix of I(0) and I(1) in VECM

text: Is there an official academic journal I could use as reference which states that a mix of I(0) and I(1) variables can be used in constructing a VECM?


title: Strange pattern in residual plot from mixed effect model

text: I've run a mixed effect model in R by using lme. The explanatory (Temp_Diff & Distance) and responsive (LF_Diff) factors are continued variables.               VARIABLE = MEAN(STD)          Temp_Diff = 0.31(0.8)          Distance = 141(115)          LF_Diff = -3.62(9.02)      Normally, residuals should be plotted against fitted values to identify violation of homogeneity. However, I got a residual plot with an obvious pattern. I can't show the picture because I am new, but the pattern shows a linear line increaed from minus to positive values. that I couldn't figure out what the causes are. I tried to drop some outlier points, but the pattern is still there. Does anyone have any experiences about that?               # ---The Code---------------------------------------     ashlf <- read.table("c://ashlf.txt", header=T);     attach(ashlf);names(ashlf);     ashlf.lme <- lme(LF_Diff~Temp_Diff+Distance, random=~1|Site/Year, data=ashlf, method="ML");     summary(ashlf.lme);     stepAIC(ashlf.lme);     plot(ashlf.lme);          ##---The output of summary(ashlf)-------------------          AIC     BIC    logLik       275.8956 285.877 -131.9478          Random effects:      Formula: ~1 | Site             (Intercept)     StdDev:    1.098288           Formula: ~1 | Year %in% Site             (Intercept) Residual     StdDev:    6.905381 1.401097          Fixed effects: LF_Diff ~ Temp_Diff + Distance                      Value Std.Error DF   t-value p-value     (Intercept)  3.607701 2.0825324 18  1.732363  0.1003     Temp_Diff   -4.273440 1.6115431 18 -2.651769  0.0162     Distance    -0.041130 0.0111402 18 -3.692059  0.0017      ![random=~1|Site/Year](http://i.stack.imgur.com/xpIRU.jpg) ![random=~1+Site|Year](http://i.stack.imgur.com/3Pa34.jpg)


title: How to write this cost equation?

text: I have a mean time to failure variable that is Weibull-distributed with shape parameter less than $1$. Whenever a failure occurs, corrective action of cost $c_c$ is performed. From time to time, preventive maintenance of cost $c_m$ is performed, which resets the time in the Weibull distribution. My question is, how do I write an annual cost equation using this? This is what I have so far. My expected number of failures in a year is $\frac{365}{\mu}$, where $\mu$ is the mean of the Weibull. Therefore, the cost for corrective action is: $$c_c * \frac{365}{\mu}$$ When I add in the preventive maintenance, it becomes something like: $$c_c * \frac{365}{\mu} + c_m * f_m$$ Where $f_m$ is the frequency of doing preventive maintenance in a year. But the problem is the frequency of preventive maintenance affects the number of times corrective action is done (in this case it increases it). So a more correct formula would be: $$c_c * fcn(f_m) + c_m * f_m$$ Where $fcn(f_m) = \frac{365}{\mu} + something$, so for every increase in $f_m$, the function value increases. Obviously, $something$ should have $f_m$ in there somewhere, but I don't know where. The idea is that every increase in frequency of preventive maintenance increases the expected number of times corrective maintenance has to be done, since the function is being reset to time zero and there is a large infant mortality in this distribution. I tried thinking about it in terms of just probability, but it's difficult. Say preventive maintenance is done some time $t_p$ after corrective maintenance. Then the probability it will have broken down in that time is $Pr(t<T=t_p)$ and the probability for the Weibull time to failure resets to time zero. But then I'm stuck, because how can I model the corrective actions in probability terms? Plus, I don't know how to convert that probability into frequency. I dont think it's as simple as $f = \frac{365}{t_p}$ because preventive maintenance happens on fixed intervals with respect to the previous preventive maintenance, not with respect to the previous corrective action, so $t_p$ always varies. And since the occurrence of one preventive maintenance resets the time to failure function, then all corrective actions are shifted. Any help or advice would be appreciated!


title: Simple extrapolation of web traffic

text: I have the following data for web traffic in Excel   * June 88,147   * July 111,839   * August 93,148   * September 93,069   * October 101,881   * November 97,345 I can do a chart with a trend line and a linear trend shows a slight upward trend, I would like the simplest way (stats is not my field!) to project / guestimate what it will be like in the future The trend line in the chart gives me `y=14.321x - 491537`, but I don't know to work that into a formula I tried `=14.321*E8-491537` (`e8` is November's value) and it gave me 902,545 which does not fit the trend at all! The goal is a rough estimation of the month by which the figures will be up to 130k+ (presuming my slight upwards trend is correct!)


title: Life expectancy of home appliances

text: I hope this is the right place where to formulate this question. In a basic course in statistics we learn how to calculate the mean life time of an electronic device. Home appliances are electronic device of course, but all that we know when we buy one is the time length of warranty. What is the relation of mean life time and warranty time in the home appliances? or better how is calculated the warranty time? What I am interested is how to calculate the probability that an appliance will brake within the warranty time.


title: Understanding the exponential distribution

text: I'm trying to wrap my head around the Exponential distribution and the meaning of its parameter. The parameter is the rate, right? So take, e.g., $$X\sim \exp(0.05)\,.$$ Now the probability of failure during the first time period is: $$P(X\le1)=1-P(X>1)=1-e^{-0.05}=0.04877\,.$$ Now I can do the math and get the correct result and so on, but cannot wrap my head around why the result should be slightly less than 5%, rather than exactly 5%. I don't get the intuition behind it.


title: Understanding survival at time function

text: This question is related to a few others (Here,Here) on the topic as I have been searching for information. Hopefully this one is sufficient. 1) I am seeing differences in the relationship between survival at time t, S(t) and the hazard at time t, h(t). In Singer and Willet page 337 (10.5) they define the estimated survival $\hat{S}(t_{j})$ = $[1-\hat{h}(t_{j})][1-\hat{h}(t_{j-1})][1-\hat{h}(t_{j-2})]....[1-\hat{h}(t_{1})]$. They also define a life table as such: ![enter image description here](http://i.stack.imgur.com/GU37f.jpg) In other sources Will Potts and others, they define $\hat{S}(t)$=$\prod_{t_{j}< t}(1-h(t_{j}))$ implying that the value of the hazard at time t is not included in the calculation for Survival at time t. Which is correct? Is the difference in how the life table is set up? 2) How would you calculate the mean time time (restricted since there is censoring at the end of the collection period). Would you sum all $\hat{S}(t)$ (including the 1 at time zero) to arrive at 7.7? You always her it is the "area under the survival curve". 3) How would you calculate the residual mean lifetime for those that survived to the end of, say year 4? Is it summing all $\hat{S}(t)$ from 0.6209 through 0.41123 divided by 0.6209 Using the fact that the formula for the mean residual survival time mrl(u) = E(T -u | T>u) is given by $\frac{1}{S(u)}\int_u^{\infty} S(t)\, dt$ and our "u" is 4 here (since we are saying T>4 in the conditional expectation)?


title: Post-hoc analysis of 3 sample prop.test

text: I am looking for the correct way of performing post-hoc analysis on the results of a 3-sample test for equality of proportions and haven't been able to find anything specific to this test. The data I am analyzing looks like this:                              pos.late neut.late neg.late         pos.early        88        12        4         neut.early       23       112        7         neg.early        58       150        5      The row names are how individuals responded in an early time period, and the column names are how they responded in a late time period. The values in each cell represent how many individuals showed that particular response over time (e.g. 88 individuals were positive in the early time period and remained positive in the late). I performed a `prop.test()` in R for each row – here is an example for the `pos.early` row:               > prop.test(c(88,23,58), c(88+23+58, 88+23+58, 88+23+58))              3-sample test for equality of proportions without         continuity correction          data:  c(88, 23, 58) out of c(88 + 23 + 58, 88 + 23 + 58, 88 + 23 + 58)     X-squared = 56.3609, df = 2, p-value = 5.773e-13     alternative hypothesis: two.sided     sample estimates:        prop 1    prop 2    prop 3      0.5207101 0.1360947 0.3431953       There is clearly a significant difference, but what post-hoc test is best used with the `prop.test()` to determine which of the `pos.early`–`pos.late`, `pos.early`–`neut.late`, `pos.early`–`neg.late` pairings are significant?


title: Why is the variance of my random effect negative?

text: I have two variables, `response` and `group`, and I fitted the model in which `group` is a random effect.               proc mixed data=myData nobound;        class group;        model response =;        random group;     run;      However, SAS returns a negative variance for the random effect... I have noticed that, unexpectedly, the variability of `response` within `group` may be larger than the overall variability of `response`. I guess this is the reason why SAS returns a negative variance. Is it your opinion too? If so, how could this be seen from a mathematical point of view? **DATA**               group   response     1   70.8     1   64.2     1   67.8     1   62     1   60      1   66.6     1   67.8     1   60     1   70     1   63.3     1   37     1   53.8     2   82.7     2   65.5     2   67.7     2   51.6     2   53.8     2   64.3     2   72.4     2   48.3     2   54.8     2   64.8     3   88.5     3   69     3   60.9     3   71.4     3   52     3   58.3     3   63.6     3   51.9     3   52.4     3   52.2     3   65.4


title: Are there any pre-generated numeric series with known statistical properties?

text: I'm trying to test a program that claims to calculate some statistical properties of a numeric sequences (like mean, median, standard deviation etc). Simple and short sequences pass the test well, but I'd like to challenge the program with either long sequences, or sequences with large and small values (to test overflows/underflows) etc. So, are there any software or online resource or other source of different numeric sequences with known statistical properties?


title: How exactly does null hypothesis testing work

text: I'm trying to gain a better understanding of how null hypothesis testing works. I have 3 questions related to the code below:   1. Am I right in saying the probability of each t value in `t_distribution` should follow the tables of critical values for t distributions found in statistical text books?   2. If yes, how using R code, can I show `t_distribution` follows the tables of critical values for t distributions found in statistical text books   3. Who exactly made the tables of critical values for t distributions found in statistical text books? # make 2 populations with (almost) zero difference               population1 <- rnorm(10000, 3, 2)     population2 <- rnorm(10000, 3, 2)      # sample 100 values from each population               library(plyr)     population1_5000_samples <- llply(1:5000, function(x) sample(population1, 100, replace = F))     population2_5000_samples <- llply(1:5000, function(x) sample(population2, 100, replace = F))      # do 5000 t-tests and extract t values               t_distribution <- laply(1:5000, function(x) t.test(population1_5000_samples[[x]], population2_5000_samples[[x]])$stat)      # distribution of t values               library(ggplot2)     ggplot(as.data.frame(t_distribution), aes(t_distribution)) + geom_histogram()      ![enter image description here](http://i.stack.imgur.com/FrODv.jpg)


title: How to obtain the model behind a simulator?

text: I am looking for an useful statistical approach or analysis tool in order to understand the data obtained from an aeroelastic simulator of wind turbine dynamics. In this case, the simulation provides data about the forces that a structure has to resist against the wind. Imagine a large structure such as a wind turbine which its base is loaded (with forces and moments) at each time step. This structural loads depend on the loads of previous instances of time (autocorrelated) and more variables such as wind intensity, angle of the attack, lift coefficients, relative speed, blade deflection of different sections of different blades, and other variables. The dataset is a time series of 12000 observations with 195 variables with many dependencies among them but some can be pointed out by an expert aeronautical engineer. She believes that some dependencies are non linear. **The main goal is to quantify the grade of dependency, in other words, what variables influence more the response variable (loads at the tower base)**. Later on I would also like to estimate the model using a data set from a different structure and see whether more influence variable are the same. I expect to see the same group of variables because the simulator’s algorithm is always the same. I thought about some approaches but I want to know your experiences in similar problems. - Generalized Linear Model or design of experiments, I have experience using them but I don't know how to capture the autocorrelation and the non-linearity. - Multivariate Time series models, I do not have too much experience, but I believe they can capture the dependencies. - What about more flexible and opaque models such as neural networks. I am a recent graduate statistician and we hardly seen such complex problems during the degree. However, I am very confident with R and I'd love to know a good approach in order to start my research at CRAN and Wikipedia. My starting point is the Time Series View Task.


title: Pruning the tree (CART) using R - different plotcp and printcp results

text: I run exactly the same analysis several times in a row. The original (not pruned tree) is the same but the results given by "plotcp" and "printcp" are different every time. So, how can I choose the best cp to prune the tree ? Thanks for response below !


title: Group effects where there is a single group

text: Suppose I have a table of subjects and measurements as follows               subject measurement     s1        1     s1        2     s1        3     s1        4     s1        5     s1        6     s1        7     s1        8     s1        9     s1       10     s2       10      Now it is unlikely that the true value of s2 is 10 as observed. It is much more likely to be closer to the mean of the s1 measurements. How do I discover the true value? I've discovered that if I construct the following table:               group subject measurement     g1  s1  1     g1  s1  2     g1  s1  3     g1  s1  4     g1  s1  5     g1  s1  6     g1  s1  7     g1  s1  8     g1  s1  9     g1  s1  10     g1  s2  10     g2  s1  1     g2  s1  2     g2  s1  3     g2  s1  4     g2  s1  5     g2  s1  6     g2  s1  7     g2  s1  8     g2  s1  9     g2  s1  10     g2  s2  10      and run the following R command                fit = lmer(measurement ~ 1 + (subject|group))  I get a plausible answer. Is it correct? Is there a better way? In R? I had to double the table up to coerce R to solve the problem as it doesn't seem to like a single group.


title: Practical difference between "sampling" and "re-sampling with replacement"

text: What is the difference between **"sampling"** and **"re-sampling with replacement"** from a practical point of view (I mean how do I code that) ? When one say: **"sampling $\tilde{x}_t^{(i)}$ from $p(x_t | x_{t-1}^{(i)})$"** , this mean that we generate $\tilde{x}_t^{(i)}$ from the conditional probability distribution $p(x_t | x_{t-1}^{(i)})$, so this is the next random value to generate from a uniform or a gaussian distribution for example, right ? Now, when one say: **"re-sampling with replacement N values $\{x_t^{(i)}; i=1,\dots,N\}$ from the set $\{\tilde{x}_t^{(i)}; i=1,\dots,N\}$ according to the associated weights/probabilities $\{w_i; i=1,\dots,N\}$"**. Does this mean that we choose each value $\tilde{x}_t^{(i)}$ according to the probability $w_i$, and repeat this until we get N values ? So how do you do that practically in a pseudo code ? **EDIT:** I understand that re-sampling here means reselecting the values from the weighted set, by choose randomly, but according to the weights. So if we had three values, with weights 0.99, 0.005 and 0.005, then the most likely outcome of resampling would be three values of the first type. But how do I code this in a general case ?


title: Intuition for Gamma-Poisson / Negative Binomial

text: I have a data set which intuitively seems Poisson-like, but it's overdispered. So I'm investigating negative binomial. From this question and this page I understand one way of viewing this is to state that we are uncertain about the density parameter $\lambda$ in Poisson, and assume $\lambda \sim \text{Gamma}(\alpha, \beta)$, in which case we get a negative binomial. I get the math here, but I don't understand why we would think that $\lambda$ is gamma-distributed. For example, the wordpress link says that when $\lambda$ represents "likelihood to default on a loan", then it would be gamma- distributed. But I don't really understand why. Is there a way I can determine if my density is gamma distributed?


title: Restricted Boltzman Machine Non-Hidden Layer Approach

text: An RBM is defined by the joint probability distribution $$p({\bf x},{\bf h})=\exp(-E({\bf x},{\bf h}))/Z$$ where $$E({\bf x},{\bf h})=-{\bf h}^TW{\bf x} - {\bf c}^T{\bf x} - {\bf b}^T{\bf h}$$ and $$Z=\sum_{{\bf x},{\bf h}}\exp(-E({\bf x},{\bf h}))$$ This gives a marginal distribution $$p({\bf x}) = \exp({\bf c}^T{\bf x} +\sum_j^H\log(1+b_j + W_j{\bf x}))/Z$$ where $H$ is the number of hidden units i.e. length of ${\bf h}$ and $W_j$ is the $j$th row of $W$. We can then see that by maximising the product of this over our training data with respect to our parameters, will tend to pick out common features in the input data due to the effect $W_j {\bf x}$ term, with the $Z({\bf b}, {\bf c}, W)$ (to make clear that $Z$ is a function of our parameters) acting as a normalisation. My question then is, why do we even bother initialising the model using the ${\bf h}$ but rather just maximise the function $$f({\bf x}) = \exp({\bf c}^T{\bf x} +\sum_j^H\log(1+b_j + W_j{\bf x}))$$ subject to $\int_{\bf x}f({\bf x})=1$. Using the hidden units allows us to use Gibbs Sampling but do we not have an analagous method for our function $f$ i.e. considering something along the lines of $W^TW{\bf x}$. Is introducing hidden units really the only way to train an RBM efficiently so that we can use contrastive divergence?


title: Normalizing error variance scale

text: Imagine we have two sets of data m1 and m2. m1 components have the dimension of [m^2/s] and m2 components are measured in [cm^2/s].   e.g.               m1={1, 2, 3, 4, 5, 6, 7, 8, 9, 10}          m2={31, 32, 33, 34, 35, 36, 37, 38, 39, 40}      Assuming that B is the proportionality constant that is estimated as a straight line fitted to the points with the least square method, it can be estimated by               NonlinearModelFit[Transpose[{10^-4 m1, m2}], a x, {a}, x,       VarianceEstimatorFunction -> (Mean[#^2] &)]          52857.14      **But, the estimated error variance is NOT non-dimensional. How can I normalize the error variance scale?** e.g.               Sqrt[NonlinearModelFit[Transpose[{10^-4 m1, m2}], a x, {a}, x,         VarianceEstimatorFunction -> (Mean[#^2] &)]["EstimatedVariance"]]          13.8873               Sqrt[NonlinearModelFit[Transpose[{m1, 10^4 m2}], a x, {a}, x,         VarianceEstimatorFunction -> (Mean[#^2] &)]["EstimatedVariance"]]          138873.      It looks quite inconsistent if I write B=52857.14 +/- # [cm^2/s] while B has no dimensions!


title: Joint distribution of two dependent variables

text: I was wondering how to (or if it is even possible) find the continuous joint distribution between two random variables $x$ and $y$ when you know the continuous marginal density distributions of both $x$ and $y$, and we know there is a correlation between the two random variables, i.e., ${\rm Cor}(x,y)$ = $\rho$ which is of unknown magnitude?


title: Different number of time-points in functional data analysis

text: I have a data set that I wish to analyze using functional data analysis methods. Data consists of repeated measures of some characteristic on a number of inviduals. I have the time of the measurements and these differ between individuals. So far, I see no problem in analyzing this, e.g. by using ready- to-use ${\tt R}$ packages. However, the number of measurements per individual differs as well for my data. How can I handle this? Can anyone point me in the direction of literature on this set-up or any packages that can handle it? What I want to do with data to start with is basically something like functional linear model (functional response and scalar predictors). I'm aware that I could smooth the data using a very light smoothing procedure (using many basis functions and no or a very small penalty parameter) and use these smoothed curves for my analysis. But is this the best way to approach the problem? It doesn't take into account that some curves are determined with more precision than are others. An additional problem to this approach is that some of the individuals only have measurements in a subinterval (say, [200, 1000]) of the full interval of the majority of the curves (say, [0, 1000]). Thus, the smooth would be very poorly estimated in those regions. Naturally, I could restrict attention to the interval [200, 1000] but this seems like a waste of data. To sum up my question, what would you do in this case? Any thoughts and hints will be much appreciated.


title: What are the primary differences between z-scores and t-scores, and are they both considered standard scores?

text: We are currently converting student test scores in this manner :               ( ScaledScore - ScaledScore Mean ) / StdDeviation ) * 15 + 100       I was referring to this calculation as a z-score, I found some information that convinced me that I should really be referring to it as a t-score as opposed to a z-score. My boss wants me to call it a "Standard Score" on our reports. Are z-scores and t-scores both considered standard scores? Is there a well known abbreviation for "Standard Score"? Can someone point me to a reference that will definitively solve this issue.


title: Clustering with asymmetrical distance measures

text: How do you cluster a feature with an asymmetrical distance measure? For example, let's say you are clustering a dataset with days of the week as a feature - the distance from Monday to Friday is not the same as the distance from Friday to Monday. How do you incorporate this into the clustering algorithm's distance measure?


title: Do non-significant correlates of a DV, which is significant, suggest anything about the effect of that DV?

text: Let's assume one has an analysis in which there are multiple correlated DVs (average correlation .46) being examined in separate univariate analyses (e.g. t-tests; insufficient df and frustration tolerance to restructure the data fro a MANOVA) in relation to a single IV. All except one of these DVs is not significantly affected by an IV manipulation. I was thinking that in general, dependence between DVs would tend to globally trend analyses to be significant or non-significant. However, when a DV with a legitimate effect is correlated with a variable which has no effect, wouldn't that correlation represent an association between them which is a source of error variance in the DV with the legitimate effect? Thus, if anything, the presence of significance in a DV associated with non-significant DVs should be taken as a sign of a true effect on that DV? In short:   1. What, if anything, do the non-significant tests tells us about the one that is significant?    2. Is there a relationship between this question and the notion of suppressor variables in regression?


title: Sampling data to have specific mean and standard deviation

text: I have a data that I want to sample such the resultant distribution of values should have specified mean and standard deviation. I can think of rejection sampling to achieve this however that seems to be overkill as rejection sampling tries to match complete distribution while I need only first two moment to be equal to specified value. There have to be easy way to do this and also multiple ways. The current data have distribution               Min.    1st Qu.     Median     Mean     3rd Qu.    Max.        0.0020  0.2120      0.3880     0.4298   0.6040     1.0000       I want it to have mean of 0.2236 and standard deviation of 0.3180463.


title: Getting started with dynamic neural networks

text: I was wondering if anybody had some good references on dynamic neural networks. The goal is to have a neural network that responds in real-time (and is governed by say a system of DEs) instead of at discrete time-steps. The prototypical application is a robot: we have a series of sensors that are constantly feeding data into the input neurons, and we have output neurons that are constantly telling the motors what to do. We want to optimize this robot to do some tasks. Is there a nice recent survey on dynamic neural nets? In particular one that discusses recurrent DNN, feedforward DNN, and DNN with and without learning. Alternatively, if continuous dynamic neural networks are not the right approach for this sort of task, then suggestions of other types of neural nets are also welcome. ## Notes I also cross posted this question to MetaOptimize.


title: Determining whether to scale up pilot study

text: I am analyzing the results of a pilot study which applied an intervention to encourage uptake of a financial relief program among distressed homeowners. The difference in outcome is 37% among the treatment group vs. 29% in the control group. At a sample size of n=400 (across both the treatment and control group), a chi-square test yields a 2-tailed p-value of 0.11. (Fisher's test yields p=0.1357). Would you recommend that the pilot study be scaled up? My inclination is to say that this is not statistically significant, so the answer would be no, but I'm unsure because my stats are fairly rusty. Thanks!


title: Fitting a line to a log-log plot

text: I have some data that I'm playing around with; for simplicity, let's suppose the data contains information on number of posts a blogger has written vs. number of people who have subscribed to that person's blog (this is just a made-up example). I want to get some rough model of the relationship between # posts vs. # subscribers, and when looking at a log-log plot, I see the following: ![enter image description here](http://i.stack.imgur.com/Jb03e.png) This looks like a rough linear relationship (on the log-log scale), and quickly checking the residuals seems to agree (no apparent pattern, no noticeable deviation from a normal distribution): ![enter image description here](http://i.stack.imgur.com/aU7Sf.png) ![enter image description here](http://i.stack.imgur.com/Zqo5G.png) So my question is: is it okay to use this linear model? I know vaguely that there are problems using linear regressions on log-log plots to estimate power law distributions, but my data isn't a power law _probability distribution_ (it's simply something that seems to roughly follow a $subscribers = A * (postings) ^ k$ model; in particular, nothing needs to sum to 1), so I'm not sure if the same critiques apply. (Perhaps I'm over-correcting at the mention of "log-log" and "linear regression" in the same sentence...) Also, all I'm really trying to do is to:   * See if there are any patterns to the blogs with positive residuals vs. blogs with negative residuals   * Suggest some rough model of how subscribers are related to number of postings.


title: What happens when I include a squared variable in my regression?

text: I start with my OLS regression: $$ y = \beta _0 + \beta_1x_1+\beta_2 D + \varepsilon $$ where D is a dummy variable, the estimates become different from zero with a low p-value. I then preform a Ramsey RESET test and find that i have some misspesification of the equation, i thus include squared x: $$ y = \beta _0 + \beta_1x_1+\beta_2x_1^2+\beta_3 D + \varepsilon $$   1. What does the squared term explain? (Non-linear increase in Y?)    2. By doing this my D estimate does not vary from zero any more, with a high p-value. How do i interpret the squared term in my equation (in general)? Edit: Improving question.


title: Where can I find an implementation of Hinton's original Boltzmann Machine?

text: I've been trying to implement the Boltzmann machine 4-2-4 encoder that appeared in A Learning Algorithm for Boltzmann Machines. but I am unable to find clear pseudocode for doing it or more specific details. I have tried "Neural networks: algorithms, applications, and programming techniques" and Fausett books so far. Where can I find an implementation of the encoder or more specific details about it?


title: Calculating L moments of a standard normal

text: As part of my research, I am having trouble calculating the 2nd L moment $\lambda_2$ and the third and fourth moment ratios $\tau_3$ and $\tau_4$ of the standard normal distribution, where L moments are defined as they are in this paper (Hosking 1990), page 3. Take $\lambda_2$ for example. $\lambda_2 = \frac{1}{2} [ EX_{2:2} - EX_{1:2}] \\ \implies \lambda_2 = \int_{-\infty}^{\infty} xF(x)f(x) dx - EX + \int_{-\infty}^{\infty} x F(x) f(x) dx \\ = \ \cdots \\ = \frac{1}{\sqrt{2 \pi}} \int_{-\infty}^{\infty} x e^{-x^2/2} dx + \frac{1}{\sqrt{2 \pi}} \int_{-\infty}^{\infty} x \ erf(\frac{x}{\sqrt{2}}) e^{-x^2/2} dx \\ = \frac{1}{\sqrt{2 \pi}} \int_{-\infty}^{\infty} x \ erf(\frac{x}{\sqrt{2}}) e^{-x^2/2} dx$ There must be an error in my calculation as this is supposed to be equal to $\frac{1}{\sqrt{ \pi}}$. $\lambda_3 = \frac{1}{3}E(X_{3:3} - 2X_{2:3}+X_{1:3}) = \int_0^1 x(F) (6F^2 - 6F + 1) dF$ $\lambda_4 = \frac{1}{4}E(X_{4:4} - 3X_{3:4}+3X_{2:4} -X_{1:4}) = \int_0^1 x(F) (20F^3 - 30F^2 + 12F - 1) dF$ $\tau_3 = \frac{\lambda_3}{\lambda_2}$ $\tau_4 = \frac{\lambda_4}{\lambda_2}$ $\tau_3$'s formula is rather involved. EDIT: I was able to figure it out. Instead of working through the integral as an ugly function of exponentials, x's and erf, work through it as a function of a power of F multiplied by the inverse erf term. The latter is analytically tractable. I am not sure the former is.


title: The role of scale parameter in GEE

text: I am learning the generalized estimating equations (GEE) and the `geepack` R package. There are some questions that I am a little confused. In a GEE-constructed model, we have $Var(Y_{it})=\phi_{it}\cdot V(\mu_{it})$, where $\phi$ is the scale parameter. We further decompose $Var(Y)$ into $V^{1/2}R(\alpha)V^{1/2}$ where $\alpha$ is the correlation parameter. **Three link-models** are specified in `geepack`, for $\mu,\phi,\alpha$, respectively. See this PDF file for details. (1) In GEE1, can I say we only need to make sure that the mean structure is correctly specified, i.e. the link model $g(\mu)=X\beta$ is correct, while it doesn't matter whether the link models for $\phi$ and $\alpha$ are correctly specified? (2) By default, the `geese.fit` function makes the scale value `scale.value = 1.0` -- does it say $\phi=1.0$? It is understandable that the default `alpha=NULL` as people can specify different correlation structures and the program will assign appropriate `alpha` values accordingly. My question is: how often people try to explicitly model the scale parameter $\phi$? (3) This question is closely related to (2) about the scale parameter $\phi$. Recall the variance function is $Var(Y_{it})=\phi_{it}\cdot V(\mu_{it})$. In the Gaussian case, we have $\phi=\sigma^2$ and $V(\mu_{it})=1$; in the binomial case, we have $\phi=1$ and $V(\mu_{it})=\mu_{it}(1-\mu_{it})$; in the Poisson case, we have $\phi=1$ and $V(\mu_{it})=\mu_{it}$. Can I say that, in the negative binomial case, $\phi=1$ and $V(\mu_{it})=\mu_{it}+\varphi\mu_{it}^2$? Here $\varphi$ is the NB2 dispersion parameter, NOT the scale parameter $\phi$ in GEE. Thank you very much!


title: Distribution function

text: Find (without using MGF) the mean and variance. $$f(x) = \exp(-kx)x^{(r-1)}k^r/(r-1)!\ \text{ for }\ x>=0$$ $$f(x) = 0\ \text{ for }\ x<0$$ $r$ positive integer, $k>0$


title: What is "feature space"?

text: What is the definition of "feature space"? For example, When reading about SVMs, I read about "mapping to feature space". When reading about CART, I read about "partitioning to feature space". I understand what's going on, especially for CART, but I think that there is some definition that I have missed. Is there a general definition of "feature space"? Is there a definition that will give me more insight into SVM kernels and/or CART?


title: Convolution of two square pulses

text: I want to convolve two square pulses. One pulse say X starts from 0 to R.Other pulse say Y starts from -R to 0.Now to convolve these two pulses I just flip the pulse X and then move it to -infinity.Now, to convolve I will pass this pulse through pulse Y. The result will be a triangular waveform which starts from -R+(1/2) and ends at 0+(1/2).Am i right?


title: For "Was this page helpful" data, should I take response rate into account?

text: To help me focus on which pages on a site to improve, I'm looking at user feedback to the question "Was this page helpful?" (Answers are "Yes" or "No".) The response rate (responses divided by unique pageviews) varies across different pages. And in fact the pages with low response rate also tend to have a higher ratio of negative ratings. This suggests to me that many users who don't find a page helpful just leave without rating it. (I have looked at other factors that could have affected response rate but the correlations are very weak.) I understand that it's problematic to generalize the result for pages that have a low response rate. But I'm reluctant to exclude these pages, since it seems that they may be the ones most in need of attention. In view of this, is it reasonable to ignore response rate per page and simply look at the positive to negative ratings ratio? Or is there a way I could account for response rate and still include these pages?


title: Random factor nested in two fixed factors

text: I have read Random effect nested under fixed effect model in R, but I have a doubt: My data is on germling survivorship, I have `Temperature` as a fixed factor (2 levels), `pH` as a fixed factor (2 levels), and I have `tank`, which I included to check for tank effect (2 tanks per treatment combination, in a total of 8 tanks). How can I make a model which nests tank within the temp/pH fixed effect? From a previous post I read, we could do: Using `aov` like this:               a0 <- aov(survivorship ~ pH*temp + Error(tank), data=d)     summary(a0)      or `lme` like this:               library(nlme)     m1 <- lme(survivorship ~ pH*temp, random=~1|tank, data=d)     anova(m1)      I have tried these two models, but I got two different answers. Can anyone explain why? Is this correct, by using `aov` like this, is `tank` nested in `pH:temp` interaction?               aov(survivorship ~ pH*temp + Error(tank)


title: How can we get the confidence (or probability) that a data-point belongs to an unknown class?

text: Using any supervised classifier, we can usually get the probability that a data point $x$ belongs to each class $y_i$, i.e. $P(y_i|x)$. However, in the case where the data-point x may belong to none of the known classes, how can we get the probability $P(?|x)$, which means that x belongs to an unknown class with probability $P(?|x)$ ?


title: How to choose predictors for regression model

text: I want to predict reaction times using several personality scores. I have 9 different personality scores. My sample consists of 23 participants. Aren't there too many predictors if I put all 9 in one model, especially regarding the small sample size? If not, how could I choose predictor? The theory does not state any predictor as more important than the others. the predictors also have low intercorrelation, i.e., multicolinearity is low even with 9 predictors in the model. One thought of mine would be to see which personality score has a significant correlation with the reaction times and only use them in a model. What is your opinion?


title: Detecting steps in time series

text: I've attached a picture of the time series I'm talking about. The top is the original series, the bottom is the differenced series. Each data point is a 5 minute average reading from a strain gauge. This strain gauge is placed on a machine. The noisy areas correspond to areas where the machine is turned on, the clean areas are when the machine is turned off. If you look at the area circled in red, there are anomalous steps in the reading that I would like to be able to detect automatically. I'm completely stumped on how I might be able to do this - any ideas? ![enter image description here](http://i.stack.imgur.com/FbGXD.png)


title: Interpreting seasonality with ACF and PACF

text: I have a dataset where empirical intuition say I should expect a weekly seasonality (i.e., the behavior in saturday and sunday is different from the rest of the week). Should this premise be true, shouldn't an autocorrelation graph give me bursts at lag multiples of 7? Here's a sample of the data:               data = TemporalData[{{{2012, 09, 28}, 19160768}, {{2012, 09, 19},          19607936}, {{2012, 09, 08}, 7867456}, {{2012, 09, 15},          11245024}, {{2012, 09, 04}, 0}, {{2012, 09, 21},          24314496}, {{2012, 09, 12}, 11233632}, {{2012, 09, 03},          9886496}, {{2012, 09, 09}, 9122272}, {{2012, 09, 24},          23103456}, {{2012, 09, 20}, 25721472}, {{2012, 09, 11},          12272160}, {{2012, 09, 25}, 21876960}, {{2012, 09, 05},          7182528}, {{2012, 09, 16}, 11754752}, {{2012, 09, 23},          23737248}, {{2012, 09, 26}, 20985984}, {{2012, 09, 10},          12123584}, {{2012, 09, 06}, 9076736}, {{2012, 09, 17},          20123328}, {{2012, 09, 18}, 20634720}, {{2012, 09, 22},          23361024}, {{2012, 09, 14}, 11804928}, {{2012, 09, 07},          9007200}, {{2012, 09, 02}, 9244192}, {{2012, 09, 13},          11335328}, {{2012, 09, 27}, 20694720}, {{2012, 10, 26},          12242112}, {{2012, 10, 15}, 10963776}, {{2012, 11, 09},          9735424}, {{2012, 10, 08}, 10078240}, {{2012, 10, 31},          10676736}, {{2012, 10, 20}, 11719840}, {{2012, 11, 05},          10475168}, {{2012, 10, 01}, 9988416}, {{2012, 10, 24},          11998688}, {{2012, 10, 12}, 10393120}, {{2012, 10, 23},          11987936}, {{2012, 10, 19}, 11165536}, {{2012, 10, 04},          9902720}, {{2012, 11, 16}, 10023648}, {{2012, 11, 21},          10047936}, {{2012, 10, 10}, 10205568}, {{2012, 11, 08},          9872832}, {{2012, 10, 21}, 12854112}, {{2012, 11, 04},          10485856}, {{2012, 10, 07}, 9565248}, {{2012, 09, 30},          9784864}, {{2012, 10, 29}, 12880064}, {{2012, 11, 10},          8945824}, {{2012, 11, 15}, 9870880}, {{2012, 09, 29},          9718080}, {{2012, 10, 18}, 10992896}, {{2012, 10, 06},          9319584}, {{2012, 11, 03}, 9077024}, {{2012, 10, 03},          10537408}, {{2012, 11, 22}, 9853216}, {{2012, 10, 11},          10191936}, {{2012, 10, 22}, 12766816}, {{2012, 11, 07},          9510624}, {{2012, 11, 14}, 9707264}, {{2012, 10, 28},          12060736}, {{2012, 11, 19}, 10946880}, {{2012, 11, 11},          9529568}, {{2012, 10, 09}, 9967680}, {{2012, 10, 17},          12093344}, {{2012, 11, 20}, 10520800}, {{2012, 10, 05},          9619136}, {{2012, 10, 25}, 11484288}, {{2012, 11, 17},          9389312}, {{2012, 10, 30}, 12078944}, {{2012, 10, 14},          9505984}, {{2012, 10, 02}, 9943648}, {{2012, 11, 24},          9458144}, {{2012, 11, 02}, 10082944}, {{2012, 11, 01},          11082912}, {{2012, 10, 13}, 9117632}, {{2012, 11, 23},          10253280}, {{2012, 11, 12}, 10240672}, {{2012, 11, 06},          9723456}, {{2012, 11, 13}, 9806880}, {{2012, 10, 16},          12368896}, {{2012, 11, 18}, 9632800}, {{2012, 10, 27}, 10606656}}]      ... and the ACF: ![enter image description here](http://i.stack.imgur.com/ARQ2A.png) ... and the PACF: ![enter image description here](http://i.stack.imgur.com/1b0Us.png)


title: continous moderator in a within subject design

text: I'm a first year PhD student and I am a bit at a loss regarding which test to use, and how to calculate power for a grant proposal. Any help will be much appreciated! I'm running all my participants in two conditions (positive/negative), so a within-subject design My dependent variable is continuous I'm looking at a moderator effect, the moderator being continuous also So I'm basically seeing if the the effect of condition on my dependent variable is moderated by a continuous variable I can't seem to figure out which test to use. Thanks in advance


title: Sampling in discrete distribution

text: At the moment I'm working with an algorithm that calculates probabilities $Z$ for a finite sequence of data points $X= \left\{ x_1,x_2,...,x_n \right\}$ using an expression like:   $p(Z_i| X_i,\theta) = \frac{1}{\phi}f(X_i,\theta)$,   where $\phi$ is a normalising constant such that $\sum_n p(Z_n|X_n) = 1$. (i.e. $f(X_i,\theta)$ gives us a number, we have a sequence of numbers and we normalise over all of them by summing and dividing). The next part of the algorithm chooses one of the $x$'s with weighted probability.   At the moment, it is possible to enumerate the sample space and calculate $\phi$, but is it possible to have some kind of MCMC approach, to avoid having to do this? So sampling a few positions, and evaluating $p(Z_i|X_i)$ only at those positions? I can see how MCMC works in the continuous case, but not sure if it's applicable here.   If it makes any difference, the values of $f(X_i,\theta)$ can have a lot of variation (e.g. 0.08, 0.70, 0.26, 64.00, 0.79, 7.11, 0.01, etc.) so areas of high probability can be separated by regions of low probability. Thanks in advance.


title: standard errors for marginal models with no clustering

text: Just a question about how robust standard errors are fitted in marginal models with GEE: I noticed that if I try to fit a model with clusters of size 1 (i.e. no clustering, therefore no intra-cluster correlation), the variance is not the same as the one I would get with a standard GLM. Here is an example with R:               set.seed (555)          X <- rbinom (100, 1, .5)     Y <- rbinom (100, 1, .5)     ID <- as.factor (1:100)          ##### standard glm #####     fit.glm <- lm (Y ~ X)     coef (summary (fit.glm))          #             Estimate Std. Error t value Pr(>|t|)     # (Intercept)    0.691     0.0651   10.61 5.77e-18     # X             -0.246     0.0971   -2.54 1.27e-02          ##### marginal model #####     library (geepack)     fit.gee <- geeglm (Y ~ X, id=ID, corstr="independence")     coef (summary (fit.gee))          #             Estimate Std.err   Wald Pr(>|W|)     # (Intercept)    0.691  0.0623 122.94   0.0000     # X             -0.246  0.0968   6.48   0.0109      I guess the answer may be that only glm are fitted with maximum likelihood but it does not seem very relevant. Any idea? Thanks in advance.


title: Why two-step dynamic panel data estimator is better then one-step?

text: Could someone please give an explanation of why two-step dynamic panel data estimator is better then one-step? Can't quiet understand it... For example, from `xtabond2` manual in STATA: "two-step estimator is asymptotically eﬃcient and robust to whatever patterns of heteroskedasticity and cross-correlation the sandwich covariance estimator models. Historically, researchers often reported one-step results as well because of downward bias in the computed standard errors in two-step. But as the next subsection explains, Windmeijer (2005) has greatly reduced this problem". So, why is two-step estimator is better then one-step estimator with robust option estimator?


title: Test for non random-walk

text: I have a set of sequential data around of around 12k points. Each is a 1 or -1. They are split something like 53%/47%. I want to test the hypothesis that this sequence comes from a random walk with a constant parameter. What sort of statistical test should I use? (I'd also be grateful for any book recommendations in this general area) Thanks very much


title: Adjust standard errors for within correlation

text: I am trying to replicate a table and in one of the notes it's written that **'standard errors are adjusted to account for the within-analyst correlation of the observations'** I am running my regressions in Matlab and obtaining normal SEs. I have been searching how to adpat my SEs to this specific situation and what's the meaning of it but I have not been successful. Can someone please help me? Thank you.


title: Weighted Standard Deviation vs. Weighted Regression Standard Error

text: This is an R question. I have $n$ observations of the variable $y$, with each observation, $i$, weighted by $w(i)$. Each $w(i)$ weight falls between zero and 1, with the sum of the weights, $m<n$. I am testing whether the weighted mean of $y$ equals zero, with $std(mean)$ denoting the weighted sample standard deviation of the mean, computed using `cov.wt` and dividing by `sqrt(m)`. At a later time, I would like to run a weighted regression of $y$ against some explanatory variables. However, as an intermediate step, I have run the regression `lm(y~1,weights=w)`, which returns an intercept equal to the weighted mean of $y$. However, the standard error of the intercept is based on $n-1$ degress of freedom and, therefore, is smaller than the standard deviation of the weighted sample standard deviation of the mean, computed using `cov.wt` and dividing by `sqrt(m)`. Is there a way to run a regression in R in which the degrees of freedom in the simple regression above would be $m-1$ rather than $n-1$? Note, I have $m$ "cases" of the variable $y$. In most cases, there is only one observation. But if there were two observations in a given case, I would give each a weight of $1/2$, if three observations in a case, each would get a weight of $1/3$, etc.


title: Testing a 2x2 contingency table: male/female, employed/unemployed

text: I major in science, and my knowledge of statistics is rather superficial. ## Problem I had to find a data set and analyze it to the best of my ability as an assignement for my statistics course. This is no longer an assignment, I just need help in interpreting why I did my analysis badly and what I should have done instead. **I used a categorical data set of employment rates in New Zealand, planning to arrange it in a 2x2 contingency table and use Pearson's chi-squared test and Fisher's exact test to test whether gender correlates with employment.** ## What I want to answer   1. Understand why I cannot use chi-squared test and Fisher's exact test for this problem and learn what I should have used instead. "Odds-ratio as a function of time", I assume? Any useful links on how do that, perfectly in R?   2. Understand the "sequential correlation" comment regarding the first part of the assignment and what exactly should I have done. ## Way to help me #1 (shorter) That's how our data looks (based on a census):                                Male     Female     Employed      1201600    1060200     Unemployed      73300      75000      I did a chi-squared test and a Fisher's exact test in R, assuming that the obtained p-value will tell me the probability of such a distribution of jobs (or one more extreme) given that the null is true (that males and females have equal chances of getting a job). I got a very small p-value, and Fisher's test gave me odds ratio of 1.16, meaning that there is a correlation, and specifically males are 16% more likely to find a job in NZ. However, according to my lecturer, I used these tests inappropriately. I didn't quite understand why, but I think he was saying that these tests assume independence, and because there's a given amount of jobs available in NZ, our samples are not independent... I'm not sure about it though (you can see his feedback quoted below). ## Way to help me #2 (longer) If you have some spare time, I would appreciate it very much if you could look at the whole assignment. I will also provide the lecturer's feedback, so if you could interpret it for me, it would be great! The assignment is very easy for a mathematician / statistician, there's only two questions there, it's just full of padding where I tried to demonstrate that I know what I'm doing, you can skip most of it. Here's the link to a PDF file with the assignment I didn't succeed in: statistics assignment.pdf. ## Lecturer's feedback > Your figure 1 exhibits sequential correlation which is the real reason why > linear regression does not work. Neither fisher's test nor chi squared is > good for your 2x2 table. This is because you want to test homogeneity, but > you are rejecting the null because of non-independence (which is not > interesting). The distinction between the two is irrelevant here (they are > asymptotically identical in any case). You could have plotted the odds ratio > as a function of time.


title: whether Y(employees injured) variation is due to X1(job function) or X2(population)

text: Here is the actual question- There are 1000 employees in a firm, and the firm has four departments namely D1, D2, D3 and D4 with 100, 200, 300, 400 employees respectively. Now, each employee is explained about what constitutes an injury(for this experiment). Then, each employee is asked whether they were injured in last year. Their responses grouped according to the departments they are in - D1= 10, D2= 20, D3= 30, D4=40, It is to be noted that each department has different job function like for example D1 is HR department, and D4 is machining department. Now, we have to determine whether no. of employees injured in each department (variation) is due to- Job function of the departments OR Population of the department? PS- I'm new to statistics, so this question can be very easy or trivial, it would definitely help me a lot if any could guide me on this. Thanks in advance. * * * This has been my thought process- Y= No. of employees injured last year. X1= Population (of employees) in department X2= Job function of the department both Y and X1 are numerical variables and X2 is categorical. If we find correlation coefficient of Y with X1 and Y with X2, it will be same. So we have to think of other way. Also, I think in this case there is collinearity between X1 and X2, which may cause one of the variable to be redundant or rather both X1 and X2 imply same thing. So, I'm kind of stuck with here. Any help regarding this question would be deeply appreciated.


title: Clarification needed about min/sim hashing + LSH

text: I have a reasonable understanding of the technique to detect similar documents consisting in first computing their minhash signatures (from their shingles, or n-grams), and then use an LSH-based algorithm to cluster them efficiently (i.e. avoid the quadratic complexity which would entail a naive pairwise exhaustive search). What I'm trying to do is to bridge three different algorithms, which I think are all related to this minhash + LSH framework, but in non-obvious ways:   1. The algorithm sketched in Section 3.4.3 of Chapter 3 of the book Mining of Massive Datasets (Rajaraman and Ullman), which seems to be the canonical description of minhashing   2. Ryan Moulton's Simple Simhashing article   3. Charikar's so-called SimHash algorihm, described in this article I find this confusing because what I assume is that although [2] uses the term "simhashing", it's actually doing minhashing in a way similar to [1], but with the crucial difference that a cluster can only be represented by a single signature (even tough multiple hash functions might be involved), while two documents have more chances of being similar with [1], because the signatures can collide in multiple "bands". [3] seems like a different beast altogether, in that the signatures are compared in terms of their Hamming distance, and the LSH technique implies multiple sorting of the signatures, instead of banding them. I also saw (somewhere else) that this last technique can incorporate a notion of weighting, which can be used to put more emphasis on certain document parts, and which seems to lack in [1] and [2]. So at last, my two questions:   * Is there a (satisfying) way in which to bridge those three algorithms?   * Is there a way to import this notion of weighting from [3] into [1]?


title: Reformulation of OLS estimators in a simple Regression with a dummy variable

text: In the Classical Regression Model i.e. $\big(E(y|x)=\alpha +\beta x$ and $Var(y|x)=\sigma^2\big)$ with only two coefficients for intercept $\alpha$ and slope $\beta$ of a dummy variable $x$, we can interpret $\alpha$ as the the mean of values for which $x=0$ and $\beta$ as the difference of the means of the data where $x=1$ and $x=0$ respectively. This makes intutitive sense, but how can I formally show that I can deduce those special expression from the standard definition: $$\hat{\alpha}=\bar y -\hat{\beta} \bar x$$ and $$\hat{\beta}=\frac{\frac{1}{n}\sum (x_i-\bar x)(y_i-\bar y)}{\frac{1}{n}\sum (x_i-\bar x)^2}.$$ I cannot reach the formulation in terms of group means... Thanks a lot, Felix


title: Test whether triplet of numbers is significantly different from $K$ other triplets

text: Given $K$ triplets $t_k:=(a_k^{(1)},a_k^{(2)},a_k^{(3)}) \in \mathbb{R}^3, k=1,...,K$ and the triplet to test against them $t=(a^{(1)},a^{(2)},a^{(3)}) \in \mathbb{R}^3$. How can one test whether $t$ differs significantly from the $K$ other triplets, in the sense that $t$ has "little probability to occur in the distribution of the $t_k$s"? That is, I want to test whether this particular triplet is somehow significantly different from the others. If this weren't triplets, e.g. if I just wanted to compare $a^{(1)}$ against $a_k^{(1)},k=1,...,K$, I would calculate the mean $\mu$ of all $a_k^{(1)},k=1,...,K$ and $a^{(1)}$ and then calculate $\alpha=\#\{k : |\mu - a_k^{(1)}| < |\mu - a_w^{(1)}|\}$. Then I could say $a^{(1)}$ is further apart from the mean as $\frac{\alpha}{K}$ other cases, i.e. reject the hypothesis that $a^{(1)}$ is just randomly distributed around the mean.


title: non-normal residuals. Significant variable

text: The question is simple. My residuals after performing a linear regression are non-normal. I am not sure whether or not one of my variables is significant. I know that I cannot use p-values. Can you suggest a non-gaussian methodology for dealing with this problem? Thanks in advance. Update: The residuals are pretty far away from non-normal (look at the scale to the qqplot on the vertical axis). ![qqplot of my residuals](http://i.imgur.com/8RlYU4c.png?1?2575) Here's a histogram of the candidate variable to be deleted. ![plot of variable](http://i.imgur.com/ZLtxmBU.png?1?6778)


title: Forecast and STL are sensitive to last few data points

text: I'm working on a time series containing hourly data for 8 days. Using the R `forecast` package, I'm trying to predict future data with `stl`. If I use only 7 days data, things work fine, but when I add the last day, the predicted values drop strongly. My questions:   1. Visually the 8 days data still show a clear periodic trend, why forecast gives that result? Looking at the data, I see that the value of the last hour of the day drops to 0 sooner than in previous days. But why forecast are so sensitive with that?   2. Any advice for this situation? Any search keyword or reference that I should know for further reading? Thank you.               xnn <- c(0.11, 0.28, 0.26, 0.27, 0, 0, 0, 0, 0, 0.19, 0.19, 0.19, 0.25,      0.35, 0.42, 0.39, 0.17, 0.19, 0.36, 0.38, 0.41, 0.35, 0.33, 0.37,      0.36, 0.33, 0.06, 0, 0.17, 0, 0, 0, 0, 0.17, 0.31, 0.23, 0.4,      0.38, 0.33, 0.41, 0.34, 0.33, 0.41, 0.33, 0.34, 0.38, 0.41, 0.41,      0.4, 0.38, 0.27, 0.17, 0, 0, 0, 0, 0, 0, 0, 0, 0.44, 0.39, 0.4,      0.4, 0.38, 0.28, 0.36, 0.38, 0.43, 0.42, 0.44, 0.45, 0.41, 0.34,      0.2, 0, 0, 0, 0, 0, 0.15, 0.17, 0.17, 0.21, 0.42, 0.46, 0.42,      0.46, 0.38, 0.34, 0.35, 0.37, 0.38, 0.36, 0.36, 0.38, 0.35, 0.36,      0.39, 0.25, 0.15, 0.15, 0.15, 0.15, 0, 0, 0.28, 0.28, 0.42, 0.43,      0.44, 0.41, 0.41, 0.37, 0.37, 0.31, 0.33, 0.32, 0.34, 0.29, 0.3,      0.28, 0.22, 0.17, 0, 0, 0, 0, 0, 0, 0, 0.17, 0.32, 0.38, 0.39,      0.4, 0.37, 0.41, 0.45, 0.39, 0.35, 0.37, 0.33, 0.38, 0.37, 0.31,      0.15, 0, 0, 0, 0, 0.15, 0.1, 0.15, 0.19, 0.26, 0.38, 0.39, 0.39,      0.38, 0.38, 0.42, 0.41, 0.4, 0.43, 0.41, 0.33, 0.28, 0.35, 0.39,      0.26, 0.31, 0, 0, 0, 0, 0, 0, 0.23, 0.3, 0.4, 0.4, 0.4, 0.43,      0.41, 0.44, 0.43, 0.43, 0.41, 0.31, 0, 0)          require(forecast)     ## hourly data     ## 8 days (24*8 = 192 data points)     dput(xnn)          ## using first 7 days, result is acceptable     x.ts <- ts(xnn[(1+24*0):(7*24)], frequency = 24)     fit <- stl(x.ts, s.window="periodic")     X11()     plot(fit)     X11()     plot(forecast(fit), type = "o")          ## include one more day, unexpected result     x.ts <- ts(xnn[(1+24*0):(8*24)], frequency = 24)     fit <- stl(x.ts, s.window="periodic")     X11()     plot(fit)     X11()     plot(forecast(fit), type = "o")      ![7days](http://i.stack.imgur.com/D3apY.png) ![8days](http://i.stack.imgur.com/oRvSk.png)


title: Fitting data to a Poisson distribution, what are the errors?

text: I have large data set of $\approx 10^6$ points where each point contains the information of a               (year, count)      of a particular event. There are many counts for a each year. Looking at the histograms of counts for a given year, I've conjectured that they follow a Poisson distribution. By generating the Poisson distribution with the maximum likelihood parameter $P(\lambda_{MLE})$, $\lambda_{MLE}=\frac{1}{n}\sum{k_i}$. I _think_ what I need to do next is compare the values of this expected distribution to my normalized histogram and compute a $\chi^2$ statistic. Wiki however warns that: > In cases where the expected value, E, is found to be small (indicating a > small underlying population probability, and/or a small number of > observations), the normal approximation of the multinomial distribution can > fail, and in such cases it is found to be more appropriate to use the > G-test, a likelihood ratio-based test statistic. This will certainly be the case for larger values of $\lambda$ at the more extreme counts. How do I calculate how "good" my data will fit to the Poisson distribution? Ideally what I would like is to show a plot of ($\lambda$ vs year) with an error bar related to this "goodness".


title: How to interpret standard errors in a Cox model?

text: I am running a multi-variate Cox regression and Stata provides the standard errors for each hazard ratio. How are theses to be interpreted? I know that I want my coefficients to be large compared to my SEs, but I don't know if the same rule applies to ratios.


title: Obtaining R pec survival patient risk percentage

text: ### Introduction I have a 300,000-row cancer dataset with around 60 variables (cancer stage, year of diagnosis, radiation therapy, histology, etc.) with a time variable ("number of months survived") and an event (alive or dead). The last two variables have complete values in the individual records. ### Survival months as outcome variable My initial goal was to create a multilayer perceptron model in WEKA given my data in order to _predict the number of months survived_ for new instances.   1. Preprocess data   2. Train the model in WEKA   3. Assess model's performance (accuracy, specificity, sensitivity)   4. Test model for new cancer records ### Patient risk as outcome variable The requirements changed thus it was changed to predicting the _patient risk of survival_ within equally-spaced time periods.   1. Divide data into:      * 24 - 47 months (2 years)     * 48 - 83 months (4 years)     * 84 - 107 months (6 years)     * 108 - 119 months (8 years)     * 120 - "up to what's available" months (10 years)   2. I will then use the function `predictSurvProb` from the package `pec` in `R` as suggested in this problem to obtain individual survival percentages for my aforemention records. The data will be divided into their own survival months bracket and respective patient risk prediction i.e. records that survived within two years will have a patient risk survival prediction percentage in two years.   3. After getting all my individual records their respective survival percentage per time period, WEKA will be used to create five models for each time period that will patient survival as the outcome variable.   4. The five models can be used to predict survival of a single record giving out five different patient survival risk  ### Problem I am still learning about `R` but I managed to apply the sample code (from the `pec` documentation for `predictSurvProb`) into my data as:               library(survival)     library(pec)     library(rms)          # fit a Cox model     coxmodel <- cph(Surv(time,vsr)~1,data=cancer,surv=TRUE)           # predicted survival probabilities can be extracted at selected time-points:     ttt <- quantile(time)          # for selected predictor values:     ndat <- data.frame(vsr=c(0,1)) # I assumed the event variable is provided here          # as follows     predictSurvProb(coxmodel,newdata=ndat,times=ttt) # has error          ## simulate some learning and some validation data     learndat <- SimSurv(100)     valdat <- SimSurv(100)          ## use the learning data to fit a Cox model     fitCox <- coxph(Surv(cancer$time,cancer$vsr)~vsr,data=cancer)          ## suppose we want to predict the survival probabilities for all patients     ## in the validation data at the following time points:     psurv <- predictSurvProb(fitCox,newdata=valdat,times=seq(24,48,72,96,120))     ## This is a matrix with survival probabilities     ## one column for each of the 5 time points     ## one row for each validation set individual      I need to obtain the patient risk calculation for 300,000 patients but the line               predictSurvProb(coxmodel,newdata=ndat,times=ttt)      shows the error               Error in .subset2(x, i, exact = exact) : subscript out of bounds      How do I solve this error?


title: Calculating power function for ANOVA

text: There are several functions in R to calculate the power of a test, for example, the pwr-functions of the `pwr` package. My question is how we get this function for ANOVA (`pwr.anova.test`)? I know that \begin{align} 1-\beta &= P(H1|H1) \\ &= P(F>F_{1-\alpha,k-1,N-k}) \\ &= 1-P(F\le F_{1-\alpha,k-1,N-k}) \\ &= 1-P(MQE/MQR\le F_{1-\alpha,k-1,N-k}) \\ ~ \\ &= 1-P\bigg(n k \bigg(\frac{^{\sum(\bar{x}_{i.}-\bar{x}_{..})^2}/_{(k-1)}} {^{\sum\sum(x_{ij}-\bar{x}_{i.})^2}/_{k(n-1)}} - f^2\!+\!f^2\bigg) \le F_{1-\alpha,k-1,N-k}\bigg) \\ ~ \\ ~ \\ &= 1-P\bigg(nk\bigg(\frac{^{\sum (\bar{x}_{i.}-\bar{x}_{..})^2}/_{(k-1)}} {^{\sum\sum(x_{ij}-\bar{x}_{i.}^2)}/_{k(n-1)}} - \frac{^{\sum(µ_i-µ)^2}/_k}{\sigma^2}\bigg)\le F_{1-\alpha,k-1,N-k}-nkf^2\bigg) \\ ~ \\ &= ???...??? \\ &= 1-P(F_{k-1,N-k}\le F_{1-\alpha,k-1,N-k}-knf^2) \\ &= \end{align}                 p.body =  pf(qf(\alpha, k-1, (n-1)*k, low=F), k-1, (n-1)k, k*n*f^2, low=F)      With $f^2$ the effect;   $knf^2$ the non-centrality parameter; and   $\sigma^2$ the error variance within groups; But what does the function which has a $F_{k-1,N-k}$ distribution look like? It has to be the (variance between groups/df) / (variance within groups/df)$-knf^2$, but how can I write these so that I see the distribution? In my calculation I did not see it. I have never found the exact way how to get to this formula. Maybe you can help.


title: How does one appropriately apply cross-validation in the context of selecting learning parameters for support vector machines?

text: The wonderful libsvm package provides a python interface and a file "easy.py" that automatically searches for learning parameters (cost & gamma) that maximize the accuracy of the classifier. Within a given candidate set of learning parameters, accuracy is operationalized by cross-validation, but I feel like this undermines the purpose of cross-validation. That is, insofar as the learning parameters themselves can be chosen an a manner that might cause an over-fit of the data, I feel like a more appropriate approach would be to apply cross validation at the level of the search itself: perform the search on a training data set and then evaluate the ultimate accuracy of SVM resulting from the finally-chosen learning parameters by evaluation within a separate testing data set. Or am I missing something here?


title: Discrepancy between stepwise and nominal logistic regression results in JMP

text: I have carried out a stepwise logistic regression in JMP. Then (using the proper button in the program window), I have chosen to build a nominal logistic regression model using (only) the variables identified by the stepwise procedure. Anyhow, comparing the summary tables of the stepwise regression and the nominal one, I have recognized that the regression coefficients are not the same, and also the p-values are not the same. There is even a variable which changes from a p-value of 0.02 to a p-value of 0.19 (much greater that 0.10, the threshold value I have chosen before stepwise procedure to retain variables in the model! How is it possible? I could use the values in the stepwise summary, but it does not contains any data allowing to build the confidence intervals. So, in suborder my question is: how can I calculate the confidence intervals using only the data reported in JMP stepwise regression summary? Edit: I have recognized just a minute ago that the differences refer to categorical variables which have yield more than one significant comparison. For example, on stepwise regression details I read variable1 is included in the model three times (and passed three times to the nominal regression procedure): A-B versus C-D-E-F-G, C-D versus E-F-G, E-F versus G. Anyhow, such variable1 is reported only one time in regression summary, which cites only the first comparison (A-B versus C-D-E-F-G). It remains a mistery for me why.


title: Need good material on multifractal analysis

text: I'm searching for some good reading material on multifractal analysis. Preferably something accessible that doesn't put the stress too much on mathematical proofs but rather on applications. As long as it gives a good review of the status of the field, the interesting results and applications, I would be happy. Also, if anything related to multifractal analysis and statistics or time series comes up, I'll take it as well. Books, papers, internetpages, videos, etc... accepted!


title: Evaluating effect sizes of interactions in multiple regression

text: I have been running 3-level multilevel models with HLM, and my main interest is in some cross-level interaction effects that I am finding. My concern is that the effect sizes of these interactions appear to be small – I am wondering whether they are really meaningful. I am turning to you to ask whether anyone could advice me on how to evaluate the size of this kind of effect. Are there any benchmarks for interaction effects in regression analyses? Would these be appropriate for multilevel models, too? I believe that it is not uncommon that the explained variance does not increase much when one adds in interaction terms in multiple regression, and that this would be even more the case for multilevel models when the variance that one is trying to explain is the level 1 variance. Is that right? Thanks in advance!


title: How to utilize the output probabilities from a probabilisitic classifier?

text: Can anybody please explain that in simple words: how can we correctly use probabilities given by some classifiers along with the predicted class values? Let us consider some implementation of a Naive Bayes classifier that works with Gaussian priors. As an output, for each test instance we have a true class label, a predicted class label, and a probability. How can we interpret this probability? Consider, the threshold is 0.5. If probability for some test instance is 0.75, it belongs to the positive class. What if the correct label for that instance is a negative class, i.e. classifier makes an error but give a high probability? Some interpret that probability as "confidence" in a non-statistical sense. Some use these probabilities as parameters of a function they try to optimize. That confuses me. For example, one owns a designer boutique and re-sells items bought cheap online. He wants to minimize costs associated with items that are not truly described (like a sweater have been worn a lot, but description says "like new"). Assume, our classifier is built on features extracted from description, seller's score, feedback, etc. We try to predict, if the item sold on ebay is truly described. If the classifier assigns it to a positive class (item is as described), one can place a higher bid. The maximum bid price one wants to pay for a truly new sweater is multiplied by the probability the sweater is as described. This is a naive model, of course. My observation is that the amount of incorrectly classified instances with high probabilities can be quite high, so it does not work for any kind of optimization. But the example I gave seems to be a popular use of probabilistic classifier. And that probability cannot be interpreted as common sense confidence, is that right?


title: How do you code three summation functions in WinBUGS?

text: I recently learned how to use winbugs. I know nothing about coding and programming. My model is not included in the set of distributions that are already available in winbugs. This is a trivariate poisson model with covariance structure: ![http://i.stack.imgur.com/HS4u7.png](http://i.stack.imgur.com/6LbEA.png) There is this `dloglike` function in winbugs to create a custom model. You just need to get the log transformation and you're good to go. Getting the log transformation is the easy part. What makes me crazy is the presence of those three summation functions. I don't know how to convert them into the Bugs language. I would like to show you my Winbugs code but it is kind of messy. Can you please help me?


title: What is the appropriate method to compare two means if data is nonnormal, sample size is different and variance seems equivalent?

text: I want to compare the means between groups `1` vs. `2` and also `2` vs. `3`. All three groups are nonnormal, as the Shapiro-Wilk test resulted in p values of 0.03, 0.04 and <0.0001 respectively. _n_ is 10, 13 and 38 respectively, and the Levene test for equal variances resulted in p values of 0.22 for `1` vs. `2` and 0.44 for `2` vs. `3` (so they seem to have equal variances). Both, the parametric _t_ -test and the nonparametric Wilcoxon test resulted in insignificant differences with p values larger than 0.5. Which test is more appropriate to this situation? The groups do not follow a normal distribution, have different sample size, but seem to be homoscedastic. Your help will be very much appreciated.


title: Interpreting the output of a neural network

text: I've implemented a neural network for prediction, and for the input data, I've used the following formula to normalize data:               Data_normalized_i= [Data_i - Min_data]/[Max_Data- Min_data]      I've some questions:   1. How do I interpret the output of my network according to my inputs?   2. Must I use the real data input to compare it with my outputs?   3. If I have to do some transformation of my outputs, how should I go about doing it? For the test error in this case, will be it calculated from the output or from the transformed outputs?


title: Recommend a graduate-level regression textbook?

text: Some friends and I (statistics grad students) would like to work through a good regression book. We're looking for   * single and multiple linear regression   * logistic regression   * model selection   * diagnostics   * splines   * penalized methods   * classifiers, trees Can anyone recommend a book (or books) that cover these topics?


title: Confidence intervals for count data of a categorical

text: I'm trying to analyze the results of a survey question. The question is multiple choice (but each respondent can only select one answer) and I want to calculate 95% confidence intervals for their responses.               A      122     B       55     C       16     D       14     E       13     ----------     Total  220      If I'm using R, can I treat each choice as a binomial (122 out of 220 chose _A_ , 98 did not choose _A_ ) and use `binconf()` to compute the interval or is there a better method for categorical data?


title: mobForest R Package

text: I have recently begún to learn about model based recursive partitioning by playing around with MOB in the party package. I came across this mobForest package but am a little baffled towards what it is actually doing. A MOB model (if I can call it that) is of the form Y ~ X1,...,Xk | Z1,...,Zl using the MOB vignette's notation. In mobForest I would have a sequence of partitions (Z1,...,Zl ;...; A1,...Al) for which evey node would have a model of the form Y ~ X1,...Xk fitted to it. Is this an ensemble of MOB models or am I completely misunderstanding this? Bootstraping the observations and sampling your variables at splits to build a sequence of trees helps tackle many problems that a single tree has. But how does that come into play here? On a side note, when I use mobForestAnalysis to do model-based random forest analysis I dont know how to use this (Model?) to predict on new observations. Can someone help me out on this?


title: Conducting correlation and one way ANOVA using data from a PostgreSQL database

text: I am trying to discern the best way to calculate a correlation and perform a one-way ANOVA on data I am taking from a PostgreSQL database.   * What tools should I use?    * Can I do this using the SQL language itself?    * Is there an easy way to export the data?


title: Fitting the moving average model

text: let v to be forecasted value for periods 1 through T and $v_{t}$ be its forecasted value at time $t$. We express $v_{t}$ as the sum of two terms, its mean at time $t$, and its deviation from the mean at time $t$, $\epsilon_{t}$. In other words, $$ v_{t}= \overline{v_{t}} + \epsilon_{t} $$ The $\overline{v_{t}}$ are chosen based on the arguments. The $\epsilon_{t}$ term is assumed to be a normally distributed random variable with mean zero and standard deviation $\sigma (ε_{t})=0.234$. The moving average formation of order q is choosen, MA(q) where q is the number of lagged terms in the moving average. We use the following moving average specification: $$\epsilon_{t} = \sum^{q}_{i=0}{\alpha_{i} \mu_{t-i}} $$ where $\mu_{t-i}$ are independently distributed standard normal random variables. To ensure that the standard deviation of $εt$ is equal to its pre- specified value, we set the $$\alpha_{i}= \frac{\sigma(\epsilon_{t})}{\sqrt{q+1}}$$ Note that $\epsilon t$ depends on $q+1$ random terms. the R-code that i have used for above model               q=31     iter=10000     for(i in 1:31){alpha_i[i] <- 0.234/(sqrt(31+1))}     err_pk <- array(0,c(iter,11))     for(i in 1:iter){err_pk[i,] <- arima.sim(list(order=c(0,0,31), ma=alpha_i), n=11, innov=rnorm(iter,0,1))}     ffe <- array(0,11)     for(i in 1:11){ffe[i] <- median(err_pk[,i])}`     plot.ts(ffe)      I am wondering that, $\alpha$ is changing through time? and also i never get the same output as in the figure 2 The end of world Population growth Nature, v412, 543 the parameter for figure in the paper are: Note: MA(30), (31 terms), $\sigma(\epsilon_{t})=0.234$, 31 initial values of $\mu=0$, 10,000 simulation Am I missing any thing?


title: R glm probit regression marginal effects

text: I am using R to replicate a study and obtain mostly the same results the author reported. At one point, however, I calculate marginal effects that seem to be unrealistically small. I would greatly appreciate if you could have a look at my reasoning and the code below and see if I am mistaken at one point or another. My sample contains 24535 observations, the dependent variable `x028bin` is a binary variable taking on the values 0 and 1, and there are furthermore 10 explaining variables. Nine of those independent variables have numeric levels, the independent variable `f025grouped` is a factor consisting of different religious denominations. I would like to run a probit regression including dummies for religious denomination and then compute marginal effects. In order to do so, I first eliminate missing values and use cross-tabs between the dependent and independent variables to verify that there are no small or 0 cells. Then I run the probit model which works fine and I also obtain reasonable results:               probit4AKIE <- glm(x028bin ~ x003 + x003squ + x025secv2 + x025terv2 + x007bin + x04chief + x011rec + a009bin + x045mod + c001bin + f025grouped, family=binomial(link="probit"), data=wvshm5red2delna, na.action=na.pass)          summary(probit4AKIE)      However, when calculating marginal effects with all variables at their means from the probit coefficients and a scale factor, the marginal effects I obtain are much too small (e.g. 2.6042e-78). The code looks like this:               ttt <- cbind(wvshm5red2delna$x003,         wvshm5red2delna$x003squ, wvshm5red2delna$x025secv2, wvshm5red2delna$x025terv2,     wvshm5red2delna$x007bin, wvshm5red2delna$x04chief, wvshm5red2delna$x011rec,         wvshm5red2delna$a009bin, wvshm5red2delna$x045mod, wvshm5red2delna$c001bin,     wvshm5red2delna$f025grouped, wvshm5red2delna$f025grouped,wvshm5red2delna$f025grouped,         wvshm5red2delna$f025grouped,wvshm5red2delna$f025grouped,wvshm5red2delna$f025grouped,     wvshm5red2delna$f025grouped,wvshm5red2delna$f025grouped,     wvshm5red2delna$f025grouped) #I put variable "f025grouped" 9 times because this variable consists of 9 levels          ttt <- as.data.frame(ttt)          xbar <- as.matrix(mean(cbind(1,ttt[1:19]))) #1:19 position of variables in dataframe ttt          betaprobit4AKIE <- probit4AKIE$coefficients          zxbar <- t(xbar) %*% betaprobit4AKIE          scalefactor <- dnorm(zxbar)          marginprobit4AKIE <- scalefactor * betaprobit4AKIE[2:20]           #(2:20 are the positions of variables in the output of the probit model 'probit4AKIE'      #(variables need to be in the same ordering as in data.frame ttt), the constant in        #the model occupies the first position)          marginprobit4AKIE #in this step I obtain values that are much too small


title: How to handle categorical dependent variable using logistic regression when one of the categories never occurs in the sample

text: I have a multinomial logistic regression model. One of the output categories is not observed in the data set that I'm using. ### Example:   * 4 different diagnoses (response variable) in the population, but in the sample, Type 3 never occurred   * 5 hormone level measurements (predictors)  ### Question   * What books/papers discuss the mathematics of handling this situation in logistic regression?


title: Difference between K-S manual test and K-S test with R?

text: On the Internet there is an example of k-s test being applied relative to distribution of number of bird varieties over different five hour periods. The observed distribution was:               a=c(0,1,1,9,4)      The expected distribution (if there is no difference between the five hours) could be:               b=c(3,3,3,3,3)      After I found the two cumulate distributions, I calculated manually, D = 0,4667 (a value that is similar with the internet). But if I try to use R, I find a different value of D:               > a=c(0,1,1,9,4)          > b=c(3,3,3,3,3)          > ks.test(a,b)          Two-sample Kolmogorov-Smirnov test          data: a and b           D = 0.6, p-value = 0.3291          alternative hypothesis: two-sided .......        * What is leading to the difference between my manual calculation and result that R gives?


title: What is this measure of error?

text: I'm reading through someone else's code for plotting the results of a psychology experiment, and (according to the code comments) they calculate the accuracy error of their behavioral paradigm as follows: $\textit{accuracy error} = \sqrt{\frac{(\textit{accuracy}) (1-\textit{accuracy})}{\textit{total trials}}}$ It's output seems to be very similar to the original accuracy. What is this? Is this some sort of multiple comparison correction? Why would they do this?


title: Bayesian estimation of Dirichlet distribution parameters

text: I want to estimate parameters of Dirichlet mixture models using Gibbs sampling and I have some questions about that:   1. Is a mixture of Dirichlet distributions equivalent to a Dirichlet process? What is their main differences if is not?   2. Also, if I want to estimate a single Dirichlet distribution's parameters, which distribution for parameters should be selected as priors in Bayesian framework? In all of the papers I found an estimation of a multinomial distribution using Dirichlet priors. I need estimation of a Dirichlet distribution using multinomial priors, perhaps. Is the posterior function also in the form of DIRICHLET(α+N) similar to the case “estimation of multinomial distribution using Dirichlet priors”? as the multiplication of probability density function for iid samples are not considered in the definition of the likelihood function. I again cannot understand why. e.g. as expressed in this paper: http://www.stat.ufl.edu/~aa/cda/bayes.pdf or http://research.microsoft.com/en-us/um/people/minka/papers/minka- multinomial.pdf * * * so thanks for your attention my data is Hyperion (a kind of hyperspectral remote sensing imagery) and i want to perform hyperspectral unmixing using mixture of Dirichlet sources and i will apply Gibbs sampling method for parameter estimation. my data is in dimension (614*512*224) which is commonly available AVIRIS sensor data for Cuprite Nevada district and is almost 200MB. also this data is available via (http://aviris.jpl.nasa.gov/data/free_data.html). unfortunately i don't know how can i sent my data. i just ask you to help me in statistical modelling tasks for my PHD thesis. i will be so grateful if you help me to solve my confusions in modelling. all the best solmaz


title: How to store checks of gradient algorithm in a matrix using R?

text: I would like to check different gradient algorithms. For example:                fr <- function(x) {   ## Rosenbrock Banana function              x1 <- x[1]              x2 <- x[2]             print(c(x1,x2))              100 * (x2 - x1 * x1)^2 + (1 - x1)^2          }     optim(c(-1.2,1),fr,method="BFGS")      prints to the screen the values at which the RBF has been evaluated. How can I store these values in a matrix ? (instead of just printing them to the screen)


title: Does correlation assume stationarity of data?

text: Inter-market analysis is a method of modeling market behavior by means of finding relationships between different markets. Often times, a correlation is computed between two markets, say S&P 500 and 30-Year US treasuries. These computations are more often than not based on price data, which is obvious to everyone that it does not fit the definition of stationary time series. Possible solutions aside (using returns instead), is the computation of correlation whose data is non-stationary even a valid statistical calculation? Would you say that such a correlation calculation is somewhat unreliable, or just plain nonsense?


title: How to optimize the k parameters in dynamic linear regression?

text: I am starting to use R's dynlm package. Currently I am just looking at the fit and eyeball which choice of lags might be the best. Is there a standard way or a strategy to determine the best k parameters for `L()`. What I often see is ridiculously high lags like k=10 in a quarterly series delivering the best fit. What could be the reason for that?


title: Different odds ratios from formula and R's fisher.test

text: In the following example               > m = matrix(c(3, 6, 5, 6), nrow=2)     > m          [,1] [,2]     [1,]    3    5     [2,]    6    6     > (OR = (3/6)/(5/6))    #1     [1] 0.6     > fisher.test(m)        #2              Fisher's Exact Test for Count Data          data:  m      p-value = 0.6699     alternative hypothesis: true odds ratio is not equal to 1      95 percent confidence interval:      0.06390055 5.07793271      sample estimates:     odds ratio       0.6155891       I calculated the odds ratio (#1) "manually", 0.600; then (#2) as one of the outputs of the Fisher's exact test, 0.616. Why didn't I get the same value?


title: Finding the Bayes estimator of $\theta$ - having trouble with likelihood calculation

text: Let $Y_n$ be the nth order statistic of a random sample of size n from a distribution with pdf $f(x|\theta) = 1/\theta$, $0<x<\theta$, zero elsewhere. Take the loss function to be $L[\theta, \delta(y)] = [\theta - \delta(y_n)]^2$. Let $\theta$ be an observed value of the random variable $\Theta$, which ahs the prior pdf $h(\theta) = \beta\alpha^{\beta}/\theta^{\beta+1}$, $\alpha < \theta < \infty$, zero elsewhere, with $\alpha > 0$, $\beta > 0$. Find hte Bayes solution $\delta(y_n)$ for a point estimate of $\theta$. I'm a little confused with finding the likelihood function. Since they are telling us that $Y_n$ is the nth statistic...does that mean that we only have one function in the likelihood function? Thanks in advance


title: What do normal residuals mean and what does this tell me about my data?

text: Pretty basic question: What does a normal distribution of residuals from a linear regression mean? In terms of, how does this reflect on my original data from the regression? I'm totally stumped, thanks guys


title: Computing percentile rank in R

text: How can I add new variable into data frame which will be percentile rank of one of the variables? I can do this in Excel easily, but I really want to do that in R. Thanks


title: Given weighted mean and raw data, derive weights to minimize error

text: Imagine I have a set of data from an experiment:   Observed Measurement 1: 5   Observed Measurement 2: 6   Observed Measurement 3: 7   Observed Measurement 4: 8   Observed Measurement 5: 9 I have a target mean of 6--I need to choose a set of weights that map the observational data to the target mean. What's the general process for deriving the weights of each measurement in the desired weighted mean? Obviously for one set of data there are infinitely many weight possibilities (particularly if there are no bounds for he possible weights), but say I had 10 sets of data, each of which used the same weights ---Weight 1 is always used on the first measurement, weight 2 is always used on the second measurement, etc-and say I wanted to minimize the use of weighting by getting each weight as close to 1.0 as possible. I have a large-n number of measurement,mean sets, well over 100, and I want to use the same weight for each and get the closest results possible. It's been a number of years since I've done any kind of statistics and I'm thinking there has to be some sort of line of best fit, least-squares, or something? If it's a process that I can automate into a little python script, more's the better. My day job mostly involves qualitative stuff these days, not quantitative, but I do have some background in math.


title: Appropriate way to estimate sample size requirements for an unknown population distribution

text: Assuming we have a sample (some observables) from a population of unknown distribution, can I assume that the right way to estimate the sample size (for estimating means, medians, etc.) is using the formula: $$n = 4s^2 / d^2$$ where $s$ is the sample standard deviation, $d$ is the margin of error (precision), and $n$ is the sample size? ... or is there some implicit assumption that the population itself must be Normal in order for this to be valid?


Which variety of Linux is best suited to work with R, SAS, SQL, Matlab? Is it Scientific Linux?


title: how to forecast daily sale using Excel

text: I am a newbie here. well I am trying to find a method or a formula to forecast meals per day, which have ̀5 meals to upload on flights, sales, wastage and passengers are what I have to consider, the old template is not completed yet , and its not so good to forecast, and I can't think of other formulas or methods to forecast , I have the sales in the past few months.. anyone can suggest me which methods can solve this problem? well I am using MS.Excel to calculate or if there is a program to suggest me it would be great :) ..thanks a lot for helping and sorry for my bad english :)


title: Intensive resources for learning R (e.g., classes, workshops)?

text: I have read the wonderful thread entitled Resources for learning R. The question asks for resources to learn R "on the cheap". I have a similar question: What is the best and quickest way to learn R, regardless of cost? I am thinking of a week-long intensive course, or similar (could even be much longer). The course would involve learning R from the very beginning and be designed for people with a background in statistics.