title: Advice on formula for factorial, blocked (possibly) repeated measures anova in R

text: Pardon my ignorance, but I really need some help figuring out how to analyze a dataset. Basically the experiment is replicated among 8 blocks, with two treatments (`a` and `b`) and a control within each block. Unfortunately (not my design) the experiment is not fully factorial (there are no a*b treatments). The experiment was run for a year, and we sampled (species richness) once every three months. I'm interested in testing for significant treatment affects, seasonal affects, and season*treatment interactions. Here are my questions:   1. Do I need to run two separate analyses for treatment `a` and `b` since the experiment was not fully factorial?   2. I'm not quite sure how to code this formula (I'm using `aov` in R). My best guess is something like this:               aov(richness~treatment_A + season + block + (season*treatment_A)       and then the same for treatment b. But I don't know if that's valid since the same plots are measure 4 times (once per season). Any advice would be _REALLY_ appreciated. As you can see, my ignorance is deep: I don't even know the proper terminology here.


title: How to interpret the results of a t-test?

text: I recently performed an experiment observing ants. For five days I measured the amount of food they ate, 25 ants in the fruits group and 25 ants in the vegetables group. Each ant was kept isolated in its own container. My hypothesis was that the ant (3 mg in body weight) would be able to eat 20 times more food than what it weighed; I came to that assumption by relating it to the amount of weight they can lift. The average for the fruits was 357 mg and the average for the vegetables was 358 mg. Keep in mind that I didn't take into account all the weight lost from the food by loss of water evaporating; I just assumed it was all eaten by the ants. For my results I get a t-value of -0.0981 ; df=223 ; standard error of difference=0.015 ; p-value of 0.9220. So since this is my first experiment dealing with t-test. I don't know how to interpret and present the results. Can someone please explain what they mean?


title: Performance of regression models

text: If we have to choose between two regression models, one has a correlation coefficient of 0.95 and the other has a correlation coefficient of 0.75. Is it always the case that the first model is to be preferred over the other if we are interested in optimizing predictive performance?


title: Cell sizes for a fractional factorial design that ignores some main effects

text: We're designing a trial to test three behavioral interventions (call them A, B, and C) for a health condition. Treatment A is an established intervention and B and C are considered add-ons to A. We want to test the following contrasts:   1. Outcomes of packages involving Treatment A vs. No treatment.   2. Outcomes of adding B, C, or B+C to A. Because we're not trying to estimate effects of B or C in the absence of A, we came up with a fractional factorial design implementing the following 5 cells:               Cell    A/B/C    included?     --------------------------     1       + + +    included     2       + + -    included     3       + - +    included     4       + - -    included     5       - + +    not included     6       - + -    not included     7       - - +    not included     8       - - -    included      **How many participants should we have in cell #8?** (the "no treatment" cell) Normally factorial designs put the same number of participants in each cell, but they also have equal numbers of cells with each treatment present vs. absent. Since we have only one cell without treatment A, should this cell be 4x the size of the others? Even if we're comfortable adjusting for unequal cell sizes, should the "no treatment" cell be substantially larger than the other cells?


title: Interpreting significance when regression coefficients are combined

text: Suppose we have the following regression, with Sun=0 for partial sun and Sun=1 for full sun:               Plant Height = a0 + a1*Bacteria + a2*Sun + a3*Bacteria*Sun      And here is the sample output from R: ![R output](http://i.stack.imgur.com/JIPvU.jpg) According to this post, for plants with full sun (Sun=1), the effects of Bacteria is: -0.08346 + 0.03368 = -0.04978, meaning that for two plants with full sun, the plant with one unit more bacteria will be 0.04978cm shorter. So my question is, **is this difference significant?** Since we get the above conclusion by combining the coefficients of Bacteria and Bacteria*Sun, but as in the output, the first coefficient has a significant p-value=0.00214 while the other coefficient for the interaction term has a large p-value=0.34887. So when we add these coefficients, how should we interpret its significance?


title: Likelihood principle: difference between weak and strong version

text: Does anyone understand the difference between weak likelihood principle and strong likelihood principle?


title: Why the division of PAC hypothesis 'gap'

text: I was checking the basic theory of PAC, and I can't understand why the division of the gap between the hypothesis and the class C that we want to learn. Link: > An individual $x$ will be classified incorrectly by the learned concept $h$ > if $x$ lays in the area between $h$ and $c$. We divide this area into 4 > strips, on top, bottom and sides of $h$. We allow these strips, > pessimistically, to overlap in the corners. In figure 1 we represent the top > strip as $t'$. If each of these strips is of area at most epsilon/4, i.e. is > contained in the strip $t$ of area $\frac {\epsilon}{4}$, then the error for > our hypothesis $h$ will be bound by epsilon. Why not do $(1-ε)^m$ instead of the $4(1-(\frac{ε}{4}))^m$ Thanks in advance!


title: Recovering true data from multiple noisy versions

text: I am trying to find if there is any way to get the true data from multiple noisy versions, but the true data has a peculiar property. **Problem Statement** Consider a matrix $F=[f_1, f_2, ... , f_n]$ where $f_i$ is `1984-dimensional`. Now, each $f_i$ can be expressed as $f_i=H^i+n_i$. Here $H^i$ is the true (unknown) data and $n_i$ is the embedded noise. The peculiar property is that each $H^i$ has high correlation with each other (in other words, all $H^i$ are similar to each other). I want to build a filter or a vector which when multiplied with $f_i$, will give me $H^i$ or approximated $H^i$ i.e. noise should reduce. **My thoughts** I understand we should place some extra assumptions on the noise type, like Gaussian etc. In my case, I don't think my noise is Gaussian. Each feature vector $f_i$ basically represents an image patch. So, noise can be random, but I am open to modelling the noise as Gaussian and see what happens. I think a machine learning approach will work better. I had two thoughs:   1. Learn a vector such that when multiplied with the data it should minimize $||w^Tf_i-w^Tf_{i+1}||$ should be minimized such that (say) $||w||\le1$.   2. Add up a subset of $f_i$ so that by adding up different $n_i$, it will approach Gaussian (?) and I assume there are methods to recover data when noise is Gaussian. I am open to any other ideas as long as I get the desired result. You can edit tags if necessary.


title: Improvement on duplicating instances

text: I have a task of Relationship extraction. There are some set of predefined relations in the corpus. I need to train classifier to recognize the type of relation or the lack of relation between every pair of nouns. I chose to use a multi-class implementation of SVM. The problem was that I have different number of examples for different relations, for example, for relation type 1 I have 1000 examples in gold standard annotation, however for relation type 2 only 50 examples. The result of classification was very good result on type 1 and very bad on type 2 on validation set. The solution was to use the same examples from the type 2 several times in my train set, such that the overall number of example of type 1 is close to the number of examples of type 2. It gives pretty good result on test set. The question is how to explain this, why it actually works. For me it looks counter-intuitive, I didn't expect that using the same examples few times could help.


title: How to analyse spatial data where the depending variable is binary

text: I have to test which factors influence game damage in fields. I mapped areas with damage and those without. It was not always possible to map 100% of a field, so there are also areas where it is unsure. Since the "unit" damage is not objective because it is not possible to determine where one ends and the next starts, I put a grid over the area and calculated for every cell independently the distance to different structures (forest, roads etc.). The resulting data look like this:               | damage |  id | dist_forest | dist_maiz | dist_roads |...     |0       |   51|           30|         20|          70|...        |0       |   51|           20|         10|          60|...        |0       |   52|           60|         10|          80|...        |0       |   52|           40|         70|          10|...        |0       |   52|           20|         60|          50|...        |1       |   53|           10|         10|          50|...        |1       |   53|           05|         20|          30|...        |1       |   54|           20|         30|          20|...        |1       |   54|           30|         20|          90|...        |1       |   54|           40|         10|          10|...      (I have about 100 individual polygons which lead to 100000 lines when resolved in square meters) I wanted to use a binary logistic regression with random effects. To solve the problem with the non-independent data I was adding the id of the damage- polygons as a random factor. The resulting model was like this:               glm <- glm(damage ~ dist_forest + dist_maiz + dist_roads + (1|cat), family=binomial(logit),data=data)      The problem ist now, that all my parameters are all highly significant. I already asked that here and it was suggested to use a special model for spatial data. > if data are (substantially) spatially dependent then standard significance > tests don't apply any way! Does anybody have any suggestions how to proceed further?


title: How can I test if a bivariate correlation is significantly different from a partial correlation?

text: I have three variables (A, T, C) with n= 1083. The Pearson correlations are as below: A*T= .170 (p<.01)   A*C= .370 (p<.01)   T*C= -.103 (p<.01) Partial correlation between A and T controlling for C = .225 (p<.01) I would like to know if the increase in A*T (r= .170) is significantly different from partial r (.225). How can I test this?


title: What kind or Result analysis can be done to analyse the characteristics of Restricted boltzmann machine?

text: I have used Mean square error value is one of the parameter to analyse.I also noticed that every time i run a design i get different Mean square error value.Not able to identify the reason . What would be other parameters used to analyse?


title: How to specify mixed ANOVA with multiple repeated measures and covariate in R

text: Working in R, how can I specify a mixed ANOVA with multiple between- and within-subjects factors in such a way that it's amenable to adding a covariate in a subsequent analysis? Also, ideally, I would like to use type 3 SS because that's what I'm used to. I originally did my analysis (without the covariate) using ezANOVA, and found a predicted 3-way interaction of two between-subjects factors (pretrain and training, below) and one within-subjects factor (section, below). A possible explanation is that the two b-s factors affect another quantity (study, below), which in turn causes the above interaction. I am hoping this explanation is NOT correct, but in order to eliminate it, I'd like to add the other quantity ('study') to my model as a covariate and see whether the original 3-way interaction is still significant. Is this an acceptable approach? Assuming yes, my issue is that I can't find a good way to add the covariate. (A) it can be added using "between_covariates" in ezANOVA, but there's a warning that this function is under development, and the results look strange, so I don't trust them without other verification. (B) Using aov, I couldn't exactly replicate my ezANOVA results _without_ the covariate, and also couldn't figure out how to add a covariate. (C) Using lm, I couldn't figure out how to get the repeated measures working. (D) Using lme from nlme, I did get a working model with the repeated measures (shown below), but the results _without_ the covariate were different from those of ezANOVA, so I don't think it's an equivalent model. I'd like to create a model equivalent to the original one, and then see whether my effect stays significant when the covariate is added. How can I do this? Reproducible code:               # This data comes from an experiment designed to investigate effects of different types of instruction on learning. Instruction is manipulated between subjects using factorial combinations of "pretraining" (3 levels) and "training" (2 levels). Learning is assessed as change in performance from a pretest, administered before pretraining & training, to a posttest, administered after. Test accuracy is my dv and test section (pre vs post) is treated as a within-subjects factor. I also have an additional within-subjects factor called "problem type", which is used to classify different types of problems on the pretest and posttest. Some problem types have more problems than others, but the number of problems of a given type is the same for pretest and posttest.          # Sample data:     nsubj = 215; nsec  = 2; nprob = 6     D = data.frame(         subjid   = rep( 1:nsubj, each=nsec*nprob ),         pretrain = rep( sample( c('a','b','c'), nsubj, replace=TRUE ), each=nsec*nprob ),         training = rep( sample( c('j','k'), nsubj, replace=TRUE ), each=nsec*nprob ),         study    = rep( sample( 1:6, nsubj, replace=TRUE ), each=nsec*nprob ),         section  = rep( rep( c( 'pretest', 'posttest' ), each=nprob ), nsubj ),         probtype = rep( c( 'v', 'w', 'x', 'y', 'z', 'z' ), nsec*nsubj ),         accuracy = sample( c( 0.0, 0.5, 1.0 ), nsubj * nsec * nprob, replace=TRUE ) )          # Model:     library( ez )     options( contrasts=c("contr.sum","contr.poly") )     ezANOVA( data=D, wid=.(subjid), dv=.(accuracy), within=.(section,probtype), between=.(pretrain,training), type=3 )          # nlme version:     library( nlme )     fit = lme( accuracy ~ section*probtype*pretrain*training, random=~1|subjid, method='REML', data=D )     anova.lme( fit, type='marginal' )     # results in similar but not identical F and p values to the original model. Denominator dfs are quite different.


title: How to compare task completion times in a crossover study

text: Ten participants have completed a task using two different methods A and B. The first five participants completed the task first using method A then B. The remaining participants completed the same task first using method B then A. The following shows the completion times (minutes) of each method and the order in which the methods were used for each participant.                                Time to complete task using                              (in minutes):     Participant ID      Method A      Method B     Used Method A Before Method B     1                      17            9                       True     2                       8            5                       True     3                      18            8                       True     4                      19           19                       True     5                      37           11                       True     6                      19           18                       False     7                       9            8                       False     8                      12           11                       False     9                      16           18                       False     10                     13           10                       False        1. Is it possible to make a valid statistical argument about the difference between the completion times of the two methods A and B?   2. Is the number of participants large enough to make a valid statistical comparison between the completion time of the two methods?   3. What statistical tests can be applied to compare the completion time of the two methods A and B?   4. Is the difference between the completion times of the two methods A and B statistically significant?


title: Self-organizing maps: fuzzy input?

text: as my first post I would like to know if there are SOM implementations (preferably R) available that accept fuzzy input. That is, I have data in which some nominal features are spread out between a number of categories. For example: feature 1 has 5 categories and an observation might have the values (which are actually probabilities) [0, 0.5, 0.25, 0.25, 0]. Thanks, Figaro


title: References for methods for calculating the confidence interval for Theil-Sen Estimator

text: For the Theil-Sen estimator I am aware of two principal methods for obtaining confidence intervals:   * The zyp R-package documentation uses:  > The confidence interval on the slope is calculated using the method defined > in (Sen, 1968). The confidence interval on the intercept is computed by > taking the standard deviation of the intercepts (sd.i), the z statistic for > the given confidence level (z), and the mean (m). The confidence level is > then: c(m - z * sd.i, m + z * sd.i)   * Wilcox (2009) proposes to use bootstrapping which makes few assumptions but can be computationally costly. Can somebody:   1. Provide a link where the method used by R is explained?   2. Tell me what assumption are made by the method used by R?   3. Link to quantitative research into the consequences/sensibility to data violating the assumptions? Of course I would be equally happy if somebody answered my questions directly rather than providing links.


title: Simple post-stratification weights in R

text: I just got my hands on the ANES (American National Election Studies) 2008 data set, and would like to do some simple analysis in R. However, I've never worked with this complex of a data set before and I've run into an issue. The survey uses oversampling and has a variable for post stratification weights. I had only the vaguest idea of what that meant, so I read the wikipedia page on it, which I understand conceptually. Unfortunately, I don't know how to manipulate R such that the post stratification weights are reflected when I do my analysis. While conceptually, the idea of oversampling didn't confuse me, the following documentation for the R "survey" package is completely unintelligible to me. I'll show what I've found so far, and I would really appreciate either an explanation of what's going on with these methods, or, if anyone knows a simpler way to apply a post-stratification weight to a data frame of variables, I'd love to here that too. So, I found the "survey" package from CRAN, and I have the manual, and, after looking through it, it seems that the most promising method is:               postStratify(design, strata, population, partial = FALSE, ...)      However, when I look at the documentation for what needs to be passed for each of these arguments, I'm completely lost. They are as follows:               design           A survey design with replicate weights          strata           A formula or data frame of post-stratifying variables          population       A table, xtabs or data.frame with population frequencies          partial          if TRUE, ignore population strata not present in the sample      None of these make a lot of sense to me, but I'm pretty sure that the design argument is supposed to be of a class also defined in this package:               svydesign(ids, probs=NULL, strata = NULL, variables = NULL, fpc=NULL,          data = NULL, nest = FALSE,          check.strata = !nest, weights=NULL,pps=FALSE,...)      If you notice, there are a ton of optional arguments here, which all seem to do similar types of things (at least to me, after reading the docs...). I'm basically at a loss for why this is so complicated in R. Am I misunderstanding things? Is there a simpler way to do this? Any help would be appreciated.


title: "weird" results from multilevel analysis

text: I am examining the effect of a group variable (with five levels) along with other predictors on achievement growth using SPSS Mixed Model. The descriptive statistics showed a consistent pattern on achievement for these five groups across all waves of data collection. That is, the average score of group C > Group E> Group B > Group D > Group A. However, the results of my analysis showed that the estimate of Group C on intercept (which was set on the last time point) was significantly lower than Group E and even lower than Group B. Is that possible? I can understand if the difference between Group C and E becomes insignificant because of other predictors in the model. However, is it possible that the addition of other predictors changes the relative standing of the groups in the model? I have rechecked and redone the analysis several times to make sure that the "weird" result was not due to my errors on recoding variables or imputing missing values. Should I be concerned about this results? Or can I say that the result was reasonable since the other predictors I added to the model may "explain away" the group difference? Any comment or advice is appreciated. Thanks!!!


title: Calculate conditional / joint probability distribution

text: I have a beta pert distribution of wall heights, with minimum 0.4m, and maximum of 1m. I also know that along any wall I have either a normal distribution of openings (they are mostly in the middle of the wall) or uniform. Now, I would like to know the probability of an opening being located at any height of any possible wall (cdf). I thought of multiplying two cdfs (the per one and say the uniform one defined between 0 and 1) but obviously it gives me non-zero probabilities when the height of the opening is greater than the wall height which is impossible. Convolution does not seem to be an option as the domain of the openings distribution will vary depending on the wall height (i.e. the pert distribution). How should I solve this? **UPDATE:** I realise what I am trying to do is calculate the probability of the opening height conditional on the wall height. I think the answer lies in the conditional formula for continuous distributions: $$ f_Y(y|X=x) = f(x,y)/f_X(x). $$ with Y being the opening height and X being the wall height. Although the problem is, I only know the wall height probability distribution so I have two unknown distributions...


title: How to sample from $\{1, 2, ..., K\}$ for $n$ random variables, each with different mass functions, in R?

text: In R, I have an $N \times K$ matrix $P$ where the $i$'th row of $P$ corresponds to a distribution on $\\{1, ..., K\\}$. Essentially, I need to sample from each row efficiently. A naive implementation is:               X = rep(0, N);     for(i in 1:N){         X[i] = sample(1:K, 1, prob = P[i, ]);     }      This is much too slow. In principle I could move this to C but I'm sure there must be an existing way of doing this. I would like something in the spirit of the following code (which does not work):               X = sample(1:K, N, replace = TRUE, prob = P)      **EDIT:** For motivation, take $N = 10000$ and $K = 100$. I have $P_1, ..., P_{5000}$ matrices all $N \times K$ and I need to sample a vector from each of them.


title: Cross-validation and R

text: Consider the following data set **train** :               z  a       1  1       0  2       0  1     1  3     0  1     1  2     1  1     0  3     0  1     1  3      with binary outcome variable **z** and predictor variable **a**. Also consider the following code:               library(caret)     set.seed(825)     fitControl <- trainControl(method= "cv", number=2, classProbs=TRUE)     svmfit1 <- train(z~a, data=train, method="svmLinear", trControl = fitControl)      What exactly is the model ` svmfit1 `? It is not using all the data to fit the model. Each time, we are leaving out 2 data points. Thus we get 5 SVM models. But which one is ` svmfit1 `? Does R just choose ` svmfit1 ` to be the model with the highest accuracy out of the five?


title: Confidence interval for a constrained fit to Gaussian-like data

text: I'm working with data from an instrument which is expected a priori to produce Gaussian (normally) distributed data: \begin{equation} G = A\exp\left(-\dfrac{(x - \mu)^2}{\sigma} \right) \end{equation} The data are normally sparse, with only about 2-3 measurements representing each $G$. In this question, I'm focussing on a single $G$, but in reality, we often have analyte signals causing overlapping $G$ that are then fitted simultaneously as described below. To fit $G$ to our measurements, we calibrate for $\mu$ and $\sigma$ a priori using reference signals, then use these calibrations to constrain all parameters except $A$. So the fit reduces to \begin{equation} G = AG_0 \end{equation} which is then fitted by least-squares minimization to determine $A$. **My main question is: what is the uncertainty ($\sigma_A$) in the fitted $A$?** * * * My initial approach was to estimate $\sigma_A$ as the RMSE of the fit. But since the RMSE is essentially the standard deviation of the residuals, this seems like an overestimate: I want the confidence interval of the fitted parameter. Can I safely use the following textbook equation for the confidence interval of a linear-regression slope in this context? (the $x_i$ are the predictor variables, $\bar{x}$ is the mean of the predictor variables, the $e_i$ are fit residuals, and 2 degrees of freedom were consumed by the 1 fitted variable and the fact that the residuals sum to zero) \begin{equation} \sigma_A = \dfrac{s_e}{\Sigma _i (x_i - \bar{x})^2 } = \dfrac{\frac{1}{n- 2}\Sigma_i e_i^2}{\Sigma _i (x_i - \bar{x})^2 } \end{equation} I think "yes": all I have done is transform my $x$ before fitting a linear parameter. I also think "no" because I'm not sure of the meaning of $x_i - \bar{x}$ in this equation -- is it specific to linear regression? * * * To give you a visualization, my data are almost as bad as these simulated data (but often the signal ontains one or two data points more): ![A gaussian fit to only 2 data](http://i.stack.imgur.com/edBWL.png) Note: I am aware that Bayesian analysis would be a better method for passing information about $\mu$ and $\sigma$ to my fit, but I am not at liberty to change the analysis software right now. I need to limit myself to an estimate of $\sigma_A$. edit: Another note that might exclude some solutions: I am analyzing thousands of measurements in bulk with no known true value.


title: Logistic regression & categorical explanatory variables

text: I'd like to perform logistic regression with some categorical explanatory variables with more categories than just binary 0/1. Is this possible and why? I am inclined to think that this would give just the wrong result because of the geometric intution, however most opinions online say that explanatory variables can be discrete/categorical (although I don't see any mentions of more categories). For 0/1 this is fine because the distances between 0 and 1 is easy. If I have even 3 categories 0/1/2, it's unclear if 1-0 is the same as 2-1. One thought I had is to do a sort of one-vs-all thing for the explanatory variables. So if I had one parameter that can on 3 categories A/B/C I remake this into 3 separate parameters: A: 0/1 B: 0/1 C: 0/1 Where 0 is not belonging to the letter class and 1 is belonging to the letter class. Background: My dependent variable should be binary. I really wanted to use logistic regression because of the probability interpretation. If you can recommend another method that can output a probability estimate (instead of only classification) that takes in categorical explanatory variables, that would also be helpful.


title: slope of non linear noisy data

text: I have data from{ Erne /Soho} for time in unit( hour )on x-axis and the intensity in unit ( 1/s.sr.cm2.Mev/n) on y axis the data so long from 1500 points to 2500 points and the curves look like Gaussian curve i need your help about :   1. Plot and find the fitting , the polyvalue and polynomial   2. Find the slope ( angle ) for the rising points from the point that beginning to rise to the peak of the curve   3. Repeat step 2 for points from peak to the end decay the data is:               x(time in hrs)  y(intensity)            X    y     0.041667    0.000171     0.125   0.000161     0.208333    0.000101     0.291667    0.000139     0.375   0.0000565     0.458333    0.000124     0.541667    0.0000669     0.625   0.0000908     0.708333    0.0000869     0.791667    0.000164     0.875   0.000137     0.958333    0.0000885     1.041667    0.000162     1.125   0.000141     1.208333    0.000125     1.291667    0.000104     1.375   0.000067     1.458333    0.0000658     1.541667    0.000135     1.625   0.0000691     1.708333    0.0000571     1.791667    0.0000895     1.875   0.000149     1.958333    0.000118     2.041667    0.000162     2.125   0.00016     2.208333    0.0000777     2.291667    0.0000882     2.375   0.000156     2.458333    0.0000706     2.541667    0.000136     2.625   0.000102     2.708333    0.000143     2.791667    0.0000921     2.875   0.0000803     2.958333    0.000106     3.041667    0.0000803     3.125   0.000117     3.208333    0.000161     3.291667    0.000183     3.375   0.0000588     3.458333    0.0000687     3.541667    0.000138     3.625   0.000101     3.708333    0.000126     3.791667    0.0000352     3.875   0.0000344     3.958333    0.0000684     4.041667    0.000102     4.125   0.000149     4.208333    0.000147     4.291667    0.0000222     4.375   0.000034     4.458333    0.0000929     4.541667    0.000092     4.625   0.00008     4.708333    0.000101     4.791667    0.000145     4.875   0.0000794     4.958333    0.000102     5.041667    0.000283     5.125   0.000341     5.208333    0.000711     5.291667    0.000892     5.375   0.001269     5.458333    0.001324     5.541667    0.001479     5.625   0.002049     5.708333    0.001915     5.791667    0.002428     5.875   0.002894     5.958333    0.003405     6.041667    0.003135     6.125   0.003133     6.208333    0.003374     6.291667    0.003532     6.375   0.003844     6.458333    0.003582     6.541667    0.004157     6.625   0.004822     6.708333    0.005595     6.791667    0.006492     6.875   0.00688     6.958333    0.007222     7.041667    0.007582     7.125   0.007961     7.208333    0.007937     7.291667    0.007869     7.375   0.007804     7.458333    0.007741     7.541667    0.007369     7.625   0.006984     7.708333    0.00662     7.791667    0.006275     7.875   0.006338     7.958333    0.006447     8.041667    0.006536     8.125   0.006673     8.208333    0.006869     8.291667    0.007084     8.375   0.007309     8.458333    0.007546     8.541667    0.007243     8.625   0.006907     8.708333    0.006596     8.791667    0.00631     8.875   0.006351     8.958333    0.006428     9.041667    0.006509     9.125   0.006592     9.208333    0.006545     9.291667    0.006488     9.375   0.006433     9.458333    0.006382     9.541667    0.006437     9.625   0.006511     9.708333    0.006594     9.791667    0.006685     9.875   0.006784     9.958333    0.006889     10.04167    0.006998     10.125  0.007113     10.20833    0.007222     10.29167    0.007335     10.375  0.007453     10.45833    0.007577     10.54167    0.007719     10.625  0.007867     10.70833    0.00802     10.79167    0.008178     10.875  0.008258     10.95833    0.008332     11.04167    0.008407     11.125  0.008483     11.20833    0.008638     11.29167    0.008805     11.375  0.008976     11.45833    0.009151     11.54167    0.008961     11.625  0.00874     11.70833    0.008531     11.79167    0.008333     11.875  0.008155     11.95833    0.007988     12.04167    0.007829     12.125  0.00768     12.20833    0.00746     12.29167    0.007242     12.375  0.007034     12.45833    0.006835     12.54167    0.006705     12.625  0.006584     12.70833    0.006466     12.79167    0.006351     12.875  0.006259     12.95833    0.006173     13.04167    0.006093     13.125  0.006017     13.20833    0.005997     13.29167    0.005985     13.375  0.005975     13.45833    0.005969     13.54167    0.005572     13.625  0.005167     13.70833    0.0048     13.79167    0.004464     13.875  0.004469     13.95833    0.004521     14.04167    0.004583     14.125  0.004694     14.20833    0.00436     14.29167    0.004653     14.375  0.004479     14.45833    0.0044     14.54167    0.004331     14.625  0.004178     14.70833    0.004028     14.79167    0.004093     14.875  0.003896     14.95833    0.003949     15.04167    0.003895     15.125  0.003843     15.20833    0.003795     15.29167    0.004257     15.375  0.00382     15.45833    0.003776     15.54167    0.003662     15.625  0.003628     15.70833    0.00356     15.79167    0.00358     15.875  0.003705     15.95833    0.00354     16.04167    0.002857     16.125  0.00322     16.20833    0.003162     16.29167    0.003105     16.375  0.003049     16.45833    0.003023     16.54167    0.002811     16.625  0.003027     16.70833    0.002789     16.79167    0.002864     16.875  0.002708     16.95833    0.002964     17.04167    0.003224     17.125  0.002653     17.20833    0.002214     17.29167    0.002505     17.375  0.002372     17.45833    0.00202     17.54167    0.002151     17.625  0.00267     17.70833    0.002533     17.79167    0.002248     17.875  0.002042     17.95833    0.002338     18.04167    0.002448     18.125  0.002201     18.20833    0.001753     18.29167    0.001924     18.375  0.001724     18.45833    0.0019     18.54167    0.00188     18.625  0.001866     18.70833    0.001861     18.79167    0.001848     18.875  0.001671     18.95833    0.002026     19.04167    0.001911     19.125  0.001747     19.20833    0.002013     19.29167    0.001459     19.375  0.001561     19.45833    0.001312     19.54167    0.001475     19.625  0.00186     19.70833    0.001447     19.79167    0.001411     19.875  0.001464     19.95833    0.001573     20.04167    0.001883     20.125  0.001447     20.20833    0.001521     20.29167    0.001419     20.375  0.00148     20.45833    0.001458     20.54167    0.001428     20.625  0.001622     20.70833    0.001325     20.79167    0.001302     20.875  0.001308     20.95833    0.001322     21.04167    0.001437     21.125  0.001255     21.20833    0.0012     21.29167    0.001379     21.375  0.001139     21.45833    0.001016     21.54167    0.001307     21.625  0.000904     21.70833    0.00098     21.79167    0.000771     21.875  0.000852     21.95833    0.001013     22.04167    0.001236     22.125  0.00102     22.20833    0.000923     22.29167    0.001115     22.375  0.001144     22.45833    0.001097     22.54167    0.000909     22.625  0.001157     22.70833    0.00095     22.79167    0.001014     22.875  0.001137     22.95833    0.00109     23.04167    0.001046     23.125  0.001032     23.20833    0.000918     23.29167    0.000785     23.375  0.000938     23.45833    0.00086     23.54167    0.000814     23.625  0.000818     23.70833    0.000865     23.79167    0.000951     23.875  0.001015     23.95833    0.000808     24.04167    0.000782     24.125  0.000609     24.20833    0.000741     24.29167    0.000756     24.375  0.00084     24.45833    0.000786     24.54167    0.000659     24.625  0.000725     24.70833    0.000603     24.79167    0.000748     24.875  0.000664     24.95833    0.000719     25.04167    0.000648     25.125  0.000688     25.20833    0.000528     25.29167    0.000644     25.375  0.000612     25.45833    0.000538     25.54167    0.000597     25.625  0.000653     25.70833    0.000598     25.79167    0.000632     25.875  0.000594     25.95833    0.000501     26.04167    0.000652     26.125  0.000739     26.20833    0.000511     26.29167    0.000611     26.375  0.00053     26.45833    0.000646     26.54167    0.000533     26.625  0.00063     26.70833    0.000551     26.79167    0.000406     26.875  0.000699     26.95833    0.000631     27.04167    0.000674     27.125  0.000469     27.20833    0.000565     27.29167    0.000397     27.375  0.000589     27.45833    0.000464     27.54167    0.000573     27.625  0.000482     27.70833    0.000569     27.79167    0.000526     27.875  0.000482     27.95833    0.000458     28.04167    0.000569     28.125  0.000515     28.20833    0.000638     28.29167    0.000516     28.375  0.000525     28.45833    0.000605     28.54167    0.000469     28.625  0.000492     28.70833    0.000262     28.79167    0.000347     28.875  0.000434     28.95833    0.000399     29.04167    0.000448     29.125  0.000473     29.20833    0.000386     29.29167    0.00043     29.375  0.000403     29.45833    0.000465     29.54167    0.000441     29.625  0.000375     29.70833    0.000403     29.79167    0.000374     29.875  0.000328     29.95833    0.000459     30.04167    0.000365     30.125  0.000403     30.20833    0.00037     30.29167    0.000587     30.375  0.000441     30.45833    0.000316     30.54167    0.000379     30.625  0.000477     30.70833    0.000425     30.79167    0.000488     30.875  0.000358     30.95833    0.000369     31.04167    0.000392     31.125  0.000399     31.20833    0.000326     31.29167    0.0004     31.375  0.000283     31.45833    0.00027     31.54167    0.000323     31.625  0.000356     31.70833    0.000415     31.79167    0.000388     31.875  0.000259     31.95833    0.000488     32.04167    0.000304     32.125  0.00041     32.20833    0.000395     32.29167    0.000344     32.375  0.000382     32.45833    0.000357     32.54167    0.000431     32.625  0.000289     32.70833    0.000346     32.79167    0.000302     32.875  0.000237     32.95833    0.000366     33.04167    0.000272     33.125  0.000205     33.20833    0.000322     33.29167    0.000299     33.375  0.000315     33.45833    0.000277     33.54167    0.000204     33.625  0.000421     33.70833    0.000304     33.79167    0.000281     33.875  0.000289     33.95833    0.000214     34.04167    0.000281     34.125  0.000266     34.20833    0.000268     34.29167    0.000267     34.375  0.000288     34.45833    0.000202     34.54167    0.000309     34.625  0.000299     34.70833    0.000393     34.79167    0.000211     34.875  0.000286     34.95833    0.000317     35.04167    0.000224     35.125  0.000233     35.20833    0.000298     35.29167    0.000234     35.375  0.000202     35.45833    0.000255     35.54167    0.000287     35.625  0.000274     35.70833    0.000319     35.79167    0.000254     35.875  0.000297     35.95833    0.000382     36.04167    0.0003     36.125  0.000263     36.20833    0.000201     36.29167    0.000243     36.375  0.00018     36.45833    0.00019     36.54167    0.000284     36.625  0.000232     36.70833    0.000224     36.79167    0.000199     36.875  0.000317     36.95833    0.00021     37.04167    0.000233     37.125  0.00019     37.20833    0.000137     37.29167    0.000127     37.375  0.000254     37.45833    0.000272     37.54167    0.000335     37.625  0.000221     37.70833    0.000234     37.79167    0.000178     37.875  0.000283     37.95833    0.000282     38.04167    0.000179     38.125  0.000136     38.20833    0.000137     38.29167    0.000147     38.375  0.000178     38.45833    0.000242     38.54167    0.000136     38.625  0.000136     38.70833    0.000316     38.79167    0.000284     38.875  0.000199     38.95833    0.000221     39.04167    0.000168     39.125  0.00021     39.20833    0.000208     39.29167    0.000271     39.375  0.000199     39.45833    0.000199     39.54167    0.000251     39.625  0.000188     39.70833    0.000231     39.79167    0.000147     39.875  0.000105     39.95833    0.000197     40.04167    0.0000835     40.125  0.00025     40.20833    0.000229     40.29167    0.000209     40.375  0.000158     40.45833    0.000147     40.54167    0.000115     40.625  0.000187     40.70833    0.000221     40.79167    0.00023     40.875  0.000229     40.95833    0.000157     41.04167    0.000241     41.125  0.000208     41.20833    0.000178     41.29167    0.000187     41.375  0.000188     41.45833    0.000187     41.54167    0.000188     41.625  0.000105     41.70833    0.000199     41.79167    0.000261     41.875  0.000155     41.95833    0.000166     42.04167    0.000146     42.125  0.000209     42.20833    0.000177     42.29167    0.000199     42.375  0.000137     42.45833    0.000168     42.54167    0.000126     42.625  0.000115     42.70833    0.0000847     42.79167    0.000167     42.875  0.000189     42.95833    0.000177     43.04167    0.0000947     43.125  0.000156     43.20833    0.000178     43.29167    0.000105     43.375  0.000193     43.45833    0.000157     43.54167    0.000136     43.625  0.000115     43.70833    0.000179     43.79167    0.000148     43.875  0.000114     43.95833    0.000127     44.04167    0.000179     44.125  0.000211     44.20833    0.00022     44.29167    0.000148     44.375  0.000126     44.45833    0.000106     44.54167    0.000096     44.625  0.000171     44.70833    0.0000954     44.79167    0.000157     44.875  0.000213     44.95833    0.000107     45.04167    0.00017     45.125  0.000136     45.20833    0.000168     45.29167    0.000162     45.375  0.000117     45.45833    0.000148     45.54167    0.000178     45.625  0.000116     45.70833    0.000136     45.79167    0.000159     45.875  0.000116     45.95833    0.000169     46.04167    0.000191     46.125  0.000104     46.20833    0.000117     46.29167    0.00016     46.375  0.000201     46.45833    0.000128     46.54167    0.000127     46.625  0.00016     46.70833    0.000139     46.79167    0.0000964     46.875  0.000126     46.95833    0.000137     47.04167    0.000159     47.125  0.000138     47.20833    0.0000862     47.29167    0.000129     47.375  0.000108     47.45833    0.0000846     47.54167    0.000171     47.625  0.00016     47.70833    0.000149     47.79167    0.00013     47.875  0.000106     47.95833    0.000181     48.04167    0.000148     48.125  0.000139     48.20833    0.000128     48.29167    0.000149     48.375  0.0000754     48.45833    0.000116     48.54167    0.000194     48.625  0.000158     48.70833    0.000106     48.79167    0.000179     48.875  0.000139     48.95833    0.000107     49.04167    0.00018     49.125  0.000116     49.20833    0.000106     49.29167    0.000161     49.375  0.000106     49.45833    0.00014     49.54167    0.0000848     49.625  0.000116     49.70833    0.000141     49.79167    0.000173     49.875  0.000117     49.95833    0.000128     50.04167    0.0000848     50.125  0.000182     50.20833    0.00014     50.29167    0.000105     50.375  0.0000863     50.45833    0.000174     50.54167    0.000213     50.625  0.000159     50.70833    0.000118     50.79167    0.000138     50.875  0.000139     50.95833    0.000191     51.04167    0.000117     51.125  0.000181     51.20833    0.000169     51.29167    0.000148     51.375  0.000138     51.45833    0.00017     51.54167    0.000203     51.625  0.000117     51.70833    0.000126     51.79167    0.00016     51.875  0.000076     51.95833    0.0000845     52.04167    0.000128     52.125  0.00015     52.20833    0.000127     52.29167    0.0000955     52.375  0.000224     52.45833    0.00013     52.54167    0.0000851     52.625  0.00012     52.70833    0.000108     52.79167    0.000116     52.875  0.000129     52.95833    0.000149     53.04167    0.0000954     53.125  0.000128     53.20833    0.000097     53.29167    0.000183     53.375  0.0000856     53.45833    0.000173     53.54167    0.00014     53.625  0.000128     53.70833    0.0000852     53.79167    0.000128     53.875  0.0000954     53.95833    0.000129     54.04167    0.0000651     54.125  0.000117     54.20833    0.0000867     54.29167    0.000164     54.375  0.0000859     54.45833    0.000194     54.54167    0.000142     54.625  0.0000988     54.70833    0.000118     54.79167    0.000176     54.875  0.000141     54.95833    0.000119     55.04167    0.000106     55.125  0.0000856     55.20833    0.000107     55.29167    0.0000965     55.375  0.000129     55.45833    0.000141     55.54167    0.000176     55.625  0.000108     55.70833    0.0000647     55.79167    0.000108     55.875  0.000143     55.95833    0.0000765     56.04167    0.000108     56.125  0.0000656     56.20833    0.0000872     56.29167    0.000129     56.375  0.000101     56.45833    0.0000883     56.54167    0.000139     56.625  0.000142     56.70833    0.000128     56.79167    0.000199     56.875  0.0000539     56.95833    0.0000549     57.04167    0.000161     57.125  0.000119     57.20833    0.000137     57.29167    0.000121     57.375  0.000177     57.45833    0.000134     57.54167    0.000153     57.625  0.000176     57.70833    0.000166     57.79167    0.00014     57.875  0.000121     57.95833    0.000108     58.04167    0.0000757     58.125  0.000164     58.20833    0.000192     58.29167    0.0000867     58.375  0.000173     58.45833    0.000129     58.54167    0.0000863     58.625  0.00015     58.70833    0.000098     58.79167    0.0000767     58.875  0.00011     58.95833    0.000111     59.04167    0.000108     59.125  0.000173     59.20833    0.000117     59.29167    0.0000442     59.375  0.00011     59.45833    0.00012     59.54167    0.00014     59.625  0.000076     59.70833    0.000195     59.79167    0.000165     59.875  0.000163     59.95833    0.000163     60.04167    0.0001     60.125  0.0000996     60.20833    0.00011     60.29167    0.0000649     60.375  0.0000764     60.45833    0.000172     60.54167    0.0000871     60.625  0.000121     60.70833    0.000207     60.79167    0.000128     60.875  0.000109     60.95833    0.000179     61.04167    0.000144     61.125  0.00012     61.20833    0.0000646     61.29167    0.0000764     61.375  0.00013     61.45833    0.00017     61.54167    0.000184     61.625  0.000153     61.70833    0.000169     61.79167    0.0000546     61.875  0.000142     61.95833    0.0000884     62.04167    0.000113     62.125  0.000153     62.20833    0.000179     62.29167    0.0000672     62.375  0.0000989     62.45833    0.00011     62.54167    0.000175     62.625  0.0000659     62.70833    0.000121     62.79167    0.000164     62.875  0.000141     62.95833    0.000146     63.04167    0.000139     63.125  0.000155     63.20833    0.000145     63.29167    0.000111     63.375  0.000165     63.45833    0.000111     63.54167    0.000133     63.625  0.0000764     63.70833    0.0000883     63.79167    0.000145     63.875  0.000124     63.95833    0.000214     64.04167    0.000089     64.125  0.000122     64.20833    0.000144     64.29167    0.000165     64.375  0.00012     64.45833    0.0000664     64.54167    0.000165     64.625  0.0000668     64.70833    0.000111     64.79167    0.000135     64.875  0.000111     64.95833    0.000122     65.04167    0.000134     65.125  0.0000989     65.20833    0.000172     65.29167    0.000167     65.375  0.000145     65.45833    0.000123     65.54167    0.000195     65.625  0.0000881     65.70833    0.0000893     65.79167    0.000146     65.875  0.000123     65.95833    0.000174     66.04167    0.000119     66.125  0.000122     66.20833    0.000151     66.29167    0.000177     66.375  0.0001     66.45833    0.0000795     66.54167    0.0000768     66.625  0.000124     66.70833    0.000176     66.79167    0.000168     66.875  0.000133     66.95833    0.0000772     67.04167    0.000112     67.125  0.0000895     67.20833    0.0000801     67.29167    0.000125     67.375  0.000112     67.45833    0.000144     67.54167    0.0000975     67.625  0.000101     67.70833    0.0000664     67.79167    0.000191     67.875  0.000119     67.95833    0.000168     68.04167    0.000213     68.125  0.0000803     68.20833    0.00013     68.29167    0.0000996     68.375  0.000101     68.45833    0.000122     68.54167    0.000068     68.625  0.000144     68.70833    0.0000908     68.79167    0.000157     68.875  0.00011     68.95833    0.0000796     69.04167    0.000115     69.125  0.000178     69.20833    0.0000911     69.29167    0.00009     69.375  0.000149     69.45833    0.000121     69.54167    0.000135     69.625  0.000113     69.70833    0.000103     69.79167    0.000145     69.875  0.0000814     69.95833    0.000111     70.04167    0.000113     70.125  0.000115     70.20833    0.000133     70.29167    0.0000561     70.375  0.0000904     70.45833    0.000124     70.54167    0.0000662     70.625  0.000101     70.70833    0.0000769     70.79167    0.0000889     70.875  0.000125     70.95833    0.000125     71.04167    0.000147     71.125  0.000132     71.20833    0.0000882     71.29167    0.000159     71.375  0.000137     71.45833    0.000135     71.54167    0.0000879     71.625  0.0000972     71.70833    0.0000781     71.79167    0.000157     71.875  0.000103     71.95833    0.000155     72.04167    0.000155     72.125  0.0000677     72.20833    0.000101     72.29167    0.000135     72.375  0.000114     72.45833    0.000211     72.54167    0.000114     72.625  0.000195     72.70833    0.000225     72.79167    0.000131     72.875  0.000121     72.95833    0.000143     73.04167    0.000101     73.125  0.0000552     73.20833    0.000101     73.29167    0.000123     73.375  0.000134     73.45833    0.0000788     73.54167    0.00011     73.625  0.0000572     73.70833    0.0000805     73.79167    0.0000909     73.875  0.000203     73.95833    0.000144     74.04167    0.000156     74.125  0.000144     74.20833    0.0000785     74.29167    0.00016     74.375  0.000101     74.45833    0.000145     74.54167    0.000155     74.625  0.0000895     74.70833    0.000125     74.79167    0.000167     74.875  0.0000559     74.95833    0.000128     75.04167    0.000101     75.125  0.0000892     75.20833    0.0000912     75.29167    0.000145     75.375  0.000112     75.45833    0.0000904     75.54167    0.000167     75.625  0.000137     75.70833    0.000103     75.79167    0.000172     75.875  0.0000915     75.95833    0.000103     76.04167    0.0000977     76.125  0.000101     76.20833    0.0000871     76.29167    0.0001     76.375  0.00008     76.45833    0.000069     76.54167    0.000116     76.625  0.0000796     76.70833    0.0000669     76.79167    0.000122     76.875  0.000145     76.95833    0.0000993     77.04167    0.0000785     77.125  0.000102     77.20833    0.000113     77.29167    0.000099     77.375  0.000141     77.45833    0.0000928     77.54167    0.000114     77.625  0.000194     77.70833    0.000219     77.79167    0.000223     77.875  0.0000917     77.95833    0.000124     78.04167    0.0000465     78.125  0.000182     78.20833    0.000147     78.29167    0.000091     78.375  0.00017     78.45833    0.000124     78.54167    0.000112     78.625  0.0000696     78.70833    0.0000942     78.79167    0.0000947     78.875  0.000136     78.95833    0.000161     79.04167    0.000138     79.125  0.000185     79.20833    0.0000931     79.29167    0.0000659     79.375  0.000109     79.45833    0.000079     79.54167    0.000101     79.625  0.000115     79.70833    0.000184     79.79167    0.000123     79.875  0.0000787     79.95833    0.000127     80.04167    0.000122     80.125  0.0000794     80.20833    0.0000657     80.29167    0.0000967     80.375  0.000173     80.45833    0.0000923     80.54167    0.0000912     80.625  0.000226     80.70833    0.000135     80.79167    0.000147     80.875  0.000123     80.95833    0.000159     81.04167    0.0000683     81.125  0.000172     81.20833    0.000161     81.29167    0.00013     81.375  0.000147     81.45833    0.000113     81.54167    0.000149     81.625  0.0000423     81.70833    0.000128     81.79167    0.000103     81.875  0.000139     81.95833    0.000191     82.04167    0.000103     82.125  0.000115     82.20833    0.0000915     82.29167    0.000147     82.375  0.0000866     82.45833    0.0000928     82.54167    0.000114     82.625  0.0000978     82.70833    0.0000872     82.79167    0.000113     82.875  0.000181     82.95833    0.000102     83.04167    0.000103     83.125  0.000154     83.20833    0.000124     83.29167    0.000154     83.375  0.000162     83.45833    0.000112     83.54167    0.000105     83.625  0.0000771     83.70833    0.000113     83.79167    0.000213     83.875  0.000185     83.95833    0.0000373     84.04167    0.00024     84.125  0.0000918     84.20833    0.0000912     84.29167    0.0000997     84.375  0.0000915     84.45833    0.0000695     84.54167    0.0000896     84.625  0.000113     84.70833    0.000113     84.79167    0.000127     84.875  0.000103     84.95833    0.000139     85.04167    0.000105     85.125  0.00021     85.20833    0.000146     85.29167    0.000198     85.375  0.000106     85.45833    0.000132     85.54167    0.0001     85.625  0.0000683     85.70833    0.0000872     85.79167    0.0000555     85.875  0.000125     85.95833    0.0000909     86.04167    0.000102     86.125  0.0000682     86.20833    0.000112     86.29167    0.0000699     86.375  0.00012     86.45833    0.0000919     86.54167    0.000142     86.625  0.000181     86.70833    0.0000685     86.79167    0.0000915     86.875  0.000112     86.95833    0.0000799     87.04167    0.0000679     87.125  0.0000343     87.20833    0.000123     87.29167    0.0000813     87.375  0.0000914     87.45833    0.000198     87.54167    0.000162     87.625  0.000205     87.70833    0.000106     87.79167    0.000153     87.875  0.00017     87.95833    0.0000789     88.04167    0.000183     88.125  0.0000906     88.20833    0.000149     88.29167    0.000106     88.375  0.0000828     88.45833    0.000206     88.54167    0.000134     88.625  0.000058     88.70833    0.000122     88.79167    0.000101     88.875  0.000133     88.95833    0.000101     89.04167    0.000141     89.125  0.000082     89.20833    0.000149     89.29167    0.000106     89.375  0.0000949     89.45833    0.000127     89.54167    0.0000907     89.625  0.00018     89.70833    0.000068     89.79167    0.000119     89.875  0.0000567     89.95833    0.000104     90.04167    0.000104     90.125  0.0000577     90.20833    0.000151     90.29167    0.000114     90.375  0.000106     90.45833    0.0000959     90.54167    0.000149     90.625  0.000171     90.70833    0.000102     90.79167    0.000128     90.875  0.000112     90.95833    0.0000881     91.04167    0.00016     91.125  0.000103     91.20833    0.000105     91.29167    0.000115     91.375  0.000118     91.45833    0.0000929     91.54167    0.000181     91.625  0.000101     91.70833    0.000101     91.79167    0.000117     91.875  0.000149     91.95833    0.0000908     92.04167    0.0000987     92.125  0.0000569     92.20833    0.000133     92.29167    0.000103     92.375  0.000106     92.45833    0.000103     92.54167    0.000115     92.625  0.0000816     92.70833    0.000137     92.79167    0.000146     92.875  0.000138     92.95833    0.000142     93.04167    0.000169     93.125  0.0000819     93.20833    0.0000943     93.29167    0.0000455     93.375  0.000114     93.45833    0.0000929     93.54167    0.000186     93.625  0.00024     93.70833    0.000114     93.79167    0.000136     93.875  0.0000706     93.95833    0.000129     94.04167    0.000115     94.125  0.000152     94.20833    0.000148     94.29167    0.000101     94.375  0.000158     94.45833    0.000159     94.54167    0.000116     94.625  0.000127     94.70833    0.000112     94.79167    0.000153     94.875  0.000207     94.95833    0.0000804     95.04167    0.000143     95.125  0.00015     95.20833    0.000105     95.29167    0.000117     95.375  0.000127     95.45833    0.000105     95.54167    0.000114     95.625  0.00014     95.70833    0.000131     95.79167    0.0000934     95.875  0.0000589     95.95833    0.0000782      Most of my data is noisy and have large fluctuation my problem is:   1. i can’t find the fitting for the curve that pass through the onset time ( the point of curve that begin to increase until the peak point and then find the angle between the line and x-axes , and find the angle of the line that pass through the peak point and the point of end decay on the curve)   2. i asked for best way in excel or matlab to write a program to do so.   3. i asked if possible separate the data to two parts : rising part and decay part and then find the slop of each line. I tried only to plot the data put i have no idea about programming to do so and what is the equation of the curve or the lines and errors.


title: How to compute an exponentially weighted covariance matrix function in R?

text: Is there a function the calculates an exponentially weighted covariance matrix in R? I recall seeing a function some time ago... I cannot turn up any thing in "R Search" or Task View pages.


title: Is it a Wishart matrix?

text: We know that an $m \times m$ random matrix $\boldsymbol{A} = \boldsymbol{H} \boldsymbol{H}^H$ is a (central) real/complex Wishart matrix with $n$ degrees of freedom and covariance matrix $\boldsymbol{\Sigma}$ $(\boldsymbol{A} \sim \mathcal{W}_m (n, \boldsymbol{\Sigma}))$, if the columns of the $m \times n$ matrix $\boldsymbol{H}$ are zero mean independent real/complex Gaussian vectors with covariance matrix $\boldsymbol{\Sigma}$. If the columns of matrix $\boldsymbol{H}$ are zero mean independent real/complex Gaussian vectors with different covariance matrix, i.e, the $i$-th column of $\boldsymbol{H}$ has covariance matrix $\boldsymbol{\Sigma}_i$ ($i=1, \cdots, n$), is matrix $\boldsymbol{A}$ a Wishart matrix, and what is its distribution? Thanks!


title: Why the infrequent use of machine learning techniques in translational biomedicine?

text: This is in followup to a previous question. here : Neural network model to predict treatment outcome and might be considered to refer to a different aspect of this question: Application of machine learning techniques in small sample clinical studies Thanks to Zach who suggested reposting. I've put some fairly serious reading in now on CART, randomForest, Neural Networks and machine learning in general, learned about WEKA and the R packages, seen and followed the Stanford engineering lectures http://www.ml- class.org/course/class/index , I'm 3 chapters into Hastie. Given the kind of data we see regularly in clinically oriented research - loads of clinical parameters + loads of biochemical parameters + pen and paper test data +/- neuroimaging data with smallish numbers, I get the feeling that I'm missing something. I'm not regularly reading about ML techniques being applied in the research literature. My question is: have I just latched on to something which is dubious and therefore regarded with justified suspicion by researching clinicians and biostatisticians who are well aware of it, or are these techniques genuinely overlooked or feared outside of "business analytics"? What keeps it "niche"?


title: Can I trust a regression if variables are autocorrelated?

text: Both variables (dependent and independent) show autocorrelation effects. Data is time-series and stationary When I run the regression residuals appear not to be correlated. My Durbin- Watson statistic is greater than upper critical value, so there is an evidence that error terms are not positively correlated. Also when I plot ACF for errors it looks like there is no correlation there and Ljung-Box statistic is smaller than critical value. Can I trust my regression output, are the t-statistics reliable?


title: What kind of statistical method do I need for comparing element concentrations?

text: I have assayed the concentrations of various metal elements in different anatomic regions of two strains of mice. Now, for each element, in each region, I want to find whether there is any difference between the two strains. The concentrations of the different metals within an animal are correlated, it is said that doing pairwise t tests will almost certainly result in nonsense P values. What kind of statistical method do I need here? Preferably, please suggest in the context of the R statistical language. Here is what I did, perhaps incorrectly (using simulated data as an example, in reality I have 7 mice in each group, 6 anatomic regions for each mouse, 13 elements for each region):               # create the data frame     > elemconc = data.frame(expand.grid(id=1:3, geno=c('exp', 'wt'), region=c('brain', 'spine'), elem=c('fe', 'cu', 'zn')), conc=rnorm(36, 10))     > elemconc       id geno region elem      conc     1   1  exp  brain   fe  8.497498     2   2  exp  brain   fe  9.280944     3   3  exp  brain   fe  9.726271     4   1   wt  brain   fe 11.556397     5   2   wt  brain   fe 10.992550     6   3   wt  brain   fe  9.711200     7   1  exp  spine   fe 11.168603     8   2  exp  spine   fe  9.331127     9   3  exp  spine   fe 11.048226     10  1   wt  spine   fe  8.480867     11  2   wt  spine   fe  8.887062     12  3   wt  spine   fe  8.329797     13  1  exp  brain   cu 10.242652     14  2  exp  brain   cu  9.865984     15  3  exp  brain   cu  9.163728     16  1   wt  brain   cu 11.099385     17  2   wt  brain   cu  9.364261     18  3   wt  brain   cu  9.718322     19  1  exp  spine   cu 10.720157     20  2  exp  spine   cu 11.505430     21  3  exp  spine   cu  9.499359     22  1   wt  spine   cu  9.855950     23  2   wt  spine   cu 10.120489     24  3   wt  spine   cu  9.526252     25  1  exp  brain   zn  9.736196     26  2  exp  brain   zn 11.938710     27  3  exp  brain   zn  9.668625     28  1   wt  brain   zn  9.961574     29  2   wt  brain   zn 10.461621     30  3   wt  brain   zn  9.873667     31  1  exp  spine   zn  9.708067     32  2  exp  spine   zn 10.109309     33  3  exp  spine   zn 10.973387     34  1   wt  spine   zn  8.406536     35  2   wt  spine   zn  7.797746     36  3   wt  spine   zn 11.127984          # use tapply to aggregate     > tapply(elemconc$conc, elemconc[c('elem', 'region')], function(x) x)        region     elem brain     spine      fe Numeric,6 Numeric,6      cu Numeric,6 Numeric,6      zn Numeric,6 Numeric,6          # check whether the order of data has been preserved after aggregation     > x['fe', 'brain']     [[1]]     [1]  8.497498  9.280944  9.726271 11.556397 10.992550  9.711200          # create an external factor for strain grouping     > tmpgeno = rep(c('exp', 'wt'), each=3)     > tmpgeno     [1] "exp" "exp" "exp" "wt"  "wt"  "wt"          # do the t test using the grouping factor     > x = tapply(elemconc$conc, elemconc[c('elem', 'region')], function(x) t.test(x~tmpgeno) )     > x        region     elem brain  spine      fe List,9 List,9      cu List,9 List,9      zn List,9 List,9


title: How would you deseasonalize this data

text: I am dealing with electricity prices, and they have 3 know "seasons":   * Prices during the day follow a pattern;   * Prices in the weekend are lower than during the week;   * In summer prices are different than in winter. The three effects interact as well, i.e. in summer the daily pattern is different than during the winter. At the end of the question the "seasons" are summarised, with the daily effects displayed for each month, full lines are weekdays, dashed lines are weekends. Now what is the best way to deseasonalize this data? I can think of a couple of options:   * Calculate with dummies the known effects, and subtract this from the original data. (This is what I have done, summary graph also attached at the bottom. In blue is the seasonal effect.);   * Work with differences, for example take the 24 hour difference to get rid of the daily effects. But can you do this to take out for example first the weekly effect and afterwards still the daily effect?   * Work with a moving average component that goes back 24 hours of even 7 days. What option do you suggest? Thanks in advance, and for those interested in the subject, this is a video explaining a very particular seasonal effect of electricity use http://www.youtube.com/watch?v=WCAzalhldg8. ![enter image description here](http://i.stack.imgur.com/wZ85Z.jpg) ![enter image description here](http://i.stack.imgur.com/Fsulm.jpg)


title: Restricted Boltzmann Machine for grayscale images

text: I know that RBM's have been used on image data for pre-training neural nets, but all I can find are RBM's on black and white images. How do you apply them to say 256 grayscale?


title: Gibbs sampling with Dirichlet Likelihood

text: I have a sequence of observations that I am representing as proportions:                X1    X2    X3    X4    X5     0.10  0.20  0.50  0.12  0.08     0.07  0.24  0.55  0.04  0.10     ...      The data are coming in as part of a more complicated hierarchical model, for which I need to perform Gibbs sampling. I have all the other conditionals identified, but I do not know how to choose a conjugate prior nor calculate the posterior for a Dirichlet likelihood. I've found this paper which reformulates the Dirichlet as an exponential family member and determines the conjugate prior with posterior calculation. However, I'm not familiar with exponential families and it's never really spelled out what the prior is in a non-exponential family form nor what the posterior should be. **Edit:** Based on this paper (equation 8), the conjugate prior for a Dirichlet is: $$ p(\vec{a}_j) = f(v, \lambda) \left[\frac{\Gamma (\Sigma_{l=1}^D a_{jl})}{\prod_{l=1}^D \Gamma(a_{jl})}\right]^v \prod_{l=1}^D e^{-\lambda_l(a_{jl}-1)} $$ where $f(v, \lambda)$ is a normalization coefficient. Since $f(v, \lambda)$ is intractable to compute, the best you can get is proportionality. I'm fine with that since I'm only interested in sampling-- it's just not clear to me how one would sample from a member of the above family.


title: BUGS with very large sample size

text: I posted this question on the JAGS help discussion, I was seeing if I could get any help here: I wish to fit a survival model using the data below (first two columns are time points, NA represents censoring, third column is number of samples). Is there a way to weight the likelihood by the sample size? If not, I could repeat each time point for given sample size but this would result in probably too larger of an overall sample.               a1=matrix(c(     0.01, 1, 463,     1, 2, 369,     9, 10, 116,     10, 11, 163,     11, 12, 149,     12, 13, 230,     12.5, NA, 150054,     11.5, NA, 146349,     10.5, NA, 118098,     9.5, NA, 41633,     8.5, NA, 30308,     7.5, NA, 25934,     6.5, NA, 29427,     2.5, NA, 37997,     1.5, NA, 40599,     0.5, NA, 43300),ncol=3,byrow=T)


title: Regression with "unidirectional" noise

text: I would like to learn a function $f: \mathbb{R}^n \rightarrow \mathbb{R}$ from data $(x_i, y_i)$ where $y_i = f(x_i) + \epsilon_i$. $\epsilon_i$ can be seen as observation noise; however, it is always smaller or equal to 0 ($\epsilon_i \leq 0$), i.e., the observations are always a lower bound on the true function value. Are there any regression approaches particularly suited for this setting? Is it required to make specific assumptions regarding the distribution of $\epsilon_i$?


title: Is possible to solve this problem with neural network?

text: I have 2 functions f(t) and g(t). I would like to find the function s(t) that minimize the error |f(s(t))-g(t)|^2 Is it possible to estimate s(t) using neural network?   I am new to the field so please try to be clear. Thanks! I am using pybrain so a solution related to pybrain would be perfect


title: How to measure the distance (or divergence - not sure) between data and a probability distribution?

text: If I have generated a set of random data and I wish to measure how well these data fit, e.g. a uniform probability distribution, what are the standard ways to do that? I am not very experienced with this sort of data processing. One thing that comes to mind is Kullback-Leibler divergence, but there you measure distance between two p.d.f.s. Here I have one p.d.f. and then a set of data I would like to compare it to, but which I can probably at best represent as a histogram. I have seen people mentioning the Kolmogorov-Smirnov statistic. What that be suitable here?


title: Covariance and variance relation

text: I collected 3D velocity data in a river. the recorded timeseries of velocity components $(u,v,w)$ contains turbulent and mean part. Reynolds Shear Stress can be estimated from $$u=u'+U$$ $$\langle u'w' \rangle=\text{cov}( u'_{t}, w'_{t} )$$ and to extract velocity fluctuation $$\langle u' \rangle=\sqrt{\text{var}( u'_{t})}$$ $$\langle w' \rangle=\sqrt{\text{var}( w'_{t})}$$ where $\langle\rangle$ shows the averaging over the time span of $t$. **My question is why $\langle u ' \rangle\langle w' \rangle$ is not equal to $\langle u' w' \rangle$** Example: $\langle u'w' \rangle=$ {55.0336, 24.3896, 22.4693, 21.9123, 31.9418, 34.791, 30.9995, -1.12979, -2.76629, -5.60678, -7.23715, -11.1097, -24.4944, -31.5994, -49.414, -92.9571, -97.5096, -91.7745, -141.235, -110.685} $\langle u' \rangle^2=$ {691.168, 438.195, 402.749, 332.111, 470.982, 565.431, 454.04, 192.986, 94.7298, 110.964, 81.9545, 128.034, 198.726, 297.308, 508.916, 910.229, 922.968, 983.665, 1311.84, 1011.1} $\langle w' \rangle^2=${46.2054, 28.1408, 43.3248, 32.9785, 32.092, 40.1115, 28.038, 13.3292, 6.74369, 7.3284, 6.98387, 11.5178, 18.7025, 29.0199, 39.9626, 60.4201, 74.1318, 80.8566, 89.0198, 78.8576}


title: How do I algorithmically determine values of T1 & T2 for canopy clustering?

text: I am trying to use canopy clustering to provide initial clusters for KMeans in mahout. Is there a way to determine / approximate the values of the distance thresholds T1 & T2 algorithmically? Right now I have T1 = 100 and T2 = 1 which doesn't seem to be doing anything good.


title: How would one visualize this data using bar charts?

text: I am using Google's visualization API for drawing charts that show files downloaded each day. I don't know if this is the best idea, but I am thinking of showing download counts on the Y-axis, weekdays (Monday, Tuesday, etc.) on the X-axis, while the legend would show the name of the files. I am facing a problem when the names/legend of files goes out of bound which makes the bars too close and confused even with color variations. So, what alternative visualization technique might be applied in this case?


title: How to compute bias-variance components of linear models using R

text: I am trying to understand the how to compute the components Variance and Bias of the total Error (e.g. http://scott.fortmann- roe.com/docs/BiasVariance.html). Assume that I have available a dataset of 1000 observations. Firstly, I split the dataset into (90-10) training and testing set and I compute the total training and test error as below:               lm1 = lm (Y ~ X1, data = data[1:900,])     lm2 = lm (Y ~ X1 + X2, data = data[1:900,])          train_error1 = sqrt(mean(residuals(lm1)^2))     train_error2 = sqrt(mean(residuals(lm2)^2))          pred_error1 = mean((predict(lm1, data = data[901:1000,]) - data[901:1000,]$Y)^2)     pred_error2 = mean((predict(lm2, data = data[901:1000,]) - data[901:1000,]$Y)^2)      In order to estimate the Variance and Bias components of the above models do I need to split the initial dataset of 1000 observations into K smaller sets and then for each of these smaller set do the above computation? If yes, what is the formula for computing these two components?


title: Alternatives to Cohen's d for non-Gaussian models

text: Cohen's d (or Hedges' g) are often used to compute effect size. They rely on the assumption of homogeneity of variance across samples however. Because of the pooling of variance that they do, I'm also in doubt that they can be useful in the situation when one want to compare groups in generalized linear models with non-Gaussian distribution of errors (Poisson for instance), possibly overdispersed for some groups. Am I right, and if yes, what should be an adequate measure of effect size in this situation ? Ideally I'm looking for a metric which could be used for both binomial and Poisson or negative binomial error distribution models.


title: How to compare proportions, dependent samples with different sample size

text: I have a dataset size=8755 that has many missing values. If I remove all the missing values, the complete dataset has size=2986. I'm trying to compare the proportion of both dataset in order to understand if those are significantly different or not. For instance compare the proportion of HIV in both groups. The problem is since it's the same dataset, although with fewer elements I don't think these are independent samples, then I'm not gonna use the z-test to compare proportions. Being dependent samples I think I need to use the NcNemar test however only 2986 samples have a match in the full dataset and this match is equal, since its the same person. Therefore I'm also don't think it's a good solution. I'm a bit lost about the best approach to tackle this problem. Any pointers would be great. Thank you!


title: Confidence intervals for cross-validated statistics

text: How does one calculate confidence intervals of cross-validated estimates? For an epidemiological paper we use cat. and cont. NRI, IDI, and difference in C index for comparison of two Cox models. The reviewer suggested showing only cross-validated estimates **and their 95% confidence intervals**. My ideas include taking the appropriate quantiles of the CV resamples, calculating the SE of those resamples and constructing Wald intervals, or bootstrapping the CI of the resamples' mean or median. But somehow these all seem phony.


title: PROC LOGISTIC, GENMOD

text: How can I get the standard error for odds ratios. I know how to obtain ORs, either by exp the estimate or by using LSMEANS or SLICE statements. Unfortunately wasnt able to find the options to have SE for ORs in the output datasets. Only the estimate SE is in there.


title: Loess through quantiles

text: I have the following data points: ![enter image description here](http://i.stack.imgur.com/d0GTF.jpg) The thin blue line is the aggregated average for each level of X (ages), the thicker blue line is a Loess fit with R. Now I have two questions regarding this Loess function:   * The number of data points per level of X is quite low at the boundaries of X. Does the Loess takes this into account? Or in other words, can I give weights to X categories depending on the number of data points?   * On the graph, you can see the dashed line is the 5% quantile. How would you draw a fit through these quantiles? Is there a way I can still take all of the data points into account? And if not, there will be a even less data at the boundaries of X, so I am afraid question 1 is even more important. Thank you! UPDATE: Apparently, there is an R package for this, being quantreg. I have however difficulties finding the optimal parameters for non parametric regression, ref the image below. Does anybody have a better idea how to model this? The green plot in the middle is the average with classic LOESS, the other curves are the 95, 75, 25 and 5 % quantile modeled with the package. Thanks again! ![enter image description here](http://i.stack.imgur.com/4omBg.jpg)


title: Self study question on Uniform Distribution

text: I am attempting a self-study question on uniform distribution but was not able to make sense of the answers provided after taking an attempt. Appreciate any help and guidiance please. **Question** > An energy company provides electricity with voltage uniformly distributed > between 220.0 volts and 250.0 volts. Any voltage level being more than 242.0 > volts is dangerous and unacceptable. What is the probability that a randomly > selected voltage level is between 245.5 and 248.5, given that it is under > the unacceptable level? **Attempt** ![enter image description here](http://i.stack.imgur.com/dNnXK.jpg)


title: Coding Categorical Variables

text: Suppose I am building a linear model in R. I will be doing standard OLS. I have 10 dummy variables (predictors) that correspond to different regions. 6 of these regions are in California, and the other 4 are in Texas. For example for my Northern California dummy variable is a 1 if the observation comes from California and 0 if not. I was thinking of creating a categorical variable with the values 1:10 with each number corresponding to a different region. Would this affect my analysis negatively in any way? I have a feeling that this would just distort the interpretation of my estimated coefficient. I also have a dummy variable for California and Texas.


title: Multivariate normal distribution log-likelihood

text: Let's say we have two 2x5 matrices A and B, and a 4x5 matrix C which is generated by concatenating A and B (C=[A;B] in MATLAB representation). Would the sum of the log-likelihoods of A and B under a 5-variate Gaussian distribution with mean $\mu$ and covariance $\Sigma$ be equal to the log- likelihood of C under the same Gaussian distribution? Since rows correspond to samples, I expect them to be equal, but I see that they are not equal when I do the calculations. Can somebody explain that? Thanks!


title: Random Forest for data imputation

text: Currently I am using Random Forest approach for Missing Values Imputation from missForest package in R. I faced the following problem: the algorithm works much longer than any other imputation approaches like KNN or SVD. Does anybody know how to overcome this problem or maybe anyone had experience with using it in other packages? I'd really appreciate any help cause I even failed to find the full description of the approach in the internet.


title: Psychology / heuristics of $\Pr(A)\le \Pr(B)$ when $A$ is a subset of $B$

text: A while ago I read about the flaws people have in assigning probabilities to events based on their beliefs and experiences. For some reason, I started thinking about that today, and realized that I don't recall the main points of the argument, and really want to understand it. I'll try to state it here as I recall it. Do you see any obvious errors in the questions or in my conclusion?               R = all rock stars     M = all rock stars who are Grand Masters (GM) in chess      What is the probability of a rock star to play chess after concert? Doesn't seem likely, right? Say, $\Pr(R) = 0.2$ What is the probability of a GM rock star to play chess after concert? A lot more likely: $\Pr(M) = 0.8$ But, $\Pr(R)$ must be $> \Pr(M)$ as $M$ is a subset of $R$, right? So one or both of the probabilities has to be wrong.   1. So, does this mean I have to modify $\Pr(R), \Pr(M)$ or both?   2. Or in order to have said $0.2$ I had to know the size of $M$?   3. Or is there a problem with how the question is stated? Also, if you can point me to a book or a paper on the subject, I'd really appreciate that.


title: Self-Study question on Probability

text: Good evening, I am working on a self-study probability question as follows: > A student is taking two courses, math and statistics. The probability the > student will pass the math course is 0.70, and the probability of passing > the statistics course is 0.80. The probability of failing both courses is > 0.1. What is the probability that the student will pass the statistics > course but not the math course? I have made several attempts to solve the questions, but my effort seems futile. My best answer was 0.25, while the answer provided was 0.2. Appreciate guidance and advice please. ![enter image description here](http://i.stack.imgur.com/HFN8N.jpg)


title: Most Powerful Test; Two-Parameter Normal Distribution

text: Is it possible to show that the two-parameter Normal distribution has monotone likelihood ratio? EDIT: This is actually part of a larger problem. We have a random sample from $\mathcal N(\mu, \sigma^2)$ and are testing: \begin{align} Ho\\!:\quad&\sigma <\sigma_o , &\mu &\in \mathbb{R} \\\ Ha\\!:\quad&\sigma = \sigma_1 , &\mu &= \mu_1 \end{align} Construct a most powerful level-$α$ test. My initial thought was to show $\mathcal N(μ,σ^2)$ has the monotone likelihood ratio (MLR), and from there we can construct a uniformly most powerful (UMP) test via the Karlin-Rubin theorem. Any suggestions on general strategy would be appreciated!


title: Probability that x is greater than some value and y is less than some value where values come from normal distribution?

text: I was attempting the following self-study question. After repeated tries, my answers were at best **0.0944** , however, the model answer seemed to be otherwise at **0.1889**. Will appreciate any pointers as I simply can't figure out even after repeated attempts. **_Question_** > According to the 2010 census in a city, university-educated employee earned > a mean annual income of USD 61,000 & standard deviation of USD4,000. If the > incomes are normally distributed and two university-educated employees are > randomly selected, what is the probability that one of the them earns more > than USD66,000 per year and the other earns less than USD66,000 per year? **_My Attempt_** ![enter image description here](http://i.stack.imgur.com/FqT0o.jpg)


title: What is an appropriate formula for residuals calculation in a model describing a single-factor repeated measures design?

text: For the repeated measures design _D.C.Montgomery_ in his "Design and Analysis of Experiments" book provides the mathematical / statistical (linear) model: $y_{ij} = \mu + \tau_i+ \beta_j + \epsilon_{ij}$, where $\tau$ is the effect of the $i$th treatment, $\beta$ is a parameter associated with the $j$th subject, and $\epsilon_{ij}$ is the independent error.   In a (document, p. 2) it is stated that to check the model's assumptions we should operate with residuals calculated as ( _A_ ): $$\hat\epsilon_{ij} = Y_{ij} - \bar Y_{i.} - \bar Y_{.j} + \bar Y.$$ The formula ( _B_ ) on Wikipedia is $$\hat\epsilon_i = Y_i - \bar Y.$$ My question is:   what formula is used to correctly calculate residuals in the model describing a single-factor repeated measures design to check its assumptions (multivariate normality, homoscedasticity, sphericity, linearity, additivity)? Thank you. p.S.   _E.Vonesh_ and _V.M.Chinchilli_ in their "Linear and Nonlinear Models for the Analysis of Repeated Measurements" monograph state that the single-sample RM design may be written as $y_{i} = \mu + \epsilon_{i}$.


title: In EM derivation why can I sum over the iid variables in the conditional expectation?

text: In EM when you take the expectation: $E[\log P(y,x \mid \theta)\mid x, \theta']$ $= \sum\limits_yP(y\mid x, \theta') \log P(y,x\mid \theta)$ I understand this but the following part I don't understand. We know that $x$ is a vector of $n$ independent identically distributed variables. Also y is a single value- the cluster x is assigned to. How is the above equivalent to: =$\sum\limits_{i=1}^{n}\sum\limits_y P_{\theta'}(y\mid x_i)log P_\theta(y,x_i)$ I realize that multiplication within a log expands to summation but how does one replace $P_{\theta'}(y\mid x)$ with $P_\theta'(y,x_i)$? Reference: slide 13 http://cs.nyu.edu/~eugenew/publications/em-talk.pdf


title: Graphical Models and Explaining Away?

text: So I have fundamental confusion about graphical models. Suppose the following graphical model is given:![enter image description here](http://i.stack.imgur.com/L0C2D.png) Now the question is do we have the following equality?: $$p(m_2|\alpha,\beta,y_3)=p(m_2|\alpha,\beta)$$ If so, why exactly and could you please introduce me a source explaining it?


title: Combining posterior distributions

text: I have 5 different posterior distributions (mcmc samples) which all estimate the same parameter beta. The 5 models are all obtained from 5 independent standardized datasets but estimate the same parameter. What I want to obtain is a single posterior distribution, combining information of my five original distribution. By instinct, following the rules of probability, I would say I can obtain a joint posterior distribution simply as the product of my 5 different posterior distributions. Is this correct, and if yes, how to do this in practice? Is there any BUGS/JAGS or R code anybody is willing to share? Thanks a lot, achaz


title: Derive probability from linear regression

text: I have a linear function that coorelates a measurement with error. As values of the measurement increase, so does the mean error. Error is normally distributed. I'd like a way to derive a probability as to how certain my measurement is. How can I do this? One thought was to transform the distribution of error at different measurements to a z-score and then find a way to derive a probability from that, but I'm not sure if that is a good approach. Any suggestions is appreciated.


title: Problems and alternatives to using a four month lag in a linear regression model

text: Let's say I'm constructing a linear model with the intention of predicting automobile sales volume. Let's say that the consumer auto purchasing cycle takes 4 months, and so we'd 'lag' each observation by four data points. So if out data was monthly, out data may end up looking like the following               month   lag_month     sales      visits       gas price     jan           feb     march     april    jan           500        50000        3.55     may      feb           550        45000        3.87      Given the lag in consumer shopping behavior, I want/'need' to lag the variables. However, I also feel like I have to lag the variables for the sake of prediction. Let's say I run a regression using this data and get the following estimates, I could input feb numbers to prediction sales four months from now, yes/no?               sales = -0.05 + 2*(visits) + .0.35*(gas)     sales = -0.05 + 2*(550) + .0.35*(3.87)  - four month lag to predict four months from now      The question I have relates to the process. This process seems incredibly poor (statistically unsound) and I'm wondering what are the problems with utilizing such an approach (four month lag). What are alternatives to lagging when the goal is prediction?


title: Question on how to normalize regression coefficient

text: Not sure if normalize is the correct word to use here, but I will try my best to illustrate what I am trying to ask. The estimator used here is least squares. Suppose you have $y=\beta_0+\beta_1x_1$, you can center it around the mean by $y=\beta_0'+\beta_1x_1'$ where $\beta_0'=\beta_0+\beta_1\bar x_1$ and $x_1'=x-\bar x$, so that $\beta_0'$ no longer has any influence on estimating $\beta_1$. By this I mean $\hat\beta_1$ in $y=\beta_1x_1'$ is equivalent to $\hat\beta_1$ in $y=\beta_0+\beta_1x_1$. We have reduced equation for easier least square calculation. How do you apply this method in general? Now I have the model $y=\beta_1e^{x_1t}+\beta_2e^{x_2t}$, I am trying to reduce it to $y=\beta_1x'$.


title: Mixed model with lmer: Variance of residuals should give the same as level 1 variance?

text: I expected that the variance of residuals from a mixed model computed by, for example, lmer should give the same as the residual variance from the summary output.               # -- Data simulation     group  <- rep(letters[1:10], each=10)     groupx <- rep(rnorm(10,0,5), each=10)     x0     <- 1:10     y0     <- rep(3*x0,10)     y      <- y0+groupx + rnorm(100,0,3)     data0  <- data.frame(grp=group,x=x0,y=y)     head(data0)          # -- library(lattice)     xyplot(y ~ x|grp, data=data0)          # -- Model     f0  <- lmer(y ~ x + (1|grp), data= data0, REML=FALSE)     f0s <- summary(f0)     f0s          # -- Compare level 1 variance (Residual):     as.data.frame(VarCorr(f0))[,c("grp","vcov")]     # -- with variance of residuals:     var(resid(f0))      In my example I got:               > as.data.frame(VarCorr(f0))[,c("grp","vcov")]            grp     vcov     1      grp 20.07274     2 Residual 10.04903          > var(resid(f0))     [1] 9.183874      I expected that `var(resid(f0))` should be exactly the same as the variance of Residuals from summary. Could someone gives me a hint what is wrong in my argumentation?


title: Sample size calculation for multilevel, longitudinal experiment

text: I'd be grateful for any advice about how to do a sample size calculation for the following design. We are interested in the impact of two variations of an intervention on hospital in-patients' perceptions of the quality of the nursing care they received. Therefore, there will be three conditions, including control. Research will be carried out in a number of hospitals; the hospitals themselves will probably not be sampled from a population: they will be hospitals that are willing to participate in the study. There will likely be a small number (say three). Within each hospital, a number of units (called wards in the UK) will again not be sampled, but they will be randomly assigned to one of the three conditions, so that there are equal numbers of units in each hospital in each arm. The interventions are at the unit level, in that they are aimed at the nurses that work in the unit. However, the outcome variable is derived from a survey of a sample of patients that were treated in the unit. New samples of patients will be surveyed at 3 monthly intervals for up to 18 months. In a pilot study, the outcome variable was the mean patient satisfaction score for each unit at each time. The statistic of interest is the interaction between experimental condition and number of months since the start of the experiment. My questions are:   1. How to calculate the required number of patients to sample in each unit at each time?   2. How to calculate the required number of units per condition in each hospital?   3. Is it appropriate to use mean scores at the unit level as the outcome measure, or would it be better to use individual patient responses as the outcome? Thanks for any help! David


title: Compute gaussian fit corresponding to a discrete variable

text: I collected data for a discrete probability distribution and I have the pairs (value, probability).               Value    Probability      90      0,0033     125      0.0204     180      0.0847     250      0.2516     355      0.4653     500      0.175     710      0.0015      I have been asked to compute a gaussian fit over that data but I am having troubles. Is it possible (and does it make sense) to fit a gaussian distribution over it? I tried to build the corresponding gaussian distribution by computing mean (`334`) and standard deviation (`100`) of my data but, of course, it does not work. Gaussian density function values are related to continuous values while the probability I have are continuous. I guess that I have too few discrete values to compute the gaussian distribution, but I am not really into statistics so I am not sure and I'd like to double check with you that I'm not missing something.


title: Breakpoints and non linear regression

text: I have done some research about breakpoints (I am not a statistician) and I found out about the breakpoint, strucchange and segmented packages in R allowing to find breakpoints assuming linear model. However, I would like to fit a periodic time series with a non linear (periodic) model, and I was wondering how I could find breakpoints for this model in R. > An example (but my model is more complicated than that) : Y ~ (a + b*t)*( > c*cos(wt) + d*sin(wt) + e*cos(2wt) + f*sin(2wt) ). I guess I am more interested in the breakpoints of the "amplitude" of my periodic function, so that I could assume a model of the form: > Y ~ (a + b*t)*f(t), with |f(t)| <= 1, where f is a periodic function to be > fitted to a non linear model, but where no breakpoints should occur. So basically, the breakpoints would only happen in the (a,b) pair of coefficients, which would be a linear regression. However, as f is unknown, this makes things harder... Thank you for your advices !


title: How one can improve correlation (higher rho values) on non normal distributions?

text: No matter what class or book I see, the normal distribution is always advocated. I have been trying to find ways to assess relationship between variables for non normal distributions, most of which are exponential. Is there any other method other than spearman correlation to deal with such things?


title: Rank-one nonnegative matrix factorization

text: For non-negative matrix factorization with Frobenius norm: $\min\limits_{U\in\mathbb{R}_+^{m\times r}, V\in\mathbb{R}_+^{r\times n}}||A-UV||_2^2$, $A\in\mathbb{R}_+^{m\times n}$ $r=1$ is a very special case because the optimal solution is given by SVD and there is no need to perform successive minimization by $V$ with fixed $U$ and vice versa. Is there anything special for $r=1$ when minimization is performed with respect to some other norm, for example, $l_1$, Kullback-Leibler divergence or any other divergence, or should it just be done by the algorithms for the general case?


title: What are the multidimensional versions of median

text: What are the multidimensional versions of the median and what are their pros and cons? I confess this doesn't have a single answer, but I think it is a useful question to ask and will be a benefit to others as well. How stable it is (i.e. how many samples one needs to get a reasonable estimate of it) is one potential, but not necessary, pro and con issue, i.e. if you know that the number of samples necessary grows exponentially with each increase in dimension making it useful in 10 dimensions but effectively useless in 200 dimensions, that would be useful to know. (I would kind of expect all of them have that, really, because there is so much "freedom" in 200 dimensions it just takes a bazillion points to nail the basic middle of a distribution down along that many directions). **Note:** I found this question after accepting an answer for this. I hadn't seen it before because I didn't realize there was both a multivariate-analysis and a multivariable tag. I am not sure why geometric median was an unpopular answer there (1 vote and near the bottom) but a popular answer here.


title: Deviance Statistic in HLM analysis

text: I am analyzing a multilevel model using the software HLM. I am trying to find out more about the **Deviance** statistic that I increasingly see being reported as a measure of the explanatory value/effect size of a model. In general, I know that HLM output reports this statistic and that the smaller it is the better. However:   * How does one use/interpret this statistic for a given analysis?   * How does one know if the deviance one obtains is at an acceptable level?


title: Alternatives to logistic regression in R

text: I would like as many algorithms that perform the same task as logistic regression. That is algorithms/models that can give a prediction to a binary response (Y) with some explanatory variable (X). I would be glad if after you name the algorithm, if you would also show how to implement it in R. Here is a code that can be updated with other models:               set.seed(55)     n <- 100     x <- c(rnorm(n), 1+rnorm(n))     y <- c(rep(0,n), rep(1,n))     r <- glm(y~x, family=binomial)     plot(y~x)     abline(lm(y~x),col='red',lty=2)     xx <- seq(min(x), max(x), length=100)     yy <- predict(r, data.frame(x=xx), type='response')     lines(xx,yy, col='blue', lwd=5, lty=2)     title(main='Logistic regression with the "glm" function')


title: Standard errors for estimates of smoothing parameters

text: My question is based on the " _forecast_ " package for R used in Forecasting with Exponential Smoothing. The State Space Approach \- Hyndman et al. 2008. I am using the `ets` function to estimate the parameters of a model. Is there a way to obtain standard errors for the estimates of the smoothing parameters?


title: Distribution with random parameter

text: I got a problem with Exercise 3.1(a). Since there's an example right before the exercise I assumed that X conditioned on M has the same distribution as in the example. I don't understand their computation of P(X=k), how did that gamma-function appear? (I'm new to the gamma-distribution) and can I use that in my solution of 3.1(a)? ![enter image description here](http://i.stack.imgur.com/8JRLU.jpg) Thank you in advance!


title: Tiny (real) datasets for giving examples in class?

text: When teaching an introductory level class, the teachers I know tend to invent some numbers and a story in order to exemplify the method they are teaching. What I would prefer is to tell a real story with real numbers. However, these stories needs to relate to a very tiny dataset, which enables manual calculations. Any suggestions for such datasets will be very welcomed. Some sample topics for the tiny datasets:   * correlation/regression (basic)   * ANOVA (1/2 ways)   * z/t tests - one/two un/paired samples   * comparisons of proportions - two/multi way tables   * So on...


title: Using priors to detect an effect? logistic Bayesian regression

text: I have designed an idea and am looking for similar approaches in other literature/areas or if I have applied the Bayesian concepts wrongly. Here is a statement of my problem:   1. I am modeling the influence of a single variable X on the odds of a death event. The base chance is 15%.   2. I am going to take the approach of using log(p/1-p)=alpha+beta*x to model the odds ( Logistic regression ).   3. I measure X in terms of 1 unit increases, so 1 increase in X will typically increase the chances of death by an additive amount. ( 1 increase in X -> +5% more risk, or +10% or +20% Note we do have 1.5,1.9.. but I was just getting across the idea of the effect of X).   4. I only have 17 data points relating X and death. That's not enough given my power analysis.   5. I construct 4 sets of priors, not based on the data, just based on an agreeable reasonable effect between the scientist and myself. I'm lacking literature here, but the researchers do feel comfortable about my approach. The first prior (normal) assumes a base chance of 15%, and an increasing additive effect of X around 5%, the second assumes a base chance of 15% and an additive effect of increasing X around 10%. ( Normal priors with certain standard deviations control the "around" part ). etc.   6. Run the bayes analysis using data I have and the 4 separate priors. It may be that one posterior based on a prior 5% effect tells me the effect is likely higher. While the 15% effect posterior tells me the effect is likely lower. I realize I may possibly get conflicting information.   7. From this I believe that the effect is somewhere between 5% and 15% and on another new sampled dataset, repeat this with the priors being near those ranges to narrow it down.   8. I can specifically quantify the judgement I used in 6, by estimating the posterior probability that the effect is more than 5%, etc. Any thoughts, criticisms? I think this approach is valid. I realize the effect of the X variable could be different ( ie not additive ). Thank you.