title: How to report log linear models of contigency tables

text: I am using log linear models (`loglm` function, library `MASS` of `R`) to evaluate if 3 variables in a 3 way contingency table are independent. I build the model of mutual independence               loglm(formula = ~A + B + C, data = test.t)      Which gives me               Statistics:                           X^2 df P(> X^2)     Likelihood Ratio 264.7872 50        0     Pearson          292.6937 50        0      From what I understand this is the LR test compares my model to the saturated model and observes that there is unexplained variance in my model and significant interactions need to be incorporated, which means I can reject my hypothesis that the 3 variables are independent. How exactly should I report them this analysis in my report? Do I need to state the LR test values, degrees of freedom and P? Is this the Pearson Chi square test in the second line? I was under the impression that the Pearson chi square test is only for 2x2 tables (the `chisq.test()` throws an error in bigger tables). Or is it the Pearson chi squared for the 2 models (my model vs the saturated model)?


title: Correlation between series (computing optimization)

text: **Objective:** Get the red line to go as high as possible. I can control the values of the green and orange lines. There are like 50 other lines that I control and can pit against the red line as well. What is a good way to determine which lines are correlated with movement in the red line? For example:   * Which lines tend to go up when the red line goes up?   * Which lines tend to go down when the red line goes up?   * Which lines seem to have an "optimal value" for the red line (anything higher or lower is correlated with the red line being not as high) Is there some kind of mathematical analysis that can be run against the series to get answers to these questions? For any programmers, am using JavaScript and Python, so any libraries in those languages would be helpful as well. **More background info** This question applies to a horizontally scaled super computing system. The red line is how much the system is getting done per minute. The green line is the database response time for a certain operation. So if the database response time that the green line represents is low when the red line is high and high when the red line is low, then I know that I need to look in to optimizing the response time of that database. The orange line represents how many threads are running in the system. Obviously, too few threads mean less will get done, but I also suspect that too many threads may bog the system down as well, so it would be cool to be able to know if there is an "optimal" value for the number of threads running in the system. Thanks so much! ![enter image description here](http://i.stack.imgur.com/71YQk.jpg)


title: Mixed model in simple english

text: Can someone please explain the intuition behind mixed models in a nutshell? Whenever I read explanations, I get overwhelmed by notation and mathematical jargon. Can someone give me a simple example or motivation to help me understand the concept in a deeper way?


title: Are Restricted Boltzmann Machines better than Stacked Auto encoders and why?

text: So I'm learning about deep learning. I first learned about stacked auto- encoders and now I'm learning about Restricted Boltzmann Machines. However non in the papers/tutorials I read I found them motivating why would one want to use RBM instead of auto-encoders. So what are the advantages of RBM over stacked auto-encoders? And when should one use RBM or auto-encoders?


title: detect line in geocoordinates

text: I have repeated samples of geocoordinates of activities in a city. In most of these samples positions will simply be random. In some samples, however, some percentage of the data will be arranged -- with some error -- on a line. Let me paint you a picture (code to generate this data below): ![Scatter with line](http://i.stack.imgur.com/Z6hPy.png) I want to be able to   * detect whether or not the sample contains a line   * estimate the two parameters of this line from the data if it does I've toyed with doing this with a Hough transform (see my previous question) but am struggling to make that work in a robust way. I have therefore tried to take another stab at the problem using a Bayesian mixture model. The code below generates model data and implements such a mixture model in PyMC, which works surprisingly well (note: most of my data will have `p=0`, the rest will have strictly positive `p`, anywhere from `0.2` to `0.5` I would guess):               import numpy as np     from matplotlib import pyplot as plt     import pymc as pm     import scipy.stats as stats          # DGP parameters     N = 300     p0 = 0.2     alpha0 = 0.3     beta0 = 0.2     sigma0 = 0.01     sigma_noise = 0.2          # DGP     x = np.random.normal(0.5,0.2,N)     y_line = np.random.normal(alpha0+beta0*x,sigma0)     line_bool = np.array(stats.bernoulli.rvs(p0, size=N))     y_rand = np.random.normal(0.5,0.2,N)     y_all = y_line*line_bool + y_rand*(1-line_bool)          # plot     fig = plt.figure(figsize=(8,8))     ax1 = fig.add_subplot(111,aspect='equal')     ax1.scatter(x[line_bool==0], y_all[line_bool==0],c='b')     ax1.scatter(x[line_bool==1], y_all[line_bool==1],c='r')          # prior for the assignment probability is a beta distribution that puts a lot of probability mass near 0     p = pm.Beta("p", 1, 3, value=0.5)     assignment = pm.Bernoulli("assignment", p, size=N)          # priors for the noise component     center_noise = pm.Normal("center_noise", 0.5, 0.4)     tau_noise = 1.0 / pm.Uniform("tau_noise", 0, 1) ** 2          # priors for the line parameters     beta_min = -10**1     beta_max = 10**1     line_pars = pm.Uniform("line_pars", beta_min, beta_max, size=2)     tau_line = 1.0 / pm.Uniform("tau_line", 0, 0.05) ** 2          # deterministic functions for the means and variances conditional on assignment     @pm.deterministic     def center_i(assignment=assignment, center_noise=center_noise, x=x, line_pars=line_pars):         center_line = line_pars[0]+line_pars[1]*x         return assignment*center_line + (1-assignment)*center_noise          @pm.deterministic     def tau_i(assignment=assignment, tau_noise=tau_noise, tau_line=tau_line):         return assignment*tau_line + (1-assignment)*tau_noise          # define the observed variables     x_obs = pm.Normal("x_obs", 0.5, 0.2, value=x, observed=True)     y_obs = pm.Normal("y_obs", center_i, tau_i, value=y_all, observed=True)          # define the model     model = pm.Model([p, assignment, center_noise, tau_noise, line_pars, tau_line])          # sample from the posterior     mcmc = pm.MCMC(model)     map_ = pm.MAP( model )     map_.fit()     mcmc.sample(50000, 20000, 3)          # look at the posterior     from pymc.Matplot import plot as mcplot     mcplot(mcmc.trace("p", 2), common_scale=False)     mcplot(mcmc.trace("line_pars", 2), common_scale=False)     mcplot(mcmc.trace("tau_noise", 2), common_scale=False)     mcplot(mcmc.trace("tau_line", 2), common_scale=False)      This works well on generated data (though I would welcome suggestions for how to improve it!), but I suspect the real data this will have to work on will be a bit messier. In particular the "background noise" (`y_rand` above) isn't nicely normal. I'm therefore wondering how I would generalize the model above:   * How would I define mixture models in which the mixture concerns not only the parameters of the distribution but the type of distribution (say, the background noise being uniformly distributed)?   * How might I make the distribution of the background noise far more general. My data being geocoordinates of activities in a city, is there way a way of taking into account that even in the course of normal activity (with `p=0`) positions may more often be arranged along a line than the simple normal model above would suggest (e.g. because they take place along a straight road)? Is there a way of "diffing out" the kind of activity that usually takes place on a Wednesday afternoon when estimating the model on data from a new Wednesday afternoon?


title: What is the distribution of edge lengths for a Euclidean minimal spanning tree?

text: Suppose that you have a Poisson point process in $\mathbb{R}^d$ with rate (density) parameter $\rho$; I particularly care about $d=2$ if there is a solution in this special case. Given a finite realization, one can construct the minimal spanning tree of the ensemble. I'd like to know (model/approximate) the distribution of the edge lengths, i.e. if I select an edge from the tree at random, what is the distribution for it's length $l$? I'm pretty sure by scaling arguments that $p_d(l \vert \rho)=f_d( \rho^{-1/d} l )$, but I don't have a good handle on how to get at the functional form.


title: Detecting abnormal points in point cloud

text: I've asked the same question at Math SE, but the suggestion is that probably this question belongs here. Given a list of point cloud in terms of $(x,y,z)$ how to determine abnormal points? The motivation is this. We need to reconstruct a terrain surface out from those point cloud, which the surveyors obtain when doing field survey. The surveyors would take an equipment and record a sufficient sample of the $x,y,z$ of a terrain. Those points will be recorded into a CAD program. The problem is that the CAD file can be corrupted from time to time with the introduction of "abnormal" points. Those points do not fit into the terrain surface generally, and tend to have erroneous $z$ value ( i.e., the $z$ value is outside of the normal range). I am aware that the definition of abnormal points is a bit loose; and I can't come up with a rigorous definition of it. However, I know what is an abnormal point when I see the drawing. Given all these constraint, is there any algorithm to detect these kinds of abnormal points?


title: Probability statement with a confidence interval

text: Suppose you construct a 95% confidence interval on a mean. Can you state that there is a 97.5% probability that the true population mean lies below the CI upper bound and that there is a 97.5% probability that the true population mean lies above the CI lower bound?


title: Finding correlation of attributes based on known decisions

text: If there are two matrix of data, rows are all a person (under 50 users). In the first dataset, columns are attributes (over 500 columns) of a user with the rating 0 for does not exist, 1 for exist. Such as 1 for married, 0 for not. The second dataset is a choice of the user. 0 is no data, -1 is does not like, 1 is likes. The user choice matrix has some missing information, the user attribute is complete. From the known user decisions, how to rank the importance of the user attributes for their decisions. I have looked into collaborative filtering, but the problem is not filtering things now, it is finding which attributes are most telling of the current dataset. My first ideas were doing a correlation of the likes/dislikes. That does show correlations. However, I can not understand what analysis can be done to rank the attributes of a user. My second thought was to use kmeans clustering for the users based on decisions, but still, it seems that this does not really account for the importance of specific attributes in correlation to the user decisions. Is this a known problem with a method of solving, or does this need to have some custom method, or is it impossible analytically?


title: I want to study the effect of the normality assumption on intervals like t-confidence interval and $\chi^2$ interval.. by using R

text: I want to study the effect of the normality assumption on intervals like t-confidence interval and $\chi^2$ interval.. (1) Generate 1000 random samples from a normal population $N$($\mu$,$\sigma$). From each of these samples, construct 95% confidence intervals for the population mean and for the population variance. Verify that out of the 1000 intervals you create, about 950 of those would actually capture the true parameter. Could someone get me started on this? (2) Design a simulation study using R to analyze the accuracy (95% capture rate) of the t- confidence interval and the $\chi^2$ interval when the population is no longer normally distributed. You should try at least five different populations (e.g. gamma, F, t, etc). Also you should experiment with different sample sizes (from small to large samples, i.e. n = 10; 20; 50; 100) to assess how the sample size affect the capture rate. Could you show me one example by using R for examining one population such as gamma ?


title: Interaction term is insignificant but the interaction plot shows signs of signficance

text: So, I ran an ANOVA on a dataset which contains over 4500 data points. Without getting into too much detail, I was wondering if there is a simple explanation to this: The p-value for my interaction term A*B is 0.145 => The interaction term is clearly insignificant. Both A and B are categorical variables with 2 levels each. However, when I plot my interaction plot, there's a clear X (i.e.) it is not parallel as one would expect. Is there something I am missing here? here's the plot- ![enter image description here](http://i.stack.imgur.com/H1qc7.jpg)


title: Values of the weights in Adaboost

text: I have implemented a simple Adaboost algorithm, using several weak classifiers, and when checking the values computed by it there are alphas with a negative value. Is that possible, or is there a bug in the implementation? Thank you.


title: How do instrumental variables address selection bias?

text: I'm wondering how an instrumental variable addresses selection bias in regression. Here's the example I'm chewing on: In _Mostly Harmless Econometrics_ , the authors discuss an IV regression relating to military service and earnings later in life. The question is, "Does serving in the military increase or decrease future earnings?" They investigate this question in the context of the Vietnam war. I understand that military service cannot be randomly assigned, and that this is a problem for causal inference. To address this issue, the researcher uses draft eligibility (as in "your draft number is called") as an instrument for actual military service. That makes sense: the Vietnam draft _randomly_ assigned young American men to the military (in theory--whether the draftees actually served touches on my question). Our other IV condition seems solid: draft eligibility and actual military service are strongly, positively correlated. Here's my question. It seems like you'd get self-selection bias: maybe richer kids can get out of serving in Vietnam, even if their draft numbers are called. (If that wasn't actually the case, let's pretend for the sake of my question). If this self-selection creates systemic bias within our sample, how does our instrumental variable address this bias? Must we narrow our scope of inference to "the types of people who couldn't escape the draft?" Or does the IV somehow salvage that part of our inference? If anybody could explain how this works, I'd be very grateful.


title: Converting a parametric survival model to a cash flow model. How do I account for aging in the population?

text: I'm building a survival model for time to failure of widgets. Other members of the team want to convert the model to a cost flow model. The basic idea is that we can use the functional form of the survival function under the parametric approach (where we assume a logistic distribution) to estimate the cohort level survival curve. Then they can calculate that if 10% of the population has failed after 1000 days, we will incur $100,000 in replacement costs. The estimation population is all widgets which have survived to 1000 days; that is the only criteria by which you are chosen as part of the population. We have the dates of all the failures, and all the widgets that have survived to the present day are considered censored. However, when we do the prediction, we will be predicting time to event for any widget that is 1000+ days old. And this is the part I'm struggling to account for. If a widget has survived for 100 additional days, this information should be included in the model somehow, but I can't figure out what rigorous approach would be. I thought I might be able to use the conditional hazard function, and exploit the relationship between the hazard and the survival function by taking the integral in a piecewise fashion. $S(t)$ is 1 for all days between 0 and 100 and then the $\exp(\int_{100}^th(t))$ thereafter, but that feels very hackerish and isn't getting me the results I expected (that might also be because the model is crap).


title: Unbalanced data or Balanced data

text: Unbalanced panels are more common in economic fields, if I want to know the behaviour of firms, what will be the differences using unbalanced data panel. Are there advantages? Does it depend on analysis's period? Or it will be better use balanced panel? Thanks!


title: What difference (if any) exists between the Response Distribution and Error Distribution in GLMs?

text: Ok, forgive my ignorance, but I keep getting confused about something at the core of GLMs. Some textbooks describe the two main parts of a GLM as the link function and the distribution of the error terms. Others describe the two main parts as the link function and the variance function, where the variance function is a description of the relationship between the mean and the variance of the response (i.e., the response distribution). But the error distribution and the response distribution seem like different things to me. If I had an equation where $Y_i = B_0 + X_i*B_1 + e$, I can see how for any given value of $X_i$ (plus constant values of $B_0$ and $B_1$), each random variable $Y_i$ would take on whatever distribution the error term had. But does that necessarily make the overall response distribution equal to that same error distribution? Do my questions even make sense?


title: How to preform non-parametric bayesian based regression (predictions) in R?

text: I am working on some **non-parametric bayesian based predictive analysis** using **R**. I have a set of data which denotes various parameters of an online transaction. Based on these parameters I want to develop a model which will provide predictions for future online transactions. The training data consist of records in this format:               transaction_id (numeric)| duration (integer)| amount | is_holiday (boolean) | status(1 or 0)                             |                   |        |                      |                             |                   |        |                      |       The problem that I am facing is that I do not know how to proceed ahead. I am do know know what are the steps that I need to follow. I looked up and found that there are few packages in R like `DPpackage` which have some functions for non-parametric bayesian modeling but there is no concrete example available about how to use it in order to perform various steps of training and testing. It would be helpful for me if someone could provide me some guidance as in which process will be better for such kind of predictive/regression analysis and how to proceed ahead, like what steps should I perform to get the training and testing done. Thanks in advance!


title: Is there a reliable recursive formula for a simple moving average (moving mean)?

text: I've tried some recursive moving average formulae (to reuse a previous output instead of summing the whole `n`-long set for every `i`) I've managed to find but none of them produces the same results as a bare moving mean does. Is there a reliable recursive formula which would produce exactly (or almost exactly) the same output as a bare moving mean?


title: How to get the actual regression value from regression analysis in R statistical package

text: this is, I think, an easy question... I did a regression analysis in R, where I wanted to check the fit of my data to a specific formula I provided... I got this to work, I see the graph and the line showing the fit of the data, but I also would like to get actual regression value obtained from this calculation, how do I do this? and How do I get the actual equation printed in the graph? * * * Hi thanks a lot for the replies. The example was very useful, but I still didnt obtain the Regression value, so I still dont know how well my data fits the equation... below Im pasting the code I used               # import the input data file "LD251-chilR.txt", for the input;      # the first column is pairwise distance and the second is LD estimate,      # need to put the file in the same folder with R program     CT251chil<-read.table("LD251-chilR.txt", sep="\t", dec=",", header=TRUE)           # your LD estimate (r2 in this case which located in column 2)     r2=CT251chil[,2]           # run nls function for getting rho estimate     nls(r2~1/(1+rho*CT251chil[,1]), start=list(rho=0.3))           # getting this rho estimate after running nls function in R     rho=0.02206872           #sort the data, no need if your data already sorted     dist<-sort(CT251chil[,1])           # put parameters in the equation as shown in my MBE paper     eq <- (((10+rho*dist)/((2+rho*dist)*(11+rho*dist)))*(1+((3+rho*dist)*(12+12*rho*dist+(rho*dist)^2)/(46*(2+rho*dist)*(11+rho*dist)))))           # plotting the graph between LD estimate and distance     plot(CT251chil[,1],CT251chil[,2], col="black", pch=20,           ylab=expression(R^2), xlab="Pairwise distance",           main="CT251-chilense", las=1)           # getting regression line     lines(dist,eq, col="black",lwd=2,lty=1)       After this point I get a graph with a line, all good. But so far I dont know how to output the actual equation on the graph get the value of R. thank you for any further help!


title: Nonparametric bootstrap confidence interval

text: I have generated a 95% confidence interval $(22,25)$ for a parameter using a nonparametric bootstrapping method. What I want to know is what value we use for the estimate of the parameter. Is it the mean of all the bootstraps, the median of all the bootstraps, the original estimate before we bootstrapped from the original dataset, or something else?


title: How do you interpret the parameters obtained from lasso logistic regression when the response is binary?

text: Are we still able to interpret the parameters in the same manner as we would in ordinary logistic regression? I'm asking this because I'm toying with the german credit data https://archive.ics.uci.edu/ml/datasets/Statlog+(German+Credit+Data) and doing a lasso logistic regression model to perform feature selection. I'm using the good_bad credit variable as the response and I'm using this model with the intent of obtaining probabilities for defaulting on a loan.


title: Multivariate model and large regression

text: I am not familiar with the concept of multivariate model and just learning about regression model. I am familiar with Autoregressive model and Moving Average. Multivariate regression model provided cues to look into what a multivariate model is but I could not follow the question as such in this thread. My problem is in the area of sparse regression which I understood as where the coefficient matrix has more number of zero entries.   1. My first question is, the coefficient array should be a vector or why is it a matrix?   2. Is sparse regression only meant for multivariate models? Can a univariate AR(100) be sparse if say 30 of the variables are non-zero lagged. Can somebody show the functional form of how multivariate model looks like and how is it different from sparse if there is any difference?   3. Any references where I can find examples of sparse regression models and books which introduce the concept.


title: Statistical significance of increase in correlation

text: I have a sample with some individuals, and also three variables measured for each individual, let's say $A$, $B$ and $C$. There is a mild correlation between $A$ and $B$ and a mild correlation between $A$ and $C$   * $r_{AB} = 0.44$, and   * $r_{AC} = 0.47$.  I want to test if this increase in correlation is significant. I have come across two ways:   1. Fisher's z-Test for comparing 2 correlation coefficients   2. Wolfe's t-Test for comparing dependent correlation coefficients I calculated the _p-values_ of each test to be, respectively, $0.53$ and $2.8\times 10^{-5}$. From my understanding, Wolfe's t-Test is most appropriate in this situation, but the huge difference in p-value is making me wonder if I did something wrong along the way... How exactly does one test whether one correlation coefficient value is significantly higher then another?


title: False discovery rate calculation in target-decoy matching context

text: A common strategy in mass spectrometry of biological molecules is to upload observed spectra to a server so that they can be matched to a LARGE database of theoretical spectra of known molecules (a.k.a. _target_ database). In order to control for false positives, a _decoy_ database consisting of incorrect/irrelevant spectra is used. I have been reading more into this subject and have come up some questions regarding the calculation of the FDR measure from this target-decoy strategy. The basic idea of the FDR value is very intuitive: $FDR = \frac{FP}{FP + TP}$ where FP and TP stands for false and true positives respectively. This makes perfect sense to me; if I'm trying to guess some peoples' names out of a phone book, and get 8 right and 2 wrong, I would have 2 _false_ out of 10 _total_ guesses, and thus my false discovery rate would be 20%. However reading this tutorial on how this is done in large scale on the servers, I got introduced to two different calculations, depending on whether or not the _target_ and _decoy_ databases are concatenated (page 2). I don't think that this is a typo as I found other occurrences * of the mysterious factor 2 in front of FP in scientific literature. However the motivation behind this is never explained (at least I couldn't find it). I would appreciate some insight on where this doubling comes from. Likewise I wonder whether or not FDR calculation this way **assumes** that the error rate for each spectra match is the same for the target database and decoy database (i.e. _assuming_ that getting 25 decoy hits _implies_ 25 target hits are also false positives). It's not really clear for me why the error rate has to be the same for the two databases. Any comments on this subject is also appreciated. * one such reference is Elias et al Nature Methods - 2, 667 - 675 (2005)


title: R: Calculating mean and standard error of mean for factors with lm() vs. direct calculation -edited

text: When dealing with data with factors R can be used to calculate the means for each group with the lm() function. This also gives the standard errors for the estimated means. But this standard error differs from what I get from a calculation by hand. Here is an example (taken from here Predicting the difference between two groups in R ) First calculate the mean with lm():                   mtcars$cyl <- factor(mtcars$cyl)         mylm <- lm(mpg ~ cyl, data = mtcars)         summary(mylm)$coef                          Estimate Std. Error   t value     Pr(>|t|)       (Intercept)  26.663636  0.9718008 27.437347 2.688358e-22       cyl6         -6.920779  1.5583482 -4.441099 1.194696e-04       cyl8        -11.563636  1.2986235 -8.904534 8.568209e-10      The intercept is the mean for the first group, the 4 cylindered cars. To get the means by direct calculation I use this:                 with(mtcars, tapply(mpg, cyl, mean))                   4        6        8          26.66364 19.74286 15.10000       To get the standard errors for the means I calculate the sample standard variation and divide by the number of observations in each group:                with(mtcars, tapply(mpg, cyl, sd)/sqrt(summary(mtcars$cyl)) )                   4         6         8         1.3597642 0.5493967 0.6842016       The direct calculation gives the same mean but the standard error is different for the 2 approaches, I had expected to get the same standard error. What is going on here? It is related to lm() fitting the mean for each group and an error term? _**Edited:_** After Svens answer (below) I can formulate my question more concise and clearly. For categorical data we can calculate the means of a variable for different groups is by using lm() without an intercept.                 mtcars$cyl <- factor(mtcars$cyl)       mylm <- lm(mpg ~ cyl, data = mtcars)       summary(mylm)$coef                Estimate Std. Error       cyl4 26.66364  0.9718008       cyl6 19.74286  1.2182168       cyl8 15.10000  0.8614094      We can compare this with an direct calculation of the means and their standard errors:                 with(mtcars, tapply(mpg, cyl, mean))                   4        6        8          26.66364 19.74286 15.10000             with(mtcars, tapply(mpg, cyl, sd)/sqrt(summary(mtcars$cyl)) )                   4         6         8         1.3597642 0.5493967 0.6842016       The means are exactly the same but the standard errors are different for these 2 methods (as Sven also notices). My question is why are they different and not the same? (when editing my question, should I delete the original text or adding my edition as I did )


title: Minimum sample size for CART

text: What would be the minimum sample size for a CART analysis? I've found this recommendation: > The original CART monograph discusses a study the authors performed working > with 215 observations and 19 predictors, where 37 records were of class 1 > and 178 of class 0. We think that this is example, with 37 examples in the > smaller class is close the smallest sample size you can usefully work with > CART. > > Recommendation: We suggest using a minimum of 100 records, with the target > variable distributed not more unbalanced than in proportions (1/3, 2/3) for > up to 30 predictors. We recommend repeated cross-validation to estimate the > out-of-sample (previously unseen data) performance. Would this minimum sample size allow dividing between training and test sample? Is it OK not to divide the sample if it is small, and only carrying out cross-validation?


title: VIF values in regression

text: Do high VIF values for a a particular variable $x$ just indicate that it is highly correlated with at least one of the other variables in the model? Does it specify which variables and how many variables that $x$ is correlated with?


title: Using and interpreting type I SOS in three-way ANOVA

text: My research question is from education: To test whether groups (university, department, gender) differ based on their scores in a test, and whether there are any interaction effects, I'll do 3-way ANOVA? Data is unbalances (unequal cells). My professor says I have to do SOS Type 1 (report parameter estimates) because of the unbalanced design. His rationale is to compare the effect of each main effect adjusted for the effects of the other? He also said he'd prefer gender to be the first, he'd want parameter estimates, and earlier he had said he prefered Cohen's d type effect sizes. My questions: 1) I put gender,university,department, than all 2way interactions, than 3-way interaction in the model. I found gender was not SS, and no SS interaction of it, what-so-ever. Fgender (1,b) = 1.83, p = .18, η2 = .006. I reported R-squared. Should I have proceeded with all other combinations? There are a lot, right? I dropped gender and redid the analysis with (2x2) [university, department, uni*dep]. 2a) I reported Funiversity (1,b) = 46.45, p < .01, η2 = .16. I looked at the estimated means, reported them, and calculated Cohen's d -0.8 ( the difference between group means divided by the square root of the mean square error). I also reported Unstandardized Buniversity = 0.48, t = 5.67, p < .01. and concluded university A were more likely to get a higher score than university B. I used an online effect size conversion calculator from eta-squared to cohen's d. It also η2 = .16 = d = 0.8. 2b) I redid the analysis with [department, university, uni*dep). I reported the exact stuff above for department. But the η2 = .13 and when I calculated Cohen's d it turned out to be 0.38. I used an online effect size conversion calculator from eta-squared to cohen's d, but it reported 0.70. What is wrong? 2c) I plotted marginal means for uni*dep. I put it there to show there is interaction. 2d) I came to the interaction. I knew it would not change in Type I SOS, so I reported Funiversity by department (1, b) = 11.80, p < .01, η2 = .04. So I said first two are practically the last is SS. A) My question is to be sure to go to my prof unembarrassed: Is this okay? What else can I do for the interaction? Can I not say uniBdepB < all others, which is obvious in the graph? B) If parameter estimates tell me the slope, how come they are not equal to estimated means? C) How can I interpret the parameter estimates for the interaction effect? Which has a value for uniAdepA but than it says the rest are redundant?


