Re: [R] 6 regions around a point

2014-12-21 Thread Chel Hee Lee

You may try this:

> get.coords.circle <- function(x0, y0, r, len, rot){
+ by <- seq(from=-pi, to=pi, length.out=len)
+ x <- r*cos(by+rot) + x0
+ y <- r*sin(by+rot) + y0
+ coords <- cbind(x,y)
+ coords[!duplicated(coords),]
+ }
> xtms <- get.coords.circle(x0=2, y0=4, r=2, len=7, rot=0)
> xtms1 <- get.coords.circle(x0=2, y0=4, r=1, len=7, rot=10)
>
> plot(xtms, type="n")
> segments(x0=xtms[1:3,1], y0=xtms[1:3,2],
+ x=xtms[4:6,1], y=xtms[4:6,2])
> text(xtms[c(4,5,6),], labels=paste("line1", 1:3, sep="_"))
> text(xtms1, labels=paste("S", 1:6, sep="_"))
> points(x=2, y=4, cex=3)
>

Please see the output.  Is this what you are looking for?  I hope this 
helps.


Chel Hee Lee

On 12/21/2014 01:25 PM, eliza botto wrote:

Dear UseRs,
A point was plotted by the following command
  plot(2,4,ylim=c(0,10),xlim=c(0,5))
how to divide the space around the plotted point into six regions each of 60 
degree as shown in the Figure 2a) in the following link 
http://infolab.usc.edu/csci599/Fall2007/papers/b-2.pdf.
Thankyou very much in advance,
Eliza.




[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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 -- To UNSUBSCRIBE and more, see
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] 6 regions around a point

2014-12-21 Thread Boris Steipe
# A general approach to "lines" on a plot is provided by segments().
# However in this special case you can use abline(). 
# You have to take care though that your aspect ratio for the
# plot is exactly 1. Therefore you have to set the asp parameter.

p <- c(2,4)
plot(p[1], p[2], xlim=c(0,5), ylim=c(0,10), xlab="", ylab="", asp=1.0)
abline(h=p[2], lty=2)
abline(p[2] - (p[1]*tan(pi/3)),  tan(pi/3), lty=2)
abline(p[2] + (p[1]*tan(pi/3)), -tan(pi/3), lty=2)

# If you need to reproduce Fig. 2a more exactly:
#  - plot your frame a bit larger, don't draw axes
#  - draw the ablines
#  - draw two arrows to symbolize the coordinate axes

p <- c(2,4)
plot(p[1], p[2], xlim=c(-0.5,10.5), ylim=c(-0.5,10.5), , xlab="", ylab="", 
axes=n, asp=1.0)
abline(h=p[2], lty=2)
abline(p[2] - (p[1]*tan(pi/3)),  tan(pi/3), lty=2)
abline(p[2] + (p[1]*tan(pi/3)), -tan(pi/3), lty=2)
arrows(0, 0, 10, 0, length=0.1)
arrows(0, 0, 0, 10, length=0.1)


Cheers,
B.

On Dec 21, 2014, at 2:25 PM, eliza botto  wrote:

> Dear UseRs,
> A point was plotted by the following command
> plot(2,4,ylim=c(0,10),xlim=c(0,5))
> how to divide the space around the plotted point into six regions each of 60 
> degree as shown in the Figure 2a) in the following link 
> http://infolab.usc.edu/csci599/Fall2007/papers/b-2.pdf.
> Thankyou very much in advance,
> Eliza.
> 
> 
> 
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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 -- To UNSUBSCRIBE and more, see
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] 6 regions around a point

2014-12-21 Thread eliza botto
Dear UseRs,
A point was plotted by the following command
 plot(2,4,ylim=c(0,10),xlim=c(0,5))
how to divide the space around the plotted point into six regions each of 60 
degree as shown in the Figure 2a) in the following link 
http://infolab.usc.edu/csci599/Fall2007/papers/b-2.pdf.
Thankyou very much in advance,
Eliza.



  
[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] Bca confidence intervals for Pearson coefficient, thanks.

2014-12-21 Thread varin sacha
Dear Professor FOX,
I really thank you lots for all your precisions. One last precision, now if I 
want tot calculate the BCa bootstrap CIs for the Cramer's V and the Eta-squared.
## Read in the data file, using headers and Tab separator> test = 
read.table(file.choose(), header = TRUE, sep = "\t")
## Check the variable names> names(test)[1] "gender"    "option"    "math.test" 
"geo.test"  "shopping"  "sports"   
## Cramer's V
> library(questionr)> tab = table(test$gender, test$option)> cramer.v(tab)      
>                    [1] 0.1490712                           
## Eta.square
> library(lsr)> test.aov = aov(math.test ~ gender, data = test)
> etaSquared(test.aov)          eta.sq eta.sq.partgender 0.1207154   0.1207154

Best Regards, looking forward to reading you once more,
Sacha
  De : John Fox 
 À : varin sacha  
Cc : r-help help  
 Envoyé le : Dimanche 21 décembre 2014 5h33
 Objet : Re: [R] Bca confidence intervals for Pearson coefficient, thanks.
   
Dear varin sacha,

I think that you misunderstand how boot() and boot.ci() work. The boot() 
function in the simplest case takes two arguments, for the data and indices 
into the data, while boot.ci() takes as its principal argument the object 
returned by boot(). All of this seems reasonably clear in ?boot and ?boot.ci.

Here's an example with different data (since as far as I can see you didn't 
supply yours):

- snip 

> library(boot)
> 
> x <- longley$Year
> y <- longley$Population
> 
> cor(cbind(x, y))
          x        y
x 1.000 0.9939528
y 0.9939528 1.000
> 
> myCor <- function(data, index){
+  cor(data[index, ])[1, 2]
+ }
> 
> set.seed(12345)
> (b <- boot(data=cbind(x, y), statistic=myCor, R=200))

ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = cbind(x, y), statistic = myCor, R = 200)


Bootstrap Statistics :
    original      bias    std. error
t1* 0.9939528 0.0008263766 0.001850004

> boot.ci(b, type="bca")
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 200 bootstrap replicates

CALL : 
boot.ci(boot.out = b, type = "bca")

Intervals : 
Level      BCa          
95%  ( 0.9895,  0.9969 )  
Calculations and Intervals on Original Scale
Warning : BCa Intervals used Extreme Quantiles
Some BCa intervals may be unstable
Warning message:
In norm.inter(t, adj.alpha) : extreme order statistics used as endpoints

--- snip ---

Note that 200 bootstrap replications are generally sufficient for bootstrap 
standard errors (a normal-theory CI would be a poor choice here, unless you 
transform the correlation coefficient), but really aren't enough for a BCa 
interval.

I hope this helps,
 John




