Re: [R] Loop Help

2012-06-17 Thread Jim Lemon

On 06/17/2012 08:50 AM, @ngel wrote:

Hello again and thank you all so much for the help. Well, I tried and did it
myself, here's what I wrote:

library(tnet)
net<- read.table("data.txt")
net<- as.tnet(net, type="longitudinal tnet")

loop<- range(net[,1])
net.static<- vector(length=as.integer(loop[2]-loop[1])+1, mode="list")
j<- 1
while(loop[1]<  loop[2]) {
net.static[[j]]<- as.static.tnet(net[net[,"t"]<=loop[1],])
degree<- degree_w(net.static[[j]])
results[[j]]<- apply(degree, 2, mean)
print(results[[j]])
loop[1]<- loop[1]+(30*30*24)
j<- j+1
}
max(loop[2]-loop[1])

and... it actually works, BUT
for every screenshot, it doesn't give the mean once, it prints it as many
times as the raws of the screenshot!
So if one screenshot has 3 lines, I get the result three times and then the
following screenshot's result as many times as the raws it has! Trying to
figure out why! Can somebody please help me?

By the way, I did tried this:

get_network_slices<-function(timestamps) {
  nslices<-length(timestamps)
  netslices<-list()
  for(slice in 1:nslices) {
  this_slice<- as.static.tnet(net[net[,"t"]<=as.POSIXlt(timestamps[slice],])
  degree_w(this_slice)
  netslices[[slice]]<-list(degree_w,apply(this_slice,2,mean))
  }
  return(netslices)
}

and got an error that this_slice is not defined!!!


Hi @ngel,
The error probably stems from me not getting the line that assigns a 
value to "this_slice" correct, so it fails and "this_slice" is not 
defined. I don't really know what that command is doing, and I have 
probably not gotten it quite right. Try just that line:


this_slice<-
 as.static.tnet(net[net[,"t"]<=as.POSIXlt(timestamps[slice],])

outside of the function, substituting a real value for 
"timestamps[slice]" and see what happens.


Jim

__
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] how to convert a vector to an upper matrix

2012-06-17 Thread Berend Hasselman

ozzi wrote
> 
> Dear Berend,
> 
> the coding is not working for more than 6. Try to plug 10 or other number.
> Would you please correct ?
> 

My mistake.

mkupper <- function(A) {
n <- length(A)
m <- (1+sqrt(1+8*n))/2
stopifnot(m==as.integer(m))
print(m)
L <- matrix(0,nrow=m,ncol=m)
L[lower.tri(L)] <- A

U <- t(L)
U
}

A <- 1:6
mkupper(A)


--
View this message in context: 
http://r.789695.n4.nabble.com/how-to-convert-a-vector-to-an-upper-matrix-tp3587715p4633630.html
Sent from the R help mailing list archive at Nabble.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] noob requesting help

2012-06-17 Thread Rui Barradas

Hello,

It's stringsAsFactors = FALSE, just one '='.

Your data structure, a list, seems to have elements, vectors, of the 
same lengths, all of them are of length 34773. You can check it with


sapply(dat, length)

It should return 4 times the value 34773. If so, there's no reason for 
David's sugestion not to work. If it does work, use dput() now, after 
converting 'dat' to a data.frame. (We still don't have a data example.)


Good luck, hope it works this time,

Rui Barradas

Em 16-06-2012 19:56, capital_P escreveu:


David Winsemius wrote



Try

dat2 <- data.frame(do.call(cbind, dat), stringsAsFactors=FALSE)


Use instead:

dat2 <- data.frame( dat, stringsAsFactors=FALSE)



Both do not seem to work:


dat2 <- data.frame(do.call(cbind, dat), stringsAsFactors=FALSE)

Error: cannot allocate vector of size 136 Kb
In addition: Warning messages:
1: In function (..., deparse.level = 1)  :
   number of rows of result is not a multiple of vector length (arg 3)
2: In structure(list(message = as.character(message), call = call),  :
   Reached total allocation of 2047Mb: see help(memory.size)
3: In structure(list(message = as.character(message), call = call),  :
   Reached total allocation of 2047Mb: see help(memory.size)




dat2 <- data.frame(dat, stringsAsFactors == FALSE)

Error in data.frame(dat, stringsAsFactors == FALSE) :
   object 'stringsAsFactors' not found

I do feel like this is going in the right direction. Thanks for your help so
far.

--
View this message in context: 
http://r.789695.n4.nabble.com/noob-requesting-help-tp4632803p4633607.html
Sent from the R help mailing list archive at Nabble.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.



__
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] visualize correlation matrix

2012-06-17 Thread Christof Kluß
Hi

is there a easy way to get something like
http://addictedtor.free.fr/graphiques/graphcode.php?graph=137

pairs(USJudgeRatings[,c(2:3,6,1,7)],
  lower.panel=panel.smooth, upper.panel=panel.cor)

but without the "lower.panel" (and without numbers and ticks at the border.)

thx
Christof

__
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] visualize correlation matrix

2012-06-17 Thread Rui Barradas

Hello,

The answer is yes, there is an easy way of getting that graphic without 
the lower.panel and axis tickmarks.



op <- par(xaxt="n", yaxt="n")
pairs(USJudgeRatings[,c(2:3,6,1,7)],
   upper.panel=panel.cor)
par(op)


where panel.cor() is the function in the link.

Hope this helps,

Rui Barradas

Em 17-06-2012 10:06, Christof Kluß escreveu:

Hi

is there a easy way to get something like
http://addictedtor.free.fr/graphiques/graphcode.php?graph=137

pairs(USJudgeRatings[,c(2:3,6,1,7)],
   lower.panel=panel.smooth, upper.panel=panel.cor)

but without the "lower.panel" (and without numbers and ticks at the border.)

thx
Christof

__
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] How to specify "newdata" in a Cox-Modell with a time dependent interaction term?

2012-06-17 Thread Jürgen Biedermann
Dear John,

Thank you very much for your help! It was very important for getting 
along further.

I found out some additional things which I want to give back to you as 
feedback (maybe other persons will have the same problems).

It's also possible to leave away the individual=T argument.
Instead you can use:

newdata.1<- 
data.frame(unique(Rossi.2[c("start","stop")]),fin=0,age=30,prio=5,Id=1,arrest.time=0)
fit <- survfit(mod.allison.5, newdata.1, id=Id) *but not* fit <- 
survfit(mod.allison.5, newdata.1, id="Id")

so you have to leave away the "" in specifying the id variable --> maybe 
the help file could be more precise here.

Additionally I found out some more things while dealing with my original 
data.

It's very important to specify the Surv object and the formula directly 
in the coxph function, *so don't do this:*

SurvO <- Surv(Rossi.2$start, Rossi.2$stop, Rossi.2$arrest.time)
form <- as.formula(SurvO ~ fin + age + age:stop + prio)
mod.allison.5<- coxph(form,
   data=Rossi.2)
mod.allison.5