title: Relativization and tackling an overdispersed dataset

text: If there is a species that is hyperabundant in a dataset and makes up 284/1638 individuals yet originated from just 1 of 10 traps. Is a quassi-poisson GLM sufficient providing residual deviance is under 1? or does this species need to be relativized or otherwise before proceeding to model selection? I really dont know how to combat this species.... Platypezid


title: Random walk out of sample forecasting

text: Probably my question is a bit stupid, but I'm having some problems in writing down in R the out-of-sample forecasting with a Random Walk. I have a multivariate time series (y) and I want to estimate the out of sample forecasting (y(t+k)-y^(t+k/t) result with a RW, for k =1,6,12. I write down this code, it works, but it seems to give me too much low errors.               residuals1 <- rep(0,58)     residuals6 <- rep(0,58)     residuals12 <- rep(0,58)          y1 <- t(y[,1])     for (i in 1:58) {              residuals1[i] <- y1[134+i+1]-y1[134+i]           residuals6[i] <- y1[134+i+6]-y1[134+i]        residuals12[i] <- y1[134+i+12]-y1[134+i]      }      I think that I'm doing some mistakes, do you have any suggestions?


title: Omitted variable bias in time series

text: This is a brief question because my lecturer mentioned it today in class but I don't quite understand. Why is omitted variable bias not a major problem in time series analysis?