John Fox, Professor
McMaster University
Hamilton, Ontario, Canada
http://socserv.mcmaster.ca/jfox/
    
On Sat, 20 Dec 2014 20:36:57 + (UTC)

> Hi to everyone,
> I am trying to calculate the Bca bootstrap confidence intervals for the 
> Pearson coefficient.
> x=Dataset$math.testy=Dataset$geo.testcor(x,y,method="pearson")[1] 0.6983799
> boot.ci(cor, conf=0.95, type=bca)Erreur dans boot.out$t0 : objet de type 
> 'closure' non indiçable
> 
> I have tried as well to calculate the Pearson coefficient using bootstrap and 
> then to calculate the Bca bootstrap CIs of the Pearson. It doesn't work 
> either. 
> boot(data = cbind(x, y), statistic = cor, R = 200)
> 
> ORDINARY NONPARAMETRIC BOOTSTRAP
> 
> 
> Call:
> boot(data = cbind(x, y), statistic = cor, R = 200)
> 
> 
> Bootstrap Statistics :
>      original    bias    std. error
> t1* -0.6243713 0.6295142  0.2506267
> t2* -0.1366533 0.1565392  0.2579134
> > boot.ci(cor, conf=0.95, type=bca)
> Erreur dans boot.out$t0 : objet de type 'closure' non indiçable
> Many thanks to tell me how to correct my R script to get the Bca CIs for my 
> Pearson coefficient. Best regards, looking forward to reading you,
> SV
> 
>     [[alternative HTML version deleted]]
> 
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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 -- To UNSUBSCRIBE and more, see
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] how update to latest version of R on mac

2014-12-21 Thread Jane Synnergren
Thanks for your response Tom, now I have successfully upgraded and my affy 
package problem seems to be is solved with this latest version of R and latest 
version of Bioconductor.

Thanks!
/Jane

---
Jane Synnergren, PhD
Systems Biology Research Center
School of Bioscience
University of Skövde
Sweden

Contact:
email: jane.synnerg...@his.se
mobile: +46 (0)708 806495
webpage: his.se/synj

Postal Address:  Visiting address:
Box 408 Kanikegränd 3A
541 28 Skövde Skövde

From: Thomas Adams mailto:tea...@gmail.com>>
Date: Sunday, December 21, 2014 4:35 PM
To: Jane Synnergren mailto:jane.synnerg...@his.se>>
Cc: R-help mailto:r-help@r-project.org>>
Subject: Re: [R] how update to latest version of R on mac

Jane,

My sincere apologies; I don't know what I was seeing/thinking -- your original 
post was correct. Namely, you should download and install: 
http://ftp.sunet.se/pub/lang/CRAN/bin/macosx/R-3.1.2-snowleopard.pkg (although, 
I am use to using the main mirror 
http://cran.r-project.org/bin/macosx/R-3.1.2-snowleopard.pkg) -- but that 
should not make a difference. The documentation does state...

R 3.1.2 binary for Mac OS X 10.6 (Snow Leopard) and higher, signed package. 
Contains R 3.1.2 framework, R.app GUI 1.65 in 64-bit for Intel Macs. The above 
file is an Installer package which can be installed by double-clicking. 
Depending on your browser, you may need to press the control key and click on 
this link to download the file.

So, this should be the way to go. I'm glad you double checked!!

Tom

On Sun, Dec 21, 2014 at 2:21 AM, Jane Synnergren 
mailto:jane.synnerg...@his.se>> wrote:
Thanks Tom!
One more quick question, regarding choosing mavericks.pkg or snowleopard.pkg.
Will mavericks work for me though I have OS X Lion 10.7.5? I thought I had to 
choose snowleopard but are still an unexperienced mac user and have really no 
idea.

Just want to dubble check before I run the installation.

/Jane
---
Jane Synnergren, PhD
Systems Biology Research Center
School of Bioscience
University of Skövde
Sweden

Contact:
email: jane.synnerg...@his.se
mobile: +46 (0)708 806495
webpage: his.se/synj

Postal Address:  Visiting address:
Box 408 Kanikegränd 3A
541 28 Skövde Skövde

From: Thomas Adams mailto:tea...@gmail.com>>
Date: Sunday, December 21, 2014 2:06 AM
To: Jane Synnergren mailto:jane.synnerg...@his.se>>
Cc: R-help mailto:r-help@r-project.org>>
Subject: Re: [R] how update to latest version of R on mac

Jane,

It's possible, since you are using Mac OS X Lion (10.7.5), you need the R 
version for that, which can be found here: 
http://cran.r-project.org/bin/macosx/old/R-3.1.1-mavericks.pkg

Cheers!
Tom

On Fri, Dec 19, 2014 at 2:56 PM, Jane Synnergren 
mailto:jane.synnerg...@his.se>> wrote:
Hi,
I get following error when trying to update the affy package to latest
version.

Error: cannot remove prior installation of package Œaffy¹

I think it is because affy requires R 3.1.1 and I have 3.0.2

I try to find a guide how I do to upgrade from 3.0.2 to latest version on
mac, but cannot find any good description. Can anyone guide me?

Do I just download R for mac (R-3.1.2-snowleopard.pkg) from
 and
dubbelclick?

What will happen with all installed packages?
I have OS X Lion 10.7.5
Which version number is stable and reliable?

I also need a guide of how to update to latest version of  bioconductor on
mac.

Descriptions for windows is available everywhere.

/Jane


---

Jane Synnergren, PhD
Systems Biology Research Center
School of Bioscience
University of Skövde
Sweden

Contact:
email: jane.synnerg...@his.se
mobile: +46 (0)708 806495
webpage: his.se/synj

Postal Address: Visiting address:

Box 408 Kanikegränd 3A
541 28 Skövde   Skövde

__
R-help@r-project.org mailing list -- To 
UNSUBSCRIBE and more, see
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 -- To UNSUBSCRIBE and more, see
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] Calculating mean, median, minimum, and maximum

2014-12-21 Thread David Winsemius

On Dec 19, 2014, at 9:35 AM, Henrik Bengtsson wrote:

