title: Relative Variable Importance in Clustering

text: In SPSS, the user can check the relative variable importance in a clustering result and produce a graph like the following: ![enter image description here](http://i.stack.imgur.com/f4UUd.png) link: http://pic.dhe.ibm.com/infocenter/spssstat/v20r0m0/index.jsp?topic=%2Fcom.ibm.spss.statistics.help%2Fidh_twostep_main.htm Then, we can identify variables dominating predictor importance or having the most impact at determining clusters. Does anyone know how importance is computed here? I would like to implement this metric in R, if it is not already there.


title: Is there a standard on-disk representation for Markov Chain models?

text: I am doing some experiments with Markov Chain models and want to know whether there is a standard representation for such models on disk in academia/industry? I want to avoid re-inventing the wheel with regards to parsing this data. Hopefully I can also validate my implementation against other implementations for the same model.


title: Get a desired percentage of censored observations in a simulation of Cox PH Model

text: This question comes from: How to simulate a Cox proportional hazards model with change point and code it in R (See answer) I want to generate a censoring variable $C=Exponential(\theta)$ that create a desired percentage of censorship on $Y$, like 33%, 55%, etc, by choosing an appropriate value for $\theta$. In order to include this censoring variable in the code I did for the previous question, I added this to the dataframe generator:               cen <- rexp(n,theta) # Censoring variable.          ycen <- pmin(Y, cen) # Censoring effect.     di <- as.numeric(Y <= cen)      ...               if (surv.df) data.frame(Y,X,ycen,di) else cbind(Y,X,ycen,di)      To do something like this (in the context of the aforementioned simulation):               coxph(Surv(ycen, di)~X)      Which $\theta$ value should I choose?


title: What is the appropriate statistical test for Presence/Absence over time

text: I am assessing the presence or absence of a number of developmental milestones in mice in their first $21$ days of life. Each day, from the day they are born until they have $21$ days, I test if pups have or not specific developmental milestone. This analysis involves $2$ groups, a control and an experimental. I want to test whether the experimental group are delayed in accomplishing these milestones compared to controls. For each test, I end up with a range of days which the presence was observed. For controls, my $n=17$, for experimental my $n=15$. So, for a particular test (for example, cliff aversion), I have observed the presence as early as day $2$ and latest as day $9$. Not sure if a Chi-square would be appropriate and if so, do I have to do a chi-square test for each day comparing both groups? Could I also perform any type of "survival" analysis that allows me to estimate the probability (what day) for each group to have the presence of a specific milestone?


title: Use of 3 different methods - is it the right approach?

text: I am writing a master's thesis on crowdfunding. For this research I set N independent variables and 3 DVs. The reason for doing so is that I want to explore the phenomenon from the aspects of all 3 DVs, as they are giving me different insights, and to crosscheck the results from all 3 models. Here is a bit more details on the three DVs and the models:   1. `Success` [dummy] - I will use logistic regression for that.    2. `Success Ratio` - Ratio variable [0-N] - I use Linear Regression and to linearize the model I use the following transformations:   $\ln(x+1)$ to keep the zeroes, as the distributions are heavily skewed. ==> $\ln(y+1)=b_0+\ln(x_1+1)b_1$..., etc.   3. `N_Backers` - Count variable [0-N] - I plan to use a Poisson GLM with log function, but it might be reasonable as well to be inline with DV2 and to use Linear Regression with transformed variables $\ln(y+1)=b_0+\ln(x_1+1)b_1$..., etc. **Question:** What is your advice on these methods, do they sound reasonable? Further, should I use linear or Poisson for the 3rd DV? My concerns are that I am not sure how the results will be interpreted using all different models. Probably the use of consistent models is the right approach.


title: Intuition for cumulative hazard function (survival analysis)

text: I'm trying to get intuition for each of the main functions in actuarial science (specifically for the Cox Proportional Hazards Model). Here's what I have so far:   * $f(x)$: starting at the start time, the probability distribution of when you will die.   * $F(x)$: just the cumulative distribution. At time $T$, what % of the population will be dead?   * $S(x)$: $1-F(x)$. At time $T$, what % of the population will be alive?   * $h(x)$: hazard function. At a given time $T$, of the people still alive, this can be used to estimate how many people will die in the next time interval, or if interval->0, 'instantaneous' death probability.   * $H(x)$: cumulative hazard. No idea. What's the idea behind combining hazard values, especially when they are continuous? If we use a discrete example with death rates across four seasons, and the hazard function is as follows:   * Starting at Spring, everyone is alive, and 20% will die   * Now in Summer, of those remaining, 50% will die   * Now in Fall, of those remaining, 75% will die   * Final season is Winter. Of those remaining, 100% will die Then the cumulative hazard is 20%, 70%, 145%, 245%?? What does that mean, and why is this useful?


title: IRT model and ceiling effect

text: I did some experiment in which tests are taken twice, pretest and posttest. I found there might be ceiling effect because the average of posttest is close to maximum test score possible. If I assume an IRT model: as the ability is getting higher (above difficulty level of the problems), the expected score distribution is skewed and never goes over the maximum score possible. So I think there might be some way to utilize the ceiling effect and skewedness of score distribution when comparing two averages of the test score. But I wonder if there's any research done already about this subject, which could be called as "comparing two groups' averages assuming IRT model considering ceiling effect"...I am thinking of simulation research possibly with some MCMC.... Any idea or advice would be welcome also!


title: Two treatments - How do I compute reliability estimates like Cronbach $\alpha$?

text: I randomized my survey, and now I have my data on SPSS (and all cleaned up) but when trying to run my Cronbach's alpha I have a problem. Because I have two treatments (sexualised and normal) all my measures are present twice in my data set. For example, “similarity” has 3 items but now I have 6 variables: 3 for the sexualized and 3 for the normal condition.               Sim1       Sim2       Sim3       aSim1       aSim2       aSim3        SPSS will not allow me to combine 2 treatments because I have no data that has answers for the 6 (obviously). **So do I calculate them separately and then get the average between the two, or do I report both Cronbach's for the 2 treatments?**


title: Exclude extraneous values in a range

text: I have a range of readings from a test (for example:0.796, 0.109, 0.11, 0.11, 0.109, 0.109, 0.109, 1.78). I want to remove the extraneous values of these so that I get a more realistic view on the more accurate test result range. If I take the average (1.78) or the median (1.09) it doesn't give me that picture. I am not able to figure out how I could mathematically have a way to achieve this, because I don't want these extraneous values to skew my results.


title: Are there useful applications of SVD that use only the smallest singular values?

text: In a number of singular value decomposition (SVD) applications, for example Latent Semantic Indexing, only the biggest singular values are used to make searches and calculate distances. Are there useful applications that drop the biggest singular values and use only the smallest ones?


title: Degenerate random variable

text: Let $X$ and $Y$ be independent $rv$ such that $XY$ is a degenerate $rv$. Can I say that individually $X$ and $Y$ are also degenerate? Why? thanks..


title: How accurate is IQR for detecting outliers

text: I'm writing a script that analyses run times of processes. I am not sure of their distribution but I want to know if a process runs "too long". So far I've been using 3 standard deviations of the last run times (n>30), but I was told that this does not provide anything useful if the data is not normal (which it does not appear to be). I found another outlier test that states: > Find the inter quartile range, which is IQR = Q3 - Q1, where Q3 is the third > quartile and Q1 is the first quartile. Then find these two numbers: > > a) Q1 - 1.5*IQR b) Q3 + 1.5*IQR > > The point is an outlier if < a or > b My data tends to be things like 2sec, 3sec, 2sec, 5sec, 300sec, 4sec, .... where 300sec is obviously an outlier. Which method is better? The IQR method or the std deviation method?


title: Help with repeated measures regression (in R)

text: I have some small datasets (example data given below) where I am interested in describing relationships between $Y_1$ and $Y_2$ across a gradient $X$ at two time points (date = "early" & "late"). As you can see from the first plot below, with $Y_1$, there appears to be a general increasing trend over across $X$ with the early data (black points), but perhaps not with the late data (brown dots). With $Y_2$, there doesn't appear to be any relationship with $X$; however, there does seem to be a general decrease in $Y_2$ in the late data relative to the early data. My question is what would be an appropriate analysis for these 2 situations? Please keep in mind that the same sites were sampled on two dates. Right now I am considering using mixed models (`lmer`, `lme` in R), but I am interested in your expert opinions. ![enter image description here](http://i.stack.imgur.com/vwaz0.png)                dat <- structure(list(Site = structure(c(1L, 1L, 2L, 2L, 3L, 3L, 4L,          4L, 5L, 5L, 6L, 6L), .Label = c("C1", "C2", "C3", "Q1", "Q2",          "Q3"), class = "factor"), date = structure(c(1L, 2L, 1L, 2L,          1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L), .Label = c("early", "late"), class = "factor"),              Y1 = c(52L, 33L, 50L, 39L, 57L, 45L, 45L, 30L, 37L, 33L,              45L, 41L), Y2 = c(10L, 4L, 7L, 4L, 7L, 4L, 9L, 3L, 8L, 4L,              8L, 4L), X = c(709.75, 709.75, 362, 362, 957.25, 957.25,              478, 478, 299.5, 299.5, 551.25, 551.25)), .Names = c("Site",          "date", "Y1", "Y2", "X"), row.names = c(NA, -12L), class = "data.frame")


title: Multicollinearity in a "population model"

text: I was talking with a colleague who told me that at the time of making logistic regression across a population did not have to worry about assumptions such as multicollinearity, because when analyzing the entire population estimated betas are "actual parameters" and since to "multicollinearity only affects the significance of the estimates" not matter in this case. This got me thinking, could it be true? If so, under what circumstances? I was reading some articles that talk about when multicollinearity is not important, but does not mention anything about population analysis. On the other hand I was reading some answers on population analysis and understand that suggest treating this type of regressions as if a large sample from a super-population analyzed. What would be the best position on this type of analysis (population) and multicollinearity (excluding the cases mentioned in the articles linked)?


title: Is there an overfitting in the data?