title: What to report from hierarchical regression?

text: Here are three regression equations: $$ \begin{align} Y &= \beta_1Z_{\text{math}}\qquad (1)\\\ Y &= \beta_2Z_{\text{math}} + \beta_3Z_{\text{science}}\qquad (2)\\\ Y &= \beta_4Z_{\text{department}} + \beta_5Z_{\text{math}} + \beta_6Z_{\text{science}} + \beta_7Z_{\text{university}} + \beta_8Z_{\text{gender}}\qquad (3) \end{align} $$ I have 3-level regression model. I want to control the effect of maths achievement alone for (maths and science achievement in model 2), and model 3 includes gender, university and department. B1math is SS in Model 1, B2science is SS in model 2, and B7university and B6science is SS in model three. I have 15, 25, and 35 R-squared estimates. I am not sure what to report, besides increase in R-squared? Should I report the standardized Beta weight of university alone (because science was controlled) or should I also report B of science? Do structure coefficient play a role in interpreting B?


title: Subsets vs pooling in regressions with interactions

text: Model: $y_{it} = \alpha_{i}+\beta_1 D_{it} + \beta_2 G_i + \beta_3 (G_{i}\times D_{it}) + \epsilon$ G is a group dummy D is a treatment y is an outcome data is a panel $\alpha$ is a fixed effect My question: if the only regressor in the model that is shared between the two groups is interacted by the group dummy, is there any sense in keeping the data pooled, or should it be subset and modeled separately?


title: Maximum Likelihood Estimation of Inverse Gamma Distribution in R or RPy

text: I am trying to fit a three parameter inverse gamma distribution to my data in either R or Python. I would like to do this using maximum likelihood estimation (MLE). The pdf of the three parameter inverse gamma is given by: ![enter image description here](http://i.stack.imgur.com/nOrUa.jpg) _Where Γ is the gamma function, ρ is the shape, α is the scale and s is the location parameter_ I haven't spotted an R package that can perform MLE to this distribution directly (if you know of one, please let me know!). So I think this leaves either:   * _(A) working out the log-likelihood function of the formula_   * _(B) transforming the data to a gamma distribution. However, this distribution only has two parameters so I'm not clear on how I would calculate the third parameter (I'm not a very mathematical person!)._ Any help on a method to use MLE to fit an inverse gamma distribution to my data would be much appreciated! Many thanks in advance.


title: How to get p-value for indirectly generated random variables?

text: This is an extension of a question I had previously asked. If the answer is a simple text-book question then feel free to reference it. The problem started out as a need to compare two experimental data sets, one demonstrating ideal conditions, and the other, demonstrating what happens when a specific event is triggered during the experiments. I end up with two sets of results, of **_non-gaussian distribution_** : Sample 1a (ideal): 32 data points $$ \mu_{1a} = 0.5505, \sigma^2_{1a}= 4.9047*10^{-5} $$ Sample 1b (non-ideal): 31 data points $$ \mu_{1b} = 0.5314*10^3, \sigma^2_{1b}= 5.4851*10^{-5} $$ While the difference is small, it is significant, and my interest revolves around that _difference in mean_. When I change the parameters of the experiment, I get another set of results: Sample 2a: 34 data points $$ \mu_{2a} = 0.5395, \sigma^2_{2a}= 2.3220*10^{-5} $$ Sample 2b: 33 data points $$ \mu_{2b} = 0.5271*10^3, \sigma^2_{2b}= 2.0283*10^{-5} $$ What I need to do is demonstrate that Sample 2, provides a different _mean difference_ then Sample 1, and I want to provide the p-value for the associated null-hypothesis. At this point, I don't mind assuming the sample difference is Gaussian (although I wouldn't for inter-sample analysis). Using Matlab, I initially planned on doing the following:               normcdf(0,abs(sol.meanLoss -  comp.meanLoss),sqrt(sol.ste^2 + comp.ste^2));      which calculates the area under the curve of a Gaussian distribution ( parameters as mean: mean difference between Sample 1 and 2 ( $(\mu_{1a}-\mu_{1b})- (\mu_{2a} -\mu_{2b})$), and std as the square root of the sum of the squares of the standard error of each sample ) from -Inf to 0. I'm not sure my answer is correct and I'm don't mind having it torn apart, but I would need something else. If the t-test is to be used, how ought I combine the Samples? Edit: note that I have calculated the standard error as $$ ste_1 = \sqrt{ \sigma^2_{1a}/N_{1a}+\sigma^2_{1b}/N_{1b}} $$


