title: Python and statistical analysis tutorials

text: I am slowly learning Python, and I am curious if anyone has come across some easy to follow tutorials or books explaining how to program simple to intermediate statistical methods, from basic descriptive tests to regression analysis and beyond. Any resource is helpful, thanks.


title: How do I visualize data where I have to compare multiple data, but one gigantic output is tanking others

text: I have 6 different data sets which I need to compare. I am able to compare five of them by plotting graphs, but when I include the sixth dataset, it tanks the other datasets because while the maximum value of the other five datasets is in thousands, for the sixth dataset, the maximum values is in billions. So I get one straight line for the first five datasets and some graph for the sixth. Therefore, I am unable to compare the first five together with the sixth. How should I visualize the data so that I can differentiate between them.


title: How do you choose a unit of analysis (level of aggregation) in a time series?

text: If you can measure a time series of observations at any level of precision in time, and your goal of the study is to identify a relationship between X and Y, is there any empirical justification for choosing a specific level of aggregation over another, or should the choice be simply taken based on theory and/or practical limitations? I have three sub-questions to this main one:   1. Is any non-random variation in X or Y within a larger level sufficient reasoning to choose a smaller level of aggregation (where non-random is any temporal pattern of the observations)?   2. Is any variation in the relationship between X and Y at a smaller level of aggregation sufficient reasoning to justify the smaller unit of analysis? If some variation is acceptable how does one decide how much variation is too much?   3. Can people cite arguments they feel are compelling/well defined for one unit of analysis over another, either for empirical reasons or for theoretical reasons? I am well aware of the modifiable area unit problem in spatial analysis (Openshaw 1984). I don't claim to be expert on the material, but all I am to think so far of it is that a smaller unit of analysis is always better, as one is less likely to commit an ecological fallacy (Robinson 1950). If one has a directly pertinent reference or answer concerning aggregation geographical units I would appreciate that answer as well.


title: Determinant perturbation approximation

text: The following problem comes from a max likelihood calculation for gaussian families, but is of independent interest. Is it possible to find a closed-form approximation for small values of $x$ for $\text{det}(B + xI)$ where I is the identity matrix and B is hermitian rank-deficient positive semidefinite?


title: What should be the optimal parameters for Random Forest classifier?

text: Currently i am using RF toolbox on MATLAB for a binary classification Problem Data Set: 50000 samples and more than 250 features So what should be the number of trees and randomly selected feature on each split to grow the trees? can any other parameter greatly affect the results?


title: Noise in the Fourier transform

text: I have a smooth but rather complex curve, sampled with a good frequency. I apply the discrete Fourier transform to it using the fast Fourier transform (FFT) algorithm and get its Fourier image. I need to find peaks on the resulting Fourier transform curve, but the image I get contains a substantial amount of noise, which is a real problem because peaks can't be clearly seen. So, my question is as follows: What are the sources of noise in the Fourier transform, and how it can be reduced? Typical plots look like this: Initial function ![Initial function](http://i.stack.imgur.com/nCb7B.png) Initial function, narrow range ![Initial function, narrow range](http://i.stack.imgur.com/vk3ld.png) Fourier transform ![Fourier transform](http://i.stack.imgur.com/3HpdW.png) Fourier transform, narrow range ![Fourier transform, narrow range](http://i.stack.imgur.com/DuwNL.png) Before transformation a constant was subtracted from the the function so that it goes to zero.


title: Weird residuals in linear regression

text: I analyze a set of multivariate measurements. It is known that several pairs of independent variables show high linear correlation. The graph below shows a scatterplot of one such pair (X and Y, upper pane), as well the residuals as a function of Y (lower left pane) and the histogram of these residuals (lower right pane) ![data example](http://i.stack.imgur.com/7L5pf.png) As one can see, there is a strange peak in the residual histogram. Many of the remaining linearly-dependent variable pairs from the same data set have a similar peak. I have double checked, and I'm sure there are no duplicate records in the data set. What might be the reason for such a behaviour? PS Please don't ask me to elaborate on the problem domain, I'm not allowed to.


title: Interpreting logistic regression output in R

text: I'm working on a multiple logistic regression in R using `glm`. The predictor variables are continuous and categorical. An extract of the summary of the model shows the following:               Coefficients:                    Estimate Std. Error z value Pr(>|z|)     (Intercept)   2.451e+00  2.439e+00   1.005   0.3150     Age           5.747e-02  3.466e-02   1.658   0.0973 .     BMI          -7.750e-02  7.090e-02  -1.093   0.2743     ...     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1      Confidence intervals:                                 2.5 %       97.5 %     (Intercept)  0.10969506 1.863217e+03     Age          0.99565783 1.142627e+00     BMI          0.80089276 1.064256e+00     ...      Odd ratios:                                Estimate Std. Error   z value Pr(>|z|)     (Intercept)  1.159642e+01  11.464683 2.7310435 1.370327     Age          1.059155e+00   1.035269 5.2491658 1.102195     B            9.254228e-01   1.073477 0.3351730 1.315670     ...      The first output shows that $Age$ is significant. However, the confidence interval for $Age$ includes the value 1 and the odds ratio for $Age$ is very close to 1. What does the significant p-value from the first output mean? Is $Age$ a predictor of the outcome or not?


title: Truncated distributions, properties

text: If one truncates a distribution, for example like the zero truncated poisson distribution: https://en.wikipedia.org/wiki/Zero-truncated_Poisson_distribution or more general: https://en.wikipedia.org/wiki/Truncated_distribution What are the properties one is using when doing this? Do the random variables just have to be iid (Independent and identically distributed). Or are there other properties one needs? Please explain someone (not too formal) why this can be done? And also: Why is p for all x in the truncated distribution rising the same factor?


title: What is behind Google Prediction API?

text: Google Prediction API is a cloud service where user can submit some training data to train some mysterious classifier and later ask it to classify incoming data, for instance to implement spam filters or predict user preferences. But what is behind the scenes?


title: Interaction between ordinal and categorical factor

text: What analyses can be used to find an interaction effect in a 2-factor design, with one ordinal and one categorical factor, with binary-valued data? Specifically, are there any types of analyses that are capable of dealing with a 2 factor design 5(ordinal) x 2(categorical), where the outcomes are either true or false? One could do a 5x2 chi square analysis, but it loses the power of the ordinality of the one factor. Alternatively, one could run independent logistic/probit regressions, but then there is the question of testing for the interaction effect. Any thoughts or suggestions that would put me in the right direction would be helpful.


title: Clustered standard errors and pooled OLS

text: Can someone please explain why clustered standard errors ensure that the error estimates of Pooled OLS results are appropriate?


title: the paired samples

text: Trace metals in drinking water affect the flavor of the water, and unusually high concentration can pose a health hazard. The following Table shows trace- metal concentrations (zinc, in mg/L) for both surface water and bottom water at six different river locations. Our aim is to see if surface water concentration (x) is predictive of bottom water concentration(y).               Location Bottom Surface      1    0.430  0.415      2    0.266  0.238      3    0.567  0.390      4    0.531  0.410      5    0.707  0.605      6    0.716  0.609       a) Draw a scatter diagram to show a possible association between the concentrations and check to see if a linear model is justified. b) Test to see if the two concentrations are independent; state your hypotheses and choice of test size.


title: Sampling of beta in bayesian regression (variable selection)

text: I am sampling a beta using a Gibbs sampling. It is variable selection model. So in different iteration of gibbs different covariates are included to the model (denoted by a variable selection indicator $\gamma$). However, this indicator have high variance ( this is because my covariates are correlated). I use: $\hat\beta = 1/n\sum \beta_{samples} $ as my estimate for further analysis. This turn out to be problematic because if I calculate the explained variance it is negative. On the other hand if I just use a single sample of $\beta_{sample}$ to find explained variance it is correct. What's wrong here? Do I need to re-estimate $\beta$ for nonzero mean $\gamma = 1/n \sum \gamma_{sample}$? This does not sound correct because this estimate does not correspond to mean of posterior of $\beta$. thanks


title: Selecting uncorrelated samples from a set of bulk data that contains correlated and dependent samples

text: i have a set of data that is generated by expensive computational model evaluations, on a total data set of 10000 samples in 40 dimensions. This sample data set is composed of different data sets, originating partly from random runs, latin hypercube DOE, radial design DOE, linear parameter studies, and a large part is based on the history data generated by several optimization runs using genetic algorithms. My thought was that a large part of the function evaluations generated during the genetic algorithms runs, could be some how used to augment them to the set of random and latin hypercube samples, in order to have a larger sample set to perform a variance based sensitivity analysis. I came up with 2 ideas, but i am an engineer, not a mathematician: 1) using the covariance matrix for the total samples matrix, trying to filter out samples until the of diagonal terms are smaller then some threshold, to avoid correlations. 2)The other idea was to make some sort of minimum distance filter to avoid areas with tightly clustered samples. Would that be sufficient? are there any tests for randomness, that i could use? The problem is that i don't know the right terminology, so maybe there exist ready to run methods for such problems, but i don't know how to find them, because i don't know their names. I am thankful for any helpful suggestions.


title: Control parameter minsplit for rpart in regression tree