> On Fri, Dec 19, 2014 at 8:48 AM, John Kane  wrote:
>> Hi,
>> It looks like you are replying to some phantom. Is your correspondent 
>> actully on R-help?
> 
> This most likely happens because OP posted the message via Nabble
> [http://r.789695.n4.nabble.com/Calculating-mean-median-minimum-and-maximum-td4700862.html]
> via "New Topic", which in turn post the message to r-help after asking
> the user to confirm:
> 
> "Mailing List Subscription Reminder
> This forum is an archive/gateway which will forward your post to the
> r-help@r-project.org mailing list.
> 
> The mailing list may require your subscription before accepting your
> post. Please note that being registered with Nabble does NOT
> automatically subscribe you to this mailing list. If you haven't
> subscribed yet, please do it now. If you aren't sure or don't
> remember, just subscribe again because there is no harm."

Some of that statement is incorrect. You will be unable to establish a new 
subscription if your email address has not changed, since it is the email 
address that acts as your  useR_name. Furthermore, if your email address 
changes, you can change the address without re-subscribing as long as you log 
in with the old email address and your password. The address at which you gain 
access to the web interface is posted in the standard footer in each message 
(unless you are using Nabble).

> 
> Since OP is not a registered user on the r-help mailing list, his/her
> post is help by r-help for moderation.  However, mtrang probably saw
> it on Nabble and replied there and since mtrang is also a registered
> user on the r-help, his/her messages reach us here on r-help, before
> OP's messages help for moderation are approved.
> 
> I'd say, this is quite annoying "feature" to everyone but Nabble
> users, particularly since the Nabble message does not include the
> previous messages (unlike email).
> 

Writing as a moderator (but not an owner);

Sometime ago (guessing Jan 2013) the automatic filters were "strengthened" in 
the ETHZ mail-server that automatically discards some Nabble postings. This was 
done because someone set up a spam bot that was posting exact duplicates of 
prior postings altered only by attaching a spam footer that had links to gawd 
knows where. These were coming in too fast to filter by the moderations staff 
so it was decided to summarily chuck them in null_dev.  The exact rules were 
never published, but I have observed that hotmail, gmail, and yahoo addresses 
will increase the probability of automatic rejection. I think that replies to 
these are sometimes held up in the moderation queue as well. I do believe that 
subscription status is also part of that algorithm.

Nabble not only makes it modestly difficult to attach prior message content but 
it also filters out the footer at the bottom of every posting so the 
subscription address and the Posting Guide address are not included. It also 
falsely claims that it is the Rhelp mailing list archive. We, of course, know 
that there is only One True Archive: https://stat.ethz.ch/pipermail/r-help/

It is possible to use Nabble (or gmail or yahoo)  as a web-front-end 
successfully and respectfully in accord with mailing list traditions of 
in-context editing and adherence to group norms. It requires an understanding 
of mailing list realities and attention to details that most Nabble users seem 
unwilling or unable to learn. They have over the years imposed considerable 
effort and time on the moderation volunteers. As a group I hope it is fair to 
say we moderators have little sympathy for too typical Nabble user who fails to 
read the Posting Guide (to which there does appear a link when users log in), 
although we do realize that a few longtime and valued contributors may use its 
facility when traveling.

-- 
David



> /Henrik
> 
>> 
>> John Kane
>> Kingston ON Canada
>> 
>> 
>>> -Original Message-
>>> From: mtr...@buffalo.edu
>>> Sent: Fri, 19 Dec 2014 08:08:29 -0800 (PST)
>>> To: r-help@r-project.org
>>> Subject: Re: [R] Calculating mean, median, minimum, and maximum
>>> 
>>> You can use the apply function which "applies" a function of your choice,
>>> and
>>> MARGIN = 2 means you want to do it columnwise:
>>> 
 apply(X = df, MARGIN=2, FUN = mean, na.rm = TRUE)
>>> Latitude Longitude   January  February March April   May
>>> June
>>>  26.9380 -109.8125  159.8454  156.4489  153.6911  150.1719  148.0885
>>> 149.2365
 apply(X = df, MARGIN=2, FUN = min, na.rm = TRUE)
>>> Latitude Longitude   January  February March April   May
>>> June
>>>   26.938  -110.688   121.204   118.713   117.293   114.398   112.357
>>> 113.910
 apply(X = df, MARGIN=2, FUN = max, na.rm = TRUE)
>>> Latitude Longitude   January  February March April   May
>>> June
>>>   26.938  -108.938   252.890   248.991   244.870   241.194   239.615
>>> 239.888
 apply(X = df, MARGIN=2, FUN 

Re: [R] MLE

2014-12-21 Thread Arne Henningsen
Dear Pari

On 21 December 2014 at 06:59, pari hesabi  wrote:
> I am trying to get the Maximum Likelihood Estimation of a parameter
> in a probability mass function. The problem is my pmf which includes
> a summation and one integral. it is not similar to other known pdfs or
> pmfs such as normal, exponential, poisson, .
> Does anybody know whether I can use the current packages(like
> Maxlik) in R for getting the MLE of the parameter?

maxLik (not Maxlik) will probably work.

> Can anybody explains me with an example?

R> library( "maxLik" )
R> ?maxLik

http://dx.doi.org/10.1007/s00180-010-0217-1

https://absalon.itslearning.com/data/ku/103018/publications/maxLik.pdf

> I would appreciate any help.

http://www.R-project.org/posting-guide.html

http://maxlik.org/

https://r-forge.r-project.org/projects/maxlik/

Best regards,
Arne

-- 
Arne Henningsen
http://www.arne-henningsen.name

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] Bca confidence intervals for Pearson coefficient, thanks.

2014-12-21 Thread John Fox
Dear Sacha,

Simply write a function that takes a data set and index vector as arguments, 
say statistic(data, index), and have it compute and return either eta^2 or V 
depending upon the application. Use the function myCor() in the previous 
example as a model.

Best,
 John 

> -Original Message-
> From: varin sacha [mailto:varinsa...@yahoo.fr]
> Sent: Sunday, December 21, 2014 3:03 PM
> To: John Fox
> Cc: r-help help
> Subject: Re: [R] Bca confidence intervals for Pearson coefficient,
> thanks.
> 
> Dear Professor FOX,
> 
> 
> I really thank you lots for all your precisions.
> One last precision, now if I want tot calculate the BCa bootstrap CIs
> for the Cramer's V and the Eta-squared.
> 
> 
> ## Read in the data file, using headers and Tab separator
> > test = read.table(file.choose(), header = TRUE, sep = "\t")
> 
> 
> ## Check the variable names
> > names(test)
> [1] "gender""option""math.test" "geo.test"  "shopping"  "sports"
> 
> 
> ## Cramer's V
> 
> > library(questionr)
> > tab = table(test$gender, test$option)
> > cramer.v(tab)
> [1] 0.1490712
> 
> ## Eta.square
> 
> > library(lsr)
> > test.aov = aov(math.test ~ gender, data = test)
> 
> > etaSquared(test.aov)
>   eta.sq eta.sq.part
> gender 0.1207154   0.1207154
> 
> 
> 
> Best Regards, looking forward to reading you once more,
> 
> Sacha
> 
> 
> 
> De : John Fox 
> À : varin sacha 
> Cc : r-help help 
> Envoyé le : Dimanche 21 décembre 2014 5h33
> Objet : Re: [R] Bca confidence intervals for Pearson coefficient,
> thanks.
> 
> 
> Dear varin sacha,
> 
> I think that you misunderstand how boot() and boot.ci() work. The boot()
> function in the simplest case takes two arguments, for the data and
> indices into the data, while boot.ci() takes as its principal argument
> the object returned by boot(). All of this seems reasonably clear in
> ?boot and ?boot.ci.
> 
> Here's an example with different data (since as far as I can see you
> didn't supply yours):
> 
> - snip 
> 
> > library(boot)
> >
> > x <- longley$Year
> > y <- longley$Population
> >
> > cor(cbind(x, y))
>   xy
> x 1.000 0.9939528
> y 0.9939528 1.000
> >
> > myCor <- function(data, index){
> +  cor(data[index, ])[1, 2]
> + }
> >
> > set.seed(12345)
> > (b <- boot(data=cbind(x, y), statistic=myCor, R=200))
> 
> ORDINARY NONPARAMETRIC BOOTSTRAP
> 
> 
> Call:
> boot(data = cbind(x, y), statistic = myCor, R = 200)
> 
> 
> Bootstrap Statistics :
> original  biasstd. error
> t1* 0.9939528 0.0008263766 0.001850004
> 
> > boot.ci(b, type="bca")
> BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
> Based on 200 bootstrap replicates
> 
> CALL :
> boot.ci(boot.out = b, type = "bca")
> 
> Intervals :
> Level  BCa
> 95%  ( 0.9895,  0.9969 )
> Calculations and Intervals on Original Scale
> Warning : BCa Intervals used Extreme Quantiles
> Some BCa intervals may be unstable
> Warning message:
> In norm.inter(t, adj.alpha) : extreme order statistics used as endpoints
> 
> --- snip ---
> 
> Note that 200 bootstrap replications are generally sufficient for
> bootstrap standard errors (a normal-theory CI would be a poor choice
> here, unless you transform the correlation coefficient), but really
> aren't enough for a BCa interval.
> 
> I hope this helps,
> John
> 
> 
> 
> 
> 
> John Fox, Professor
> McMaster University
> Hamilton, Ontario, Canada
> http://socserv.mcmaster.ca/jfox/
> 
> On Sat, 20 Dec 2014 20:36:57 + (UTC)
> varin sacha  wrote:
> > Hi to everyone,
> > I am trying to calculate the Bca bootstrap confidence intervals for
> the Pearson coefficient.
> > x=Dataset$math.testy=Dataset$geo.testcor(x,y,method="pearson")[1]
> 0.6983799
> > boot.ci(cor, conf=0.95, type=bca)Erreur dans boot.out$t0 : objet de
> type 'closure' non indiçable
> >
> > I have tried as well to calculate the Pearson coefficient using
> bootstrap and then to calculate the Bca bootstrap CIs of the Pearson. It
> doesn't work either.
> > boot(data = cbind(x, y), statistic = cor, R = 200)
> >
> > ORDINARY NONPARAMETRIC BOOTSTRAP
> >
> >
> > Call:
> > boot(data = cbind(x, y), statistic = cor, R = 200)
> >
> >
> > Bootstrap Statistics :
> >  originalbiasstd. error
> > t1* -0.6243713 0.6295142  0.2506267
> > t2* -0.1366533 0.1565392  0.2579134
> > > boot.ci(cor, conf=0.95, type=bca)
> > Erreur dans boot.out$t0 : objet de type 'closure' non indiçable
> > Many thanks to tell me how to correct my R script to get the Bca CIs
> for my Pearson coefficient. Best regards, looking forward to reading
> you,
> > SV
> 
> >
> > [[alternative HTML version deleted]]
> >
> > __
> > R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> > 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, minima

Re: [R] How to create a time series object with time only (no date)

2014-12-21 Thread Gabor Grothendieck
On Sun, Dec 21, 2014 at 12:09 AM, ce  wrote:
>
> Dear all,
>
> I want to create a time series object from 00:00:00 to 23:59:00 without dates 
> ?
> I can't figure it out with xts ?
>

This uses zoo, rather than xts:

library(zoo)
library(chron)

tt <- seq(times("00:00:00"), times("23:59:00"), by = times("00:01:00"))
dat <- 1:1440
z <- zoo(dat, tt)


-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] list of lists, is this element empty