title: About correlation of ordinal variables having different number of categories and about correlation of mixed type of variables

text: I have three questions:   1. For data with different response categories (likert type or ordinal) what type of correlation should be calculated? I know for ordinal variables **Spearman's rank correlation** can be used. Am I right?    2. Do I need to order the categories and recode the responses according to their importance to calculate these correlations? I think I should clarify my question. Suppose that I have two variables that I believe responsible behind profitability of a company. I have the following variables: For variable 'regular staff meetings', responses are collected as: 1='yes' and 2='no' For variable 'proportion of staffs given formal performance appraisals', responses are collected as: 1=10%, 2=20%,..., 10=100% Now as intuitively I believe that the more regular staff meetings are held or performance appraisals are given the better will be the profit of that company, so should we recode for variable 'regular staff meetings' 'yes'=2 and 'no'=1 that responses with 'yes' receives the first rank while spearman's rank correlation is calculated? Similarly for variable 'proportion of staffs given formal performance appraisals' responses with 10 (100% staffs receiving performance appraisal) will receive the first rank while spearman's rank correlation is calculated. Otherwise if I don't recode the responses from the variable 'regular staff meetings' then the ranks the responses receive becomes of opposite direction as responses with 'yes'=1 receives the second rank but it was supposed to increase profitability!   3. What if I want correlation between mixed type of variables (both **interval** and **ordinal** )? Is it possible to calculate? Thanks for your patience and upcoming kind clarification. :)


title: Difference between SOM and Hopfield

text: I am learning about SOM and also about Hopfield networks. I yet don't understand why I would use `SOM` over `Hopfield` and vice versa. What kind of problems I should resolve with SOM and which with Hopfield?


title: BIC or AIC to determine the optimal number of clusters in a scale-free graph?

text: I am currently trying to partition a scale-free ("big") graph (around 20k vertices, 500k edges) into appropriate sub-graphs. Having derived the Laplacian of the graph, I tried running an approach based on the spectral gap and Fiedler-vector, however, not really unexpectedly, ended up with vertex valuations (i.e. components of the corresponding eigenvector) being close to zero for a majority of the nodes. Clearly, there is no obvious cut in the graph. Nevertheless, even if it's just for the sake of showing that several methods fail on graphs following the spectral features of the one I am working on, I would like to further explore spectral clustering approaches - some of which require a fixed k denoting the number of partitions. I am aware of the use of the BIC and AIC with respect to k-means-clustering. What would interest me is, whether these criteria are also used in the realm of spectral graph clustering? Is there any rationale allowing to establish a link between the spectra of graphs and model selection criteria like the BIC and AIC? Any input is much appreciated! * * * Additions: So, I have run a few tests. I have tried RSB with the median for the cutoff value c. I used high-evidence (low false positive rate, possibly high false negative rate) cluster-data to validate against (roughly ~250 non-overlapping groups), in a rather "poor man's" fashion, so nothing fancy at all. The initial cut already affected more than 235 clusters, even though many of them are actually rather small (we are talking about an avg. of around 75 here). I tried deviating from the median by the MAD (towards the valuation with the highest absolute value) which resulted in a bad performance as well. After some further tries, I ended up choosing the 1st- or 3rd-quartile of the valuation distribution, which allowed some small and rather trivial cuts. Nevertheless, the spectral gap never looked promising and the characteristic valuation simply horrible. For computing them I used ARPACK (IRLM), so I expect the results to be considerably accurate in double precision. Here's a plot of the characteristic valuation (log2, just quick and dirty) after the first 2 iterations (which both yielded 2 clusters of roughly 36 nodes each) - the core seems to be too dense. ![run0](http://666kb.com/i/c5lyabals2hkrir2u.gif) ![run1](http://666kb.com/i/c5lyax7pf6uega93q.gif) I thought about at least buying Fan Chung's more recent book on spectral clustering (spectral clustering), since I liked reading through the previous work (at least the first two chapters). They were dry to the bone, but nevertheless quite informative.


title: Is the absolute value of the difference between two Poisson distributions a Poisson distribution?

text: My question is basically as follows: what is the distribution of the absolute value of the Skellam distribution. Thanks!


title: Is there an "initial condition" for ARMA model?

text: 1. ARMA model is a stochastic version of recursive relation. For deterministic recursive relations, we solve them and need initial conditions to fully get the solution. So I wonder what is an "initial condition" for ARMA model like?   2. In Introduction to Time Series and Forecasting, by Peter J. Brockwell, Richard A. Davis, they define an ARMA model to be causal, if its output process can be represented as MA($\infty$) (see equation (3.1.5)). In equation (3.1.7), they gave the explicit formulation of the MA($\infty$) representation. I don't find they specify or use any "initial" condition in the definition of causality and in the derivation of the MA($\infty$) representation. This is quite different from that we need initial conditions to fully get the solution to a deterministic recursive relation. So how shall we understand the MA($\infty$) representation of an ARMA output process? Thanks!


title: How are ratings generally distributed?

text: Do user ratings of products — games, films, books, restaurants, etc. — tend to follow any particular distributions? If so, what form do they commonly take? Are there public data sources that show these distributions (e.g. for Amazon books, Netflix movies, etc.)?


title: What statistical test should I use for this situation?

text: I've got three groups A,B,C. Each has a different number of records and I have calculated the mean for each group. Now I want to know if these means are equal or not? What statistical test should I be doing?               Group   No.ofrec  Mean     A       100       10     B       200        5     C       300       20


title: Sales prediction of Amazon Sales Rank

text: I understand that Zipf's Rank-size distribution would allow me to estimate the sales of a book given the amazon sales rank. How can I fit the distribution without knowing book sales leading to a rank. Sales Rank is calculated with time since last sale in mind, a decay function.


title: Distribution of a MA(1) process

text: Suppose I have this MA(1) model: $y_t = \mu + \epsilon_t + \theta \epsilon_{t-1}$ with $\epsilon_t \sim \mathcal{N}(0,\sigma^2)$ The marginal distribution of $y_t$ for all $t$ is $\mathcal{N}(\mu,\sigma^2(1 + \theta^2))$ but the classical theory says that since the $y_t$ variables are not independent, the likelihood of the sample is not the product of the marginal densities. One has to use conditional densities instead (given $\epsilon_0 = 0$): $y_1 \sim \mathcal{N}(\mu,\sigma^2)$   $y_2 \sim \mathcal{N}(\mu + \theta \epsilon_1,\sigma^2)$   $...$   $y_t \sim \mathcal{N}(\mu + \theta \epsilon_{t-1},\sigma^2)$ And the likelihood is the product of these conditional densities, right? Now, I wrote a small program in R simulating a MA(1) process and displaying the distribution of the generated $y_t$ variables:               n = 100000     mu = 5     theta = 3     sigma = 2          y = rep(1, n)     eps.1 = 0          for (t in 1:n)     {         eps = rnorm(1, mean = 0, sd = sigma)         y[t] = mu + eps + theta * eps.1         eps.1 = eps     }          hist(y, probability = TRUE)     curve(dnorm(x, mean = mu, sd = sigma * sqrt(1 + theta^2)), col = "blue", add = TRUE)      And I find out that the joint distribution of $y_t$ is $\mathcal{N}(\mu,\sigma^2(1 + \theta^2))$. How come?


title: Probability that a sum of potential numbers is greater than some value

text: Say I am about to receive 5 cash prizes and I have the probability of receiving each cash prize. Let's denote a set of cash prizes with $k$. So, below is the set of cash prizes and the set of corresponding probabilities: $$k=\\{65,25,30,54,30\\}$$ $$p(k)=\\{0.8,0.5,0.25,0.2,0.4\\}$$ Hence, the expected value of my cash prizes is $$E[k]=65*0.8+25*0.5+30*0.25+54*0.2+30*0.4=94.8$$ So, my question is: what is the probability that the sum of the cash prizes will be greater than 85\$? 90$?. **Attempt** : I did not come up with anything meaningful. I do not know if the question is valid in the first place. Thanks for your help.


title: Examples of machine learning in personal computers

text: Are there any examples of machine learning in our PC's? The only one that I know is Windows Speech Recognition application. I do not mean any add-ons. I mean any basic applications that are available in most PC's. **EDIT:** Any reference would be helpful.


title: What test/method is suited for analysing and modeling answers from a survey?

text: People got a list of different documents and were asked if they find things important or not and if they think the documents would be easy to get or not. People belong to different groups and have different professions. For answering they had three choices:   0 = not important   2 = moderately important   4 = very important and the same for easiness:   0 = not difficult   2 = moderately difficult   4 = very difficult   _no answer_ was also possible, so there can be also missing values. My question is: How could the data be analysed with statistical methods to find out   * how profession and group influence what is important for the people? (= effect and interaction of those 2 variables on the "importance-value")   * and for which documents there are statistically significant differences what is important for them Expected result (simple example): for a salesman the product data sheet is very important, someone working in service will need the documentation more often etc. and marketing people often need press releases. What kind of test or analysis method would be appropriate to analyse this kind of data? For sure, this is no designed experiment, so the number of people per group and per profession is not always the same. I assume this is kind of a standard problem - where can I find a good tutorial for such kind of analysis?


title: Longitudinal / panel data in component failure analysis

text: We just had a courses on panel data in econometrics, and longitudinal data in classic statistics, And I was wondering if something similar exists in engineering, for example: You have 50 identical cars and you ride with them until a particular component fails. You use the data of the car (forces on the car, accelerations, impacts, etc...) as predictors and you use both the time dimension and the cross sectional dimensions to predict the chance that a component fails on another car Or you have 50 wind turbines, and after a couple of years you see that a component fails. Again you want to know the chance that these components fails on another turbine, given the power the turbine produced, the turbulation it sees, etc... Given that a failure is binary data, would it be smart to model something like a GEE model on this data set? Thank you.


