[R] Computing Median for Subset of Data
From my larger data set I created a subset of it by using: subset_1 - subset(timeuse, IndepTrans = 1, Physical = 1) where my larger data set is timeuse and the smaller subset is subset_1. The subset was conditioned on IndepTrans equaling 1 in the data and Physical equaling 1 as well. I want to be able to compute the median of a variable first for the larger data set timeuse then for the subset file subset_1. How do I identify to R which data set I'm wanting the median computed for? I've tried many possibilities but for some reason can't figure it out. Thanks, Matt. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Computing Median for Subset of Data
On Jun 2, 2013, at 08:08 , Matt Stati wrote: From my larger data set I created a subset of it by using: subset_1 - subset(timeuse, IndepTrans = 1, Physical = 1) That's not going to work. Try subset(timeuse, (IndepTrans == 1) (Physical == 1)) (Whenever you do things like this, run summary(subset_1) to check.) -- Peter Dalgaard, Professor, Center for Statistics, Copenhagen Business School Solbjerg Plads 3, 2000 Frederiksberg, Denmark Phone: (+45)38153501 Email: pd@cbs.dk Priv: pda...@gmail.com __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Computing Median for Subset of Data
On 02-06-2013, at 08:08, Matt Stati mattst...@yahoo.com wrote: From my larger data set I created a subset of it by using: subset_1 - subset(timeuse, IndepTrans = 1, Physical = 1) where my larger data set is timeuse and the smaller subset is subset_1. The subset was conditioned on IndepTrans equaling 1 in the data and Physical equaling 1 as well. I want to be able to compute the median of a variable first for the larger data set timeuse then for the subset file subset_1. How do I identify to R which data set I'm wanting the median computed for? I've tried many possibilities but for some reason can't figure it out. ?with with(timeuse, median(…)) with(subset_1, median(…)) Berend __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] object not found
Hi Jim
It is defined and initialised:
line 69 in the differential equation: dmCaT_soma_AB - (mCaTx(v_soma) -
mCaT_soma_AB)/taumCaT(v_soma)
and initialised in line 83
dmCaT_soma_AB=0
However, l just noticed the mistake. It needs to be initialised without the
d, i.e. mCaT_soma_AB=0 and not dmCaT_soma_AB
Regards
Jannetta
Regards
Jannetta
On 2 June 2013 02:14, jim holtman jholt...@gmail.com wrote:
I just did a search through the source and never did find a place that
'mCaT_soma_AB' was defined. How can you expect to use it if it has not
been defined.
On Sat, Jun 1, 2013 at 8:50 PM, Jannetta Steyn janne...@henning.orgwrote:
Hi All
I am trying to implement a neural model in R from a paper but I have run
into a bit of a problem. The model is getting so complex that I can't
really test it in bits any more. At the moment I get this error:
Error in eval(expr, envir, enclos) : object 'mCaT_soma_AB' not found
I think I have been staring at it too long because I can't see what I have
done wrong. Can anyone perhaps spot the error I have made. The full code
is
posted below.
Many thanks
Jannetta
# TODO: Add comment
#
# Author: a9912577
###
library(deSolve)
ST - function(time, init, parms) {
with(as.list(c(init, parms)),{
#functions to calculate activation m and inactivation h of the currents
mNax - function(v) 1/(1+exp(-(v+24.7)/5.29));
taumNa - function(v) 1.32 - (1.26/(1+exp(-v+120)/25));
hNax - function(v) 1/(1+exp((v+48.9)/5.18));
tauhNa - function(v) (0.67/(1+exp(-(v+62.9)/10))) *
(1.5+1/(1+exp(v+34.9)/3.6));
mCaTx - function(v) 1/(1+exp(-(v+25)/7.2))
taumCaT - function(v) 55-(49.5/(1+exp(-(v+58)/17)))
hCaTx - function(v) 1/(1+exp((v+36)/7))
tauhCaT_AB -function(v) 87.5-(75/(1+exp(-(v+50)/16.9)))
tauhCaT_PD -function(v) 350-(76/(1+exp(-(v+50)/16.9)))
mCaSx - function(v) 1/(1+exp(-(v+22)/8.5))
taumCaS - function(v) 16-(13.1/(1+exp(-(v+25.1)/26.5)))
mNapx - function(v) 1 / (1+exp(-(v+26.8)/8.2))
taumNap - function(v) 19.8-(10.7/(1+exp(-(v+26.5)/86.)))
hNapx - function(v) 1/1+exp((v+48.5)/4.8)
tauhNap - function(v) 666-(379/(1+exp(-(v+33.6)/11.7)))
mhx -function(v) 1/(1+exp((v+70)/6))
taumh - function(v) 272+(1499/(1+exp(-(v+42.2)/8.73)))
mKx - function(v) 1/(1+exp(-(v+14.2)/11.8));
taumK - function(v) 7.2-(6.4/(1+exp(-(v+28.3)/19.2)))
# AB soma
iCaT_soma_AB - gCaT_soma_AB * mCaT_soma_AB ^ 3 * hCaT_soma_AB * (v_soma -
ECaT_soma_AB)
iCaS_soma_AB - gCaS_soma_AB * mCaS_soma_AB ^ 3 * (v_soma - ECaS_soma_AB)
iNap_soma_AB - gNap_soma_AB * mNap_soma_AB ^ 3 * hNap_soma_AB * (v_soma -
ENap_soma_AB)
ih_soma_AB - gh_soma_AB * mh_soma_AB ^ 3 * hh_soma_AB * (v_soma -
Eh_soma_AB)
iK_soma_AB - gK_soma_AB * mK_soma_AB ^ 4 * mK_soma_AB * (v_soma -
EK_soma_AB)
iKCa_soma_AB - gKCa_soma_AB * mKCa_soma_AB ^ 4 * (v_soma - EKCa_soma_AB)
# Total current for Calcium
totalICa - iCaT_soma_AB + iCaS_soma_AB
# Differential equations
dCaConc_soma - (-F_AB * totalICa - CaConc_soma + C0_AB)/tauCa_AB
mKCax_AB - function(v, CaConc_soma)
(CaConc_soma/(CaConc_soma+30))*(1/(1+exp(-(v+51)/4)))
mKCax_PD - function(v, CaConc_soma)
(CaConc_soma/(CaConc_soma+30))*(1/(1+epx(-(v+51)/8)))
taumKCa -function(v) 90.3 - (75.09 / (1+exp(-(v+46)/22.7)))
mAx - function(v) 1/(1+exp(-(v+27)/8.7))
taumA - function(v) 11.6-(10.4/(1+exp(-(v+32.9)/15.2)))
hAx - function(v) 1 / (1+exp((v+46.9)/4.9))
tauhA - function(v) 38.6 - (29.2/(1+exp(-(v+38.9)/26.5)))
mProcx - function(v) 1 / (1+exp(-(v+12)/3.05))
taumProc - 0.5
# Currents as product of maximal conducatance(g), activation(m) and
inactivation(h)
# Driving force (v-E) where E is the reversal potential of the particular
ion
# AB axon
iNa_axon_AB - gNa_axon_AB * mNa_axon ^ 3 * hNa_axon * (v - ENa_axon_AB)
iK_axon_AB - gK_axon_AB * mK_axon ^ 4 * (v - EK_axon_AB)
iLeak_axon_AB - gLeak_axon_AB * (v - ELeak_axon_AB)
dv - (0 - iNa_axon_AB - iK_axon_AB - iLeak_axon_AB) / C_axon_AB
dmNa_axon_AB - (mNax(v) - mNa_axon_AB)/taumNa(v)
dhNa_axon_AB - (hNax(v) - hNa_axon_AB)/tauhNa(v)
dmK_axon_AB - (mKx(v) - mK_axon_AB)/taumK(v)
dv_soma - (I - iCaT_soma_AB - iCaS_soma_AB - iNap_soma_AB - ih_soma_AB -
iK_soma_AB - iKCa_soma_AB)
dmCaT_soma_AB - (mCaTx(v_soma) - mCaT_soma_AB)/taumCaT(v_soma)
dhCaT_soma_AB - (hCaTx(v_soma) - hCaT_soma_AB)/tauhCaT_AB(v_soma)
dmCaS_soma_AB - (mCaSx(v_soma) - mCaS_soma_AB)/taumCaS(v_soma)
dmNap_soma_AB - (mNapx(v_soma) - mNap_soma_AB)/taumNap(v_soma)
dhNap_soma_AB - (hNapx(v_soma) - hNap_soma_AB)/tauhNap(v_soma)
dmh_soma_AB - (mhx(v_soma) - mh_soma_AB)/taumh(v_soma)
dmK_soma_AB - (mKx(v_soma) - mK_soma_AB)/taumK(v_soma)
dmKCa_soma_AB - (mKCax_AB(v_soma,CaConc_soma) -
mKCa_soma_AB)/taumKCa(v_soma)
list(c(dv,dv_soma,dmNa_axon_AB, dhNa_axon_AB, dmK_axon_AB, dmCaT_soma_AB,
dhCaT_soma_AB, dmCaS_soma_AB, dmNap_soma_AB, dhNap_soma_AB, dmh_soma_AB,
dmK_soma_AB, dmKCa_soma_AB,dCaConc_soma))
})}
## Set initial state
Re: [R] object not found
but that is not the variable that is not found. mCaT_soma_AB Sent from my iPad On Jun 2, 2013, at 4:12, Jannetta Steyn janne...@henning.org wrote: mCaT_soma_AB __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] metafor mixed-effect models: testing different predictor combinations for lowest AIC
At 22:12 01/06/2013, emhanlon wrote: Good afternoon, I am trying to use a mixed-effect model to find the best predictors of wildlife response (fledgling success, offspring growth rate, etc...) to human disturbance. I am using a Log Response Ratio as my measure of effect size along with study variance, and testing to see how: habitat, diet, mass, generation time, sociality, etc... influence how a species responds to disturbance. I have 9 total predictors, but a sample size that would only allow for 3, at most, in a model. Does anyone know of any code that could allow me to automatically test all combinations of 1,2, or 3 predictors to get the model with the lowest overall AIC? You could use the facility that mods can be a matrix and use combn to generate all the possible combinations of 9 taken 1, 2, 3 at a time. Each time you call ram.uni just save whatever you want to use as your criterion. I am not sure how good an idea that is and I would advise looking at all the solutions near (for some meaning of near) the apparently best one. You might be surprised how close they are which rather casts doubt on the wisdom of selecting variables. #Code I had been using in different combinations manually: res - rma(LRR, est_var, mods = cbind(Gentime, Sociality, Season), data = data) #Also, if the best model has more than one predictor, how can you get specific mean effect sizes for different levels of the #variable, especially if it is categorical (non-ranked)? Thank you so much for your help, Edward -- View this message in context: http://r.789695.n4.nabble.com/metafor-mixed-effect-models-testing-different-predictor-combinations-for-lowest-AIC-tp4668486.html Sent from the R help mailing list archive at Nabble.com. Michael Dewey i...@aghmed.fsnet.co.uk http://www.aghmed.fsnet.co.uk/home.html __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Windows Vista computer will not use changes I made to R console file
Le samedi 01 juin 2013 à 15:28 -0400, Julie Royster a écrit : My newer laptop died so I installed R on an older Vista laptop. I did follow the online instructions about giving myself full control under the R properties security tab. I made changes in the R console file under the etc folder to use bold type in color Black, but it is not doing what I ask. I am still getting blue type that is not bold. This previously worked great on the other computers I have used. Can anyone tell me what to do? THANKS!!! Julie in Raleigh NC Maybe you have changes saved to a Rconsole file in your home folder? What happens if you make changes from Rgui and save them? Regards [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Multivariate EWMA covariance estimator?
Hi,
since I want to calculate the VaR of a portfolio consiting of 4 assets
(returns saved into eonreturn,henkelreturn and so on) I have to
estimate the covariance matrix. I do not want to take the rectangular
version with equal weights, but the exponentially weighted moving
average in a multivariate version. I want to estimate a covariance
matrix at every time point t. Then I want to comput the VaR at this
time point t. Afterwards, I will look at the exceedances and do a
backtest.
I tried to implement it as follows (data attached):
lambda-0.9
summe2-0
dummy2-0
covestiexpo-list(NA)
meanvalues-NA
for(i in 101:length(eonreturn)){
meanvalues-matrix(c(mean(eonreturn[(i-100):(i-1)]),mean(henkelreturn[(i-100):(i-1)]),mean(siemensreturn[(i-100):(i-1)]),mean(adidasreturn[(i-100):(i-1)])),4)
for(a in 1:100){
dummy2-lambda^(a-1)*t(datamatrix[(i-a),]-t(meanvalues))%*%(datamatrix[(i-a),]-t(meanvalues))
summe2-summe2+dummy2
}
covestiexpo[[i]]-(1-lambda)*summe2
}
So the covestieexpo[[101]] would be the covariance estimate for the
101th day, taking into account the last 100 observations. Now, the
problem is, that there seems to be something wrong, since the
covariance estimates are cleraly wrong, they seem to be too big. At
the beginning, compared to the normal covariance estimate the
difference is as follows:
covestiexpo[[101]]
[,1][,2][,3][,4]
[1,] 0.004559042 0.002346775 0.004379735 0.003068916
[2,] 0.002346775 0.001978469 0.002536891 0.001909276
[3,] 0.004379735 0.002536891 0.005531590 0.003259803
[4,] 0.003068916 0.001909276 0.003259803 0.003140198
compared to cov(datamatrix[1:100,])
[,1] [,2] [,3][,4]
[1,] 0.0018118239 0.0007432779 0.0015301070 0.001119120
[2,] 0.0007432779 0.0008355960 0.0009281029 0.000754449
[3,] 0.0015301070 0.0009281029 0.0021073171 0.001269626
[4,] 0.0011191199 0.0007544490 0.0012696257 0.001325716
So already here, it is obvious, that something is not correct, if I
look at a period far ahead:
covestiexpo[[1200]]
[,1] [,2] [,3] [,4]
[1,] 0.5312575 0.1939061 0.3419379 0.2475233
[2,] 0.1939061 0.3204951 0.2303478 0.2022423
[3,] 0.3419379 0.2303478 0.5288435 0.2943051
[4,] 0.2475233 0.2022423 0.2943051 0.4599648
you can see, that the values are way too large, so where is my mistake?
Thanks a lot for your help!
--
Neumann, Conrad
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] Robust GAM that covers Gaussian Distributions
Dear All, I was looking the r-archives and crantastic...for a package that has a robust approach to generalized additive models. I found two packages robustgam and rgam but their implemented functions cover only binomial and poisson distributions (pls correct me if I am wrong). I would greatly appreciate if anyone could share with us other packages or robust approaches of general additive modeling that might have a better performance with small data sets (n = 50 -100 records). Thank you very much all for reading this message. I am hoping and looking forward to receiving your reply. Sincerely, Christos Giannoulis [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Windows Vista computer will not use changes I made to R console file
On 02.06.2013 13:49, Milan Bouchet-Valat wrote: Le samedi 01 juin 2013 à 15:28 -0400, Julie Royster a écrit : My newer laptop died so I installed R on an older Vista laptop. I did follow the online instructions about giving myself full control under the R properties security tab. I made changes in the R console file under the etc folder to use bold type in color Black, but it is not doing what I ask. I am still getting blue type that is not bold. This previously worked great on the other computers I have used. Can anyone tell me what to do? THANKS!!! Julie in Raleigh NC Maybe you have changes saved to a Rconsole file in your home folder? What happens if you make changes from Rgui and save them? Right, if there is a user level file, it will take precedence. See ?Startup There was a bug in an older version of R (cannot remember which one) not respecting all the setting. So as always, we also need the R version you are talking about. Best, UWe Ligges Regards [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Multivariate EWMA covariance estimator?
On 02-06-2013, at 15:17, Neuman Co neumanc...@gmail.com wrote:
Hi,
since I want to calculate the VaR of a portfolio consiting of 4 assets
(returns saved into eonreturn,henkelreturn and so on) I have to
estimate the covariance matrix. I do not want to take the rectangular
version with equal weights, but the exponentially weighted moving
average in a multivariate version. I want to estimate a covariance
matrix at every time point t. Then I want to comput the VaR at this
time point t. Afterwards, I will look at the exceedances and do a
backtest.
I tried to implement it as follows (data attached):
lambda-0.9
summe2-0
dummy2-0
covestiexpo-list(NA)
meanvalues-NA
for(i in 101:length(eonreturn)){
meanvalues-matrix(c(mean(eonreturn[(i-100):(i-1)]),mean(henkelreturn[(i-100):(i-1)]),mean(siemensreturn[(i-100):(i-1)]),mean(adidasreturn[(i-100):(i-1)])),4)
for(a in 1:100){
dummy2-lambda^(a-1)*t(datamatrix[(i-a),]-t(meanvalues))%*%(datamatrix[(i-a),]-t(meanvalues))
summe2-summe2+dummy2
}
covestiexpo[[i]]-(1-lambda)*summe2
}
So the covestieexpo[[101]] would be the covariance estimate for the
101th day, taking into account the last 100 observations. Now, the
problem is, that there seems to be something wrong, since the
covariance estimates are cleraly wrong, they seem to be too big. At
the beginning, compared to the normal covariance estimate the
difference is as follows:
covestiexpo[[101]]
[,1][,2][,3][,4]
[1,] 0.004559042 0.002346775 0.004379735 0.003068916
[2,] 0.002346775 0.001978469 0.002536891 0.001909276
[3,] 0.004379735 0.002536891 0.005531590 0.003259803
[4,] 0.003068916 0.001909276 0.003259803 0.003140198
compared to cov(datamatrix[1:100,])
[,1] [,2] [,3][,4]
[1,] 0.0018118239 0.0007432779 0.0015301070 0.001119120
[2,] 0.0007432779 0.0008355960 0.0009281029 0.000754449
[3,] 0.0015301070 0.0009281029 0.0021073171 0.001269626
[4,] 0.0011191199 0.0007544490 0.0012696257 0.001325716
So already here, it is obvious, that something is not correct, if I
look at a period far ahead:
covestiexpo[[1200]]
[,1] [,2] [,3] [,4]
[1,] 0.5312575 0.1939061 0.3419379 0.2475233
[2,] 0.1939061 0.3204951 0.2303478 0.2022423
[3,] 0.3419379 0.2303478 0.5288435 0.2943051
[4,] 0.2475233 0.2022423 0.2943051 0.4599648
you can see, that the values are way too large, so where is my mistake?
Without actual data this is an unverifiable statement.
But you probably have to move the statement
summe2 - 0
to inside the i-forloop just before the a-forloop.
summe2 - 0
for(a in 1:100){
…
so that summe2 is initialized to 0 every time you use it as an accumulator in
the a-forloop.
Furthermore there is no need to initialize dummy2. It gets overwritten
continuously.
Berend
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Regularized Discriminant Analysis scores, anyone?
On 02.06.2013 05:01, Matthew Fagan wrote: Hi all, I am attempting to do Regularized Discriminant Analysis (RDA) on a large dataset, and I want to extract the RDA discriminant score matrix. But the predict function in the klaR package, unlike the predict function for LDA in the MASS package, doesn't seem to give me an option to extract the scores. Any suggestions? There are no such scores: same as for qda, you do not follow the Fisher idea of the linear discriminant components any more: Your space is now partitioned by ellipsoid like structures based on the estimation of the inner-class covariance matrices. rda as implemented in klaR (see the reference given on the help page) is a regularization that helps to overcome problems when estimating non-singular covariance matrices for the separate classes. i have already tried (and failed; ran out of 16 GB of memory) to do this with the rda package: don't know why, but the klaR package seems to be much more efficient with memory. I have included an example below: The rda package provides a completely different regularization technique, see the reference given on the help page. Best, Uwe Ligges library(klaR) library(MASS) data(iris) x - rda(Species ~ ., data = iris, gamma = 0.05, lambda = 0.2) rda1-predict(x, iris[, 1:4]) str(rda1) # This gets you an object with posterior probabilities and classes, but no discriminant scores! # if you run lda y - lda(Species ~ ., data = iris) lda1-predict(y, iris[, 1:4]) str(lda1) head(lda1$x) # gets you the discriminant scores for the LDA. But how to do this for RDA? # curiously, the QDA function in MASS has this same problem, although you can get around it using the rrcov package. Regards, and thank very much for any help, Matt __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] measuring distances between colours?
Sorry about the bug. How embarrassing. Especially because I've learned over
the years to trust my gut feelings when something doesn't feel quite right,
and when I was testing the function, I remember thinking surely there a
better matching named color than 'magenta'.
Thanks for the fix.
Kevin
On Sat, Jun 1, 2013 at 11:30 AM, John Fox j...@mcmaster.ca wrote:
Hi Michael,
This has become a bit of a comedy of errors.
The bug is in Kevin Wright's code, which I adapted, and you too in your
version, which uses local() rather than function() to produce the closure.
The matrix which.col contains character data, as a consequence of binding
the minimum squared distances to colour names, and thus the comparison
cols.near[2,] near^2 doesn't work properly when, ironically, the distance
is small enough so that it's rendered in scientific notation.
Converting to numeric appears to work:
rgb2col2 - local({
+ all.names - colors()
+ all.hsv - rgb2hsv(col2rgb(all.names))
+ find.near - function(x.hsv) {
+ # return the nearest R color name and distance
+ sq.dist - colSums((all.hsv - x.hsv)^2)
+ rbind(all.names[which.min(sq.dist)], min(sq.dist))
+ }
+ function(cols.hex, near=.25){
+ cols.hsv - rgb2hsv(col2rgb(cols.hex))
+ cols.near - apply(cols.hsv, 2, find.near)
+ ifelse(as.numeric(cols.near[2,]) = near^2, cols.near[1,],
cols.hex)
+ }
+ })
rgb2col2(c(#010101, #EE, #AA, #00AA00, #AA,
+ #00, #AA00AA, #00))
[1] black gray93darkred green4blue4
darkgoldenrod
[7] darkmagenta cyan4
The same bug is in the code that I just posted using Lab colours, so (for
posterity) here's a fixed version of that, using local():
rgb2col - local({
+ all.names - colors()
+ all.lab - t(convertColor(t(col2rgb(all.names)), from = sRGB,
+ to = Lab, scale.in = 255))
+ find.near - function(x.lab) {
+ sq.dist - colSums((all.lab - x.lab)^2)
+ rbind(all.names[which.min(sq.dist)], min(sq.dist))
+ }
+ function(cols.hex, near = 2.3) {
+ cols.lab - t(convertColor(t(col2rgb(cols.hex)), from = sRGB,
+ to = Lab, scale.in = 255))
+ cols.near - apply(cols.lab, 2, find.near)
+ ifelse(as.numeric(cols.near[2, ]) near^2, cols.near[1, ],
toupper(cols.hex))
+ }
+ })
rgb2col(c(#010101, #EE, #AA, #00AA00, #AA,
#00, #AA00AA, #00))
[1] black gray93 #AA #00AA00 #AA #00
[7] #AA00AA #00
rgb2col(c(#010101, #EE, #AA, #00AA00, #AA,
#00, #AA00AA, #00), near=15)
[1] black gray93firebrick3limegreen
[5] blue4 #00 darkmagenta lightseagreen
So with Lab colours, setting near to the JND of 2.3 leaves many of these
colours unmatched. I experimented a bit, and using 15 (as above) produces
matches that appear reasonably close to me.
I used squared distances to avoid taking the square-roots of all the
distances. Since the criterion for near colours, which is on the distance
scale, is squared to make the comparison, this shouldn't be problematic.
I hope that finally this will be a satisfactory solution.
Best,
John
-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-bounces@r-
project.org] On Behalf Of Michael Friendly
Sent: Saturday, June 01, 2013 11:33 AM
To: John Fox
Cc: 'r-help'; 'Martin Maechler'
Subject: Re: [R] measuring distances between colours?
Just a quick note: The following two versions of your function don't
give the same results. I'm not sure why, and also not sure why the
criterion for 'near' should be expressed in squared distance.
# version 1
rgb2col - local({
hex2dec - function(hexnums) {
# suggestion of Eik Vettorazzi
sapply(strtoi(hexnums, 16L), function(x) x %/% 256^(2:0) %%
256)
}
findMatch - function(dec.col) {
sq.dist - colSums((hsv - dec.col)^2)
rbind(which.min(sq.dist), min(sq.dist))
}
colors - colors()
hsv - rgb2hsv(col2rgb(colors))
function(cols, near=0.25) {
cols - sub(^#, , toupper(cols))
dec.cols - rgb2hsv(hex2dec(cols))
which.col - apply(dec.cols, 2, findMatch)
matches - colors[which.col[1, ]]
unmatched - which.col[2, ] near^2
matches[unmatched] - paste(#, cols[unmatched], sep=)
matches
}
})
# version 2
rgb2col2 - local({
all.names - colors()
all.hsv - rgb2hsv(col2rgb(all.names))
find.near - function(x.hsv) {
# return the nearest R color name and distance
sq.dist - colSums((all.hsv - x.hsv)^2)
rbind(all.names[which.min(sq.dist)], min(sq.dist))
}
function(cols.hex, near=.25){
cols.hsv - rgb2hsv(col2rgb(cols.hex))
cols.near
Re: [R] Official way to set/retrieve options in packages?
What would be an example of setting, saving, and re-loading an option to a
user's .Rprofile -- and would this be a no-no in a CRAN package?
--j
On Sat, Jun 1, 2013 at 4:57 PM, Prof Brian Ripley rip...@stats.ox.ac.ukwrote:
On 01/06/2013 22:44, Anthony Damico wrote:
hope this helps.. :)
# define an object `x`
x - list( any value here , 10 )
# set `myoption` to that object
options( myoption = x )
# retrieve it later (perhaps within a function elsewhere in the
package)
( y - getOption( myoption ) )
it's nice to name your options `mypackage.myoption` so users know what
package the option is associated with in case they type `options()`
here's the `.onLoad` function in the R survey package. notice how the
options are only set *if* they don't already exist--
But a nicer convention is that used by most packages in R itself: if the
option is not set, the function using it assumes a suitable default.
That would make sense for all the FALSE defaults below.
Note though that this is not 'persistent': users have to set options in
their startup files (see ?Startup). There is no official location to
store package configurations. Users generally dislike software saving
settings in their own file space so it seems very much preferable to use
the standard R mechanisms (.Rprofile etc).
survey:::.onLoad
function (...)
{
if (is.null(getOption(survey.lonely.psu)))
options(survey.lonely.psu = fail)
if (is.null(getOption(survey.ultimate.cluster)))
options(survey.ultimate.cluster = FALSE)
if (is.null(getOption(survey.want.obsolete)))
options(survey.want.obsolete = FALSE)
if (is.null(getOption(survey.adjust.domain.lonely)))
options(survey.adjust.domain.lonely = FALSE)
if (is.null(getOption(survey.drop.replicates)))
options(survey.drop.replicates = TRUE)
if (is.null(getOption(survey.multicore)))
options(survey.multicore = FALSE)
if (is.null(getOption(survey.replicates.mse)))
options(survey.replicates.mse = FALSE)
}
environment: namespace:survey
On Sat, Jun 1, 2013 at 4:01 PM, Jonathan Greenberg j...@illinois.edu
wrote:
R-helpers:
Say I'm developing a package that has a set of user-definable options
that
I would like to be persistent across R-invocations (they are saved
someplace). Of course, I can create a little text file to be
written/read,
but I was wondering if there is an officially sanctioned way to do
this?
I see there is an options() and getOptions() function, but I'm
unclear how
I would use this in my own package to create/save new options for my
particular package. Cheers!
--j
--
Jonathan A. Greenberg, PhD
Assistant Professor
Global Environmental Analysis and Remote Sensing (GEARS) Laboratory
Department of Geography and Geographic Information Science
University of Illinois at Urbana-Champaign
607 South Mathews Avenue, MC 150
Urbana, IL 61801
Phone: 217-300-1924
http://www.geog.illinois.edu/~jgrn/
AIM: jgrn307, MSN: jgrn...@hotmail.com, Gchat: jgrn307, Skype: jgrn3007
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
--
Brian D. Ripley, rip...@stats.ox.ac.uk
Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel: +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UKFax: +44 1865 272595
--
Jonathan A. Greenberg, PhD
Assistant Professor
Global Environmental Analysis and Remote Sensing (GEARS) Laboratory
Department of Geography and Geographic Information Science
University of Illinois at Urbana-Champaign
607 South Mathews Avenue, MC 150
Urbana, IL 61801
Phone: 217-300-1924
http://www.geog.illinois.edu/~jgrn/
AIM: jgrn307, MSN: jgrn...@hotmail.com, Gchat: jgrn307, Skype: jgrn3007
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Official way to set/retrieve options in packages?
Have you read
?Startup
which explains startup procedures, including reading of .Rprofile.
I cannot speak for Brian Ripley, but I can tell you that I would
prefer that packages did not mess with my .Rprofile files.
OTOH, I have no objections if packages suggest that I set certain
options in them or even give me lines of code that I can add if I
choose ( or give me an option to call a function to do so).
I understand that this is not the automatic procedure you seem to be seeking...
Cheers,
Bert
On Sun, Jun 2, 2013 at 9:06 AM, Jonathan Greenberg j...@illinois.edu wrote:
What would be an example of setting, saving, and re-loading an option to a
user's .Rprofile -- and would this be a no-no in a CRAN package?
--j
On Sat, Jun 1, 2013 at 4:57 PM, Prof Brian Ripley
rip...@stats.ox.ac.ukwrote:
On 01/06/2013 22:44, Anthony Damico wrote:
hope this helps.. :)
# define an object `x`
x - list( any value here , 10 )
# set `myoption` to that object
options( myoption = x )
# retrieve it later (perhaps within a function elsewhere in the
package)
( y - getOption( myoption ) )
it's nice to name your options `mypackage.myoption` so users know what
package the option is associated with in case they type `options()`
here's the `.onLoad` function in the R survey package. notice how the
options are only set *if* they don't already exist--
But a nicer convention is that used by most packages in R itself: if the
option is not set, the function using it assumes a suitable default.
That would make sense for all the FALSE defaults below.
Note though that this is not 'persistent': users have to set options in
their startup files (see ?Startup). There is no official location to
store package configurations. Users generally dislike software saving
settings in their own file space so it seems very much preferable to use
the standard R mechanisms (.Rprofile etc).
survey:::.onLoad
function (...)
{
if (is.null(getOption(survey.lonely.psu)))
options(survey.lonely.psu = fail)
if (is.null(getOption(survey.ultimate.cluster)))
options(survey.ultimate.cluster = FALSE)
if (is.null(getOption(survey.want.obsolete)))
options(survey.want.obsolete = FALSE)
if (is.null(getOption(survey.adjust.domain.lonely)))
options(survey.adjust.domain.lonely = FALSE)
if (is.null(getOption(survey.drop.replicates)))
options(survey.drop.replicates = TRUE)
if (is.null(getOption(survey.multicore)))
options(survey.multicore = FALSE)
if (is.null(getOption(survey.replicates.mse)))
options(survey.replicates.mse = FALSE)
}
environment: namespace:survey
On Sat, Jun 1, 2013 at 4:01 PM, Jonathan Greenberg j...@illinois.edu
wrote:
R-helpers:
Say I'm developing a package that has a set of user-definable options
that
I would like to be persistent across R-invocations (they are saved
someplace). Of course, I can create a little text file to be
written/read,
but I was wondering if there is an officially sanctioned way to do
this?
I see there is an options() and getOptions() function, but I'm
unclear how
I would use this in my own package to create/save new options for my
particular package. Cheers!
--j
--
Jonathan A. Greenberg, PhD
Assistant Professor
Global Environmental Analysis and Remote Sensing (GEARS) Laboratory
Department of Geography and Geographic Information Science
University of Illinois at Urbana-Champaign
607 South Mathews Avenue, MC 150
Urbana, IL 61801
Phone: 217-300-1924
http://www.geog.illinois.edu/~jgrn/
AIM: jgrn307, MSN: jgrn...@hotmail.com, Gchat: jgrn307, Skype: jgrn3007
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
--
Brian D. Ripley, rip...@stats.ox.ac.uk
Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel: +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UKFax: +44 1865 272595
--
Jonathan A. Greenberg, PhD
Assistant Professor
Global Environmental Analysis and Remote Sensing (GEARS) Laboratory
Department of Geography and Geographic Information Science
University of Illinois at Urbana-Champaign
607 South Mathews Avenue, MC 150
Urbana, IL 61801
Phone:
Re: [R] Multivariate EWMA covariance estimator?
Thanks a lot for your answer, one more question:
I now use 100 values, so not infinity values. That means I cut some
values off, so the weights will not sum up to one. With which factor
do I have to multiply the (1-lambda)*summe2 to rescale it? So that I
do not always underestimate the variance anymore?
2013/6/2 Berend Hasselman b...@xs4all.nl:
On 02-06-2013, at 15:17, Neuman Co neumanc...@gmail.com wrote:
Hi,
since I want to calculate the VaR of a portfolio consiting of 4 assets
(returns saved into eonreturn,henkelreturn and so on) I have to
estimate the covariance matrix. I do not want to take the rectangular
version with equal weights, but the exponentially weighted moving
average in a multivariate version. I want to estimate a covariance
matrix at every time point t. Then I want to comput the VaR at this
time point t. Afterwards, I will look at the exceedances and do a
backtest.
I tried to implement it as follows (data attached):
lambda-0.9
summe2-0
dummy2-0
covestiexpo-list(NA)
meanvalues-NA
for(i in 101:length(eonreturn)){
meanvalues-matrix(c(mean(eonreturn[(i-100):(i-1)]),mean(henkelreturn[(i-100):(i-1)]),mean(siemensreturn[(i-100):(i-1)]),mean(adidasreturn[(i-100):(i-1)])),4)
for(a in 1:100){
dummy2-lambda^(a-1)*t(datamatrix[(i-a),]-t(meanvalues))%*%(datamatrix[(i-a),]-t(meanvalues))
summe2-summe2+dummy2
}
covestiexpo[[i]]-(1-lambda)*summe2
}
So the covestieexpo[[101]] would be the covariance estimate for the
101th day, taking into account the last 100 observations. Now, the
problem is, that there seems to be something wrong, since the
covariance estimates are cleraly wrong, they seem to be too big. At
the beginning, compared to the normal covariance estimate the
difference is as follows:
covestiexpo[[101]]
[,1][,2][,3][,4]
[1,] 0.004559042 0.002346775 0.004379735 0.003068916
[2,] 0.002346775 0.001978469 0.002536891 0.001909276
[3,] 0.004379735 0.002536891 0.005531590 0.003259803
[4,] 0.003068916 0.001909276 0.003259803 0.003140198
compared to cov(datamatrix[1:100,])
[,1] [,2] [,3][,4]
[1,] 0.0018118239 0.0007432779 0.0015301070 0.001119120
[2,] 0.0007432779 0.0008355960 0.0009281029 0.000754449
[3,] 0.0015301070 0.0009281029 0.0021073171 0.001269626
[4,] 0.0011191199 0.0007544490 0.0012696257 0.001325716
So already here, it is obvious, that something is not correct, if I
look at a period far ahead:
covestiexpo[[1200]]
[,1] [,2] [,3] [,4]
[1,] 0.5312575 0.1939061 0.3419379 0.2475233
[2,] 0.1939061 0.3204951 0.2303478 0.2022423
[3,] 0.3419379 0.2303478 0.5288435 0.2943051
[4,] 0.2475233 0.2022423 0.2943051 0.4599648
you can see, that the values are way too large, so where is my mistake?
Without actual data this is an unverifiable statement.
But you probably have to move the statement
summe2 - 0
to inside the i-forloop just before the a-forloop.
summe2 - 0
for(a in 1:100){
…
so that summe2 is initialized to 0 every time you use it as an accumulator in
the a-forloop.
Furthermore there is no need to initialize dummy2. It gets overwritten
continuously.
Berend
--
Neumann, Conrad
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Multivariate EWMA covariance estimator?
On 02-06-2013, at 19:03, Neuman Co neumanc...@gmail.com wrote:
Thanks a lot for your answer, one more question:
I now use 100 values, so not infinity values. That means I cut some
values off, so the weights will not sum up to one. With which factor
do I have to multiply the (1-lambda)*summe2 to rescale it? So that I
do not always underestimate the variance anymore?
I don't know but maybe something like this
1/sum(lambda^((1:100)-1))/(1-lambda)
which in your case is 1.27
Berend
2013/6/2 Berend Hasselman b...@xs4all.nl:
On 02-06-2013, at 15:17, Neuman Co neumanc...@gmail.com wrote:
Hi,
since I want to calculate the VaR of a portfolio consiting of 4 assets
(returns saved into eonreturn,henkelreturn and so on) I have to
estimate the covariance matrix. I do not want to take the rectangular
version with equal weights, but the exponentially weighted moving
average in a multivariate version. I want to estimate a covariance
matrix at every time point t. Then I want to comput the VaR at this
time point t. Afterwards, I will look at the exceedances and do a
backtest.
I tried to implement it as follows (data attached):
lambda-0.9
summe2-0
dummy2-0
covestiexpo-list(NA)
meanvalues-NA
for(i in 101:length(eonreturn)){
meanvalues-matrix(c(mean(eonreturn[(i-100):(i-1)]),mean(henkelreturn[(i-100):(i-1)]),mean(siemensreturn[(i-100):(i-1)]),mean(adidasreturn[(i-100):(i-1)])),4)
for(a in 1:100){
dummy2-lambda^(a-1)*t(datamatrix[(i-a),]-t(meanvalues))%*%(datamatrix[(i-a),]-t(meanvalues))
summe2-summe2+dummy2
}
covestiexpo[[i]]-(1-lambda)*summe2
}
So the covestieexpo[[101]] would be the covariance estimate for the
101th day, taking into account the last 100 observations. Now, the
problem is, that there seems to be something wrong, since the
covariance estimates are cleraly wrong, they seem to be too big. At
the beginning, compared to the normal covariance estimate the
difference is as follows:
covestiexpo[[101]]
[,1][,2][,3][,4]
[1,] 0.004559042 0.002346775 0.004379735 0.003068916
[2,] 0.002346775 0.001978469 0.002536891 0.001909276
[3,] 0.004379735 0.002536891 0.005531590 0.003259803
[4,] 0.003068916 0.001909276 0.003259803 0.003140198
compared to cov(datamatrix[1:100,])
[,1] [,2] [,3][,4]
[1,] 0.0018118239 0.0007432779 0.0015301070 0.001119120
[2,] 0.0007432779 0.0008355960 0.0009281029 0.000754449
[3,] 0.0015301070 0.0009281029 0.0021073171 0.001269626
[4,] 0.0011191199 0.0007544490 0.0012696257 0.001325716
So already here, it is obvious, that something is not correct, if I
look at a period far ahead:
covestiexpo[[1200]]
[,1] [,2] [,3] [,4]
[1,] 0.5312575 0.1939061 0.3419379 0.2475233
[2,] 0.1939061 0.3204951 0.2303478 0.2022423
[3,] 0.3419379 0.2303478 0.5288435 0.2943051
[4,] 0.2475233 0.2022423 0.2943051 0.4599648
you can see, that the values are way too large, so where is my mistake?
Without actual data this is an unverifiable statement.
But you probably have to move the statement
summe2 - 0
to inside the i-forloop just before the a-forloop.
summe2 - 0
for(a in 1:100){
…
so that summe2 is initialized to 0 every time you use it as an accumulator
in the a-forloop.
Furthermore there is no need to initialize dummy2. It gets overwritten
continuously.
Berend
--
Neumann, Conrad
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Multivariate EWMA covariance estimator?
Again, a big thanks for your answer.
On this webpage:
http://financetrainingcourse.com/education/2010/03/master-class-calculating-value-at-risk-var-final-steps/
I found, that I have to rescale by dividing the weights calculated in
Step B2 by 1-?n
The ? is the lambda, since the webpage cannot display it, I also
found it on another webpage, therefore, I changed my code to the
following:
dummy2-lambda^(a-1)/(1-lambda^100)*t(datamatrix[(i-a),]-t(meanvalues))%*%(datamatrix[(i-a),]-t(meanvalues))
Do you think this is correct?
One further question: You also told me, that I do not have to
initialize my dummy2, what does this mean? I wrote dummy2-0 because I
have to create this variable before using it for the loop?
2013/6/2 Berend Hasselman b...@xs4all.nl:
On 02-06-2013, at 19:03, Neuman Co neumanc...@gmail.com wrote:
Thanks a lot for your answer, one more question:
I now use 100 values, so not infinity values. That means I cut some
values off, so the weights will not sum up to one. With which factor
do I have to multiply the (1-lambda)*summe2 to rescale it? So that I
do not always underestimate the variance anymore?
I don't know but maybe something like this
1/sum(lambda^((1:100)-1))/(1-lambda)
which in your case is 1.27
Berend
2013/6/2 Berend Hasselman b...@xs4all.nl:
On 02-06-2013, at 15:17, Neuman Co neumanc...@gmail.com wrote:
Hi,
since I want to calculate the VaR of a portfolio consiting of 4 assets
(returns saved into eonreturn,henkelreturn and so on) I have to
estimate the covariance matrix. I do not want to take the rectangular
version with equal weights, but the exponentially weighted moving
average in a multivariate version. I want to estimate a covariance
matrix at every time point t. Then I want to comput the VaR at this
time point t. Afterwards, I will look at the exceedances and do a
backtest.
I tried to implement it as follows (data attached):
lambda-0.9
summe2-0
dummy2-0
covestiexpo-list(NA)
meanvalues-NA
for(i in 101:length(eonreturn)){
meanvalues-matrix(c(mean(eonreturn[(i-100):(i-1)]),mean(henkelreturn[(i-100):(i-1)]),mean(siemensreturn[(i-100):(i-1)]),mean(adidasreturn[(i-100):(i-1)])),4)
for(a in 1:100){
dummy2-lambda^(a-1)*t(datamatrix[(i-a),]-t(meanvalues))%*%(datamatrix[(i-a),]-t(meanvalues))
summe2-summe2+dummy2
}
covestiexpo[[i]]-(1-lambda)*summe2
}
So the covestieexpo[[101]] would be the covariance estimate for the
101th day, taking into account the last 100 observations. Now, the
problem is, that there seems to be something wrong, since the
covariance estimates are cleraly wrong, they seem to be too big. At
the beginning, compared to the normal covariance estimate the
difference is as follows:
covestiexpo[[101]]
[,1][,2][,3][,4]
[1,] 0.004559042 0.002346775 0.004379735 0.003068916
[2,] 0.002346775 0.001978469 0.002536891 0.001909276
[3,] 0.004379735 0.002536891 0.005531590 0.003259803
[4,] 0.003068916 0.001909276 0.003259803 0.003140198
compared to cov(datamatrix[1:100,])
[,1] [,2] [,3][,4]
[1,] 0.0018118239 0.0007432779 0.0015301070 0.001119120
[2,] 0.0007432779 0.0008355960 0.0009281029 0.000754449
[3,] 0.0015301070 0.0009281029 0.0021073171 0.001269626
[4,] 0.0011191199 0.0007544490 0.0012696257 0.001325716
So already here, it is obvious, that something is not correct, if I
look at a period far ahead:
covestiexpo[[1200]]
[,1] [,2] [,3] [,4]
[1,] 0.5312575 0.1939061 0.3419379 0.2475233
[2,] 0.1939061 0.3204951 0.2303478 0.2022423
[3,] 0.3419379 0.2303478 0.5288435 0.2943051
[4,] 0.2475233 0.2022423 0.2943051 0.4599648
you can see, that the values are way too large, so where is my mistake?
Without actual data this is an unverifiable statement.
But you probably have to move the statement
summe2 - 0
to inside the i-forloop just before the a-forloop.
summe2 - 0
for(a in 1:100){
…
so that summe2 is initialized to 0 every time you use it as an accumulator
in the a-forloop.
Furthermore there is no need to initialize dummy2. It gets overwritten
continuously.
Berend
--
Neumann, Conrad
--
Neumann, Conrad
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Multivariate EWMA covariance estimator?
On 02-06-2013, at 19:45, Neuman Co neumanc...@gmail.com wrote:
Again, a big thanks for your answer.
On this webpage:
http://financetrainingcourse.com/education/2010/03/master-class-calculating-value-at-risk-var-final-steps/
I found, that I have to rescale by dividing the weights calculated in
Step B2 by 1-?n
The ? is the lambda, since the webpage cannot display it, I also
found it on another webpage, therefore, I changed my code to the
following:
dummy2-lambda^(a-1)/(1-lambda^100)*t(datamatrix[(i-a),]-t(meanvalues))%*%(datamatrix[(i-a),]-t(meanvalues))
Do you think this is correct?
I'm not going to do homework.
You have to scale the weights. So you need to scale lambda^(a-1) by their sum
for a=1:100 if I understand correctly.
One further question: You also told me, that I do not have to
initialize my dummy2, what does this mean?
Exactly that: you don't have to initialize dummy2.
You are assigning a value to dummy2 for each value of a.
You are thus either creating a new object dummy2 or overwriting any existing
object dummy2 and thus destroying a previous value.
I wrote dummy2-0 because I
have to create this variable before using it for the loop?
You are not using it on the righthand side of an expression. You are assigning
to it.
You only need to initialize dummy2 if you have an expression involving dummy2
on the righthand side e.g. dummy2 - dummy2 + ……
Berend
2013/6/2 Berend Hasselman b...@xs4all.nl:
On 02-06-2013, at 19:03, Neuman Co neumanc...@gmail.com wrote:
Thanks a lot for your answer, one more question:
I now use 100 values, so not infinity values. That means I cut some
values off, so the weights will not sum up to one. With which factor
do I have to multiply the (1-lambda)*summe2 to rescale it? So that I
do not always underestimate the variance anymore?
I don't know but maybe something like this
1/sum(lambda^((1:100)-1))/(1-lambda)
which in your case is 1.27
Berend
2013/6/2 Berend Hasselman b...@xs4all.nl:
On 02-06-2013, at 15:17, Neuman Co neumanc...@gmail.com wrote:
Hi,
since I want to calculate the VaR of a portfolio consiting of 4 assets
(returns saved into eonreturn,henkelreturn and so on) I have to
estimate the covariance matrix. I do not want to take the rectangular
version with equal weights, but the exponentially weighted moving
average in a multivariate version. I want to estimate a covariance
matrix at every time point t. Then I want to comput the VaR at this
time point t. Afterwards, I will look at the exceedances and do a
backtest.
I tried to implement it as follows (data attached):
lambda-0.9
summe2-0
dummy2-0
covestiexpo-list(NA)
meanvalues-NA
for(i in 101:length(eonreturn)){
meanvalues-matrix(c(mean(eonreturn[(i-100):(i-1)]),mean(henkelreturn[(i-100):(i-1)]),mean(siemensreturn[(i-100):(i-1)]),mean(adidasreturn[(i-100):(i-1)])),4)
for(a in 1:100){
dummy2-lambda^(a-1)*t(datamatrix[(i-a),]-t(meanvalues))%*%(datamatrix[(i-a),]-t(meanvalues))
summe2-summe2+dummy2
}
covestiexpo[[i]]-(1-lambda)*summe2
}
So the covestieexpo[[101]] would be the covariance estimate for the
101th day, taking into account the last 100 observations. Now, the
problem is, that there seems to be something wrong, since the
covariance estimates are cleraly wrong, they seem to be too big. At
the beginning, compared to the normal covariance estimate the
difference is as follows:
covestiexpo[[101]]
[,1][,2][,3][,4]
[1,] 0.004559042 0.002346775 0.004379735 0.003068916
[2,] 0.002346775 0.001978469 0.002536891 0.001909276
[3,] 0.004379735 0.002536891 0.005531590 0.003259803
[4,] 0.003068916 0.001909276 0.003259803 0.003140198
compared to cov(datamatrix[1:100,])
[,1] [,2] [,3][,4]
[1,] 0.0018118239 0.0007432779 0.0015301070 0.001119120
[2,] 0.0007432779 0.0008355960 0.0009281029 0.000754449
[3,] 0.0015301070 0.0009281029 0.0021073171 0.001269626
[4,] 0.0011191199 0.0007544490 0.0012696257 0.001325716
So already here, it is obvious, that something is not correct, if I
look at a period far ahead:
covestiexpo[[1200]]
[,1] [,2] [,3] [,4]
[1,] 0.5312575 0.1939061 0.3419379 0.2475233
[2,] 0.1939061 0.3204951 0.2303478 0.2022423
[3,] 0.3419379 0.2303478 0.5288435 0.2943051
[4,] 0.2475233 0.2022423 0.2943051 0.4599648
you can see, that the values are way too large, so where is my mistake?
Without actual data this is an unverifiable statement.
But you probably have to move the statement
summe2 - 0
to inside the i-forloop just before the a-forloop.
summe2 - 0
for(a in 1:100){
…
so that summe2 is initialized to 0 every time you use it as an accumulator
in the a-forloop.
Furthermore there is no need to initialize dummy2. It gets overwritten
continuously.
Berend
--
Neumann, Conrad
--
Neumann, Conrad
Re: [R] measuring distances between colours?
Dear Kevin,
When computer code is bug free, we'll probably all be out of business. Thank
you for improving my original code.
Best,
John
-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-bounces@r-
project.org] On Behalf Of Kevin Wright
Sent: Sunday, June 02, 2013 10:43 AM
To: John Fox
Cc: r-help; Michael Friendly; Martin Maechler
Subject: Re: [R] measuring distances between colours?
Sorry about the bug. How embarrassing. Especially because I've learned
over
the years to trust my gut feelings when something doesn't feel quite
right,
and when I was testing the function, I remember thinking surely there
a
better matching named color than 'magenta'.
Thanks for the fix.
Kevin
On Sat, Jun 1, 2013 at 11:30 AM, John Fox j...@mcmaster.ca wrote:
Hi Michael,
This has become a bit of a comedy of errors.
The bug is in Kevin Wright's code, which I adapted, and you too in
your
version, which uses local() rather than function() to produce the
closure.
The matrix which.col contains character data, as a consequence of
binding
the minimum squared distances to colour names, and thus the
comparison
cols.near[2,] near^2 doesn't work properly when, ironically, the
distance
is small enough so that it's rendered in scientific notation.
Converting to numeric appears to work:
rgb2col2 - local({
+ all.names - colors()
+ all.hsv - rgb2hsv(col2rgb(all.names))
+ find.near - function(x.hsv) {
+ # return the nearest R color name and distance
+ sq.dist - colSums((all.hsv - x.hsv)^2)
+ rbind(all.names[which.min(sq.dist)], min(sq.dist))
+ }
+ function(cols.hex, near=.25){
+ cols.hsv - rgb2hsv(col2rgb(cols.hex))
+ cols.near - apply(cols.hsv, 2, find.near)
+ ifelse(as.numeric(cols.near[2,]) = near^2, cols.near[1,],
cols.hex)
+ }
+ })
rgb2col2(c(#010101, #EE, #AA, #00AA00, #AA,
+ #00, #AA00AA, #00))
[1] black gray93darkred green4
blue4
darkgoldenrod
[7] darkmagenta cyan4
The same bug is in the code that I just posted using Lab colours, so
(for
posterity) here's a fixed version of that, using local():
rgb2col - local({
+ all.names - colors()
+ all.lab - t(convertColor(t(col2rgb(all.names)), from = sRGB,
+ to = Lab, scale.in = 255))
+ find.near - function(x.lab) {
+ sq.dist - colSums((all.lab - x.lab)^2)
+ rbind(all.names[which.min(sq.dist)], min(sq.dist))
+ }
+ function(cols.hex, near = 2.3) {
+ cols.lab - t(convertColor(t(col2rgb(cols.hex)), from =
sRGB,
+ to = Lab, scale.in = 255))
+ cols.near - apply(cols.lab, 2, find.near)
+ ifelse(as.numeric(cols.near[2, ]) near^2, cols.near[1, ],
toupper(cols.hex))
+ }
+ })
rgb2col(c(#010101, #EE, #AA, #00AA00, #AA,
#00, #AA00AA, #00))
[1] black gray93 #AA #00AA00 #AA #00
[7] #AA00AA #00
rgb2col(c(#010101, #EE, #AA, #00AA00, #AA,
#00, #AA00AA, #00), near=15)
[1] black gray93firebrick3limegreen
[5] blue4 #00 darkmagenta lightseagreen
So with Lab colours, setting near to the JND of 2.3 leaves many of
these
colours unmatched. I experimented a bit, and using 15 (as above)
produces
matches that appear reasonably close to me.
I used squared distances to avoid taking the square-roots of all the
distances. Since the criterion for near colours, which is on the
distance
scale, is squared to make the comparison, this shouldn't be
problematic.
I hope that finally this will be a satisfactory solution.
Best,
John
-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-bounces@r-
project.org] On Behalf Of Michael Friendly
Sent: Saturday, June 01, 2013 11:33 AM
To: John Fox
Cc: 'r-help'; 'Martin Maechler'
Subject: Re: [R] measuring distances between colours?
Just a quick note: The following two versions of your function
don't
give the same results. I'm not sure why, and also not sure why the
criterion for 'near' should be expressed in squared distance.
# version 1
rgb2col - local({
hex2dec - function(hexnums) {
# suggestion of Eik Vettorazzi
sapply(strtoi(hexnums, 16L), function(x) x %/% 256^(2:0)
%%
256)
}
findMatch - function(dec.col) {
sq.dist - colSums((hsv - dec.col)^2)
rbind(which.min(sq.dist), min(sq.dist))
}
colors - colors()
hsv - rgb2hsv(col2rgb(colors))
function(cols, near=0.25) {
cols - sub(^#, , toupper(cols))
dec.cols - rgb2hsv(hex2dec(cols))
which.col - apply(dec.cols, 2, findMatch)
matches - colors[which.col[1, ]]
unmatched -
[R] useR meetup group in Munich, Friday 7th: Web Service Frameworks with R
Dear all, I would like to invite Munich (Germany) area R users for our second meeting: 7th June 2013. The group is aimed to bring together practitioners from industry and academia in order to exchange knowledge and experience in solving data analysis statistical problems by using R. More information about the group at: http://www.meetup.com/munich-useR-group/ Our second meeting will host two talks about Web Service Frameworks with R (shiny and RevolutionDeployR). Meet you in Munich Markus __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Scale package change comma defaults?
I was using the comma() funtion in the scales page and was wondering how to chage a the defaults. comm(1000) gives me 1,000 which is what I usually want but how would I change the output to 1.000. I had thought that I could simply do comm(1000, big.mark = .) but I am getting Error in format.default(x, ..., big.mark = ,, scientific = FALSE, trim = TRUE) : formal argument big.mark matched by multiple actual arguments And since I'm here I might as well ask if there is a way to keep a couple fo decemal points rather than rounding to the first integer. Thanks John Kane Kingston ON Canada sessionInfo() R version 3.0.1 (2013-05-16) Platform: i686-pc-linux-gnu (32-bit) locale: [1] LC_CTYPE=en_CA.UTF-8 LC_NUMERIC=C LC_TIME=en_CA.UTF-8 [4] LC_COLLATE=en_CA.UTF-8 LC_MONETARY=en_CA.UTF-8 LC_MESSAGES=en_CA.UTF-8 [7] LC_PAPER=C LC_NAME=C LC_ADDRESS=C [10] LC_TELEPHONE=C LC_MEASUREMENT=en_CA.UTF-8 LC_IDENTIFICATION=C attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] scales_0.2.3 FREE 3D MARINE AQUARIUM SCREENSAVER - Watch dolphins, sharks orcas on your desktop! __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Scale package change comma defaults?
John Kane jrkrideau at inbox.com writes:
I was using the comma() funtion in the scales page and was
wondering how to chage a the defaults. comm(1000) gives me 1,000
which is what I usually want but how would I change the output to
1.000.
Since the comma() function is just a wrapper for format():
scales::comma
function (x, ...)
{
format(x, ..., big.mark = ,, scientific = FALSE, trim = TRUE)
}
Why not just define your own function that does what you want?
myComma - function (x, big.mark=,, ...)
{
format(x, ..., big.mark = big.mark, scientific = FALSE, trim = TRUE)
}
myComma(1000,.)
I had thought that I could simply do
comm(1000, big.mark = .)
but I am getting
Error in format.default(x, ..., big.mark = ,, scientific = FALSE, trim =
TRUE) :
formal argument big.mark matched by multiple actual arguments
And since I'm here I might as well ask if there is a way
to keep a couple fo decemal points rather than rounding to the first integer.
Not quite sure what you mean here; see ?format for more details.
format(1200,big.mark=.,nsmall=2)
might be what you want, but if you're going to use . for big.mark
then you might want:
options(OutDec=,)
[1] 1.200,00
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] X-axis in meandist {vegan}
I am using vegan package for estimating dissimilarity within and between clusters, and plotting them. The following code gives me a dendrogram without X-axis. How would I add X axis to the plot? Below dataEnv is a matrix of my environmental variables. dataFact$Source_by_Type has the levels that delineate clusters. dataEnvDist - vegdist(dataEnv, method=bray) ## bray curtis dissimilarity meanDist - meandist(dataEnvDist, dataFact$Source_by_Type) plot(meanDist, ylab=Mean dissimilarity) Thanks in advance. Kumar Mainali -- Section of Integrative Biology University of Texas at Austin Austin, Texas 78712, USA [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Regularized Discriminant Analysis scores, anyone?
Thank you Dr. Ligges, i very much appreciate the quick reply. i wondered if that was the case, based on the math as I (poorly) understood it. However i remain confused. page 107 from the rrcov package PDF makes me think I can derive LDA-style discriminant scores for a QDA: library(rrcov) data(iris) qda1-QdaClassic(x=iris[,1:4], grouping=iris[,5]) pred_qda-predict(qda1, iris[,1:4]) head(pred_qda@x) plotdat-pred_qda@x plot(plotdat[,1], plotdat[,2]) plot(plotdat[,2], plotdat[,3]) pred_qda$x looks like QDA discriminant scores. No doubt you are right, but if you have a moment, I'd love to know what these scores are and what they summarize. In addition, I have run into this nice set of lengthy R code to manually calculate discriminant scores for a QDA: https://cs.uwaterloo.ca/~a2curtis/courses/2005/ML-classification.pdf None of this means i can calculate discriminant scores for a RDA, of course, but QDA is my back-up choice. Bottom line: am i am completely misinterpreting what I am seeing here, mathematically? Or is this just the result of different ways of implementing QDA in R? Regards, and thanks again, Matt On 6/2/2013 10:39 AM, Uwe Ligges wrote: On 02.06.2013 05:01, Matthew Fagan wrote: Hi all, I am attempting to do Regularized Discriminant Analysis (RDA) on a large dataset, and I want to extract the RDA discriminant score matrix. But the predict function in the klaR package, unlike the predict function for LDA in the MASS package, doesn't seem to give me an option to extract the scores. Any suggestions? There are no such scores: same as for qda, you do not follow the Fisher idea of the linear discriminant components any more: Your space is now partitioned by ellipsoid like structures based on the estimation of the inner-class covariance matrices. rda as implemented in klaR (see the reference given on the help page) is a regularization that helps to overcome problems when estimating non-singular covariance matrices for the separate classes. i have already tried (and failed; ran out of 16 GB of memory) to do this with the rda package: don't know why, but the klaR package seems to be much more efficient with memory. I have included an example below: The rda package provides a completely different regularization technique, see the reference given on the help page. Best, Uwe Ligges library(klaR) library(MASS) data(iris) x - rda(Species ~ ., data = iris, gamma = 0.05, lambda = 0.2) rda1-predict(x, iris[, 1:4]) str(rda1) # This gets you an object with posterior probabilities and classes, but no discriminant scores! # if you run lda y - lda(Species ~ ., data = iris) lda1-predict(y, iris[, 1:4]) str(lda1) head(lda1$x) # gets you the discriminant scores for the LDA. But how to do this for RDA? # curiously, the QDA function in MASS has this same problem, although you can get around it using the rrcov package. Regards, and thank very much for any help, Matt __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Matthew Fagan Columbia University Department of Ecology, Evolution, and Environmental Biology 512-569-1417 (cell/home) (212) 854-9987 (office) (212) 854-8188 (fax) __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Problem with package 'mi' for multiple imputation
Dear R-Community, I have just started using R (3.0.0) for my statistics and I find it the best programm for running my regressions. However, as I am still a junior user I soon came to the limit of my technical knowhow. Here's my situation: My dataset consists of 30 observations and 100 variables (numeric and factors). Some of the observations have missing values and as a result these observations will be omitted when running the regressions, which in turn reduces the sample. In order to counteract this, I have installed the package 'mi' for applying the method of multiple imputation. To reduce the structural problems in the data, I have corrected the variable types (i.e. numeric to factors and vice versa) and eliminated one observation, as it was an outlier and led to collinearity amongst some variables. To control for the problem of contrasts can be applied only to factors with 2 or more levels all factorvariables with only 1 level where excluded. Here's the information matrix: names include order number.mis all.mis type collinear 1 InvestmentYear Yes 1 1 No positive-continuous No 2 InvestmentSize Yes NA 0 No positive-continuous No 3 InvestmentSyndication Yes NA 0 No binary No 4 InvestmentRole Yes NA 0 No binary No 5 InvestmentExitYearAgreed Yes 2 7 No nonnegative No 6 InvestmentExitYearActual Yes 3 3 No nonnegative No 7 InvestmentExitStrategyAgreed Yes 4 5 No unordered-categorical No 8 InvestmentExitStrategyActual Yes 5 3 No unordered-categorical No 9 InvestmentAmountRL Yes 6 7 No continuous No 10 Sector Yes NA 0 No unordered-categorical No 11 BusinessStage Yes NA 0 No unordered-categorical No 12 Origin Yes 7 2 No unordered-categorical No 13 ManagementFoundingMember Yes NA 0 No binary No 14 ManagementFormerPC Yes NA 0 No binary No 15 ManagementFormerPCQuantity Yes NA 0 No ordered-categorical No 16 ManagementFormerPCSuccess Yes 8 1 No ordered-categorical No 17 ManagementCompositionCEO Yes 9 2 No ordered-categorical No 18 ManagementCompositionTechnical Yes 10 2 No ordered-categorical No 19 ManagementCompositionFinancial Yes 11 2 No binary No 20 ManagementCompositionScientific Yes 12 2 No ordered-categorical No 21 ManagementCompositionMarketing Yes 13 2 No fixed No 22 ManagementCompositionSales Yes 14 2 No fixed No 23 ManagementExperienceCEO Yes 15 6 No nonnegative No 24 ManagementExperienceTechnical Yes 16 7 No nonnegative No 25 ManagementExperienceFinancial Yes 17 6 No nonnegative No 26 ManagementExperienceScientific Yes 18 7 No nonnegative No 27 ManagementExperienceMarketing Yes 19 6 No binary No 28 ManagementExperienceSales Yes 20 6 No fixed No 29 ManagementExperienceStartup Yes 21 6 No ordered-categorical No 30 ManagementJointExperience Yes 22 6 No ordered-categorical No 31 ManagementDegreesNSBachelors Yes 23 2 No binary No 32 ManagementDegreesNSMasters Yes 24 2 No binary No 33 ManagementDegreesNSPhD Yes 25 2 No ordered-categorical No 34 ManagementDegreesNSProf Yes 26 2 No ordered-categorical No 35 ManagementDegreesEBachelors Yes 27 2 No fixed No 36 ManagementDegreesEMasters Yes 28 2 No binary No 37 ManagementDegreesEPhD Yes 29 2 No binary No 38 ManagementDegreesEProf Yes 30 2 No fixed No 39 ManagementAverageCompensation Yes 31 5 No nonnegative No 40 ManagementStockCompensation Yes NA 0 No unordered-categorical No 41 BoDMembers Yes NA 0 No nonnegative No 42 BoDCompsitionFounder Yes 32 1 No ordered-categorical No 43 BoDCompsitionInvestor Yes 33 1 No nonnegative No 44 BoDCompsitionExecutive Yes 34 1 No ordered-categorical No 45 BoDCompsitionOther Yes NA 0 No ordered-categorical No 46 ScientificBoardMembers Yes 35 6 No nonnegative No 47 FundsTotal Yes NA 0 No positive-continuous No 48 FundsResearch Yes NA 0 No positive-continuous No 49 FundsManufacturing Yes NA 0 No ordered-categorical No 50 FundsMarketing Yes NA 0 No fixed No 51 FundsOther Yes NA 0 No binary No 52 FundsCommittedCorporate Yes NA 0 No nonnegative No 53 FundsCommittedCorporateAmount Yes 36 3 No nonnegative No 54 FundsCommittedPrivate Yes NA 0 No nonnegative No 55 FundsCommittedPrivateAmount Yes 37 6 No nonnegative No 56 FundsCommittedFounder Yes 38 1 No ordered-categorical No 57 FundsCommittedFounderAmount Yes 39 2 No ordered-categorical No 58 FundsCommittedNone Yes NA 0 No binary No 59 GovernmentalAward Yes 40 2 No binary No 60 GovernmentalAwardAmount Yes 41 9 No nonnegative No 61 NonGovernmentalAward Yes 42 3 No binary No 62 NonGovernmentalAwardAmount Yes 43 3 No ordered-categorical No 63 Label Yes 44 2 No binary No
[R] Strange behaviour of R graphics copied to PowerPoint
Hello, I am using R to create graphics, especially to plot time series charts. These charts are then copied as metafiles (for best quality) to a PowerPoint presentation and then saved to PDF (via the Save As dialog). Attached is two pictures. The first picture shows how my chart looks like in the R Graphics window, and the second picture shows how the chart becomes after saving it to PDF. http://r.789695.n4.nabble.com/file/n4668522/R.png http://r.789695.n4.nabble.com/file/n4668522/Rppt.png As you can see. After saving the metafile to PDF via PowerPoint, some straight lines appears (it seems like all of the lines has the same origin in the upper left corner and ends somewhere on the times series line). This happens in both plot() and ggplot(). The problem appears more often when using daily data in my time series. With monthly data the problem don't exist. Have anyone experienced this before? Do you think the problem is related to R or to Powerpoint? Thanks all, E [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Conditional Forecast with MSBVAR
Hi there! I'm trying to use the R package MSBVAR to generate a conditional forecast from a bayesian VAR. The condition that I am imposing for the forecast is observed values for all but one of the variables. I then want to forecast the path of that one variable, given observed values for all others. I have no problem importing the data and estimating the VAR, however, when I try and use the command hc.forecast I get the following errors: Error in solve.default((hstar.update), (crossprod(X.update, Y.update) + : system is computationally singular: reciprocal condition number = 3.6882e-17 Error in chol.default(vcv.Bh) : the leading minor of order 88 is not positive definite This is how I am using the command: conditional-hc.forecast(flatvar,condition,20,burnin=3000,gibbs=5000) where flat var is a bayesian VAR with a flat prior, condition is a matrix that specifies the condition mentioned above (observed values, and 0's for the variable I want to forecast), 20 is the number of steps I want to forecast. Any help is very much appreciated. I'm new to this so if I am leaving out information please let me know! [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] X-axis in meandist {vegan}
Dear Kumar Mainali, On 03/06/2013, at 00:01 AM, Kumar Mainali wrote: I am using vegan package for estimating dissimilarity within and between clusters, and plotting them. The following code gives me a dendrogram without X-axis. How would I add X axis to the plot? Below dataEnv is a matrix of my environmental variables. dataFact$Source_by_Type has the levels that delineate clusters. dataEnvDist - vegdist(dataEnv, method=bray) ## bray curtis dissimilarity meanDist - meandist(dataEnvDist, dataFact$Source_by_Type) plot(meanDist, ylab=Mean dissimilarity) What should the X-axis show? You can add X-axis with command axis(1), but it hardly shows anything useful. Cheers, Jari Oksanen -- Jari Oksanen, Dept Biology, Univ Oulu, 90014 Finland __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] 3.0.1 update and compiler package
Hello, would it be possible to make a decision whether this forwarded to the list or declined. Its now 2 days. Thanks Christoph On 31/05/13 14:33, Christoph Knapp wrote: Hi, I recently updated to R 3.0.1. I'm running linux ubuntu 12.04. I realized that I have to update all the installed packages so I run update.packages(checkBuilt=TRUE) as described here http://www.r-bloggers.com/r-3-0-0-is-released-whats-new-and-how-to-upgrade/ . The first thing it did was telling me that it will not update several packages. When it was finished I used the warnings function to have a look at what did not work. See the list below. warnings() Warning messages: 1: In install.packages(update[instlib == l, Package], l, ... : installation of package ‘boot’ had non-zero exit status 2: In install.packages(update[instlib == l, Package], l, ... : installation of package ‘cluster’ had non-zero exit status 3: In install.packages(update[instlib == l, Package], l, ... : installation of package ‘foreign’ had non-zero exit status 4: In install.packages(update[instlib == l, Package], l, ... : installation of package ‘KernSmooth’ had non-zero exit status 5: In install.packages(update[instlib == l, Package], l, ... : installation of package ‘MASS’ had non-zero exit status 6: In install.packages(update[instlib == l, Package], l, ... : installation of package ‘Matrix’ had non-zero exit status 7: In install.packages(update[instlib == l, Package], l, ... : installation of package ‘nlme’ had non-zero exit status 8: In install.packages(update[instlib == l, Package], l, ... : installation of package ‘nnet’ had non-zero exit status 9: In install.packages(update[instlib == l, Package], l, ... : installation of package ‘rpart’ had non-zero exit status 10: In install.packages(update[instlib == l, Package], l, ... : installation of package ‘spatial’ had non-zero exit status 11: In install.packages(update[instlib == l, Package], l, ... : installation of package ‘survival’ had non-zero exit status 12: In install.packages(update[instlib == l, Package], l, ... : installation of package ‘class’ had non-zero exit status 13: In install.packages(update[instlib == l, Package], l, ... : installation of package ‘epiR’ had non-zero exit status 14: In install.packages(update[instlib == l, Package], l, ... : installation of package ‘gmodels’ had non-zero exit status 15: In install.packages(update[instlib == l, Package], l, ... : installation of package ‘gplots’ had non-zero exit status 16: In install.packages(update[instlib == l, Package], l, ... : installation of package ‘mgcv’ had non-zero exit status 17: In install.packages(update[instlib == l, Package], l, ... : installation of package ‘gregmisc’ had non-zero exit status I tried to reinstall them manually but this always failed because of a package dependency to the compiler package. Now, if I try to install the compiler package it tells me. install.packages(compiler) Installing package into ‘/home/christoph/R/x86_64-pc-linux-gnu-library/3.0’ (as ‘lib’ is unspecified) Warning message: package ‘compiler’ is not available (for R version 3.0.1) The last line also came up all the time when the packages were updated Doing a bit of research does not deliver much only that the compiler package was included into R at version 2.13.0 (http://dirk.eddelbuettel.com/blog/2011/04/12/). Most of those packages which do not work any more are pretty important for some of my scripts and I would not even know what packages replace the packages above. Would anyone know how to fix this? Regards __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] 3.0.1 update and compiler package
Hi, How did you upgraded your version of R? From source or from a Linux package? Regards, Pascal On 05/31/2013 11:33 AM, Christoph Knapp wrote: Hi, I recently updated to R 3.0.1. I'm running linux ubuntu 12.04. I realized that I have to update all the installed packages so I run update.packages(checkBuilt=TRUE) as described here http://www.r-bloggers.com/r-3-0-0-is-released-whats-new-and-how-to-upgrade/ . The first thing it did was telling me that it will not update several packages. When it was finished I used the warnings function to have a look at what did not work. See the list below. warnings() Warning messages: 1: In install.packages(update[instlib == l, Package], l, ... : installation of package ‘boot’ had non-zero exit status 2: In install.packages(update[instlib == l, Package], l, ... : installation of package ‘cluster’ had non-zero exit status 3: In install.packages(update[instlib == l, Package], l, ... : installation of package ‘foreign’ had non-zero exit status 4: In install.packages(update[instlib == l, Package], l, ... : installation of package ‘KernSmooth’ had non-zero exit status 5: In install.packages(update[instlib == l, Package], l, ... : installation of package ‘MASS’ had non-zero exit status 6: In install.packages(update[instlib == l, Package], l, ... : installation of package ‘Matrix’ had non-zero exit status 7: In install.packages(update[instlib == l, Package], l, ... : installation of package ‘nlme’ had non-zero exit status 8: In install.packages(update[instlib == l, Package], l, ... : installation of package ‘nnet’ had non-zero exit status 9: In install.packages(update[instlib == l, Package], l, ... : installation of package ‘rpart’ had non-zero exit status 10: In install.packages(update[instlib == l, Package], l, ... : installation of package ‘spatial’ had non-zero exit status 11: In install.packages(update[instlib == l, Package], l, ... : installation of package ‘survival’ had non-zero exit status 12: In install.packages(update[instlib == l, Package], l, ... : installation of package ‘class’ had non-zero exit status 13: In install.packages(update[instlib == l, Package], l, ... : installation of package ‘epiR’ had non-zero exit status 14: In install.packages(update[instlib == l, Package], l, ... : installation of package ‘gmodels’ had non-zero exit status 15: In install.packages(update[instlib == l, Package], l, ... : installation of package ‘gplots’ had non-zero exit status 16: In install.packages(update[instlib == l, Package], l, ... : installation of package ‘mgcv’ had non-zero exit status 17: In install.packages(update[instlib == l, Package], l, ... : installation of package ‘gregmisc’ had non-zero exit status I tried to reinstall them manually but this always failed because of a package dependency to the compiler package. Now, if I try to install the compiler package it tells me. install.packages(compiler) Installing package into ‘/home/christoph/R/x86_64-pc-linux-gnu-library/3.0’ (as ‘lib’ is unspecified) Warning message: package ‘compiler’ is not available (for R version 3.0.1) The last line also came up all the time when the packages were updated Doing a bit of research does not deliver much only that the compiler package was included into R at version 2.13.0 (http://dirk.eddelbuettel.com/blog/2011/04/12/). Most of those packages which do not work any more are pretty important for some of my scripts and I would not even know what packages replace the packages above. Would anyone know how to fix this? Regards __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