newdata.1<- 
data.frame(unique(Rossi.2[c("start","stop")]),fin=0,age=30,prio=5,arrest.time=0)
fit <- survfit(mod.allison.5, newdata.1, individual=T)
 > Fehler in model.frame.default(data = newdata.1, formula = SurvO ~ fin 
+  :
   Variablenlängen sind unterschiedlich (gefunden für 'fin')
--> Error, length of variables are different

I think the error message is not optimal here, anyway.

The problem obviously consists in the survfit-function not recognizing 
the names of the variables anymore.
It seems completely logic, but anyway, it took me many hours.

Thanks again,
Jürgen



- Korrespondenz --

Betreff: Re: How to specify "newdata" in a Cox-Modell with a time 
dependent interaction term?
Von: John Fox 
An: J?rgen Biedermann 
Datum: 16.06.2012 17:22
> Dear Jurgen,
>
>> fit<- survfit(mod.allison.5, newdata.1, individual=TRUE)
>> fit
> Call: survfit(formula = mod.allison.5, newdata = newdata.1, individual = TRUE)
>
> records   n.max n.start  events  median 0.95LCL 0.95UCL
>19809 432 432 114  NA  NA  NA
>
> plot(fit)  # shows one survival curve
>
> See ?survfit.coxph.
>
> Best,
>   John
>
>
> On Sat, 16 Jun 2012 16:23:41 +0200
>   Jürgen Biedermann  wrote:
>> Dear John,
>>
>> Thank you very much for the quick answer. I will use the mentioned citation 
>> of your book in the future.
>>
>> However:
>> The Id argument is necessary for the function to know that the different 
>> rows belong to one individual.
>>
>> Or how the Help-File says:
>>
>> id:
>> optional variable name of subject identifiers. If this is present, then each 
>> group
>> of rows with the same subject id represents the covariate path through time 
>> of
>> a single subject, and the result will contain one curve per subject. If the 
>> coxph
>> fit had strata then that must also be specified in newdata. If missing, then 
>> each
>> individual row of newdata is presumed to represent a distinct subject and 
>> there
>> will be nrow(newdata) times the number of strata curves in the result (one 
>> for
>> each strata/subject combination).
>>
>> I tried your proposal anyway, but, as expected, fifty different survival 
>> curves appeared.
>>
>> So, there is still no solution.
>>
>> Best regards
>> Jürgen
>>
>> -- 
>> ---
>> Jürgen Biedermann
>> Bergmannstraße 3
>> 10961 Berlin-Kreuzberg
>> Mobil: +49 176 247 54 354
>> Home: +49 30 940 53 044
>> e-mail: juergen.biederm...@gmail.com
>>
>>
>> - Korrespondenz --
>>
>> Betreff: Re: How to specify "newdata" in a Cox-Modell with a time dependent 
>> interaction term?
>> Von: John Fox
>> An: J?rgen Biedermann
>> Datum: 16.06.2012 15:55
>>> Dear Jürgen,
>>>
>>> All the values of your Id variable are equal to 1; how can that make sense? 
>>> Removing the argument Id=id may get you what you want.
>>>
>>> By the way, the document to which you refer was an appendix to a 2002 book 
>>> that has been superseded by Fox and Weisberg, An R Companion to Applied 
>>> Regression, Second Edition (Sage, 2011). Appendices for that book are 
>>> available 
>>> at.
>>>
>>> I hope this helps,
>>>John
>>>
>>> 
>>> John Fox
>>> Sen. William McMaster Prof. of Social Statistics
>>> Department of Sociology
>>> McMaster University
>>> Hamilton, Ontario, Canada
>>> http://socserv.mcmaster.ca/jfox/
>>> 
>>> On Sat, 16 Jun 2012 15:04:48 +0200
>>>Jürgen Biedermann   wrote:
 Dear Mr. Therneau, Mr. Fox, or to whoever, who has some time...

 I don't find a solution to use the "survfit" function (package: survival)  
 for a defined pattern of covariates with a Cox-Model including a time 
 dependent interaction term. Somehow the definition of my "newdata" 
 argument seems to be erroneous.
 I already googled the problem, found many persons having the same or a 
 similar problem, but still

Re: [R] count data without NA in certain time intervals and plot it

2012-06-17 Thread Tagmarie
Thank you Arun for your time!
Your idea is maybe only the first step to what I want but it was
nevertheless a new tool for me and interessing to learn. 

I added a "week"-column to your data set: 
dattrial<-data.frame(a=c(1,NA,rnorm(4,10)), Week=c(3,3,3,4,4,4))

I am looking for a way to count the number of rows for each week which do
contain data (without NA).
In the next step I want to create a graph which shows the week on the x-axis
and the counted number of data for each week on the y-axis. 

Thank you! 

--
View this message in context: 
http://r.789695.n4.nabble.com/count-data-without-NA-in-certain-time-intervals-and-plot-it-tp4633611p4633635.html
Sent from the R help mailing list archive at Nabble.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] count data without NA in certain time intervals and plot it

2012-06-17 Thread Rui Barradas

Hello,

I've seen your reply to arun's reply and gave it a try.
Since arun's code included more than one column, I've added another in 
one of the examples.


# Example 1
dattrial1 <- data.frame(a=c(1,NA,rnorm(4,10)), Week=c(3,3,3,4,4,4))

d1 <- split(dattrial1, dattrial1$Week)
count <- sapply(d1, function(x) sum(!is.na(x$a)))
count


# Example 2
dattrial2 <- data.frame(a=c(1,NA,rnorm(4,10)), b=c(1,2,NA,3,4,6), 
Week=c(3,3,3,4,4,4))


d2 <- split(dattrial2, dattrial2$Week)
count <- sapply(d2, function(x){
yes <- apply(x, 1, function(y) all(!is.na(y)))
sum(yes)
})
count


# Works for both examples
plot(names(count), count, type="b", col="red", pch=16)


Hope this helps,

Rui Barradas

Em 16-06-2012 21:11, Tagmarie escreveu:

Hello,
I'm quite new to R and still spend hours trying to figure out single things
so I hope nobody rolls his eyes over my question.

I have a data set over time and converted it to the POSTIXct format. I added
a column in the data set for the week and the month.

I try to get a plot which shows the weeks on the x-axis and the number of
datasets without NAs on the y-axis. That doesn't sound too difficult but I
can't figure it out.

Does anybody have an idea?

--
View this message in context: 
http://r.789695.n4.nabble.com/count-data-without-NA-in-certain-time-intervals-and-plot-it-tp4633611.html
Sent from the R help mailing list archive at Nabble.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.



__
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] how to convert a vector to an upper matrix

2012-06-17 Thread ozzi
Dear Berend,

the coding is not working for more than 6. Try to plug 10 or other number.
Would you please correct ?

--
View this message in context: 
http://r.789695.n4.nabble.com/how-to-convert-a-vector-to-an-upper-matrix-tp3587715p4633628.html
Sent from the R help mailing list archive at Nabble.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] count data without NA in certain time intervals and plot it

2012-06-17 Thread Tagmarie
Great! That works! 
Thank you Rui! 
I would have spent days (which I don't have left before handing my report
in) getting there by myself! 
Have a great rest-weekend!

--
View this message in context: 
http://r.789695.n4.nabble.com/count-data-without-NA-in-certain-time-intervals-and-plot-it-tp4633611p4633638.html
Sent from the R help mailing list archive at Nabble.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] Temporal disaggregation

2012-06-17 Thread stef salvez
you are right Jeff and sorry for this

I will try to explain what I want.

I have the following dataset

dat <- data.frame("country" = c(rep(1,4)),
"date" = c("23/11/08","28/12/08","25/01/09","22/02/09"),
"price" = c(2,3,4,5))

Normally,  prices are observed every  4 weeks (28 days). Each
observation that it is published is a 4-week average. In some cases,
though, I have jumps (35 days-see from "23/11/08" to "28/12/08").  So
some prices that are published are 5-week average

I want to interpolate the 4-week average prices to monthly prices
so as to have

dat <- data.frame("country" = c(rep(1,4)),
"date" = c("11/08","12/08","01/09","02/09"),
"price" = c(3,2,1,4))

where the "new" "price" = c(3,2,1,4) will contain the interpolated
prices. So by starting what I have , that is,  - "price" = c(2,3,4,5))
-I want to calculate monthly average prices via interpolation.

I hope to have been more explicit.

thank you and apologies for one more time,