title: How to interprete results from Poisson random coefficients analysis?

text: I have obtained the following results in R from a random coefficient Poisson analysis. Generalized linear mixed model fit by maximum likelihood ['glmerMod'] Family: poisson ( log ) Formula: frequency ~ 1 + insgen + ageveh + make + area + (1 | ID) Data: Panel                 AIC       BIC    logLik   deviance       1099.9670 1134.9262 -541.9835 1083.9670 Random effects: Groups Name Variance Std.Dev. ID (Intercept) 1.551e-11 3.939e-06 Number of obs: 584, groups: ID, 584 Fixed effects:                         Estimate    Std. Error   z value   Pr(>|z|)        (Intercept) -22.98292 8432.07738 -0.003 0.9978 insgenM 0.02616 0.08806 0.297 0.7664 ageveho 0.05889 0.08586 0.686 0.4928 make -0.10447 0.04126 -2.532 0.0113 * area1 23.68571 8432.07738 0.003 0.9978 area2 23.85969 8432.07738 0.003 0.9977 area3 23.77374 8432.07738 0.003 0.9978 * * * Signif. codes: 0 ‘ ** _’ 0.001 ‘** ’ 0.01 ‘_’ 0.05 ‘.’ 0.1 ‘ ’ 1 Correlation of Fixed Effects:                   (Intr)  insgnM  ageveh  make    area1   area2       insgenM 0.000 ageveho 0.000 -0.037 make 0.000 0.071 0.108 area1 -1.000 0.000 0.000 0.000 area2 -1.000 0.000 0.000 0.000 1.000 area3 -1.000 0.000 0.000 0.000 1.000 1.000 I have included predictors:gender (Male,female), area which has 4 levels (area 1,2,3,4) and vehicle age (New and old) and make of car.


On which data should the lift be calculated i.e. Training set or Test set and why? What does the lift value 115% mean


title: Choosing the right forecasting technique

text: I'm currently attempting to forecast visitor data for stores. My dataset includes daily visitor totals of three years. Note that the dataset isn't complete (stores can be closed for a few days, etc). I've augmented the dataset with weather data per day and holiday information (I have about 18 features per datapoint, including weekday,month,holiday,weather-data). My computer science background doesn't include a lot of statistical methodology, so I have a hard time figuring out which techniques I can use/should look into. My focus right now is on 1-ahead forecasting with a simple neural network in pylearn2, but the predictions that I'm getting out of it are highly inaccurate. I would like to evaluate some other techniques for comparison and evaluation. I could use some hints in the right directions, so I'm wondering: Which techniques could I look into for forecasting these numbers? And how could I simplify the problem to increase my succes?


title: Distribution of a Mean and Variance

text: Say we have observations $x_1 \dots x_n$ and we have some sort of Bayesian framework where we would like to estimate a distribution for the mean $\mu$ of our observations and the variance $\sigma^2$ of our observations (these distributions will then be updated based on further observations). I have 2 questions: 1) By the central limit theorem, for large $n$, $\hat{\mu} \sim N(\mu, \frac{\sigma^2}{n})$. I was therefore thinking of using this as my prior distribution for $\mu$. If I was going to do the same thing for the variance, what would the limiting distribution be? 2) Would it perhaps be better to choose many different size $k$ subsamples of my observations, calculate the sample mean in each case and use this as an empirical prior distribution for $\mu$ (and do the same thing for the variance)? P.S. Any links to further information on the topic would be appreciated.


title: Problem defining ARIMA order

text: This is a long post so I hope you can bear with me, and please correct me where I'm wrong. _My goal is to produce a daily forecast based on 3 or 4 weeks of historical data._ The data is 15 minute data of the local load of one of a transformer lines. Im having trouble finding the model order of a seasonal ARIMA process. Consider the electricity demand time series: ![Original Time Series](http://i.share.pho.to/80d86574_l.png) When the first 3 weeks are taken as subset and differenced the follwing ACF/PACF plots are computed: ![Subset](http://i.share.pho.to/5c165aef_l.png) ![First difference](http://i.share.pho.to/b7300cc2_l.png) ![Seasonal and first difference](http://i.share.pho.to/570c5397_l.png) This looks like the series is kinda stationary. But the seasonality could also be weekly (see Seasonal difference week and second order differences [here]http://share.pho.to/3owoq, what do you think?) So let's conclude that the model take shape of : $$ ARIMA(p,1,q)(P,1,Q)_{96} $$ In the last figure a distinct spike at lag 96 indicates a seasonal MA(1) component (maybe AR(1) could be aswell as there is a distinct spike in PACF aswell). The spikes at lag 1:4 indicate a MA(4) component which corresponds with exponential decay in the PACF with a bit of imagination. Thus the initial model manually selected could be: $$ ARIMA(0,1,4)(0,1,1)_{96} $$ with               Series: x      ARIMA(0,1,4)(0,1,1)[96]               Coefficients:         ma1      ma2      ma3      ma4     sma1         -0.2187  -0.2233  -0.0996  -0.0983  -0.9796         s.e.   0.0231   0.0234   0.0257   0.0251   0.0804              sigma^2 estimated as 364612:  log likelihood=-15138.91         **AIC=30289.82   AICc=30289.87   BIC=30323.18**      The auto.arima function computes the following model (with stepwise and approximation on TRUE, else it takes to long to converge): $$ ARIMA(1,1,1)(2,0,2)_{96} $$ with               Series: x      ARIMA(1,1,1)(2,0,2)[96]               Coefficients:         ar1      ma1    sar1    sar2     sma1     sma2         0.7607  -1.0010  0.4834  0.4979  -0.3369  -0.4168         s.e.  0.0163   0.0001  0.0033  0.0116   0.0216   0.0255              sigma^2 estimated as 406766:  log likelihood=-15872.02         **AIC=31744.99   AICc=31745.05   BIC=31784.25**      Which means no seasonal differencing is applied. Here are the residuals of the both models. The Ljung Box statistic give a very small p value, indicating that there is still some autocorrelation present (? correct me if im wrong). ## Forecasting Thus to determine which is better, an out-sample accuracy test is then the best. So for both models an forecast is made 24 hours ahead which is compared with each other. The results are: ![auto.arima](http://i.share.pho.to/5d1dd934_l.png) ![manual model](http://i.share.pho.to/7ca69c97_l.png) Auto:                                     ME     RMSE      MAE       MPE      MAPE      MASE        ACF1 Theil's U     Training set   -2.586653 606.3188 439.1367 -1.284165  7.599403 0.4914563 -0.01219792        NA     Test set     -330.144797 896.6998 754.0080 -7.749675 13.268985 0.8438420  0.70219229  1.617834      Manual                                      ME     RMSE      MAE        MPE      MAPE      MASE         ACF1 Theil's U     Training set 2.456596e-03 589.1267 435.6571 -0.7815229  7.509774 0.4875621 -0.002034122        NA     Test set     2.878919e+02 919.7398 696.0593  3.4756363 10.317420 0.7789892  0.731013599  1.281764      ## Questions As you can think of this is an analysis on the first three weeks of a dataset. I'm struggling in my mind with the following questions:   1. How do I select the best ARIMA model (by trying all different orders and checking the best MASE/MAPE/MSE? where the selection of performance measurement can be a discussion in it's own..)   2. If I generate a new model and forecast for every new day forecast (as in online forecasting), do I need to take the yearly trend into account and how? (as in such a small subset my guess would be that the trend is neglible)   3. Would you expect that the model order stays the same throughout the dataset, i.e. when taking another subset will that give me the same model?   4. What is a good way, within this method to cope with holidays? Or is ARIMAX with external holiday dummies needed for this?   5. Do I need to use Fourier series approach to try models with `seasonality=672` as discussed in Long seasonal periods ?   6. If so would this be like `fit<-Arima(timeseries,order=c(0,1,4), xreg=fourier(1:n,4,672)` (where the function fourier is as defined in Hyndman's blog post)   7. Are initial P and Q components included with the fourier series? Most theoretical knowledge obtained from FPP, great stuff! Before advising on using exponential smoothing or (dynamic)linear regression this is also being worked on to compare. ## Data https://www.dropbox.com/sh/mzx61sskya5ze6x/Zq3A7Q6htH/trafo.txt # Code               data<-read.csv("file", sep=";")     load<-data[,3]      I removed the few zero values with week before values               stepback<-672     load[is.na(load)] <- 0 # Assumed no 0's in first 672 values!     idx <- which(load == 0)     idx <- idx[which(idx>stepback)]      load[idx] <- load[idx-stepback]           ED<-ts(load,start=0, end=c(760,96),frequency=96)     x<-window(ED,start=0, end=c(20,96))      It is also possible to post a reproducible example but this will make the post even longer, but possible if needed. So if there is anything I should provide please let me know, thanks in advance.


title: What kind of analysis one should do on a course

text: A course problem given by some unknown in a discussion forum is the following: A student has been given data of age at death of novels (67 data points), age at death of poems (32 data points) and age at death of nonfiction (24 data points). Then lecturer asked to use graphical and numerical methods to describe the data without saying what kind of information he wants on the data. What should a student do? Ask the lecturer what is really on his mind? Do some plots of the data he or she has been taught? Should the course exercises be well defined or is it common that one can pass the courses only if he or she makes an educated guess at what the lecturer wants to be contained in the answer?


title: R: How to create a function from a model?

text: I am using an automatic model selection procedure, "step". The model of depart (the biggest possible) is a polynom, say of the 4th degree.               Depart<-lm(y~x+I(x^2)+I(x^3)+I(x^4))          Final<-step(Depart)      I need to transform the Final model to a corresponding function. How can i do this?


title: Response surface of a particular discontinuous function

text: I have a function `IR` that depends on several (maybe 100) `input random variables`. I know `IR varies between 0 and 1`. Therefore it has two discontinuities, one at `IR=0` and the other at `IR=1`. I also know that the `pdf of IR` takes non-zero probability values at 0 and 1 and is continuous between 0 and 1. That is, the pdf of `IR` is a mixture pdf with discrete probabilities at 0 and 1 and a continuous distribution between 0 and 1. **My 1st question is** : If my training set and testing set are large enough to cover approximately the entire range of the input random variables and the accuracy of the regression of my response surface is very good (95% or higher for example), _if I change the distributions of the input random variables (but not their minimum and maximum values) in what conditions can I still use the previously obtained response surface of IR to get estimates of the pdf of`IR`?_ **My 2nd question is** : Since `IR` takes continuous values between 0 and 1, if I use SVM or random forest methods to get a response surface by regression do I need to discretize IR in `n` features? If I'm right, this means I have to divide the IR interval by 1000 and so IR would take values 0, 0.001, 0.002,...., 0.999, 1. _However, if my testing or training set contain a value equal to 0.0005 this result would be classified in what feature?_ **My 3rd question is** : _Is it preferable to define`IR` with a continuous function, from -infinite to +infinite, to obtain better response surface accuracies, possibly using other types of reponse surface methods and not only SVM or random forest methods?_ Thanks Joao


