Re: [R] Alternatives for explicit for() loops
Thanks a lot Jim and Boris for replying.
Sent from my iPhone
> On Nov 9, 2015, at 1:13 AM, jim holtman wrote:
>
> You need to take a close look at the function incomb that you are creating.
> I see what appears to be a constant value
> ("*(gamma((1/beta)+1))*((alpha)^(-(1/beta)))") being computed that you might
> only have to compute once before the function. You are also referencing many
> variables (m, LED, j, ...) that are not being passed in that are in the
> global environment; it is probably better to pass them in to the function. I
> am not sure what 'pbapply' is doing for you since I see this is new to the
> code that you first sent out.
>
> I would be good it you told us what the function is trying to do; you are
> showing us how you want to do it, not what you want to do. Are there other
> ways of doing it? If speed is your problem, then consider the "Rcpp" package
> and write the function is C++ which might be faster, but again, take a look
> at what you are doing to see if there are other ways. I don't have time to
> dig into the code, since there is a lack of comments, to understand why you
> are using, e.g., 'choose', 'prod', etc.). There are probably a lot of ways
> of speeding up the code, if could tell us what you want to accomplish.
>
>
> Jim Holtman
> Data Munger Guru
>
> What is the problem that you are trying to solve?
> Tell me what you want to do, not how you want to do it.
>
>> On Sun, Nov 8, 2015 at 4:48 AM, Maram SAlem
>> wrote:
>> Thanks all for replying.
>>
>> In fact I've used the the Rprof() function and found out that the incomb()
>> function (in my code above) takes about 80% of the time, but I didn't
>> figure out which part of the function is causing the delay. So I thought
>> that this may be due to the for() loops.
>> I MUST run this code for rather large values of n and m, so is there any way
>> that can help me do that without having to wait for more than three days to
>> reach an output. N.B. I'll have to repeat these runs for may be 70 or 80
>> times , and this means HUGE time
>>
>> I'd appreciate any sort of help.
>> Thanks in advance.
>>
>> Maram Salem
>>
>>> On 6 November 2015 at 16:54, jim holtman wrote:
>>> If you have code that is running for a long time, then take a small case
>>> that only runs for 5-10 minutes and turn on the RProfiler so that you can
>>> see where you are spending your time. In most cases, it is probably not
>>> the 'for' loops that are causing the problem, but some function/calculation
>>> you are doing within the loop that is consuming the time, and until you
>>> determine what section of code that is, is it hard to tell exactly what the
>>> problem is, much less the solution.
>>>
>>>
>>> Jim Holtman
>>> Data Munger Guru
>>>
>>> What is the problem that you are trying to solve?
>>> Tell me what you want to do, not how you want to do it.
>>>
On Wed, Nov 4, 2015 at 9:09 AM, Maram SAlem
wrote:
Hi Jim,
Thanks a lot for replying.
In fact I'm trying to run a simulation study that enables me to calculate
the Bayes risk of a sampling plan selected from progressively type-II
censored Weibull model. One of the steps involves evaluating the expected
test time, which is a rather complicated formula that involves nested
multiple summations where the counters of the summation signs are
dependent, that's why I thought of I should create the incomb() function
inside the loop, or may be I didn't figure out how to relate its arguments
to the ones inside the loop had I created it outside it. I'm trying to
create a matrix of all the possible combinations involved in the
summations and then use the apply() function on each row of that matrix.
The problem is that the code I wrote works perfectly well for rather small
values of the sample size,n, and the censoring number, m (for example,
n=8,m=4),but when n and m are increased (say, n=25,m=15) the code keeps on
running for days with no output. That's why I thoug!
ht I should try to avoid explicit loops as much as possible, so I did my best
in this regard but still the code takes too long to execute,(more than three
days), thus, i believe there must be something wrong.
Here's the full code:
library(pbapply)
f1 <- function(n, m) {
stopifnot(n > m)
r0 <- t(diff(combn(n-1, m-1)) - 1L)
r1 <- rep(seq(from=0, len=n-m+1), choose( seq(to=m-2, by=-1,
len=n-m+1), m-2))
cbind(r0[, ncol(r0):1, drop=FALSE], r1, deparse.level=0)
}
simpfun<- function (x,n,m,p,alpha,beta)
{
a<-factorial(n-m)/(prod((factorial(x)))*(factorial((n-m)- sum(x
b <- ((m-1):1)
c<- a*((p)^(sum(x)))*((1-p)^(((m-1)*(n-m))- sum(x%*%(as.matrix(b)
d <- n - cumsum(x) - (1:(m-1))
e<- n*(prod(d))*c
LD<-list()
for (i in 1:(m-1)) {
LD[[i]]<-seq(
Re: [R] Error in rep.int() invalid 'times' value
Thanks a lot Petr for ur reply and advice. Hope I 'd be able to minmize time as
much as possible.
Regards,
Maram Salem
Sent from my iPhone
> On Oct 21, 2015, at 9:28 AM, PIKAL Petr wrote:
>
> Hi
>
> Several options
>
> Coding in C+ or other similar language (you seem to be more familiar with
> them)
>
> Using debug and find out how your function behaves in steps
>
> Using function with smaller m, n with Rprof to see where the time is spent.
>
> I believe that coding in R like it was C+ is the way to hell. There is
> nothing wrong with cycles however if you compute something which can be
> computed easily by vectorised approach you loose efficiency.
>
> e.g.
>
> you compute this
> m<-5
>> for ( i in 1:m-1)
> +{
> +b[i]<- (m-i)
> +}
>> b
> [1] 4 3 2 1
>
> but you can achieve it by
>
> b <- ((m-1):1)
>
> Time comparison:
>
>> m<-1e5
>> system.time(for ( i in 1:m-1) {b[i]<- (m-i)})
> user system elapsed
> 11.610.00 12.11
>> system.time(bb <- ((m-1):1))
> user system elapsed
> 0 0 0
>> all.equal(b,bb)
> [1] TRUE
>
> So the time gain is huge only in this small computation.
>
> Cheers
> Petr
>
>> -Original Message-
>> From: R-help [mailto:r-help-boun...@r-project.org] On Behalf Of Maram
>> SAlem
>> Sent: Tuesday, October 20, 2015 6:33 PM
>> To: William Dunlap
>> Cc: r-help@r-project.org
>> Subject: Re: [R] Error in rep.int() invalid 'times' value
>>
>> Yes Indeed William. f1() works perfectly well and took only 30 secs to
>> execute f1(25,15), but I wonder if there is anyway to speed up the
>> execution of the rest of my code (almost seven hours now) ?
>>
>> Thanks for helping.
>>
>> Maram Salem
>>
>>> On 20 October 2015 at 18:11, William Dunlap wrote:
>>>
>>> f0 is essentially your original code put into a function, so expect
>> it
>>> to fail in the same way your original code did.
>>> f1 should give the same answer as f0, but it should use less memory
>>> and time.
>>> Bill Dunlap
>>> TIBCO Software
>>> wdunlap tibco.com
>>>
>>>
>>> On Tue, Oct 20, 2015 at 2:05 AM, Maram SAlem
>>>
>>> wrote:
Thanks William. I've tried the first code, ( with f0() ), but still
for n=25, m=15 , I got this:
> s<-f0(25,15)
Error in rep.int(rep.int(seq_len(nx), rep.int(rep.fac, nx)), orep)
>> :
invalid 'times' value
In addition: Warning message:
In rep.int(rep.int(seq_len(nx), rep.int(rep.fac, nx)), orep) :
NAs introduced by coercion to integer range
I don't know if this is related to the memory limits of my laptop,
or it doesn't have to do with the memory.
Any help on how to fix this error will be greatly appreciated.
Thanks All.
Maram Salem
On 15 October 2015 at 17:52, William Dunlap
>> wrote:
>
> Doing enumerative combinatorics with rejection methods rarely
>> works
> well. Try mapping your problem to the problem of choosing
> m-1 items from n-1. E.g., your code was
>
> f0 <- function(n, m) {
> stopifnot(n > m)
> D<-matrix(0,nrow=n-m+1,ncol=m-1)
> for (i in 1:m-1){
> D[,i]<-seq(0,n-m,1)
> }
> ED <- do.call(`expand.grid`,as.data.frame(D))
> ED<-unname(as.matrix(ED))
> lk<-which(rowSums(ED)<=(n-m))
> ED[lk,]
> }
>
> and I think the following does the same thing in much less space
>> by
> transforming the output of combn().
>
> f1 <- function(n, m) {
> stopifnot(n > m)
> r0 <- t(diff(combn(n-1, m-1)) - 1L)
> r1 <- rep(seq(from=0, len=n-m+1), choose( seq(to=m-2, by=-1,
> len=n-m+1), m-2))
> cbind(r0[, ncol(r0):1, drop=FALSE], r1, deparse.level=0) }
>
> The code for adding the last column is a bit clumsy and could
> probably
>>> be
> improved. Both f0 and f1 could also be cleaned up to work for
>> m<=2.
>
> See Feller vol. 1 or Benjamin's "Proofs that (really) count" for
> more on this sort of thing.
>
>
>
> Bill Dunlap
> TIBCO Software
> wdunlap tibco.com
>
> On Thu, Oct 15, 2015 at 7:45 AM, Maram SAlem
> >>>
> wrote:
>>
>> Dear All,
>>
>> I'm trying to do a simple task (which is in fact a tiny part of a
>>> larger
>> code).
>>
>> I want to create a matrix, D, each of its columns is a sequence
>> from 0
>>> to
>> (n-m), by 1. Then, using D, I want to create another matrix ED,
>> whose rows represent all the possible combinations of the
>> elements
>> of the columns
>>> of
>> D. Then from ED, I'll select only the rows whose sum is less than
>> or equal to (n-m), which will be called the matrix s. I used the
>> following code:
>>
>>> n=5
>>> m=3
>>> D<-matrix(0,nrow=n-m+1,ncol=m-1) for (i in 1:m-1)
>> + {
>> + D[,i]<-seq(0,n-m,1)
>> + }
>>> ED <- do.call(`expand.grid`,as.data.frame(D))
>>> ED<-as.matrix(ED)
>>
>>> lk<-which(row
Re: [R] Error in rep.int() invalid 'times' value
Thanks a lot for helping William. Will check the reference as well.
Sent from my iPhone
> On Oct 15, 2015, at 5:52 PM, William Dunlap wrote:
>
> Doing enumerative combinatorics with rejection methods rarely
> works well. Try mapping your problem to the problem of choosing
> m-1 items from n-1. E.g., your code was
> f0 <- function(n, m) {
>stopifnot(n > m)
>D<-matrix(0,nrow=n-m+1,ncol=m-1)
>for (i in 1:m-1){
> D[,i]<-seq(0,n-m,1)
>}
>ED <- do.call(`expand.grid`,as.data.frame(D))
>ED<-unname(as.matrix(ED))
>lk<-which(rowSums(ED)<=(n-m))
>ED[lk,]
> }
> and I think the following does the same thing in much less space by
> transforming the output of combn().
> f1 <- function(n, m) {
>stopifnot(n > m)
>r0 <- t(diff(combn(n-1, m-1)) - 1L)
>r1 <- rep(seq(from=0, len=n-m+1), choose( seq(to=m-2, by=-1, len=n-m+1),
> m-2))
>cbind(r0[, ncol(r0):1, drop=FALSE], r1, deparse.level=0)
> }
> The code for adding the last column is a bit clumsy and could probably be
> improved. Both f0 and f1 could also be cleaned up to work for m<=2.
>
> See Feller vol. 1 or Benjamin's "Proofs that (really) count" for more on this
> sort of thing.
>
>
>
> Bill Dunlap
> TIBCO Software
> wdunlap tibco.com
>
>> On Thu, Oct 15, 2015 at 7:45 AM, Maram SAlem
>> wrote:
>> Dear All,
>>
>> I'm trying to do a simple task (which is in fact a tiny part of a larger
>> code).
>>
>> I want to create a matrix, D, each of its columns is a sequence from 0 to
>> (n-m), by 1. Then, using D, I want to create another matrix ED, whose rows
>> represent all the possible combinations of the elements of the columns of
>> D. Then from ED, I'll select only the rows whose sum is less than or equal
>> to (n-m), which will be called the matrix s. I used the following code:
>>
>> > n=5
>> > m=3
>> > D<-matrix(0,nrow=n-m+1,ncol=m-1)
>> > for (i in 1:m-1)
>> + {
>> + D[,i]<-seq(0,n-m,1)
>> + }
>> > ED <- do.call(`expand.grid`,as.data.frame(D))
>> > ED<-as.matrix(ED)
>>
>> > lk<-which(rowSums(ED)<=(n-m))
>>
>> > s<-ED[lk,]
>>
>>
>> This works perfectly well. But for rather larger values of n and m (which
>> are not so large actually), the number of all possible combinations of the
>> columns of D gets extremely large giving me this error (for n=25, m=15):
>>
>> > ED <- do.call(`expand.grid`,as.data.frame(D))
>> Error in rep.int(rep.int(seq_len(nx), rep.int(rep.fac, nx)), orep) :
>> invalid 'times' value
>> In addition: Warning message:
>> In rep.int(rep.int(seq_len(nx), rep.int(rep.fac, nx)), orep) :
>> NAs introduced by coercion to integer range
>>
>>
>> Any help or suggestions will be greatly appreciated.
>>
>> Thanks,
>>
>> Maram Salem
>>
>> [[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] Last msg not sent to the list
Thanks a lot Ivan and Henrik. Maram Sent from my iPhone > On Oct 14, 2015, at 3:30 PM, Henrik Bengtsson > wrote: > > In addition, > > if you go to https://stat.ethz.ch/mailman/listinfo/r-help (which is in > the footer of every R-help message), you'll find a link to 'R-help > Archives' (https://stat.ethz.ch/pipermail/r-help/). On that latter > page, you'll see all messages that have been sent out to the list. > That will allow you to make sure your message went out. > > /Henrik > > On Wed, Oct 14, 2015 at 3:54 AM, Ivan Calandra > wrote: >> Maram, >> >> I have received both of your e-mails on this topic, so they made it to the >> list. >> There is the option " Receive your own posts to the list?" on the membership >> configuration website (https://stat.ethz.ch/mailman/options/r-help/). If it >> is checked to "no", that would explain why you didn't receive your own >> posts. >> >> As to why nobody answered, no idea. Try again? >> >> HTH, >> Ivan* >> * >> >> -- >> Ivan Calandra, PhD >> University of Reims Champagne-Ardenne >> GEGENAA - EA 3795 >> CREA - 2 esplanade Roland Garros >> 51100 Reims, France >> +33(0)3 26 77 36 89 >> ivan.calan...@univ-reims.fr >> https://www.researchgate.net/profile/Ivan_Calandra >> >> Le 14/10/15 12:38, marammagdysa...@gmail.com a écrit : >> >>> Dear All, >>> >>> My last mail entitled: "using the apply() family to evaluate nested >>> functions with common arguments" to the r-help list didn't reach me though >>> I've sent it 2 days ago. I've included my suggested code and asked about >>> some details to make it work. In addition, I haven't received any feedback >>> from the r-help that may be that mail had something wrong or needs some >>> approval or ... , as the ones I used to receive when I first applied to the >>> list. >>> Any idea why is that?! >>> Thanks. >>> Maram Salem >>> >>> Sent from my iPhone >>> __ >>> 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-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] Last msg not sent to the list
Dear All, My last mail entitled: "using the apply() family to evaluate nested functions with common arguments" to the r-help list didn't reach me though I've sent it 2 days ago. I've included my suggested code and asked about some details to make it work. In addition, I haven't received any feedback from the r-help that may be that mail had something wrong or needs some approval or ... , as the ones I used to receive when I first applied to the list. Any idea why is that?! Thanks. Maram Salem Sent from my iPhone __ 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] Fwd: Warning message with maxLik()
Sent from my iPhone
Begin forwarded message:
> From: Maram SAlem
> Date: July 21, 2015 at 11:40:56 PM GMT+2
> To: Arne Henningsen
> Cc: "r-help@r-project.org"
> Subject: Re: [R] Warning message with maxLik()
>
> Dear Arne,
>
> The elements of the theta vector are indeed strictly positive. I've just
> tried to use instead : lamda = log (theta), which means that theta = exp
> (lamda), so as to get rid of the log() function that appears in the
> log-likelihood and is causing the 50 warnings, but still the estimates I got
> for lamda and then those I got for theta (using theta=exp(lamda)) are
> irrelvant and their standard errors are infinite, which means that therer is
> still a problem that I can't yet figure out.
>
> Thanks,
> Maram
>
>> On 18 July 2015 at 08:01, Arne Henningsen wrote:
>> Dear Maram
>>
>> - Please do not start a new thread for the same issue but reply to
>> previous messages in this thread [1].
>>
>> - Please read my previous responses [1] more carefully, e.g. to use
>> "theta <- exp( param )" which guarantees that all elements of "theta"
>> are always positive.
>>
>> [1]
>> http://r.789695.n4.nabble.com/NaN-produced-from-log-with-positive-input-td4709463.html
>>
>> Best regards,
>> Arne
>>
>>
>>
>> 2015-07-18 2:46 GMT+02:00 Maram SAlem :
>> > Dear All,
>> > I'm trying to get the MLe for a certain distribution using maxLik ()
>> > function. I wrote the log-likelihood function as follows:
>> > theta <-vector(mode = "numeric", length = 3)
>> > r<- 17
>> > n <-30
>> >
>> > T<-c(7.048,0.743,2.404,1.374,2.233,1.52,23.531,5.182,4.502,1.362,1.15,1.86,1.692,11.659,1.631,2.212,5.451)
>> > C<-
>> > c(0.562,5.69,12.603,3.999,6.156,4.004,5.248,4.878,7.122,17.069,23.996,1.538,7.792)
>> > # The loglik. func.
>> > loglik <- function(param) {
>> > theta[1]<- param[1]
>> > theta[2]<- param[2]
>> > theta[3]<- param[3]
>> >
>> > l<-(r*log(theta[3]))+(r*log(theta[1]+theta[2]))+(n*theta[3]*log(theta[1]))+(n*theta[3]*log(theta[2]))+
>> > (-1*(theta[3]+1))*sum(log((T*(theta[1]+theta[2]))+(theta[1]*theta[2])))+
>> > (-1*theta[3]*sum(log((C*(theta[1]+theta[2]))+(theta[1]*theta[2]
>> > return(l)
>> > }
>> >
>> > then, I evaluated it at theta<- c(40,50,2)
>> >
>> > v<-loglik(param=theta)
>> > v
>> > [1] -56.66653
>> >
>> > I used this same log-likelihood function, once with analytic gradient and
>> > another time with numerical one, with the maxLik function, and in both
>> > cases I got the same 50 warning messages and an MLE which is completely
>> > unrealistic as per my applied example.
>> >
>> > a <- maxLik(loglik, gradlik, hesslik, start=c(40,50,2))
>> >
>> > where gradlik and hesslik are the analytic gradient and Hessian matrix,
>> > respectively, given by:
>> >
>> > U <- vector(mode="numeric",length=3)
>> > gradlik<-function(param = theta,n, T,C)
>> > {
>> > U <- vector(mode="numeric",length=3)
>> > theta[1] <- param[1]
>> > theta[2] <- param[2]
>> > theta[3] <- param[3]
>> > r<- 17
>> > n <-30
>> > T<-c(7.048,0.743,2.404,1.374,2.233,1.52,23.531,5.182,4.502,1.362,1.15,1.86,1.692,11.659,1.631,2.212,5.451)
>> > C<-
>> > c(0.562,5.69,12.603,3.999,6.156,4.004,5.248,4.878,7.122,17.069,23.996,1.538,7.792)
>> > U[1]<- (r/(theta[1]+theta[2]))+((n*theta[3])/theta[1])+(
>> > -1*(theta[3]+1))*sum((T+theta[2])/((theta[1]+theta[2])*T+(theta[1]*theta[2])))+
>> > (-1*(theta[3]))*sum((C+theta[2])/((theta[1]+theta[2])*C+(theta[1]*theta[2])))
>> > U[2]<-(r/(theta[1]+theta[2]))+((n*theta[3])/theta[2])+
>> > (-1*(theta[3]+1))*sum((T+theta[1])/((theta[1]+theta[2])*T+(theta[1]*theta[2])))+
>> > (-1*(theta[3]))*sum((C+theta[1])/((theta[1]+theta[2])*C+(theta[1]*theta[2])))
>> > U[3]<-(r/theta[3])+(n*log(theta[1]*theta[2]))+
>> > (-1)*sum(log((T*(theta[1]+theta[2]))+(theta[1]*theta[2])))+(-1)*sum(log((C*(theta[1]+theta[2]))+(theta[1]*theta[2])))
>> > return(U)
>> > }
>> > hesslik<-function(param=theta,n,T,C)
>> > {
>> > theta[1] <- param[1]
>> > theta[2] <- param[2]
>> > theta[3] <- param[3]
>> > r<- 17
>> > n <-30
>> > T<-c(7.048,0.743,2.404,1.374,2.233,1.52,23.531,5.182,4.502,1.362,1.15,1.86,1.692,11.659,1.631,2.212,5.451)
>> > C<-
>> > c(0.562,5.69,12.603,3.999,6.156,4.004,5.248,4.878,7.122,17.069,23.996,1.538,7.792)
>> > G<- matrix(nrow=3,ncol=3)
>> > G[1,1]<-((-1*r)/((theta[1]+theta[2])^2))+((-1*n*theta[3])/(theta[1])^2)+
>> > (theta[3]+1)*sum(((T+theta[2])/((theta[1]+theta[2])*T+(theta[1]*theta[2])))^2)+(
>> > theta[3])*sum(((C+theta[2])/((theta[1]+theta[2])*C+(theta[1]*theta[2])))^2)
>> > G[1,2]<-((-1*r)/((theta[1]+theta[2])^2))+
>> > (theta[3]+1)*sum(((T)/((theta[1]+theta[2])*T+(theta[1]*theta[2])))^2)+
>> > (theta[3])*sum(((C)/((theta[1]+theta[2])*C+(theta[1]*theta[2])))^2)
>> > G[2,1]<-G[1,2]
>> > G[1,3]<-(n/theta[1])+(-1)*sum(
>> > (T+theta[2])/((theta[1]+theta[2])*T+(theta[1]*theta[2])))+(-1)*sum((C+theta[2])/((theta[1]+theta[2])*C+(theta[1]*theta[2])))
>> > G[3,1]<-G[1,3]
>> > G[2,2]<-((-1*r)/((theta[1]+theta[2])^2))+((-1*n*theta[3])/(theta[2])^2)+
>> > (theta[3]+1)*sum(((T+theta[1