On 6/17/12, Jeff Newmiller  wrote:
> stop repeating yourself. and stop asking us to do your work.
>
> Try reading the posting guide. Give us an example of the output you DO want.
> Show us the code you already have. Use the dput function to give us example
> data to work with. Learn to use the str function so you know what data you
> are really working with. Ask specific questions.
> ---
> 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.
>
>
>
> stef salvez  wrote:
>
>>Dear R users,
>>I have a panel data set (in MS excel)  on prices across countries and
>>time
>>
>>countrytime price
>> 1 "23/11/08"2
>>1   "28/12/08"   3
>>1"25/01/09"   4
>>1   "22/02/09"   5
>>1"29/03/09"  6
>>1  "26/04/09"   32
>>1  "24/05/09"   23
>>1  "28/06/09"   32
>>2   "26/10/08"45
>>2  "23/11/08" 46
>>2  "21/12/08"   90
>>2  "18/01/09"54
>>2  "15/02/09" 65
>>2   "16/03/09"   77
>>2  "12/04/09"7
>>2   "10/05/09"   6
>>
>>
>>As you can see,
>>
>>1)the start and end date of the time series for countries 1 and 2 are
>>different. For example, for country 1 the time series begins on
>>"23/11/08" while for country 2 the time series begins on "26-10-2008”.
>>
>>2)My data on prices are available every 28 days (or equivalently every
>>4
>>weeks). So, each observation is a 4-week average. But in some cases I
>>have jumps (35 days or 29 days instead of
>>28 days). For example from the above table we have such jumps: from
>>"28/12/08" to "28/12/08" , from 22/02/09" to "29/03/09", etc
>>
>>My goal is to have a unified sequence of dates across countries.
>>Otherwise I will not be able to do further data/econometric analysis,
>>Unless you have different suggestion, I want  to take what I have and
>>calculate monthly average prices and also report how many prices those
>>averages are based on. I suppose that I will still have gaps and may
>>well need to interpolate.
>>
>>Put differently, I want to interpolate the 4-week average prices to
>>monthly average prices.
>>The problem is also the jumps where I have 5 weeks in some cases and I
>>want to find the monthly average of it.
>> I do not want something like
>>
>> country  yearmon avg.price freq
>>11 Nov 2008 21
>>21 Dec 2008 31
>>31 Jan 2009 41
>>41 Feb 2009 51
>>51 Mar 2009 61
>>61 Apr 2009321
>>71 May 2009231
>>81 Jun 2009321
>>92 Oct 2008451
>>10   2 Nov 2008461
>>11   2 Dec 2008901
>>12   2 Jan 2009541
>>13   2 Feb 2009651
>>14   2 Mar 2009771
>>15   2 Apr 2009 71
>>16   2 May 2009 61
>>
>>
>>Please, I would be grateful to you if you could provide the exact code
>>for doing this
>>
>>thank you
>>
>>__
>>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-

Re: [R] noob requesting help

2012-06-17 Thread capital_P

Rui Barradas wrote
> 
> It's stringsAsFactors = FALSE, just one '='.
> 
> sapply(dat, length)
> 
> It should return 4 times the value 34773. 
> 
> use dput() 
> 

it worked!

> dat2 <- data.frame(dat, stringsAsFactors = FALSE)

> sapply(dat2, length)
device_info_serial   hour   time 
 34773  34773  34773 
tripID 
 34773 

> dput(head(dat2, 20))
structure(list(device_info_serial = c(121L, 121L, 121L, 121L, 
121L, 121L, 121L, 121L, 121L, 121L, 121L, 121L, 121L, 121L, 121L, 
121L, 121L, 121L, 121L, 121L), hour = c(10L, 11L, 11L, 11L, 11L, 
11L, 11L, 11L, 11L, 11L, 11L, 11L, 12L, 12L, 12L, 12L, 12L, 12L, 
12L, 13L), time = structure(c(1242896364, 1242896671, 1242897111, 
1242897399, 1242897686, 1242897978, 1242898281, 1242898589, 1242898916, 
1242899253, 1242899542, 1242899835, 1242901122, 1242901902, 1242902199, 
1242902489, 1242902807, 1242903096, 1242903402, 1242903708), class =
c("POSIXct", 
"POSIXt"), tzone = ""), tripID = c(3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L)), .Names =
c("device_info_serial", 
"hour", "time", "tripID"), row.names = c(NA, 20L), class = "data.frame")
> 

> departures <- lapply(split(dat2, list(dat2$device_info_serial,
> dat2$tripID)), function(x) x[x$time == min(x$time),])

Now I have a new problem:

I need to make a histogram of the hours (in this case, their time of
departure), for the whole population and for each bird seperately.

Before, I used this: 

hist(active$hour, breaks = 24, xlim=c(0,24))

and:

b121 <- active[which(active$device_info_serial == 121),]
b130 <- active[which(active$device_info_serial == 130),]
b132 <- active[which(active$device_info_serial == 132),]
b133 <- active[which(active$device_info_serial == 133),]
etc.

But now there are no departures$device_info_serial or departures$hour,
because departures is a list, not a data frame.

I've tried:

> dep <- data.frame(departures, stringsAsFactors = FALSE)
Error in data.frame(`121.1` = list(device_info_serial = integer(0), hour =
integer(0),  : 
  arguments imply differing number of rows: 0, 1, 2

> dep <- data.frame(departures)
Error in data.frame(`121.1` = list(device_info_serial = integer(0), hour =
integer(0),  : 
  arguments imply differing number of rows: 0, 1, 2

--
View this message in context: 
http://r.789695.n4.nabble.com/noob-requesting-help-tp4632803p4633637.html
Sent from the R help mailing list archive at Nabble.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] Temporal disaggregation

2012-06-17 Thread stef salvez
or put differently, I want to change the data frequency of the time
series to monthly

thanks

On 6/17/12, stef salvez  wrote:
> you are right Jeff and sorry for this
>
> I will try to explain what I want.
>
> I have the following dataset
>
> dat <- data.frame("country" = c(rep(1,4)),
> "date" = c("23/11/08","28/12/08","25/01/09","22/02/09"),
> "price" = c(2,3,4,5))
>
> Normally,  prices are observed every  4 weeks (28 days). Each
> observation that it is published is a 4-week average. In some cases,
> though, I have jumps (35 days-see from "23/11/08" to "28/12/08").  So
> some prices that are published are 5-week average
>
> I want to interpolate the 4-week average prices to monthly prices
> so as to have
>
> dat <- data.frame("country" = c(rep(1,4)),
> "date" = c("11/08","12/08","01/09","02/09"),
> "price" = c(3,2,1,4))
>
> where the "new" "price" = c(3,2,1,4) will contain the interpolated
> prices. So by starting what I have , that is,  - "price" = c(2,3,4,5))
> -I want to calculate monthly average prices via interpolation.
>
> I hope to have been more explicit.
>
> thank you and apologies for one more time,
>
>
>
>
>
> On 6/17/12, Jeff Newmiller  wrote:
>> stop repeating yourself. and stop asking us to do your work.
>>
>> Try reading the posting guide. Give us an example of the output you DO
>> want.
>> Show us the code you already have. Use the dput function to give us
>> example
>> data to work with. Learn to use the str function so you know what data
>> you
>> are really working with. Ask specific questions.
>> ---
>> 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.
>>
>>
>>
>> stef salvez  wrote:
>>
>>>Dear R users,
>>>I have a panel data set (in MS excel)  on prices across countries and
>>>time
>>>
>>>countrytime price
>>> 1 "23/11/08"2
>>>1   "28/12/08"   3
>>>1"25/01/09"   4
>>>1   "22/02/09"   5
>>>1"29/03/09"  6
>>>1  "26/04/09"   32
>>>1  "24/05/09"   23
>>>1  "28/06/09"   32
>>>2   "26/10/08"45
>>>2  "23/11/08" 46
>>>2  "21/12/08"   90
>>>2  "18/01/09"54
>>>2  "15/02/09" 65
>>>2   "16/03/09"   77
>>>2  "12/04/09"7
>>>2   "10/05/09"   6
>>>
>>>
>>>As you can see,
>>>
>>>1)the start and end date of the time series for countries 1 and 2 are
>>>different. For example, for country 1 the time series begins on
>>>"23/11/08" while for country 2 the time series begins on "26-10-2008”.
>>>
>>>2)My data on prices are available every 28 days (or equivalently every
>>>4
>>>weeks). So, each observation is a 4-week average. But in some cases I
>>>have jumps (35 days or 29 days instead of
>>>28 days). For example from the above table we have such jumps: from
>>>"28/12/08" to "28/12/08" , from 22/02/09" to "29/03/09", etc
>>>
>>>My goal is to have a unified sequence of dates across countries.
>>>Otherwise I will not be able to do further data/econometric analysis,
>>>Unless you have different suggestion, I want  to take what I have and
>>>calculate monthly average prices and also report how many prices those
>>>averages are based on. I suppose that I will still have gaps and may
>>>well need to interpolate.
>>>
>>>Put differently, I want to interpolate the 4-week average prices to
>>>monthly average prices.
>>>The problem is also the jumps where I have 5 weeks in some cases and I
>>>want to find the monthly average of it.
>>> I do not want something like
>>>
>>> country  yearmon avg.price freq
>>>11 Nov 2008 21
>>>21 Dec 2008 31
>>>31 Jan 2009 41
>>>41 Feb 2009 51
>>>51 Mar 2009 61
>>>61 Apr 2009321
>>>71 May 2009231
>>>81 Jun 2009321
>>>92 Oct 2008451
>>>10   2 Nov 2008461
>>>11   2 Dec 2008901
>>>12   2 Jan 2009541
>>>13   2 Feb 2009651
>>>14   2 Mar 2009771
>>>15   2 Apr 2009 71
>>>16   2 May 2009 61
>>>
>>>
>>>Please, I would be grateful to you if you could provide the exact code
>>>for doing this
>>>
>>>thank you
>>>
>>>__