title: How do we know which variables (X1,X2,....,X6) are affecting the outcomes (y)?

text: The data I am using are sales person data. The data contain six independent variables (six vulnerabilities) scored 1-20 and the output is win/loss of the deal. My question is how do we know which variable or combination of variables are affecting the output i.e. win/loss, and what should be the methods or step to come up with my objective?


title: Outliers for Normalization: Is it important?

text: I'm new to statistics and have to do some research analysis using SPSS. I need to know, what is the purpose of removing outliers? I want to normalize my set of data obtained through a national survey, but there are a few variables that do not show a bell curve even if I have transformed my data. log, ln, sqrt etc. I have tried everything. What else can be done? I have to do multiple regression analysis on the data to predict future data.


title: Variance-Covariance matrix interpretation

text: Assume we have a linear model `Model1` and `vcov(Model1)` gives the following matrix:                            (Intercept)    latitude  sea.distance   altitude     (Intercept)    28.898100 -23.6439000  -34.1523000  0.50790600     latitude      -23.643900  19.7032500   28.4602500 -0.42471450     sea.distance  -34.152300  28.4602500   42.4714500 -0.62612550     altitude        0.507906  -0.4247145   -0.6261255  0.00928242      For this example, what does this matrix actually display? What assumptions we can safely make for our model and it's independent variables?


title: How to deal with outlier in the experiment

text: My question is related to this one. I'm working on an experiment involving transferring data between different locations all over the world via the Internet. As the Internet can be busy, sometimes it takes longer to transfer the same amount of data than the other times. The difference can be significant, e.g. for the 1mb of data, it normally takes 2 seconds to transfer, but occasionally (around 1 out of 20 times), it can take up to 20 seconds or more. I am comparing different setups for the same problem. For each one, I perform multiple times and get the average result, i.e. time to solve the problem. However, as outlier can happen at one setup but not on another, my comparison is affected. So, at the moment, for each setup, I perform M times but only get the average of N of the smallest results (N < M). I also plan to state and justify my method in the paper as well. So, my question is whether my method is acceptable and if it makes my paper less reliable. I'm also interested in how experimental outlier is handled in general. I'm working in CS but answers from other fields are welcome.


title: A clinical test incorrectly identified 50/150 false positives and 23/150 false negatives.Calculate the number of true positives and true negatives

text: Was wondering if anybody could explain how I can calculate the answer. Would be much appreciated.


title: Using R's prop.test function to compare subsets of a population

text: I'm new to CV, so please excuse any faux pas. I've collected sample data from a population and would like to discern whether or not the experimental group outperforms the control group. Let's say I have 10,000 observations within the experimental group and saw 4500 successes. Let's also say that I observed a 47% success rate in my control group (so if my control group is 5% of the total sample size, that would be 235 successes out of 500 trials). Now, am I right in my usage and interpretation of the following?:               prop.test(4500,10000,.47,alt="less")          data:  4500 out of 10000, null probability 0.47     X-squared = 15.9776, df = 1, p-value = 3.205e-05     alternative hypothesis: true p is less than 0.47     95 percent confidence interval:      0.0000000 0.4582455     sample estimates:        p      0.45      Since the p-value is far less than .05, I can safely reject the null hypothesis that that the experimental group will perform as well or better than the control group. Also, I know that 95% of the time, the probability of success for the experimental group will fall between 0% and ~45.8%. Also, is the `prop.test` function a good indicator of sufficient sample size? If not, what is? Thanks!


title: Discovering Values with Forecast

text: So, I have numbers that I know are going to be between 1.01 and 1000. Each value corresponds to 1 second. Is this is a time-series? I've been searching all over, but I can't figure out if these values represent an exponential curve, a linear curve (I think not), or something else. And witch forecasting model to do what I want. I want to find the next 5-10 values. I think you can help me with this, but this is another issue, or not. The file is linked to another program (not this file) to capture the values. These values are going to be captured second by second. It's a dynamic graphic, not a static one. This is the values: 8,6 8,6 8,6 8,6 8,6 8,6 8,6 8,6 8,6 8,6 8,4 8 8 8 8 8,6 8,6 8,2 8,2 8,2 8 8 8 7,8 8 8 7,6 7,4 7,4 7,4 7,4 7,2 7 7 8 8 8 7,8 7,8 8,4 8,4 8,2 8,4 8,4 8,2 8,2 8,2 7,2 7,6 7,8 7 6,8 7 7 7,2 7,4 7,4 7,2 7,6 7,4 7,8 7,4 6,4 7,4 7,2 7,2 6,4 6,2 6,4 6,4 6,2 6,2 6,2 6,2 6 5,6 5,5 6,6 6,6 6,8 6,8 7 7,4 6,8 7,8 7,8 7,8 7 7 7,8 6,2 7,2 6,2 6,2 6,4 6,8 6,6 7 5,4 5,1 6,6 5,8 7,4 9 9 8,4 9 11 11 11 12 12 12 25 25 16 30 20 40 40 40 95 95 95 100 100 150 280 280 280 990 990 990 990 990 990 990 990 990 990 990 990


title: Calculating relative importance of independent variables in Linear and Logistic regression

text: After fitting Logistic/Linear regression model, we get estimates of parameters showing importance of attribute in predicting dependent variable and can be considered as weight of that independent variable. My question is how to calculate relative importance of attributes in linear and logistic regression. Till now I, used to take mod of each coefficient and map sum of coefficients to 100 (i.e. calculate percentage of weight of attribute). Is it the right way.


title: How should I interpret GAP statistic?

text: I used GAP statistic to estimate k clusters in R. However I'm not sure if I interpret it well. ![enter image description here](http://i.stack.imgur.com/jTY0S.png) From the plot above I assume that I should use 3 clusters. ![enter image description here](http://i.stack.imgur.com/HTia8.png) From the second plot I should choose 6 clusters. Is it correct interpretation of GAP statistic? I would be grateful for any explanation.


title: Determining trend significance in a time series

text: I have some time series data and want to test for the existence of and estimate the parameters of a linear trend in a dependent variable w.r.t. time, i.e. time is my independent variable. The time points cannot be considered IID under the null of no trend. Specifically, the error terms for points sampled near each other in time are positively correlated. Error terms for samples obtained at sufficiently different times can be considered IID for all practical purposes. I do not have a well-specified model of how the error terms are correlated for points close to each other in time. All I know from domain knowledge is that they are positively correlated to some degree or another. Other than this issue, I believe the assumptions of ordinarly least squares linear regression (homoskedasticity, linearity, normally distributed error terms) are met. Modulo the correlated error term issue, OLS would solve my problem. I am a complete novice at dealing with time series data. Is there any "standard" way to proceed in these circumstances?


title: Finding the mode of a function by MCMC sampling

text: When trying to find the mode of a nonnegative function $f$ (i.e. maximize the function), one way to do it is to sampling the function viewed as an unnormalized density of some distribution via MCMC. Suppose we have had a sufficiently long sequence of samples via this method, I was wondering how to determine the mode from the samples? Specifically, as we know that a part taken from the end of the sequence will be assumed to be approximately subject to the distribution corresponding to the target function $f$. As far as I have thought already, there might be two choices for the estimation of the mode of the function:   1. take the last sample in the original sequence,    2. take a small subsequence from the end of the original sequence and evaluate the function f on every sample in the subsequence and pick the one with the maximal function value.  I saw the first choice from my class note, but the second was my first thought before looking at the note. So I was wondering what choice for the mode of the target function might be reasonable or better and why? Possible to have any references on this? You don't have to limit your scope to the two I just mentioned. Thanks and regards!


title: Appropriate test for multivariate experiment result with unknown distributions

text: I have an experiment with any possible (reasonable) number of parameters (independent variables). I run the experiment several times for each possible combination of my variables. The data I get will be generally numeric. However I know nothing (and any assumptions are difficult) about the distribution of my data. What I am interested in is a measure of how well do my parameters predict the data I get. Which statistic should I use? How do I calculate it (by hand, a link to a tutorial would be very sweet)? ### Edit I am trying to solve this as generally as possible (hence the slightly non- specific description) for a piece of software I'm working on. To make it a bit more clear a bit of an example: I have these parameters:               decay: 0.1 | 0.2 | 0.3     particles: 10 | 100     velocity: 30 | 70      This gives 12 combinations (3 * 2 * 2) and I'll measure my dependent variable (say temperature) five times for each combinations. Thus my final dataset will have 60 measurements of temperature. Now suppose that temperature was in fact given by: $t = K(0.3v + 0.6d + \varepsilon)$ where $t$ is temperature, $K$ is some constant, $v$ is velocity, $d$ is decay and $\varepsilon$ is some sort of random effect. Particles is completely unrelated to the measured temperature. Now I'd like to perform a test that would tell me that velocity has a ~0.3 effect, decay ~0.6 and particles ~0 effect. However I may have more or less variables and more or less measurements.


title: Statistical classification of text

text: I'm a programmer without statistical background, and I'm currently looking at different classification methods for a large number of different documents that I want to classify into pre-defined categories. I've been reading about kNN, SVM and NN. However, I have some trouble getting started. What resources do you recommend? I do know single variable and multi variable calculus quite well, so my math should be strong enough. I also own Bishop's book on Neural Networks, but it has proven to be a bit dense as an introduction.


title: Bayesian references

text: Any good places to start getting into Bayesian statistics? I'm a graduate student in social sciences, with a decent amount of stats classes under my belt, but I'm far from fluent. Any references would be much appreciated. Thanks!


title: Recommendations for learning probability and Bayesian statistics?

text: I have been very interested lately in learning Bayesian Statistics, but I have only a little bit of background in the frequentist statistics, only one term at University. Some of the books that I have seen are extremely mathematically oriented like:   * Bayesian Theory, José Bernardo. Probability and Statistics, Degroot   (this is actually a nice book, but it jumps too much between chapters in the examples, making the process of reading pretty sequential in all the book) The topics that I would like to learn are:   * Bayesian inference.   * Generation of independent samples from distributions.   * Monte Carlo integration, importance sampling.   * Posterior distribution with numerical quadrature or Laplace expansion.   * MCMC methods: Gibbs and Metropolis-Hastings sampling.   * Auxiliary variable methods in MCMC.   * EM algorithm.   * Multi-model inference.   * MCMC theory. I have tried to find something in Coursera but nothing. What books or online courses do you recommend?


