title: The relationship between the gamma distribution and the normal distribution

text: I recently found it necessary to derive a pdf for the square of a normal random variable with mean 0. For whatever reason, I chose not to normalise the variance beforehand. If I did this correctly then this pdf is as follows: $$ N^2(x; \sigma^2) = \frac{1}{\sigma \sqrt{2 \pi} \sqrt{x}} e^{\frac{-x}{2\sigma^2}} $$ I noticed this was in fact just a parametrisation of a gamma distribution: $$ N^2(x; \sigma^2) = \operatorname{Gamma}(x; \frac{1}{2}, 2 \sigma^2) $$ And then, from the fact the sum of two gammas (with the same scale parameter) equals another gamma, it follows that the gamma is equivalent to the sum of $k$ squared normal random variables. $$ N^2_\Sigma(x; k, \sigma^2) = \operatorname{Gamma}(x; \frac{k}{2}, 2 \sigma^2) $$ This was a bit surprising to me. Even though I knew the $\chi^2$ distribution -- a distribution of the sum of squared _standard_ normal RVs -- was a special case of the gamma, I didn't realise the gamma was essentially just a generalisation allowing for the sum of normal random variables of _any_ variance. This also leads to other characterisations I had not come across before, such as the exponential distribution being equivalent to the sum of two squared normal distributions. This is all somewhat mysterious to me. Is the normal distribution fundamental to the derivation of the gamma distribution, in the manner I outlined above? Most resources I checked make no mention that the two distributions are intrinsically related like this, or even for that matter describe how the gamma is derived. This makes me think some lower-level truth is at play that I have simply highlighted in a convoluted way?


title: Comparing two opposite ranked sets of data

text: Looking for the best formula to compare two sets of numbers **Set #1 - Rank: Smaller numbers indicate higher "rank" (Ex: 1 is the best, 1,000,000 is the worst)** **Set #2 - Number of Items: The more items in a category, the worse off the category is** This formula is essentially for profitability, and is relating to Amazon Kindle Books Here's some background about how Amazon ranks things (set #1) For each category, there's roughly 10 best sellers that can be pulled, sometimes less if the category is smaller. For each best selling book, it's assigned a SalesRank (#1 is the best-selling book)   1. I counted how many best sellers there were for a category, then took each of their SalesRanks and averaged them (is there a better way?)   2. I then counted how many books were in each category. So, what I'm looking for is a formula (preferably spitting out on a 1-100 scale) that finds categories with high sales ranks, but a low number of books in the category Thanks in advance for any help from brilliant minds, and let me know if you need more info!


title: How to define silhouette for one cluster?

text: I want to compare two clustering algorithms. I took data that the first algorithm gathered in one cluster. The second algorithm gave 3 clusters for the same points. In order to compare the results, I tried to compute the two different silhouette for these algorithms, but I realized that in the first case, since there was only one cluster, there is no other cluster to compute the dissimilarities to. How can silhouette be defined in the case of one cluster? If it is not possible to define it, what would be a good way to compare the results of my two clustering methods?


title: Open source Java library for statistics at the level offered by a graduate statistics course