2014-12-21 Thread William Dunlap
Your 'x' has length 2, so x[[3]] cannot be calculated ('subscript out of
bounds' is what I get).  You can check for this with length(x)<3.

In general, you want to be more precise: 'does not have a value', 'is
NULL', and 'is empty' are not synonymous.  I'm not sure what 'does not have
a value' means to you.  NULL is a value with a certain type, 'NULL', and
length, 0.  seq_len(0) is empty, but not NULL, since it has type 'integer'.
 c(1,NA) contains a 'missing value'.  Can you explain what what you are
trying to check for, giving some context?



Bill Dunlap
TIBCO Software
wdunlap tibco.com

On Sat, Dec 20, 2014 at 7:58 AM, Ragia Ibrahim  wrote:
>
> Hello,
> Kindly I have a list of lists as follow
> x
> [[1]]
> [1] 7
>
> [[2]]
> [1] 3 4 5
>
> as showen x[[3]] does not have a value and it has NULL, how can I check on
> this
> how to test if x[[3]] is empty.
>
> thanks in advance
> Ragia
>
> [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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 -- To UNSUBSCRIBE and more, see
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] Subject: short-sale constraint with nonpositive-definite matrix in portfolio optimization

2014-12-21 Thread Ravi Varadhan
Hi,

You can try a projected gradient approach, which is implemented in the spg() 
function in the "BB" package. You have to provide a projection function which 
will take an infeasible matrix as input and will give a feasible (i.e. 
positive-definite) matrix as output. This kind of matrix projection is quite 
easy to do, and in fact, there are functions in R to do this (e.g., see 
posdefify() function in "sfsmisc" package).


There is a recent review article on spectral projected gradient algorithm in J 
Statistical Software.

http://www.jstatsoft.org/v60/i03


Hope this is helpful,

Ravi

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] where are the NEWS for R <3.0.0?

2014-12-21 Thread Ben Marwick
Thanks, yes I found them in my local installation, just as you said. In 
case anyone else wants to get them programatically, here's what I did 
(on windows 7, R version 3.1.1):


# navigate to local R installation dir
setwd(system.file())
setwd('../../doc')

# get all news files
news_files <- list.files( pattern = "NEWS", full.names = TRUE, recursive 
= TRUE)


Thanks again,

Ben

On 20/12/2014 6:22 PM, Duncan Murdoch wrote:

On 19/12/2014, 5:48 PM, Ben Marwick wrote:

I'm looking for the NEWS files for versions of R before 3.0.0, does
anyone know where they are? I'm interested in getting the names of
people who have been acknowledged when changes are made to R.

At the bottom of http://cran.r-project.org/src/base/NEWS.html it says

"Older news can be found in text format in files NEWS.0, NEWS.1 and
NEWS.2 in the ‘doc’ directory. News in HTML format for R versions from
2.10.0 to 2.15.3 is in NEWS.2.html."

The text includes these URLs:

http://cran.r-project.org/src/NEWS.0
http://cran.r-project.org/src/NEWS.1
http://cran.r-project.org/src/NEWS.2
http://cran.r-project.org/src/base/NEWS.2.html

All of these go to 404s, and when src is replaced by doc, as suggested
by the text, they also go to 404s. There are no NEWS files in here
either: http://cran.r-project.org/doc/html/ So there are a bunch of
broken links here.


If you read these from within the HTML help system in R, the links
aren't broken.  CRAN copied some of the files from there, but not all.

You can also find the files in the RHOME/doc directory (and subdirectories).

Duncan Murdoch



I found a bunch of PDFs listed here http://cran.r-project.org/doc/Rnews/
dating from 2008-2001. Are these only source of the earlier NEWS files?

They're less convenient because the NEWS section is not a single plain
text file, but a section of the PDF document. I'd be most grateful to
know where I can find html or txt files of the earlier NEWS files. If
the PDFs are all there are, it would be good to know that for sure.

The other method I was considering was scraping bugzilla for names of
people who submitted bugs that resulted in a fix
(https://bugs.r-project.org/bugzilla3/buglist.cgi?product=R&query_format=advanced&resolution=FIXED).
But text mining the NEWS files would be quicker for a first approximation.

thanks,

Ben





__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] mediation analysis in R using the method described by Lange T., et al. in the paper "A simple unified approach for estimating natural direct and direct effects

2014-12-21 Thread Ruzan Udumyan
Dear All,

I need to do  a mediation analysis with survival data (using Cox model).

   - The independent variable is categorical (with 3 levels, coded as 1, 2,
   3).
   - Mediator variable is a 1-10 score variable (derived from a continuous
   variable).

I must confess, I am a Stata user. However, could not find an appropriate
command in Stata for mediation analysis with Cox model, and was very happy
to find a way of doing it in R.
I am running into problems though, which I do not know how to solve. More
specifically, I get an error message from RStudio after running the syntax
mentioned below.

I am kindly asking you to help me fix the problem.

With many thanks and my best regards,

Ruzan


Error in model.frame.default(formula = Surv(time, event) ~ factor(stress) +  :
  variable lengths differ (found for 'factor(stress)')
 In addition: Warning messages:1: In checkwz(wz, M = M, trace = trace,
wzepsilon = control$wzepsilon) :
  31 elements replaced by 1.819e-122: In checkwz(wz, M = M, trace =
trace, wzepsilon = control$wzepsilon) :
  187 elements replaced by 1.819e-123: In checkwz(wz, M = M, trace =
trace, wzepsilon = control$wzepsilon) :
  1492 elements replaced by 1.819e-124: In checkwz(wz, M = M, trace =
trace, wzepsilon = control$wzepsilon) :
  2081 elements replaced by 1.819e-125: In checkwz(wz, M = M, trace =
trace, wzepsilon = control$wzepsilon) :
  2081 elements replaced by 1.819e-126: In checkwz(wz, M = M, trace =
trace, wzepsilon = control$wzepsilon) :
  4783 elements replaced by 1.819e-127: In checkwz(wz, M = M, trace =
trace, wzepsilon = control$wzepsilon) :
  8309 elements replaced by 1.819e-128: In checkwz(wz, M = M, trace =
trace, wzepsilon = control$wzepsilon) :
  8642 elements replaced by 1.819e-129: In checkwz(wz, M = M, trace =
trace, wzepsilon = control$wzepsilon) :
  8642 elements replaced by 1.819e-1210: In checkwz(wz, M = M, trace =
trace, wzepsilon = control$wzepsilon) :
  8642 elements replaced by 1.819e-12

>
Here is the syntax:

doEffectDecomp = function(myData)
 {
   #step 1
   myData$stresstemp <- myData$stress
   library(VGAM)
   fitM <-vglm(phys1 ~ factor(stresstemp)+ C, data = myData, family =
multinomial())  # C stands for confounders
   #summary(fitM)
 #Step 2. Next we construct new variables id and star
 N <- nrow(myData)
 myData$id <- 1:N # construct id variable
 levelsOfstress <- unique(myData$stress)
 myData1 <- myData
 myData2 <- myData
 myData3 <- myData
 myData1$star <- levelsOfstress[1]
myData2$star <- levelsOfstress[2]
   myData3$star <- levelsOfstress[3]
   newMyData <- rbind(myData1, myData2, myData3)

#Step 3. compute weights
 newMyData$stresstemp <- newMyData$stress
tempDir <- as.matrix(predict(fitM,type = "response",
newdata=newMyData))[cbind(1:(3*N),newMyData$phys1)]
   newMyData$stresstemp <- newMyData$star
   tempIndir <- as.matrix(predict(fitM,type = "response",
newdata=newMyData))[cbind(1:(3*N),newMyData$phys1)]
   newMyData$weightM <- tempIndir/tempDir
   #hist(newMyData$weightM)
   #Step 4. Finally the MSM model for direct and indirect effects can be
estimated.
# weighted Cox model
library(survival)
   cox <- coxph(Surv(time,event) ~factor(stress) + factor(star) +  C,
method="efron", data=newMyData, weights=newMyData$weightM)

  summary(cox)

  # return value: TE; DE; IE
   # Return value: Estimates for total, direct, indirect effect
   TE = exp(sum(coef(cox)[c('stressTRUE', 'starTRUE')])) # I am not
sure this is correct for the categorical exposure with 3 levels
   DE = exp(unname(coef(cox)['stressTRUE']))
   IE = exp(sum(coef(cox)['starTRUE']))
   PM = log(IE) / log(TE)
return(c(exp(coef(cox)), TE=TE, DE=DE, IE=IE, PM=PM))
 }



doEffectDecomp(myData)

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] how update to latest version of R on mac

2014-12-21 Thread Jane Synnergren
Thanks Tom!
One more quick question, regarding choosing mavericks.pkg or snowleopard.pkg.
Will mavericks work for me though I have OS X Lion 10.7.5? I thought I had to 
choose snowleopard but are still an unexperienced mac user and have really no 
idea.

Just want to dubble check before I run the installation.

/Jane
---
Jane Synnergren, PhD
Systems Biology Research Center
School of Bioscience
University of Skövde
Sweden

Contact:
email: jane.synnerg...@his.se
mobile: +46 (0)708 806495
webpage: his.se/synj