[R] How to save my result with Loops

2012-06-17 Thread Cassie
Hi everyone, i'm not very good with R, and i know that my question is quite
simple.
I generated this vector for a sample of 100 and i want to know how to save
these as 
ez1, ez1, ez3...ez100.
I guess i have to do it by  "for", but i didn't succeed..
i put here my output

ezt
  [1] -0.05470228 -1.76356668 -1.68843463 -0.77283767  0.25077826
-0.65236833  1.05230056 -0.48900254 -1.03857310 -0.30378185  0.86384470 
1.35341065 -0.55714575
 [14] -0.31853348  0.48801617  1.98541092 -2.52174364  0.17085272 
0.63973657  1.07164064 -0.35393600  0.61847775  1.18536018  1.40576222 
0.88830675 -0.84600747
 [27] -0.80453071 -2.90893372  0.47658777 -0.74312903  0.51442825 
0.22481451 -0.33104807  0.63688003  1.75841693  0.34841081 -1.18367323
-0.80041424  0.14039505
 [40]  0.78883365  1.50628498  0.81714749 -1.50882506 -1.11457623
-0.11722597 -0.23157740 -0.50457286 -0.76091885  0.89006377 -0.33416350
-2.27666710 -0.71381032
 [53] -2.12997125 -0.88065165  0.49617389  1.06932037 -0.88977000
-1.11977792  0.33601881  0.50490347 -1.46019653  1.30649432 -1.51794234 
0.26189545  0.86140985
 [66]  1.31436956  2.63651459  0.77599055  0.51114807  0.39171407
-0.07196456 -0.60086534  0.70217126  1.31626039 -0.67222186 -0.17144670 
0.32945944 -0.36040439
 [79] -1.56054165 -2.73620770  0.73432896  0.37370346  0.20326237 
0.72310549 -1.65473285  0.03052296  2.18543356 -1.80562305  1.32452946 
0.41418175 -1.82803429
 [92] -2.42929015  0.71443886 -1.13828906  2.30465334 -0.72840312 
0.49937916  0.66331583 -0.05890777 -0.98056989
> 

thank you very much

--
View this message in context: 
http://r.789695.n4.nabble.com/How-to-save-my-result-with-Loops-tp4633640.html
Sent from the R help mailing list archive at Nabble.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] How to save my result with Loops

2012-06-17 Thread R. Michael Weylandt
Short answer: don't. R is built around vectors and lists. Trying to
put individual variates in their own "scalar" variables is almost
always a bad idea.

Long answer: you can do it with assign() but I won't tell you how in
light of my short answer.

Michael

On Sun, Jun 17, 2012 at 8:30 AM, Cassie  wrote:
> Hi everyone, i'm not very good with R, and i know that my question is quite
> simple.
> I generated this vector for a sample of 100 and i want to know how to save
> these as
> ez1, ez1, ez3...ez100.
> I guess i have to do it by  "for", but i didn't succeed..
> i put here my output
>
> ezt
>  [1] -0.05470228 -1.76356668 -1.68843463 -0.77283767  0.25077826
> -0.65236833  1.05230056 -0.48900254 -1.03857310 -0.30378185  0.86384470
> 1.35341065 -0.55714575
>  [14] -0.31853348  0.48801617  1.98541092 -2.52174364  0.17085272
> 0.63973657  1.07164064 -0.35393600  0.61847775  1.18536018  1.40576222
> 0.88830675 -0.84600747
>  [27] -0.80453071 -2.90893372  0.47658777 -0.74312903  0.51442825
> 0.22481451 -0.33104807  0.63688003  1.75841693  0.34841081 -1.18367323
> -0.80041424  0.14039505
>  [40]  0.78883365  1.50628498  0.81714749 -1.50882506 -1.11457623
> -0.11722597 -0.23157740 -0.50457286 -0.76091885  0.89006377 -0.33416350
> -2.27666710 -0.71381032
>  [53] -2.12997125 -0.88065165  0.49617389  1.06932037 -0.88977000
> -1.11977792  0.33601881  0.50490347 -1.46019653  1.30649432 -1.51794234
> 0.26189545  0.86140985
>  [66]  1.31436956  2.63651459  0.77599055  0.51114807  0.39171407
> -0.07196456 -0.60086534  0.70217126  1.31626039 -0.67222186 -0.17144670
> 0.32945944 -0.36040439
>  [79] -1.56054165 -2.73620770  0.73432896  0.37370346  0.20326237
> 0.72310549 -1.65473285  0.03052296  2.18543356 -1.80562305  1.32452946
> 0.41418175 -1.82803429
>  [92] -2.42929015  0.71443886 -1.13828906  2.30465334 -0.72840312
> 0.49937916  0.66331583 -0.05890777 -0.98056989
>>
>
> thank you very much
>
> --
> View this message in context: 
> http://r.789695.n4.nabble.com/How-to-save-my-result-with-Loops-tp4633640.html
> Sent from the R help mailing list archive at Nabble.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.

__
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] How to save my result with Loops

2012-06-17 Thread John Kane
It is not clear exactly what you want.  Are you saying that you want to create 
an object out of each number in that vector or do you just want to name each 
element in the vector?

It is best to supply the code that you used so that readers can see what you 
are doing.  

Also when you are providing data please use something like dput()

Let's say your vector is called vec.  You can type dput(vec) and get a data 
structure that you can copy and paste into your email and which a reader can 
directly paste into R.


John Kane
Kingston ON Canada


> -Original Message-
> From: letamo...@hotmail.com
> Sent: Sun, 17 Jun 2012 06:30:13 -0700 (PDT)
> To: r-help@r-project.org
> Subject: [R] How to save my result with Loops
> 
> Hi everyone, i'm not very good with R, and i know that my question is
> quite
> simple.
> I generated this vector for a sample of 100 and i want to know how to
> save
> these as
> ez1, ez1, ez3...ez100.
> I guess i have to do it by  "for", but i didn't succeed..
> i put here my output
> 
> ezt
>   [1] -0.05470228 -1.76356668 -1.68843463 -0.77283767  0.25077826
> -0.65236833  1.05230056 -0.48900254 -1.03857310 -0.30378185  0.86384470
> 1.35341065 -0.55714575
>  [14] -0.31853348  0.48801617  1.98541092 -2.52174364  0.17085272
> 0.63973657  1.07164064 -0.35393600  0.61847775  1.18536018  1.40576222
> 0.88830675 -0.84600747
>  [27] -0.80453071 -2.90893372  0.47658777 -0.74312903  0.51442825
> 0.22481451 -0.33104807  0.63688003  1.75841693  0.34841081 -1.18367323
> -0.80041424  0.14039505
>  [40]  0.78883365  1.50628498  0.81714749 -1.50882506 -1.11457623
> -0.11722597 -0.23157740 -0.50457286 -0.76091885  0.89006377 -0.33416350
> -2.27666710 -0.71381032
>  [53] -2.12997125 -0.88065165  0.49617389  1.06932037 -0.88977000
> -1.11977792  0.33601881  0.50490347 -1.46019653  1.30649432 -1.51794234
> 0.26189545  0.86140985
>  [66]  1.31436956  2.63651459  0.77599055  0.51114807  0.39171407
> -0.07196456 -0.60086534  0.70217126  1.31626039 -0.67222186 -0.17144670
> 0.32945944 -0.36040439
>  [79] -1.56054165 -2.73620770  0.73432896  0.37370346  0.20326237
> 0.72310549 -1.65473285  0.03052296  2.18543356 -1.80562305  1.32452946
> 0.41418175 -1.82803429
>  [92] -2.42929015  0.71443886 -1.13828906  2.30465334 -0.72840312
> 0.49937916  0.66331583 -0.05890777 -0.98056989
>> 
> 
> thank you very much
> 
> --
> View this message in context:
> http://r.789695.n4.nabble.com/How-to-save-my-result-with-Loops-tp4633640.html
> Sent from the R help mailing list archive at Nabble.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.


