Thanks a lot Petr.
I did exactly what you did and found that f1() works for n=25 and m=15. But
I just wanted to figure out how f1() works, as I used its output for this
larger code and It has been running for almost 5 hours now.
s<-f1(25,15)
simpfun<- function (x,n,m,p,alpha,beta)
{
a<-factorial(n-m)/(prod((factorial(x)))*(factorial((n-m)- sum(x))))
b<- vector(mode = "numeric", length = m-1)
for ( i in 1:m-1)
{
b[i]<- (m-i)
}
c<- a*((p)^(sum(x)))*((1-p)^(((m-1)*(n-m))-sum(b%*%x)))
d <-vector(mode = "numeric", length = m-1)
for (i in 1:m-1)
{
d[i]<- n- (sum(x[(1):(i)])) - i
}
e<- n*(prod(d))*c
LD<-list()
for (i in 1:(m-1))
{
LD[[i]]<-seq(0,x[i],1)
}
LD[[m]]<-seq(0,(n-m-sum(x)),1)
LED<-expand.grid (LD)
LED<-as.matrix(LED)
store1<-vector(mode= "numeric",length=nrow(LED))
h<- vector(mode= "numeric", length= (m-1))
lm<- vector(mode= "numeric", length= (m-1))
for (j in 1:length(store1) )
{
incomb<-function(x,alpha,beta) {
g<-((-1)^(sum(LED[j,])))*(gamma((1/beta)+1))*((alpha)^(-(1/beta)))
for (i in 1:(m-1))
{
h[i]<- choose(x[i],LED[j,i])
}
ik<-prod(h)*choose((n-m-sum(x)),LED[j,m])
for (i in 1:(m-1))
{
lm[i]<-(sum(LED[j,1:i])) + i
}
plm<-prod(lm)
gil<-g*ik/(plm)
hlm<-vector(mode= "numeric",length=(sum(LED[j,])+(m-1)))
dsa<-length(hlm)
for (i in 1:dsa)
{
ppp<- sum(LED[j,])+(m-1)
hlm[i]<-
(choose(ppp,i))*((-1)^(i))*((i+1)^((-1)*((1/beta)+1)))
}
shl<-gil*(sum(hlm)+1)
return (shl)
}
store1[j]<-incomb(x,alpha=0.2,beta=2)
}
val1<- sum(store1)*e
return(val1)
}
va<-apply(s,1,simpfun,n=25,m=15,p=0.3,alpha=0.2,beta=2)
EXP<-sum(va)
I don't know if something is wrong or this is normal, but I've used R
before with looping codes but it never took too long.
Any suggestions please?
Thanks in advance.
Maram Salem
On 20 October 2015 at 15:27, PIKAL Petr <petr.pi...@precheza.cz> wrote:
> Hi
>
>
>
> Why are you confused?
>
>
>
> > f0(5,3)
>
> [,1] [,2]
>
> [1,] 0 0
>
> [2,] 1 0
>
> [3,] 2 0
>
> [4,] 0 1
>
> [5,] 1 1
>
> [6,] 0 2
>
> > f1(5,3)
>
> [,1] [,2]
>
> [1,] 0 0
>
> [2,] 1 0
>
> [3,] 2 0
>
> [4,] 0 1
>
> [5,] 1 1
>
> [6,] 0 2
>
> >
>
>
>
> seems to give the same result.
>
>
>
> > res0<-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
>
>
>
> > res1<-f1(25,15)
>
>
>
> So f1 works for required m,n.
>
>
>
> William just used different approach to get the same result, which is
> indeed quite common in R.
>
>
>
> Cheers
>
> Petr
>
>
>
> *From:* Maram SAlem [mailto:marammagdysa...@gmail.com]
> *Sent:* Tuesday, October 20, 2015 1:50 PM
> *To:* PIKAL Petr
> *Cc:* r-help@r-project.org
>
> *Subject:* Re: [R] Error in rep.int() invalid 'times' value
>
>
>
> Thanks for the hint Petr.
>
>
>
> I'm just a little bit confused with the function f1(). could you please
> help me and insert comments within f1() to be able to relate it with f0()?
>
>
>
> Thanks a lot.
>
>
>
> Maram Salem
>
>
>
> On 20 October 2015 at 11:29, PIKAL Petr <petr.pi...@precheza.cz> wrote:
>
> Hi
>
> What about using the second function f1 which William suggested?
>
> Cheers
> Petr
>
>
>
> > -----Original Message-----
> > From: R-help [mailto:r-help-boun...@r-project.org] On Behalf Of Maram
> > SAlem
> > Sent: Tuesday, October 20, 2015 11:06 AM
> > To: William Dunlap; r-help@r-project.org
> > Subject: Re: [R] Error in rep.int() invalid 'times' value
> >
> > 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 <wdun...@tibco.com> 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
> > > <marammagdysa...@gmail.com>
> > > 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
> > >>
>
>
>
> ------------------------------
> Tento e-mail a jakékoliv k němu připojené dokumenty jsou důvěrné a jsou
> určeny pouze jeho adresátům.
> Jestliže jste obdržel(a) tento e-mail omylem, informujte laskavě
> neprodleně jeho odesílatele. Obsah tohoto emailu i s přílohami a jeho kopie
> vymažte ze svého systému.
> Nejste-li zamýšleným adresátem tohoto emailu, nejste oprávněni tento email
> jakkoliv užívat, rozšiřovat, kopírovat či zveřejňovat.
> Odesílatel e-mailu neodpovídá za eventuální škodu způsobenou modifikacemi
> či zpožděním přenosu e-mailu.
>
> V případě, že je tento e-mail součástí obchodního jednání:
> - vyhrazuje si odesílatel právo ukončit kdykoliv jednání o uzavření
> smlouvy, a to z jakéhokoliv důvodu i bez uvedení důvodu.
> - a obsahuje-li nabídku, je adresát oprávněn nabídku bezodkladně přijmout;
> Odesílatel tohoto e-mailu (nabídky) vylučuje přijetí nabídky ze strany
> příjemce s dodatkem či odchylkou.
> - trvá odesílatel na tom, že příslušná smlouva je uzavřena teprve
> výslovným dosažením shody na všech jejích náležitostech.
> - odesílatel tohoto emailu informuje, že není oprávněn uzavírat za
> společnost žádné smlouvy s výjimkou případů, kdy k tomu byl písemně zmocněn
> nebo písemně pověřen a takové pověření nebo plná moc byly adresátovi tohoto
> emailu případně osobě, kterou adresát zastupuje, předloženy nebo jejich
> existence je adresátovi či osobě jím zastoupené známá.
>
> This e-mail and any documents attached to it may be confidential and are
> intended only for its intended recipients.
> If you received this e-mail by mistake, please immediately inform its
> sender. Delete the contents of this e-mail with all attachments and its
> copies from your system.
> If you are not the intended recipient of this e-mail, you are not
> authorized to use, disseminate, copy or disclose this e-mail in any manner.
> The sender of this e-mail shall not be liable for any possible damage
> caused by modifications of the e-mail or by delay with transfer of the
> email.
>
> In case that this e-mail forms part of business dealings:
> - the sender reserves the right to end negotiations about entering into a
> contract in any time, for any reason, and without stating any reasoning.
> - if the e-mail contains an offer, the recipient is entitled to
> immediately accept such offer; The sender of this e-mail (offer) excludes
> any acceptance of the offer on the part of the recipient containing any
> amendment or variation.
> - the sender insists on that the respective contract is concluded only
> upon an express mutual agreement on all its aspects.
> - the sender of this e-mail informs that he/she is not authorized to enter
> into any contracts on behalf of the company except for cases in which
> he/she is expressly authorized to do so in writing, and such authorization
> or power of attorney is submitted to the recipient or the person
> represented by the recipient, or the existence of such authorization is
> known to the recipient of the person represented by the recipient.
>
[[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.
