title: auto.arima and prediction

text: I'm really new to ARIMA methods and am trying to forecast electricity load. I've integrated: electricity load, temperature, weekday (dummy), public holidays, and school holidays. My model tries to perform a non seasonal ARIMA with linear regression for each hour of the day. Here is my code for an example of one of the 24 hours (6 AM):               # ElecLoad contains hourly loads and other data for 2005 and 2006 (=2*365*24 entries):     # 1. Electricity load in MW     # 2. day of weak: sunday=0, monday=1, etc      # 3. Hour of the day 0 -> 23     # 4. Public Holiday: 1 if Public Holiday, 0 otherwise     # 5. Scool vacation: 1 if no scool     # 6. Temperature in °F          # Create the weak matrix = dumy variables for the weakdays     weakmatrix<-model.matrix(~as.factor(ElecLoad[,2]))     #Remove intercept     weakmatrix<-weakmatrix[,-1]          #Generate FullTable     FullTable<-cbind(load=ElecLoad[,1], weakmatrix, ElecLoad[,4],                      ElecLoad[,3],ElecLoad[,5],ElecLoad[,5]^2, ElecLoad[,6])     colnames(FullTable)<-c("Load","mon","tue","wed","thu","fri","sat",                            "ScoolHol","PubHol","Temp","Temp2","Hour")          #Create the xreg = substed for a specific hour of the day (column 12 = Hour)     xreg<-subset(FullTable[,2:11], FullTable[,12] == 7)          #Create the Load time serie, also a subset of the full table     LoadTs<-ts(subset(FullTable[,1], FullTable[,12] == 7),start=1,frequency=1)          #Launch of auto.arima     ArimaLoad<-auto.arima(LoadTs, xreg=xreg, lambda=0)      When I try to forecast with the same 2 years data as `xreg`, here is my output               plot(forecast(ArimaLoad,xreg=xreg), include=0)      ![enter image description here](http://i.stack.imgur.com/MpNeH.png) While when I try to plot the fitted it looks identical to my original Load               plot(fitted(ArimaLoad))      ![enter image description here](http://i.stack.imgur.com/zsw81.png) I don't understand why the `prediction()` is so much different than the `fitted()` with the same `xreg` matrix. Is this a normal behaviour, how can I improve my model to better fit with the real situation? * * * Thank you so much for your support. I'm not sure I understood everything from what you propose. You mean that I should build a first model to forecast the daily average load (I prefer the average than the sum because due to DST, some days don't have 24 hours...). This model would be deterministic, but I don't see what kind of model you're thinking off? Is a multilinear regression ok? I prefer to consider the log(load) to make the different parameters multiplicative which I think is better fit to the reality. Then I should have 24 hourly models, taking the daily average then split with a sort seasonal effect? Should I use somewhere an ARIMA model? I'm not convince of considering the month as having an effect, in my opinion there is no reason that consumption is more important in January than August except if we consider the Temperature and Holidays effects. The hour of the day is related to the activity that's the reason why I'm considering the specific model for each hour. The same way each day of the weak is different. I've tried a multilinear regression for the same hour (7:00 AM) and the result looks not so bad.               #Create the frame.data     Load<-subset(FullTable[,1], FullTable[,12] == 7)     FullData<-cbind(LogLoad=log(Load), xreg)     FrameData<-data.frame(FullData)          # multilinear regression     mlin<-lm(LogLoad ~ mon+tue+wed+thu+fri+sat+ScoolHol+PubHol+Temp+`Temp2`, FrameData)     plot(exp(mlin$model$LogLoad), type="l",col="blue")     lines(exp(fitted(mlin)), col="red")      ![enter image description here](http://i.stack.imgur.com/MpNeH.png) fitted() in red which is now exactly the same as predict() if I re-use the same data entry (2005-2006) and looks not so far from the original load in blue (no so bad for a simple model). I still don't fully understand why it did not work with ARIMA as it also takes into consideration multilinear regression. Now my "simple" model already takes into account several parameters, like the temperature, the holiday, the school vacations the day of the weak and the hour of the day (local time, not UCT). How can I improve my model further more? How can I make sure that the parameters are invariant? Is there a specific method?


title: What can be inferred from this data?