Postal Address:  Visiting address:
Box 408 Kanikegränd 3A
541 28 Skövde Skövde

From: Thomas Adams mailto:tea...@gmail.com>>
Date: Sunday, December 21, 2014 2:06 AM
To: Jane Synnergren mailto:jane.synnerg...@his.se>>
Cc: R-help mailto:r-help@r-project.org>>
Subject: Re: [R] how update to latest version of R on mac

Jane,

It's possible, since you are using Mac OS X Lion (10.7.5), you need the R 
version for that, which can be found here: 
http://cran.r-project.org/bin/macosx/old/R-3.1.1-mavericks.pkg

Cheers!
Tom

On Fri, Dec 19, 2014 at 2:56 PM, Jane Synnergren 
mailto:jane.synnerg...@his.se>> wrote:
Hi,
I get following error when trying to update the affy package to latest
version.

Error: cannot remove prior installation of package Œaffy¹

I think it is because affy requires R 3.1.1 and I have 3.0.2

I try to find a guide how I do to upgrade from 3.0.2 to latest version on
mac, but cannot find any good description. Can anyone guide me?

Do I just download R for mac (R-3.1.2-snowleopard.pkg) from
 and
dubbelclick?

What will happen with all installed packages?
I have OS X Lion 10.7.5
Which version number is stable and reliable?

I also need a guide of how to update to latest version of  bioconductor on
mac.

Descriptions for windows is available everywhere.

/Jane


---

Jane Synnergren, PhD
Systems Biology Research Center
School of Bioscience
University of Skövde
Sweden

Contact:
email: jane.synnerg...@his.se
mobile: +46 (0)708 806495
webpage: his.se/synj

Postal Address: Visiting address:

Box 408 Kanikegränd 3A
541 28 Skövde   Skövde

__
R-help@r-project.org mailing list -- To 
UNSUBSCRIBE and more, see
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 -- To UNSUBSCRIBE and more, see
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] MLE

2014-12-21 Thread pari hesabi
Dear Sir/Madam
I am trying to get the Maximum Likelihood Estimation of a parameter in a 
probability mass function. The problem is my pmf which includes a summation and 
one integral. it is not similar to other known pdfs or pmfs such as normal, 
exponential, poisson, .
Does anybody know whether I can use the current packages(like Maxlik) in R for 
getting the MLE of the parameter?
Can anybody explains me with an example?
I would appreciate any help.
Best Regards,
Pari  
__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] how update to latest version of R on mac

2014-12-21 Thread Thomas Adams
Jane,

My sincere apologies; I don't know what I was seeing/thinking -- your
original post was correct. Namely, you should download and install:
http://ftp.sunet.se/pub/lang/CRAN/bin/macosx/R-3.1.2-snowleopard.pkg
(although, I am use to using the main mirror
http://cran.r-project.org/bin/macosx/R-3.1.2-snowleopard.pkg) -- but that
should not make a difference. The documentation does state...

*R 3.1.2* binary for Mac OS X 10.6 (Snow Leopard) *and higher*, signed
package. Contains R 3.1.2 framework, R.app GUI 1.65 in 64-bit for Intel
Macs. The above file is an Installer package which can be installed by
double-clicking. Depending on your browser, you may need to press the
control key and click on this link to download the file.

So, this should be the way to go. I'm glad you double checked!!

Tom

On Sun, Dec 21, 2014 at 2:21 AM, Jane Synnergren 
wrote:

>   Thanks Tom!
> One more quick question, regarding choosing mavericks.pkg or
> snowleopard.pkg.
> Will mavericks work for me though I have OS X Lion 10.7.5? I thought I had
> to choose snowleopard but are still an unexperienced mac user and have
> really no idea.
>
>  Just want to dubble check before I run the installation.
>
>  /Jane
>
> ---
> Jane Synnergren, PhD
> Systems Biology Research Center
> School of Bioscience
> University of Skövde
> Sweden
>
>  Contact:
> email: jane.synnerg...@his.se
> mobile: +46 (0)708 806495
> webpage: his.se/synj
>
>  Postal Address:  Visiting address:
>  Box 408 Kanikegränd 3A
> 541 28 Skövde Skövde
>
>   From: Thomas Adams 
> Date: Sunday, December 21, 2014 2:06 AM
> To: Jane Synnergren 
> Cc: R-help 
> Subject: Re: [R] how update to latest version of R on mac
>
>Jane,
>
>  It's possible, since you are using Mac OS X Lion (10.7.5), you need the R
> version for that, which can be found here:
> http://cran.r-project.org/bin/macosx/old/R-3.1.1-mavericks.pkg
>
>  Cheers!
>  Tom
>
> On Fri, Dec 19, 2014 at 2:56 PM, Jane Synnergren 
> wrote:
>>
>> Hi,
>> I get following error when trying to update the affy package to latest
>> version.
>>
>> Error: cannot remove prior installation of package Œaffy¹
>>
>> I think it is because affy requires R 3.1.1 and I have 3.0.2
>>
>> I try to find a guide how I do to upgrade from 3.0.2 to latest version on
>> mac, but cannot find any good description. Can anyone guide me?
>>
>> Do I just download R for mac (R-3.1.2-snowleopard.pkg) from
>> 
>> and
>> dubbelclick?
>>
>> What will happen with all installed packages?
>> I have OS X Lion 10.7.5
>> Which version number is stable and reliable?
>>
>> I also need a guide of how to update to latest version of  bioconductor on
>> mac.
>>
>> Descriptions for windows is available everywhere.
>>
>> /Jane
>>
>>
>>
>> ---
>> 
>> Jane Synnergren, PhD
>> Systems Biology Research Center
>> School of Bioscience
>> University of Skövde
>> Sweden
>>
>> Contact:
>> email: jane.synnerg...@his.se
>> mobile: +46 (0)708 806495
>> webpage: his.se/synj
>>
>> Postal Address: Visiting address:
>>
>> Box 408 Kanikegränd 3A
>> 541 28 Skövde   Skövde
>>
>> __
>> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
>> 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 -- To UNSUBSCRIBE and more, see
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] R - Aggregate 3-Hourly Block Data into Weekly (Melt)

2014-12-21 Thread Shouro Dasgupta
Apologies for re-posting. Any suggestions for clipping/subsetting panel
data (ID=US counties, Time=Week)? I would like start the data to start on
the first Monday and end on the last Sunday for the time period 2026-2045.
Thanks and apologies again.

Sincerely,

Shouro

On Fri, Dec 19, 2014 at 11:01 AM, Shouro Dasgupta  wrote:

> Thank you very much for your reply. I really appreciate it. I apologize
> for the HTML version, I have made modifications and replied to your
> questions/comments below. Thanks again
>
>
>
> tmp1 <- structure(list(FIPS = c(1001L, 1003L, 1005L), X2026.01.01.1 = c(
> 285.5533142,
>
>   285.5533142, 286.2481079), X2026.01.01.2 = c(283.4977112, 283.4977112,
>
>   285.0860291), X2026.01.01.3 = c(281.9733887, 281.9733887, 284.1548767
>
>   ), X2026.01.01.4 = c(280.0234985, 280.0234985, 282.6075745),
>
>   X2026.01.01.5 = c(278.7125854, 278.7125854, 281.2553711),
>
>   X2026.01.01.6 = c(278.5204773, 278.5204773, 280.6148071),
> X2026.01.01.7 = c(282.3938, 282.3938, 283.1096), X2026.01.01.8 = c(
> 285.9133, 285.9133, 286.1951) .Names = c("FIPS", "X2026.01.01.1",
> "X2026.01.01.2", "X2026.01.01.3", "X2026.01.01.4",
>
>   "X2026.01.01.5", "X2026.01.01.6", "X2026.01.01.7", "X2026.01.01.8" ),
> class = "data.frame", row.names = c(NA, -3L))
>
>
>
> *Looks like this*
>
> FIPS
>
> X2026.01.01.1
>
> X2026.01.01.2
>
> X2026.01.01.3
>
> X2026.01.01.4
>
> X2026.01.01.5
>
> X2026.01.01.6
>
> X2026.01.01.7
>
> X2026.01.01.8
>
> 1001
>
> 285.5533
>
> 283.4977
>
> 281.9734
>
> 280.0235
>
> 278.7126
>
> 278.5205
>
> 282.3938
>
> 285.9133
>
> 1003
>
> 285.5533
>
> 283.4977
>
> 281.9734
>
> 280.0235
>
> 278.7126
>
> 278.5205
>
> 282.3938
>
> 285.9133
>
> 1005
>
> 286.2481
>
> 285.086
>
> 284.1549
>
> 282.6076
>
> 281.2554
>
> 280.6148
>
> 283.1096
>
> 286.1951
>
>
>
> For X2026.01.01.1 represents Year=2026, Month=01, Day=01, Hour block=1.
>
> I have extracted the data by FIPS code and reshaped the yearly data files
> using melt();
>
>
>
> for (i in filelist) {
>
>   tmp1 <- as.data.table(read.csv(i,header=T, sep=","))
>
>   tmp2 <- melt(tmp1, id="FIPS")
>
>   tmp2$year <- as.numeric(substr(tmp2$variable,2,5))
>
>   tmp2$month <- as.numeric(substr(tmp2$variable,7,8))
>
>   tmp2$day <- as.numeric(substr(tmp2$variable,10,11))
>
> }
>
>
> I have added date string and weekdays using the following codes:
>
> *Date Variable*
>
> tmp2$date <- with(tmp2, ymd(sprintf('%04d%02d%02d', year, month,
> day)))
>
>
>
> *Day Variable*
>
> tmp2$day <- weekdays(as.Date(tmp2$date))
>
>
>
> *Question 1: *Apologies for clipping the data and not showing all the
> hour blocks. I have included a full 8-hour block now. For each year, I have
> 3-hour blocks for every day.
>
>
>
> *Question 2: *I have two time periods for each GCM; 2026-2045 and
> 2081-2100. There are occasions when days would have 7 hour blocks instead
> of 8, it could be a data reporting issue from the models.
>
> *Reply to Comment 1: *The data has been downscaled to the US from gridded
> data; Resolution: T42 in atm. 1/3~1ºlat. x 1ºlon. tripolar grids in ocn. So
> Daylight Savings Time could well be an issue. I will look into it.
>
>
>
> *Reply to Comment 2: * Thank you for the suggestion. I realize that using
> rep() is risky, however, I have assign week numbers to each FIPS code (ID)
> for each year. I was thinking of something similar to this:
>
>
>
> weeks<-rep(seq(1,52,1),each=(unique(tmp2$FIPS)**8*)
>
>
>
> Any alternative code will be highly appreciated. My first goal is to
> subset/clip the data to begin on the first Monday and end on the last
> Sunday of each year.
>
> On Fri, Dec 19, 2014 at 9:37 AM, Jeff Newmiller 
> wrote:
>>
>> Thank you for attempting to convey your problem clearly using example
>> code... but your use of HTML email has very nearly undone all your
>> efforts.  Also, use of "dput" to make an R-readable block of data is more
>> reliable than read.table to get the data into our R sessions quickly.
>>
>> First question: you say these are three-hour results, but there are only
>> six per day in your example.
>>
>> Second question: you say the "number of blocks vary between 7 and 8", but
>> your example data does not illustrate that problem. (If there were 8 blocks
>> per day the three-hour statement would make more sense.)
>>
>> Comment: You have not mentioned timezone information... this information
>> looks ripe for GMT, but if that is a bad assumption then daylight savings
>> might account for some variations in blocks per day.
>>
>> Comment: I think your plan of using rep to identify the week numbers is
>> risky. I would recommend using a Date or POSIXt type to make the timestamps
>> computable, and then find the date corresponding to the beginning of the
>> week that the timestamp falls into. Then aggregate grouping on those time
>> values. Unfortunately, the specific way you go about identifying the
>> beginning of week may depend on the timezone information.
>>
>>
>> On Thu, 18 Dec 2014, Shouro Dasgupta wrote:

Re: [R] How to create a time series object with time only (no date)

2014-12-21 Thread Joshua Ulrich
On Sun, Dec 21, 2014 at 7:59 AM, ce  wrote:
>
> Thanks Joshua,
>
> Would you kindly explain if I have an xts array with different dates how I 
> change all dates to 1970-01-01 without touching the time ? I tried with 
> indexFormat without success. indexFormat(s) <- "1970-01-01 %H:%M:%S" . when I 
> plot a graph it still shows original dates.

You can't just change how the data are printed.  You have to actually
change the underlying data.  Here's one example of how you could do
that:

x <- xts(1:5, .POSIXct(1:5+86400*1:5, tz="UTC"))
index(x) <- as.POSIXct(paste("1970-01-01",
  format(index(x), "%H:%M:%S")), tz="UTC")

Note that you should ensure your timezone is UTC, GMT, or any timezone
that doesn't have daylight saving time.  Otherwise you might have
instances in your data where certain hours either do not exist or
exist twice.

> ce
>
>
> -Original Message-
> From: "Joshua Ulrich" [josh.m.ulr...@gmail.com]
> Date: 12/21/2014 12:53 AM
> To: "ce" 
> CC: "R-Help" 
> Subject: Re: [R] How to create a time series object with time only (no date)
>
> On Dec 20, 2014 11:11 PM, "ce"  wrote:
>>
>>
>> Dear all,
>>
>> I want to create a time series object from 00:00:00 to 23:59:00 without 
>> dates ?
>> I cant figure it out with xts ?
>>
> You cant create an xts object without a date in the index. If the date doesnt 
> matter, you can just set it to 1970-01-01 (or any other day).
>> ce
>>
>> __
>> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
>> 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.
>
>



-- 
Joshua Ulrich  |  about.me/joshuaulrich
FOSS Trading  |  www.fosstrading.com

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] How to create a time series object with time only (no date)

2014-12-21 Thread ce

Thanks Joshua,

Would you kindly explain if I have an xts array with different dates how I 
change all dates to 1970-01-01 without touching the time ? I tried with 
indexFormat without success. indexFormat(s) <- "1970-01-01 %H:%M:%S" . when I 
plot a graph it still shows original dates. 
ce


-Original Message-
From: "Joshua Ulrich" [josh.m.ulr...@gmail.com]
Date: 12/21/2014 12:53 AM
To: "ce" 
CC: "R-Help" 
Subject: Re: [R] How to create a time series object with time only (no date)

On Dec 20, 2014 11:11 PM, "ce"  wrote:
>
>
> Dear all,
>
> I want to create a time series object from 00:00:00 to 23:59:00 without dates 
> ?
> I cant figure it out with xts ?
>
You cant create an xts object without a date in the index. If the date doesnt 
matter, you can just set it to 1970-01-01 (or any other day).
> ce
>
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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 -- To UNSUBSCRIBE and more, see
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] non negativity constraints if else function