title: How does one interpret a Bland-Altman plot?

text: In Plain English, how does one interpret a Bland-Altman plot? What are the advantages of using a Bland-Altman plot over other methods of comparing two different measurement methods?


title: What statistical test is appropriate when comparing heart rates from different devices?

text: I am checking the accuracy of one device that finds the heart rate against the other device that gives correct measurements. Specifically, I have performed an experiment on a group of people. I used both device 1 and device 2 simultaneously on the each person and collected 7 values. Having 10 participants, I have in total 70 values for each device. I placed the data from all participants together and I have two columns, one for each device, of 70 values. Now, I need to compare the two set of values to find whether they are significantly different or not. I am not sure what test to use here. If I use the t-test, which one should I use? Paired, Two Samples assuming equal Variances or Two samples assuming unequal variances? I am a bit confused.


title: Calculating pointwise mutual information between two strings

text: I have a dataset consisting 5000 sentences. I need to calculate PMI between 3-gram and 5-grams in this dataset. For example: > The 5-gram is: $x_1$$x_2$$x_3$$x_4$$x_5$ > > And the 3-gram is: $x_2$$x_3$$x_4$ How can I calculate PMI($x_2$$x_3$$x_4$, $x_1$$x_2$$x_3$$x_4$$x_5$) in this dataset? What is the exact formula? As far as I know, in order to calculate PMI(y,z), it is needed to keep track of these counts in the dataset: > Count(y,z) --> The number of co-occurrences of y and z > > Count(i,z) --> The number of occurrences of z > > Count(y,i) --> The number of occurrences of y > > N --> The sample size The final formula is: > $PMI(y,z) = \frac{Count(y,z)N}{Count(i,z)Count(y,i)}$ In the context of my problem, these counts are listed as below: > $Count_1$($x_2$$x_3$$x_4$, $x_1$$x_2$$x_3$$x_4$$x_5$) : All occurrences of > exact 5-gram " $x_1$$x_2$$x_3$$x_4$$x_5$ " > > $Count_2$(i, $x_1$$x_2$$x_3$$x_4$$x_5$) : All occurrences of exact 5-gram > "$x_1$$x_2$$x_3$$x_4$$x_5$ " !!!!! > > $Count_3$($x_2$$x_3$$x_4$, i) : All occurences of 5-grams in the form of " _ > $x_2$$x_3$$x_4$ _ ", that is the first and last words are subsituted with > all the possible words in the dataset. > > $N$ : All the 3-grams. !!!!! My problem is with $Count_2$ and $N$. As you see, $Count_2$ is equal to $Count_1$. Is it sensible? And I'm not sure about the way I have counted $N$.


title: Fisher's exact test, very small sample size, RxC table, why not work

text: I apply Fisher's exact test on very very small sample size. Each location in the table would have about 1 or 2 samples. The fisher's method has no room to enumerate the the tables, and returns me a very bad result. So is there any other improved exact test for very tiny sample size data?


title: Error in an article with a conditional probability?

text: I’ve recently read the article "Visual Tracking of Human Visitors under Variable-Lighting Conditions for a Responsive Audio Art Installation," A. Godbehere, A. Matsukawa, K. Goldberg, American Control Conference, Montreal, June 2012. On page 4 it says: > Making use of Assumption I-C5, we let $p(f|F) = 1-p(f|B)$ [$f$ is the feature calculated and $F$ stands for foreground, $B$ for background]. Assumption I-C5 is given on page 3. The quoted sentence is a non sense whatever assumption one makes. $p(f|F)$ is then replaced in the following relation, from earlier on the same page (which is just Bayes’ rule): > $$p(B|f) = \frac{p(f|B)p(B)}{p(f|B)p(B) + p(f|F)p(F)}$$ and the final Boolean result comes from a threshold on the above formula, where $p(F)$ and $p(B)$ are adjustable parameters (now it's really $p(F) = 1- p(B)$). I suspect that the authors noticed that this approach worked and accepted it.   Maybe the right way to proceed would have been to fix a constant p(f|F), as a uniform distribution models a situation where we have no way to know it, and so it will become a new adjustable parameter. Like in the assumption made by the authors, we have an increasing function, $p(B|f)$, of variable $p(f|B)$ with a few constants. It is equivalent to the above approach but with a different choice of parameters and a different threshold. Any thoughts about it?


title: Binomial probability in JMP

text: I'm trying to do a problem in JMP, but I'm new to the program and don't understand it completely. I think I'm close, but need to be pointed in the right direction. Here is the problem: > A standard test for extrasensory perception (ESP) asks subjects to identify > which of four shapes (circle, square, diamond, or waves) is on the front of > a card, viewed by the experimenter but not the subject. Suppose that a > subject takes a test with 20 of these cards (each card is equally likely to > display one of the four shapes). In this investigation, we will focus on the > number of correct identifications made in 20 trials. > > a. Use the binomial distribution) to find the probability that a guessing > subject would get 7 or more correct. How do you get the probability for 7 or more using the binomial probability in JMP? Here is a screen shot of what I've done so far in JMP. The values I used were 20 for N, and .25 for pi: ![enter image description here](http://i.stack.imgur.com/WhBCe.jpg) Does what I have look correct so far? Where do I go from here?


title: Having trouble understanding netflix RBM