FREE ONLINE PHOTOSHARING - Share your photos online with your friends and 
family!
Visit http://www.inbox.com/photosharing to find out more!

__
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] How to save my result with Loops

2012-06-17 Thread Bert Gunter
On Sun, Jun 17, 2012 at 6:30 AM, Cassie  wrote:

> Hi everyone, i'm not very good with R, and i know that my question is quite
> simple.
>

So make an effort to educate yourself instead of badgering this list.
Please read "An Intro to R" -- which ships with R (help.start() or Help
Menu --> Manuals to find it) . -- or some other beginning R tutorial of
your choice. With a small  initial investment of work, enlightenment will
ensue, happiness will reign, your self-esteem will be enhanced, and world
peace will break out. ... Well, maybe not the last... But you'll be able to
turn on the light and won't have to grope around in the dark hoping someone
will shine a flashlight.

Cheers,
Bert


> I generated this vector for a sample of 100 and i want to know how to save
> these as
> ez1, ez1, ez3...ez100.
> I guess i have to do it by  "for", but i didn't succeed..
> i put here my output
>
> ezt
>  [1] -0.05470228 -1.76356668 -1.68843463 -0.77283767  0.25077826
> -0.65236833  1.05230056 -0.48900254 -1.03857310 -0.30378185  0.86384470
> 1.35341065 -0.55714575
>  [14] -0.31853348  0.48801617  1.98541092 -2.52174364  0.17085272
> 0.63973657  1.07164064 -0.35393600  0.61847775  1.18536018  1.40576222
> 0.88830675 -0.84600747
>  [27] -0.80453071 -2.90893372  0.47658777 -0.74312903  0.51442825
> 0.22481451 -0.33104807  0.63688003  1.75841693  0.34841081 -1.18367323
> -0.80041424  0.14039505
>  [40]  0.78883365  1.50628498  0.81714749 -1.50882506 -1.11457623
> -0.11722597 -0.23157740 -0.50457286 -0.76091885  0.89006377 -0.33416350
> -2.27666710 -0.71381032
>  [53] -2.12997125 -0.88065165  0.49617389  1.06932037 -0.88977000
> -1.11977792  0.33601881  0.50490347 -1.46019653  1.30649432 -1.51794234
> 0.26189545  0.86140985
>  [66]  1.31436956  2.63651459  0.77599055  0.51114807  0.39171407
> -0.07196456 -0.60086534  0.70217126  1.31626039 -0.67222186 -0.17144670
> 0.32945944 -0.36040439
>  [79] -1.56054165 -2.73620770  0.73432896  0.37370346  0.20326237
> 0.72310549 -1.65473285  0.03052296  2.18543356 -1.80562305  1.32452946
> 0.41418175 -1.82803429
>  [92] -2.42929015  0.71443886 -1.13828906  2.30465334 -0.72840312
> 0.49937916  0.66331583 -0.05890777 -0.98056989
> >
>
> thank you very much
>
> --
> View this message in context:
> http://r.789695.n4.nabble.com/How-to-save-my-result-with-Loops-tp4633640.html
> Sent from the R help mailing list archive at Nabble.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.
>



-- 

Bert Gunter
Genentech Nonclinical Biostatistics

Internal Contact Info:
Phone: 467-7374
Website:
http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm

[[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] FW: Multivariate Normal and loops

2012-06-17 Thread R. Michael Weylandt
It's requested to cc the list on our back-and-forths so that i) you
can get help faster (I was moving all yesterday and it's a Hallmark
holiday in the US today) since others may respond; ii) others can
benefit by reading the archives.



On Sun, Jun 17, 2012 at 4:34 AM, Giulia Motta  wrote:
>
>
> the problem is that i want to do something like
> http://www.eco.uc3m.es/~jgonzalo/teaching/PhDTimeSeries/gonzalojec94.pdf in
> page 13
> with theta= 0.5, sigma = 0.5, a1=0 a2=-1 beta= 1
> and so i have to give a name (ewt, ezt) at the erratic component of my
> multivariate normal
> 
> From: letamo...@hotmail.com
> To:
> Subject: RE: [R] Multivariate Normal and loops
> Date: Sun, 17 Jun 2012 09:03:25 +
>
>
> I'm very grateful to you!
> But i'm wondering how i can do this, creating a matrix in which i can
> numerate every element?
> I don't know if i'm clear (probably not).
> Thank you very much
>

Unfortunately, I'm not totally clear as to what you are looking for
here. You can add row and column names (labels) with code like

rownames(A) <- c("a","cat","dog")

and then subscript appropriately:

A["cat",]

will pull out all columns (the empty index after the comma) in the row
labelled "cat". This is all documented in the manual I suggested to
you.

If there's something else you mean, try being a little clearer (using
examples from other programming languages if you know them is always
good) and we'll see what we can do.

Best,
Michael

>
>> From: michael.weyla...@gmail.com
>> Date: Fri, 15 Jun 2012 10:34:05 -0500
>> Subject: Re: [R] Multivariate Normal and loops
>> To: letamo...@hotmail.com
>> CC: r-help@r-project.org
>>
>> I don't think you need loops at all actually :-) [The main thing to
>> working with R vs. C is "vectorization" -- embrace it now and life
>> will be much better]
>>
>> You could do something like this:
>>
>> library(MASS) # Provides a multivariate normal distribution
>>
>> mvrnorm(100, c(1,5), matrix(c(2,1,1,2), ncol = 2))
>>
>> This gives you 100 samples of a multivariate normal with mean vector
>> (1,5) and covariance matrix
>>
>> (2, 1,
>> 1, 2)
>>
>> Feel free to tweak as necessary and type ?mvrnorm to see more
>> documentation.
>>
>> Hopefully this gets you started.
>>
>> For more on vectorization, see the "Introduction to R" that comes with
>> R. You can get it by typing help.start() at your command prompt.
>>
>> Hope this helps,
>> Michael
>>
>> On Fri, Jun 15, 2012 at 9:26 AM, Cassie  wrote:
>> > Hi,
>> > i'm not english and i'm not very familiar to R, and i'm asking if you
>> > can
>> > help me.
>> > I'm wondering how to create a multivariate normal an then repeat this
>> > for a
>> > sample of T=1000, and the save this result.
>> > Thank you very much for your helping
>> >
>> > --
>> > View this message in context:
>> > http://r.789695.n4.nabble.com/Multivariate-Normal-and-loops-tp4633504.html
>> > Sent from the R help mailing list archive at Nabble.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.

__
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] strang behaviour of mice package

2012-06-17 Thread Weidong Gu
Hi,

In my case, that was due to the existence of highly correlated
variables. see loggedEvents in the help file of mice.

Weidong Gu

On Sat, Jun 16, 2012 at 3:33 PM, b04877054  wrote:
> I have the same problem. How did you end up solving yours?
>
> --
> View this message in context: 
> http://r.789695.n4.nabble.com/strang-behaviour-of-mice-package-tp3679224p4633609.html
> Sent from the R help mailing list archive at Nabble.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.

__
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] Multi-threads in R

2012-06-17 Thread Gary Dong
Dear R users,

I'm wonder if there is a easy way to make R use multi-CPUs on my computer.
My computer has four CPUs but R uses only one. Thanks.

Gary