2014-12-21 Thread Jeff Newmiller
Please keep the conversation on the list, Esra... I don't do personal 
tutoring online, and others may answer your questions sooner and better 
than I can.


Some comments below:

On Sat, 20 Dec 2014, Esra Ulasan wrote:


Hello,
Thank you very much for your concern. I have these codes for portfolio 
optimization solved by Lagrange function.
I want to put a non-negativity constraint on weights. If the element of weight 
vector is negative remove it from the vector and recalculate the other elements 
of weight vector again after removing negative ones. It is going to be a 
non-negativity constraint on weights. I am trying to do it but I could not 
solve it.
Thank you very much.
Esra

optimization<-function(returns) {


# what is "x"? From context, perhaps a matrix?
# also, you are abusing positional parameters by not putting x first.. 
# read ?colMeans
# also, "mean" is the name of a VERY common R function... using it as a 
# variable name is very confusing at best



 mean <- colMeans(na.rm=FALSE,x)


# what is "asset.names"?


 names(mean) <-assets.names


# 100L is a simpler way to write an integer, though it does not look 
# to me like M really needs to be an integer so much



 M <- as.integer(100)#nuber of ports on the eff.front.
 S <- cov(x)
 Rmax<- 0.01282409#max monthly return value
 Thetah <- solve(S)


# what is p?


 u <- rep(1,p)
 a<- matrix(rep(0,4), nrow=2)
 a[1,1] <- t(u)%*%Thetah%*%u
 a[1,2] <- t(mean)%*%Thetah%*%u
 a[2,1] <- a[1,2]
 a[2,2] <- t(mean)%*%Thetah%*%mean


# you might be interested in the ?det function


 d <- a[1,1]*a[2,2]-a[1,2]*a[1,2]
 f <- (Thetah%*%(a[2,2]*u-a[1,2]*mean))/d
 g <- (Thetah%*%(-a[1,2]*u+a[1,1]*mean))/d
 r <- seq(0, Rmax, length=M)
 w <- matrix((rep(0, p*M)), nrow=p)

 for(i in 1:M) { w[,i] = f+r[i]*g#portfolio weights
 if (w[,i] <0) {w[,i]=0} else {w[,i]=w[,i]}
   }
 rownames(w)=assets.names


# I think the following two lines can entirely replace the for loop 
#above


w <- f + r * g

 w=ifelse(w<0, 0, w)



# I am not familiar with this type of modification of matrices


 #removing zero rows from weight matrix
 row_sub = apply(w, 1, function(row) all(row !=0 ))
 w2=w[row_sub,]

 s <- sqrt( a[1,1]*((r - a[1,2]/a[1,1])^2)/d + 1/a[1,1] ) # variance of the 
frontier
 ss <- sqrt(diag(S))
 names(ss) <- assets.names
 minp <- c(sqrt(1/a[1,1]), a[1,2]/a[1,1])# risk-return values of MVP
 wminp <- f + (a[1,2]/a[1,1])*g  #weights(allocation of 
assets) of MVP
 names(wminp) <- assets.names
 sharpe <- c((sqrt(1/a[1,1])) / a[1,2]/a[1,1])  #maximum sharpe ratio
 tanp <- c(sqrt(a[2,2])/a[1,2], a[2,2]/a[1,2])  #risk-return values of 
tangency portfolio
 wtanp <- f+(a[2,2]/a[1,2])*g   #weight (allocation of 
assets) of tangency portfolio



# Not sure what you are hoping to achieve by concatenating some scalars 
# and vectors together for return from the function.



return(c(s, r, ss, mean, minp, tanp, wminp, wtanp, w))



# seem to be missing a brace "}"

Have you read the help file for the "optim" function? It includes some 
options for constrained optimization that might help you.


(no more comments)



20 Ara 2014 tarihinde 20:41 saatinde, Jeff Newmiller  
?unlar? yazd?:


"if weight element is negative set it to zero, else recalculate the weights 
again"... this seems like the only way out of this loop is to calculate negative 
weights... and since you set them to zero at that point they will all be zero when you 
are done. Somehow I doubt that is what you intended to say.

On the other hand, if you had provided us with example input data and result 
data then we would not have to guess what you want.

I suspect that loops are not really what you want at all... but I could be 
wrong.

---
Jeff NewmillerThe .   .  Go Live...
DCN:Basics: ##.#.   ##.#.  Live Go...
 Live:   OO#.. Dead: OO#..  Playing
Research Engineer (Solar/BatteriesO.O#.   #.O#.  with
/Software/Embedded Controllers)   .OO#.   .OO#.  rocks...1k
---
Sent from my phone. Please excuse my brevity.

On December 19, 2014 8:26:55 PM PST, Esra Ulasan  wrote:

Hello,

I have tried the solve the non-negativity constraint "if else function"
in R. But I have done something wrong because it still gives the same
solution. I want that, if weight element is negative set it to zero,
else recalculate the weights again. These are the codes:

for(i in 1:M){ w[,i] = f+r[i]*g  #portfolio weights
  for(i in 1:M){
if (w <0){w=0}else{w=w}
  }
}
If you help me I would be happy
Thank you
[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mai