text: Suppose you ask 100 people, would you purchase this item for $1000? You get responses that fall into one of three categories, Yes, No or Maybe. Assuming you get a non-trivial number of each replies, can it be estimated what price would give an equal number of Yes's and No's (with the rest being Maybe's). How would one model this? Could it be determined with some assumptions? What about with asking the question a second time with a different price? This is hypothetical and not homework, and I can't seem to determine if there is anything that can be inferred.


title: Combining Deterministic and Random Unbiased Estimators

text: I am trying to compute an expectation $E[f(X;\theta,n)]$ where $\theta$ and $n$ are known parameters. I have an easy-to-compute deterministic function $\tilde{f}(\theta,n)$ that provides an approximation to $E[f(X;\theta,n)]$ that gets better as $n\rightarrow\infty$ (but always underestimates the true value). I can also sample directly from $f(X;\theta,n)$ and compute the Monte Carlo estimate $\hat{f}$, which is unbiased but can have a large variance for some $\theta$ and $n$. I am wondering if there are similar problems that might yield methods for estimating $E[f(X;\theta,n)]$ which combine these two estimates. Since we know $\tilde{f}(\theta,n)\rightarrow E[f(X;\theta,n)]$ as $n$ increases, it seems like it might also be possible to compute the estimators for a range of values of $n$ and use the asymptotic behavior to somehow constrain the fitted values. The only approach I can think of is to use a Bayesian estimate of $E[f(X;\theta,n)]$ by taking the sample mean of the Monte Carlo samples and incorporating a prior having mean $\tilde{f}$. This seems like an OK start, but at the very least I'd like to be able to incorporate the knowledge that $E[f]>\tilde{f}$. It also seems quite likely that a more sophisticated approach exists for a similar problem, but I'm just not sure where to start in searching the literature.


title: Regime-switching regression

text: I know that dependence of $x$ and $y$ in my time series data changes in time. Suppose we have monthly time series from 1990 to 2011. If we independently regress $x$ and $y$ for 1990 to 2000 and from 2000 to 2011, the results are quite different. There is low dependence in the former sample and strong in the latter. Although I approximately know the moment of structural shift, I wish to have a statistical model, which allows reveal this moment more precisely.   I know about regime-switching regression models, but they only provide different parameters and do not test the moment of structural shift.   Is there any regime-switching regression model for time series, which helps to identify the moment of time, when switch happens? I know about markov- switching autoregression models (MSVAR), but in my problem I need to apply simple regression without autoregression terms. Are there other more simple and suitable approaches?


title: Cross-correlation of massive arrays

text: I have 16 1D arrays with approximately 10-11 million double-precision elements each. I need to execute a cross-correlation across them, i.e., 1 with 2, 1 with 3, ..., 1 with 16, 2 with 3, 2 with 4, ..., 2 with 16, and so on. This cannot be done efficiently on my MacBook Intel Core 2 duo 2.4 GHz, with 4GB of RAM. My question is, what is the typical approach, if not brute force (faster processor, more RAM) that people use to overcome this problem, or problems like it? I asked this question on Stack Overflow a few minutes ago, but I realized that there might be more people here who have dealt with this particular problem.


title: Resources for machine learning applied to fault detection

text: I am looking for a book or for a collection of examples where machine learning is used for fault detection.   Any suggestions?


title: Finding Abnormally Low Ticket Prices

text: I'm creating a program where my goal is to analyse lists of ticket prices and locate tickets that are priced abnormally low in relation to the spread of prices. My preliminary model does the following to the price list in this order..(examples are at the end of my writeup - they include the modifications made after running through the model)   1. Removes any tickets that are priced over 999.00 . Sometimes tickets will be priced intentionally at 9999.99 to ensure they do not sell - so i don't want to use them in my calculations.   2. Removes any tickets that are priced 1.5 standard deviations above the new mean. I found that unrealistically high priced tickets were unnecessarily influencing the mean.    3. Recalculates the mean and SD and pulls tickets that are 2.9 standard deviations BELOW the mean (I'm calling this the "Price Floor". This results in some success, however it also produces Price Floors that are negative in some cases. I would like to improve my model to be able to take any price list and develop "Price Floors" that are more realistic. My questions are:   1. Am I going down the right path with my current model?   2. Would you change any aspects of my current model?   3. Are there other methods available that may be more effective than standard deviation? Below are two samples of the data I'm working with using the model I described above. In Set 1, the Price Floor is 36.48, thus producing one result of 36.1. In Set 2, the Price Floor is -3.45. While I would not expect a result from Set 2, I feel that my model is inadequate due to the unrealistic Price Floor. Thank you so much for any help you can provide. **Set 1** (Price Floor: 36.48)               [0] => 36.1     [1] => 41.6     [2] => 43.8     [3] => 44.9     [4] => 45     [5] => 46     [6] => 46     [7] => 46     [8] => 46     [9] => 46     [10] => 48.2     [11] => 49.03     [12] => 49.3     [13] => 49.3     [14] => 49.3     [15] => 50.4     [16] => 51.5     [17] => 51.5     [18] => 51.5     [19] => 51.5     [20] => 51.5     [21] => 52.05     [22] => 52.47     [23] => 52.6     [24] => 52.6     [25] => 52.6     [26] => 52.6     [27] => 52.6     [28] => 52.6     [29] => 53.7     [30] => 53.7     [31] => 53.7     [32] => 53.7     [33] => 53.7     [34] => 53.7     [35] => 54.8     [36] => 54.8     [37] => 54.8     [38] => 55.9     [39] => 55.9     [40] => 55.9     [41] => 55.9     [42] => 55.9     [43] => 55.9     [44] => 55.9     [45] => 55.9     [46] => 55.9     [47] => 55.9     [48] => 55.9     [49] => 57     [50] => 57     [51] => 57     [52] => 57     [53] => 57     [54] => 57     [55] => 59.2     [56] => 59.2     [57] => 59.2     [58] => 59.2     [59] => 59.2     [60] => 60.3     [61] => 61.4     [62] => 62.5     [63] => 64     [64] => 64.15     [65] => 64.15     [66] => 64.15     [67] => 65     [68] => 65     [69] => 66     [70] => 67.45     [71] => 67.45      **Set 2** (Price Floor -3.45)               [0] => 46     [1] => 58.1     [2] => 58.38     [3] => 59.2     [4] => 61.4     [5] => 62.5     [6] => 64.7     [7] => 65.25     [8] => 65.25     [9] => 65.25     [10] => 65.25     [11] => 66.35     [12] => 66.9     [13] => 67.99     [14] => 67.99     [15] => 67.99     [16] => 67.99     [17] => 67.99     [18] => 68.55     [19] => 70.2     [20] => 71.85     [21] => 71.85     [22] => 72.4     [23] => 72.4     [24] => 73.5     [25] => 73.5     [26] => 73.5     [27] => 74.83     [28] => 75.7     [29] => 79     [30] => 79     [31] => 79.55     [32] => 80.1     [33] => 84.5     [34] => 84.5     [35] => 84.5     [36] => 85.05     [37] => 87.8     [38] => 89.45     [39] => 93.3     [40] => 95.5     [41] => 97.15     [42] => 98.76     [43] => 98.8     [44] => 99.66     [45] => 106.5     [46] => 108.7     [47] => 112     [48] => 112     [49] => 119.15     [50] => 119.15     [51] => 122.45     [52] => 123     [53] => 124.65     [54] => 125.75     [55] => 125.75     [56] => 128.5     [57] => 128.5     [58] => 133.45     [59] => 134     [60] => 135.65     [61] => 136.2     [62] => 139.5     [63] => 139.5     [64] => 139.5     [65] => 140.6     [66] => 149.95     [67] => 150.5     [68] => 154.9     [69] => 160.4     [70] => 161.5     [71] => 164.25     [72] => 164.8     [73] => 165.9     [74] => 167     [75] => 167     [76] => 167     [77] => 167.55     [78] => 169.75     [79] => 171.4     [80] => 172.5


title: Lagging over a grouped time series

text: I have a few tens of thousands of observations that are in a time series but grouped by locations. For example:               location date     observationA observationB     ---------------------------------------      A       1-2010   22           12      A       2-2010   26           15      A       3-2010   45           16      A       4-2010   46           27      B       1-2010   167          48      B       2-2010   134          56      B       3-2010   201          53      B       4-2010   207          42      I want to see if month _x_ 's `observationA` has any linear relationship with month _x_ +1's `observationB`. I did some research and found a `zoo` function, but it doesn't appear to have a way to limit the lag by group. So if I used zoo and lagged `observationB` by 1 row, I'd end up with the location A's last `observationB` as location B's first `observationB`. I'd rather have the first `observationB` of any location to be `NA` or some other obvious value to indicate "don't touch this row". I guess what I'm getting at is whether there's a built-in way of doing this in R? If not, I imagine I can get this done with a standard loop construct. Or do I even need to manipulate the data?


title: Difference between neural net weight decay and learning rate

text: In the context of neural networks, what is the difference between the learning rate and weight decay?


title: How to run survival analysis on big dataset?

text: I am recently involved in a project that needs to analyze the survival time of objects. Therefore, I plan to use the rms package to build a Cox model. The problem is, since the dataset I have is so big (about 450,000 instances, and each instances is associated with 9 covariables), that the R environment is fail to handle this. Anyone suggestions?


title: Effects of ignoring the sampling strategy

text: Does the use of calibration/post-stratification imply that the sampling strategy/stratification/individual selection probabilities can be ignored altogether? In a survey reweighting exercise, I have population totals for a complete cross-classification of four attributes: country, age, sex, and work status. **EDIT** : The sampling scheme is stratified by country (about the same number of respondents per country), and within each country stratified by age and sex (very roughly uniform). The work status is not part of the sampling stratification. My task is to produce generic weights that can be used for identifying interactions between variables (stratification-response and response- response). Estimation of population means or totals might be another application. (At least that's what I understand. To be honest, I don't really know the application -- how much does it matter?) **Question** : How will the calibrated weights be affected if information on the sampling strategy is added to the calibration procedure, e.g. in the form of prior weights? Is it a problem when for some combinations of categories there are only very few respondents? Question related to the practical exercise: Choice of calibration totals


title: Curing noncoverage with hot-deck imputation?

text: Is it advisable to use hot-deck imputation to allow for poststratification in the presence of empty post-strata? In a survey reweighting exercise, I have population totals for a complete cross-classification of four attributes: country, age, sex, and work status. However, the survey suffers from noncoverage: For a few category combinations, there is no corresponding sample in the survey. This means that simple post- stratification cannot be applied directly. For the missing category combinations, I could use survey records that are "similar" to the missing ones (hot-deck imputation), and then apply simple post-stratification. I could even repeat this $k$ times and average over the $k$ sets of weights. **Question** : Is this considered "better" or "worse" than collapsing post-strata or using GREG estimation? Why? (Pointers to literature are welcome, too.) The following questions are related: Choice of calibration totals, Effects of ignoring the sampling strategy


title: Difference between survival analysis and classification?

text: Since the time variable can be treated as a normal feature in classification, why not using more powerful classification methods (such as, C4.5, SVM) to predict the occurrence of an event? Why lots of people still use the classic but old Cox model? In case of the right-censoring data, since the time would change for an instances, so I think same object with different time values could be treated as two different instances in classification. Is this OK? Is there are some highly-cited paper on this topic? Thank you!


title: PCA loadings interpretation

text: Regarding the final output of PCA:   1. Using rotation, I'll get a loading for my 50 components, so the equation would be $pc_1=0.05*a_1+0.02*a_2+\dots+0.04*a_{50}$. Now, can I derive the value of $pc_1$ by using my last observed return for $a_1,a_2,\dots,a_{50}$?    2. As I looked at my data, I have found at least 10 principal components are required to explain 80% of the variance, so using same logic as above, I'll get $pc_1,pc_2,\dots,pc_{10}$. Now the regression would be $y=w_1*pc_1+w_2*pc_2+...+w_{10}*pc_{10}+e $. In each of them $w_1,w_2,\dots,w_{10}$ can be their respective eigenvalue and the remaining term is an error term?


title: Participant identifier for a field study in a sensitive work environment

text: We are doing an evaluation in different organizations and need to get data from workers about their work behavior. As this information is something that could get workers in trouble, privacy is an issue, so we go for anonymized questionnaires. To connect the data over time, we are looking for a participant identifier. Based on prior experiences, asking for: > Kindly write down your participant code. This code consists of the first > letter of your place of birth, the first letter of your own first name, the > first letter of the first name of your mother, and your own day of birth. has the disadvantage that some participants think that their colleagues or superiors can guess the code and there are always a few participants who do not know one parent. Also, participants make mistakes with this code. So my question is -- is there a way to ask for a code that allows us to connect data over time (i.e., it is unique for a person), yet is also accepted by participants? **Update:** I had a look at the answers here (thank you very much) -- one caveat, we are organizing the evaluation, but there are different researchers (separate institutions) visiting the sites. Meaning that keeping a list of the participants real names and their codes does not work (if one institution kept it, it would create a bottleneck). Ideally the code is something that is solely on the participants side.


title: What are some good exploratory analysis and diagnostic plots for count data?

text: Does anyone know of good reference material on exploratory analysis and diagnostic plots for count data?


title: Summarize time series regression coefficients after giving 1 coefficient a fixed number (from a distribution)

text: In a time series regression model I am working on, our intervention happened just several months after an big influential external event. The full time series model, therefore, underestimates the external event, but puts more weights on our intervention effect. However, we know that this is not quite correct from other researchs and domain experts. So the current remedy is to give an estimated effect to the external event associated with a given distribution, say the Normal distribution. After simulating for $N$ times, we plug this simulated number as a fixed value in the model, the rest of the coefficients are still estimated from the R software (using "arima" function). Therefore they result in $N$ models, or less. Based on this approach, the distribution of some coefficients are not normally distributed. My question is:   1. Should I still choose the model based on the model selection criteria, or should I summarize each coefficient as a distribution?    2. There is a parameter originally not statistically significant, but with the expected sign. I kept it there as a confounding because it could change the estimate of the intervention parameter by 10%. If we summarize coefficients as distributions, the standard deviation of this distribution is much smaller than the standard error of the coefficient in the full model. And the 95% confidence range does not contain zero as well. Therefore how should I treat it for the reporting purpose, which goes back to Question 1, how should I summarize the models due to this method? It will be really appreciated if someone can provide some suggestion.


title: Estimating a sparse inverse covariance matrix with known sparsity

text: The inverse of the covariance matrix for a distribution can be a good value for the mass matrix of a Hamiltonian monte carlo distribution. If the distribution in question is the posterior of a Bayesian graphical model, many or most of the variables will be conditionally independent of each other. Thus, the inverse of the covariance matrix will have zeros for those pairs of variables. You can figure out the sparsity pattern of the matrix by examining the model graph. I would like to efficiently estimate the inverse of the covariance matrix for a distribution given some IID samples from the distribution and assuming that we know its sparsity pattern. It appears this problem is well studied for the more general case where we're trying to infer the sparsity pattern. For example, this paper is well cited and draws upon a number of other papers. However, I don't know of papers that deal with the case where the sparsity pattern is known. Any pointers?


title: Software to compute p-values and Cohen's kappa for 10 by 10 contingency table?

text: I have a 10 by 10 contingency table with most of my cells = 0, and my total cell count = 70. Is there a statistical package or online calculator which I can use to obtain p-values and Cohen's kappa?


title: Need some help clarifying my procedure using kmeans on multiple batches of data

text: I am working on an online image feature recognition program(BOW histograms) that gets objects in a live cam and extracts the SIFT features. After getting a bunch of pictures, I get the kmeans of the batch and create a global cluster. Next, I repeat the step all over and get more pics to create a new separate batch of clusters. Ideally though, I want to limit the number of clusters that I have and clusters all images at once, but this is impossible since there is an infinite amount of new data coming in. I'm not quite sure how to approach this. Here is how my code works right now: > 1.Image is taken from live video feed, once enough pictures are saved, get > kmeans of sift features.(get 200 clusters) > > 2.Repeat step 1, a new batch of live feed pictures, get kmeans again. > Combine the kmeans vectors with the previous kmeans like :[A B],(total 400 > clusters) You can see that this is bad, because I quickly get too much clusters, and each batch of clusters will definitely have overlaps with another batch. What I want: > 1.Image taken from live video feed, once pics are saved, get kmeans(200 > clusters) > > 2.Repeat step 1, get kmeans again, which updates the the previous clusters. > (still 200 clusters, but some clusters have been updated or changed) Nothing that I've seen could accommodate that, unless I'm just not understanding them correctly. many implementations, such as data-streaming, seemingly only allows you to make the cluster in one go. After all the data are streamed, it doesn't seem so easy to update the existing cluster. I would assume though that this should be possible, but none of the code I found seems to back that up. I would appreciate any help on clearing this up. Thanks!


title: Projecting the lifetime value of a subscription customer

text: We are looking for a way to predict the lifetime value of a subscription customer for a new subscription based company. The problem is that we won't know what the lifetime value of a customer is until the end of their lifetime (when a customer cancels their subscription). We want to be able to predict, for the average customer, what their value is based on 2 1/2 months of data. Can someone suggest a simple, but mostly reliable method for predicting this value? I am not a statistician or math guru, so I would appreciate a simple answer in easy to understand terms. For the sake of simplicity, assume that the price a customer pays per month is constant at $20 and the average cancelations are 10% (that is, if 1,000 joined, 100 have canceled their subscription). I found this article: http://lsvp.com/2010/07/19/how-to-estimate-lifetime-value/ It uses this model to estimate a similar business: https://docs.google.com/spreadsheet/ccc?key=0AnVV-5PRmDIKdHcwQWlvd3F3MDFodXNvamcxZjZ3VlE&hl=en#gid=0 I'm tempted to use this model but concerned that 2 1/2 months of data won't work with this model, which uses 12 months of data.


title: SVM mathematical background

text: I try to get the basic understanding behing SVM algorithm, however I have a problem with basic mathematics. I follow the lecture Support Vector Machine. > Suppose the two classes can be separated by a hyperplane:$(w \cdot x) + b = > 0$ Acoording to wikipedia, hyperplane is defined as $n(r-r_0)=0$, does it mean that $b=-w \cdot r_0$? I tried to consider 2-dimensional case, when $w \cdot x +b =0$ is a line, but it's completely doesn't make sense, $b=-w \cdot x$, where $b$ should be a constant, how can I generalize it to a 2 dimensional case. In addition, why $w$ is actually orthogonal to the plus and minus plane? I appreciate your help.


title: Numerical solution of varying coefficients ODE

text: I have a set of observed raw data and use 2nd order ODE to fit the data $$y''+b_1(t)y'+b_0(t)y = 0$$ The $b_1(t)$ and $b_0(t)$ are time-dependent and I use principal differential analysis(PDA) (R-package: **fda** , function: **pda.fd** )to get the estimate of $\hat{b_1}(t)$ and $\hat{b_0}(t)$ . To check the validity of the estimates $\hat{b_1}(t)$ and $\hat{b_0}(t)$, I use collocation method (R-package **bvpSolve** , function: **bvpcol** ) to get the numerical solution of the ODE $$y''+\hat{b_1}(t)y'+\hat{b_0}(t)y = 0$$ Then compare the solution from ODE above with the smoothing curve fitting of the raw data. My question is that my numerical solution from bvpcol can caputure the shape of the fitting curve but not for the value of the function. There are different in term of some constant multiples. ![enter image description here](http://i.stack.imgur.com/NwTy0.jpg) See the figure of my output. The **gray dot** is my raw data, the **red line** is Fourier expansion of the raw data, the **green line** is numerical solution of bvpcol function and the **blue line** the green-line/1.62. We can see the green line can capture the shape but with values that are constant times of fourier expansion. I fit several other data and have similar situation but different constant. I am wondering it is the problem of numerical solution of ODE or some other reasons and how to solve this problem to get a good accordance between numerical solution(green) and true Fourier expansion? Any help and idea is appreciated! Here is a raw data and code: RData is here               library(fda)     library(bvpSolve)     # load the data     load('y.RData')     tvec = 1:length(y)     tvec = (tvec-min(tvec))/(max(tvec)-min(tvec))     # create basis     fbasis = create.fourier.basis(c(0,1),nbasis=47)     bbasis = create.bspline.basis(c(0,1),norder=8,nbasis=47)     bfdPar = fdPar(bbasis)     yfd = smooth.basis(tvec,y,fbasis)$fd         yfdlist = list(yfd)         bwtlist = rep(list(bfdPar),2)         # PDA fit         bwt = pda.fd(yfdlist,bwtlist)$bwtlist     # output of estimated coefficients     beta0.fd<-bwt[[1]]$fd          beta1.fd<-bwt[[2]]$fd      # define the vary-coef function in terms of t     fbeta0<-function(t)eval.fd(t,beta0.fd)     fbeta1<-function(t)eval.fd(t,beta1.fd)     # define 2nd order ODE     fun2 <- function(t,y,pars) {       with(as.list(c(y,pars)),{         beta0 = pars[[1]];         beta1 = pars[[2]];         dy1 = y[2]         dy2 = -beta1(t)*y[2]-beta0(t)*y[1]         return(list(c(dy1,dy2)))       })     }     # BVP      yinit<-c(p1[1],NA)     yend<-c(p1[length(p1)],NA)      t<-seq(tvec[1],tvec[length(tvec)],0.005)     col<-bvpcol(yini=yinit,yend=yend,x=t,func=fun2,parms=c(fbeta0,fbeta1),atol=1e-5,islin=T)     # plot output     plot(col[,1],col[,2],col='green',type='l')     points(tvec,p1,col='darkgray')     lines(yfd,col='red',lwd=2)     lines(col[,1],col[,2],col='green',type='l')     lines(col[,1],col[,2]/1.62,col='blue',type='l',lwd=2,lty=4)     legend('topleft',col=c('green','darkgray','red','blue'),            legend=c('ODE solution','raw data','basis curve fitting','ODE solution/1.62'),lty=1)      **An updated question: are $\hat{b_1}(t)$ and $\hat{b_0}(t)$ unique?** As I change the basis(say bspline,monomial) for expansion of $b_0(t)$ and $b_1(t)$, usually I have very different estimated curves, but the numerical solution of ODE (green line) are very close. That makes me believe different functional form of $\hat{b_1}(t)$ and $\hat{b_0}(t)$ can contribute the same solution. I started a new post with much more detail here.


title: Concepts of mixed effects in statistics and econometrics, how to cope with them?

text: When I search for the definition about fixed-effects, random-effects or mixed- effects model here or elsewhere on the internet, there are a lot of discrepances. My first exposure to linear mixed-effects model was in longitudinal data analysis in Biostatistics. The definition is clear to me that the fixed-effect is the population-averaged effect, and random-effects is the subject-specific effect. Then the mixed-effects model is the model that contains both fixed-effects and random-effects. The mixed-effects model is usually the random effects model because it contains at least one fixed- effects parameters. Like time slope, you have one mean slope for all individuals in the data, and random-effects are those subject-specific slope deviating from the mean slope. However in Econometrics, the deveoplments of fixed-effects and random-effects models have distinct definitions, which is whether heterogeneity correlates or not with the error term. Some statistical tests were developed to test whether fixed-effects or random-effects model should be used. There are lot of social science analyses adopting the Econometric approach as well. Therefore when I read the discussions about fixed-effects, random-effects or mixed-effects models posted by people from different areas, they always confuse me. Even though sometimes the mathematical defintions are similar, the modelling process and consideration behind it are quite different. I hope there are some general discussions on Statistics and Econometrics about their discrepance in defintions or concepts rather than methodologies or algorithms used.


title: Help identifying an error function

text: I've recently come across an old Excel sheet that is used to help assess whether someone has collected enough samples when performing a time study (essentially trying to figure out the distribution of time required to complete a task). In the data column (labeled "time"), the person performing the time study enters the length of time it took to complete a task, with each row representing one trial. The sheet contains the following formula, labeled "Current % Error": \begin{equation}\frac{t \sigma}{\mu \sqrt{n-1}}\end{equation} where $t$ is a t-score which varies based on sample size, $\mu$ is sample average, and $\sigma$ is sample standard deviation. My question is, **what algorithm is this**? I know that we can measure confidence intervals using the similar equation $\mu \pm \frac{t \sigma}{\sqrt{n-1}}$, but I'm not sure what's being calculated when we divide by the mean. What does this correspond to, if anything? ![Figure](http://i.stack.imgur.com/tpk8V.png)


title: How to use prior probability in inferencing from HMM for activity recognition?

text: I am interested in modelling human activities using sensor data with HMMs and would like to incorporate prior knowledge during inference. The normal procedure is to model K different activities with K separate HMMs. To test an unknown sequence, compute its likelihood from each of the HMMs and the HMM with maximum value is assigned as the class label. This is all done under the assumption that priors over HMM are uniform. A problem can arise when one of the class is rare or unusual and its prior probability may be very low in comparison to other classes and therefore the uniform priors may not be a good assumption. Therefore, I am interested in posterior probability and not just the likelihood to capture the combined effect. My observations are continuous (features extracted from sensor) and not discrete values. My questions are:   1. Can inference be done using a bayesian network type approach that include multiplication of prior with likelihood?    2. In my case the prior will be the count of activities available per HMM. Can that be estimated using a dirichlet prior to avoid zero-count problem for rare class (assuming I approximate an HMM for a rare class). Does that make sense?   3. The multivariate observation data is approximated using single gaussian (and not mixtures), in that case likelihood will be gaussian?, can it be mixed with dirichlet prior to compute posterior probability? or the likelihood is still multinomial as it represents K different outcomes from K different HMMs? Sorry if I have mixed with some of the basic concepts, I am new and I seek guidance to move further.


title: How to simulate an unreplicated factorial design?

text: I am dealing with a unreplicated factorial design. I have some illustrative examples but I need to simulate some unreplicated factorial designs. I do not how and what to use. Can $R$ handle this? For example, I would like to analyse a $2^{4}$ factorial design (factors are A, B, C and D) with only one run and 15 contrasts. I have a single column for response. I would like to compare some methods in the literature to see which method detects active effects better. Thus, I set the active effects to have the same magnitude of $1.5\sigma$ and I would like to generate $100$ response vectors using errors that are i.i.d. with $\mathcal N(0 ,1)$. My true model has four active effects and I would like to simulate $100$ response vectors using this true model $y=3+1.5A+1.5B+1.5C+1.5BC$. But I do not know how to generate data like this using R. ![An example from Montgomery for $2^4$ unreplicated factorial design](http://i.stack.imgur.com/e7F8i.png) * * * Thanks gung for your reply. I just wrote a simple code before I saw your answer here. I think, I need to build up a bit more R knowledge. Anyway, here it is: For the analysis of unreplicated factorial designs with $k$ factors and $p=2^{k}-1$ factorial effects (the main effects and interactions), the following model is generally used \begin{equation} y=\sum\limits_{i=0}^{p}x_{i}\beta_{i}+\varepsilon_{i} \end{equation} So, Firstly I introduced my sign table for $2^{4}$ and $\beta$ coefficients of so-called active effects. Sign table consists of rows (runs) and columns (contrasts with general mean). ![enter image description here](http://i.stack.imgur.com/e7F8i.png) And then, I created my regression equation with magnitudes of active effects and zeros of remaining inactive effects. My simulated model, for example, was $y=3+1.5A+1.5B+1.5C+1.5BC$. ![enter image description here](http://i.stack.imgur.com/9N4GA.png) And then, I run the code below               x=read.csv("sign2.txt", header=TRUE)     sign= as.matrix(x)     is.matrix(sign)          y=read.csv("beta2.txt", header=TRUE)     beta= as.matrix(y)     is.matrix(beta)          signt=t(sign)          bs=t(beta %*% signt)          epsilon=matrix( rnorm(16*1,mean=0,sd=1), 16, 1)           response=bs+epsilon      However, unfortunately, it's for one simulation. I will put a loop command to run the simulation n-times.


title: Conditional/normalized multiple imputation

text: Suppose I have a data set with a certain outcome $Y$, covariates $X$, and a certain status variable $Z$, which can take a finite (small) number of values, say 1, 2 and 3. Any of these variables may be missing in the data set, so I want to multiply impute my data. On top of the imputations from the model $Y|\\{X,Z\\}$, I want to obtain normalized imputations $Y|\\{X,Z=1\\}$ -- i.e., force one of my predictor variables to be set to a specific level. The context is somewhat similar to that of the BMI-liars exercise in Sec. 7.3 in Stef van Buuren's FIMD book. The status $Z$ corresponds to different sources of measurements on $Y$, and I suspect that the status $Z=1$ is the most accurate, so I want to get a feeling of what outcomes on $Y$ would have been if everybody were measured using the source $Z=1$. The difference though is that I don't have any parallel measures, like his self-reported and instrumented BMI. So what I need, computationally, is run the burn-ins, calibrate the imputation model(s), and in the last iteration for $Y$, substitute $Z=1$ instead of its actual or predicted values. There may be a way to create a passive variable that is constant $=1$, but then it would be dropped from the imputation equation as collinear with the intercept term. If I just create a copy of $Y$ and make it missing for $Z \neq 1$, and put $Y$ and $Z$ as predictors, then I get a perfect prediction with singular matrices, so that's a no-go, either. Any ideas how that can be implemented using reasonably standard packages? I would like to use Stata or R for this.


title: R: Extremely large bandwidth, npregbw function in np package

text: Using npregbw function in package np, I get the extremely large bandwidth values for local linear smoothing.               library(np)     bandwidth <- npregbw(xdat = xdat, ydat = ydat, regtype="ll", bwmethod="cv.ls")$bw      I have 4 datasets (2 simulations, 2 real examples) and for some of them I get (seems like) the right bandwidth values (0.1 - 0.5) for the others the enormous values (10^5 - 10^6). The sample size of the datasets is the same, the structure very similar. I ran into the same problem as the other users posting the question to this issue. There is no available answer. Here: http://thr3ads.net/r-help/2010/09/1036675-local-linear-and-local-constant- kernel-regression-with-np Or here: https://stat.ethz.ch/pipermail/r-help/2011-September/289538.html Any similar experiences? Any thoughts how to fix it?


title: How to create index from correlated variables?

text: My question is to create an index to measure a manager's earnings forecast ability. I use earnings forecast accuracy (dollars difference) and the earnings forecast horizon (days difference) to capture the ability. However, those two variables condition on each other. How should I weight them and create index? The smaller dollar difference and the longer horizon will reflect higher forecast ability.


title: How to specify the unit information prior

text: Continuing on from this question and this question re BIC and its approximation to the Bayes factor with a unit information prior (Kass & Wasserman, 1995), I'm trying to quantify this relationship as a stepping stone into Bayesian stats. So far, my calculation of the BIC approximation of the Bayes factor (based upon my impression of Wagenmakers 2007) is linearly related to my Bayes factor that is calculated from my interpretation of the unit information prior using the INLA package in R. Good start! However, my BIC Bayes factor is ~ 3 times smaller than the Bayes factor calculated with INLA and I'm not sure why. The prior I've used in the " **inla** " function is N(0, 1/(variance * n)) and this seems to me the likely place where I'm out. I'm not sure how I got the multiply by n in the formula, but it appears to work... roughly. Kass and Wasserman have N(0, variance / n) which when converted to precision would be N(0, n / variance), but this gives me a less good relationship. Help based on other Bayesian packages is also welcome. _EDIT_               *Deleted code, see below answer instead*      _EDIT_ So I'm pretty sure I've figured out the one sample case. I would still appreciate help for the two sample case and the regression case (which I'll start working on now).


title: Cross validation - feature information outside the fold

text: i am building a model that predicts likelihood of the presence of an infectious disease across 10000 villages. There is spatial dimension - each villages is surrounded by e.g. 8 other villages, and presence of disease in a nearby village increased the likelihood of disease in the village in question. I created a feature that sums the number of cases in the 8 nearby villages. However - initially - i created this value on the entire dataset. With regards to e.g. 10 fold cross validation - is this "wrong"? Am i bringing in information from "outside the fold" and hence bias the performance? Should i run the feature creation within each fold, even those some of the nearby villages will be outside the fold? (also - computationally this might become very heavy...) Thoughts welcome - much appreciated


title: Wavelets with forecast model

text: I am new to wavelets decomposition technique and am trying to use it with time series model. What I use now is that   1. use discrete wavelet transform. I use Daubechies with 4 as mother wavelet with 3 levels (obtained 4 waves).   2. reconstruct those wavelets and then use time series model e.g., ARIMA on thoes 4 waves. My question is   1. Do we have to reconstruct wavelets everytime before forecast ?   2. What can we use the coefficient for ?


title: Software implementation of Support Vector Machines

text: I don't have a background of computer science, but i have my major in mathematics. I am interested in Support Vector Machines. I went through the theory as well as some practical examples. My problem is that I understood all the concepts, like separating hyperplanes, kernels and the importance of the parameters C and gamma. But what I would like to know is, if we implement this algorithm of Support Vector Machines in software, then what outputs the software will give us and how do we implement those outputs?. What are the important things to note down?


title: How to test the effect of a grouping variable with a non-linear model?

text: I have a question regarding the use of a grouping variable in a non-linear model. Since the nls() function does not allow for factor variables, I have been struggling to figure out if one can test the effect of a factor on the model fit. I have included an example below where I want to fit a "seasonalized von Bertalanffy" growth model to different growth treatments (most commonly applied to fish growth). I would like to test the effect of the lake where the fish grew as well as the food given (just an artificial example). I am familiar with a workaround to this problem - applying an F-test comparing models fit to pooled data vs. separate fits as outlined by Chen et al. (1992) (ARSS - "Analysis of residual sum of squares"). In other words, for the example below, does the fitting of two models significantly reduce sum of the squared residuals (in this example, yes): ![enter image description here](http://i.stack.imgur.com/ho09y.png) I imagine there is a simpler way to do this in R using nlme(), but I am running into problems. First of all, by using a grouping variable, the degrees of freedom is higher than I obtain with my fitting of separate models. Second, I am unable to nest grouping variables - I don't see where my problem is. Any help using nlme or other methods is greatly appreciated. Below is code for my artificial example:               ###seasonalized von Bertalanffy growth model     soVBGF <- function(S.inf, k, age, age.0, age.s, c){         S.inf * (1-exp(-k*((age-age.0)+(c*sin(2*pi*(age-age.s))/2*pi)-(c*sin(2*pi*(age.0-age.s))/2*pi))))     }          ###Make artificial data     food <- c("corn", "corn", "wheat", "wheat")     lake <- c("king", "queen", "king", "queen")          #cornking, cornqueen, wheatking, wheatqueen     S.inf <- c(140, 140, 130, 130)     k <- c(0.5, 0.6, 0.8, 0.9)     age.0 <- c(-0.1, -0.05, -0.12, -0.052)     age.s <- c(0.5, 0.5, 0.5, 0.5)     cs <- c(0.05, 0.1, 0.05, 0.1)          PARS <- data.frame(food=food, lake=lake, S.inf=S.inf, k=k, age.0=age.0, age.s=age.s, c=cs)          #make data     set.seed(3)     db <- c()     PCH <- NaN*seq(4)     COL <- NaN*seq(4)     for(i in seq(4)){         age <- runif(min=0.2, max=5, 100)         age <- age[order(age)]         size <- soVBGF(PARS$S.inf[i], PARS$k[i], age, PARS$age.0[i], PARS$age.s[i], PARS$c[i]) + rnorm(length(age), sd=3)         PCH[i] <- c(1,2)[which(levels(PARS$food) == PARS$food[i])]         COL[i] <- c(2,3)[which(levels(PARS$lake) == PARS$lake[i])]         db <- rbind(db, data.frame(age=age, size=size, food=PARS$food[i], lake=PARS$lake[i], pch=PCH[i], col=COL[i]))     }          #visualize data     plot(db$size ~ db$age, col=db$col, pch=db$pch)     legend("bottomright", legend=paste(PARS$food, PARS$lake), col=COL, pch=PCH)               ###fit growth model     library(nlme)          starting.values <- c(S.inf=140, k=0.5, c=0.1, age.0=0, age.s=0)          #fit to pooled data ("small model")     fit0 <- nls(size ~ soVBGF(S.inf, k, age, age.0, age.s, c),        data=db,       start=starting.values     )     summary(fit0)          #fit to each lake separatly ("large model")     fit.king <- nls(size ~ soVBGF(S.inf, k, age, age.0, age.s, c),        data=db,       start=starting.values,       subset=db$lake=="king"     )     summary(fit.king)          fit.queen <- nls(size ~ soVBGF(S.inf, k, age, age.0, age.s, c),        data=db,       start=starting.values,       subset=db$lake=="queen"     )     summary(fit.queen)               #analysis of residual sum of squares (F-test)     resid.small <- resid(fit0)     resid.big <- c(resid(fit.king),resid(fit.queen))     df.small <- summary(fit0)$df     df.big <- summary(fit.king)$df+summary(fit.queen)$df          F.value <- ((sum(resid.small^2)-sum(resid.big^2))/(df.big[1]-df.small[1])) / (sum(resid.big^2)/(df.big[2]))     P.value <- pf(F.value , (df.big[1]-df.small[1]), df.big[2], lower.tail = FALSE)     F.value; P.value               ###plot models     plot(db$size ~ db$age, col=db$col, pch=db$pch)     legend("bottomright", legend=paste(PARS$food, PARS$lake), col=COL, pch=PCH)     legend("topleft", legend=c("soVGBF pooled", "soVGBF king", "soVGBF queen"), col=c(1,2,3), lwd=2)          #plot "small" model (pooled data)     tmp <- data.frame(age=seq(min(db$age), max(db$age),,100))     pred <- predict(fit0, tmp)     lines(tmp$age, pred, col=1, lwd=2)          #plot "large" model (seperate fits)     tmp <- data.frame(age=seq(min(db$age), max(db$age),,100), lake="king")     pred <- predict(fit.king, tmp)     lines(tmp$age, pred, col=2, lwd=2)     tmp <- data.frame(age=seq(min(db$age), max(db$age),,100), lake="queen")     pred <- predict(fit.queen, tmp)     lines(tmp$age, pred, col=3, lwd=2)                    ###Can this be done in one step using a grouping variable?     #with "lake" as grouping variable     starting.values <- c(S.inf=140, k=0.5, c=0.1, age.0=0, age.s=0)     fit1 <- nlme(model = size ~ soVBGF(S.inf, k, age, age.0, age.s, c),        data=db,       fixed = S.inf + k + c + age.0 + age.s ~ 1,       group = ~ lake,       start=starting.values     )     summary(fit1)          #similar residuals to the seperatly fitted models     sum(resid(fit.king)^2+resid(fit.queen)^2)     sum(resid(fit1)^2)          #but different degrees of freedom? (10 vs. 21?)     summary(fit.king)$df+summary(fit.queen)$df     AIC(fit1, fit0)               ###I would also like to nest my grouping factors. This doesn't work...     #with "lake" and "food" as grouping variables     starting.values <- c(S.inf=140, k=0.5, c=0.1, age.0=0, age.s=0)     fit2 <- nlme(model = size ~ soVBGF(S.inf, k, age, age.0, age.s, c),        data=db,       fixed = S.inf + k + c + age.0 + age.s ~ 1,       group = ~ lake/food,       start=starting.values     )      Reference: Chen, Y., Jackson, D.A. and Harvey, H.H., 1992. A comparison of von Bertalanffy and polynomial functions in modelling fish growth data. 49, 6: 1228-1235.


title: Mixed Model in SAS

text: I want to analyse the data from an experiment. The participants get different acoustic signals and judge them. The score is a high number, if they feel comfortable with the signal in the other case it is a lower number. I have 5 different signals and 15 participants. The experiment was performed 5 times for each subject, so I've get 65 values. I want to use a mixed model. I get different results from SPSS and SAS and cannot find the error. It could be that I used PROC MIXED wrong.               proc mixed data=ein1;     class Signal ProNr;     model se = Signal;     random ProNr;     repeated /subject=ProNr;     run;     quit;      The result is:                                                   Die Prozedur MIXED                                        Covariance Parameter Estimates                                      Kov.Parm     Subjekt    Schätzwert                                      ProNr                            0                                 Residual     ProNr          112855                                                 Anpassungsstatistiken                                    -2 Res Log-Likelihood           881.1                               AIC (kleiner ist besser)        883.1                               AICC (kleiner ist besser)       883.2                               BIC (kleiner ist besser)        883.7                                             Typ 3 Tests der festen Effekte                                         Zähler            Nenner               Effekt       Freiheitsgrade    Freiheitsgrade    F-Statistik    Pr > F                    Signal                     4                48           1.00    0.4169      Did I choose the right parameters? I think that the result says that Signal isn't significant (0.4169 > 0.05). Is this the right interpretation of the results? I used in SPSS an anova with repeated measurements. SPSS corrected the results an the result with greenhouse geisser correction is 0.013. Why this results are so different? **EDIT1:** The data looks like:               se   Pro_Nr signal     3.4799  1   1     3.3682  2   1     4.3217  3   1     2.9976  4   1     5.5861  5   1     6.7242  6   1     2.6379  7   1     1.9341  8   1      Pro_nr contains values from 1 to 15, every subject get one special number. signal contains numbers from 1 to 5. Each signal gets its own number. **Edit2:** Proc glm isn't the right procedure to analyse this data, because no group variable can be defined. I have to use proc mixed for this problem. **Edit3:** Is it possible to compare the results of an anova and a mixed model? Did I get the same results with both analyses? New code:               PROC MIXED DATA=ein1;     CLASS Frequ ProNr;     MODEL se = Frequ;     REPEATED / subject=ProNr type=cs;     RUN;     QUIT;      Did I formulated the model (and so the program code) right to answer the question, wheater se depends on frequ (under the random and repeated effect of subjects)? Do you have some favourite literature to mixed models?


title: Is there a way to test a probabilistic prediction for an uncontrollable experiment?

text: This is something that I've wondered about mostly in connection with weather forecasting, but I would like to ask the general question because it comes up from time to time in other contexts. Suppose I have a process, denoted $F$, which, given some initial state $x$, produces an outcome $y$ from among a discrete set of possible outcomes: $F(x) \in\\{y_1,\ldots,y_k\\}$. The process is nondeterministic, so that running it multiple times with the same initial state $x$ will not always produce the same result. Also, suppose I have a theoretical model for the process which, given $x$, predicts a probability distribution over the various outcomes $y_1, \ldots, y_k$. I would like to test the hypothesis that the model is "correct." Or more precisely, given the following hypothesis: > For all initial states $x$, the distribution of results produced by the > process converges to the probability distribution predicted by the model as > the number of runs starting from state $x$ approaches infinity. I would like to know at what significance level I can reject this hypothesis given a set of data points, each consisting of an initial state $X_i$ and the corresponding outcome $Y_i$. If all the initial states $X_i$ are the same, then apparently I can do this by using a Pearson's $\chi^2$ test or Fisher's exact test. (Someone please correct me if this is wrong) If there are several _different_ states among the $\\{X_i\\}$, then I can group the data by initial state and run Pearson's or Fisher's test on the group where $X_i=x_1$, and on the group where $X_i=x_2$, etc. all separately. In this way I could determine a significance level for the related hypothesis > For _one specific_ initial state $x_j$, the distribution of results... _(as > above)_ separately for each $x_j$. But I'm not interested in a bunch of individual significance levels for different initial states. I would like to have _one_ statistical test that combines _all_ my data and gives me _one_ significance level for the model as a whole. This is essential because, in the extreme case, all the $X_i$ are distinct, so if I group my data by initial state, each group has only one element - which makes certainly Pearson's test, and to a large extent Fisher's as well, basically useless. So, to sum up: **is there a statistical test that will give me a significance level at which I can reject the hypothesis listed at the top of the question?** If not in general, are there assumptions that make it possible? (e.g. does it matter if $x_j\in\mathbb{R}^n\forall x_j$ for some known, fixed, finite $n$?) Or is it known that this is impossible, and in that case, what's the reason? In case it wasn't obvious, here's the canonical example I had in mind: $x$ would be the current state of the atmosphere (temperature, pressure, humidity, etc.) collected at weather stations around the US, plus satellite data, plus other stuff; the prediction would be something like P(rain)=0.6, P(cloudy)=0.3, P(sunny)=0.1 at some specified place and future time. There are a few previously asked questions I've found that have relevant answers, especially this one and this one. Between them, there are several methods listed for _comparing_ probabilistic predictions:   * Brier score   * Scoring rules   * Logistic regression   * ROC curve But as far as I can tell (and I may well be wrong), all these provide comparative measures, not an actual statistical significance - in other words, the result of a given test (like the Brier score) is only useful for comparing against results of the same test on a different theoretical model. * * * I'm not particularly experienced with statistics, so if there are advanced concepts involved in answering this, it'd be great to have some references for further reading.


title: How to sample from a Wishart distribution?

text: From Wikipedia, we know that $n$, the degrees of freedom, should be larger than $p-1$ where $p$ is the dimension of the scale matrix. Also, from the bottom part of the same article, we see "Bartlett decomposition" to sample Wishart distribution where the diagonal element is generated from $\chi^{2}(n - i + 1)$. However, if you type the following commands in Matlab 2011b:               Sigma = eye(2)     nu = 1     wishrnd(Sigma,nu)      It will return something! Here, `nu` (degrees of freedom) is no longer greater than the dimension of the scale matrix, right? Also, if you look close in to the implementation of `wishrnd` by issuing `type wishrnd`, we see the following line:               a = diag(sqrt(chi2rnd(df-(0:n-1))));      Clearly, the implementation generates $\chi^{2}(n-i)$ instead of $\chi^{2}(n-i+1)$. Something incompatible here?


title: Given a known population size and rate of occurrence, how do I calculate a reasonable sample size?

text: I'm working with automating some data-retrieval, and now need to do some verification of the results. The system has automatically flagged 6% of the records, and I would like to determine a reasonable size of the remaining records to be manually verified, which will give me some assurance that other errors, if present, are identified. Here's the quick numbers:   * 2,502 total records   * 152 (6%) known errors Let's say I want an 80% chance of spotting errors in the remaining ~2,350 records. How many should I request be reviewed? Is there an easy formula to calculate it with different thresholds? * * * Update: Thanks for the input, everyone. My goal is to determine a rough number of "good" records which should be manually reviewed to determine if there are still errors. These are all results which were labeled as "high confidence" by the automated system we used, so in theory we expect them to be 100% accurate. I'm trying to be diligent in not fully trusting the automated process, but also conscientious of the cost of reviewing these. I would estimate that review will take approximately 1 minute per result; so reviewing all 2500 would be 40 hours of labor.


title: Ridge regression in Matlab produces different results than direct calculation

text: I am trying to run Ridge regression in Matlab (`ridge`): $L$ is some matrix, $x$ is some random vector and $y=Lx+\alpha n$ is another vector. I expect $\mbox{ridge}(y,L,\alpha) $ to return the same result as: $(LL^T +\alpha^2 I)^{-1}L^T y$, but the latter is significantly better. I can't understand the problem as I think this is exactly what `ridge()` does. I even tried with $\alpha=1$. Appreciate any help.


title: Runs test for randomness for k elements

text: Say I have a sequence of length $n$ of $k$ different elements, s.t. $n > k, k > 2$. I want to show this sequence is random, based on the number of runs. I found this paper by AM Mood which explains how to find the expected values and variances of the number of runs of each type, but I'm unsure how to show that a sequence is random from this data. I'm new to hypothesis testing and statistics in general, so any help would be appreciated.


title: Choosing attributes for clustering/classification

text: The situation is as follows. There are 400 examples in the training set and 200 discrete classes (each class has two examples). There are a few thousand attributes. When I run dimensionality reduction to 2D or 3D, I would like to see (optimally) 200 clusters of 2 points each - one cluster for each class. However in practice that's not the case. The question is, how do you go about choosing the set of attributes that will give optimal results with clustering/classification (whether with dimensionality reduction or without it).


title: Can I fit a mixed model with subjects that only have 1 observation?

text: I have a very large dataset where I have repeated measurements over time for individual locations. Some locations might have 10 data points and some locations have only 1 data point. I fit a mixed model and use locations as random effects. My question is can I still use the location that only has 1 data point (since you can't make a regression line with just 1 data) or should I exclude those locations?


title: Significance testing of slopes with replicates

text: I've made replicate measurements of a parameter (fluorescence) which is expected to increase with time, and I'm having a hard time understanding how to test for significance of the slope of the parameter vs. time (using a linear least squares model). Intuitively, it makes sense to me to combine all time-fluorescence pairs for each replicate into a grand dataset, and look at the significance of the slope of that fit. However, the intercept of the parameter is expected to differ between replicates due to an instrumental artifact. It's not clear to me whether it is legitimate to combine the data for each replicate due to the expected difference in intercepts. Many thanks for any help. to put this in R terms:               t1 <- 1:10     f1 <- rnorm(10)     t2 <- 1:10     f2 <- rnorm(10)     t3 <- 1:10     f3 <- rnorm(10)      Is it legitimate to look at the significance of the slope of               model <- lm(c(f1, f2, f3) ~ c(t1, t2, t3))      Or should I be doing something different?


title: Steps to figure out a posterior distribution when it might be simple enough to have an analytic form?

text: > **Possible Duplicate:**   >  Steps to figure out a posterior distribution when it might be simple enough > to have an analytic form? I must admit this is an embarrassing question. I do a lot of programming for statistics work, but rarely do I ever deal with analytical formulas for things (it's usually all just based on algorithms for processing the data). In this case, I am trying to compute a Bayesian estimate of some coefficients for an autoregression, with 11 data samples: $$ Y_{i} = \mu + \alpha\cdot{}Y_{i-1} + \epsilon_{i} $$ where $\epsilon_{i}$ is Gaussian with mean 0 and variance $\sigma_{e}^{2}$ The prior distribution on the vector $(\mu, \alpha)^{t}$ is Gaussian with mean $(0,0)$ and a diagonal covariance matrix with diagonal entries equal to $\sigma_{p}^{2}$. Based on the autoregression formula, this means that the distribution of the data points (the $Y_{i}$) are normal with mean $\mu + \alpha\cdot{}Y_{i-1}$ and variance $\sigma_{e}^{2}$. Thus, the density for all of the data points $(Y)$ jointly (assuming independence, which is fine for the program I am writing), would be: $$ p(Y \quad | (\mu, \alpha)^{t}) = \prod_{i=2}^{11}\frac{1}{\sqrt{2\pi\sigma_{e}^{2}}}\exp{\frac{-(Y_{i} - \mu - \alpha\cdot{}Y_{i-1})^{2}}{2\sigma_{e}^{2}}}.$$ By Bayes' theorem, we can take the product of the above density with the prior density, and then we'll just need the normalizing constant. My hunch is that this should work out to be a Gaussian distribution, so we can worry about the normalizing constant at the end rather than explicitly calculating it with integrals over $\mu$ and $\alpha$. This is the part I am having trouble with. How do I compute the multiplication of the prior density (which is multivariate) and this product of univariate data densities? The posterior needs to be purely a density of $\mu$ and $\alpha$, but I cannot see how you'll get that out of such a product. Any pointers are really helpful, even if you just point me in the right direction and then I need to go and do the messy algebra (which is what I've already attempted several times). As a starting point, here is the form of the numerator from Bayes' rule: $$ \frac{1}{(2\pi\sigma_{e}^{2})^{5}\cdot{}2\pi\sigma_{p}^{2}} \exp{\biggl [ \frac{1}{2\sigma_{e}^{2}}\sum_{i=2}^{11}(Y_{i} - \mu - \alpha\cdot{}Y_{i-1})^{2} - \frac{\mu^{2}}{2\sigma_{p}^{2}} - \frac{\alpha^{2}}{2\sigma_{p}^{2}} \biggr ] }.$$ The issue is how to see that this reduces down to a Gaussian density of $(\mu, \alpha)^{t}$.


title: Weighing contrasts differently for unbalanced designs

text: I'm using a mixed effect multilevel model for a crossed design in R. By design, one of my conditions has half the number of observations as the other condition (I am not talking about missing values here). Of interest is a contrast comparing the two. Here's the questions:   * Should I weigh the contrasts differently taking into account the number of observations?   * If yes, should more weight be given to the condition with larger number of observations since it has more bearing on the data or should I do the reverse to balance out the effect of number of observations? Thanks, Bonnie


title: Normal distribution approximation using central limit theorem

text: I have a variable with 150 observations and it is non normal. Wouldn't CLT be useful in this context. What does large n really mean in CLT?


title: Adding between subject factor reduces p value of main effect?

text: I'm investigating the effect of tumor location on mood (depression / anxiety) and also seeing if time has an effect at two intervals (`t0` and `t1`). I had a data set with `Time` as a repeated measurement variable, `Horizontal location of tumor` as an independent variable (categorical) and `Mood score` as a dependent variable. In my first analysis `Time` turned out to be significant so I thought great, I may have got something here. `Horizontal location` was not significant (i.e., regardless of horizontal location of the tumor time still had an effect). Next, I added another between-subject variable, namely `Vertical location of tumor`, to the analysis. However, if I now run this three-way factorial repeated measure design `time` doesn't turn out to be a significant factor on mood anymore! I feel this is weird because it was significant in a two way and now it isn't in a three way? (Horizontal location and vertical location both don't turn out significant.) **So my question boils down to: how come that adding an additional between factor variable makes the main effect of time turn non-significant (increases the p-value)?** Am I allowed to run two two-way interaction with `Horizontal` / `Vertical` location separated instead of one big three-way or does that confound the results?


title: Reversing IV and DV in multiple regression

text: I am interested in a relationship between adolescents' thoughts about cigarettes (e.g., how many times have you thought about cigarettes in the past day?) and cigarette use (e.g., how many cigarettes did you smoke in the past month?), where both measures are continuous. With a longitudinal study (e.g., three-wave study), it would be possible to examine if cigarette thought would predict cigarette use, cigarette use would predict cigarette thoughts, or both (i.e., a reciprocal relationship). However, if I conduct a cross-section study, where I would measure the above two variables and some demographic variables (e.g., gender, ethnicity, etc.), could I do more than a correlational analysis?? For example, would it be reasonable to conduct two multiple regression analyses with cigarette thought as DV, cigarette use, gender, and ethnicity as IVs in one model, and then switch DV and one of IV (i.g., cigarette use <-> cigarette thoughts) in a next model? I am aware that you should stop at a correlational analysis when you do not have a specific hypothesis about two variables; however, I feel that even within a context of a cross-sectional study, it may be fine to run two multiple regression analyses to examine different predictable relationships between two variables as long as you do not talk about a 'causal relationship.' Is this reasoning wrong? Any inputs would be greatly appreciated.


title: Is the union of two reservoir samples a random sample?

text: I have two reservoir samples (with replacement) of a distribution which are obtained independently. If I just produce the union of these two samples, will the resulting sample also be considered as a random sample?


title: Ideas for outputting a prediction equation for Random Forests

text: I've read through the following posts that answered the question I was going to ask: Use Random Forest model to make predictions from sensor data Decision tree for output prediction Here's what I've done so far: I compared Logistic Regression to Random Forests and RF outperformed Logistic. Now the medical researchers I work with want to turn my RF results into a medical diagnostic tool. For example: If you are an Asian Male between 25 and 35, have Vitamin D below xx and Blood Pressure above xx, you have a 76% chance of developing disease xxx. However, RF doesn't lend itself to simple mathematical equations (see above links). So here's my question: what ideas do you all have for using RF to develop a diagnostic tool (without having to export hundreds of trees). Here's a few of my ideas:   1. Use RF for variable selection, then use Logistic (using all possible interactions) to make the diagnostic equation.   2. Somehow aggregate the RF forest into one "mega-tree," that somehow averages the node splits across trees.    3. Similar to #2 and #1, use RF to select variables (say m variables total), then build hundreds of classification trees, all of which uses every m variable, then pick the best single tree.  Any other ideas? Also, doing #1 is easy, but any ideas on how to implement #2 and #3?


title: How do I calculate the Bayes error of a multivariate normal Bayesian classifier?

text: I have a 4 dimensional feature and each of them are independent normal distributions. I want to calculate the bayesian error associated with this classifier. The covariance matrix and the mean have been given for multivariate normal distribution. I understand that I need to compute the error integral for the 4 dimensional space. However When I tried to compute the integral in mathematica, it throws up an error. I was wondering if there are any approximations or a workaround to compute the integral given all the above data?


title: Probability of a user participating in a thread

text: If I am given the count of users participating in a forum thread on a continuous number of days, how can I can find the probability that a number of users will participate in the same thread on another given day in the future?


title: Introduction to Kalman filters

text: What are good introductory books on Kalman filters? I like lots of examples and practical techniques, and less theory.


title: Interpreting Granger causality test's results

text: I'm trying to educate myself on Granger Causality. I've read the posts on this site and several good articles online. I also came across a very helpful tool, the Bivariate Granger Causality - Free Statistics Calculator, that allows you to enter your time series and calculate the Granger Stats. Below, is the output from the sample data included on the site. I have also taken a crack at interpreting the results. My Questions:   * Is my interpretation directionally correct?    * What key insights have I overlooked?    * Also what is the meaning and interpretation of the CCF charts? (I'm assuming CCF is cross correlation.)  Here are the results and plots that I have interpreted:               Summary of computational transaction     Raw Input   view raw input (R code)     Raw Output  view raw output of R engine     Computing time  2 seconds     R Server    'Herman Ole Andreas Wold' @ wold.wessa.net          Granger Causality Test: Y = f(X)     Model   Res.DF  Diff. DF    F   p-value     Complete model  356              Reduced model   357 -1  17.9144959720894    2.94360540545316e-05          Granger Causality Test: X = f(Y)     Model   Res.DF  Diff. DF    F   p-value     Complete model  356              Reduced model   357 -1  0.0929541667364279  0.760632773377753      ![enter image description here](http://i.stack.imgur.com/R6jem.png) ![enter image description here](http://i.stack.imgur.com/WpScL.png) ![enter image description here](http://i.stack.imgur.com/RUG3Z.png) My interpretation:   * Test was based upon 357 data points and was performed with a lag value of 1   * The p-value of 0.0000294 means I can reject the null hypothesis that x does not cause y for the Y = f(x).   * The p-value of .76 allows me to accept the null for X = f(Y)   * The fact that first hypothesis was rejected and second accepted is a good thing    * I'm a little rusty on my F-test so I don't really have anything to say on this for now.   * I'm also not sure how to interpret the CCF graph. I really appreciate it if any of you who are well versed with Granger- causality could let me know if I'm interpeting this correctly and also fill in some of the blanks. Thanks for your help.


title: 52 variables after backward variable selection on logistic regression on 160 variable at beginning, whether it is illusion or good modeling

text: Thank you for you guys effort on building such a nice community and I learned lot by reading and asking. Background about my logistic regression:   1. Target Variable: Binary target about volatility of market. "1" for trading days with top 15% percentile volatility day which indicates the extreme day on market. "0" for rest of normal and kind of calm day for market.    2. Predictor Variable: frequency of 160 trading related key words from famous trading online forum. For example, the word “trader” in 2008-01-01 we found it occurred for 50 times and we found in 2008-01-01 we have totally 500 posts on that forum. Therefore, we kind of normalize the appearance of "trader" by dividing 50 of "trader" by 500 of "number of posts". We got 0.1 here for 2008-01-01. We have 160 words dictionary and the training data is from 2008-01-01 to 2011-12-31. Each row has 160 column with number indicating the frequency of specific words. That is all. Nothing variable of price, volatility here. Kind of pure sentiment data based on frequency counting.   3. The Goal of Modelling: I think it is quite obvious that we would want to see whether we could have the sentiment words to be good predictor on whether tomorrow would be a EXTRAME VOLATILE DAY?   4. Software I use: SAS Enterprese Miner 3 in SAS 9.1.3. Here is definition and details of modeling:   A. I select logistic regression   B. No data partition because I only have 1000days as my training data so I select to use complete data which means no sampling process and only training data.   C. I select the backward as we want to see what are those significant predictor out of 160, here is some shortage if using forward.   D. There is no validation data set so I chose the Crossvalidation-error as my ultimate criteria. Significant level for variable to stay is 0.05.   5. Result: thankful for the sas easy usage, I got the following graph directly from the the result and modeling manager: First one is the result: as it shows how many variable left: ![http://i.stack.imgur.com/jGCsT.jpg](http://i.stack.imgur.com/Av8FL.jpg) wow, so many variables... scaring. Second one about the response chart which indicates how much we improve the prediction compared to pure random pick or guessing about tomorrow, here, in my case, is about 15% chance you guess it right about whether tomorrow is a volatile day. Response chart looks nice and smooth, right? ![http://i.stack.imgur.com/XoCes.jpg](http://i.stack.imgur.com/LE39t.jpg) WORRIING: TOO MANY VARIABLE LEFT HERE. KIND OF SCARING. HOW DO I IMPROVE MODELLING AND INTERPRET THE RESULT?


title: Difficulty in understanding Hidden Markov Model for syntax parsing using Viterbi algorithm

text: I intend to apply Kevin Murphy's Hidden markov model (HMM) toolbox. I have a set of production rules(arbitrary) $A_0 \to AB [p=1]$, $A\to aC [p=1]$, $B\to bbC [p=0.5]$, $B\to b [p=0.5]$ where $A_0$ is the start symbol, A, B, C are the non terminals, and a and b are the terminals and the probabilities are put in brackets. The observation is {aabc}. I do not know how to use the Viterbi parsing algorithm and cannot understand what all values in the parameters to use in the toolbox. Can anyone explain with a simple example how to use the program so as to select a maximum likelihood parse and how to train? The objective is to classify a language or a text string. Any other example apart from language would also help. Primary question is how to feed in all these data in the toolbox so as to train the HMM and how to use it for the task since evrything in it uses rand function, how do I customize it?   1. Number of states Q=4 ? ($A_0$, A, B, C)   2. Number of observation symbols of length T=4; number of sequences nex=1?


title: R - Confused on Residual Terminology

text: * Root mean square error   * residual sum of squares   * residual standard error   * mean squared error   * test error I thought I used to understand these terms but the more I do statistic problems the more I have gotten myself confused where I second guess myself. I would like some re-assurance & a concrete example I can find the equations easily enough online but I am having trouble getting a 'explain like I'm 5' explanation of these terms so I can crystallize in my head the differences and how one leads to another. If anyone can take this code below and point out how I would calculate each one of these terms I would appreciate it. R code would be great.. Using this example below:               summary(lm(mpg~hp, data=mtcars))      Show me in R code how to find:               rmse = ____     rss = ____     residual_standard_error = ______  # i know its there but need understanding     mean_squared_error = _______     test_error = ________      Bonus points for explaining like i'm 5 the differences/similarities between these. example:               rmse = squareroot(mss)


title: Are these participants repeatedly measured?

text: I have Likert items for a number of students learning Spanish and English at the same time. I really need to compare the means for each (non-normal) variable to see which language is voted higher for variables such as "I watch TV in this language". Now, these participants are clearly not matched pairs, they are the same people being tested under two conditions (Spanish or English). Repeated measures refers to time differences. Is the time idea essential or can I call my participants repeated measures? If not, what are they, and what sort of test is indicated for comparing the means? To anyone who responds, many thanks - I am a newbie


title: Limiting distribution of the first order statistic of a general distribution

text: Let $Z_i,Z_2,\ldots$ be IID Random Variables with density $f$. Suppose that $P(Z_i>0)=1$ and that $\lambda=\lim_{x \to 0+} f(x)>0$. How can I show that $X_n=n \times \min\\{Z_i\\}$ has a limiting exponential distribution with mean$1/\lambda$? I know that the CDF of the first order Statistic is for $t>0$ $$F_{min \\{ Z_i \\}} (t)= 1- \left[ 1-F(t) \right]^{n}$$ Convergence in distribution requires that $F_n (t) \to F(t)$ for some distribution function $F$ but I cannot quite take the limit of the above, with $t$ replaced by $x/n$, as the distribution is not known. I guess I have to insert the second piece of information somewhere but I do not recognise where. Could you please give me a hint or two? Thank you.


title: GSVM-RU Implementation

text: Are there any implementations (free or otherwise) of GSVM-RU (granular support vector machines - repetitive under-sampling) as described by SVMs Modeling for Highly Imbalanced Classification ?


title: What kind of analysis gives you the statement "If you DONT reach X amount by time T, then your chances go down by P percentage?

text: I am trying to model growth for data I have regarding downloads of applications. I would like to make a statement, if you "DONT reach X amount of downloads by time T, then your chances of reaching 15 Million users go down by 75%" Is there some sort of analysis we can do for cutoffs? Thanks!


title: What is the intuition of invertible process in time series?

text: I'm reading a book on time series and I started scratching my head in the following part: ![enter image description here](http://i.stack.imgur.com/iLgJx.png) ![enter image description here](http://i.stack.imgur.com/IkKZF.png) Could someone explain the intuition for me? I couldn't get it from this text. Why do we need the process to be invertible? What is the big picture here? Thank you for any help. I'm new on this stuff so if you could be kind to use student-level terms when explaining this :)


title: Factor intercorrelations and model design

text: Having done an exploratory factor analysis (SPSS), what exactly can I infer from the factor intercorrelation table that is output after the factors are extracted/rotated? If two factors are relatively highly correlated (let's say the .4-.6 range or higher), does this mean that when constructing my model graphically, I could put a bidirectional line between the two latent variables representing the factors? Is that the sort of thing the intercorrelation table represents? If not, what part of the output/what further analysis are you supposed to use to justify putting those bi- or unidireectional lines between factors? I am a newbie on a deadline so help would be greatly appreciated!


title: How to calculate the HHG (Heller Heller Gorfine) distance in R

text: In looking at this question and investigating some recent developments in measuring correlation, I came across the HHG (Heller Heller Gorfine) test. Heller et al. promote it as superior to the recently MIC test that made a splash in 2011. I've found their `R` code and a variant but I can't see how to get a distance measure that ranges between zero and one. How can I calculate that distance? Here's a reproducible example:               ## get the code     download.file("http://www.math.tau.ac.il/~ruheller/Software/HHG2x2_0.1-1.tar.gz", "HHG2x2_0.1-1.tar.gz")     install.packages("HHG2x2_0.1-1.tar.gz", repos = NULL, type="source")     library(HHG2x2)     writeChar(con="myHHG.R", getURL("https://raw.github.com/andrewdyates/HHG_R/master/R/myHHG.R", ssl.verifypeer = FALSE))     source("myHHG.R")          ## built-in example     X = datagenCircle(50)     Dx = as.matrix(dist((X[1,]),diag=TRUE,upper=TRUE))     Dy = as.matrix(dist((X[2,]),diag=TRUE,upper=TRUE))     myHHG(Dx,Dy); pvHHG(Dx,Dy)      Which returns:               $sum_chisquared     [1] 5515.762          $sum_lr     [1] 3028.207          $max_chisquared     [1] 20.02963          $max_lr     [1] 11.30457          $pv     [1] 0.00229954          $output_monte     [1] 10001          $A_threshold     [1] 2.985682          $B_threshold     [1] -4.553877      How can I get a distance metric from this package? Or have I misunderstood what the HHG is about? This seems to be the main source for the HHG test: Heller, R., Y. Heller, et al. (2012). "A consistent multivariate test of association based on ranks of distances." Biometrika. 99 (4) doi:10.1093/biomet/ass070 http://arxiv.org/abs/1201.3522


title: CCA/KCCA for more than two views

text: Canonical Correlation Analysis (CCA) _(and itskernel equivalent (KCCA))_ can be used to find linear _(nonlinear)_ relationships between two aligned multivariate datasets (or views). Is there a way to extend this to more than two views? I suppose one way would be to apply CCA/KCCA recursively in a tree structure, but this appears to be rather inefficient. Is there a single optimisation that can achieve this in one step? Or are there any alternative methods that do the same thing?


title: Data on institutional quality/strength?

text: I'd like to perform some statistical methods on data describing the quality/strength of institutions for countries in the world. I've been searching for this data on the World Bank, OECD, United Nations and IMF. I must have missed this index, as I've heard about people using it from a respectable source like those mentioned previously. Does anyone know about such index from a reliable source?


title: How does a Support Vector Machine (SVM) work?

text: How does a Support Vector Machine (SVM) work, and what differentiates it from other linear classifiers, such as the Linear Perceptron, Linear Discriminant Analysis, or Logistic Regression? * (* _I'm thinking in terms of the underlying motivations for the algorithm, optimisation strategies, generalisation capabilities, and run-time complexity_ )


title: Kernel matrix normalisation

text: Originally posted in stats.SE but never got an answer so reposting here. Is it ever a _bad_ idea to normalise the kernel matrix? By this I mean the method described on page 113 of Shawe-Taylor & Cristianini's "Kernel Methods for Pattern Analysis" (matlab notation):               % original kernel matrix stored in variable K     % output uses the same variable K     % D is a diagonal matrix storing the inverse of the norms     D = diag(1./sqrt(diag(K)));     K = D * K * D;      I'm specifically thinking of text mining, where the kernel is formed from the $tf-idf$ matrix $T$ ( _i.e._ $K = T T'$). The reason I'd like to normalise is that I'm thinking of combining several of such kernels ( _e.g._ from bigrams, trigrams etc.) but I'm getting some really weird results when I do. Summing the kernels without normalisation results in one or the other dominating. (I could use a weighted sum - such as using $1/||K||_F$ where $||\cdot||_F$ is the Frobenius norm but this seems a bit of a hack - this is unsupervised learning).


title: What's a good approach to estimate the probability of word frequencies?

text: I have a document corpus and I want to estimate the probability of occurrence of a certain word $w$. Simply calculating the frequencies and use such a number as an estimation is not a good choice. Is there any work on this topic describing a better approach?


title: How to detect variables discriminating a sample from the rest of the samples?

text: I have a traditionally structured data set where rows are observations and columns are variables. There are only a few observations but comparatively more variables. The observations are regions of a country and the variables are characteristics ( _i.e._ , population size, mean income, number of males, number of females, etc.). I need to find out which characteristics can distinguish one region from the rest. For example, I need to find out if population size distinguishes region 1 most from regions 2, 3, 4, etc. Which statistical methods are available to perform this type of analysis? Your help is appreciated. Just in case it is not clear, my data looks something like the following.                    A, B, C, D     R1   8, 9, 5, 4     R2   5, 8, 4, 5     R3   7, 9, 7, 4      Where R1 = region 1, R2 = region 2, R3 = region 3, A = population, B = average income, C = number of males, D = number of females. My queries are something like this:   * how is R1 different from the set {R2, R3}?    * how is R2 different from the set {R1, R3}?    * how is R3 different from {R1, R2}?  The expected answers look something like this:   * A make R1 different from {R2, R3}?    * {A, D} make R2 different from {R1, R3}?


title: Interpreting a regression modeled on twice differenced data

text: I have built a OLS model with data that was twice differenced. As I understand (and maybe I'm wrong) the coefficients (betas) can be applied to the original undifferenced data to provide Y at that time.... However, when I apply my model to the real undifferenced data, I'm way off and I suspect it's the intercept, or something


title: I want to compare two treatment that were given at 5 time-points

text: I have two groups. One received treatment A and a second received treatment B. The score (dependent variable) was taken at 5 different time points. Can I run a repeated measures ANOVA ? I am only interested if to know if treatment A is different (better or not) than treatment B


title: Learning on huge datasets

text: Basically, there are two common ways to learn against huge datasets (when you're confronted by time/space restrictions):   1. Cheating :) - use just a "manageable" subset for training. The loss of accuracy may be negligible because of the law of diminishing returns - the predictive performance of the model often flattens out long before all the training data is incorporated into it.   2. Parallel computing - split the problem into smaller parts and solve each one on a separate machine/processor. You need a parallel version of the algorithm though, but good news is that a lot of common algorithms are naturally parallel: nearest-neighbor, decision trees, etc. Are there other methods? Is there any rule of thumb when to use each? What are the drawbacks of each approach?


title: Iris dataset how to solve

text: I am new to Machine Learning. For Iris data set problem we can solve the problem with Multinomial logistic regression and as well as neural network. Which would give better performance regarding with cost and error?


title: F-measure for document clustering evaluation - NaN

text: I'm developing the Java application for text document clustering, and I'm researching some evaluation methods. I implemented F-measure (http://en.wikipedia.org/wiki/F1_score), but I have a problem - the returned value is NaN. It happens where a cluster doesn't contain any data from a specific category - precision and recall are equal to zero. How should I handle this situation - F-measure in that case should be zero as well? I will be very grateful for any advice.


title: Population and sample with very different means: why is that so?

text: ![http://postimg.org/image/vppji4o3r/](http://i.stack.imgur.com/WbjNG.png) I have a problem with determining the effect of sample size on the shape of a prior probability distribution. The parameter values are on the $x$ and the prob on the $y$ axis. The right curve consists of individuals (smaller sample size), while the left consists of colonies (with larger sample size). I just don't know why the position of the mean changes. I expected that the larger sample should have peak more (central limit) but not that the mean would shift. Could somebody tell me whether there is any mechanism which would result in the behaviour seen on the picture?


title: Which test to use to compare vitamin D levels

text: I'm working on research into the vitamin D levels of professional and amateur soccer players, and relationship with factors such as skin color, vitamin D-intake etc. I don't know which test to use. I'll give you an idea of the situation:   * I have a continuous variable: the measured vitamin D levels   * A categorial variable: the skin color (either black or white ) I've read about the One-way ANOVA, but it requires three options for the categorical variable. I used the ANOVA test, the p-value was 0.005. Is it possible to use this test with just two options for the categorical variable? Or is there another test I could use? And if it's possible to use the ANOVA, does the value F mean anything, or just the p-value?


title: Theoretical upper bounds of classification accuracy?

text: I'm looking for theoretical upper bounds of classification accuracy. Please let me know if you are familiar with results like the following. **The setup below is a general one, but please share results you know that apply in particular settings as well. Thanks!** Let $Y$ be the binary variable to classify and $X = (X_1,...,X_K)$ be $K$ explanatory variables. We know the joint distribution (in population) of $(y,X) \sim \mu$ and want to build a model of $\hat{y} = f(X)$. Then for the class of XXX functions $f()$, there will be an upper bound to the maximum classification accuracy: $P(\hat{y} = y) \le g(\mu(y,X))$?


title: Inverse transformation sampling for mixture distribution of two normal distributions

text: I am confused by the special way required to use inverse method in the following problem, Here is the problem: > Consider a mixture distribution of two normal distributions, where the > desired PDF $f(x)$ is given by: > > $f(x) = r\, f_a(x) + (1 − r)\, f_b(x)$,where $f_a$ and $f_b$ are normal PDFs > with means $a$ and $b$, respectively (standard deviation is 1 for both). > Using two uniform random variables $u_1$ and $u_2$, explain how we can use > the inversion method to sample from $f(x)$. Note, the `qnorm` command in R > may be helpful here. My confusion is from "two uniform random variables $u_1$ and $u_2$". My thought is that we find out the cdf, $F(x)$ (which can be obtained via `pnorm()` in R), and then we can use some numerical method (such as Newton- Raphson) to generate $x\sim f(x)$, so here it only needs one uniform distribution and does not need `qnorm()`. What's wrong with my method? Does the problem suggest a better method?


title: How does a generalized linear mixed model estimate means and how does this differ from calculating means by hand?

text: I've recently read a paper which used generalised linear mixed models to estimate mean annual and monthly values for the response variable in the model. The response variable was a normally distributed continuous variable and year and month were used as predictors in the model. The annual and monthly mean were then calculated for the response variable.   1. How is the annual and monthly mean estimated for the response variable in a generalised linear mixed model?    2. How would the mean calculated by a generalised linear model differ to a mean calculated by hand?


title: How to combine features extracted by PCA, LDA and LBP?

text: What I'm thinking is to combine PCA features, LDA features and LBP features together to get a higher accuracy, since I think the three features are all kind of histogram vectors and when we decide the label of a test image, we compare the histogram of the test image to the histograms of existing images. So is it possible to combine the feature vectors together? However I think the feature histogram of LBP is very long. If we use 8-by-8 grids then the histogram will be $8\times 8\times 256 = 16384$ long. And the eigenface weight vectors are of the length of how many eigenface components we choose, say 80. So I'm not sure whether we can do this. I'll keep thinking on this and any help is highly appreciated.


title: What is the advantage of median polish over the median?

text: I still can't figure out the advantage of median polish over a regular median for the summarisation of probe to make probesets in microarray analysis. As I understand from here the idea is to use the median of the different arrays to summarise the data along with the median of the different probes: The use of median polish for feature selection I've been making some simulations to try to see what might happen. If we imagine the following is a simulation microarray dataset, where the columns are the probes we wish to summarise, and the rows are the different microarray chips:               library(preprocessCore)     >  y <- matrix(10+rnorm(100),20,5)     > y                [,1]      [,2]      [,3]      [,4]      [,5]      [1,]  9.334358  9.993648  8.551274 10.109988 10.243317      [2,] 11.448786  8.908376 11.536720 10.679236 10.831209      [3,] 11.979813  8.726539  9.740086 10.103683  9.349783      [4,] 10.552108 12.772855 10.484486  9.362849  9.426693      [5,] 10.056581  9.890734  9.907472 10.283063  9.909602      [6,] 12.187766 10.290644  8.770036 11.241425 12.856710      [7,] 11.071675  9.932000 11.761954 10.806470 11.013961      [8,]  9.560737  9.234133 11.307681  8.672639  9.570637      [9,]  8.952978  8.549438  9.962865  8.527808  8.421271     [10,]  8.853584 10.117102 10.040929  9.551693  9.880730     [11,]  8.794862 11.276158  7.579099  9.167762  8.863161     [12,] 10.923659  9.873338  9.081718  9.501927 11.956930     [13,] 10.150289  8.472951  9.367948  9.376648  9.963847     [14,] 10.271810 10.738851 11.253192  8.169423 10.286973     [15,] 10.424398 10.356835 10.004031 10.790331  9.922300     [16,] 10.494939  7.793422 10.182820  9.597307  8.726760     [17,] 11.447273 10.366120 11.620400  9.698011 10.706059     [18,]  9.597951 11.161659  9.923795  8.462029 10.945888     [19,]  8.796656 10.962160 10.204844  7.944251  9.332819     [20,]  9.462135 10.621933  9.697430 11.112817  9.905340      So we can summarise the columns ussing the median (which treats each probeset (that is, column) separately, and also using median polish (which takes into account the columns and the rows, more detail in http://www.stats.ox.ac.uk/pub/MASS4/VR4stat.pdf               > colSummarizeMedian(y)     $Estimates     [1] 10.211050 10.055375  9.983448  9.574500  9.915951          $StdErrors     [1] NA NA NA NA NA          > colSummarizeMedianpolish(y)     $Estimates     [1] 9.975078 9.898037 9.907471 9.788084 9.907471          $StdErrors     [1] NA NA NA NA NA      Then I add an "outlier" by changing one of the arrays and recalculate. As we would expect the median stays the same, but the median polish doesn't. However, I don't see why the median polish is a better answer here than the median answer? I don't see why we are borrowing information from the other columns, and what we gain from this. It seems especially true since most of the time some of the columns will be chips from differential treatments. Perhaps I have missed something rather fundamental here, but any help would be much appreciated, I am at somewhat of a mental block. Also, if people think this is better suited to a more bioinformatics q+a site, please let me know!               > y[20,] <- y[20,] + rnorm(5, sd=5)     > colSummarizeMedian(y)     $Estimates     [1] 10.348104  9.990152  9.983448  9.526810  9.915951          $StdErrors     [1] NA NA NA NA NA          > colSummarizeMedianpolish(y)     $Estimates     [1] 10.030423  9.879213  9.933043  9.713632  9.893364          $StdErrors     [1] NA NA NA NA NA