text: I am taking a graduate course in Applied Statistics that uses the following textbook (to give you a feel for the level of the material being covered): Statistical Concepts and Methods, by G. K. Bhattacharyya and R. A. Johnson. The Professor requires us to use SAS for the homeworks. My question is that: is there a Java library(ies), that can be used instead of SAS for problems typically seen in such classes. I am currently trying to make do with Apache Math Commons and though I am impressed with the library (it's ease of use and understandability) it seems to lack even simple things such as the ability to draw histograms (thinking of combining it with a charting library). I have looked at Colt, but my initial interest died down pretty quickly. Would appreciate any input -- and I've looked at similar questions on Stackoverflow but have not found anything compelling. NOTE: I am aware of R, SciPy and Octave and java libraries that make calls to them -- I am looking for a Java native library or set of libraries that can together provide the features I'm looking for. NOTE: The topics covered in such a class typically include: one-samle and two- sample tests and confidence intervals for means and medians, descriptive statistics, goodness-of-fit tests, one- and two-way ANOVA, simultaneous inference, testing variances, regression analysis, and categorical data analysis.


title: Identify different periods of variance in a time series

text: I have a time series $x_t$ which may go through different phases of volatility. One example might be some stock that has high variance from 9 AM to 11 AM, low variance from 11 AM to 2 PM, and then high variance again afterwards. Is there a way to identify these different periods of variance? I am thinking of taking a sliding window of length $L$, computing the variance in that window and running change point detection on that, but I think this requires me to know the distribution of the estimated variance. Another idea I had was to take the variance of the windows and try to fit a Markov chain to it, but I don't know beforehand how many states there should be. Sorry, the motivating example is not just in finance. I would like to have a way to model "risk" in a time series, which might be risk associated with a financial asset. Another example is in wind energy production - energy produced from wind turbines is highly volatile and unpredictable, but it would be nice to have a learning algorithm that could identify what days of the month is variance of wind energy production high and when is it low


title: Effects of multiple within-subjects factors when not all combinations of levels are available

text: My question: how can I analyze the effects of multiple within-subjects factors when my design does NOT include the all combinations of their factor levels, as would be required by a standard factorial ANOVA? Here's my design: Each participant is trained to discriminate 2 categories in each of 2 domains. In one of the domains, they receive treatment A, and in the other, treatment B. Which treatment goes with which domain is selected randomly for each participant, as is which domain is shown first, so all 4 combinations of domain orders and treatment orders are possible and equally likely. My DV is discrimination performance, which comes from a test administered right after treatment. There is one test for the first treatment / domain and another for the second. There are actually 3 categories, rather than 2, in each domain, but a given participant only sees 2 of them. Which 2 they see is determined by another factor, called 'distance', taking levels 'near' and 'far'. If they are in the near condition, they get the 2 categories designated 'near' for both of the 2 domains, and likewise for far. So, ignoring the random order of presentation of the 2 domains and the 2 treatments, I have 2 within-subject factors (domain, treatment) and 1 between- subject factor (distance). The problem is that I don't have the factorial combinations of the 2 within-subject factors: a given participant only sees 2 of the 4 possible combinations, e.g. domain 1 with treatment A and domain 2 with treatment B. So, if I plug the above factors into a standard ANOVA, it won't work because some combinations of the within-subjects factors don't appear. Finally my question. What's the right model to analyze this data? Someone mentioned to me that this is a "nested design", but what I found on the net regarding this term didn't seem very relevant. Any help would be appreciated!


title: How to calculate the pairwise LD for the given data?

text: I have the following data, which is the output from the MS Hudson software.               segsites: 6      positions: 0.1256 0.3122 0.3218 0.4970 0.5951 0.7943      001010      110101      010100      001010      010100       I want to make an R function to calculate the R-Squared across pairs separated by <10% (the difference between positions of SNPs must be < 0.10) of the simulated genomic region. How would I go about doing this?


title: Is single precision floating point so bad?

text: I've been looking at some of the packages from the High perf task view dealing with GPU computations, and given that most GPU seem to be an order of magnitude stronger at performing single precision arithmetics than DP ones, i was wondering:   1. Why none of the packages gives more control to the user on the type of precision required ? I can see many applications in statistics where SP arithmetics (i.e. number coded with 7 digit accuracy) are good enough for practical use (if i am overestimating the gains involved, let me know).    2. Is python more flexible on this ? If so why ? I don't see why the absence of a 'single' type in R would make including such an option (together with a warning) in say GPUtools or magma impossible (though i'll be happy to be shown wrong). PS: i'm specifically thinking of applications were the numbers are already dimension-wise scaled and centered (so that Chebychev's inequality is binding) dimension-wise.


title: Plotting standard errors

text: This is a very basic questions however I have not been able to find an answer. When you plot your regular bar plot with your standard errors of the means for comparison do you plot the standard errors and means from the raw data ?   Or do you plot the ones predicted by the model you are fitting. Thank you in advance for your help. I usually use the latter please advice. Alfred


title: Have I correctly specified my model in lmer?

text: I have scoured lots of help sites and am still confused about how to specify more complicated nested terms in a mixed model as well. I am also confused as the use of `:` and `/` and `|` in specifying interactions and nesting with random factors using `lmer()` in the `lme4` package in `R`. For the purpose of this question, let's assume I have accurately portrayed my data with this standard statistical model: $$ Y_{ijk} = u + \text{station}_i + \text{tow}_{j(i)} + \text{day}_k + (\text{station}\times \text{day})_{ik} + (\text{tow}\times\text{day})_{j(i)k} $$ `station` is fixed, `tow` and `day` are random. `Tow` is (implicitly) nested within `station`. In other words, I'm hoping that my model includes Station(i,fixed), Tow(j,random,implicitly nested within Station), Day(k,random), and interaction between Tow and Day, and the interaction between Day and Station. I have consulted with a statistician to create my model and at this time believe it to be representative of my data, but will also add in a description of my data for those who are interested at the bottom of my post so as not to clutter. So far what I've been able to piece together is the following in `lmer`:               lmer(y ~ station + (1|station:tow) + (1|Day) + (1|station:day) + (1|tow:day),           data=my.data)      Does this accurately depict my statistical model? Any suggestions for how to improve my code if it does not read correctly? **I've bolded the specific terms I'm having difficulty specifying in my lmer formula** **#1. tow nested within station when tow is random and station is fixed**   I'm confused, however about differentiating between nested and interaction terms that are random using `:` and `/` . In my above example, I have `(1|station:tow)` in which I'm hoping reads tow nested within station. I've read conflicting comments on various sites whether or not I should be using `:` or `/` here within the random `(1|...)` format of `lmer`. **#2. The interaction between station and day when station is fixed and day is random**   I then have `(1|station:day)` but this time I'm hoping it reads the interaction between station and day. It seems like I could use station*day to account for the individual effects of station and day as well as their interaction (rather than including each of the three terms separately as I do above), but I don't see how to specify this when one is fixed and the other is random. Would `station*(1|day)` do that? **#3. The interaction between tow and day (both random) when tow is nested in station (fixed)** Then lastly, I have `(1|tow:day)` which I'm hoping reads the interaction of `tow` and `day`, but I'm wondering if I need to specify again that tow is nested (implicitly) in station? I am new to both `R` and `lmer` and statistical modeling and greatly appreciate the trouble of thorough explanations in any responses to my questions if possible.$$$$ More details on my data: I am asking whether concentrations of plankton vary across a physical front in the nearshore ocean. I have three stations, inshore, within, and offshore of this front. Station is thus fixed. At each station, I take three replicate plankton tows (from which I sort, count, and get a concentration in terms of # of bugs per meter cubed of water). Tow is random: in three tows I hope to account for the general variability in plankton at that particular station. Tow is intrinsically nested in station as each tow does not have a unique ID (123,123,123 is the ID for tows at each station). I then did this on multiple, independent days with a new front that had formed. I think I can think of Day as a blocking factor? Day is random as repeating this on multiple independent front days is attempting to capture variability from day to day and be representative of all days where this front is present. I want to know about the interaction terms to see if Tows change in variability from day to day and if stations always yield similar data or does it depend on the day? Again, thank you for your time and help, I appreciate it!


title: Determining uncertainty propagation for an overdetermined set of equations

text: I have this problem of the following format. \begin{align*} x_L &= X\left(\frac{dx_L}{dX}\right) + Y\left(\frac{dx_L}{dY}\right) + Z\left(\frac{dx_L}{dZ}\right)\\ &\\ x_R &= X\left(\frac{dx_R}{dX}\right) + Y\left(\frac{dx_R}{dY}\right) + Z\left(\frac{dx_R}{dZ}\right)\\ &\\ y_L &= X\left(\frac{dy_L}{dX}\right) + Y\left(\frac{dy_L}{dY}\right) + Z\left(\frac{dy_L}{dZ}\right)\\ &\\ y_R &= X\left(\frac{dy_R}{dX}\right) + Y\left(\frac{dy_R}{dY}\right) + Z\left(\frac{dy_R}{dZ}\right)\\ \end{align*} where all of the left hand sides ($x_L$, $x_R$, $y_L$, and $y_R$ are known, and have known uncertainty values) and while each multiplicative term is known ($dx_L/dX$, $dx_R/dX$, etc) their uncertainties are difficult to approximate (numbers come from calibration), and will probably have to be assumed negligibly small. $X$, $Y$, and $Z$ are unknowns. As you can see I have 4 equations and three unknowns, so the system is overdetermined. I am not sure how to propagate uncertainty in the left hand terms to the values $X$, $Y$ and $Z$. I have used a method outlined in this thread for a properly specified set of equations. Propagation of uncertainty through a linear system of equations Except i take 4 separate uncertainties, for each of the four combinations of three equations (1,2,3) (1,2,4) (1,3,4) (2,3,4). Depending on which 3 equations are used, the uncertainties can vary by two orders of magnitude. I have almost no ideas on how to combined these, since the actual solutions to the overdetermined set takes a least squares fit between the four solutions. Any help would be appreciated. For those interested, what we actually have here is a grid of relative displacements in image space to real space for stereo particle image velocimetry measurements to determine 3d velocity profiles. $L$ and $R$ indicate displacements measured in the Left and Right cameras respectively. ($x_L$) is horizontal image displacement in pixels on the left camera, ($y_L$) is the vertical image displacement in pixels on the Left camera etc. Edit: to Whomever made my formatting pretty, Thanks!


title: How to calculate E(Y) for a given X in a mixed logistic regression?

text: I made a mixed regression model and I would like to calculate the probability of $Y=1$ for a given value of the predictor for each estimation of the random effect. I am not sure, given it is a logit link function (binary dependent variable), if I can use that: $$\frac{\exp(B_0 + B_1X + u)}{(1+\exp(B_0 + B_1X +u))}$$ where $u$ is the estimate of the random effect. Do anyone has an idea?


title: Estimating population correlation based on current data and a previous study

text: Assume:   * A previous study looking at the relationship between $X$ and $Y$ obtained a correlation of $r = 0.50$ using a sample of $n = 100$. The raw data is not available.   * The current study also looking at the relationship between $X$ and $Y$ obtained a correlation of $r = 0.45$ with $n = 50$. How would you do the following tasks:   1. Give your best estimate of the population correlation between $X$ and $Y$ assuming the two studies are estimating the same correlation.   2. Give your best estimate of the population correlation between $X$ and $Y$ in the current study assuming that the first study is slightly different than the current study (e.g., it used a different measurement procedure, a different type of sample, etc.). Of course the weight given to the previous study would depend on perceived similarity with the second study. Thus, are there standard ways of quantifying similarity between studies in such calculations?


title: How to verify a distribution is normalized?

text: I have a question in which it asks to verify whether if the Uniform distribution (${\rm Uniform}(a,b)$) is normalized.   1. For one, what does it mean for any distribution to be normalized?    2. And two, how do we go about verifying whether a distribution is normalized or not?  I understand by computing $$ \frac{X-\text{mean}}{\text{sd}} $$ we get normalized _data_ , but here it's asking to verify whether a _distribution_ is normalized or not.


title: Diagnostic plots for count regression

text: What diagnostic plots (and perhaps formal tests) do you find most informative for regressions where the outcome is a count variable? I'm especially interested in Poisson and negative binomial models, as well as zero-inflated and hurdle counterparts of each. Most of the sources I've found simply plot the residuals vs. fitted values without discussion of what these plots "should" look like. Wisdom and references greatly appreciated. The back story on why I'm asking this, if it's relevant, is my other question. Related discussions:   * Interpreting residual diagnostic plots for glm models?   * Assumptions of generalized linear models   * GLMs - Diagnostics and Which Family


title: Binary classification using radial basis kernel SVM with a single feature

text: Is there any interpretation (graphical or otherwise) of a radial basis kernel SVM being trained with a single feature? I can visualize the effect in 2 dimensions (the result being a separation boundary that is curved rather than a linear line. (e.g http://en.wikipedia.org/wiki/File:Kernel_Machine.png). I'm having trouble thinking of what this would be like if your original data only had a single feature. What would the boundary line look like for this case?


title: Reclassifications based on altered probability distributions

text: Background. Using non-negative matrix factorisation on mutation data, I classified mutations in samples   based into clusters defined by signatures. The mutations in each sample are broken down into the three signatures, where each signature is defined by the probability distribution of seeing particular trinucleotide contexts being enriched in it, having assumed that all trinucleotides occur at equal frequency (i.e, 1/64)               Signature 1 = 1/64.p1(a) + 1/64.p2(b) + 1/64.p3(c) + 1/64.p4(d) + 1/64.p5(d)     Signature 2 = 1/64.p6(a) + 1/64.p7(b) + 1/64.p8(c) + 1/64.p9(d) + 1/64.p10(d)     Signature 3 = 1/64.p11(a) + 1/64.p12(b) + 1/64.p13(c) + 1/64.p14(d) + 1/64.p15(d)      Each sample's mutations are broken down into Signature 1,Signature 2 and Signature 3 based   on those probability distributions. However, factoring in the actual frequencies they occur in the human exome with leads to   the probability distribution for each signature changing such that they are represented by...               Signature 1 = x1.p1(a) + x2.p2(b) + x3.p3(c) + x4.p4(d) + x5.p5(d)     Signature 2 = x1.p6(a) + x2.p7(b) + x3.p8(c) + x4.p9(d) + x5.p10(d)     Signature 3 = x1.p11(a) + x2.p12(b) + x3.p13(c) + x4.p14(d) + x5.p15(d)      where x1, x2 , x3, x4 and x5 are the weights for the normalised probabilities of a,b,c,d,e. Is there any way of recalculating the breakdown of mutations in each sample based on the altered distributions using the information from the original distributions and breakdowns? Alternatively, how valid is it to recalculate the mutation counts using weighting for the difference between the two probability distributions before running non-negative matrix factorisation? Cheers.


title: How to interpret Weka Logistic Regression output?

text: Please help interpret results of logistic regression produced by `weka.classifiers.functions.Logistic` from the WEKA library. I use numeric data from WEKA examples:               @relation weather          @attribute outlook {sunny, overcast, rainy}     @attribute temperature real     @attribute humidity real     @attribute windy {TRUE, FALSE}     @attribute play {yes, no}          @data     sunny,85,85,FALSE,no     sunny,80,90,TRUE,no     overcast,83,86,FALSE,yes     rainy,70,96,FALSE,yes     rainy,68,80,FALSE,yes     rainy,65,70,TRUE,no     overcast,64,65,TRUE,yes     sunny,72,95,FALSE,no     sunny,69,70,FALSE,yes     rainy,75,80,FALSE,yes     sunny,75,70,TRUE,yes     overcast,72,90,TRUE,yes     overcast,81,75,FALSE,yes     rainy,71,91,TRUE,no      To create the logistic regression model I use the following command:               java -cp WEKA_INS/weka.jar weka.classifiers.functions.Logistic -t WEKA_INS/data/weather.numeric.arff -T WEKA_INS/data/weather.numeric.arff -d ./weather.numeric.model.arff      Here are what the three arguments mean:               -t <name of training file> : Sets training file.     -T <name of test file> : Sets test file.      -d <name of output file> : Sets model output file.      Running the above command produced the following output:               Logistic Regression with ridge parameter of 1.0E-8     Coefficients...                   Class     Variable                    yes     ===============================     outlook=sunny           -6.4257     outlook=overcast        13.5922     outlook=rainy           -5.6562     temperature             -0.0776     humidity                -0.1556     windy                    3.7317     Intercept                22.234          Odds Ratios...                   Class     Variable                    yes     ===============================     outlook=sunny            0.0016     outlook=overcast    799848.4264     outlook=rainy            0.0035     temperature              0.9254     humidity                 0.8559     windy                   41.7508               Time taken to build model: 0.05 seconds     Time taken to test model on training data: 0 seconds          === Error on training data ===     Correctly Classified Instances          11               78.5714 %     Incorrectly Classified Instances         3               21.4286 %     Kappa statistic                          0.5532     Mean absolute error                      0.2066     Root mean squared error                  0.3273     Relative absolute error                 44.4963 %     Root relative squared error             68.2597 %     Total Number of Instances               14               === Confusion Matrix ===      a b   <-- classified as      7 2 | a = yes      1 4 | b = no      Questions:   1. First section of the report:              // Coefficients...                        Class     Variable                    yes     ===============================     outlook=sunny           -6.4257     outlook=overcast        13.5922     outlook=rainy           -5.6562     temperature             -0.0776     humidity                -0.1556     windy                    3.7317     Intercept                22.234          * Do I understand right that `Coefficients` are in fact weights that are applied to each attribute before adding them together to produce the value of class attribute `play` equal to `yes`?   2. Second section of the report:              // Odds Ratios...                        Class     Variable                    yes     ===============================     outlook=sunny            0.0016     outlook=overcast    799848.4264     outlook=rainy            0.0035     temperature              0.9254     humidity                 0.8559     windy                   41.7508          * What is the meaning of "Odds Ratios"?     * Do they all also relate to class attribute `play` equal to `yes`?     * Why is the value `outlook=overcast` so much bigger than the value of `outlook=sunny`?   3. Confusion matrix              === Confusion Matrix ===      a b   <-- classified as      7 2 | a = yes      1 4 | b = no          * What is the meaning of "Confusion Matrix"? Thanks a lot for your help!


title: Running a directional (one-tailed) Mann-Whitney U test in python

text: As described in this thread, the implementation of _Mann-Whitney U_ provided by scipy doesn't seem to return the direction of the effect (i.e., it's always one-tailed). I other words, it returns if the two populations differ significantly or not but not which one is greater. **Might there be a different python implementation that can be used for both two-tailed and one-tailed testing?** A possible solution is to wrap around R's MWU implementation like in this example, but I'd rather not have that extra dependency on R and `rpy2`.


title: Normal distribution for discrete values

text: I understand that a normal distribution is a continuous distribution, but I came across this question in my test set that baffles me. **Question:**   The toy cars produced in a factory have an average weight of 20 lb. The standard deviation of the individual toy cars is 3 lb. We know the weights of each toy car follow a normal distribution. Compute the proportion of 200 toy cars that weigh 14 lb each.


title: Randomized Results Sampling

text: I'm taking a look at randomized results and have a question about an end case. Specifically, suppose my scheme is as follows: 1) I want to ask question Q_A. The probability of answering Q_A in the affirmative is p. 2) I have innocuous question Q_I. 3) I have everyone in my sample flip a coin. If it comes up heads, they answer Q_A truthfully. If it comes up tails, they flip the coin AGAIN and if on the second flip it comes up heads, they answer "Yes" and if it comes up tails, they answer "NO". 4) I'm interested in estimating p. To do this I though of the following: $$P(Yes) = P(Yes| Coin 1= H) *P*(Heads) + P(Yes|Coin1 = T)*P(Tails) \Rightarrow $$ $$P(Yes) = p*(1/2) + (1/2)*(1/2) \Rightarrow $$ $$p = 2*(P(Yes)-1/4)$$ However, that seems to behave badly at certain points. For instance, suppose 9/10 of the class raised there hand. I then get p=1.3 which is impossible. I tried plotting the likelihood but it seems to go from 0 at -1/4 and maxes out at 3/4. Thank you for your help!


title: Expected value of posterior vs. success probability

text: ## Context Suppose I have two models, $H_1$ and $H_2$ for which I know the prior probabilities $p(H_1)$ and $p(H_2)$. Furthermore, I know the class-conditional distributions $p(x|H_1)$ and $p(x|H_2)$ of a random variable $X \in \mathbb{R}^n$. I get to observe a realization of $X$, call it $x_0$. The posterior distribution of the models given the observation I just made can be obtained via Bayes' formula: $$p(H_i|x_0) = \frac{p(x_0|H_i)p(H_i)}{p(x_0)}$$ And the Bayes-optimal decision rule is the (very intuitive) MAP rule: $$\hat{H} = \operatorname*{arg\,max}_{H_i} p(H_i|x_0)$$ Now, let's say I'm interested in how certain I am when I'm deciding for $H_1$. The probability that I'm correctly identifying $H_1$ when indeed $X$ is generated from model $H_1$ is: $$P_{suc,H_1} = \int_{R_1}p(x|H_1)dx,$$ where $R_1$ is the decision region associated to $H_1$, i.e. $$R_1 = \{x \in \mathbb{R}^n : \operatorname*{arg\,max}_{H_i} p(H_i|x_0) = H_1\}$$ ## Problem Now that I've introduced the terminology & notation, let me come to the actual question ;) In my case, $R_1$ has a complicated shape (it's a pretty much arbitrary convex polytope), and I'm just not able to perform the integration in order to obtain $P_{suc,H_1}$. I'm looking for alternative ways to get this quantity. My intuition is as follows: **the probability of success is related (equal?) to the expected value of the posterior under the hypothesis $H_1$**. Mathematically: $$P_{suc, H_1} \longleftrightarrow \mathbb{E}_{p(x|H_1)}[p(H_1|x)]$$ or: $$\int_{R_1}1p(x|H_1)dx \longleftrightarrow \int_{\mathbb{R}^n}p(H_1|x)p(x|H_1)dx$$ Actually, I've already managed to prove that the two are not equal (in the left integral, replace $1$ with $p(H_1|x) + p(H_2|x)$; in the right integral, use $\mathbb{R}^n = R_1 + R_2$; some terms cancel out and you can analyze what's left over.) ## Questions I guess I really have two questions:   1. What's the interpretation of the quantity $\int_{\mathbb{R}^n}p(H_1|x)p(x|H_1)dx$? Intuitively, it should relate to $P_{suc,1}$, but I'm not sure how / why.   2. Is there any other way I could compute the error probability by integrations over some simpler regions (even if the integrand becomes more complicated)?


title: How to estimate a dynamic Tobit model

text: I have data which correspond to a corner solution. The Tobit-model seems to be adequate for this data. However, I also wants to control for a baseline variable (t-1) and unobserved heterogenity. This is normally not permissable in the Tobit model. I understand that Wooldridge have solved the problem in this paper: http://download.clib.psu.ac.th/datawebclib/e_resource/trial_database/WileyInterScienceCD/pdf/JAE/JAE_4.pdf (page 49-50) My problem is that I am weak in mathematics. Thus, I don't understand what Wooldridge proposes in the paper just mentioned. Can somebody please explain to me (in plain language), step by step, the process which I should follow in order to build a dynamic Tobit model.


title: Coefficient of determination of a orthogonal regression

text: ## Context I must find the general form of a linear equation. The $X$ and $Y$ values are the location coordinates of touches on a screen. I want to find the best fit line, described by equation $mX + nY + p = 0$ , so I need to find $m, n, p$. I want to be able to make regression of horizontal lines as well as vertical lines, or diagonal lines. I decided to make **Orthogonal Distance Regression** (ODR) instead of Ordinary Least Square (OLS) one. I am using SciPy library (python langage), and I've written below function to make the job. It is using Single Value Decomposition (SVD). I've tested it, and it is giving excellent results. Now, I want to get the coefficient of determination $R^2$. ## Coefficient of determination According to Wikipedia, the definition is $$ R^2 = 1 - \frac{SS_{res}}{SS_{tot}} $$ where $SS_{res}$ is the sum of squares of residuals, and $SS_{tot}$ is the total sum of squares. To apply it in my case I considered: $$ SS_{res} = \sum_id((x_i,y_i);(\hat{x}, \hat{y}))^2 $$ where $ d((x_i,y_i);(\hat{x}, \hat{y}))$ is the distance of point $(x_i,y_i)$ to the best fit line. and $$ SS_{tot} = \sum_id((x_i,y_i);(\overline{x}, \overline{y}))^2 $$ where $\overline{x}$ is the mean of all the $x_i$ and $\overline{y}$ is the mean of all the $y_i$. ## My question Is it the correct way to compute the coefficient of determination? I am asking because it looks optimistic, in the sense that it is very often closed to 1.


title: On the Fly Time Series modeling

text: I'm dealing with a system which monitors and records a time series (half hourly) which I plan to use to build a double seasonal time model (if possible using something that already exists, such as tbats from R's forecast). My issues is related to the fact that I don't plan to build this model very often. Let's say I have three weeks of data which allowed me to take into account both the daily and weekly seasonalities. Two months later, it's a Tuesday, it's 10AM and I want to make a forecast for 12AM. Do I need to use all the data from the moment I made the model to this current date, or is there a way I can use a shorter period (such as the current day data) ? Speed is paramount and if I need to make a forecast every half an hour, I would really like only using a reduced period of time to make the forecast. I'm basically looking for a method for dealing with growing time series as most literature deals with a static data sets.


title: Is there a Discriminative RBM implementation for Python?

text: First, sorry for my poor english.. I'm trying to implement Discriminative Restricted Boltzmann Machine (paper), but I'm having some trouble in computing `p(y|x)` because its value explodes while computing the product of `n_hidden` `1+exp(o_yj(x))` values. Is there any way to solve this problem? Whatever, what I need is just an implementation of DRBM, I'm looking for some python codes which implemented DRBM.


title: dimension reduction of discrete numerical data

text: I have a bivariate discrete numerical dataset and would like to reduce its dimensions to a single variable. A 9 x 8 table of counts of the (x,y) data values is:                    0  1  2  3  4  5  6 7        -----------------------     8 | 10 38 58 79 50 19 15 1     7 | 22 61 67 86 43 19  7 0     6 | 44 70 44 44 19  6  2 0     5 | 45 33 29  8  5  2  1 0     4 | 62 15  8  4  1  0  1 0     3 | 12  6  7  2  2  1  0 0     2 |  1  2  1  0  0  0  0 0     1 |  3  0  0  0  0  0  0 0     0 | 10  2  2  0  1  0  0 0      Something that seems reasonable is to fit a curve from (0,0) to (7,8) and project the data onto this [in the study (0,0) is least knowledgeable, (7,8) most knowledgeable]. The reduced variable would be the distance from (0,0) to (7,8) after projection. Does anyone know how to do this? The X variable is the sum of number of facts correctly recalled after reviewing some materials. The Y variable is a sum of correct responses to TRUE/FALSE questions relating to those materials. Another example of these type of data is the following                     0   1   2  3          -------------       8 | 5  82 143 37       7 |21 139 163 23       6 |23 122  86 15       5 |17  63  29  6       4 |21  29   8  1       3 | 6  13   3  0       2 | 3   4   2  0       1 | 2   0   0  0       0 | 4   0   0  0      in this case the X variable has values 0 to 3. Suggestions for other methods also appreciated.


title: Proving increasing function using expectation

text: Given $c>0$ and define a function $f:(0,\infty)\mapsto \mathbb{R}$ by \begin{equation} f(\sigma)=\frac1{\sqrt {2\pi}\sigma}\int_{-\infty}^\infty\max\{-c,\min\{x,c\}\}^2e^{-\frac{x^2}{2\sigma^2}}. \end{equation} We can see $f$ as the second moment of $X_c=\max\{-c,\min\{X,c\}\}$ where $X\sim N(0,\sigma^2)$. I could prove that $f$ is increasing using derivative which is very tedious. I'd like to see $f$ is increasing using expectation approach. Could anyone help? Thanks in advance


title: Logistic Regression with time slices

text: I'm using logistic regression to predict student retention in an online course. I have a data of student interactions within a web platform of an online course. The course spans 6 weeks, with new lecture resources and new assignments uploaded at the beginning of each week. The weekly assignments due at the end of each week. Students can watch lecture videos, view/write forum posts, and do assignments. Students can drop out of the course any time (i.e. no longer interacting within the course platform). Each week, I want to predict the likelihood of a student staying within the course in the next week ( **stay in the next week=1, out in the next week=0** ). The predictors are the number of times the student watch the lecture videos ( **video_views** ), the number of posts the student read ( **posts_read** ), the number of posts the student wrote ( **post_written** ) and the student's score of this week assignment ( **score** ). I'm thinking of building 6 models using logistic regression, for each week. But I also want to make a connection between, say, week 6's model and week 5's model. Can such connection be shown if I use percentage of cumulative values for each predictors (e.g. cumulative number of assignment score out of total assignment score for the whole course) instead of weekly values? Should I just build a single model with the course week ( **course_week** ) as another predictor? This is my first time using logistic regression, and I'm not sure about putting a time variable in a logistic regression model.


title: How to plot non-uniform bar chart using GNUplot?

text: I want to plot using GNUplot a bunch of data that is in this form:               123 130 5      ... where the first number is the x starting value, the second is the x ending value, and the 3rd number is the y value. Ideally I would plot these in an accumulative/stacked manner. Is there any easy way to do this with GNUplot?


title: Based only on these sensitivity and specificity values, what is the best decision method?

text: If I have the following `sensitivity` and `specificity` values, what is the best decision we can say in this case?               sensitivity     specificity     -----------     -----------       66.3             74.7       87.2             65.9       56.4             76.4       79.5             94.3      Thanks.


title: Binary or Multinomial Logistic Regression in SPSS: Interpretation and Reference Categories

text: I am trying to analyze my data using Multinomial Logistic Regression whereby my dependent variable is a clinical outcome (sick vs healthy) and 1 independent variables (Factors) are in several categories. The problem I have is trying to figure out how I can set one of the category as a reference group in SPSS. I have done the inverse analysis by switching dependent and factor variables so I can set the reference but after now that I think about it, I realize it doesn't make much sense. I also compared the OR values from MLR to that of a 2x2 analysis but it is vastly different. I also tried Binary Logistic Regression and created dummy variables for each category but I didn't have sensible values either. Any advice? _EDIT:_ SPSS command and output [I have set the first category (1) in my independent variable as reference]               NOMREG Group (BASE=FIRST ORDER=ASCENDING) BY Status     /CRITERIA CIN(95) DELTA(0) MXITER(100) MXSTEP(5) CHKSEP(20) LCONVERGE(0)          PCONVERGE(0.000001) SINGULAR(0.00000001)      /MODEL      /STEPWISE=PIN(.05) POUT(0.1) MINEFFECT(0) RULE(SINGLE) ENTRYMETHOD(LR)         REMOVALMETHOD(LR)      /INTERCEPT=EXCLUDE      /PRINT=PARAMETER SUMMARY LRT CPS STEP MFI.      ![output](http://i.stack.imgur.com/Ay75U.png)


title: Structural Equation Models - building with a large sample, and testing with a small one?

text: I'm planning a study using Structural Equation Modelling to test different accounts of language learning. I would like to study a group of children with language difficulties. However, they are very hard to recruit and I'm aware that these models typically require at least 200 participants. One possibility would be (1) Test SEMs on a large cohort of typically developing children - who are much easier to recruit (n = 500). (NB I'm imagining 3 or 4 different models each with about 3 or 4 latent variables). (2) Apply these models to a smaller cohort of language impaired individuals (n = 80) I'm imagining that in the second stage one would keep the coefficients from stage one, and would use some kind of procedure to manipulate the values of the latent (exogenous) variables to maximise model fit. Research questions would be   * Which model best describes the data for the typically developing children?   * Does this model also best describe the data for the language impaired children? I'm wondering if this is statistically feasible and whether anyone has tried this approach before. Ta.


title: Which test do I run to check inter-rater agreement for multiple ordinal variables?

text: I'd like to check the inter-rater agreement between the six raters for the six ordinal variables scored for each of the 24 score sheets passed to each rater (6 ordinal variables per sheet, 24 sheets per rater, 6 raters in this study 6 X 24 X 6). Which inter-rater agreement test should I run on this ordinal data? Do I analyze all the six variables across the six raters' 24 sets of data together (i.e. is there an all-encompassing correlation test? 6 variables X 24 score sheets X 6 raters)? Or should I check the level of agreement among the six raters for all 24 score sheets per variable to focus more on agreement per variable (1 variable X 24 score sheets X 6 raters)? Thank you for your time and consideration in reading this question? I really appreciate your feedback!


title: Ordinal Data analysis

text: I have to analyze the results of a survey. The questions regard different motivators that the interviewees encounter ("I do it to help my colleagues" as a motivator for instance); all the answers are ordinal (1=strongly disagree -> 5=strongly agree) and the aim of this study is to see which factor (or set of factors) (nationality, age, site they work in, gender, role in the organization, center they belong) is influencing the answer. I have looked in a lot of websites and lectures (ordered logit/probit model, discrete dependent variables, Likert scale, correlation, regressions, and so on) but it seems that my study could match a lot of ways to analyze these data. So far I only have 15 interview results, which is not a lot: we expect around 100 answers at the end. I would like to have causal inference: their current behavior (motivations) should help us to draw a conclusion (causal effect, does their nationality influence their motivation factors?) regarding the factors of motivations. Could someone tell me which is the best way to treat these results?


title: Excel: which products are most frequently ordered together? (clustering question)

text: I'd like to recruit your help in coming up with an Excel-based method to analyse a set of raw ordering data where each item is on its own row. So, in the data below, order 111 contains two part numbers, ABC and DEF, in quantities of 2 and 3, respectively.               Order Number | Part Number | Quantity |     111          | ABC         |     2    |     111          | DEF         |     3    |     222          | ABC         |     1    |     333          | ABC         |     4    |     333          | DEF         |     4    |      You can think of ABC as peanut butter and DEF as jelly. Usually they're ordered in a 1:1 ratio, but someone may be a heavy peanut butter user and only order peanut butter (as in order 222). The part number array is large, let's say in the several hundreds (peanut butter, jelly, apple butter, fluffernutter, lemon curd, whatever). I'd like to be able to answer:   * For a given part number (ABC, for example) what is the most commonly ordered partner part number (peanut butter, ABC, is usually paired with jelly, DEF, 83% of the time, apple butter, 3%, by itself 14%) and   * What is the ratio of the queried part number against its partner? (peanut butter:jelly = 0.8 because jelly goes faster usually)


title: Feature selection from wavelet transformation in R

text: I am new to wavelets. Currently, I am developing a prediction model using time series data. I am using the `wavelets` package in R. I am taking part of the time series, doing a wavelet transformation, and trying to predict the final outcome.   1. I am doing a discrete wavelet decomposition. How should I select the required features/coefficients out of those coefficients?   2. Should I do a PCA to extract the features?   3. I am getting a lot of coefficients estimated to be zero. Should I remove those coefficients?  Any pointers would be helpful.


title: linearity of a time series

text: I am currently trying to correct forecast data using Kalman filter (python). I do not know where to start. I wanted to know how can I do a test to Know if my time series is linear or non linear? Is there a test?


title: Validation of model with relatively sparse data

text: I am not a statistician, so please forgive my lack of precision with terminology in the following background and question. I have done a fair amount of reading and research about the following, but feel as though I may not be headed in the right direction. I have a mathematical model for which I would like to estimate some unknown parameters in the form of distributions. I have created a hierarchical Bayesian model and am able to get estimates of the parameter distributions at both an individual and population level. This seems to be working fine. To conduct this estimation, I have been using all of the available data I can find, which constitutes about 10 studies. The number of samples in each dataset is 'small' (about 5-10 points) and the number of points in each dataset is different. My concern is that I am using all of the data for 'calibration' and am not presently setting any aside for validation. I have read about cross validation and various k-fold techniques for optimal parameter estimation. However, I am unclear on whether this type of technique makes sense for my situation. All of the materials I have read describe how these techniques work for single parameters (not distributions) and for data where each set is of approximately the same size. Do you have any recommendations about ways of validating the type of model I have described and identifying the 'best' set of parameter distributions? Also, any pointers to relevant (introductory) papers or books would be appreciated. Thank you very much for your help.


title: How to define Confidence Intervals from a distribution of 1000 t-statistics?

text: I have a vector of length 300 containing some kind of values (say, scores of a math test). The distribution is not normal. I want to test if the average score of a small group (30 samples) is significantly different from the mean of all the remaining 470 individuals. I performed a Welch's t-test obtaining a t-statistic. In order to validate my results I performed a bootstrapping: I computed 1000 t-tests from a random creation of subgroups (again,30 vs. 470) generating a distribution of 1000 t-statistics. My question is, how can I define the Confidence Intervals for my t-statistics distribution and check if the t-statistics of the original subgroup is outside the CI?


title: Sampling brownian motions

text: I wish to sample standard linear Brownian motions on the interval $[0,1]$. I do this by dividing the interval into $n$ equal sub-intervals, deciding $B(0)=0$, and letting $B\left(\frac{k}{n}\right)=B\left(\frac{k-1}{n}\right)+\mathcal{N}\left(0,\frac{\sigma^2}{n}\right)$ for $k\ge 1$, after deciding for some $\sigma>0$. I do this $s$ times. My question is as follows: how large should $n$ and $s$ be so I'll feel comfortable enough to say that the sampled Brownian motion represent, in some sense, the real distribution of Brownian motions? In fact, there are two questions here:   * How large should $n$ be so that the random walks should resemble, in some sense, Brownian motions?   * How large should $s$ be so that the distribution of the sampled fractional Brownian motions would resemble, in some sense, the real distribution of such fractional Brownian motions?


title: Questions about text mining code in R

text: I am trying to do text mining, and seeing the code, I have the whole picture about what it's trying to do about the text. But the problem is on some specific part of code, I don't know why the format is this way, and what the parameters present. So do you have some suggestions about references or books about R language so that I can check what is this function used for and the interpretation of parameter in this function? Below are several questions about doing text mining.   1. Why we need two `[ ]` in the third line when turn the `termDocumentMatrix` into a matrix?              cand=c("Romney","Obama")     tdm<-list(name=cand,tdm=s.tdm)     #s.tdm is TermDocumentMatrix of a text.     tdm.dm<-t(data.matrix(tdm[["tdm"]]))      2.               filepath<-"C:/e"     cor.score<-if(length(grep("http|html",filepath))){cor.score<-Corpus(URISource(filepath))}else{score.cor <- generateSpeechDocCorpus(filepath)}      This sentence is trying to see if the filepath is URL or not, I understand using "grep" to check if filepath has string "http" or "html", but why we need sentence "length" outside grep? I am confused. AND for the last term in the code:               generateSpeechDocCorpus(filepath),      I can also use               Corpus(DirSource(directory=filepath,encoding="ANSI"))      to achieve the same purpose. So what is the difference between generateSpeechDocCorpus and Corpus?


title: How do I determine/compute a cutoff for chance level/not chance level?

text: I am in the process of describing my research design for my dissertation and ran into a roadblock. In my design, I am converting 20 y/n responses from 190 participants to two dichotomous groups: 1). Chance level and 2). Above chance level. If chance level is 50% or 10 responses correct/incorrect, how do I determine what would not be chance? How do I determine the cutoff for placing an individual into the above chance level group? I have heard a few things such as 25% over chance is no longer chance. But even here, I am not sure how to calculate it...is 50% the base or is 100% the base? I could really use your help with this matter. Thank you for your consideration. ====================================================================================== Hi Glen, I asked a question about probabilities and cut-offs, to which you replied with a great answer. I appreciate the table you provided as it gives me a choice. However, now my dilemma is that I would like to know the source so I can cite it. My mentor wants to know my source! Skootz


title: CLT for sum of products of Gaussian variables

text: I have two multivariate Gaussians $\mathbf{x} \sim \mathcal{N}(\mu_x, \Sigma_x^2)$ and $\mathbf{y} \sim \mathcal{N}(\mu_y, \Sigma_y^2)$. I assume that the covariances are diagonal. I am interested in modelling the distribution of their product $a = \mathbf{x}^T\mathbf{y}$, ideally as another Gaussian. Since it is a sum of products and can be rewritten as $a = \sum_i x_iy_i$. Does the central limit theorem hold in this case?


title: Multiplication and division of values with geometric standard deviation

text: What is the geometric standard deviation of a value, which is the result of dividing two independent values, each of which has its own geometric standard deviation ? It is a frequent situation in science, that the the signal of the machines is in linear relationship to the logarithm of the analysed quantity, so it is appropriate to use geometric standard deviation to characterize the the error of the measured quantity. (I calculated geometric mean and geometric standard deviation from triplicate measurements for each value). Now I want to divide two such quantity values (one by another) to obtain a relative concentration. What will be the geometric standard deviation of the relative concentration ? (I assume, that the two original values are independent - I do not want to deal with covariances etc.) * * * **EDIT:** In desperation, I came up with a theory, but would appreciate a confirmation from a matematician, because I am a molecular biologist. I wanted to do x=a/b. I took my triplicate measurements and made all possible combinations of values, which are a1/b1, a1/b2, a1/b3, a2/b1, a2/b2, a2/b3, a3/b1, a3/b2, a3/b3. Then I calculated geometric mean and geometric standard deviation, as if these 9 results were 9 independent measurements. **IS THIS STEP CORRECT ???** Another approach: I assume, that "lognormal distribution, geometric means, multiplication/division of geometric means, logarithms of geometric standard deviations" are all in the same analogy to "normal distribution, arithmetic means, addition/substraction of arithmetic means, arithmetic standard deviations". For my "x=a/b" , it leads me to these formulas: ![enter image description here](http://i.stack.imgur.com/qF6hF.jpg) The resulting equation gives me the **same results** , as when I calculate from those 9 independent combinations described above. I calculate in Excel, which shows many decimal numbers, so it demonstrates very nicely, that the results of both calculations are identical. **Can someone confirm the formulas are correct, please ?**


title: How can I explain this correlation? (two trendlines with same orientation)

text: I obtained a graph by comparing two values (chemical elements) and I obtained two general trend lines parallel with each other (one from the black circles/triangles/etc., and the other from the grey ones). How can be explained in terms of graph observations / statistical correlation? There is a bimodal correlation? ![enter image description here](http://i.stack.imgur.com/HmbJ8.png)


title: Repeated measure two-way ANOVA (non-parametric)

text: Non-parametric ANOVA – a hot topic that is unanswered. There are many questions on this topic online. However, they all seem to end in a debate and no definite answer or clear explanation (that I can relate to my data set). The only good answers I can find seem to relate to one treatment, that was repeated over time. How do I handle two treatments? More importantly, should this be handled as a linear model? I cannot meet assumptions of homogeneity, normality, equal sample size, and transformations were not helpful. Data setup:   * Soil cores growing plants.   * Subject – about 33 samples. A core will consist of a treatment combo ( _n_ = 5 to 6; repeats)   * `Treatment01` - clay, sand, loamy sand   * `Treatment02` - water, fertiliser   * Time - an sample is analysed from each core sample, every week, for 3 months.   * Variable - Nutrients (mg/kg soil), moisture, pH, mass of soil etc. For simplicity, let's do carbon. Here are my failed attempts on R. This problem rests in the statistical theory. `Anova(lm(Carbon ~ Treatment01 * Treatment02, data), type = "3")`   or `model1 <\- lm(Carbon ~ Treatment01 * Treatment02)` Then something like...               print(lsmeans(model1, list(pairwise ~ Treatment01)), adjust = c("tukey"))     print(lsmeans(model1, list(pairwise ~ Treatment02)), adjust = c("tukey"))     print(lsmeans(model1, list(pairwise ~ Treatment01 | Treatment02)), adjust = c("tukey"))     print(lsmeans(model1, list(pairwise ~ Treatment02 | Treatment01)), adjust = c("tukey"))      Any suggestions? What are the disadvantages in pulling the data set apart to do a series of Friedman or Kruskal–Wallis tests (depending on how much I pull this data set apart?)


title: Why can't I use the variance of the sample average in the Central Limit Theorem for the weak-stationary process?

text: Under mild conditions $\dfrac{\bar{X}-\mu}{\sqrt{\sigma^2/n}}$ approaches the standard normal (where $\sigma^2$ is the process variance, not the marginal variance $\sigma^2_x$). Why is the denominator not the standard error for the sample mean of dependent data, specifically $\sqrt{\frac{\sigma_x^2}{n}\times [1 + 2*\sum_{i=1}^{n-1}(1-i/n)*\rho_i]}$ ?


title: Mixed effects modelling; what to do when model is over-specified?

text: I'm trying to use mixed-effects modelling to analyse some data. There are a number of variables that I need to specify within the model, two of which are between-participants (`x1` and `x2`) and two of which are repeated-measures (`z1` and `z2`). I'm interested in individual differences both generally and in relation to the repeated-measure variables. From what I understand, according to Barr's (2012) 'Keep it maximal' paper, the `lmer` syntax should be:               y ~ x1*x2*z1*z2+ (1+(z1*z2)|ID)      When I try to run this in `R`, it fails, I think the error message is telling me that I've over-specified the model:               Error in checkZdims(reTrms$Ztlist, n = n, control, allow.n = FALSE) :      number of observations (=289) <= number of random effects (=292) for term ((z1 * z2) | ID); the random-effects parameters and the residual variance (or scale parameter) are probably unidentifiable      Is there a way to maximally specify my model, or do I need to restrict my random effects term? If so, what would the new term be?


title: How to plot a learning curve based on a sequence of date stamped successes and failures?

text: **I have** data on results of surgical operations:   * date of operation (in YYYY-MM-DD format)   * result of operation (0 for successful operation, 1 for unsuccessful operation)  Date is usually used only for case ordering. **I need to**   * plot a learning curve. (I had found that the resulting curve called **cumulative average**.) x should be the number of cases y - probability of failure (usually decreases).   * get the data on current performance quality  **The question is:**   * What is the way to achieve this in R (with ggplot2 if possible)?   * Is $\dfrac{successful}{total}*100\%$ will be sufficient statistic for this case or something more advanced is usually required?


title: Quantiles from the combination of normal distributions

text: I have information on the distributions of anthropometric dimensions (like shoulder span) for children of different ages. For each age and dimension, I have mean, standard deviation. (I also have eight quantiles, but I don't think I'll be able to get what I want from them.) For each dimension, I would like to estimate particular quantiles of the length distribution. If I assume that each of the dimensions is normally distributed, I can do this with the means and standard deviations. Is there a pretty formula that I can use to get the value associated with a particular quantile of the distribution? The reverse is quite easy: For a particular value, get the area to the right of the value for each of the normal distributions (ages). Sum the results and divide by the number of distributions. **Update** : Here's the same question in graphical form. Assume that each of the colored distributions is normally distributed. ![The same question in graphical form](http://i.stack.imgur.com/nP07r.png) Also, I obviously can just try a bunch of different lengths and keep changing them until I get one that's close enough to the desired quantile for my precision. I'm wondering if there is a better way than this. And if this is the right approach, is there a name for it?


title: Support vector regression for multivariate time series prediction

text: Has anyone attempted time series prediction using support vector regression? I understand support vector machines and partially understand support vector regression, but I don't understand how they can be used to model time series, especially multivariate time series. I've tried to read a few papers, but they are too high level. Can anyone explain in lay terms how they would work, especially in relation to multivariate time series? Thanks in advance. EDIT:To elaborate a bit, let me try to explain with a stock price example. Say we have stock prices for N days. Then, for each day we could construct a feature vector, which, in a simple case, could be be the previous day's price and the current day's price. The response for each feature vector would be the next day's price. Thus, given yesterday's price and today's price the objective would be to predict the next days price. What I don't understand is, say we have six months training data, how would you give greater emphasis to the more recent feature vectors?


title: Caculate Type II error for Binomial distribution

text: Basically this is a question if I got the general idea behind the Type II Error right. I’ll set up a little reproducible example: Given $p_0=0.6$ with an $\alpha$-level of $\alpha=0.1$. As Binomial tests are conservative with respect to $\alpha$ we first calculate (using `R`) the actual $\alpha$ used: ## Calculating the actual $\alpha$-level: $P(X< c_1)= 0.0123$               for(i in 1:10){ # lower bound         if(sum(dbinom(0:10, 10, 0.6)[1:(i+1)])<0.05){             c <- sum(dbinom(0:10, 10, 0.6)[1:(i+1)])         }         else break         print(c)     }      $P(X>c_2)=0.0464$               for(i in 10:1){ # upper bound         if(sum(dbinom(0:10, 10, 0.6)[11:i])<0.05){             d <- sum(dbinom(0:10, 10, 0.6)[11:i])         }         else break         print(d)     }      Hence we have an actual $\alpha$ of $\alpha=P(X<c_1)+P(X>c_2)=0.0123+0.464=0.587$. ## Calculating the critical values $c_1$ and $c_2$: The critical values $c_1$ and $c_2$ are: $c_1=3$               qbinom(0.0123, 10, 0.5)      $c_1=8$               qbinom(1-0.0464, 10, 0.5)      ## Calculating Type II Error for $p_1=0.5$ and $p_2=0.7$: $\beta_{0.5}=0.817$               pbinom(8, 10, 0.5) - pbinom(3, 10, 0.5)      $\beta_{0.7}=0.84$               pbinom(8, 10, 0.7) - pbinom(3, 10, 0.7)      ## Logic behind this: The Type II Error $\beta$ with respect to a specific value of a parameter $\vartheta\in H_1$ is defined as $\beta_\vartheta = P(c_1\le X\le c_2|\vartheta\in H_1)$. The critical constants $c_1$ and $c_2$ are calculated under $H_0$ but the probability of $X$ being greater than or equal $c_1$ and less than or equal to $c_2$ is calculated under $H_1$. Is this correct?


title: Sample session for analysis of conjoint data in R

text: While I have successfully used R for Cluster Analysis and Correspondence Analysis, I'm new to the use of R-project for Choice-Based Conjoint. Per recommendations in this post, I'm looking at the following [R] Packages for purposes of analyzing choice-based conjoint data:   * ChoiceModelR   * mlogit   * AlgDesign for constructing choice sets   * prefmod for analysing paired comparison data   * conf.design for constructing factorial designs For reference purposes, I'm curious whether there is a sample R-project session showing the use of one or more of these packages to:   * Specify attributes and attributes levels   * Analyze choice-based conjoint data to list out respondent-level utility values   * Simulate market share Has anyone seen a sample R-project session of this nature? Or can someone who has created such a session post it?


title: Errors-in-variables regression: is it valid to pool data from three sites?

text: I recently had a client come to me to do a bootstrap analysis because an FDA reviewer said that their errors-in-variables regression was invalid because when pooling data from sites the analysis include pooling data from three sites where two sites included some samples that were the same. **BACKGROUND** The client had a new assaying method they wanted to show was "equivalent" to an existing approved method. Their approach was to compare the results of both methods applied to the same samples. Three sites were used to do the testing. Errors-in-variables (Deming regression) was applied to the data at each site. The idea is that if the regression showed the slope parameter to be close to 1 and the intercept near 0 this would show that the two assaying techniques gave nearly the same results and hence the new method should be approved. At site 1 they had 45 samples giving them 45 paired observations. Site 2 had 40 samples and site 3, 43 samples. They did three separate Deming regressions (assuming a ratio of 1 for the measurement errors for the two methods). So the algorithm minimized the sum of squared perpendicular distances. Separate regressions were done at each site and a pooled regression was also done using all the data from all three sites. In their submission the client pointed out that some of the samples used at sites 1 and 2 were the same. In the review the FDA reviewer said that the Deming regression was invalid because common samples were used which causes "interference" that invalidates the assumptions of the model. They requested that a bootstrap adjustment be applied to the Deming results to take account of this interference. At that point since the client did not know how to do the bootstrap I was brought in. The term interference was strange and I was not sure exactly what the reviewer was getting at. I assumed that the point really was that the because the pooled data had common samples there would be correlation for the common samples and hence the model error terms would not all be independent. **THE CLIENT'S ANALYSIS** The three separate regressions were very similar. Each had slope parameters close to 1 and intercepts near 0. The 95% confidence interval contained 1 and 0 for the slope and intercept respectively in each case. The main difference was a slightly higher residual variance at site 3. Furthermore they compared this to the results from doing OLS and found them to be very similar (in only one case did the confidence interval for the slope based on OLS not contain 1). In the case where the OLS CI for the slope did not contain 1 the upper bound of the interval was something like 0.99. With the results being so similar at all three sites pooling the site data seemed reasonable. The client did a pooled Deming regression which also lead to similar results. Given these results I wrote a report for the client disputing the claim that the regressions were invalid. My argument is that because there are similar measurement errors in both variables the client was right to use Deming regression as a way to show agreement / disagreement. The individual site regressions had no problems of correlated errors because no samples were repeated within a given site. Pooling data to get tighter confidence intervals. The pooling which included use of the the common samples twice might produce positive correlation between the residuals for those common samples which would mean that the confidence intervals for the regression parameters would be too narrow (the estimated model residual variance biased on the low side). This difficulty could be remedied by simply pooling the data with the common samples from site 1 say left out. Also the three individual site models do not have the problem and are valid. This seems to me to provide strong evidence of agreement even without the pooling. Furthermore the measurements were taken independently at sites 1 and 2 for the common sites. So I think that even the pooled analysis using all the data is valid because the measurement errors for a sample at site 1 are not correlated with the measurement errors in the corresponding sample at site 2. This really just amounts to repeating a point in the design space which should not be a problem. It does not create correlation / "interference". In my report I wrote that a bootstrap analysis was unnecessary because there is no correlation to adjust for. The three site models were valid (no possible "interference" within sites) and a pooled analysis could be done removing the common samples at site 1 when doing the pooling. Such a pooled analysis could not have an interference problem. A bootstrap adjustment would not be necessary because there is no bias to adjust for. **CONCLUSION** The client agreed with my analysis but was afraid to take it to the FDA. They want me to do the bootstrap adjustment anyway. **MY QUESTIONS** A) Do you agree with (1) My analysis of the client's results and (2) my argument that the bootstrap is unnecessary. B) Given that I have to bootstrap the Deming regression are there any procedures SAS or R that are available for me to do the Deming regression on the bootstrap samples? EDIT: Given the suggestion of Bill Huber I plan to look at bounds on the errors-in-variables regression by regression both y on x and x on y. We already know that for one version of OLS the answer is essentially the same as errors-in-variables when the two error variances are assumed to be equal. If this is true for the other regression then I think that will show that the Deming regression gives an appropriate solution. Do you agree? In order to meet the client's request I need to do the requested bootstrap analysis that was vaguely defined. Ethically I think it would be wrong to just provide the bootstrap because it doesn't really solve the client's real problem, which is to justify their assay measuring procedure. So I will give them both analyses and request at least that they tell the FDA that in addition to do the bootstrap I did inverse regression and bounded the Deming regressions which I think is more appropriate. Also I think that analysis will show that their method is equivalent to the reference and the Deming regression is therefore adequate also. I plan to use the R program that @whuber suggested in his answer to enable me to bootstrap the Deming regression. I am not very familiar with R but I think I can do it. I have R installed along with R Studio. Will that make it easy enough for a novice like me? Also I have SAS and am more comfortable programming in SAS. So if anyone knows a way to do this in SAS I would appreciate knowing about it.


title: Can statistical significance be obtained in an underpowered study?

text: Excuse my ignorance. I tried to look up the questions but couldn't find it. I had a discussion today whether it is possible that a study result will have a p<0.05 yet that would be caused by chance (i.e. a possible underpowered study - if I got all the terms correctly). My thoughts were that if there is a statistically significant difference (proven by p=0.03) then the study is automatically powered enough to detect that difference. Thank you!


title: Using n-grams to find data that does not 'belong'

text: Recently I posted a question over in CS.SE dealing with methods of classifying data. Essentially the problem is that I have a collection of strings (100's of thousands). Most of these strings are the same "type"-- in my initial case they are street addresses. However, there are a few strings in this collection that don't belong (city names, zip codes, or even completely unrelated nonsense like a number e.g. '53231232'). My goal is to identify these strings that don't belong. One suggestion on CS.SE was to use n-gram analysis. For addresses especially this seemed like it could work, considering addresses tend to follow similar patterns (start with a number, end in St/Rd etc.). Just to get me started I counted bi-gram frequencies from my collection of strings. So I ended up with something like this (first 5 rows sorted by frequency descending):               ST #          104958     STREET #      89884     RD #          83532     AVE #         82628     ROAD #        65682      There were a total of 1.5 million distinct bigrams, and about 7.6 million total. My problem now is that I want to pull each string and evaluate it against these results. For example if I take "1 Microsoft Way" I can compare the frequencies of each bigram:               # 1               7409     1 Microsoft       1     Microsoft Way     1     Way #             16363      While '1 Microsoft' and 'Microsoft Way' may only appear once, '# 1' and 'Way #' appear a relatively large number of times and thus I would hope this is classified as an address. Similarly, if a zip code is examined: 98052               # 98052     1     98052 #     1      There's nothing to suggest this is an address and I'd hope it would be classified as 'not an address'. My problem is how can I "combine" the frequencies I've found for individual addresses into a meaningful statistic? I've though about computing z-scores for each bi-gram in the sample address, but this still leaves me with 1 z-score per bi-gram and I don't know how I would combine these. I've though about maybe finding the average z-score but I'm pretty sure that has no statistical basis. I'm also not sure how the fact that the original corpus contains the erroneous data influences what statistical tests I can and cannot use (if it does at all). Essentially what I want to do is have a probability that "1 Microsoft Way" is an address based on the bigrams it contains and the frequency of those bigrams in the larger corpus. So say, "1 Microsoft Way" gets a .90. At the same time, I'd hope "98052" gets a very low probability, say .01. Can anyone suggest any statistical tests or other methodologies that will allow me to achieve this?


title: Do random forest variable importance measures take into account the interactions?

text: Do random forest measures of variable importance (mean change of accuracy, mean change of Gini index) take the interactions into account? I think I know how we come up with the variable importance plot (by permuting each of the predictors), and it doesn't seem that random forest captures the interaction. Does anybody have another point of view? Thanks.


title: The relationship between the number of support vectors and the number of features

text: I ran an SVM against a given data set, and made the following observation: If I change the number of features for building the classifier, the number of resulting support vectors will also be changed. I would like to know how to explain this kind of scenario.


title: Panel regression with proportion as both dependent and independent

text: I need to do a panel regression where my data is proportions for both dependent and one of the independent variables. I do not have the underlying data for the proportion so I cannot expand the data and use logistic regression. Any ideas on how to handle this? Poisson regression maybe?


title: SPSS version 19 omits one covariate (of 8) in a binary logistic regression

text: When I run a binary logistic regression using SPSS 19, the program can analyse seven covariates, but consistently omits the eighth. This occurs regardless of which covariate it is. Is this a limitation of this version of the software, or can someone suggest what I may be doing wrong?


title: Error "system is computationally singular" when running a glm

text: I'm using the robustbase package to run a glm estimation. However when I do it, I get the following error:               Error in solve.default(crossprod(X, DiagB * X)/nobs, EEq) :        system is computationally singular: reciprocal condition number = 1.66807e-16      What does this mean/indicate? And how can I debug it? PS. If you need anything (the formula/specification or data) to answer, I'll gladly provide it.


title: Analysis of proportions over time

text: My knowledge of statistics is limited and I am looking for resources to read on the matter if possible. Anyways, I am currently trying to estimate a confidence interval for a proportion over time. The specific example pertains to a trader's Win % as he make more trades. My current idea is to take the first $n$ sample data to calculate a proportion, then see how that changes over the amount of trades. I want to ultimately calculate the variability of the win %. My confidence interval would be calculated as such: $$p \pm \sqrt[2]{\frac{1}{n} * p(1-p)}$$ I know this is clearly the wrong approach to look into but I would appreciate it if I am pointed to in the right direction.


title: Is there a one-tailed test for difference in data variance in a mixed model?

text: I would like to test if the variation of my response variable is higher for one level of an explanatory variable than for the other. However I am not interested in the total variation, but only variation towards higher values of the response variable. My model (using the package `nlme` in R): `lme(y~factor + covariate, random=~1|Site)` I should specify that I do not want to look at differences in the residual variance but in the variation of the observed values.


title: Regression through PCA

text: I am working with principal component analysis (PCA) to have a regression model. Let's say we have 3 variables that we get every month, and so far we gathered 10 months of data. We have built the PCA on the raw data matrix.               p = prcomp(na.omit(rdatat), center=TRUE,scale=TRUE)          loadings=p$rotation[]          pca=cbind(p$x[,1], p$x[,2], p$x[,3])      Now we can put `p$x[,1]` as input of the regression and the result will be in the form of $Y=$ intercept$+ K\,$Component1, where $K$ is a coefficient from regression. How can I build Component1 the month after in order to predict $Y$? I have the raw values of the 3 variables I use in PCA, but how can I use them?


title: Transform Data to Desired Mean and Standard Deviation

text: I am looking for a method to transform my dataset from its current mean and standard deviation to a target mean and a target standard deviation. Basically, I want to shrink/expand the dispersion and scale all numbers to a mean. It doesn't work to do two separate linear transformations, one for standard deviation, and then one for mean. What method should I use? Could the solution possibly be applied to an example where a point 1.02 in a dataset with SD .4 and a mean 0.88 is transformed when I adjust the mean of the dataset to 0.5 and the SD to 0.1667? What is the new value of the point?


title: Can chi square be used to compare proportions?

text: I've read that the chi square test is useful to see if a sample is significantly different from a set of expected values. For example, here is a table of results of a survey regarding people's favourite colours (n=15+13+10+17=55 total respondents):               red,blue,green,yellow          15,13,10,17      A chi square test can tell me if this sample is significantly different from the null hypothesis of equal probability of people liking each colour. Question: Can the test be run on the proportions of total respondents who like a certain colour? Like below:               red,blue,green,yellow          0.273,0.236,0.182,0.309      Where, of course, 0.273+0.236+0.182+0.309=1. If the chi square test is not suitable in this case, what test would be? Thanks! Edit: I tried @Roman Luštrik answer below, and got the following output, why am I not getting a p-value and why does R say "Chi-squared approximation may be incorrect"?               > chisq.test(c(0,0,0,8,6,2,0,0),p = c(0.406197174,0.088746395,0.025193306,0.42041479,0.03192905,0.018328576,0.009190708,0))              Chi-squared test for given probabilities          data:  c(0, 0, 0, 8, 6, 2, 0, 0)      X-squared = NaN, df = 7, p-value = NA          Warning message:     In chisq.test(c(0, 0, 0, 8, 6, 2, 0, 0), p = c(0.406197174, 0.088746395,  :       Chi-squared approximation may be incorrect


title: Likert Scale Analysis

text: What are some interesting techniques that can be used to analyze Likert data? As a frame of reference, I have constructed a survey with about 50 items which are intended to assess the attitudes of the survey participants towards two very broad categories. Each item is on a 5 point scale. The survey will be offered multiple times, thus pre and post testing is possible, although not necessarily the aim of the survey.


title: How to calculate pvalue for the Spearman correlation test

text: I'm trying to code up the spearman correlation coefficient test of R in C++.               > x = runif(100)     > y = runif(100)     > cor.test(x,y, method="spearman")     Spearman's rank correlation rho          data:  x and y      S = 165182, p-value = 0.9306     alternative hypothesis: true rho is not equal to 0      sample estimates:         rho          0.008808881      I manage to get S and rho, but I cannot find any documentation to go from S to the pvalue. Is there any GSL distribution I can use for this? BTW, I don't want to test significance using Student's t distribution or permutations. Thanks!


title: What is an 'atom' and what are 'atomic weights'?

text: I have come across the following statement: > A notable feature of the Hierarchical Dirichlet Process is that all > Dirichlet Processes' $G_j$ share the same set of atoms and only the atom > weights differ. This is a result of the almost sure discreteness of the top- > level DP. What is meant by atom and atom weights? Googling gets to me Wiki articles about measure theory and I'm not able to understand them and also not sure if they are what I'm looking for. Could anyone provide a simple explanation of their meaning?


title: Annual time series count data where the dependent variable is a count averaging 3,000 and no zeros

text: I need your assistance on time series count data. I got some annual time series data I want to run, however the dependent variable is a count (number of deaths) while the independent variables are either continuous or discrete. Ideally count data is modeled using standard Poisson or NB model (mostly its usually cross sectional data) or panel Poisson or panel NB for panel data. For time series the use of the standard models of Poisson or NB will not take care of serial correlation according to literature I have reviewed by Cameron and Trivedi (1999), Quddus (2008) and others (see attached documents). They therefore recommend the use of other versions of the models eg INAR(1) models. However, none of the literature explains how exactly this is to be done in any of the soft wares. Stata 13 does not avail that option neither does Eviews 5 (softwares that I have access to). How then is it to be executed? Or which command or adjustment to the standard models are to be applied when running the data? Or is it just a matter of including a lag of Y as an explanatory variable and executing Poisson? And will that get me a value for the binomial thinning operator? I have tested for stationarity and found that my dependent count variable is I(0) while the independent variables are either I(0), I(1) or I(2). How then do I proceed if I cannot use standard Poisson or NB and neither can I use other time series models such as ARIMA because of the normality assumption which allow for negative values? Your quick guidance will be highly appreciated.


title: Intuitive understanding of regularization

text: As I understand from this link regularization used to reduce overfitting of model.   1. Is overfitting bad when we have really a lot of data?    2. I don't understand why "very large weights fit the training data very well"?   3. Is regularization always needed?


title: Ljung Box test in R

text: I'm new to stats. I want to get clarity about Ljung Box test. Its a omnibus test. I'm using R to test autocorrelations of a time series. Using below command               Box.test(LogReturns, lag = 10, type = c("Ljung-Box"))      To test autocorrelations upto lag 10, can I infer with above single command or I need to run it in loop from lag = 1 to 10? R documentation is not clear to me. Any suggestion? Thanks


title: Compared results of 2 goodness-of-fit chi-square test

text: I am conducting `Goodness-of-fit` test on 2 cities using chi square.                                 | Males    Females                       | --------------     Test 1 : Boston   | 50000    10000     say chi square = 29.45 with df=1      Similiarly,                                 | Males     Females                       | --------------     Test 2 : NYC      | 500000    20000     say chi square = 455.45 with df=1          Can I say since chi-square for NYC was much higher than that of Boston, the 'unknown' relationship between city and sex ratio is stronger in NYC than Boston?


title: Transition from pure mathematics to statistics via operational research

text: I am currently studying an M.Sc. in pure mathematics (my dissertation focuses on category/set theory)). However recently I have decided that I wish to transition to a much more applied field. During my undergrad I studied almost all pure courses, with only one statistics module. I have enrolled in an M.Sc. in operational research for next year where I will be studying a mix of operational research and statistics courses. All of this is in the UK. My question is, should I choose to, will it be possible/realistic for me to apply to do a Ph.D. in some sort of statistics? Not only in terms of actually getting accepted but also getting funding/actually being able to do this? Thanks for any help


title: Does Greenhouse-Geisser correction influence the effect size estimation, the statistical significance treshold, or both of them?

text: I discussed recently with one of my colleagues about the Greenhouse-Geisser correction for non-sphericity. I told him first that this correction is very conservative and that it would be more efficient to perform an alternative analysis not relying on the sphericity assumption. He answered that, in his case, the effects were statistically significant after applying the G-G correction and that consequently he did not have to care about the conservative nature of his analysis. As a reply, I said that since G-G correction influences the statistical significance, this correction might maybe influence the bias and efficiency of the effect size estimations. However, this first reply was intuitive and after having thought about that later, I think now that I was probably wrong. Indeed, the correction might affect only the statistical significance treshold as if we worked with a higher alpha than the conventional .05. However, while writing my question I just remember of formula allowing to compute effect size estimation using the F statistic and its df. So maybe I was also wrong in the way I re-thought to the question. The G-G correction might affect both the statistical significance treshold and the effect size estimation. What is your opinion about that? Thanks in advance for your answers. EDIT: I just tried directly to see if the _df_ adjustement following G-G correction leads to a dramatical change in the partial eta squared value. I used the following formula to compute the partial eta squared from the _F_ statistic and its _df_ : partial eta squared = _F_ /( _F_ + ( _df2_ /*df1*)) I used the following values _F_ (2, 10) = 12.534 before G-G correction (from this webpage on sphericity) and I found a partial eta squared = .71483974. The non-sphericity correction is epsilon = .638. So the corrected _df1_ = 2*.638 = 1.276 and the corrected _df2_ = 10*.638 = 6.380. So after correction, we can used the following values _F_ (1.276, 6.380) = 12.534 to compute a partial eta squared = .71483974. So we find exactly the same estimation of effect size. Why is it the case? When we look closer to the formula partial eta squared = _F_ /( _F_ + ( _df2_ /*df1*)), we can see that it is the ratio of the two _df_ which is considered to compute the effect size. So since the G-G correction consists in multiplying each _df_ by epsilon we obtain the following formula if we try to use the corrected _df_ : partial eta squared = _F_ /( _F_ + ( _df2_ *epsilon/ _df1_ *epsilon)). So here we can see that the two epsilon in the formula can be removed and that the formula remains the same regardeless of the adjustment of the _df_. In conclusion, it seems from this example that G-G correction does not affect the point estimate of the effect size. Its influence on statistical significance seems to be mainly, if not exclusively, a matter of statistical significance treshold adjustment. Maybe this observation is trivial for everyone else but maybe it will help someone. Nevertheless, my problem is not entirely solved with this illustration. Indeed, it seems that G-G correction does not affect point estimate. But what about confidence intervals? Intuitively I would say that if such a correction affect the comutation of the p-value it would also influence the compuation of confidence intervals according to the link between the two. In the webpage abovementioned, it is stated that after applying the G-G correction we have _p_ = .009 compared with the _p = .002_ observed before applying the correction. So if there are no mistakes in this webpage, the G-G correction might affect statistical significance throught an increase of the p-value. Consistently we should observed an increase in the width of the confidence interval for the effect size. But does the G-G correction also decrease the statistical significant treshold? Is there someone who knows something about that? Thanks again


title: Are these statements about p-values correct?

text: I've seen the following equivalent statements about p-values:   * _"It is the probability of wrongly rejecting the null hypothesis if it is in fact true."_ [1]   * _"The P value or calculated probability is the estimated probability of rejecting the null hypothesis (H0) of a study question when that hypothesis is true."_ [2] Some people I know insist these are wrong, and say these are equivalent to type I error. But the type I error is in fact determined by the significance level. Are these right or wrong? [1] http://www.stats.gla.ac.uk/steps/glossary/hypothesis_testing.html#pvalue [2] http://www.statsdirect.com/help/default.htm#basics/p_values.htm


title: How to fit a mixture of Gamma distributions to the PMF of a discrete distribution?

text: I have a PMF of some discrete distribution that has been numerically computed.   Note that I do **not** have any _samples_ to work with here, so techniques like Maximum-Likelihood and Expectation-Maximization don't apply. I only have the PMF of the discrete distribution itself, which is simply a long, nonnegative vector whose components sum to 1. The discrete distribution looks reasonably well-modeled as a mixture of _N_ gamma distributions ( _N_ is known). What's a reasonable way to go about fitting the mixture to it? The only way I can think of is to hand-code my own coordinate- or gradient- descent algorithm, but it seems too much effort (both on my part and in terms of the amount of computation necessary). Is there a better way? (While not necessary, a SciPy or MATLAB/Octave example could be extremely helpful. I'm hoping for a method I can code myself in a language like C++, but I realize that might not be practical, so I'm interested in other approaches as well.) * * * Some example data as requested, in case it helps:               0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0.00000378547917956     0.000254067618914     0.00156479688482     0.0044414187977     0.00881560818165     0.0150067507346     0.0250934783012     0.0364480196843     0.0446846535887     0.0481736403324     0.0473452833494     0.0436535252283     0.0387132874982     0.0337816696454     0.0295972032879     0.0267001698978     0.0279827988189     0.0362165748226     0.0486471989886     0.0602335084185     0.0672898116143     0.0684179033513     0.0641605623675     0.0561730620339     0.046373159578     0.0363756352803     0.0272651661276     0.0196061621026     0.0135630878889     0.009043229377     0.00581909570544     0.00361712103199     0.00217346975615     0.0012631914959     0.000710408159256     0.000386756023978     0.000203892354011     0.00010411799478     0.0000515137355117     0.0000246996739328     0.0000114793542013     0.00000517223941954     0.00000225965576783     0.000000957343112229     0.00000039337449953     0.000000156784391692     0.0000000606172441131     0.0000000227362899619     0.00000000827373414225     0.00000000292123458756     0.00000000100077057752     0.000000000332677663195     0.00000000010731171507     0.0000000000335905747662     0.0000000000102032826632     0.00000000000300759417371     0.000000000000860422844084     0.000000000000238586927992     0.0000000000000640598685209     0.0000000000000160982338571     0.00000000000000432986979604     0.0000000000000008881784197     0     0


title: Penalized Bayesian quantile regression with LASSO and adaptive LASSO penalty

text: I have some questions about penalized Bayesian quantile regression with LASSO and adaptive LASSO penalty:   1. Would you give me a detailed outline for the formula, especially as used in Bayesian inference?    2. What package in R is easy to use and suitable for this topic?   3. Would you give me R code for running inference Bayesian quantile regression with LASSO penalty?


title: Why does a fixed-effect OLS need unique time elements?

text: The `plm` function of the `plm` library in R is giving me grief over having duplicate time-id couples, even when I'm running a model that I don't think should need a time variable at all (see reproducible example below). I can think of three possibilities:   1. My understanding of fixed effects regression is wrong, and they really do require unique time indices (or time indices at all!).   2. plm() is just being overly-finicky here and should relax this requirement.   3. The particular estimation technique that plm() uses--the within transformation--requires time indices, even though the order doesn't seem to matter and the less computationally-efficient version (including dummies in a straight-up OLS model) doesn't need them. Any thoughts?               set.seed(1)     n <- 1000     test <- data.frame( grp = as.factor(rep( letters, (n/length(letters))+1 ))[seq(n)], x = runif(n), z = runif(n) )     test$y <- with( test, 2*x + 3*z + rnorm(n) )     lm( y ~ x + z, data = test )     lm( y ~ x + z + grp, data = test )          require(plm)     # Model fails if I don't specify a time index, despite effect = "individual"     plm( y ~ x + z, data = test, model = "within", effect="individual", index = "grp" )      # Create time variable and add it to the index but still specify individual FE not time FE also     library(plyr)     test <- ddply( test, .(grp), function(dat) transform( dat, t = seq(nrow(dat)) ) )     # Now plm() works; note coefficients clearly include the fixed effects, as they match the lm() version above     plm( y ~ x + z, data = test, model = "within", effect="individual", index = c("grp","t") )      # Scramble time variables and show they don't matter as long as they're unique within a cluster     test <- ddply( test, .(grp), function(dat) transform( dat, t = sample(t) ) )     plm( y ~ x + z, data = test, model = "within", effect="individual", index = c("grp","t") )      # Add a duplicate time entry and show that it causes plm() to fail     test[ 2, "t" ] <- test[ 1, "t" ]      plm( y ~ x + z, data = test, model = "within", effect="individual", index = c("grp","t") )       **Why this matters** I'm trying to bootstrap my model, and when I do the requirement that the index-time pairs be unique is causing headaches which seem unnecessary if (2) is true.


title: Loss function that relates ROPE with HDI?

text: In Doing Bayesian Data Analysis (link to the book) and Bayesian Estimation Supersedes the t-Test, J. Kruschke proposes using the following criterion to reject or accept the null hypothesis in a location experiment:   * Accept null hypothesis if the 95% HDI falls completely inside the ROPE   * Reject the null hypothesis if the 95% HDI falls completely outisde the ROPE Presumably ROPE is given by a desired effect size (correct me if I am wrong). Is there a way to connect the above decision rule with the optimization of a posterior expected loss? If so, what would is the loss function that corresponds to this criterion?


title: Variable Construct (Dimension) Insertion (Representation) in SmartPLS

text: I have a question please. I am working on SmartPLS. The question is regarding variable construct representation. Any variable could have construct including dimensions. For example, the entrepreneurial orientation construct is : risk- taking, pro-activeness, and innovativeness. All the three dimensions are under the entrepreneurial orientation variable. The question is whether I need to represent each dimension alone as a separate latent variable in SmartPLS, or I can represent them (dimensions) all with a single latent variable? Thank you,


title: Can nonstationarity be told from the autocorrelation function?

text: Here "stationarity" means the first and second moments don't change over time. From a page of Time Series: Theory and Methods, by Peter J. Brockwell, Richard A. Davis > In this chapter we shall examine the problem of selecting an appropriate > model for a given set of observations $\{X_t t = 1, ..., n > \}$. If the > data (a) exhibits no apparent deviations from stationarity and (b) **has a > rapidly decreasing autocorrelation function** , we shall seek a suitable > ARMA process to represent the mean-corrected data. If not, then we shall > first look fo r a transformation of the data which generates a new series > with the properties (a) and (b). From a page of Time Series: Theory and Methods, by Peter J. Brockwell, Richard A. Davis > **Trend and seasonality** are usually detected by inspecting the graph of > the (possibly transformed) series. However they **are also characterized by > sample autocorrelation functions which are slowly decaying and nearly > periodic respectively**. From wikipedia > non-stationarity is often indicated by an autocorrelation plot with very > slow decay. **My questions:** Why should a stationary ARMA process have a rapidly decreasing autocovariance function? Should a stationary AR process also have a rapidly decreasing autocovariance function? Why is non-stationarity often indicated by an autocorrelation plot with very slow decay and values near 1 at small lags? Why are trend and seasonality characterized by autocorrelation functions that are slowly decaying and nearly periodic, respectively?


title: EM Algorithm to maximize the censored likelihood

text: I want to implement the EM algorithm to maximize the following censored likelihood when the data is normal: ![enter image description here](http://i.stack.imgur.com/4TQbt.png) The complete data likelihood is: ![enter image description here](http://i.stack.imgur.com/sVDuZ.jpg) resulting in the expected complete- data log-likelihood ![enter image description here](http://i.stack.imgur.com/aG8eX.png) where Zi's are missing and distributed from Normal(theta,1). The EM sequence is: ![enter image description here](http://i.stack.imgur.com/c7qpM.png) My R-code to achieve the maximum of theta is:                #initial values:     ybar=0;sd_y=1     a=1;n=30;m=20       theta=0       dif=iter=1          while (dif>10^(-4)){     theta=c(theta,m*ybar/n+(n-m)*(theta[iter]+     dnorm(a-theta[iter])/(1-pnorm(a-theta[iter])))/n)     iter=iter+1     diff(theta[iter:(iter+1)])}      The Algorithm is not working. Any helps would be appreciated. Thanks!