[[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] noob requesting help

2012-06-17 Thread Rui Barradas

Hello,

Try

names(departures) <- NULL
dep <- data.frame(departures[ !sapply(departures, is.null) ])

You now have data.frame 'dep'.

Rui Barradas

Em 17-06-2012 12:19, capital_P escreveu:


Rui Barradas wrote


It's stringsAsFactors = FALSE, just one '='.

sapply(dat, length)

It should return 4 times the value 34773.

use dput()



it worked!


dat2 <- data.frame(dat, stringsAsFactors = FALSE)



sapply(dat2, length)

device_info_serial   hour   time
  34773  34773  34773
 tripID
  34773


dput(head(dat2, 20))

structure(list(device_info_serial = c(121L, 121L, 121L, 121L,
121L, 121L, 121L, 121L, 121L, 121L, 121L, 121L, 121L, 121L, 121L,
121L, 121L, 121L, 121L, 121L), hour = c(10L, 11L, 11L, 11L, 11L,
11L, 11L, 11L, 11L, 11L, 11L, 11L, 12L, 12L, 12L, 12L, 12L, 12L,
12L, 13L), time = structure(c(1242896364, 1242896671, 1242897111,
1242897399, 1242897686, 1242897978, 1242898281, 1242898589, 1242898916,
1242899253, 1242899542, 1242899835, 1242901122, 1242901902, 1242902199,
1242902489, 1242902807, 1242903096, 1242903402, 1242903708), class =
c("POSIXct",
"POSIXt"), tzone = ""), tripID = c(3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L)), .Names =
c("device_info_serial",
"hour", "time", "tripID"), row.names = c(NA, 20L), class = "data.frame")





departures <- lapply(split(dat2, list(dat2$device_info_serial,
dat2$tripID)), function(x) x[x$time == min(x$time),])


Now I have a new problem:

I need to make a histogram of the hours (in this case, their time of
departure), for the whole population and for each bird seperately.

Before, I used this:

hist(active$hour, breaks = 24, xlim=c(0,24))

and:

b121 <- active[which(active$device_info_serial == 121),]
b130 <- active[which(active$device_info_serial == 130),]
b132 <- active[which(active$device_info_serial == 132),]
b133 <- active[which(active$device_info_serial == 133),]
etc.

But now there are no departures$device_info_serial or departures$hour,
because departures is a list, not a data frame.

I've tried:


dep <- data.frame(departures, stringsAsFactors = FALSE)

Error in data.frame(`121.1` = list(device_info_serial = integer(0), hour =
integer(0),  :
   arguments imply differing number of rows: 0, 1, 2


dep <- data.frame(departures)

Error in data.frame(`121.1` = list(device_info_serial = integer(0), hour =
integer(0),  :
   arguments imply differing number of rows: 0, 1, 2

--
View this message in context: 
http://r.789695.n4.nabble.com/noob-requesting-help-tp4632803p4633637.html
Sent from the R help mailing list archive at Nabble.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.



__
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] Multi-threads in R

2012-06-17 Thread R. Michael Weylandt
Take a look at the parallel package which ships with all current versions of R. 

Michael

On Jun 17, 2012, at 11:39 AM, Gary Dong  wrote:

> Dear R users,
> 
> I'm wonder if there is a easy way to make R use multi-CPUs on my computer.
> My computer has four CPUs but R uses only one. Thanks.
> 
> Gary
> 
>[[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.


Re: [R] Multi-threads in R

2012-06-17 Thread Uwe Ligges



On 17.06.2012 19:04, R. Michael Weylandt  wrote:

Take a look at the parallel package which ships with all current versions of R.


which is the right answer for the body of the message. For the subject 
line: Take a look at multi-threaded BLAS which can be used with R.


Uwe Ligges


Michael

On Jun 17, 2012, at 11:39 AM, Gary Dong  wrote:


Dear R users,

I'm wonder if there is a easy way to make R use multi-CPUs on my computer.
My computer has four CPUs but R uses only one. Thanks.

Gary

[[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] Efficient distance calculation on big matrix

2012-06-17 Thread Stefan Evert

> I'm working on analyzing a large data set, lets asume that
> dim(Data)=c(1000,8700). I want to calculate the canberra distance
> between the columns of this matrix, and using a toy example ('test' is
> a matrix filled with random numbers 0-1):
> 
>> system.time(d<-as.matrix(dist(t(test), method = "canberra", diag = FALSE, 
>> upper = FALSE, p = 2)))
>user   system  elapsed
> 1417.7133.219 1421.144
> The system.time results also confuse me a bit, since 99% of the time
> is not system time but user time. What does that mean?

User time is the time that R spends working on your problem; system time refers 
to tasks done by the operating system, e.g. disk access, managing locks and, 
most importantly, swapping when you run out of RAM.  With multi-threading, 
system time can be much larger than the time that has actually elapsed.

> 
> Is there any way to calculate the distance which would take less time?

Well, one thing you can do is to get a faster computer. :-)  The command above 
takes only 670 seconds on my MacBook Pro (without multi-threading).

Calculating a distance matrix is an expensive computation.  In your example, R 
has to carry out (8700 * 8700 * 1000) / 2 = 37.8 billion floating point 
divisions.  With approx. 27 clock cycles per division (according to tables I've 
found on the Web), this takes at least 340 seconds even on a 3GHz CPU (and 
ignoring memory access, addition/subtraction, loops, etc.).

You can shave off some of the time if you implement the distance calculation in 
C, inline the code to avoid callback functions in loops, operate on columns of 
the matrix directly (which should be more cache-friendly than rows) and don't 
check for NA's, NaN's and other degenerate cases.

I've done just that in my experimental R package "wordspace", which isn't on 
CRAN yet:

> library(wordspace)
> A <- matrix(runif(8.7e6), 1000, 8700)

> system.time(d1 <- as.matrix(dist(t(A), method="canberra")))
   user  system elapsed 
669.207   2.724 669.305 

> system.time(d2 <- dist.matrix(A, method="canberra", byrow=FALSE))
   user  system elapsed 
250.534   0.784 250.301 

> all(d1 == d2)
[1] TRUE

If you aren't tied to Canberra distance, you can use a less expensive metric 
such as the Manhattan distance for an additional, more substantial speed boost:

> system.time(d3 <- dist.matrix(A, method="manhattan", byrow=FALSE))
   user  system elapsed 
 42.488   0.999  43.569 

This is still single-threaded, so you can run multiple of these calculations in 
parallel depending on how many cores your server has.

Hope this helps,
Stefan


PS: In case you'd like to give it a try yourself and aren't daunted by a 
complete lack of documentation:

svn checkout svn://scm.r-forge.r-project.org/svnroot/wordspace/pkg


[ ev...@linglit.tu-darmstadt.de | http://purl.org/stefan.evert ]

__
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] noob requesting help

2012-06-17 Thread David Winsemius


On Jun 17, 2012, at 7:19 AM, capital_P wrote:



Rui Barradas wrote


It's stringsAsFactors = FALSE, just one '='.

sapply(dat, length)

It should return 4 times the value 34773.

use dput()



it worked!


dat2 <- data.frame(dat, stringsAsFactors = FALSE)



sapply(dat2, length)

device_info_serial   hour   time
34773  34773  34773
   tripID
34773


dput(head(dat2, 20))

structure(list(device_info_serial = c(121L, 121L, 121L, 121L,
121L, 121L, 121L, 121L, 121L, 121L, 121L, 121L, 121L, 121L, 121L,
121L, 121L, 121L, 121L, 121L), hour = c(10L, 11L, 11L, 11L, 11L,
11L, 11L, 11L, 11L, 11L, 11L, 11L, 12L, 12L, 12L, 12L, 12L, 12L,
12L, 13L), time = structure(c(1242896364, 1242896671, 1242897111,
1242897399, 1242897686, 1242897978, 1242898281, 1242898589,  
1242898916,
1242899253, 1242899542, 1242899835, 1242901122, 1242901902,  
1242902199,

1242902489, 1242902807, 1242903096, 1242903402, 1242903708), class =
c("POSIXct",
"POSIXt"), tzone = ""), tripID = c(3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L)), .Names =
c("device_info_serial",
"hour", "time", "tripID"), row.names = c(NA, 20L), class =  
"data.frame")





departures <- lapply(split(dat2, list(dat2$device_info_serial,
dat2$tripID)), function(x) x[x$time == min(x$time),])


That's going to remove a lot of rows.



Now I have a new problem:

I need to make a histogram of the hours (in this case, their time of
departure), for the whole population and for each bird seperately.


Since you are not telling us what columns correspond to those  
concepts, I'm going to use the terms you are. Something along the  
lines of:


require(lattice)
histogram( ~ hours | bird,  data=dat2, nint =24 , xlim=c(0,24))
histogram( ~ hours ,data=dat2,  nint =24 , xlim=c(0,24))




Before, I used this:

hist(active$hour, breaks = 24, xlim=c(0,24))

and:

b121 <- active[which(active$device_info_serial == 121),]
b130 <- active[which(active$device_info_serial == 130),]
b132 <- active[which(active$device_info_serial == 132),]
b133 <- active[which(active$device_info_serial == 133),]
etc.

But now there are no departures$device_info_serial or departures$hour,
because departures is a list, not a data frame.



The result on that example was just a list with one element, a  
dataframe.


> str(departures)
List of 1
 $ 121.3:'data.frame':  1 obs. of  4 variables:
  ..$ device_info_serial: int 121
  ..$ hour  : int 10
  ..$ time  : POSIXct[1:1], format: "2009-05-21 04:59:24"
  ..$ tripID: int 3

So to get what you might have expected

> dep <- departures[[1]]
> dep
  device_info_serial hourtime tripID
1121   10 2009-05-21 04:59:24  3


But since you didn't actually describe what you did want (or if you  
did you left upthread), its only a guess.



I've tried:


dep <- data.frame(departures, stringsAsFactors = FALSE)
Error in data.frame(`121.1` = list(device_info_serial = integer(0),  
hour =

integer(0),  :
 arguments imply differing number of rows: 0, 1, 2


Learn to use str().




dep <- data.frame(departures)
Error in data.frame(`121.1` = list(device_info_serial = integer(0),  
hour =

integer(0),  :
 arguments imply differing number of rows: 0, 1, 2

--
View this message in context: 
http://r.789695.n4.nabble.com/noob-requesting-help-tp4632803p4633637.html
Sent from the R help mailing list archive at Nabble.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.


David Winsemius, MD
West Hartford, CT

__
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] how to delete a matrix column

2012-06-17 Thread 李红旺
> a<-c(1,4)
> a
[1] 1 4
> b<-a*5
> b
[1] 5 20

a is a very long vector , how can i get c(1:5,4:20)? i do not want to
use a loop.

thanks very much!

__
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] how to delete a matrix column

2012-06-17 Thread Jeff Newmiller
"a" looks pretty short to me, just two elements.

You should spend some (more) time reading the "Introduction to R" document that 
is supplied with the software... particularly the parts on indexing.
---
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.

"李红旺"  wrote:

>> a<-c(1,4)
>> a
>[1] 1 4
>> b<-a*5
>> b
>[1] 5 20
>
>a is a very long vector , how can i get c(1:5,4:20)? i do not want to
>use a loop.
>
>thanks very much!
>
>__
>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] how to delete a matrix column

2012-06-17 Thread Petr Savicky
On Mon, Jun 18, 2012 at 01:24:21AM +0800, 李红旺 wrote:
> > a<-c(1,4)
> > a
> [1] 1 4
> > b<-a*5
> > b
> [1] 5 20
> 
> a is a very long vector , how can i get c(1:5,4:20)? i do not want to
> use a loop.

Hi.

How large are the numbers in a? If max(a) is not too large, try
the following.

  a <- c(1, 4, 3, 1)
  ma <- max(a)
  b <- matrix(a, nrow=4*ma+1, ncol=length(a), byrow=TRUE) +
   matrix(0:(4*ma), nrow=4*ma+1, ncol=length(a))
  out <- c(b[b <= rep(5*a, each=4*ma+1)])
  out

  [1]  1  2  3  4  5  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20  3  4  5
  [26]  6  7  8  9 10 11 12 13 14 15  1  2  3  4  5

Hope this helps.

Petr Savicky.

__
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] time-series classification with k-nn

2012-06-17 Thread Pierre Antoine DuBoDeNa
Hello,

i am trying to apply k-nn classification for my time-series, however the
euclidean distance is not the best choice as the features i use are not all
normalized (others have values form 0-1 others are negative etc.) and also
it doesn't do any feature evaluation and give different weights to features.

Any idea if any package or code could be helpful? to learn distance metric
for k-nn ?

Also any way to compare with the "classic" distance metrics, like DTW etc.
? I am interested in learning the thresholds of similarity..

Best,
PA

[[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] Multi-threads in R

2012-06-17 Thread Gary Dong
Thanks for all replied.

I read the introduction of R parallel. Is it for loops only?

Gary


On Sun, Jun 17, 2012 at 10:04 AM, R. Michael Weylandt <
michael.weyla...@gmail.com>  wrote:

> Take a look at the parallel package which ships with all current versions
> of R.
>
> Michael
>
> On Jun 17, 2012, at 11:39 AM, Gary Dong  wrote:
>
> > Dear R users,
> >
> > I'm wonder if there is a easy way to make R use multi-CPUs on my
> computer.
> > My computer has four CPUs but R uses only one. Thanks.
> >
> > Gary
> >
> >[[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.


Re: [R] how to delete a matrix column

2012-06-17 Thread Ted Harding
On 17-Jun-2012 17:24:21 ÀîºìÍú wrote:
>> a<-c(1,4)
>> a
> [1] 1 4
>> b<-a*5
>> b
> [1] 5 20
> 
> a is a very long vector , how can i get c(1:5,4:20)?
> i do not want to use a loop.
> 
> thanks very much!

Your question is not quite clear, but I think you mean:

  a is a very long vector, and you want a vector with elements

  a[1],a[2],a[3],a[4],a[5],a[4],a[5],a[6],...,a[20]

If that is the case, then a[c(1:5,4:20)] will do it. Example:

a <- (1:100)
a[c(1:5,4:20)]
# [1]  1  2  3  4  5  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20

Does that help?

If not, try to re-formulate the question!
Ted.

-
E-Mail: (Ted Harding) 
Date: 17-Jun-2012  Time: 20:58:01
This message was sent by XFMail

__
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] count data without NA in certain time intervals and plot it

2012-06-17 Thread arun
Hi,

Sorry, I didn't understand your question in the first post.  I saw Rui's reply 
and your reply that it is solved.

I have another solution if it helps you.


dattrial<-data.frame(a=c(1,NA,rnorm(4,10)), Week=c(3,3,3,4,4,4))
dattrial_wk3<-subset(dattrial,Week==3)
dattrial_wk4<-subset(dattrial,Week==4)
count1<-colSums(!is.na(dattrial_wk3))
count2<-colSums(!is.na(dattrial_wk4))
dattrialnew<-data.frame(rbind(count1[1],count2[2]),Week=(rle(dattrial$Week)$values))
 plot(dattrialnew$Week,dattrialnew$a,type="l",col="blue",pch=14,xlab="Week",ylab="Count")


A.K.




- Original Message -
From: Tagmarie 
To: r-help@r-project.org
Cc: 
Sent: Sunday, June 17, 2012 5:40 AM
Subject: Re: [R] count data without NA in certain time intervals and plot it

Thank you Arun for your time!
Your idea is maybe only the first step to what I want but it was
nevertheless a new tool for me and interessing to learn. 

I added a "week"-column to your data set: 
dattrial<-data.frame(a=c(1,NA,rnorm(4,10)), Week=c(3,3,3,4,4,4))

I am looking for a way to count the number of rows for each week which do
contain data (without NA).
In the next step I want to create a graph which shows the week on the x-axis
and the counted number of data for each week on the y-axis. 

Thank you! 

--
View this message in context: 
http://r.789695.n4.nabble.com/count-data-without-NA-in-certain-time-intervals-and-plot-it-tp4633611p4633635.html
Sent from the R help mailing list archive at Nabble.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.


__
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] visualize correlation matrix

2012-06-17 Thread Kevin Wright
You can also try the corrgram package.

Kevin


On Sun, Jun 17, 2012 at 4:06 AM, Christof Kluß  wrote:
> Hi
>
> is there a easy way to get something like
> http://addictedtor.free.fr/graphiques/graphcode.php?graph=137
>
> pairs(USJudgeRatings[,c(2:3,6,1,7)],
>  lower.panel=panel.smooth, upper.panel=panel.cor)
>
> but without the "lower.panel" (and without numbers and ticks at the border.)
>
> thx
> Christof
>
> __
> 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.



-- 
Kevin Wright

__
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] Multi-threads in R

2012-06-17 Thread R. Michael Weylandt
I would argue (somewhat emphatically) that the parallel facilities you
are looking at are absolutely not for `for` loops. `for` loops are a
control structure native to imperative programming and, as such, are
inherently stateful. This provides many advantages, but easy
parallelization is absolutely not one of them. Rather consider R as a
functional language and the appropriate "control structure" is the
Map(), most famously of MapReduce, but also inherent in R's *apply()
family of functions. The careful reader will note that here is where
parallelization is exposed. By being state-free, there's no inherent
sequentiality and, consequently, no barriers to parallel computing.

For a somewhat concrete example consider what it would take to do a
state-dependent simulation (like a gibbs sampler) in parallel: it is,
to my mind, quite non-trivial. Compare this to the question of doing
column sums of a matrix in parallel -- effectively trivial. The
parallelization afforded by R is of that second sort.

A red herring in this discussion is the foreach() formalism, but a
moment's thought suggests to me that foreach() is much closer to Map()
than to for().

Best,
Michael

On Sun, Jun 17, 2012 at 1:23 PM, Gary Dong  wrote:
> Thanks for all replied.
>
> I read the introduction of R parallel. Is it for loops only?
>
> Gary
>
>
> On Sun, Jun 17, 2012 at 10:04 AM, R. Michael Weylandt
>   wrote:
>>
>> Take a look at the parallel package which ships with all current versions
>> of R.
>>
>> Michael
>>
>> On Jun 17, 2012, at 11:39 AM, Gary Dong  wrote:
>>
>> > Dear R users,
>> >
>> > I'm wonder if there is a easy way to make R use multi-CPUs on my
>> > computer.
>> > My computer has four CPUs but R uses only one. Thanks.
>> >
>> > Gary
>> >
>> >    [[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.


Re: [R] Multi-threads in R

2012-06-17 Thread Fabrice Tourre
I think multicore is one of answer if you can write your function in to lapply.

On Mon, Jun 18, 2012 at 12:14 PM, R. Michael Weylandt
 wrote:
> I would argue (somewhat emphatically) that the parallel facilities you
> are looking at are absolutely not for `for` loops. `for` loops are a
> control structure native to imperative programming and, as such, are
> inherently stateful. This provides many advantages, but easy
> parallelization is absolutely not one of them. Rather consider R as a
> functional language and the appropriate "control structure" is the
> Map(), most famously of MapReduce, but also inherent in R's *apply()
> family of functions. The careful reader will note that here is where
> parallelization is exposed. By being state-free, there's no inherent
> sequentiality and, consequently, no barriers to parallel computing.
>
> For a somewhat concrete example consider what it would take to do a
> state-dependent simulation (like a gibbs sampler) in parallel: it is,
> to my mind, quite non-trivial. Compare this to the question of doing
> column sums of a matrix in parallel -- effectively trivial. The
> parallelization afforded by R is of that second sort.
>
> A red herring in this discussion is the foreach() formalism, but a
> moment's thought suggests to me that foreach() is much closer to Map()
> than to for().
>
> Best,
> Michael
>
> On Sun, Jun 17, 2012 at 1:23 PM, Gary Dong  wrote:
>> Thanks for all replied.
>>
>> I read the introduction of R parallel. Is it for loops only?
>>
>> Gary
>>
>>
>> On Sun, Jun 17, 2012 at 10:04 AM, R. Michael Weylandt
>>   wrote:
>>>
>>> Take a look at the parallel package which ships with all current versions
>>> of R.
>>>
>>> Michael
>>>
>>> On Jun 17, 2012, at 11:39 AM, Gary Dong  wrote:
>>>
>>> > Dear R users,
>>> >
>>> > I'm wonder if there is a easy way to make R use multi-CPUs on my
>>> > computer.
>>> > My computer has four CPUs but R uses only one. Thanks.
>>> >
>>> > Gary
>>> >
>>> >    [[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] Cholesky decomposition error

2012-06-17 Thread nataraj
Thanks Kjetil for your detailed code for log-Cholesky but my situation is like 
optimization of the variables inside the matrix.

Let me change your own example matrix to explain the problem.

 A martrix is  [,1] [,2][,3]
[1,]3*Vs+Vn1*Vs 1*Vs
[2,]1*Vs 3*Vs+Vn1*Vs
[3,]1*Vs 1*Vs   3*Vs+Vn

So I have to optimize two variables "Vs" and "Vn" using the maximum likelihood 
function, which need the determinant and inverse of the matrix A to compute the 
value.

Isn't it logically correct that I have to seed some initial values for the Vs 
and Vn and use it for predicting the value of the matrix and check for its 
positive definite condition, if the condition fulfill, then I have to compute 
the determinant value and inverse of the matrix in order to use them in the 
maximum likelihood function and the optimization iteration to be carried out 
further until parameters converges.

I sincerely believe that you have points to implements log-Cholesky in my 
situation, but I am simply not able to see where is the Choleksy decomposition 
to be implemented in this my workflow.

Regards,
B.Nataraj




-Original Message-
From: Kjetil Halvorsen [mailto:kjetilbrinchmannhalvor...@gmail.com]
Sent: Sunday, June 17, 2012 4:10 AM
To: Nataraj B (ORLL-Biotech)
Cc: r-help@r-project.org
Subject: Re: [R] Cholesky decomposition error

see below.

On Fri, Jun 15, 2012 at 11:53 PM,   wrote:
> Dear Mr.Kjetil,
>
> Thanks for your comment. You have already pointed me the article in reply to 
> one of my earlier post to this list and I am following the paper. Now I am 
> checking for condition for positive definiteness for original matrix using a 
> simple script (got from earlier posting in the list) and if the condition 
> passed then the optimization function perform decomposition of the matrix 
> using chol() function.
>
> I believe all 5 parameterization techniques listed in the paper weigh each 
> other based on its performance and not on accuracy, please correct me if I am 
> wrong. Since I am novice, I just want to start with the easiest method 
> "Cholesky" for my purpose.
>
> Are you recommending log-Cholesky function , if so could you please tell me 
> the name of the function which does log-Cholesky decomposition of a matrix.

I dont think there is a prewritten function for log cholesy, but it is
very easy to write!

> makelogchol
function(A) {#A should be a square posdef matrix
m <- dim(A)[1]
uind <- upper.tri(A)
R <- chol(A)
diag <- log(diag(R))
upper <- R[uind]
return(list(diag=diag, upper=upper))
}

and its "inverse" --- to get the cholesky factor back:

> makemat
function(diag, upper) {
m <- length(diag)
mu <- m*(m-1)/2
if (length(upper)!=mu) stop("incompatible lengths")
A <- matrix(0.0, m, m)
ind <- upper.tri(A)
A[ind] <- upper
diag(A) <- exp(diag)
A
}

and an example of its use:

Lets say A is the posdef matrix where you will start your
optimization, so you need
the log cholesy parameters to feed to optim:

> A
 [,1] [,2] [,3]
[1,]311
[2,]131
[3,]113
> start <- makelogchol(A)
> start
$diag
[1] 0.5493061 0.4904146 0.4581454

$upper
[1] 0.5773503 0.5773503 0.4082483

Then to get the cholesky factors back you call makemat:

> R <- makemat(start$diag,  start$upper)
> R
 [,1]  [,2]  [,3]
[1,] 1.732051 0.5773503 0.5773503
[2,] 0.00 1.6329932 0.4082483
[3,] 0.00 0.000 1.5811388
> t(R) %*% R
 [,1] [,2] [,3]
[1,]311
[2,]131
[3,]113
>

Kjetil







>
> Regards,
> B.Nataraj
>
>
>
> -Original Message-
> From: Kjetil Halvorsen [mailto:kjetilbrinchmannhalvor...@gmail.com]
> Sent: Friday, June 15, 2012 11:09 PM
> To: Nataraj B (ORLL-Biotech)
> Cc: gunter.ber...@gene.com; r-help@r-project.org
> Subject: Re: [R] Cholesky decomposition error
>
> see inline.
>
> On Fri, Jun 15, 2012 at 4:33 AM,   wrote:
>> Thanks for your reply. I am sorry and I am bit hurried up to say before 
>> doing a proper due diligence, I have found out that during the optimization 
>> the variables tend to vary the values of the matrix , the function report 
>> error at some point (in particular iteration step) when the matrix become 
>> non-decomposable due to not a positive definiteness.
>
> ¿What are you trying to optimize? If the objective function depends on
> a matrix argument
> which has to be a positive definite function, you must parametrize the
> matrix such that the matrix
> inside the optimizer always is positive definite. So if your positive
> definite matrix is A, then, for example, represent it as
> its cholesky decomposition A= L L^T where L is lower triangular with
> positive diagonal. Here the stricly
> upper diagonal part varies freely, but the diagonal not, so represent
> the diagonal as exp( l_i)
> where now the l_i varies freely. This is called the log-Cholesky
> parametrization. For other ideas along this lines,