text: I processed rpart( ) on the same dataset. One did not use the control parameter "minsplit", but the other one did. I do not understand why I got the different first node in two processes. My understanding is that the "minsplit" will only extend the branches of the tree, but it should not change the structure of the tree indicated by the significant predictors. Please correct me if my understanding is not right and explain the reason of two different results. Thanks! 1st method - without "minsplit": > rt <\- rpart(SO2~Temp+Manuf+Pop+Wind+Precip+Days, data=usair, > method="anova") par(xpd=NA) plot(rt) text(rt, use.n=TRUE, all=TRUE) ![enter > image description here](http://i.stack.imgur.com/77J4V.png) 2nd method - with "minsplit":               rt <- rpart(SO2~Temp+Manuf+Pop+Wind+Precip+Days, minsplit=10, data=usair, method="anova")     par(xpd=NA)     plot(rt)     text(rt, use.n=TRUE, all=TRUE)      ![enter image description here](http://i.stack.imgur.com/rncB4.png)


title: Estimating the mean with least median of squares

text: I have a set of real numbers $x_{1}...x_{n}$ and would like to estimate their mean (let's call it y) so that $median(x_{i}-y)^{2}$ is minimal. Is there an algorithm and a correctness proof for finding it? I understand I need to look at segments of the type $[x_{i},x_{i+n/2}]$ and check for the shortest one. But how can I prove that? Thank you


title: How to calculate "individual precision" of a data set?

text: Looking at the Wikipedia article on _Conjugate prior_ , in the section on _Continuous distributions_ , We see the first two entries discuss "Normal with known variance" and "Normal with known precision". The "Interpretation of hyperparameters" column refers to the "sum of individual precisions". While I understand precision is the reciprocal of variance, I was under the impression variance was a characteristic of the sample population and not of individual data points. Is My understanding incorrect? If so, how do I calculate the "individual precision"?


title: Martingale and deterministic functions

text: Suppose: $u_t \sim N(0,1) \ iid.$, $X_t = g(X_{t-1}) \cdot u_t$ whereas $g(X)$ can be any deterministic function. Is this sufficient to define a martingale? So does it hold: $E(X_t|X_{t-1}, \ldots , X_1) = E(g(X_{t-1}) \cdot u_t|X_{t-1}, \ldots , X_1) =^{\ pi.} E(g(X_{t-1})|X_{t-1}, \ldots , X_1) \cdot E(u_t|X_{t-1}, \ldots , X_1) = E(g(X_{t-1})|X_{t-1}, \ldots , X_1) \cdot 0 = 0$ It seems pretty obvious but I think I missed some assumptions for $g(X)$ or am I totally wrong?


title: Studying whether a policy had an effect outside the normal fluctuations in a time series, with time series cross sectional data

text: I'd really appreciate a point in the right direction of what models I should study, and the best modelling approach to evaluate a sentencing policy. I certainly am not asking for a full solution, but it would save me a huge amount of time if I knew where to start. I have experience of time series models with economic time series (including concepts like nonstationarity and cointegration), but no experience at all with time series cross sectional data. **Background** Sentence lengths vary through time for many reasons such as changes in judicial attitudes, changes in the severity of offending coming before courts, public perceptions of sentencing etc. One important reason why sentencing levels change is the release of new/revised sentencing guidelines, which are documents that advise judges on appropriate sentence lengths. I would like to study whether a new sentencing guideline had an effect on sentence lengths. That is, following the release of a sentencing guideline, did sentencing lengths change more than is consistent with the normal month to month fluctuations in sentencing that happen irrespective of changes to guidelines? **Dataset** I have a dataset that records, for each sentence passed, the sentence length, the date of sentence but also lots of additional information that can be used as controls, such as the aggravating and mitigating factors present in each case. The dataset runs from a year before the guideline was introduced to around a year after its introduction. For any particular offence, I may have around 5,000 records a year (total n=10,000). I also have data that may allow me to control for some common shocks which may hit the system, such as a change in judicial attitudes. e.g., if judges become more penal. This may be possible by using data on average sentence lengths for offences _not_ affected by new sentencing guidelines (the guideline of interest in on assault, but I also have data on, say, burglary). As far as I can tell, charts of average sentence lengths appear almost nonstationary - they wander up and down but do not seem drawn to any long run average or trend . Intuition tells me sentence lengths cannot wander indefinitely up or down. However, over the time horizons I'm studying, it may be reasonable to assume they can wander without bounds. **Attempts so far** I have grouped sentences into months and calculated in each month average sentence lengths. I have then estimated simple ARDL models. This allows me to generate a forecast of the time series beginning just prior to the introduction of the new guideline. I can also plot error bands around this forecast to generate a 'fan chart'. If the outturn of average sentence lengths strays outside the 5th and 95th percentiles of the forecast, then I conclude that the sentencing guidelines probably had an effect. I have included average sentence lengths for other offences (e.g. burglary) in these models in the hope that this captures some of the common shocks hitting the system (e.g. a general increase in penality). However, as yet, I haven't included the controls I have for the aggravating and mitigating factors in each case, because I don't know how to best combine this cross sectional data with the time series data. By doing so, I would hope to "narrow" the level of uncertainty reflected in the fan chart because I would effectively be explaining away some of the sources of variability in sentencing. I would therefore have greater power in distinguishing whether the new guideline had an affect. I feel like a way forward may be to run cross sectional regressions of the aggravating and mitigating factors on sentence length, save the residuals, and use the monthly averages of the residuals rather than average sentence lengths in the time series work. There must be a 'best way' to do this kind of thing, but I'm unsure of what models are most appropriate. Any ideas or help is greatly appreciated - thank you. Apologies for such a lengthy question. Robin


title: What's the formula of normalized correlation?

text: I read a paper which used normalized correlation to evaluate the distance between two vectors. But I searched on the Internet and found little about normalized correlations, but I still got some clues. In one paper the formula for normalized correlation is given as follows: $$dc(y_1,y_2) = \frac{y_1^T y_2}{|y_1||y_2|}$$ I was confused that it's just the cosine similarity formula! So what on earth is the formula of normalized correlation? Thanks in advance! PS: I found some explanation in this article The Corrected Normalized Correlation Coefficient: A Novel Way Of Matching Score Calculation for LDA- Based Face Verification $$ \delta_{cos}(y,\bar{y_{j}}) = \frac {y^{T} \bar{y_{j}}} { \|y\|\|\bar{y_{j}}\|} $$ In statistics, the expression above is often referred to as the normalized correlation coefficient and is used for measuring the extent to which two samples, in our case the vectors $y$ and $\bar{y_{j}}$ are linearly related. When the absolute value of the normalized correlation coefficient equals one, then there exists a linear relation between the two samples, while on the other hand, when the value of the normalized correlation coefficient equals zero, then the two samples have no linear relation. Generally, the higher the absolute value of the coefficient, the stronger the linear relation between the two samples. Based on this fact the absolute value of the normalized correlation coefficient is commonly employed for computing the matching score between the input vector $y$ in the client template $\bar{y_{j}}$


title: If $Z$ is standard normal, and $Z^2$ is chi-squared, is $Z/(\sqrt{Z^2/n})$ a t distribution?

text: If $Z$ is $\mathcal N(0, 1)$, and $Z^2$ is $\chi^2_1$, is $\frac Z {\sqrt{Z^2/n}}$ a $t$ distribution with $n=1$ degree of freedom?


title: Interval Censored Survival Analysis

text: I am conducting a study on neurological diseased patients. I am looking into possible factors related to time-to-dementia. Costs associated with neurological testing allow us to assess all patients in the prospective cohort on an annual basis. We have (can afford) 4 fixed follow up visits for all patients - baseline, 1-year, 2-year, 3-year. We assess dementia at each visit. My data look roughly as follows:   * Baseline: 0/140 patients with dementia.   * 1-year: 10/120 patients with dementia.   * 2-year: 20/100 patients with dementia.   * 3-year: 30/80 patients with dementia.  A lot of censoring is occurring year-to-year. Further, I don't know when patients actually develop/convert to dementia. Just know it occurred between:   * (0,1]   * (1,2]   * (2,3] A sensible looking reference on "fixed interval" survival data seemed to suggest that I could analyze these data (after manipulating the dataset) using GLM (binomial distribution, and complementary log-log link). Also compared to Cox PH using Efron method to account for ties. Thoughts on this approach? http://www.ics.uci.edu/~vqnguyen/stat255/Lecture13.pdf If I want a "regression" approach to investigate time to dementia what might be possible alternatives to the above model(s)?


title: Statistics Modelling Question

text: Here is the situation: I want to sell my house. The price I'm getting from people who want to by my house follows i.i.d. with $X_n \sim \text{Uniform}(0,1)$, where $X_n$ is the highest offer on the $n$th day. I spend $C<0.1$ for advertising my house. My overall profit is $X_n-C_n$ when I sell the house on the $n$th day with price $X_n$. What kind of price should I consider to sell the house? I have the below so far: $E(X) = 0.5$ First day return = $X_1-C$ Second day return = $X_2-2C$ ... the earlier $n$(th day) is, the higher the return, so $n$ should be $<20$ (if $C=0.05$) $X_n$ is decreasing as $n$ increases, so the value should be $< 0.5$? I am thinking that it is not possible to calculate the exact value? So the answer should be expressed in term of $C$?


title: How to interpret p-value of Kolmogorov-Smirnov test (python)?

text: I have Two samples that I want to test (using python) if they are drawn from the same distribution. To do that I use the statistical function ks_2samp from scipy.stats. It returns 2 values and I find difficulties how to interpret them. Help please!


title: Comparing a variable based on category in R

text: I have a data like this:               ds <- structure(list(return = c(-1.11777941176471, -0.883516129032258,            3.71416129032258, 0.028483870967742, -2.05184375, -0.26984375,            5.07122058823529, 15.54915625, 0.557730769230769, 3.40277358490566           ), factor_5 = c(0L, 1L,  1L, 0L, 0L, 1L, 1L, 1L, 0L, 0L)),            .Names = c("return", "factor_5"), row.names = c(NA, 10L), class = "data.frame")      These are nothing but stock returns and their classification. Now I need to test which returns are better - the returns of stock belonging to category 1 or 0. The mean, skewness and all other simple statistics are more or less the same and the density plots are also almost same. But now I need to say which category yields better returns. Is there any other test that I can carry out to do this?


title: Parallel lines on residual vs fitted plot

text: I have a multiple regression problem, which I tried to solve using simple multiple regression:               model1 <- lm(Y ~ X1 + X2 + X3 + X4 + X5, data=data)      This seems to be explaining the 85% of variance (according to R-squared) which seems pretty good. However what worries me is the weird looking Residuals vs Fitted plot, see below: ![enter image description here](http://i.stack.imgur.com/Y2UBn.png) I suspect the reason why we have such parallel lines is because the Y value has only 10 unique values corresponding to about 160 of X values. Perhaps I should use a different type of regression in this case? **Edit** : I've seen in the following paper a similar behavior. Note it's a one-page only paper so when you preview it you can read it all. I think it explains pretty well why I observe this behavior but I'm still not sure if any other regression would work better here? **Edit2:** The closest example to our case I can think of is the change in interest rates. FED announces new interest rate every few months (we don't know when and how often). In the meantime we gather our independent variables on the daily basis (such as daily inflation rate, stock market data, etc.). As a result we will have a situation where we can have many measurements for one interest rate.


title: Confusion between sample and population

text: > A team has conducted a study in which an introductory class was randomly > divided into two groups. Group 1 was administered a dosage of alcohol. Group > 2 was given an equvialent dosage of marijuana. Fifteen minutes after > administration, both groups were asked to solve a puzzle. The experimenter > recorded the amount of time in seconds it took for each subject to solve the > puzzle. Now, in this case, will the given data for each group be considered as sample or population?


title: aggregating multiple shortest-path networks

text: I have a dataset of 500 observations (firms) and for each firm I know how they are directly and indirectly related based on some characteristics. I can therefore construct a network of connectivity across all stocks. For each firm I compute the shortest path to all other firms. I end up with 500 networks of degree of separation. For instance firm A is directly related to B and C and indirectly connected to D and E through B. Firm B and C belong to the first set of degree separation and D and E to the second degree of separation. How can I aggregate all these 500 networks into one? Is there a mathematical / statistical methodology where I can aggregate all these networks? For instance, I would construct a network with a node (say A) representing a group of firms that are all directly connected and this node A is connected to another node, B. Firms in node B are also firms that directly connected but what connects B and A is that firms in A and B are often connected at the second degree of separation.


title: Life table probabilities

text: I am using 2007's Social Security Life Table: http://www.ssa.gov/oact/STATS/table4c6.html I would like to use the table to generate a number for each scenario I ran. This number represents the remaining life (in years) scenario n will live. The way I plan to do this is to first input an age 'x' let say age 60. I get the corresponding number on column 'Number of Lives' which is 85,227. (This simply reveals that out of 100K people, 85,227 survived) Next, I randomly select one integer from 1:100K and if it is less than 85,227, then this scenario will continue on to next year else this case will die off. I repeat this process, at each age get the corresponding number from column 'Number of Lives' and so on.. I want to know if such method is sound from a probability perspective. Does it violate any axioms. Do I still preserve the underlying population distribution of the table? Any better ways of doing it? Thanks,


title: Time Series Econometrics: VECM with multiple cointegration vectors

text: Has anyone gone through some papers using Vector Error Correction Models in causality applications with more than one cointegration vectors, say two. I guess there will be more than one ECM terms. How to assess the endogeneity of the left hand variables if t-stats on different ECM coefficients yield different (conflicting) results. For left hand side variables to be endogenous should we have both ECM coefficients negative and significant? Javed Iqbal


title: Multiple hypothesis testing correction with Benjamini-Hochberg, p-values or q-values?

text: Given a list of p-values generated from independent tests, sorted in ascending order, one can use the Benjamini-Hochberg procedure for multiple testing correction. For each p-value, the Benjamini-Hochberg procedure allows you to calculate the False Discovery Rate (FDR) for each of the p-values. That is, at each "position" in the sorted list of p-values, it will tell you what proportion of those are likely to be false rejections of the null hypothesis. My question is, are these FDR values to be referred to as " **q-values** ", or as " **corrected p-values** ", or as something else entirely? **EDIT 2010-07-12:** I would like to more fully describe the correction procedure we are using. First, we sort the test results in increasing order by their un-corrected original p-value. Then, we iterate over the list, calculating what I have been interpreting as "the FDR expected if we were to reject the null hypothesis for this and all tests prior in the list," using the B-H correction, with an alpha equal to the observed, un-corrected p-value for the respective iteration. We then take, as what we've been calling our "q-value", the maximum of the previously corrected value (FDR at iteration i - 1) or the current value (at i), to preserve monotonicity. Below is some Python code which represents this procedure:               def calc_benjamini_hochberg_corrections(p_values, num_total_tests):         """         Calculates the Benjamini-Hochberg correction for multiple hypothesis         testing from a list of p-values *sorted in ascending order*.              See         http://en.wikipedia.org/wiki/False_discovery_rate#Independent_tests         for more detail on the theory behind the correction.              **NOTE:** This is a generator, not a function. It will yield values         until all calculations have completed.              :Parameters:         - `p_values`: a list or iterable of p-values sorted in ascending           order         - `num_total_tests`: the total number of tests (p-values)              """         prev_bh_value = 0         for i, p_value in enumerate(p_values):             bh_value = p_value * num_total_tests / (i + 1)             # Sometimes this correction can give values greater than 1,             # so we set those values at 1             bh_value = min(bh_value, 1)                  # To preserve monotonicity in the values, we take the             # maximum of the previous value or this one, so that we             # don't yield a value less than the previous.             bh_value = max(bh_value, prev_bh_value)             prev_bh_value = bh_value             yield bh_value


title: Correct formatting (capitalization, italicization, hyphenation) of "p-value"?

text: I realize this is pedantic and trite, but as a researcher in a field outside of statistics, with limited formal education in statistics, I always wonder if I'm writing "p-value" correctly. Specifically:   1. Is the "p" supposed to be capitalized?   2. Is the "p" supposed to be italicized? (Or in mathematical font, in TeX?)   3. Is there supposed to be a hyphen between "p" and "value"?   4. Alternatively, is there no "proper" way of writing "p-value" at all, and any dolt will understand what I mean if I just place "p" next to "value" in some permutation of these options?


title: How to find text blocks in a scanned document?

text: I am trying to detect text in a scanned document by examining variations in the lightness of the scan collapsed vertically. Here's a sample of the input I would receive, with the lightness plot of each vertical pixel strip superimposed: ![Example](http://i.stack.imgur.com/coUW6.jpg) Note: I've applied a Gaussian smoothing function to the data ~ 10 times, but it seems to be pretty wiggly to begin with. It is easy to see that the left margin is _really_ wiggly (i.e., has many extrema). **Problem** : I want to generate a set of critical points of the image. I've resorted to computing the number of extrema of the function within an interval (using the derivative and its proximity to zero) and dividing that by the length of the interval, but that isn't easy on the computer. (I use Python, and I couldn't find many low-pass filters for the data.) Thanks!


title: Nonlinear models which are hard to estimate

text: Genetic algorithms are avoided in econometry literature as often as possible, but still sometimes they are inevitable. The question is: **Which well known models are the most difficult to estimate using conventional algorithms?** (By Conventional algorithms I mean Gauss-Newton method, Levenberg–Marquardt algorithm, and so on) **Motivation:** I want to test some heuristic methods, and I need some benchmark to be sure that this particular model is really hard to estimate.


title: Multivariate orthogonal regression in R

text: I have a project in which I need to perform orthogonal regression in a multivariate space. For the univariate case, I've found Teetor's R Cookbook suggests using principle components:               reg_orth = prcomp( ~ y + x, data=ds[train,])     reg_orth_slope = reg_orth$rotation[2,1]/reg_orth$rotation[1,1]     reg_orth_int = reg_orth$center[2]- reg_orth_slope * reg_orth$center[1]     reg_orth_pred = ds[test,'y'] * reg_orth_slope + reg_orth_int      This is fine and seems to work great, but I actually have 5 independent variables that I want to use to predict my independent variable. I tried to consider what is done to compute the slope and intercept, above, but with 5 independent variables, I'm not sure what to do. I would have thought orthogonal regresssion would be quite common, but I'm not finding any information on it in the multivariate case.


title: Intuitive reason why the Fisher Information of Binomial is inversely proportional to $p(1-p)$

text: It confuses/blows my mind that the Binomial has variance proportional to $p(1-p)$. Equivalently, the Fisher information is proportional to $\frac{1}{p(1-p)}$. What is the reason for this? Why is the Fisher Information minimized at $p=0.5$? That is, why is inference most difficult at $p=0.5$? ### Context: I'm working on a sample size calculator, and the formula for $N$, the sample size needed, is an increasing factor of $p(1-p)$, the result of a variance estimation in the derivation.


title: Using Ito's lemma for Brownian motion

text: I am a little confused by Ito's lemma. I reviewed the basic application for geometric brownian motion. I'm now trying to apply it to a different functional form to make myself better. My mind leaps to the following: \begin{equation} dx=x (\mu dt+(\sigma+\delta^t) dz) \end{equation} \begin{equation*} \delta<1 \end{equation*} \begin{equation*} dz=\epsilon \sqrt{dt}, \epsilon = N(0,1) \end{equation*} \begin{equation*} F=ln(x) \end{equation*} I tried to put my work below: \begin{equation*} dF \approx \frac{dF}{dt}dt+\frac{dF}{dx}dx+\frac{1}{2}\frac{d^2F}{d^2dx}dx^2 \end{equation*} \begin{equation*} dF \approx 0 +\frac{1}{x}x(\mu dt+(\sigma+\delta^t)dz)-\frac{1}{2}\frac{1}{x^2}x^2(\mu dt+(\sigma+\delta^t)dz)^2 \end{equation*} \begin{equation*} dF \approx \mu dt+(\sigma+\delta^t)dz-\frac{1}{2} (\mu^2 dt^2+(\sigma+\delta^t)^2 dz^2+ 2 \mu (\sigma+\delta^t)dz dt) \end{equation*} \begin{equation} dF \approx \mu dt+(\sigma+\delta^t)dz-\frac{1}{2}(\sigma+\delta^t)^2dz^2 \end{equation} I'm wondering, is it true that df/dt=0 now that the functional form is different? There's clearly a t inside x at this point, so F(x(t)). Previously, df/dt was zero, and I thought it was because Brownian motion isn't time- differentiable. But here, I think that it df/dt might be $\frac{1}{x}*x*\delta^t*ln(\delta)*dz$ by the chain rule. Is df/dt still 0 for this functional form? Why/Why not?


title: Marginal calculation in linear gaussian model

text: If we have a linear gaussian model where we have, $P(Z) = \mathcal{N}(\mu_z,\sigma^2_z)$ and $P(X|Z) = \mathcal{N}(w_o+zw,\sigma^2_x)$ How do we calculate the marginal $P(X) = \int_Z P(Z)P(X|Z)$ ?


title: Statistic model for removing bad data

text: Let's say I have a (ads live days, revenue) data set.   The data set shows how much revenues each ads generates during the days it is live.     * ads1 generates 100 dollars during the 5 days when it is live.    * ads2 generates 200 dollars during the 10 days when it is live.    * ads3 generates 300 dollars during the 15 days when it is live. ...   so I have (100,5) , (200,10) , (300,15) .....   Imagine I have 100 ads in the data set and I can only launch 80 of them next month.   I have to decide which ads perform better so I can continue launching next month. What is the statistic model?   I was thinking using linear regression, to keep the nodes that have less variance .... but I am not sure if that is correct.


title: Likelihood function of DSGE model using Kalman filter

text: In Frank Schorfheide's class notes on likelihood functions of DSGE models, he expresses the value of the likelihood function for a given vector of parameters $\theta$, and time series $Y^T$ as: $$p(Y^{T}|\theta)=(2\pi)^{-nT/2}\left(\prod_{t=1}^{T}\left|F_{t|t-1}\right|\right)^{-1/2}exp\\{-\frac{1}{2}\sum_{t=1}^{T}v_{t}F_{t|t-1}v_t\prime\\}$$ where $v_t$ is the innovation in $y$ $$v_t=y_t-\hat{y}_{t|t-1}$$ and the marginal distribution of $y_t$ is $$y_t|Y^{t-1}\sim\mathcal{N}\left(\hat{y}_{t|t-1}, F_{t|t-1}\right)$$ I've just got a few questions about what these terms look like. First, does anyone have an idea what $n$ is in the exponent of the first term in the first equation? I think it might be a misprint, but I'm not sure. Second, what does $F_{t|t-1}$ look like? For an $n\times 1$ vector $y$ I'm picturing an $n\times n$ matrix, but what would the values of $F_{j,k}$ be equal to? I'm picturing the covariance between $y_{t,j}$ and $\hat{y}_{t|t-1,k}$ - is that correct? Lastly, I'm assuming from the results of my code that the value the likelihood function returns is a scalar, but it doesn't look like the formula produces one -- for an $n\times 1$ vector $y$, wouldn't the second term in the first equation be $n\times n$? Or do you think it's meant to be the determinant of $F$?


title: Sampling from bivariate distribution with known density using MCMC

text: I tried to simulate from a bivariate density $p(x,y)$ using Metropolis algorithms in R and had no luck. The density can be expressed as $p(y|x)p(x)$, where $p(x)$ is Singh-Maddala distribution $p(x)=\dfrac{aq x^{a-1}}{b^a (1 + (\frac{x}{b})^a)^{1+q}}$ with parameters $a$, $q$, $b$, and $p(y|x)$ is log-normal with log-mean as a fraction of $x$, and log-sd a constant. To test whether my sample is the one I want, I looked at the marginal density of $x$, which should be $p(x)$. I tried different Metropolis algorithms from R packages MCMCpack, mcmc and dream. I discarded burn-in, used thinning, used samples with size up to million, but the resulting marginal density was never the one I supplied. Here is the final edition of my code I used:               logvrls <- function(x,el,sdlog,a,scl,q.arg) {         if(x[2]>0) {              dlnorm(x[1],meanlog=el*log(x[2]),sdlog=sdlog,log=TRUE)+              dsinmad(x[2],a=a,scale=scl,q.arg=q.arg,log=TRUE)         }         else -Inf         }          a <- 1.35     q <- 3.3     scale <- 10/gamma(1 + 1/a)/gamma(q - 1/a)*  gamma(q)           Initvrls <- function(pars,nseq,meanlog,sdlog,a,scale,q) {         cbind(rlnorm(nseq,meanlog,sdlog),rsinmad(nseq,a,scale,q))     }          library(dream)     aa <- dream(logvrls,             func.type="logposterior.density",             pars=list(c(0,Inf),c(0,Inf)),             FUN.pars=list(el=0.2,sdlog=0.2,a=a,scl=scale,q.arg=q),             INIT=Initvrls,             INIT.pars=list(meanlog=1,sdlog=0.1,a=a,scale=scale,q=q),             control=list(nseq=3,thin.t=10)             )      I've settled on dream package, since it samples until the convergence. I've tested whether I have the correct results in three ways. Using KS statistic, comparing quantiles, and estimating the parameters of Singh-Maddala distribution with maximum likelihood from the resulting sample:               ks.test(as.numeric(aa$Seq[[2]][,2]),psinmad,a=a,scale=scale,q.arg=q)          lsinmad <- function(x,sample)         sum(dsinmad(sample,a=x[1],scale=x[2],q.arg=x[3],log=TRUE))      optim(c(2,20,2),lsinmad,method="BFGS",sample=aa$Seq[[1]][,2])           qq <- eq(0.025,.975,by=0.025)         tst <- cbind(qq,                   sapply(aa$Seq,function(l)round(quantile(l[,2],qq),3)),                   round(qsinmad(qq,a,scale,q),3))      colnames(tst) <- c("Quantile","S1","S2","S3","True")           library(ggplot2)      qplot(x=Quantile,y=value,            data=melt(data.frame(tst),id=1),             colour=variable,group=variable,geom="line")      When I look at the results of these comparisons, KS statistic almost always rejects the null hypothesis that sample is from Singh-Maddala distribution with supplied parameters. Maximum likelihood estimated parameters sometimes comes close to its true values, but usually too far out of comfort zone, to accept that sampling procedure was succesfull. Ditto for the quantiles, empirical quantiles are not too far, but too far away. My question is what I am doing wrong? My own hypotheses:   1. MCMC is not appropriate for this type of sampling   2. MCMC cannot converge, due to theoretical reasons (the distribution function does not satisfy required properties, whatever they are)   3. I do not use the Metropolis algorithm correctly   4. My distribution tests are not correct, since I do not have independent sample.


title: Yates' continuity correction in confidence interval returned by prop.test

text: This is the confidence interval estimated by prop.test               n <- 600; x <- 276; p <- 0.40     prop.test(x, n, p, alternative="two.sided", conf.level=0.95, correct=T)     95 percent confidence interval:      0.4196787 0.5008409       I tried to reproduce it, reading the code under prop.test. Here is a simplified way to get those two limits               ESTIMATE <- x/n     YATES <- 0.5     conf.level <- 0.95     z <- qnorm((1 + conf.level)/2)     YATES <- min(YATES, abs(x - n * p))      z22n <- z^2/(2 * n)     p.c <- ESTIMATE + YATES/n     (p.c + z22n + z * sqrt(p.c * (1 - p.c)/n + z22n/(2 * n)))/(1 + 2 * z22n)     [1] 0.5008409     p.c <- ESTIMATE - YATES/n     (p.c + z22n - z * sqrt(p.c * (1 - p.c)/n + z22n/(2 * n)))/(1 + 2 * z22n)     [1] 0.4196787      Can you explain to me why the underlying probability of success (p) is used in line 5? or maybe you could suggest where can I find more info about this YATES correction that affects the ESTIMATE. Thank you


title: How to run chi-squared test on imputed data

text: I have a survey data set with missing values and I generated 10 multiple imputations in which the missing values were imputed. There are several categorical variables in the data sets and I'd like to see if there's any association among these variables using chi-squared test. However, I've been searching how to run this analysis in either SAS or Stata but didn't find a solution. It seems that mi estimate doesn't support svy: tab in Stata. I can run chi-squared test on each imputed data set using in SAS using proc surveyfreq, however, I don't know how to combine the results. Can anyone please guide me through the steps on how to test the associations between categorical variables? Thank you in advance! Jin


title: Should I test for all the OLS assumptions for a pooled OLS from panel data?

text: As the title says, I need to perform a Pooled OLS, a Fixed effects and a Random effects analysis. In the case of a normal OLS, one should test for normality, collinearity, homoscedasticity, linearity, etc. I have been following the steps described here, but I am not so sure if I should do that in my case. My model is like the one described here. ![enter image description here](http://i.stack.imgur.com/6Loo7.png) Should I follow the same tests as in a normal OLS regression? If so, should I include the time dummy in the tests?


title: Guidelines on assessing how data is gathered

text: I will be working with a dataset sample which I will use for optimization purposes. However, before I start, I want to briefly analyze data in order to get a better sense and identify possible directions to explore. My problem is I don't know how the data was sampled, which sampling method was used or whether the same method was used for all parts of the data, and I'd like to get a sense about those. One of my ideas was to form subsamples from the sample and perform hypothesis tests but I don't know if this should be the way to go. I am aware that this is a vague question, but I was wondering if you could point out references/lectures/papers or comments in order to help me have a good start.


title: Is a heat-map of gene expression more informative if Z-scores are used instead of actual expression measurement values?

text: I have a heat-map of gene expression measurements (log2-transformed microarray signals, after inter-microarray data normalization, etc.) that I am using to illustrate the expression of 72 genes ('rows' of the heat-map) which I had identified as differentially expressed among different sub-groups of the 60 samples ('columns' of the heat-map, ordered by sub-groups) of my study. The ranges of gene measurements are within the 1-12 range (e.g., 4-8 for gene X, 2-10 for gene Y, and so on). It is a two-color heat-map, with the brightest green, black, and brightest red colors of the color scale used for values 1, 4 and 12, respectively. A reviewer has commented that the heat-map will be more informative if Z-scores of the gene expression measurements are used instead. I don't get this because to me it seems that the heat-map will be less informative; Z-scoring will reduce the dimensionality of the data as one can no longer compare one gene to another for a given sample. Can anyone comment on this? Thanks. An image showing the current and reviewer-proposed heat-maps can be seen here: http://i.imgur.com/a2hmT.png


title: Appropriateness of calculating scale means based on available non-missing responses (i.e., person-mean imputation)

text: This question came up in a consulting context, and I was interested in your thoughts. ### Context One strategy for dealing with occasional missing data when calculating scale means looks like this in the language of SPSS:               COMPUTE depmean =  mean.4(dep1, dep2, dep3, dep4, dep5, dep6).     EXECUTE.      I.e., calculate the mean of a psychological scale such as depression by taking the mean of six items. If a participant has four or more non-missing items, return the mean of the non-missing items. If the participant has three or fewer non-missing items, return missing. Of course the number of items in the scale and the threshold number items for calculating the mean can vary. ### Question   * In general, under what conditions, would you see this method of dealing with missing data to be appropriate?   * If you perceive it to be inappropriate, what alternative procedure would you recommend?


title: R: statistical test

text: I am working on 4 different species of tomatoes. From the data I had, I looked at the occurrence of a particular "event" in certain intervals of their genome (this interval is identical in all 4 plants) and I have a file for each of the species with their probability of occurrence. The file looks something like this: > ch01:a1-b1   >  p = 0.45   >  ch01:a2-b2   >  p = 0.005   >  ...   >  ... I have 4 such files (one for each species) where the chromosomal locations/intervals are identical and the probability value differs. Now I would like to find out those intervals where the occurrence of this event stands out. What do I mean by that? If I find, for example, over a particular interval, the event is more significant in plant/file 1 and 3, but not the others, then this interval is of interest to me. I would like to find all such intervals where "the occurrence of event is not all `equally` significant/insignificant". My question is how should I go about this? Which parametric/non-parametric test can I use? Can I use the test on each interval separately? What are the usual issues to look for and overcome.


title: How do I remove all but one specific duplicate record in an R data frame?

text: I have a data frame that contains some duplicate ids. I want to remove records with duplicate ids, keeping only the row with the maximum value. So for structured like this (other variables not shown):               id var_1     1 2     1 4     2 1     2 3     3 5     4 2      I want to generate this:               id var_1     1 4     2 3     3 5     4 2      I know about unique() and duplicated(), but I can't figure out how to incorporate the maximization rule...


title: Correct test for proportions

text: I have N objects which can have either 1 or both properties, A & B. I want to know if having property A raises the chance of having property B. My data is: ![data](http://i.imgur.com/4ILrUtG.png) Where S means they have the property, and N means they dont (ie, 263 objects have A & B, while 188 have B but not A). What is the correct test to see whether having property A raises the chance of having property B? I guess I can just use a 2 proportion test , but I feel like theres somethnig more specialized to this problem... Thanks!


title: Bootstrapping with a small number of observations

text: Let's say I've collected a _small_ number (N) of observations about a hypothesis that I'd like to test. I could use the bootstrap method to produce a sample distribution for the mean result of N observations, but I'm concerned that this model could break down when N gets very small, introducing error into the sample distribution itself. So my question is, how can I determine what the minimum N is that I need for reasonable results; or more quantitatively, how is N tied to the sampling error as N->0? _Update:_ I am coming to understand that the minimum value for N will vary based on the nature of the underlying data. So, in this case what meta- observations can I make to help me determine this? I don't know the true underlying distribution, or else I wouldn't need the bootstrapping...


title: interpret Alleffects() from effects package

text: I have this example logit model where some of the variables are factors but I'm not too sure how to interpret the effects. If I understand logit models correctly the coefficients that we get from the fitted model are the change in log-odds per unit change in the explanatory variable holding everything else constant. If take the `exp()` of the coefficients then I have the odds. I'm interested in the impact over the probability of some cathegorical variable to be " _male_ " or " _female_ " for instance. Could you please help me understand this? if `allEffects()` is not what I'm looking for, could you please let me know how could I get them?               library(effects)          titanic <- glm(survived ~ passengerClass + sex + age,data=Titanic,                    family=binomial)     titanic.all <- allEffects(titanic, typical=median,                                given.values=c(passengerClass2nd=1/3,                               passengerClass3rd=1/3, sexmale=0.5))     plot(titanic.all,          ticks=list(at=c(.01, .05, seq(.1, .9, by=.2), .95, .99)),          ask=FALSE)      EDIT: I don't think it is a duplicate. I'm interested in the output of the package `effects`, in particular in the output of the function `allEffects`. I found one document with the following notes "Notice that the print method for the object returned by allEffects reports tables of the effects, which, by default, are on the scale of the response variable, for a logit model, on the probability scale"


title: Problem on interpretation of data

text: Following Critical Phenomena in Natural Sciences of Didier Sornette, I am plotting the maximum value among $N$ variables. Doing this I can see the tail behaviour of a distribution (page 19 of C.P.N.S. of Sornette). Calling $\lambda$ the maximum value among the first $N$ variables, in the case the tail is of the form $f(x)\sim e^{-x/a}$ then $\lambda\sim a \ln N$. In the case the tail of the distribution is of the form of power law $f(x)\sim x^{-(a+1)}$ then $\lambda\sim N^{1/a}$. Now I have many samples of data coming from the same physical phenomenon, but they show two different behaviour. Ones with tail of power law type and ones with exponential tail type. Does this means that are presents two different mechanisms generating the same physical phenomenon?


title: Compare different statistical models that forecast the same result

text: **Background:** I'm assessing the future condition of a product where we have 18,000 units in the inventory. The product can be assessed as Poor, Fair, Good, or Excellent, based on a `condition_index` rating that we determine. For example, if for one product, the index reading is below 20, it's considered poor, 20 to 50 is fair, 50 to 70 is good, and 70+ is considered excellent It takes eight years to go through a full cycle of condition data for all 18,000 units. With these tests, we take each product in the inventory and provide a `condition_index` manually with a reader. On the other hand, we can use deterioration curves that say, if the product is x years old then it is in y condition. The way we're reporting at the end of the year is the following way:               Condition   Frequency     Poor        10% (1800)     Fair        20% (3600)     Good        30% (5400)     Excellent   40% (7200)      Here are how we perform the models: **Model 1:**   * Starting at product with `ID = 1`, we measure the `condition_index` directly from the product.    * It takes 8 years to cycle through all the products as we can only test about 2,500 a year.    * At the end of the year, we average all the ratings from all years up according to the table above, and determine the `condition_index` that way. **Model 2:**   * We use the yield curves based on the age of current products. The yield curves go from 100% to a 0 where 100% represents a product with no damage and 0 represents a product that is fully damage and unusable.   * For example, if we replace, rehab, or replace `product_id = 5` with a new one, then its age goes to 100%.    * Based on the cycle of changing and maintaining the products, we determine the condition index that we have. so based on the age of a product, we determine a condition_index and then report on the whole inventory that way. Our previous modelers have used Model 2 as our "official" condition reporting method; citing that consistency of data is not reliable in Model 1 as well as model 2 provides for a consistent replicable methodology. I'm not sure that I fully agree as I think parts fo the data from Model 1 should be used since it is actual collected product data and even better if we can combine it with model 2. **What are some of the ways I can compare method 1 to method 2 to determine which is the proper one to use?** **EDIT 1** Note that the current results show that model 1 is predicting that we have 20% in poor vs 10% in poor when using model 2


title: Symmetrical regression

text: Ordinary least squares regression is clearly not symmetrical! Regressing $Y$ on $X$ is not the same as regressing $X$ on $Y$. What kind of regression is symmetrical? I have 5 variables ($A$, $B$, $C$, $D$, $E$) and I want to find out the true relationship among the variables. $$ W_a * A + W_b * B + W_c * C + W_d * D + W_e * E + W_0 = 0 $$ where ($W_a$, $W_b$, $W_c$, $W_d$, $W_e$ and $W_0$ are constants). Is Deming regression symmetrical?


title: Mixed model in R with paired design and covariates

text: I make habitat suitability models for animal species. The purpose of my research is to investigate the accuracy of different models. I clearly have a nested design:   * accuracy_measure -> response variable   * 2 model types (model_type) -> fixed effect   * 230 species (species) -> random effect   * 10 replicates/species (replicate) -> random effect   * 10 subreplicate/replicate -> observation So I have: $10\times 10\times 230$ observations/model, identified as `speciesID`-`replicate`-`subreplicate` (`species_ID` ranging 1:230, `replicate` 1:10 and `subreplicate` 1:10) One could think about such mixed-effect model:               my.model <- lme(fixed = accuracy_measure~model_type, random = ~1|species/replicate)      I do not expand here into model simplification, yet... but here are my questions:   1. my replicates & subreplicates are paired in the sense that they come from the same split of the data. As an example for species X, the 20 observations of `replicate1` of model A and B (10 for A and 10 for B) are linked by a same data split which is likely to influence my accuracy measure. In the same way the 2 observations `replicate1`-`subreplicate2` of model A and B (1 for A and 1 for B) are also linked. Is there a way to introduce such pairing in a mixed effect model?   2. several continuous covariates, attributes of species (number of points for modelization, size in km2 of the range) may influence the measures of accuracy and I may be interested in investigating those effects. How could I include a covariate in such model? How should I strucutre it given that there are species covariates (high order of nesting)? I hope that my questions are relevant!


title: Calculate tail probabilities from density() call in R

text: This question concerns how to implement the following problem in R.               x = rnorm(1000)     hist(x,freq=FALSE)     lines(density(x))      How would you calculate the upper (or lower) tail probability for a given cutoff (e.g. +1) given the density estimate above? NOTE: the following solution isn't good enough. I need the calculation based on the smoothed density curve, i.e. an integral of the curve, not the empirical histogram.               sum(x>1)/length(x)      Also please do not suggest the use any of the standard `pnorm` functions because they are only correct if the underlying distribution is correctly specified. Thanks!


title: Doubt in derivative of logarithm

text: This seem to silly but I wanted to confirm if the derivative of the log- likelihood $\hskip 2 pt l(x_i)$. The derivative of $$\frac{d (\sum_{i=1}^{M} log(x_i))}{dx} = \frac{1}{x_i} \sum_{i=1}^{M} \frac{1}{x_i}$$ Is this correct? UPDATE of the Question : the pdf is $$p(x) = dP(x)/dx = d*x^{d-1}$$ where $P(x) = x^d$ $0<x<1$ Then the log-likelihood of x is $$ln L(x) = ln[p(x_1).p(x_2).\ldots.p(x_M)] = M*ln(d) + (d-1)\sum_{i=1}^{M} log(x_i)$$ Taking derivative of $ln(x)$ w.r.t to x or $x_i$ so as to find an estimate of $x$. The purpose is to find the minimum x's and the minimum value where found would correspond to the parameter. How do I proceed?


title: How do you use in sample error estimator in regression?

text: In their book The Elements of Statistical Learning, Friedman and coll. talk about the in sample error $Err_{in}$ (p. 229): $$Err_{in} = \frac{1}{n}\sum_{i=1}^n\mathbb{E}\left[\left(Y_i-\hat{f}\left( x_i\right)\right)^2\right] $$ there exists an unbiased estimator of this error in case of linear regression which is $$Err_{in}=\text{training_error} + 2\frac{d}{n}\hat{\sigma_{\epsilon}}^2$$ with $d$ the dimension, $n$ the number of observations and $\hat{\sigma_{\epsilon}}^2$ the unbiased estimator of the residuals. My question is how do you use this criterion ?


title: Cloud computing platforms for machine learning

text: I've got a small list of companies that provide a platform for running R, python, or octave scripts on clusters built on top of amazon EC2. Are there other names I should add?   1. Cloudnumbers   2. Opani   3. crdata


title: Heatmap for co-localization of factors in a set of clusters?

text: I want to represent the co-occurrence of factor groups within clusters in a heat map that reflects the colocalization frequency of each pair of factors in clusters (yellow for more frequently colocalized, red means less). After trying different things, I came up with the code below. Is this a sensible way of representing this data set?               set.seed(1)     x = c(       paste("cluster-",sample(c(1:30000),800000, replace=TRUE),sep=""),       paste("cluster-",sample(c(1:300),100000, replace=TRUE),sep=""),       paste("cluster-",sample(c(600:900),100000, replace=TRUE),sep="")          )     y = c(       paste("factor-",sample(c(letters[1:19]),800000, replace=TRUE),sep=""),       paste("factor-",sample(c(letters[1:3]),100000, replace=TRUE),sep=""),       paste("factor-",sample(c(letters[6:9]),100000, replace=TRUE),sep="")          )     d = data.frame(x,y)     t = table(d)     dat = rbind(t[,])     dats = dat/apply(dat,1,sum)     cdats = cor(dats)     heatmap(cdats)      PS: also, I would like to add a correlation score legend of the bottom right corner that shows the colour code for correlation from 1 to 0. How can I do that?


title: How to test for a break in a time series cycle

text: I've been scratching my head over this issue and would appreciate some help. I have a time series from 1920-2011 which I've used a Baxter Kings filter on to detrend. I would like to test whether the series has fundamentally changed between periods 1950-1979 and 1980 to 2009. Spectral analysis shows periodicity is increasing in frequency, but I would like to go at this in more ways than just that. Spec analysis also says that coherence is unity at each frequency, which seems suspect to me. If anyone has any ideas, please throw them out there, I would be so appreciative!


title: What is the latent class in pLSA

text: I want to implement Probabilistic Latent Semantic Analysis(pLSA) in Python. I have searched many times but couldn't find a simple tutorial.               terms     d1  d2  d3  d4     t1        1   1   0   1     t2        1   1   0   1     t3        1   1   1   1     t4        0   1   0   1      First we need to get P ( w , d ) = ∑ c P ( c ) P ( d | c ) P ( w | c ) = P ( d ) ∑ c P ( c | d ) P ( w | c ) What is this latent class(c) used in the algorithm? How to find it? Please provide a simple explanation using the example dataset?


title: Random number generation distributed like a translated weibull from uniform random generator

text: If $X$ is uniformly distributed on $(0,1)$, then the random variable $ \lambda(-\ln(1-X))^{1/k}\ $, is Weibull distributed with parameters $k$ and $\lambda$. With this, I can get random numbers distributed weibull from a uniform random number generator. but if I have a translated Weibull distribution $ f(x, k, \lambda, \theta) = \frac{k}{\lambda} (\frac{x-\theta}{\lambda})^{k-1} \exp^{-(\frac{x-\theta}{\lambda})^k} $ how is the transformation to generate random numbers distributed like a translated weibull from a generator of uniform random numbers ? ### EDIT the question it's not about the density function, it about of the transformation. How I can generate random numbers distributed like a translated weibull ? (using a uniform random number generator)


title: Reconstruction error of PCA when space dimensionality is larger than number of points

text: I've got a question and I have done several experiments in R, yet couldn't figure out why. The question is for a data set of N*D, N for number of data points and D for dimension, the maximum number of principal component is max(N, D). Then I discovered a fact that, when N is larger than D, by using D principal components, any old data(data within the data set) could be reconstructed absolutely the same as it was, so can any new data, the result of `((new_data - mean) %*% loadings %*% t(loadings) + mean) ### code in R` is the same as new_data. However, when D is larger than N, by using N principal components, reconstruction of data within the data set is still precise, that of the new data differs from the originals. I just can't figure out why. Can you please explain or give some hint?


title: Autocorrelation of nonstationary IMA(1,1) process

text: I knwo that the autocorrelation in MA(1) process varies between -.5 and +.5, if we consider d(t)=c+e(t)−θ⋅e(t−1),then for positive values of Theta, autocorrelation is negative and for negative values of Theta autocorrelation is positive. Now I'm wondering if the autocorrelation of IMA(1,1) process also follows the same pattern as MA(1)? If the process is IMA(1,1) as follows d(t)=c+d(t-1)+e(t)−θ⋅e(t−1), whether for negative values autocorrelation is positive and for positive values autocorrelation is negatve? how can I calculate the range of autocorrelation for IMA(1,1) process?


title: Calculating probabilities of an nth step transition matrix for discrete time markov chains

text: "Let $\\{X_n, n \geq 0\\}$ be a DTMC with state space $S = \\{1, 2, 3, 4, 5\\}$ and the following transition probability matrix: $$ P = \begin{pmatrix} 0.1 & 0.0 & 0.2 & 0.3 & 0.4 \\\ 0.0 & 0.6 & 0.0 & 0.4 & 0.0 \\\ 0.2 & 0.0 & 0.0 & 0.4 & 0.4 \\\ 0.0 & 0.4 & 0.0 & 0.5 & 0.1 \\\ 0.6 & 0.0 & 0.3 & 0.1 & 0.0 \end{pmatrix} $$ with the initial distribution as: $$ a = \begin{pmatrix} 0.5 & 0 & 0 & 0 & 0.5 \end{pmatrix}. $$ Compute $P(X_2 = 2, X_4 = 5)$" Also, how do I do this? I'm not even sure what it means. Is this asking me to work out the probability of $X_4 = 5$ given that I have $X_2 = 2$? EDIT: The matrix is correct. From the previous part of the question (that I didn't post), $P(X_2 = 2) = 0.08$. In the answers, my lecturer uses the Chapman Kolmogorov equations. Does this make the answer any clearer? EDIT 2: Edited to bump to the top as I have commented on whuber's post with a couple of questions


title: What does #iter mean in libsvm when using RBF as a kernel?

text: I have a dataset with 60,000 training examples. I tried the following command under windows and with libsvm:               svm-train.exe -t 2 -g 0.07 -c 1 images.train      and it started giving multiple outputs for example:               .....*..*     optimization finished, #iter = 7090     nu = 0.117216     obj = -667.069565, rho = -0.055205     nSV = 4702, nBSV = 26     ....*..*     optimization finished, #iter = 6485     nu = 0.100224     obj = -565.615150, rho = -0.198917     nSV = 4342, nBSV = 29     ..*.*     optimization finished, #iter = 3671     nu = 0.027754     obj = -176.698407, rho = -0.838783     nSV = 2560, nBSV = 52     ....*..*     optimization finished, #iter = 6192     nu = 0.104976     obj = -604.912230, rho = -0.284038     nSV = 4109, nBSV = 96     ......^C      I stopped it because I didn't understand what was going on. I searched online in the libsvm website however I did not find any explanation on why I get multiple outputs and what these outputs mean exactly. I tried running the heart_scale and I only got one output:               svm-train.exe -t 2 -g 0.07 -c 1 heart_scale     *     optimization finished, #iter = 134     nu = 0.433785     obj = -101.855060, rho = 0.426412     nSV = 130, nBSV = 107     Total nSV = 130      so does the amount of training samples have anything to do with this and if so what's the exact correlation? thanks Thank you in advance


title: How to assess heterogeneity between fitted, correlated parameter values?

text: I have four datasets, and for each one, I fit the same 3-parameter model. I'm interested in making a statement about the heterogeneity of each parameter value of the model, across the four datasets (and also comparing it to parameters from a fit to the pooled dataset). My first thought was to use Cochrane's $Q$ and $I^2$ to assess heterogeneity for each parameter, but the three parameters aren't independent - they are intrinsically linked in the model itself (essentially, the model is a probit function (2 parameters) using a generalized mean (1 parameter) as input, so the generalized mean parameter is going to be highly correlated with one of the probit parameters). In some sense, I'm trying to compare the results of a traditional meta- analysis (combining the three parameter values from the 4 dataset fits) with a pooled analysis approach (the model parameter values fit from the combined dataset). Any insight on what I'm missing would be much appreciated!


title: expectation of an exponential function

text: What is the expectation of an exponential function: $$\mathbb{E}[\exp(A x)] = \exp((1/2) A^2)\,?$$ I am struggling to find references that shows this, can anyone help me please? I am assuming Gaussian distribution. A is a constant and x is a random variable that is gaussian distributed.


title: Number of significant figures in a chi square fit

text: I'm working on a scientific paper in which goodness-of-fit and a parameter are estimated by a $\chi^2$ minimization, and I'd like to know whether there is a general rule for how many significant figures one should quote for the value of the $\chi^2$ and its resulting P-value. In particular I'm fitting one parameter from data in 31 bins, each with its contents in the gaussian regime. I'd guess the precision of my minimization algorithm is not the limiting factor (my code spits out $\chi^2$ = 38.8898 / 30, so the P-value is something like 0.1282) for signifcant figures, but what is? Is it the number of degrees of freedom? The degree to which the gaussian assumption is invalid in the least-populated bin(s)? Most papers I see typically quote to tenths in $\chi^2$ (so 38.9 / 30 in my case) and hundredths in P-value (0.13 for me), but I'm not sure if that's always reasonable. P.S.- I'd like to apologize if any terminology I used is funny: in physics we use annoyingly non-standard terms when discussing statistics.


title: Denoising of microarray data using PCA

text: I'm looking for references and comments regarding the validaty of the following method for data denoising, which I found while reading a code doing analysis of some gene expression dataset. The dataset consists of columns $x_1, ..., x_n$ of length $m$ (expression levels for $n$ genes observed in m samples). Someone with the knowledge of the dataset said that if we look at the top 10% columns with the maximum variance, we find that those columns have the maximum variance due to some artifacts (noise) in the measurement of expression level of the corresponding genes. In addition, we know that the rest 90% of columns are either not affected by any kind of noise, or are affected by some noise but to a much lesser degree than the top 10% columns. Now, in the code that I'm examining the following method is used to remove the noise from the dataset. They calculate the principal components $y_1, ..., y_n$ for variables $x_1, ..., x_n$. They took $y_1$ (the leading principal components) and (this is my guess) assume that it mostly captures the variance caused by the artifacts described above. Then they transform the data (all $n$ columns) using the following rule: $x_i = x_i - (\text{projection of } x_i \text{ onto } y_1)$. That is, from each column they remove the component that is collinear with $y_1$, and keep the component that is orthogonal to $y_1$. Can anybody please provide any references for this method or comment on its applicability in this case?


title: Using an aggregated outcome variable as contextual in multilevel analysis

text: A main hypothesis in our study is that people’s decisions (our dichotomous outcome variable) are influenced by how people in their social context have decided in the same matter. In a multilevel analysis perspective, is it feasible to use as a level 2, or contextual variable, the aggregated (mean) of the outcome variable? KML


title: Top-Down or Bottom-Up Approach for demand forecasting

text: I have 5000 SKUs which all of them are highly positive autocorrelated, to get the item level forecast for all5000 SKUs (disaggregate forecast) which approach can provide more accurate forecasts, BU or TD, SES is forecasting method? and why? BU approach: we do the forecast for all 500 SKUs directly TD approach: first we sum up the demand for all SKUs, then we do the forecast for aggregate demand, finally according to the weight of each SKU, we will obtain the disaggregate forecast for each SKU!(f(i,t)=pi*Ft) Thank you in advance!


title: Is it possible to compute integrals in JAGS?

text: Is it possible to specify a model that uses an integral in the JAGS dialect of BUGS? (I believe this is possible in *BUGS, but I use a Mac, on which it seems JAGS is better supported.)


title: What is the standard deviation of the sum of THREE correlated random variables?

text: Basically I want to know how to extend the following formula to the three variable case: $$\mbox{var}(aX+bY) = a^2\mbox{var}(X)+ b^2\mbox{var}(Y) + 2ab \sqrt{ \mbox{var}(X) \mbox{var}(Y)} \mbox{corr}(X,Y)$$ How can I calculate $\mbox{var}(aX + bY + cZ) $? The answer to the two variable case is explained here. It's possible that the general case is covered too, but I didn't quite follow the matrix algebra and wondered if there is an explicit solution for the three variable case.


title: What advantages do "internally studentized residuals" offer over raw estimated residuals in terms of diagnosing potential influential datapoints?

text: The reason I ask this is because it seems that internally studentized residuals seem to have the same pattern as raw estimated residuals. It would be great if someone could offer an explanation.


title: Best test to use

text: Suppose a group of people are exposed to some chemical for $20$ minutes. Their heart rate is measured before being exposed. The exposure to the chemical seems to cause their heart rate to speed up. What test can be used to precisely find the point where the heart rate differs significantly from baseline? A paired t-test?


title: Two highly correlated variables where both correlate with a third: Correlation and Causation

text: My question is certainly quite basic for statisticians! Let's suppose `Var1` and `Var2` are highly correlated with a poor $R^2$. `Var3` is another variable that we will use in a regression (standard linear model (Gaussian error distribution, OLS estimator)) as a response variable. `Var1` and `Var2` are the explanatory variables of this regression where the interaction effect is not computed. Let's call this regression `regression1` Is it possible that `Var1` affects `Var3` (highly significantly) but not `Var2` (not significant at all)? Let's assume that the causation is the following: `Var1` causes `Var2` and `Var3` Is it plausible that both `Var2` and `Var1` are associated with highly significant p-value in `regression1`? * * * For example: In a given species of fish, there is a very significant relation between the size of the fishes and the color of the fishes. But the linear model that gives this very significant p-value does not explain much of the total variance in fish color and size. When performing a regression of depth on fish size and fish color (without the interaction term), I get two highly significant p-values (for the two explanatory variables). Can I infer that both the size of the fishes and the color of the fishes "influence separately" the depth at which the fishes are? Or maybe depth has an effect on the color of the fishes (because the diet is different for example) and the color of the fishes influenced the size of the fishes (because fishes that are bright need be big to escape predators or something like that). In such case, the size of the fishes would be associated to a significant p-value (in the regression of depth) only because it is correlated to the other explanatory variable which is the color of the fishes.


title: What are pure random sampling and orthogonal sampling?

text: I'm working on a project to calculate the Area of Mandelbrot Set by Monte Carlo method. I implemented my code by generating standard random numbers and throwing them into the area. I used simple sampling, and I also implemented Latin Hypercube Sampling. I should also implement my project with _Pure Sampling_ and _Orthogonal Sampling_ method. I read Wikipedia but it is not clear for me what is exactly Orthogonal Sampling and how I should implement it. I implement Latin Hypercube Sampling by dividing the area and generating random arrays for each sub-area and counting. What about Orthogonal? That would be great if someone help me to understand the concept at least.


title: Book on data modeling

text: I'm looking for reference to get a better understanding and overview of the methods for data modeling (mostly related to business questions like finding categories in data and setting of scoring functions for forecasting). After some internet research I couldn't single out a canonical reference. I have some very basic knowledge of statistics and I believe the topics that could be most useful for me would be SVM, Decision Trees, Clustering, Covariance Matrix, Correlations, Multivariate analysis ... just to name some keywords for the direction. Also I'd like to have some foundation to eventually dive into other statistical and data mining topics. While the reference should introduce the concepts, I'd like it to be fast- paced (not chatty, not too many examples) and mathematically thorough while concentrating on applications. Of course there are specialized books on each topic, but I'd rather learn the basics of all of these topic quickly. So ideally one or two books for all of them together. I cannot think of a different way to explain :) Can you suggest good books? Or maybe just help me get some orientation in the vast field of statistics so that I know what to look for.


title: Is the Wilcoxon paired rank test appropriate for this data?

text: I have two vectors of normalized data values representing two paired conditions ("methylation" and "expression"). A scatterplot of my data looks like this: ![enter image description here](http://i.stack.imgur.com/oq9dK.png) I'd like to know if I am using the Wilcoxon signed-rank test correctly in order to answer some questions. Specifically, I would like to identify those data points whose difference differs significantly from the median. The p-value that `wilcox.test()` reports in R suggests that the alternative hypothesis is true (that the median of the difference of the two conditions is not zero, or that H0 fails):               > wilcox.test(mtx$methylation_difference_normalized, mtx$expression_difference_normalized, alternative = "two.sided", paired = TRUE, exact = FALSE, correct = TRUE, conf.int = TRUE, conf.level = 0.95)     V = 57049625, p-value < 2.2e-16     alternative hypothesis: true location shift is not equal to 0     95 percent confidence interval:      0.2722975 0.3382661     sample estimates:     (pseudo)median      0.3250247      Assuming this is the case, is the confidence interval that is reported usable for identifying points in the scatterplot which are significant? Given this result, for example, can I simply take the (absolute?) difference of any pair of observations, and if that difference falls outside this interval, do I identify that data point as statistically significant? Failing that, is there another test I should be using instead to accomplish the same task?


title: How Gradient boosting can be more interpretable than CART?

text: I found this document which compare some learning methods and I don't understand this table : ![enter image description here](http://i.stack.imgur.com/UEat0.png) Gradient boosting has a better intepratability score than CART. How is it possible ? I thought gradient boosting was an ensemble method, so we have to take a look at all the trees in order to understand how the classification is built. Dou you think there is a mistake here ?


title: Do I discuss the lack of interaction if I already discussed lack of main effects?

text: I conducted an experiment with 2 independent group variables (gender and language - language indicates either native speakers or 2nd language speakers) and ran a 2 way ANOVA. There were no main effects, and no interaction. My hypotheses had been that there would be main effects and interaction. So, I rejected all 3 hypotheses. I discussed what the reason may have been for having no main effects. Do I separately need to discuss the reasons for no interaction? Also, even through the interaction did not reach significance, there was a 'trend' for males to perform better in a second language, while for females to perform better in their first language. Do I need to comment on that if I already mentioned it in results section?


title: Perform a Kruskal-Wallis test in Stata

text: Could anyone help me perform a Kruskal-Wallis test in Stata? I am not familiar with Kruskal-Wallis test, but I have to perform a nonparametric test, since the normality hypothesis is rejected, to confirm the results of a parametric test already performed (F-test). I have got the data already ranked in groups though I am not understanding the results. chi-squared = 505.166 with 179 d.f. probability = 0.0001 chi-squared with ties = 505.195 with 179 d.f. probability = 0.0001


title: Spread-Level Plot versus Power Transformation Functions in R

text: I'm having trouble interpreting the results from the Spread-Level Plot function in R (car package). The documentation says: > PowerTransformation   >  spread-stabilizing power transformation, calculated as 1 - slope of the > line fit to the plot. This is not explicit enough for me. Should this transformation be applied to every variable in the regression? For example, assume I have an lm object given by:               myFit <- lm(y ~ x1 + x2)      Then I use Spread-Level Plot:               slp(myFit)      If the 'suggested power transformation' is 0.5, then does that imply a homoscedastic model could be fit using one of the following?               refitA <- lm(sqrt(y) ~ sqrt(x1) + sqrt(x2))     refitB <- lm(sqrt(y) ~ x1 + x2)     refitC <- lm(sqrt(y) ~ sqrt(x1 + x2))      If I understand correct, refitA would be the suggested model to approximate homoscedasticity. On the other hand, if I _only_ want to transform the LHS, I would use the `powerTransform` function (also from car package). i.e., an "estimated transform parameter" of 0.5 from the powerTransform function would imply that refitB is homoscedastic. Is this correct? Thanks!


title: Problem with Finding Likelihood: Bayesian

text: I am really unfamiliar with Bayesian methods particularly parameter estimation. Suppose I have a test to find a parameter, theta which is the number of packaged bag for retail sale that could contain a hidden toy. I consider the following values of theta to be possible: `{0, 1, 2, 3, 4, 5}`. I am interested in the following two hypotheses:               H0: theta == (no toys)      H1: theta =/= 0 (there are toys in a bag)      To test these hypotheses, I theoretically perform an experiment. Instead of testing all five of the bags (avoid damaging the packaging), I pick one of the bags at random and test for a toy. Using a testing machine, it is designed to say 1 if there are toys found and 0 if no toys were found. When a bag was tested, the result (the data) was:               x1 = 0       That is, no toys were found in the first bag that was tested. I've used parameter estimation theory using prior distributions. But I am struggling with finding the likelihood. Am I supposed to do a Beta distribution to do so? Or binomial, in which there is 0 successes in 1 trial for theta?


title: What kind of stochastic processes can Ljung–Box_test apply to or assume?

text: Ljung-Box test tests zero autocorrelation in a time series. I wonder what type of stochastic processes the test assumes on the time series? For example, whenever autocorrelation exists? By assumption, I mean the one that holds under both the null and alternative hypotheses. So the null will be the assumption plus zero autocorrelation. (Therefore the alternative will be the assumption plus nonzero autocorrelation.) Thanks! :)


title: Matlab classify function for multi-class classification

text: What kind of Machine Learning method does the function "classify" in Matlab use for the multi-class classification? Is it SVM? If so, how does it use the classifier for multiple classes?


title: How does Pearson's Chi Squared Statistic approximate a Chi Squared Distribution

text: So if Pearson's Chi Squared Statistic is given for a $1 \times N$ table, then its form is: $$\sum_{i=1}^n\frac{(O_i - E_i)^2}{E_i}$$ Then this approximates $\chi_{n-1}^2$, the Chi-Squared Distribution with $n-1$ degrees of freedom, as the sample size $N$ gets larger. What I don't understand is how this asymptotic approximation works. I feel like the $E_i$'s in the denominators should be replaced with $\frac{s_i^2}{n_i}$. Since that would give you $\chi_n^2 = \sum_{i=1}^nZ_i^2$, for $Z_i\sim n(0,1)$. But of course this has $n$ degrees of freedom, not $n-1$, so clearly something else is going on.


title: Simulating unconditional Gaussian markov random field

text: I am trying to sample the Gaussian markov random field or say multivariate gaussian distribution with some spatial correlation given by the precision matrix Q. Here is the algorithm that I am using               Compute the Cholesky factorization Q=LL'     Sample z ~ N(0,I)     Solve L'x = z      where x is the samples obtained. I have a grid of 20x20. And I have precision matrix such that any point in the spatial grid is related to its four neighbors. So if my Q matrix is of size 400x400, then the first row which is for the first point, the Q matrix is nonzero for columns 1,2 and 21 where 1 denotes the first element itself, 2 denotes the element right of it and column 21 denotes the element below it. So, I got the samples x from it and I plotted them in a grid of 20x20. However, I couldn't see much relation. May be I am wrongly interpreting the plot. So anyone with any suggestions? ![enter image description here](http://i.stack.imgur.com/Baoki.jpg) Further, I tried to increase the window size to 5x5Here is my code. However, when I plot the generated data, they don't seem to be spatially correlated or it's getting worse.               %Initialize the grid     row = 20;     column = 20;          N = row*column;     linearMatrix = zeros(row*column,3);     delta = 0.0001;          index=1;     for i=1:row         for j=1:column             linearMatrix(index,1) = i;             linearMatrix(index,2) = j;             linearMatrix(index,3) = index;             index = index+1;         end     end          sparseIndexX = zeros(N,1);     sparseIndexY = zeros(N,1);     sparseIndexValue = zeros(N,1);          index = 1;     for j=1:N         x = linearMatrix(j,1);         y = linearMatrix(j,2);         for k=1:N             x1 = linearMatrix(k,1);             y1 = linearMatrix(k,2);                  if(j == k)                 %diagonal element                 sparseIndexX(index,1) = j;                 sparseIndexY(index,1) = k;                           % For the four corner elements                 if j == 1 || j == N || (x == 1 && y == column) || (x == row && y == 1)                     sparseIndexValue(index,1) = 8;                 elseif (x == 1 && y == 2)  || (x == 1 && y == column-1) || (x == 2 && y == 1)  || (x == 2 && y == column) || + ...                    (x == row-1 && y == 1) ||  (x == row-1 && y == column) || (x == row && y == 2) || (x == row && y == column-1)                     % For the elements at the edge but not the corner                     sparseIndexValue(index,1) = 11;                 elseif (x==2 && y==2) || (x==2 && y==column-1) || (x==row-1 && y==2) || (x==row-1 && y==column-1)                     sparseIndexValue(index,1) = 15;                 elseif (x==2 && y>2) || (x==2 && y<column-1) || (x==row-1 && y>2) || (x==row-1 && y<column-1) || (x>2 && y==2) || (x<row-1 && y==2) || (x>2 && y==column-1) || (x<row-1 && y==column-1)                     sparseIndexValue(index,1) = 19;                 elseif x > 2 && y > 2 && x < row-1 && y < column-1                     % For the elements in the middle                     sparseIndexValue(index,1) = 24;                      else                     sparseIndexValue(index,1) = 14;                 end                      index = index + 1;             else                 %if the elements are adjacent                  if (abs(x-x1) <= 2 && abs(y-y1) == 0) || (abs(x-x1) == 0 && abs(y-y1) <= 2) || (abs(x-x1) <= 2 && abs(y-y1) <= 2)                     sparseIndexX(index,1) = j;                     sparseIndexY(index,1) = k;                     sparseIndexValue(index,1) = -1;                     index = index + 1;                 end             end         end     end          %Generate the precision matrix for spatial correlation     QSpatial = sparse(sparseIndexX', sparseIndexY', sparseIndexValue', N, N);          Q=QSpatial;          %Make the matrix strictly diagonally dominant so that it is positive definite and     %invertible     for i=1:N         Q(i,i) = Q(i,i) + delta;     end          v=zeros(N,1);          %Get the cholesky decomposition     %Step 1     L=chol(Q,'lower');          %Sample from the normal distribution     %Step 2     z=randn(1,N);     z=z';          %Back substitution method to generate the sample     %Step 3     for i=N:-1:1         total = 0;         if i~= N             for j=(i+1):N                 total=total + L(j,i) * v(j,1);             end         end         v(i,1) = 1/L(i,i) * (z(i,1) - total);     end          %Step 4     x=v;               x=reshape(x,20,20);     x=x';     imagesc(x);     colorbar;        This is the plot of the samples ![enter image description here](http://i.stack.imgur.com/BBB8k.jpg) Does anyone have any suggestion?


title: Why do linear regression and ANOVA give different $p$-value in case of considering interaction between variable?

text: I was trying to fit one time-series data (without replicates) using regression model. The data looks like follows:               > xx.2               value time treat         1  8.788269    1     0         2  7.964719    6     0         3  8.204051   12     0         4  9.041368   24     0         5  8.181555   48     0         6  8.041419   96     0         7  7.992336  144     0         8  7.948658    1     1         9  8.090211    6     1         10 8.031459   12     1         11 8.118308   24     1         12 7.699051   48     1         13 7.537120   96     1         14 7.268570  144     1      Because of lack of replicates, I treat the time as continuous variable. Column "treat" shows the case and control data, respectively. First, I fit the the model "value = time*treat" with "lm" in `R`:               summary(lm(value~time*treat,data=xx.2))          Call:     lm(formula = value ~ time * treat, data = xx.2)          Residuals:          Min       1Q   Median       3Q      Max      -0.50627 -0.12345  0.00296  0.04124  0.63785           Coefficients:                  Estimate Std. Error t value Pr(>|t|)         (Intercept)  8.493476   0.156345  54.325 1.08e-13 ***     time        -0.003748   0.002277  -1.646   0.1307         treat       -0.411271   0.221106  -1.860   0.0925 .       time:treat  -0.001938   0.003220  -0.602   0.5606          The pvalue of time and treat is not significant. While with anova, I got different results:                summary(aov(value~time*treat,data=xx.2))                 Df Sum Sq Mean Sq F value Pr(>F)       time         1 0.7726  0.7726   8.586 0.0150 *     treat        1 0.8852  0.8852   9.837 0.0106 *     time:treat   1 0.0326  0.0326   0.362 0.5606       Residuals   10 0.8998  0.0900                       The pvalue for time and treat changed. With linear regression, if I am right, it means the time and treat has no significant influence on value, but with ANOVA, it means time and treat has significant influence on value. Could someone explain to me why there is difference in these two methods, and which one to use?


title: Fitting Lognormal Distribution in WinBugs

text: Let's suppose that some data points follow $Lognormal(\mu,\sigma^2)$ and both parameters are unknown . My goal is to obtain the posterior distribution by assigning conjugate prior distributions on both $\mu$ and $\sigma^2$,. How this can be done in WinBugs?


title: How to explain what density is and the interpretation of the curve's height to non-statisticians

text: I tend to use histograms of continuous variables adding estimated density curves in order to compare several charts easily. However, I find difficulties when I try to explain what density is and the interpretation of the curve's height to non-statisticians. The stalwart Wikipedia provide us a good definition: > In probability theory, a probability density function (pdf), or density of a > continuous random variable, is a function that describes the relative > likelihood for this random variable to take on a given value. Do you think the concept "relative likelihood" would be clear enough for non- statisticians? What is your approach in such situations?


title: fitting the tail of a distribution in a regression tree

text: I have 3 integer valued time series $a_t$, $b_t$ and $y_t$ with $k$ observations. I want to fit $y_t$ with the 2 first, and for that purpose I use a regression tree like this:   * test all combinations of $a_t,\ldots a_{t-k}$, $b_t,\ldots b_{t-k}$, for each value between 0 and $N$. For each combination, take $\bar{y}_s$ the mean value on the subset as an estimator,   * select the splitting variable/value that minimize $\sum_{i\in s}(y_i - \bar{y}_s)^2$ I stop the recursion when the tree is deep enough or if I don't have enough point to compute $\bar{y}_s$. The distribution being very peaky and the dataset big, some combination of $a_t,b_t,a_{t-1},b_{t-1},\ldots$ are extremely frequent and some of them appear just a few times. While $\bar{y}_s$ seems to be a good estimator on the training set, it degrades very quickly on the tail of the distribution on the test set (but it is good on the test set for the center of the distribution). How to improve my estimator in my regression tree to better fit the tail of the distribution and keep a constant performance on the tail of the test set too ?


title: Convex cost functions and Fisher consistency

text: In Section 4.2 of the paper by Boucheron et al., the authors argue that the minimizer $f^*$ of the cost functional $$A(f) = \mathbb{E}\\{\phi(-Yf(X))\\}$$ is such that the classifier $g^*$, constructed from $f^*$ by $$ g^*(x) = \begin{cases} 1 & \text{if } f^*(x)\geq 0 \\\ -1 & \text{otherwise} \end{cases}$$ is the Bayes classifier. Note that the cost function $\phi$ is a non-negative, strictly convex, and non-decreasing function. In order to prove the claim, the authors argue that $\phi'$ is strictly increasing. Where is that coming from?


title: Negative binomial regression with $\beta>1$

text: In analyzing my non-independent count data with a negative binomial regression, I am finding that one of my variables (a proportion variable) consistently is showing a beta weight > 1. As a matter of fact it's over 3. I read that beta weights can be greater than one if there is suppression going on in a model with multiple predictors, but even when I put my proportion variable in the model by itself the value is still over three. I am currently running the analyses in SPSS and have used the same syntax I have always used for a negative binomial regression and this is the first time this has ever happened. Does anyone have any ideas about what may be going on?


title: R: truncated regression with upper and lower bounds

text: I'm looking for a package to do truncated regression in R. I found the truncreg and truncSP packages, but they seem to allow either a lower bound or an upper bound, but not both at the same time. All my response values are bound between 0 and 1, so I need to have both bounds taken into account. Any help? Thanks a lot!


title: Wilson Score CI for difference of two proportions

text: I am wondering somebody here may provide a formula or source for calculating Wilson score confidence interval of difference of two proportions (namely CI for p1-p2). It seems most available web sources are for one sample proportion. Thanks in advance


title: Specification of a random slopes model with two grouping factors

text: Suppose people take a 10-item exam. For each item, $k = 1,2,…,10$, exactly one rater assigns a score to exactly one person, with the constraint that no person sees the same rater twice. There are $i = 1,2,…,5000$ people and $j = 1,2,…,500$ raters. So people and raters are partially crossed, with most person-rater pairs having 0 data points, and the rest having exactly 1 data point. Update: items were randomly assigned to positions. Suppose I want to estimate a random slopes model where scores vary linearly with item position and the variability in scores is due to person intercept variability, person slope variability, rater intercept variability, rater slope variability, and residual error variability – in R, `lmer(rating ~ itemPosition + (itemPosition | personID) + (itemPosition | raterID))` Firstly, is there a better way to write the model down in non-matrix form than what I have below? Because there are two grouping factors, I am confused about the indices. If I let $q$ index data points, I want something like $Y_q = \beta_0 + \beta_1 ItemPosition_q + p_{0i} + p_{1i} ItemPosition_q + r_{0j} + r_{1j} ItemPosition_q + e_q$ where $p_{0i} \sim N(0, \sigma_{p0}^2)$ denotes a person-level random intercept, $r_{0j} \sim N(0, \sigma_{r0}^2)$ denotes a rater-level random intercept, $p_{1i} \sim N(0, \sigma_{p1}^2)$ denotes a person-level random slope, $r_{1i} \sim N(0, \sigma_{r1}^2)$ denotes a rater-level randoms slope, $e_q \sim N(0, \sigma^2)$ is the residual error, the person random effects are uncorrelated with the rater random effects, the random effects are uncorrelated with the error term, $cov(p_{0i}, p_{1i}) = \sigma_{p01}$, and $cov(r_{0i}, r_{1i}) = \sigma_{r01}$. More specifically, can I replace the index $q$ with $ijk$? Secondly, do I have the correct expression for $$Var(p_{0i} + p_{1i} ItemPosition_q + {r_0j} + r_{1j} ItemPosition_q + e_q) ?$$ I get $(\sigma_{0p}^2 + 2\sigma_{p01} ItemPosition_q + \sigma_{1p}^2 ItemPosition_q^2)$ + $(\sigma_{0r}^2 + 2\sigma_{r01} ItemPosition_q + \sigma_{1r}^2 ItemPosition_q^2)$ + $\sigma^2.$ I want to use this as the denominator to calculate the proportion of variance due to persons, for example. Lastly, any general comments about the appropriateness of the model are also welcome but not at all necessary (as I haven’t provided any information about the data).