text: According to this paper (pdf), the energy function of the restricted Boltzmann machine (RBM) is defined as: ![RBM Energy Function](http://s24.postimg.org/54331he3p/2014_02_05_01_18_29.png) and the paper shows that the conditional probability of softmax unit is: ![Conditional Prob](http://s22.postimg.org/nik01545t/2014_02_04_21_29_08.png) I'm having difficulty understanding how this conditional probability can be derived from the energy function using the Boltzmann distribution. Can anyone help me understand how to derive $p(v=1|h)$ and $p(h=1|v)$ from the given energy function?


title: Finding the state for Markov chain in R based on a generated sample

text: I am looking for a R code of discrete Markov chains. I have done the following:               set.seed(123)          p_11 <- 0.0468      #transition probability         p_12 <- 1 - p_11    #transition probability         p_21 <- 0.3232      #transition probability         p_22 <- 1- p_21     #transition probability          num_samples      <- 5    #number of samples     num_time_periods <- 21   # number of years          sample_u <- array(0,c(num_samples,num_time_periods))          for(i in 2:num_time_periods){         sample_u[,i] <- runif(num_samples, min=0, max=1)}            sample_u          ifelse(sample_u <= p_11, rho_1 <- 1 , rho_1 <- 2 )      I know that I can simplify the transition matrix using the following commands:               trans.mat <- matrix(c( 0.0468 , 0.9532, 0.3232, 0.6768), 2, 2, byrow = TRUE)      How do I simulate a Markov chain from the `sample_u` according to that transition matrix? The condition is: If `u<= p_11`, then `rho_1 =1` (in state 1) , if `u >p_11`, then `rho_1=2` (in state 2). I know that my code is wrong because it is not referring to the previous state. I think I need to correct the `ifelse` command but I'm not sure how to do this. I want to find the state for each period in a matrix form.


title: Standard deviation of recall, precision, accuracy and F-score?

text: I'm not sure how to correctly calculate the **standard deviation** of recall, precision, accuracy and F-score. For example:               TP  FP  FN  TN  total   recall   precision    F-Score    accuracy                                40  10  15  50  115     0.727      0.800       0.762      0.783      TP = True Positive, FN = False Negative, etc. _(Do I use the formula for standard deviation of a binomial distribution: $\sigma = \sqrt{np(1-p)}$? OR $\sigma = \sqrt{\frac{p(1-p)}{n}}$ with $n$ = total and $p$ = recall, precision, F-Score, or accuracy.)_ How can i calculate the standard deviation? I'm not even sure if it makes sense but i was asked to check this.


title: How do I fit a non-linear regression to only one factor level on my scatter plot?

text: My data is structured like so: 'data.frame': 50 obs. of 3 variables:   * Project : Factor w/ 2 levels "A","B": 1 1 1 1 1    * x: int 2 2 2 2 6 4 4 4 6 4 ...    * y: num 0.622 0.425 0.363 0.344 0.346 ... I have 'attached' the data and plotted both levels in my scatter plot using:               plot(x,y,pch=as.numeric(Project))     TSF<-x[order(x)]     SD<-y[order(x)]      And fitted a non-linear regression to the data               nls_fit <- nls(SD ~ a - (b*TSF)+ ((c*TSF)^(2)), start = list(a = 0.34, b = 0.017,      +     c = 0.0003))     lines(TSF, predict(nls_fit), col = "red")      This works well... But how can I fit this equation to only the factor "A" of the data? I'm very new at this so if you have the time an inclination it would be great if you could describe what any code you write does. Thanks Kris


title: Different representations of the one-tailed/one-side hypothesis test

text: > **Possible Duplicate:**   >  What hypothesis for A is better than B? Consider the Pearson correlation estimate $r$. I want to do a hypothesis test. $$\text{Option 1}$$ $H_0: \text{corr}(X,Y) \leq 0$ $H_1: \text{corr}(X,Y) > 0$ $$\text{Option 2}$$ $H_0: \text{corr}(X,Y) = 0$ $H_1: \text{corr}(X,Y) > 0$ * * * I'm wondering if there's any difference between the two in terms of application (and theory)? I want to do $\text{Option 1}$ but the only information I can find is on $\text{Option 2}$.


title: What is a good way to measure how well a set of data fits to a set of functions

text: I got a set of data sets as reference data, each comprised of 14 non- independent variables that is to be used as a basis for developing a method to compute data of that kind with fewer computational resources. So basically a set of points in 14-dimensional space. What I need to do is parametrize my new method so that the output data, again points in 14-dimensional space, fits the reference data as well as possible. I've done a Principal component analysis, which reduced my 14-dimensional point cloud to a set of 5 linear functions that contain 98% of the variance of my reference data. What I'm looking for now is a way to generate some kind of measure that tells me how well the new data fits to the reference data, either using the original data set or the PCA-derived linear functions. And this is where my knowledge ends right now, so I would be grateful for pointers towards good methods to derive such a "how well does the new data fit the reference data" measure of any kind that I can use as benchmark for my parameters.


title: Minimum sample size for PCA when the main goal is to estimate the first or second principal component?

text: If I have a dataset with $n$ observations and $p$ variables (dimensions), and generally $n$ is small (n=12-16), and $p$ may range from small (p = 4-10) or perhaps much larger (p= 30-50). I remember learning in MVA that $n >> p$, but it seems like this may not be so. Note that for my purposes I am rarely interested in any principal components past PC2.   * What are the rules of thumb for minimum sample size when PCA is OK to use, and when it is not?   * Is it ever OK to use the 1st few PCs even if n=p or n less than p?    * Are there any references on this?   * Does it matter if your main goal is to use PC1 and possibly PC2 either:      1. simply graphically, or      2. as synthetic variable then used in regression?


title: Writing up findings from path analysis

text: I have obtained my results from a path analysis (using MPlus) and I´m not sure how go about writing a results section in APA style. The model is showing a nice fit to the data, my hypotheses are supported, and I have a clear idea of what my result mean, but I just don´t know WHAT statistics are conventionally reported in this kind of analyses and HOW to present them. I have found a book, Kline (2010), and a couple of papers, McDonald & Ho (2002), and Boomsma (2000), but again, I feel unsure as to whether these may be good (or practical) references to do what I need. I feel a written example of a results section presenting a path analysis may help me, but I haven´t found any so far. (This work is for my MSc dissertation, so space constrains for writing the results should not be a problem, although I do plan writing a manuscript version to potentially publish it. If you need additional information that may help with your answer please let me know). Could anyone please point me out in the right direction? Thanks!


title: Matrix factorization vs random walk with restart for recommender systems

text: Suppose I want to handle "friend recommendation" problem on a large social network graph. I came across random-walk-with-restart as one technique used. I was thinking of using matrix factorization as another approach. Is there any comparative study between the two approaches?


title: Which loss function for?

text: I want to model some data with a "Boosted Regression Trees" in order to do so I have to define a "loss function". As a response variable I have count data. So I assume that a "poisson distribution". However my response variable are not poisson distributed. I read somewhere that I for count data also a gaussian or laplace distribution could take into consideration. Can someone give me some advice and maybe explain the difference between all three loss functions? Greetings


title: Autoregressive moving average or feed-forward neural network

text: When we have to make a forecast, the books tell us that the main method is the autoregressive moving average model. In my opinion there is another big tool, the feed forward neural network (FFNN). So I think that we could use two main tools:   * Autoregressive moving average   * Feed forward neural network Of course there must be differences, but I am not an expert. Who, having sufficient experience in these two methods, can explain to me the differences between these two methods on making predictions?


title: Predict only the first N principal components in a PCA analysis

text: I'm using R to analyze a very large dataset. I conduct a PCA on one dataset,               PCA <- prcomp(formula = ~., data = train, scale = T, na.action=na.exclude)      and then want to apply the PCA on another dataset,               test_rot <- data.frame(predict(PCA,test,na.action=na.omit))      This works however it is fairly memory intensive and I'm really only interested in the first N principal components (say: 50). Is there a way to only predict the first 50 PCs in that second line of code?


title: Difference between regression and approximation tasks

text: Do you know any references (books/articles) which give a systematic description of the relation between approximation and regression tasks? I know that, although these terms originate from different worlds - functional theory and statistics, they are commonly used interchangeably but I really am interested in precise definitions and differences.


title: Mixed ANOVA: small and unbalanced samples

text: I have to analyze two samples ($n_1=10, n_2=18$) in a design in which there is a between-subject factor (`Groups`: 2 levels) and a within-subjects factor (`Stimuli`: 2 levels). I'm wondering if, in such situations, a mixed (between/within) ANOVA is justified, given the small and unbalanced samples size. The rule of thumb for ANOVA is at least 30 subjects for cell, I know, but often I read papers in which ANOVA is used for smaller samples. **Could anyone suggests me the best practice for this situation? Can such analysis in this context be defended in a publication process?**


title: when to weight frequencies and variables, and how to determine category boundaries

text: Should income categories of a sample intercept survey be weighted to match the income proportions of the population of a city? Also, what guidelines should be used when deciding on category breaks between income groups? When looking at responses to various questions between income groups do the frequency tables for the questions need to be weighted by the proportion of the income group in the overall survey?


title: Discriminative vs. Generative Models

text: This has be asked before, but I still have not grasped it completely. I know that generative models model the feature distribution and that this includes modelling the P(x|y) and P(y), which are not required if we are trying to classify (find P(y|x)). **Question:** Many text books say that it is easier to include features in discriminative models which is rarely explained. Also they mention that discriminative models allow overlapping features (features that are interdependent). Could anybody explain what this means and why is it true, or guide me to the place to read? I see that it is possible to include features in generative models as well, and I can't see why it should be easier or more efficient in discriminative models.


title: Comparing Results in larger vs smaller populations

text: We have ten schools with different-size populations, ranging from 1052 at the smallest school to 4833 at the largest school. Each school has a goal of achieving 60 percent of its enrollees showing improvement on tests. What are the statistical advantages or disadvantages, if any, for the larger or smaller schools?


title: Plotting the coverage of the confidence interval as a function of sample size using Monte Carlo in R

text: I have a distribution $x\sim\text{exponential}(\lambda)$ and I know that $\hat{\lambda}=\bar{x}$ and Fisher's information to be $\lambda^{-2}$. Now, I need to use `R` to generate the data from this distribution with $\hat{\lambda}=1$, and plot the coverage of the CI of the sample as a function of sample size for different sample sizes using 1000 Monte Carlo reps. All I have so far for `R` is:                x=rexp(25,1)     > x      [1] 0.14576099 0.78689809 0.56063641 1.70545364 0.40507393 0.22229612      [7] 2.15959513 0.19580686 2.14139285 0.30020825 1.04155197 0.11201061     [13] 0.11249431 4.12707736 0.35781892 0.81850649 0.06649549 2.35163544     [19] 0.23481295 0.54353838 0.34712200 0.71332001 0.55924596 0.50507862     [25] 0.41171776     > tcrit<-qt(0.975,24)     > tcrit     [1] 2.063899     > mu<-mean(x)     > mu     [1] 0.8370219     > sig<-sqrt(var(x))     > sig     [1] 0.9603002     > hist(x)      I don't know how to create these CI's in `R`.


title: Bayesian fitting - multiplying two probabilities with differing orders of magnitude

text: I am fitting a model to data using Bayesian inference. This is my first time of using this method. My posterior is $P = P_{prior} + P_{photometric} + P_{spectroscopic}$. Value of $P$ is negative and is approaching 0 while it converges. The problem is that $P_{spectroscopic}$ and $P_{photometric}$ are computed from different types of data and values $P_{spectroscopic}$ and $P_{photometric}$ are of different orders of magnitude: $P_{spectroscopic}\sim-100$ and $P_{photometric}\sim-10$ so $P_{spectroscopic}$ takes over the convergence and doesn't affect $P_{photometric}$. My question is: should usually all the probabilities be of the same order of magnitude in such cases (or soon converge to comparable values)? Because now one term dominates and the other is basically ignored and the convergence is awful (or it does not converge at all). What should I do in such cases? [Edit Aug 25]: Corrected the obvious mistake with the addition in the equation for $P$; included latex


title: How to estimate correlation matrix from largest eigenvalues?

text: I'm trying to estimate a correlation matrix from the 5 largest eigenvalues and associated eigenvectors of the sample correlation matrix. My problem is that the output from the following Matlab code results in much higher average correlation (0.8 v 0.4) for the financial data I'm looking at.               [V,D] = eig(InputCorr);     eigvals = diag(D);     eigvals(1:(end-5)) = 0;     eigvals = eigvals*size(InputCorr,1)/sum(eigvals);     BB = bsxfun(@times,V,eigvals')* V';     T = 1 ./ sqrt(diag(BB));     OutputCorr = 0.5*(BB+BB') .* (T*T');       Is such a marked increase in average correlation to be expected or am I doing something wrong? Update: The $A=λ_{1}v_{1}v_{1}' \+ ... + λ_{n}v_{n}v_{n}'$ decomposition suggested by @Jonathan Lisic is very useful. If I truncate it to just the first five terms, then the off-diagonal elements of the matrix are what I expect (similar to the corresponding elements of the original matrix and similar average correlation), but the diagonal elements are well below 1. So the increase in correlation is caused by rescaling the matrix to get the diagonal elements to equal to one. I am actually only interested in the off- diagonal elements, is it kosher to use these directly without doing any rescaling?


title: Allow data to dictate the priors and then run the model using these priors? (e.g., data-driven priors from same data set)

text: It is my understanding that we should not be allowing the same data set we are analyzing to drive/define what the prior distributions look like in a Bayesian analysis. Specifically, it is inappropriate to define prior distributions for a Bayesian analysis based on summary statistics from the same data set that you are then going to use the priors to help fit a model to. Does anyone happen to know of resources that specifically discuss this as being inappropriate? I need some citations for this issue.


title: Are there problems with arbitrary application of bootstrap?

text: Suppose I have a statistics (say a price index) and I want to obtain standard errors for it. I have heard that blind application of bootstrap may not be a good practice. If true 1- What could go wrong if I just apply nonparametric bootstrap and obtain standard errors for my index? 2- If it is OK, are the rsulting standard errors robust to heteroskedasticity, serial correlation etc.? Thank you


title: Open source classification algorithms, preferably in C++

text: I am in search of open source classification algorithms. I am working on a computer vision project that uses classification for scene recognition. I wish to bench test a range of machine learning algorithms for the classification stage. The implementations must be on open source licenses and preferably written in C++. All and any help would be highly appreciated. Thanks in advance.


title: Predicting a baseline response using ARIMA forecasts

text: I have a modelling dilemma. I am creating a model that attempts to predict demand (leads not sales) based upon the correlation to advertising spend. We know that without advertising spend, demand is driven by seasonality. So our models include seasonal factors like month of the year and even day of the week. If I were building a regular linear regression model, I would fit a linear regression model to a training dataset, to get estimates of the coefficients of the seasonal factors and advertising spend to demand. In order to get an estimate of future baseline demand, I would forecast demand using all the coefficients from the model and then I would estimate a baseline by setting adspend equal to zero. For ARIMA models, there are additional factors such as AR and MA terms. Would I estimate my baseline the same way by just setting the coefficient on advertising spend equal to zero? Thanks for any thoughts.


title: How to get the centroid of one cluster?

text: I have obtained several clusters of pedestrian trajectories by **hierarchical clustering**. Now for each cluster, I wish to merge all the pedestrian trajectories into one "centroid trajectory". I am currently implementing the hierarchical clustering with the Python scipy.cluster.hierarchy. Is there any off-the-shelf function to get the centroid of one cluster? If there isn't, general suggestions on merging methods are also greatly welcomed.