text: In my training set, I have 4026 instances, two classes and 104 attributes. The training error (error on the entire training set) is about 14% and 10 fold CV error on the training set is about 27%. Does this mean that my model is overfitting the training data? **EDIT (More Results):**   * Total Instance 4026 (+: 951, -:3075)   * Model Random forest (default settings in Weka)   * Training Error:14%, 10 Fold CV Error: 27%, Training F-score (+):67%, Training F-score (-): 91%, 10 Fold CV F-score (+): 31%, 10 Fold CV F-score (-): 83%   * After Random Undersampling, Total Instance 1902 (+ and -:951 each)   * After Random Undersampling, Training Error:16%, 10 Fold CV Error: 36%, Training F-score (+):85%, Training F-score (-): 84%, 10 Fold CV F-score (+): 65%, 10 Fold CV F-score (-): 63%


title: Difference Estimate & Confidence Interval for Paired Mann-Whitney Test

text: I have two paired samples and I am running a paired Mann-Whitney test to determine whether they are different. In Matlab, this takes the form of p = signrank(x,y), which tests the null hypothesis that x-y comes from a population of median zero. Unfortunately, Matlab doesnt give difference estimates for x-y or confidence intervals for this estimate. In the past, to estimate the shift estimate for nonpaired samples, I used the Hodges Lehman estimate (see Difference of 'centers' of 2 non-normal samples with Mann-Whitney test & http://en.wikipedia.org/wiki/Hodges%E2%80%93Lehmann_estimator). This is essentially the median of all pairwise differences. Here, since the samples are paired, I imagine that the HL estimate is incorrent. Intuitively, I was thinking of using the median of all paired differences. In either case, my question is how to get the difference estimate for x-y and the confidence interval for that estimate. Below is a picture of my data: The top histogram shows both samples (correct & incorrect), the bottom histogram shows the pairwise differences. A boxplot of the data is on the right. The estimated p-Value for the paired MW test is above & is <.05. Thus the data is significantly difference with correct > incorrect (even though it admittedly is hard to see on the histogram. Any suggestions on how to make it more visualisable are greatly appreciated!) ![data](http://i.imgur.com/dV04A2V.png)


title: Random intercepts model - one measurement per subject

text: I have repeated measurements on several hundred subjects. I plan to fit a random intercepts model. Most subjects have 3 to 4 measurements but about 20% only have one measurement. For the subjects with one measurement, should they share the same intercept or should they each have their own intercept as well? I'm not sure if there is enough information for the latter. I plan on doing the modelling in BUGS where I can control the intercepts for each subject.


title: Conducting IPD meta-analysis when studies use different measures

text: I'm interested in conducting an individual patient data (IPD) meta-analysis in which the main outcome of interest is a correlation coefficient between and IV and DV. The IV is measured slightly differently across studies, and conversion to make them all match might not be possible. Clearly, I can still do a two-stage IPD meta-analysis based on correlation coefficients calculated for each study individually. I would also like to conduct a one-stage IPD analysis for a few reasons: to more flexibly model an important covariate, to assess whether the relationship between IV and DV is linear, and to model aggregation bias (per Stewart et al. 2012). Given the heterogeneity in the IV measures, is this still possible?


title: Estimate absolute residuals in Stata

text: I have to regress the following equation in Stata: $$R_{i,t} = \delta{i} + \sum_i \alpha_i*W_i + \sum_i \beta_i * R_{i,t−n} + \mu_{i,t}$$ I'm ok with this regression. My problem comes out when I have to regress the second regression of the model I'm studying, that is: $$|μ_{i,t}| = \omega_i + ϕiM_{i,t} + \sum_i \rho_i*|μ_{i,t−j}| + \phi_iiA_{i,t} + \eta_{i,t} $$ In this regression I regress the residual I found in the first regression on its past value. How can I get the value of absolute residual? Is exist a particular formula by computing? Or do I have to consider simply the absolute value of the number $\mu_{i,t}$ I find in the first regression?


title: Design of matrix of contrasts in R

text: I am doing some post-hoc comparisons (in lme4, but here I'll just present a simple linear model), and I am having a hard time making sure that I am building the right matrix of contrasts to test differences between some combinations of factors. **Here is my model:**               dataset <- data.frame(response=rnorm(1200), factorX=rep(c("1", "2"), 600), factor2=rep(c("A", "B", "C"), 400))          model <- lm(response ~ factorX*factor2,data = dataset)      In my dataset, I measured a response variable depending on two factors: factorX (which has 2 levels: 1 and 2, 1 being the reference level) and factor2 (which has 3 levels: A, B and C, A being the reference level). **I get this model output:**               Coefficients:                       Estimate Std. Error t value Pr(>|t|)     (Intercept)        0.02918    0.07183   0.406    0.685     factorX2          -0.03780    0.10158  -0.372    0.710     factor2B           0.02027    0.10158   0.200    0.842     factor2C          -0.10972    0.10158  -1.080    0.280     factorX2:factor2B  0.01409    0.14366   0.098    0.922     factorX2:factor2C  0.08925    0.14366   0.621    0.535      **Question** : I want to know if there is a significant difference between some combinations of factor levels, so I build the following contrasts, and wanted to know if they corresponded to the right comparisons: `c(0, 0, 0, 0, 0, 1)` => factorX1_factor2A vs. factorX2_factor2C `c(0, 0, 1, 0, -1, 0)` => factorX1_factor2B vs. factorX2_factor2B `c(0, -1, -1, 0, 0, 1)` => factorX2_factor2B vs. factorX2_factor2C `c(0, 0, 0, 1, 0, 0)` => factorX1_factor2A vs. factorX1_factor2C `c(0, 1, 0, 0, 0, 0)` => factorX1_factor2A vs. factorX2_factor2A Is it correct?! Thanks!


title: Illustrative datasets and analysis for multilevel modelling

text: I recently took an introductory course on multilevel modelling. Most of the datasets and examples we used were from the social sciences. I've just got a 2 week internship in a biostatistics department, where they want me to start a project concerning the variation at hospital level of patient outcomes for an emergency condition that has a high mortality rate, both between hospitals and over a 5 year time span. I'm starting the internship next week, and I was hoping to find a book or online resource where a similar analysis (preferably with R, Stata or MLwiN) has been done, preferably one where they provide their data sets for the reader. Any links would be most welcome. Edit: I will be working with a dataset detailing all recorded aspects of the patient's in-hospital care. The main outcome of interest is death within 30 days of admission.


title: Assessing quality of similarity measure

text: I'll try to compute different similarity measures on some data. The result is always a similarity matrix on the data. I want to evaluate the quality of the matrix by checking if the values for positive and negative similarities have a wide spread so that one can set an threshold which is very general to all data. Means, postive and negative similarities should not be close to each other. I have no way to check precision and recall because I'm not in the postition to select threshold yet. I first want to check the optimal similarity measure. Is there a standard procedure for such a task?


title: Estimating adjusted rate ratios for categorical independent variable measured only as population percents

text: I have counts of people with a condition, by age (5-year groups), sex, and area (100-200 areas). For each area, I have:   * estimated population by age (5-year groups) and sex,   * estimated population percent living in city, small town or rural (sum to 100%),   * estimated population in high, middle and low socioeconomic status (sum to 100%), and   * estimated population proportion treated (via food or water supply fortification). I know nothing about individual-level correlation between rurality and socioeconomic status (or even if it's the same across areas). I expect proportion treated is near 100% in cities, lower in small towns, much lower in rural areas. Researcher wants adjusted rate (or odds) ratios for having the condition (and p-values) for rural vs cities, low vs high SES, untreated vs treated, young vs older. Can I produce valid estimates? If so, how should I model it?


title: If I repeat every sample observation in a linear regression model and rerun the regression how would the result be affected?

text: Say I have N observations, possibly multiple factors and I repeat each observation twice (or M times) how would a regression on this new set of size NM compare to a regression on just the original observations?


title: PCA followed by LDA: dimension reduction strategy

text: I have a high dimensional dataset (n*p: 30 * 100) which I want to use as an testing dataset to build a two group classifier (LDA or QDA). I've read that you can do PCA to do an dimension reduction of your dataset to select to most important features. But I'm a bit confused what you use exactly as the input to build the classifier. I'm familiar with PCA using SVD and what it means. Consider following situation: - I do a SVD of my dataset. - I look at the scores of the first couple of principal components. - When I assign my scores a label indicating from which group they come, I see that the 3th PC best separates my 2 groups. (allthough it only explains 7% of the total variance) What do I do next? 1. I take the 3th PC transform to the original parameter space (scores * loadings * scale + mean) and build my classifier 2. I look at the loadings in the 3th PC and try to decide which parameters in my original parameter space are important and build a classifier only using these. 3. ... Option 2 seems the most sensible in my opinion but I'm not entirely sure. Also If I see that only the 3th PC is important to explain the variance in my two groups, can I forget about the first two PC in my further analysis? Many thanks in advance.


title: Solving the Kolmogorov forward equation for transition probabilities

text: Let $\lambda \mu > 0$ and let $X$ be a Markov chain on $\{1,2\}$ with generators $$ Q = \begin{pmatrix} -\mu & \mu \\ \lambda & -\lambda \end{pmatrix}$$ Write down the forward equations and solve them fro the transition probabilities, $p_{i,j}(t), i,j = 1, 2$. The answers say: "We have that $$p_{11}' = -\mu p_{11} + \lambda p_{12}, \, \, \, p_{22}' = -\lambda p_{22} + \mu p_{21}$$ where $p_{12} = 1 - p_{11}$ and $p_{21} = 1 - p_{22}$. Solve these subject to $p_{ij}(t) = \delta_{ij}$, the Kronecker delta, to obtain that the matrix $P_t = \{p_{ij}(t)\}$." The bits I don't understand are these 1) How do they get those equations in the first place? For example, why have they added $\lambda p_{12}$ when the coefficient for $p_{12} = \mu$? 2) How do you solve these equations? What do you do with that Kronecker delta to solve them?


title: Will Multivariate Gaussian classifier work for text classification?

text: So far i have evaluated mn Bayes and Bernoulli, so my question is if i take the counts of the words of each document and use them for assigning the document to the particular class will it work with Multivariate Gaussian classifier (Bayes with Gaussian model)?


title: SVMs for nested/prioritzed text document classes.

text: Let's say A, B, C, D, E, F are classes of text documents (i.e. sentences). Usually they occur distinctively, say: one document is precisely C. But sometimes one document is C AND A, at the same time. Because A is more important than C, the document class is defined . So at least two questions arise: _1) How this classification would be made properly with an SVM?_ _2) Is this possible in R?_ Note: I already tried to do this with the standard approach (linear classification, 10000 features, e1071 SVM implementation). The classes distribution isn't constant (F is much more frequent than A, B, or C). Sadly here the most most important classes (A, B, C) are wrongly classified as F. _3) Is there a possibility to set a priority for an standard SVM?_


title: factorGraph approximations: Splitting variables

text: I am reading Chris Bishop's chapter on Expectation Propagation and there is the bit on how approximating factor graphs as shown in the images. ![enter image description here](http://i.stack.imgur.com/wkAVG.png) So, the original joint distribution can be written as: $$ p(x) = f_a(x_1, x_2) \; f_b(x_2, x_3) \; f_c(x_2, x_4) $$ Now, the assumption is that each of these factors also factorise according to the variables they depend on, so we have an approximating distribution as: $$ q(x) = \widetilde{f}_{a1}(x_1) \; \widetilde{f}_{a2}(x_2) \; \widetilde{f}_{b2}(x_2) \; \widetilde{f}_{c2}(x_2) \; \widetilde{f}_{b3}(x_3) \; \widetilde{f}_{c4}(x_4) $$ Now, just let us consider the original factor $f_a$. It has been split into $\widetilde{f}_{a1}(x_1)$ and $\widetilde{f}_{a2}(x_2)$. So, I am guessing the original factor which depended on both $x_1$ and $x_2$ can be approximated as a product of two unary variables i.e. $f_a(x_1, x_2) \approx \widetilde{f}_{a1}(x_1) \widetilde{f}_{a2}(x_2)$. Now imagine that the original factor $f_a$ represented some joint probability as $P(x_1|x_2) P(x_2)$. Now imagine this $P(x_1|x_2)$ is a multivariate normal and $P(x_2)$ is a gamma distributed. They both represent priors in the joint distribution. Now, how does this relate to the approximating factors $\widetilde{f}_{a1}(x_1)$ and $\widetilde{f}_{a2}$? Can I give them whatever distributions I find suitable? For example, can I make $\widetilde{f}_{a1}(x_1)$ a multivariate normal and $\widetilde{f}_{a2}$ as a Gamma distribution? What about any other arbitrary distribution? Is it unto the user discretion to choose the appropriate factor forms? The reason I ask this is that ultimately the full approximate joint distribution has a form typically a multidimensional normal distribution and there must be some constraints on what can be chosen for these factor messages? This is on page 514-515 of the Bishop book. he does not speak much about the form of these factor messages? I did ask a similar question before but it did not get much attention. I tried to make some progress on it since then but so far have been unsuccessful.


title: How to create a correction factor based on measures from a known material?

text: I work with medical imaging and have the following issue: -- every time I scan a patient I obtain images of a calibration phantom alone, followed by images of subject + calibration phantom in the same image Although the calibration phantom has a stable composition, measures of its density vary from scan to scan because of   * equipment variation   * subject size   * other unknown factors Because of these fluctuations, ideally I should apply a correction factor to the measures I perform on a subject. How can this be achieved with the data I am collecting?


title: Calculate quality estimates for distribution of real data

text: Imagine I have measured heights of humans (in cm) and show the data in a histogram:               # Generate sample data     sampleSize = 100000     sampleData = round(rnorm(n = sampleSize, mean = 175, sd = 14))          # Draw histogram of sample     h = hist(sampleData, breaks = max(sampleData) - min(sampleData), plot=T)      Now I want to know how realistic each size in the distribution is. For example I found a human with the size of 100cm and I want to know if this individual is overrepresented in my sample. **I want to give a measure for the uncertainty for each height class.** How can I do that? My first idea was to bootstrap the measured data, calculate the standard deviation for each class among the bootstrapped samples and give a relative standard-deviation (scaled by 1000) for each class.               ######################################################################     # Calculate the standard Deviation for each sampled size of a bootstrap measurement     cycles = 100000          # initialize results-matrix     result = matrix(0, nrow = length(h$counts), ncol = cycles)          # sample original data     for (i in 1:cycles) {         tmp = hist(sample(sampleData, replace = T, size = sampleSize), breaks = h$breaks, plot=F)         result[,i] = tmp$counts         print(paste("cycle:", i))     }          # calculate standard deviations     for (i in 1:length(h$counts)) {         h$sdboot[i] = sd(result[i,])     }          ######################################################################     pdf("boot.pdf")         # Plotting         #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++         plot(h)         # Bootstrap: standarddeviations         polygon(x=c(rev(h$mids), h$mids), y= c(rev(h$counts + h$sdboot), h$counts - h$sdboot), col="#ff000022")              #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++         plot(h)         # Bootstrap: standarddeviations         polygon(             x=c(rev(h$mids), h$mids),             y= c(rev(h$counts + 1000 * (h$sdboot / h$counts)), h$counts - 1000 * (h$sdboot / h$counts)),             col="#00ff0022"         )              #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++         plot(x=h$mids, y=h$sdboot/h$counts)          dev.off()      The resulting graphs (2 and 3) show the increasing uncertainty with lower numbers of counts in the class, but I'm not sure how informative that is. It does not really answer the question " **How likely is it, that you miss members of a certain height class in a measurement with the same sampleSize?** " or " **How likely is it, that the value is of the measured size just by chance** " One problem might be, that I do not have a model for my real data (For human height one could assume a normal distribution). Is it somehow possible to give a measure for the uncertainty of my real data?


title: Cluster analysis or regression?

text: I have a table with each row representing single printer model, its features, and price. I want to know how price is formed based on these features. What should I start with? Multiple regression, so I could cut off insignificant features? Cluster analysis to define small clusters with equal price? What are the ways to do the task?


title: Classification with increasing number of classes

text: I have dataset with large number of labelled data (say that there are _k_ classes). I have also another, much smaller dataset with unlabelled data that I want also to be labelled. The problem is that in the second dataset, number of classes need not to be the same as number of classes in first dataset (more precisely, the there may be more classes). In other words, some objects from unlabelled dataset can be classified to one of _k_ possible classes since they are "too close of them" but some objects should be clasified to new classes. Real world example: In some applications, number of classes can increase with time. For example, number of known species increases with time, so when I see any unknown plant in nature, biologist can either to classify it or to say that this is new discovery (new class). Is there any technique that deals with this type of classification? Any help is appreciated.


title: How to define in R the most important variables?

text: I have a data set with my target gene and more than thousand transcriptional factors somehow correlated with this gene. There is data of these factors in more than 70 variable conditions. What I'm looking to do is use CV to find most useful variables which correlate with my target gene. How could I do that in R please?


title: Interaction depth parameter in GBM

text: In the GBM package one is supposed to be able to provide interaction.depth>2, which means higher-order interactions between features. However, the resulting trees (as seen by pretty.gbm.tree) never show such interactions (and indeed - each row corresponds to just a single feature). I'm not even sure anymore that the package actually supports depth>2.... Does anybody have any idea?


title: Longitudinal panel dataset: Consequences of missing values

text: I am analyzing a longitudinal panel dataset using OLS. The data spans around 40 years, but for some variables data was unavailable for certain categories. In most cases, the data for given category stretches not as far back as for others (or ends earlier), but there are also a few cases where a category has an entire blank series. Unfortunately there is no way that I can retrieve the missing data. Hence, I will have to estimate my models with gaps in my dataset. I am wondering whether there are consequences for my estimation results that I should be aware of. I'm thinking here statistical biases or issues regarding the validity of my results. Suggestions would be very much appreciated.


title: Guessing the appropriate 'distribution' of an observation to perform some form of sensitivity analysis

text: I wish to run some simulations on the 'estimates' of an event provided by experts (or historical data if available). For example: The technology manager of company X says that probably 300-1000 computers in the organization (out of say 100,000 in total) have Windows XP that needs to be phased out. Now this is an estimate, he's not sure of the number. Let's take 2 scenarios:   1. It's worth spending time/money to sample and estimate (i.e., value of information is worth the effort)   2. It's not worth the effort to perform any form of sampling and let's go with the estimate. In either case, I wish to 'simulate' this model (with other similar) variables. You may think of it is wanting to run Monte-Carlo Simulations, but I don't want to put a specific solution to the idea, for the sake of discussion. The question is then this - given sampling is done, can you 'guess' the appropriate probability distribution (I think this is doable, if I remember correctly from statistics, not sure though). And if sampling is NOT performed, how best to guess it? Note: If my question is 'fuzzy' I'll surely revise it as I receive some comments/answers to make it more clear. The example above may not be the best one but nevertheless let's go with it for the sake of example. If a better example is suggested, I'll surely add it to the question


title: Unbalanced Haar Wavelets and R package unbalhaar

text: I have found the article by Fryzlewicz (2007) and his R package unbalhaar. I do not understand well the outputs of his functions and in particular of the `best.unbal.haar.bu()` function. Here is the documentation of it : best.unbal.haar.bu. From his article, I have understood that the `$smooth` argument is the constant factor, and that `$detail[1,]` is the location of the coefficients and `$detail[3,]` is the value of the coefficients. But I do not see what could be `$detail[2,]` called "weights" in the documentation... Does anyone know about this work and know what it is ? Thank you for your help. PS : I have not posted this question on stackoverflow because I think it is more about the maths than about the coding thing.


title: Finding moments for a theoretical density function

text: I am working on finding higher order moments for a given theoretical function, to be used in modelling of daily log-returns. The PDF is, $f_r(x) =$ $\begin{cases} \quad \frac{1}{2}ae^{a\left(x-\mu\right)} & \text{if x } < \mu \\ \quad \frac{1}{2}be^{-b\left(x-\mu\right)} & \text{if x } \geq \mu \\ \end{cases} $ I believe the mean and variance are correct as, E[X] =$-\frac{1}{2 a}+\frac{1}{2 b}+\mu$ and Var[X] = $\frac{3}{4 a^2}+\frac{1}{2 a b}+\frac{3}{4 b^2}$ In computing skewness I have for $\mu_3 = \acute{\mu_3} - 3\acute{\mu_1}\acute{\mu_2} +2\acute{\mu_1}^3 $, $= -\frac{7}{4 a^3}-\frac{3}{4 a^2 b}+\frac{3}{4 a b^2}+\frac{7}{4 b^3}$. And for the kurtosis, $\mu_4 = \acute{\mu_4} - 4\acute{\mu_1}\acute{\mu_3} +6\acute{\mu_1}^2\acute{\mu_2} -3\acute{\mu_1}^4 $, $= \frac{117}{16 a^4}+\frac{15}{4 a^3 b}+\frac{15}{8 a^2 b^2}+\frac{15}{4 a b^3}+\frac{117}{16 b^4}$ My first questions, are the $\mu_3$ and $\mu_4$ the correct moments for Skew and Kurt? When I am trying to estimate the parameters to fit the theoretical and sample moments, it does not seem so.. I am happy to provide more details etc if deemed necessary to clarify the issue. Thanks.


title: The prisoner paradox

text: I am given an exercise, and I can't quite figure it out. > ### The Prisoner Paradox > > Three prisoners in solitary confinement, A, B and C, have been sentenced to > death on the same day but, because there is a national holiday, the governor > decides that one will be granted a pardon. The prisoners are informed of > this but told that they will not know which one of them is to be spared > until the day scheduled for the executions. > > Prisoner A says to the jailer “I already know that at least one the other > two prisoners will be executed, so if you tell me the name of one who will > be executed, you won’t have given me any information about my own > execution”. > > The jailer accepts this and tells him that C will definitely die. > > A then reasons “Before I knew C was to be executed I had a 1 in 3 chance of > receiving a pardon. Now I know that either B or myself will be pardoned the > odds have improved to 1 in 2.”. > > But the jailer points out “You could have reached a similar conclusion if I > had said B will die, and I was bound to answer either B or C, so why did you > need to ask?”. > > What are A’s chances of receiving a pardon and why? Construct an explanation > that would convince others that you are right. > > You could tackle this by Bayes theorem, by drawing a belief network, or by > common sense. Whichever approach you choose should deepen your understanding > of the deceptively simple concept of conditional probability. Here's my analysis: This looks like the the Monty Hall problem, but not quite. If A says `I change my place with B` after he is told C will die, he has 2/3 chances to be saved. If he doesn't, then I would say his chances are 1/3 to live, like when you don't change your choice in the Monty Hall problem. But at the same time, he is in a group of 2 guys, and one should die, so it is tempting to say that his chances are 1/2. So the paradox is still here, how would you approach this. Also, I have no idea how i could make a belief network about this, so i'm interested to see that.


title: PCA on returns or levels?

text: I am doing a PCA on a few economic indices (which I standardized) (24 timeseries). When doing the PCA on the values/level of the index, I get rather similar loadings on each index, which I think is good, as no one should dominate. However when I run it on the returns of these standardized indices, the loadings are much more extreme on 1 or 2 and almost none for the rest. Also when doing PCA on both levels and returns, I only get non negligible loadings on the levels! What is the interpretation of this, how are these different results possible? Additional info: I am doing it in matlab with the princomp function


title: What is the difference between generalized estimating equations and GLMM?

text: I'm running a GEE on 3-level unbalanced data, using a logit link. How does this differ (in terms of the conclusions I can draw and the meaning of the coefficients) from a GLM with mixed effects (GLMM) and logit link? More detail: The observations are single bernoulli trials. They are grouped clustered into classrooms and schools. Using R. Casewise omission of NAs. 6 predictors also interaction terms. (I am not flipping children to see if they land heads-up.) I'm inclined to exponentiate the coefficients to odds-ratios. Does this have the same meaning in both? There is something lurking in the back of my mind about "marginal means" in GEE models. I need that bit explained to me. Thanks.


title: Narrow confidence interval -- higher accuracy?

text: I have two questions about confidence intervals: Apparently a narrow confidence interval implies that there is a smaller chance of obtaining an observation within that interval, therefore, our accuracy is higher. Also a 95% confidence interval is narrower than a 99% confidence interval which is wider. The 99% confidence interval is more accurate than the 95%. Can someone give a simple explanation that could help me understand this difference between accuracy and narrowness?


title: Goodness of fit for orthogonal least square regression

text: We have a set of XY coordinates which we have fitted with an orthogonal least square regression model. We have a vector of residuals for each point and the fitted line. How do we assess goodness of fit quantitatively? We're currently using Matlab but are open to other approaches.


title: Alternatives to Pearson correlation

text: We are trying to predict values of variable A having N other variables. What we do, we calculate Pearson correlation between A and each of the other N variables, for last M values, using fixed M. We use variable with largest correlation coefficient as predictor. This scheme works fine when analyzing N variables during 24 month and decreasing between month 25-36.   What could come as improvement of Pearson correlation in this context? Variables are bi-valued.


title: Scrambling and correlation in low discrepancy sequences (Halton/Sobol)

text: I am currently working on a project where I generate random values using low discrepancy / quasi-random point sets, such as Halton and Sobol point sets. These are essentially $d$-dimensional vectors that mimic a $d$-dimensional uniform(0,1) variables, but have a better spread. In theory, they are supposed to help reduce the variance of my estimates in another part of the project. Unfortunately, I've been running into issues working with them and much of the literature on them is dense. I was therefore hoping to get some insight from someone who has experience with them, or at least figure out a way to empirically assess what is going on: **If you have worked with them:**   * What exactly is scrambling? And what effect does it have on the stream of points that are generated? In particular, is there an effect when the dimension of the points that are generated increases?   * Why is it that if I generate two streams of Sobol points with MatousekAffineOwen scrambling, I get two different streams of points. Why is this not the case when I use reverse-radix scrambling with Halton points? Are there other scrambling methods that exist for these point sets - and if so, is there a MATLAB implementation of them?  **If you have not worked with them:**   * Say I have $n$ sequences $S_1, S_2, \ldots,S_n$ of supposedly random numbers, what type of statistics should I use to show that they are not correlated to each other? And what number $n$ would I need to prove that my result is statistically significant? Also, how could I do the same thing if I had $n$ sequences $S_1, S_2, \ldots,S_n$ of $d$-dimensional random $[0,1]$ vectors? * * * **Follow-Up Questions on Cardinal's Answer**   1. Theoretically speaking, can we pair any scrambling method with any low discrepany sequence? MATLAB only allows me to apply reverse-radix scrambling on Halton sequences, and am wondering whether that's simply an implementation issue or a compability issue.    2. I am looking for a way that will allow me to generate two (t,m,s) nets that are uncorrelated with each another. Will MatouseAffineOwen allow me to do this? How about if I used a deterministic scrambling algorithm and simply decided to choose every 'kth' value where k was a prime?


title: Building background for machine learning for CS student

text: I am a CS graduate student and I am starting to get really interested in Machine Learning (and Predictive Analytics). I have started working on a text classification project with a professor to learn the field but I am realizing that my background is pretty weak. The professor is too busy to teach me the basics, so I have to do it on my own. I have done some calculus, stat 101 and linear algebra in the distant past, but I do not remember too much. But I am pretty sure if I get a good book, I can pick things up pretty quickly. As of now I have started off my Machine Learning studies by reading Alpaydin's Machine Learning but it is pretty challenging given the big holes in my background. I am looking for suggestions on books on major topics (Linear Algebra, Statistics, Probability, Optimization, and anything else that might be relevant) that would help me ramp up relatively quickly, given my background. I am looking at doing it in multiple iterations. First, get what I need to do be able to do get some practical work done (like my current classification project), and then go back and read books to get deeper understanding, and so on. Please note that at this point I am only looking at formal treatment of the subjects ( **textbooks only** ). **EDIT 1** : I think I may have found stats and probability book (All of Statistics: A Concise Course in Statistical Inference by Larry Wasserman). Just need to figure out how to put the remaining pieces together. I welcome any further ideas in this regard.


title: Computational geometry vs perceptron for shattering

text: Given two point sets, $B$ consisting of blue points, and $R$ of red points, on the plane. The problem is to formulate a theoretical model to compare the average runtime of the Computational Geometric (CG) Algorithm and the  Perceptron Learning (PL)  for shattering two sets of points by a line. ( Briefly, CG tests if there exists a line, $ax+by = c$ such that the system of equations $ax_r+by_r\leq c; ax_b+by_b \geq c$ for all $(x_b,y_b) \in B$ and $(x_r,y_r) \in R$ is consistent; if so, it finds one such.) AFAIK, only empirical models are known thus far, where we take a large number of such point sets and test out the performance of CG and PL. While I can't see any clear solution to this, I am looking for one along the following lines: Choose 2 random integers, $b$ and $r$ in the range $[1,n]$ where $n$ is a large integer, to appear in the answer as a parameter. Choose $b$ blue points and $r$ red points within the unit square. Compute the average runtime performance of CG (essentially determining the complexity of determining consistency of linear equations on a random set of values) and PL (essentially determining the number of weight changes on a random data set) by a multi-fold integration. Further, average across all possible choices of $b$ and $r$ in the range $[1,n]$. Can you offer more insight or a better formulation?


title: How to obtain confidence limits of predicted values in ARIMA?

text: How can one obtain confidence limits of predicted values in ARIMA?


title: Using kurtosis to assess significance of components from independent component analysis

text: In PCA eigenvalues determine the order of components. In ICA I am using kurtosis to obtain the ordering. What are some accepted methods to assess the number, (given I have the order) of components that are singificant apart from prior knowledge about the signal?


title: Divide and Conquer

text: Satellite-Routed sensor system has attracted attention in next generation network. In this system large number of sensor terminals send sensed data to the monitoring station located in different area via a satellite. The problem of this system that large number of sensor terminals trying to transmit their sensed data simultaneously to satellite (collision will be occur). Suppose that each channel in the satellite are serving 100 sensors, and one sensor can communicate with this channel at a time. At a given time suppose 4 sensors have data to send for the satellite channel so collision will be occur, in this manner the satellite divide the 100 sensors into two groups and see in the first group how many sensors have data to send if one sensor this is ok if more than one sensor again the satellite divide the 50 sensors into two groups and so on. So the satellite decide division when more than one sensor has data to send at one time ( this divide and conquer need a lot of time in division ). So if you have an efficient way how to divide these sensors into groups I will be thankful for you ? so I need a modified way of divide and conquer to fastest the division process ? I need a better way in this way than divide and conquer All the best


title: How to properly analyze distance from a reference?

text: I'm measuring distances of various samples from a reference point. The distance is defined as a non-negative number, where $d=0$ means that the test case is identical to the reference. My general question is: Given a set of "typical" distances, what is the proper way to tell whether a given $d_1$ "too large", compared to the "typical"? In my particular case the distance distribution is shown on the following graph ![enter image description here](http://i.stack.imgur.com/BcFYJ.png) I failed to transform these data to anything symmetrical, so I can't use normal approximation. Any suggestions?


title: How to design an 8-run experiment in 5 factors?

text: I'm trying to design an experiment where I measure a variable as a function of 5 two-level factors, labelled A, B, C, D and E. I'm trying to understand how to best design this experiment so I can conduct it in 8 runs. I've tried to follow the guidance given in Box, Hunter & Hunter, and found two experimental $2^{5-2}$ designs that seem to contradict each others. One is given in the table p. 272, where D=AB and E=AC. That would yield the following design:               ABCDE     +++++     ++-+-     +-+-+     +----     -++--     -+--+     --++-     ---++      The other one is the example whose table is given on p. 236, defined by D=BC and E=ABC:               ABCDE     +++++     ++---     +-+--     +--++     -+++-     -+--+     --+-+     ---++      On what criteria should one choose one design over another? Sorry if this sounds like a newbie's question, but I'd really like to understand this issue.


title: Simulate ARIMA by hand

text: I was working on ARIMA in R and I am trying not to use library `forecast` as much as possible. I have a code for finding the best ARIMA model, but it is showing some warning messages. Here is my function:               best.aic<-1e8     n<-length(x.ts)     for(p in 0:maxord[1]) for(d in 0:maxord[2]) for(q in 0:maxord[3])     {                fit<-arima(x.ts,order=c(p,d,q))         fit.aic<--2*fit$loglik+(log(n)+1)*length(fit$coef)         if(fit.aic<best.aic)         {             best.aic<-fit.aic             best.fit<-fit             best.model<-c(p,d,q)         }     }     list(best.aic,best.fit,best.model)      It gives the same warning 4 times:               In arima(x.ts, order = c(p, d, q)) :     possible convergence problem: optim gave code=1      Also, I will appreciate any help concerning simulation of ARIMA by hand, not by `arima.sim` function. Thank you.


title: Number of levels in categorical predictors/combining categories in logistic regression

text: I have an independent variable “State in the US” and it has around 16 different levels (PA, NY, NJ, etc). Now I have combined the states which appear infrequently into "others" category. But I still have around 12 levels which are equally populated. I have a very large data set (1 million observations) and no zero cell issues with the variables. Questions:   * I would like to know if there was a limit on the number of levels a categorical variable can have when used with logistic regression. If yes, why is it the case?   * And what is the best way to bin categorical variables like Zip code or occupation code, apart from binning similar levels together?


title: Trouble writing likelihood

text: I'm having trouble writing the likelihood for a homework question: Suppose $X=Z(Y+\theta)+(1-Z)(\theta-Y)$ and $X_1, X_2, X_3 \overset{iid}{\sim} f_X(x)$ where $Z\sim \rm{Bernoulli}(0.5)$ and $Y\sim \rm{exponential}(1)$ and $Z$ and $Y$ are independent. Let $x=(7.2, 9.3, 8.6)$ be the observed data. Can you guys help me (give hints, ideas) on how to write the likelihood (or log-likelhood function)? So far I have: $L(\theta;x) = f_X(x|\theta) = \prod_{i=1}^3 f_X(x|\theta) = \prod_{i=1}^3 f_Z(??|\theta)f_Y(???|\theta)$ but I can't figure out what to do really in this strange situation.


title: Logistic Regression/Naive Bayes with highly correlated data

text: **Background:** We work with data from sports event, more accurately with data about the spectators of sports events: how many people are being violent, what kind of event is this, etc. We have quite a lot of data from the past few years and we try to find the "right" number of security people we need to minimize the violence while keeping some kind of "budget". **Aim:** we want to forecast for a given set of explanatory variables (weather, type of sports event, location, etc.) the expected violence (low, medium, high) depending on the number of "security guards" that are present at the game. **Problem:** The historical data is of course highly correlated: the number of security guards is in some measure proportional to the violence, and also very related to the other variables (as they probably used some safety experts to evaluate the danger). Using the non-collinearity assumption to make a naive bayes seems wrong. > **Question:** What is a correct way to forecast the violence depending on > the number of guards present during the event? **My guess:** I should "discretize" the number of guards in 3-4 "bins" (e.g. few, some, many, a lot) to remove some of the correlation, and use the corresponding training sets to forecast the output violence depending on the input variables. But I lose a lot of information by training only on a subset of my data.


title: Inference for the skeptical (but not math-averse) reader

text: I just watched a lecture on statistical inference ("comparing proportions and means"), part of an intro to stats online course. The material made as little sense to me as it always does (by now I must have seen this stuff dozens of times, spread out over the last three decades). > I'm looking for a book on "basic Stats-101" (point estimation, estimate > assessment, statistical inference, hypothesis testing, study design) that > takes seriously the problem of convincing a skeptical reader... Below I give some examples of the _type_ of question that the author I'm searching for would take seriously and know how to address convincingly. But first let me take a minute to stress that in this post **_I'm not asking these questions._** _Please, do not answer them!_ I give them just as examples, and by way of "litmus test" (for the type of author of searching for).   1. If a "proportion" is simply the _mean_ of a Boolean variable (i.e. one that takes only the values 0 and 1), why are different procedures taught for doing statistical inference with "proportions" and with "means"?   2. If the normal distribution is so robust that assuming normality gives good results even in cases when that data is not quite normally distributed, and if the t-distribution is so normal-looking, why all the fuss about using the t-distribution instead of the normal?   3. What _exactly_ are "degrees of freedom", and why do we worry about them?   4. What does it mean to speak of the "true" value of a parameter, considering that we are just using distributions that happen to _look_ similar to the data?   5. How come "exploratory data analysis" is a good thing, while "data snooping" is an evil thing? As I've said, I'm put off by the _attitude_ that's implied by a neglect of such questions. It's not the "epistemological stance" that I want to see in someone who's teaching me something. I'm looking for authors who respect the reader's skepticism and rationality, and who know how to address them (without necessarily going off into pages and pages of formalisms and technicalities). I realize that this is a tall order, and maybe especially so when it comes to statistics. Therefore, I don't expect that many authors will have succeeded at it. But at the moment I'd be content with finding just _one_. Let me add that I'm not math-averse. On the contrary, I love math. (I'm comfortable with analysis [aka "advanced calculus"], linear algebra, probability theory, even basic measure theory.) That said, my interest at the moment is in "applied", "practical", "everyday", "real-world" statistics (as opposed to theoretical niceties). (But I don't want a cookbook either!) FWIW, I have read the first few chapters of _Data analysis using regression and multilevel/hierarchical models_ by Gelman and Hill, and I like the authors' tone. Their focus is practical, but go into the theory when needed. They also often step back, and assess standard practices critically, and offer candid opinions that appeal to a skeptical reader's commonsense. Unfortunately, these authors have not written a book devoted to the subject I'm asking about in this post ("Stats 101" stuff, as described above). I'm also aware that one of these authors (Gelman) co-authored the highly-regarded _Bayesian data analysis_ , but, again, this is not what I'm looking for at the moment. EDIT: Dikran Marsupial raises the following objection: > I don't think there is necessarily anything wrong with neglecting questions, > there comes a point where addressing every question detracts from the > exposition of the basic concepts which is often more important (especially > in a stats 101 book!). I agree with that. It would be more accurate for me to say that I'm looking for a "second look at basic stats." In fact, with this as my motivation, I looked at the textbooks used at graduate courses on inference (say), and found that they _too_ neglected questions like the ones I've listed. If anything, they seemed even less inclined to delve into such questions (so that they can focus on matters like the conditions for some convergence-or-other of this-or- that...). The problem is that the more advanced books are addressed to a _radically different_ population of readers, one where the "skepticism of the outsider" has been drastically depleted. IOW, those who are taking graduate-level statistics are _past the point of being bothered by the questions that bother me_. They're not skeptical about any of this stuff anymore. (How did they get over the skepticism hump? Maybe some were never too critical in the first place, especially if they learned their stats fairly early on--I know I was not a particularly critical freshman myself, for example, though I did not take stats then. Others may have had teachers who filled in where their textbooks fell short. Some may have been clever enough to figure out the answers to such questions for themselves. Who knows.)


title: What is the difference between the Monte Carlo (MC) and Monte Carlo Markov Chain (MCMC) method?

text: The goal of both methods seems to be to derive an estimate of a posterior/target distribution. If a process model exists which links some input parameters (which are themselves uncertain and can be described by a PDF) to an output parameter through a model equation or other computations, why would one choose one method over the other? Would both be applicable? Can one make a statement on the benefit of one method over the other with respect to the number of required draws/simulation runs in order to reach a sufficiently good approximation of the target PDF?


title: Optimizing "twenty questions"-like decision tree queries

text: I'm looking for an algorithm which can optimize the selection of queries used to build a decision tree. However, unlike most decision trees, I am constrained to ask the same set of questions for every case. I'll couch this in the form of the parlor game Twenty Questions . One person thinks of something (e.g. an elephant) and the other needs to guess the object based on a set of binary questions. For examples, "is it bigger than a breadbox?" and "could I put it in my mouth?". In principle, a perfect decision tree could select from up to 2**20 ~ 1 million different answers. This assumes that there is no correlation between one question and the next. (There are very few things bigger than a breadbox that I can put into my mouth. Cotton candy is likely one of them.) As a point of strategy, and quoting from the above Wikipedia link: > Careful selection of questions can greatly improve the odds of the > questioner winning the game. For example, a question such as "Does it > involve technology for communications, entertainment or work?" can allow the > questioner to cover a broad range of areas using a single question that can > be answered with a simple "yes" or "no". If the answerer responds with > "yes," the questioner can use the next question to narrow down the answer; > if the answerer responds with "no," the questioner has successfully > eliminated a number of possibilities for the answer. The goal of this sort of broad question is to increase the selection/rejection point to 50% of the universe under consideration. This gives the questioner the most information. **Question 1: Given a set of several million single- topic questions with low discrimination, how does one construct a more informative multi-topic question which is the union of several of those low- level questions?** Also, as I mentioned, I have an extra constraint. I always have to ask the same set of questions. I can't use feedback from a previous question to drive the next. In my case, I'm building an inverted index as part of a chemical substructure search. There are some 3 million unique features in my target set of about 50 million records. These features include "contains a ring of size 5" and "contains at least 4 carbon atoms." About 1/2 of these features only exist in a single record. When a substructure query comes in, I can identify its features, get the corresponding inverted indices, and get the intersection. This can help speed up the search a lot. For legacy reasons, I am limited to storing 1024 indices. This should be enough to distinguish between 2**1024 different molecules, but unfortunately most features aren't that discriminatory, either because they are rare or because they are too well correlated with other features. **Question 2: How can I take this into account when selecting my list of pre-determined queries?** Algorithms like C4.5 for making decision trees are well known. However, C4.5 picks a single feature for each decision. I would like to enrich discrimination by merging multiple questions into one. ("Does it have at least 16 carbon atoms or at least 5 oxygen atoms"?) It must also have a cost model because it takes time to evaluate some of the more complicated features, and that time might not be worthwhile. I don't know where to go from here. I can't come up with a good heuristic myself. My hope is Twenty Questions is well-known, so there's existing work related to this topic that people can refer me to.


title: Introducing observation errors in jags code

text: I want to introduce observations errors around my data in Jags, but I face some trouble coding it without having a double definition error on node Y3 So far I have :               log_Y1_real[y]~dnorm(log(Y1[y]),prec_obs_Y1)     #Y1 is one variable     log_Y2_real[y]~dnorm(log(Y2[y]),prec_obs_Y2)     #Y2 is another variable     Y1_real[y]<-exp(log_Y1_real[y])     Y2_real[y]<-exp(log_Y2_real[y])          Y3_real<-log(Y1_real/Y2_real)     # My model is fit on Y3,defined as log(Y1/Y2), which should integrate obs error from Y1 and Y2               Y3_real~dnorm(mu_Y3,prec_Y3)           #then the model..     muY3~model....      Thanks a lot for your insights !


title: CoreNLP SemanticGraph - search for edges with specific lemmas

text: I'm using Stanford CoreNLP's dependency parser, and wondering how to make a generic search for SemanticEdge(s) with specific head lemma, dependent lemma, and lexical relationship. For example, if I have an actual dependency like this:   * dobj(discover-4, insights-6) How do I search for it using lemmas instead of the literal word and the index? Basically I want to be able to pattern match the parts of the dependency graph using generic rules...


title: Downsampling stationary time series data, effect on variance

text: Suppose I have stationary time series data, like this: ![enter image description here](http://i.stack.imgur.com/N7Oua.png) The time series has Gaussian noise around a true mean. If I then take windows of N samples across the series and average them to generate a new serie, How are the variance of the old and new series related?


I read the article "Two-way imputation: A Bayesian method for estimating missing scores in tests and questionnaires, and an accurate approximation" by Van Ginkel et al. (available here) about a Bayesian two-way imputation method for estimating missing scores in tests and questionnaires. Does anyone have suggestions about how to implement this? Is there a package or function (in e.g. R or SPSS) available that can deal with this method?


title: Finding the right parameter for a student's $t$-distribution

text: Suppose I have a homogeneous time series $x_n$, where the values arrive each second. For $n\ge 60$, let $\sigma_n$ be the standard deviation of the values $x_{n-59},\ldots,x_n$. Under the assumption that $x_{n+60}\sim\sigma_n\mathcal{T}(\nu)+\varepsilon_n$, where $\mathcal{T}(\nu)$ is the student's $t$-distribution with $\nu$ degrees of freedom ($2<\nu\le 59$ is constant) and $\varepsilon_n$ is an error term which should be minimized as much as possible, how can I estimate $\nu$? Formulated otherwise: I wish to find $\nu$ such that the vector $\left\{\varepsilon_n\right\}_{n\ge 0}$ will be as small as possible (in some norm). Are there any methods for achieving this? The length of that data on which I wish to fit the distribution may be arbitrarily long. I'd expect that if the data is normally distributed, the fitting would return 59. In fact, I don't even know how to approach such a problem. The error terms are random themselves, so in what sense can I "minimize" them?


title: Smoothing function for displaying stacked lines without the smoothing introducing crossings

text: I'm displaying time-series data as a "stacked line" or "stacked area" chart. (E.g. with percentage data, data points at 10%, 20% and 30% are displayed at 10%, 30% and 60% on the chart.) Unsmoothed lines through the data points obviously never cross (and never go below the minimum data value or above the maximum total; think 0% and 100% again for instance). I'm looking for a smoothing function that will preserve these properties; that is, the smoothed lines should _never cross_ (including between data points), should _never exceed the maximum/minimum values given by the data points_ , and (ideally) should _only touch at data points_. Does such a beast exist? As is probably apparent, I'm not a statistician. The smoothing here is only for visual appeal, so solutions that are theoretically less than ideal are still definitely on the cards.


title: What is a block in experimental design?

text: I have two questions about the notion of block in experimental design : (1) What is the difference between a block and a factor ? (2) I tried to read some books but something is not clear: it seems that the authors always assume that there is no interaction between the "block factor' and other factors. Is it right, and if it is, why ?


title: How is the proof that the Quartimax/Varimax-rotation converges?

text: Empirically the quartimax-/varimax-rotation has proven useful and it was always converging in my applications. But from my readings years ago I have a vague remembering, that the fact of a proof of convergence has been mentioned but I've never seen the actual proof. Perhaps this is not difficult and may be done here, or does someone have a reference? What I'm really after is the proof for the "inverse" of the quartimax: where the rotation is defined for the _minimizing_ (instead of the maximizing) of the criterion; I hope I can use the quartimax-convergence-proof directly or at least have a path how to proceed with such a proof on my own... * * * [update] Hmm, I've got some hint to articles of ten Berge(1995, 60(3)) and of Jennrich(2001, 66(2)) (both appeared in "Psychometrika" that couple of years ago), which both seem to deal with _problems_ of the convergence of that rotation - so this seems not to be trivial (or as easy as I hoped it would possibly be)... (Unfortunately, I shall not have access to our library before monday)


title: How do I report error from imbalanced data in a random forest algorithm?

text: I have built what I think is a very good predictive model using randomforest. The initial dataset was imbalanced for the outcome 2:1, so I randomly resampled the dataset to balance it, then trimmed the predictors down to 20 or so and managed to get the sensitivity and specificity of the model up to the 90s based on 10-fold cross validation. Can I report that? Do I not have to test it on the imbalanced original dataset? I'm kind of afraid of these results as they look a bit too good, although they did take some man-hours to tune the machine to within an inch of its life. I've seen such things reported in the biomedical literature without a separately sourced validation dataset. Is 10-fold cross-validation "enough"? Essentially I want to make sure I haven't "cheated". Have I inflated the precision by resampling and if so what should I do about it? See below the WEKA buffer output.               === Classifier model (full training set) ===          Random forest of 200 trees, each constructed while considering 5 random features.     Out of bag error: 0.0199                    Time taken to build model: 0.24 seconds          === Stratified cross-validation ===     === Summary ===          Correctly Classified Instances         147               97.351  %     Incorrectly Classified Instances         4                2.649  %     Kappa statistic                          0.947      Mean absolute error                      0.0531     Root mean squared error                  0.1419     Relative absolute error                 10.6329 %     Root relative squared error             28.384  %     Coverage of cases (0.95 level)         100      %     Mean rel. region size (0.95 level)      59.2715 %     Total Number of Instances              151               === Detailed Accuracy By Class ===                         TP Rate   FP Rate   Precision   Recall  F-Measure   ROC Area  Class                      0.973     0.026      0.973     0.973     0.973      0.998    FALSE                      0.974     0.027      0.974     0.974     0.974      0.998    TRUE     Weighted Avg.    0.974     0.027      0.974     0.974     0.974      0.998          === Confusion Matrix ===                 a   |  b   <-- classified as           71    | 2 |  a = FALSE            2    | 76 |  b = TRUE


title: The meaning of convergence in Variational Inference?

text: My friend and I are discussing about the convergence of Variational Inference, especial for Expectation Propagation method. After running some loops, the likelihood of my graphical model can be plotted like the below figure (the x axis is the number of iteration and the y axis is the likelihood value). When the likelihood is convergent(after 15 iterations), the distributions of variables in my probabilistic graphical model are also convergent. My friend says that we can use the value of variable that maximize the likelihood. I do not agree because I think the convergent value should be chosen. If we use the value that maximize the likelihood, why do we need the convergence of inference algorithm? What is the meaning of convergence in Variational Inference? However, I do not have any theory evidence for my opinion. Would you please give me some hints / proof / paper / book to reject or support my point? ![Likelihood over the loop](http://i.stack.imgur.com/QGr7j.jpg)


title: Probability of an SRS having a mean between two values

text: Can anybody help me out with the following question? The distribution of the age at which all married males got married is right skewed with mean = 22.9 years and standard deviation = 1.5 years. What is the probability that an SRS of 75 married males would have a sample mean between 22.8 and 23.0?


title: Calculating and Normalizing ngram relevancy scores from free text extraction

text: I currently look for a set of ngrams in many sets of documents to establish a relevancy score for each set - _eg. I look for the n-gram "adhesive tape" in ~1M sets of 1-500 documents._ The values I have at my disposal are:   * the total documents in each set   * the total documents each ngram was found in within each set   * the total number of times each ngram was found in each document in each set   * the total number of times each ngram was found in each set I've tried using raw occurrences of each ngram, and number of documents each ngram occurs in, and various combo calculations. None of these are really establishing very good relevancy. The problems I'm facing ( among others I'm probably not thinking of ) are:   * different ngrams occur more often than other ngrams - _eg. "adhesive" will simply be used more in english text than "adhesive tape"._   * the number of documents in each set will vary   * the amount of text in each document will vary **What are some statistical methods I could try for this type of relevancy extraction from free text? I have a limited background in this type of semantic analysis.**


title: Whitening data before regression, should I whiten the response variable too?

text: I have some data X where the samples are not independent (they're correlated with each other), and I'm trying to do a regression of some continuous variable y on X. This sample correlation could confound my analysis. Hence, I thought of using ZCA whitening to correct for this inter-sample covariance in X $$ X_{\mbox{zca}} = U D^{-1} U^T X, $$ based on SVD $$ X=UDV^T. $$ In addition to correcting the input data, should I also correct the response variable y (using the same singular values/vectors used to transform X) $$ y_{\mbox{zca}}=UD^{-1}U^Ty $$ and then do the regression of y on X?


title: How to deal with quasi-complete separation in a logistic GLMM?

text: **Update** : Since I now know that my problem is called **quasi-complete separation** I updated the question to reflect this (thanks to Aaron). * * * I have a dataset from an experiment in which 29 human participants (factor `code`) worked on a set of trials and the `response` was either 1 or 0. In addition, we manipulated the materials so that we had three crossed factors, `p.validity` (valid versus invalid), `type` (affirmation versus denial), and `counterexamples` (few versus many):               d.binom <- read.table("http://pastebin.com/raw.php?i=0yDpEri8")     str(d.binom)     ## 'data.frame':   464 obs. of  5 variables:     ##      $ code           : Factor w/ 29 levels "A04C","A14G",..: 1 1 1 1 1 1 1 1 1 1 ...     ##      $ response       : int  1 1 1 1 0 1 1 1 1 1 ...     ##      $ counterexamples: Factor w/ 2 levels "few","many": 2 2 1 1 2 2 2 2 1 1 ...     ##      $ type           : Factor w/ 2 levels "affirmation",..: 1 2 1 2 1 2 1 2 1 2 ...     ##      $ p.validity     : Factor w/ 2 levels "invalid","valid": 1 1 2 2 1 1 2 2 1 1 ...      Overall there is only a small number of 0s:               mean(d.binom$response)     ## [1] 0.9504      One hypothesis is that there is an effect of `validity`, however, preliminary analysis suggests there might be an effect of `counterexamples`. As I have dependent data (each participant worked on all trials) I would like to use a GLMM on the data. Unfortunately, `counterexamples`quasi-completely separate the data (at least for one level):               with(d.binom, table(response, counterexamples))     ##         counterexamples     ## response few many     ##        0   1   22     ##        1 231  210      This is also reflected in the model:               require(lme4)     options(contrasts=c('contr.sum', 'contr.poly'))               m2 <- glmer(response ~ type * p.validity * counterexamples + (1|code),                  data = d.binom, family = binomial)     summary(m2)     ## [output truncated]     ## Fixed effects:     ##                                      Estimate Std. Error z value Pr(>|z|)     ##   (Intercept)                            9.42     831.02    0.01     0.99     ##   type1                                 -1.97     831.02    0.00     1.00     ##   p.validity1                            1.78     831.02    0.00     1.00     ##   counterexamples1                       7.02     831.02    0.01     0.99     ##   type1:p.validity1                      1.97     831.02    0.00     1.00     ##   type1:counterexamples1                -2.16     831.02    0.00     1.00     ##   p.validity1:counterexamples1           2.35     831.02    0.00     1.00     ##   type1:p.validity1:counterexamples1     2.16     831.02    0.00     1.00      The standard errors for the parameters are simply insane. As my final goal is to assess whether or not certain effects are significant, standard errors are not totally unimportant.   * **How can I deal with the quasi complete separation?** What I want is to obtain estimates from which I can judge whether or not a certain effect is significant or not (e.g., using `PRmodcomp` from package `pkrtest`, but this is another step not described here). Approaches using other packages are fine as well.


title: How to calculate the population survival probability in a two stage randomized clinical trial?

text: I am trying to reproduce a paper on two stage randomized clinical trials. In these clinical trials, a patient is randomized to either of the maintenance therapies $B_1$ or $B_2$ upon their remission after getting the induction therapy $A_1$. The remission probability is taken to be 0.80 and the 2nd stage randomization probability (to treatment $B$) is taken to be 0.50. They have generated survival times $(t)$ for patients getting induction therapy $A_1$ by calculating the following quantities.               #Generating remission status:     r<-rbinom(n,1,p1)             #Generating t0 from exponential with mean 182.5 days:     t0<-rep(0,n)     for(i in 1:n){     t0[i]<-ifelse(r[i]==0,rexp(1,1/(365*0.5)),0)     }          #Generating tr from exponential with mean 365 days:     tr<-rep(0,n)     for(i in 1:n){     tr[i]<-ifelse(r[i]==1,rexp(1,1/365),0)     }          #Generating B treatment indicator:     z<-rbinom(n,1,p3)             #Generating t1 from exponential with mean 365 days:     t1<-rep(0,n)     for(i in 1:n){     t1[i]<-ifelse(z[i]==1,rexp(1,1/365),0)     }          #Generating t2 from exponential with mean 547.5 days:     t2<-rep(0,n)     for(i in 1:n){     t2[i]<-ifelse(z[i]==0,rexp(1,1/365),0)     }          #Generating survival time t(i):     t<-rep(0,n)     for(i in 1:n){     t[i]<-(1-r[i])*t0[i]+r[i]*(tr[i]+z[i]*t1[i]+(1-z[i])*t2[i])     }      Now, I can't understand how do I calculate the population survival probability for treatment policy $A_1B_1$. Say, the assessing time is 700. I think I need to calculate the mean/rate from the sum of two Exponentials $tr$ and $t1$ by incorporating the probability of remission and the probability of second randomization. But I am not getting properly how to do it. For your convenience you can see this link.


title: Dummy's, significant and not significant

text: I have in my model two dummy's, for a categorical variable with 3 categories (native, EU-immigrant and non-EU-immigrant). The reference category is 'native'. These are merely control variables. Now I find a significant coefficient for one dummy, and not for the other. See last 2 rows here. Note these are log odds.                   L1PRED  S.E.     Response    R             Fixed Part           cons    -1.775  0.080     female  0.431   0.042     age 0.016   0.001     FA  1.053   0.052     TV  1.778   0.112     second  -0.277  0.052     third   -0.269  0.070     fourth  -0.298  0.093     fifth   -0.163  0.102     EU_imm  -0.011  0.092     nEU_imm 0.953   0.439      I was wondering what the right thing to do here is? Because any control variable that has no effect is usually left out of the model. But I cannot do that here, I think? Would it be most common to recode 'native' and 'EU-immigrants' into one category, since there is no significant effect on the outcome?


title: Compare means of two groups with a variable that has multiple sub-group

text: I was unable to understand how to compare multiple mean values in each group. Within the subgroup the mean values influence each other (If one changes, the other changes automatically). I collected travel distance and the mode of transportation choice at two different locations. There are 32 people in each location with "x" different modes. A person can choose one mode and record their travel distance. How should I compare to see if the travel distances traveled by different modes of these two groups are significantly different? If it is categorical data, I know I have to use chi square, but does the t-test work here? I believe the t-test is to compare two groups with one mean in each group. But I have 2x separate mean values (one for the each mode and two locations).


title: What is the name of this perceptron-like classifier?

text: I wanted to find a variant of the perceptron which works for non-separable data, so I tried using $f(x)=\mathrm{\tanh}(x)$ instead of the hard threshold function and finding a $w$ that minimises the function $$\frac{1}{2} \sum_i (y_i - f(w \cdot x_i))^2$$ where $y_i \in \{-1,1\}$ are the classes of the training examples and $x_i$ are the features. I minimised it by using "batch" gradient descent and it seems to work nicely and gives me what I want. My question is: is this the same as **Adaline**? I found a few references to Adaline, but I can't tell whether Adaline refers to the classifer itself, or also the algorithm which is used to train it, or even the physical machine which originally implemented that algorithm (by the way, I am equally confused about the use of "perceptron"). Am I justified in saying that I am using Adaline?


title: Mean survival time of a Weibull distribution

text: I'm trying to calculate the mean survival time of a Weibull distribution, and am getting what feels like an errant estimate of the mean--and each source I look up for how to calculate the mean gives a slightly different formulation. Going from Klein and Moeschberger, the mean is: $$\frac{\Gamma(1+1/\alpha)}{\lambda^{1/\alpha}}$$ As I understand it, $\alpha$ is the shape parameter, and $\lambda$ is $\exp(\beta_0...\beta_ k)$ for a model with $k$ terms. Is this correct? Thus, for a model with two $\beta$ terms, $\beta_0 = 2.18$ and $\beta_1 = 0.66$ along with an $\alpha$ of 0.88, is the mean survival time, as evaluated by R for a non-integer Gamma function as follows:               b0 <- 2.18     b1 <- 0.66     alpha <- 0.88          mean <- (gamma(1+(1/alpha)))/(exp(2.18+0.66)**(1/alpha))      Correct? That produces a mean of 0.0423, which seems...small, even for such a skewed distribution.


title: Bootstrapping unbalanced clustered data (non-parametric bootstrap)

text: I am trying to figure out how to simulate bootstrap samples from a dataset with unbalanced clusters. The approach I would like to adopt is **non- parametric pairs bootstrap** , which easily allows to maintain the dependence structure of the clusters. Suppose for a moment that data were balanced (e.g., 500 mothers, each with 2 children). The two-level simulation algorithm with `B` iterations would be: For $b = 1,\dots, B$,   * Sample 500 mothers with replacement.   * Sample 2 children _without_ replacement. Hence, both the clusters internal composition is maintained unaltered with respect to the initial sample and the final sample size is equal to the one of the original dataset ($N = 1000$). Now, suppose that some mothers has 3 children. This implies that by adopting the above strategy the final simulated sample in general will not be composed by 1000 observations. To your knowledge, are there statistical issues in this second case? If so, how would you proceed? After having read on Davidson book [1] that the unbalanced clusters case would require more advanced techniques I made an extensive bibliographic research, but I found little or nothing about it in terms of simulation algorithms. > [1] Davidson, A. C., Hinkley, D. V. (1997). _Bootstrap methods and their > applications_. Cambridge University Press.


title: Logistic Regression - Bayesian Approach - Assessing Classification Precision

text: I have recently begun to read about bayesian statistics and I am playing around with the R2WinBUGS package. I'm trying to fit a logistic regression to the spam data (that can be found on the webpage of the elements of statistical learning) using R and WinBUGS. My approach was to first divide the data into 80% training and 20% testing sets. I can fit the model using the 80% set but I dont know how to write WinBugs code to predict on new observations (say the 20% test set) and I wonder if this approach to study model/classification precisicion make sense in a Bayesian Approach?


title: Estimating comparative success of different brochures

text: ## The Real World Problem One of my clients is preparing to send a direct mailer to their subscribed user list, and this statistical challenge came up. Their marketing team has 3 different brochures, and wants to know which brochure gets the highest response rate. They would also like to know if sending the mailer with a hand-written address, on a thick envelope, improves results compared to a normal envelope. Let's assume the following:   * For each brochure $b_i$ ($i = 1,2,3$), a person receiving that brochure _who actually opens it and reads it_ will respond with probability $r_i$, where $r_i$ is the true response rate for that brochure   * The thick, high-quality envelopes have a true open rate of $o_{thick}$ while normal envelopes have an open rate of $o_{normal}$   * From previous mailings, we expect actual observed response rates will be between approximately 1% and 5%.  ## Our Goals We want to find the best brochure while sending the fewest number of mailers. We also want to estimate the two open rates. Upon gathering empirical response rates from actual sent mailers, if the true difference between response rates $r_i$ is greater than half a percent, we should be able to detect that difference as statisitcally significant with $p < .05$ ## My thoughts so far We randomly assign users to each of the 3 brochures, such that $N$ users receive each brochure. We want to know what $N$ we need to achieve our desired sensitivity in detecting differences in response rates. Assuming the worst case, we need to be able to detect a difference between true rates of 1% and 1.5%. The SD for this difference is $\sqrt{\frac{(.01*.99) + (.015*.985)}{N}}$. Setting twice that quantity (2 standard deviations gives us 95% confidence) equal to .005 (our desired half a percent) leads to the solution $N = 3948$. ## Questions   * Is this the optimal design or can we do better?   * Is my calculation of $N$ correct? Finally, what is the best way to estimate $o_{normal}$ and $o_{thick}$, or simply the difference between the two? My idea was to randomly assign half of each brochure group to each type of envelope. Within each brochure group, the observed response rates would be the product of the open rates and the $r_i$. This would complicate my calculation of $N$ above, since really I should have used this product in my calculation. My answer would then depend on an estimate of the average open rate -- $\frac{o_{normal} + o_{thick}}{2}$ -- which I'd have to guess at. Also, I am not sure how to determine the distribution of the difference between $o_{normal}$ and $o_{thick}$, since we now have three different estimates of that difference, each of which depends on a different $r_i$, each of which we have only empirical estimates of, empirical estimates which themselves depend on our guess at the average open rate. Thanks very much for any help with this.


title: How to control the cost of misclassification in Random Forests?

text: **Is it possible to control the cost of misclassification in the R packagerandomForest?** In my own work false negatives (e.g., missing in error that a person may have a disease) are far more costly than false positives. The package rpart allows the user to control misclassification costs by specifying a loss matrix to weight misclassifications differently. Does anything similar exist for `randomForest`? Should I, for instance, use the `classwt` option to control the Gini criterion?


title: Using Fourier analysis to check if data is oscillating

text: A simulation gives the population numbers for every species in the domain per frame. These vary over time and can be quite noisy, for example we might have: Species A: `1, 1, 2, 3, 4, 5, 5, 4, 3, 2, 1` Species B: `5, 4, 3, 3, 2, 2, 2, 1, 2, 3, 4` Fourier analysis gives: Species A               [1]  31.0000000+0i -11.3435375+0i   0.7284459+0i  -0.4486906+0i     [5]   0.6467467+0i   0.4170355+0i   0.4170355+0i   0.6467467+0i     [9]  -0.4486906+0i   0.7284459+0i -11.3435375+0i      Species B               [1] 31.0000000+0.0000000i  7.7083901-1.7455710i  2.5347604+1.2715540i     [4]  0.0859162+0.3689912i  1.8966739-0.8223734i -0.2257406+0.1538824i     [7] -0.2257406-0.1538824i  1.8966739+0.8223734i  0.0859162-0.3689912i     [10]  2.5347604-1.2715540i  7.7083901+1.7455710i      How can I use this to give an idea if they are oscillating? I wish to somehow use this in a genetic algorithm as the fitness function to see how close each dataset is to oscillating.


title: Correlate heuristic metrics with future true performance

text: I am faced with a problem, that I'm pretty sure is a statistical one, but me taking 1 course in probability followed by 1 course in statistics back in university did not prepare me to adequately face it. With that in mind, while a R-based sample code is always helpful and appreciated, what I'm really after is some directions to learn about the problem at large. At the moment I'm at a stage where I'm not really sure what keywords to plug into google :) So here is the set up (fully computerized simulation, even though objects might suggest physical entities).   * I have a bunch of robots   * The robots share power supply (this is the point of optimization)   * Robots are not identical (further complicates optimization)   * These robots perform various tasks   * The tasks have value assigned to them. Once the task is complete by a robot I can tell how much benefit that brought, but I can't know it ahead of time.   * I have a series of heuristics (H1, H2, ..., Hn) that try to approximate the real performance of a task before it is completed. I use these to decide which task a given robot should perform, obviously these are inaccurate. From simulations **ignoring the shared power** supply I have:   * Series representing the true performance of each robot (by letting each robot just do a bunch of tasks, exhaustive enumeration is impossible, so a sampling of possible tasks is used instead)   * Series representing the total system performance (all individual robot performances combined)   * Series corresponding to each heuristic function (a perfect heuristic would be identical to the true performance series, but I don't have one like that) So what I want to establish are the following 2 things:   1. Decide which heuristic functions are better at predicting future performance based on how they did in the simulation window (could be individual functions or combination). The problem here is that it's not a simple combination, since H1 can have values in [0, 1] and H2 in thousands.   2. Decide power allocations % for each robot. In real use I'd like to make sure that I use the guiding heuristics (from #1) that ensures that each robot brings the most benefit. I intend to achieve this by insuring that the robot with the most accurate heuristics and the highest predictions of future performance gets the most power. But how much more power? How should I split it? This is the ultimate problem I need to solve.


title: Condition for Law of Large Numbers, Monte Carlo

text: In some lecture notes I am reading, there is the following; Consider $X_{1},...,X_{n}$, each with pdf $g$ (the instrumental distribution). Our aim is to estimate $E_{f}[h(X)]$ where $h(X)$ is some measurable function and $f$ a pdf (the target distribution). Then: $\large\frac{1}{n}\displaystyle\sum\limits_{i=1}^n\frac{f(X_{i})h(X_{i})}{g(X_{i})} \rightarrow E_{g}\left[\frac{f(X_{i})h(X_{i})}{g(X_{i})}\right] = E_{f}[h(X)]$ as $n\rightarrow \infty$, provided $E_{g}\left[\left|\large\frac{f(X_{i})h(X_{i})}{g(X_{i})}\right|\right] < \infty$ and supp($g$) $\supset$ supp($f.h$). Now, it is the condition on the last line I am not happy with. I understand that the expectation given must be finite to be able to apply the law of large numbers. However, I really do not understand how the last statement comes into play. Is it to ensure that the expectation is finite? Is it to ensure that the series above converges? It makes sense to me that we should want supp($g$) $\supset$ supp($f$), since $g$ is the pdf we are using to aquire information from the target distribution $f$. Thanks for your insight and thoughts.


title: What's so 'moment' about 'moments' of a probability distribution?

text: I KNOW what moments are and how to calculate them and how to use the moment generating function for getting higher order moments. Yes, I know the math. Now that I need to get my statistics knowledge lubricated for work, I thought I might as well ask this question - it's been nagging me for about a few years and back in college no professor knew the answer or would just dismiss the question (honestly). So what does the word "moment" mean in this case? Why this choice of word? It doesn't sound intuitive to me (or I never heard it that way back in college :) Come to think of it I am equally curious with its usage in "moment of inertia" ;) but let's not focus on that for now. So what does a "moment" of a distribution mean and what does it seek to do and why THAT word! :) Why does any one care about moments? At this moment I am feeling otherwise about that moment ;) PS: Yes, I've probably asked a similar question on variance but I do value intuitive understanding over 'look in the book to find out' :)


title: $p(D)$ in Bayesian Statistics

text: Say I have the following obvious Bayesian computation: ![enter image description here](http://i.stack.imgur.com/1Bjsx.png) where $\theta$ is a model parameter that we try to infer and $D$ is **observed data**. I have always understood $p(D)$ to relate to knowledge that we have about the data $D$, but this concept has always been somewhat abstract. Many texts for example, simply ignore $p(D)$. Moreover, quite often, $p(D)$ is referred to as a **normalizing constant**. But if this is the case, **why write $p(D)$**? Is it **always** the case that $p(D)$ is a constant?


title: Scripts for Mixed-Effects Models in S and S-PLUS

text: I couldn't find the scripts for Mixed-Effects Models in S and S-PLUS. The only place is http://stat.bell-labs.com which has code for first two chapters. I'd highly appreciate if someone point out the place to get the scripts for this book. Thanks in advance for your time and cooperation.


title: Is median fairer than mean?

text: I recently read the advice that you should generally use median not mean to eliminate outliers. Example: The following article http://www.amazon.com/Forensic-Science-Introduction-Scientific- Investigative/product-reviews/1420064932/ has 16 reviews at the moment:               review= c(5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 4, 4, 3, 2, 1, 1)     summary(review)  ## "ordinary" summary          Min. 1st Qu.  Median    Mean 3rd Qu.    Max.      1.000   3.750   5.000   4.062   5.000   5.000       Because they use **Mean** the article gets 4 stars but if they would use **Median** it would get 5 stars. Isn't the median a 'fairer' judge? * * * An experiment shows that median's error is always bigger than mean. **Is median worse?** `               library(foreach)          #the overall population of bookjudgments     n<-5     p<-0.5     expected.value<-n*p     peoplesbelieve <-rbinom(10^6,n, p)          #16 ratings made for 100 books     ratings <- foreach(i=1:100, .combine=rbind) %do% sample(peoplesbelieve,16)     stat <- foreach(i=1:100, .combine=rbind) %do% c(mean=mean(ratings[i,]), median=median(ratings[i,]))          #which mean square error is bigger? Mean's or Median's?     meansqrterror.mean<-mean((stat[,"mean"]-expected.value)^2)     meansqrterror.median<-mean((stat[,"median"]-expected.value)^2)          res<-paste("mean MSE",meansqrterror.mean)     res<-paste(res, "| median MSE", meansqrterror.median)     print(res)      `


title: Repeated measures design with measurements from different groups of animals

text: In a repeated measured design we measure a particular variable at different time points from the same subjects. In animal experiments, if animals are sacrificed at every time point to measure a variable, theoretically the measurements at every time point are from different animals, though the measured variable is same. Can such experiment setting be called repeated measure design ?


title: Upper bounds for the copula density?

text: The Fréchet–Hoeffding upper bound applies to the copula distribution function and it is given by $$C(u_1,...,u_d)\leq \min\{u_1,..,u_d\}.$$ Is there a similar (in the sense that it depends on the marginal densities) upper bound for the copula density $c(u_1,...,u_d)$ instead of the CDF? Any reference would be greatly appreciated.