Re: [R] The code itself disappears after starting to execute the for loop
Thanks a lot Bert , will check out your suggestions. I've unchecked the buffer output option in GUI but still have the same problem. Thanks for your time and concern. Maram Salem Sent from my iPhone > On Nov 23, 2016, at 5:55 PM, Bert Gunter <bgunter.4...@gmail.com> wrote: > > In addition to Jim's comments, which you have not yet satisfactorily > addressed (buffering in GUI??), > > 1. Show your code! > > 2. Show ouput of sessionInfo() > > 3. Upgrade to the latest R version maybe > > 4. Perhaps write to package maintainer (see ?maintainer) if nothing or > no one helps. > > Cheers, > Bert > > > Bert Gunter > > "The trouble with having an open mind is that people keep coming along > and sticking things into it." > -- Opus (aka Berkeley Breathed in his "Bloom County" comic strip ) > > >> On Tue, Nov 22, 2016 at 10:05 AM, Maram SAlem <marammagdysa...@gmail.com> >> wrote: >> Thanks for helping Jim. >> >> I'm actually using the pbapply function together with the print function >> within a loop. In earlier versions, the progress bar and the output of the >> print function used to appear after each iteration of the loop. But with the >> 3.3.1. Version nothing appears, instead the console turns white and the >> cursor turns blue ( busy) and I know nothing about the progress of the >> running code. >> >> I just want to see the bar and the output of the print function as I used >> to, any help? >> >> Thanks in advance. >> Maram Salem >> >> >> >> Sent from my iPhone >> >>> On Nov 3, 2016, at 8:30 PM, jim holtman <jholt...@gmail.com> wrote: >>> >>> A little more information would help. How exactly are out creating the >>> output to the console? Are you using 'print', 'cat' or something else? Do >>> you have buffered output checked on the GUI (you probably don't want it >>> checked or you output will be delayed till the buffer is full -- this might >>> be the cause of your problem. >>> >>> >>> 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 Thu, Nov 3, 2016 at 1:55 PM, Maram SAlem <marammagdysa...@gmail.com> >>>> wrote: >>>> Hi all, >>>> >>>> I've a question concerning the R 3.3.1 version. I have a long code that I >>>> used to run on versions earlier to the 3.3.1 version, and when I copied >>>> the code to the R console, I can still see the code while the loop is >>>> executing , along with the output printed after each iteration of the loop. >>>> >>>> Now, on the 3.3.1 version, after I copy the code to the console, it >>>> disappears and I only see the printed output of only one iteration at a >>>> time, that is, after the first iteration the printed output disappears ( >>>> though it's only 6 lines, just giving me some guidance, not a long output). >>>> This is causing me some problems, so I don't know if there is a general >>>> option for R that enables me to still see the code and the output of all >>>> the iterations till the loop is over, as was the case with earlier R >>>> versions. >>>> >>>> I didn't include the code as it's a long one. >>>> >>>> Thanks a lot in advance, >>>> >>>> Maram >>>> >>>> >>>> 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. >>> >> >>[[alternative HTML version deleted]] >> >> __ >> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see >> https://stat.ethz.ch/mailman/listinfo/r-help >> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html >> and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] The code itself disappears after starting to execute the for loop
Thanks for helping Jim. I'm actually using the pbapply function together with the print function within a loop. In earlier versions, the progress bar and the output of the print function used to appear after each iteration of the loop. But with the 3.3.1. Version nothing appears, instead the console turns white and the cursor turns blue ( busy) and I know nothing about the progress of the running code. I just want to see the bar and the output of the print function as I used to, any help? Thanks in advance. Maram Salem Sent from my iPhone > On Nov 3, 2016, at 8:30 PM, jim holtman <jholt...@gmail.com> wrote: > > A little more information would help. How exactly are out creating the > output to the console? Are you using 'print', 'cat' or something else? Do > you have buffered output checked on the GUI (you probably don't want it > checked or you output will be delayed till the buffer is full -- this might > be the cause of your problem. > > > 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 Thu, Nov 3, 2016 at 1:55 PM, Maram SAlem <marammagdysa...@gmail.com> >> wrote: >> Hi all, >> >> I've a question concerning the R 3.3.1 version. I have a long code that I >> used to run on versions earlier to the 3.3.1 version, and when I copied the >> code to the R console, I can still see the code while the loop is executing >> , along with the output printed after each iteration of the loop. >> >> Now, on the 3.3.1 version, after I copy the code to the console, it >> disappears and I only see the printed output of only one iteration at a >> time, that is, after the first iteration the printed output disappears ( >> though it's only 6 lines, just giving me some guidance, not a long output). >> This is causing me some problems, so I don't know if there is a general >> option for R that enables me to still see the code and the output of all the >> iterations till the loop is over, as was the case with earlier R versions. >> >> I didn't include the code as it's a long one. >> >> Thanks a lot in advance, >> >> Maram >> >> >> 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. > [[alternative HTML version deleted]] __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] The code itself disappears after starting to execute the for loop
Hi all, I've a question concerning the R 3.3.1 version. I have a long code that I used to run on versions earlier to the 3.3.1 version, and when I copied the code to the R console, I can still see the code while the loop is executing , along with the output printed after each iteration of the loop. Now, on the 3.3.1 version, after I copy the code to the console, it disappears and I only see the printed output of only one iteration at a time, that is, after the first iteration the printed output disappears ( though it's only 6 lines, just giving me some guidance, not a long output). This is causing me some problems, so I don't know if there is a general option for R that enables me to still see the code and the output of all the iterations till the loop is over, as was the case with earlier R versions. I didn't include the code as it's a long one. Thanks a lot in advance, Maram 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.
Re: [R] NaNs produced as a returned value for a function
Thanks for helping Dunlap.
I just don't get why does the function bettarg() executes well if I
evaluate it (by hand) step-by step . nut for the same input values , if I
tried to execute bettarg() it will give me the error
Error in if (y <= accept.prob) { : missing value where TRUE/FALSE needed
Although both y and accept.prob have values and are not missing.
Any help would be appreciated as this error is causing my entire(large and
complicated) code to stop execution.
Maram Salem
On 14 February 2016 at 22:08, William Dunlap <wdun...@tibco.com> wrote:
> You can do things like
>while ( !is.nan( r <- randomFunction(x) )) {}
># r will be a non-NaN value of randomFunction(x) now
> or
>for(i in seq_len(1000)) {
> if (!is.nan( r <- randomFunction(x))) {
> break
> }
> if (i == 1000) {
> stop("no good values of randomFunction(x) in 1000 tries")
> }
>}
> to avoid infinite loops when randomFunction [almost] always produces a NaN.
>
> You may be able to avoid some NaNs by replacing log(b^n) with n*log(b) and
> log(prod(x^b)) with sum(b*log(x)). That would avoid unneeded Inf values,
> which
> can lead to NaN down the line (Inf-Inf -> NaN, Inf/Inf -> NaN).
>
>
> Bill Dunlap
> TIBCO Software
> wdunlap tibco.com
>
> On Sun, Feb 14, 2016 at 9:22 AM, Maram SAlem <marammagdysa...@gmail.com>
> wrote:
>
>> Hi all,
>>
>> I'm trying to write 2 functions(as a part of a larger code) to evaluate a
>> certain equation. The function is :
>>
>> X= c (0.3893094 2.0962311 2.6007558 3.0761810 3.3246947 3.3917976
>> 4.1981546 6.8826140 12.3128013 15.5588470)
>> R=c (0 1 0 0 0 1 1 1 0 1)
>>
>> alpha.update=function(X, R, alpha.curr, beta.curr=1 ,m=10,
>> hyp=c(3,15,6,22.5))
>>
>> {
>>
>> o<-numeric(m)
>>
>> for (i in 1:m) {
>>
>>o[i]<- (1+R[i])*((X[i])^(beta.curr))
>>
>> }
>>
>>sh<-sum(o) + hyp[2] + (hyp[4]* beta.curr)
>>
>>rg<-rgamma(1, shape= m+hyp[1]+hyp[3] , rate = sh )
>>
>>return(rg)
>>
>>}
>>
>>
>> alpha.curr<- alpha.update(X, R, alpha.curr=0.2, beta.curr=1 ,m, hyp)
>>
>>
>> bettarg<- function(X, R, alpha.curr, beta.curr=1 ,m=10,
>> hyp=c(3,15,6,22.5))
>>
>>{
>>
>>o<-numeric(m)
>>
>>for (i in 1:m) {
>>
>>o[i]<- (1+R[i])*((X[i])^( beta.curr))
>>
>> }
>>
>> logbt<- log(beta.curr ^(m+hyp[3]-1)) + log(prod((X)^( beta.curr
>> -1)))
>> + (-1*alpha.curr *(sum(o) + (hyp[4]* beta.curr)))
>>
>>
>>
>> bt<- exp(logbt)
>>
>> return(bt)
>>
>>}
>>
>>
>> The problem is that the function bettarg() sometimes produces NaN, and
>> this stops the evaluation of my equation, so how can I force it to ignore
>> the NaNs produced and repeat the evaluation again till it prroduces a
>> number?
>>
>>
>> Thanks in advance,
>>
>>
>> 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.
[R] NaNs produced as a returned value for a function
Hi all,
I'm trying to write 2 functions(as a part of a larger code) to evaluate a
certain equation. The function is :
X= c (0.3893094 2.0962311 2.6007558 3.0761810 3.3246947 3.3917976
4.1981546 6.8826140 12.3128013 15.5588470)
R=c (0 1 0 0 0 1 1 1 0 1)
alpha.update=function(X, R, alpha.curr, beta.curr=1 ,m=10,
hyp=c(3,15,6,22.5))
{
o<-numeric(m)
for (i in 1:m) {
o[i]<- (1+R[i])*((X[i])^(beta.curr))
}
sh<-sum(o) + hyp[2] + (hyp[4]* beta.curr)
rg<-rgamma(1, shape= m+hyp[1]+hyp[3] , rate = sh )
return(rg)
}
alpha.curr<- alpha.update(X, R, alpha.curr=0.2, beta.curr=1 ,m, hyp)
bettarg<- function(X, R, alpha.curr, beta.curr=1 ,m=10,
hyp=c(3,15,6,22.5))
{
o<-numeric(m)
for (i in 1:m) {
o[i]<- (1+R[i])*((X[i])^( beta.curr))
}
logbt<- log(beta.curr ^(m+hyp[3]-1)) + log(prod((X)^( beta.curr -1)))
+ (-1*alpha.curr *(sum(o) + (hyp[4]* beta.curr)))
bt<- exp(logbt)
return(bt)
}
The problem is that the function bettarg() sometimes produces NaN, and
this stops the evaluation of my equation, so how can I force it to ignore
the NaNs produced and repeat the evaluation again till it prroduces a
number?
Thanks in advance,
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.
Re: [R] missing value where TRUE/FALSE needed ERROR
Thanks for helping thanoon.
I used the if() {
} else {
.
}
and it worked better
On 9 February 2016 at 15:56, thanoon younis <thanoon.youni...@gmail.com>
wrote:
> Hi,
>
> Try to check the function only
> if (runif(1) <= accept.prob) beta.cand else beta.curr
> something is missing here. Why you take only runif(1).
>
> Regards
>
> On 9 February 2016 at 16:01, Maram SAlem <marammagdysa...@gmail.com>
> wrote:
>
>> Hi all,
>>
>> I'm trying to write a function to implement a Metropolis-within-Gibbs
>> algorithm for two parameters.I'm including a naive version here so as to
>> be
>> able to spot the error I got. So I first generate the vectors, X and R,
>> that will help to start the algorithm using (for example):
>>
>> n=8; m=5; p=0.1; t=0.9 ; JH=10;
>>
>> R <- numeric(m)
>>
>> W <- numeric(m)
>>
>> V <- numeric(m)
>>
>> U <- numeric(m)
>>
>> X <- numeric(m)
>>
>> Bay.alpha<- numeric (JH)
>>
>> Bay.beta<- numeric (JH)
>>
>> Bay.Surv <- numeric (JH)
>>
>> hyp=c(3,15,6,22.5)
>>
>> theta<-c(0.2,2)
>>
>>
>> alpha.curr<-theta[1]
>>
>> beta.curr<- theta[2]
>>
>>
>>
>> R[1]<-rbinom(1, n-m, p)
>>
>> for (i in 2:m-1) {
>>
>> R[i]<-rbinom(1,n-m-sum(R[1:i-1]),p)
>>
>> }
>>
>> R[m]<-n-m-sum(R[1:m-1])
>>
>> W<-runif(m, min = 0, max = 1)
>>
>> for (i in 1:m){
>>
>> V[i]<-W[i]^(1/(i+sum(R[(m-i+1):m])))
>>
>> }
>>
>> for (i in 1:m){
>>
>> U[i]<- 1- prod(V[(m-i+1):m])
>>
>> }
>>
>> for (i in 1:m){
>>
>> X[i]<- ((-1/theta[1])*log(1-U[i]))^(1/theta[2])
>>
>> }
>>
>>
>> Then, I defined three functions 1- alpha.update() for updating alpha
>> (Gibbs
>> step) 2- bettarg(), for the target distribution of beta 3- beta.update()
>> for updating beta using a Metropolis Hastings technique.
>>
>>
>>alpha.update=function(X, R, alpha.curr, beta.curr ,m, hyp)
>>
>> {
>>
>> o<-numeric(m)
>>
>> for (i in 1:m) {
>>
>>o[i]<- (1+R[i])*((X[i])^(beta.curr))
>>
>> }
>>
>>sh<-sum(o) + hyp[2] + (hyp[4]* beta.curr)
>>
>>rg<-rgamma(1, shape= m+hyp[1]+hyp[3] , rate = sh )
>>
>>return(rg)
>>
>>}
>>
>>
>>
>>bettarg<- function(X, R, alpha.curr, beta.curr ,m, hyp)
>>
>>{
>>
>>o<-numeric(m)
>>
>>for (i in 1:m) {
>>
>>o[i]<- (1+R[i])*((X[i])^( beta.curr))
>>
>> }
>>
>> bt<- beta.curr ^(m+hyp[3]-1) * prod((X)^( beta.curr -1)) *
>> exp(-1*alpha.curr *(sum(o) + (hyp[4]* beta.curr)))
>>
>> return(bt)
>>
>>}
>>
>>beta.update<- function (X, R, alpha.curr, beta.curr ,m, hyp,
>> cand.sd)
>>
>>{
>>
>>beta.cand<- rnorm(1, mean = beta.curr, sd = cand.sd)
>>
>> AB<- bettarg (X, R, alpha.curr, beta.curr = beta.cand ,m, hyp)
>>
>> CD<- bettarg (X, R, alpha.curr, beta.curr ,m, hyp)
>>
>> accept.prob <- AB/CD
>>
>> if (runif(1) <= accept.prob)
>>
>> beta.cand
>>
>> else beta.curr
>>
>> }
>>
>>
>> Then I started a tiny chain using ten iterations but got this ERROR:
>>
>>
>> for (k in 1:JH)
>>
>> {
>>
>> alpha.curr<- alpha.update(X, R, alpha.curr, beta.curr ,m, hyp)
>>
>> Bay.alpha[k]<-alpha.curr
>>
>> beta.curr<- beta.update (X, R, alpha.curr, beta.curr ,m, hyp, cand.sd=2)
>>
>> Bay.beta[k]<- beta.curr
>>
>> Bay.Surv [k]<- exp (-1*Bay.alpha[k] * (t)^Bay.beta[k])
>>
>> }
>>
>>
>> Error in if (runif(1) <= accept.prob) beta.cand else beta.curr :
>>
>> missing value where TRUE/FALSE needed
>>
>> In addition: Warning message:
>>
>> In rgamma(1, shape = m + hyp[1] + hyp[3], rate = sh) : NAs produced
>>
>>
>> I ran the code several times for only one iteration but everything was
>> fine
>> with no errors or warnings, so I don't know from where does the missing
>> value/ NA come from?
>>
>> In addit
[R] missing value where TRUE/FALSE needed ERROR
Hi all,
I'm trying to write a function to implement a Metropolis-within-Gibbs
algorithm for two parameters.I'm including a naive version here so as to be
able to spot the error I got. So I first generate the vectors, X and R,
that will help to start the algorithm using (for example):
n=8; m=5; p=0.1; t=0.9 ; JH=10;
R <- numeric(m)
W <- numeric(m)
V <- numeric(m)
U <- numeric(m)
X <- numeric(m)
Bay.alpha<- numeric (JH)
Bay.beta<- numeric (JH)
Bay.Surv <- numeric (JH)
hyp=c(3,15,6,22.5)
theta<-c(0.2,2)
alpha.curr<-theta[1]
beta.curr<- theta[2]
R[1]<-rbinom(1, n-m, p)
for (i in 2:m-1) {
R[i]<-rbinom(1,n-m-sum(R[1:i-1]),p)
}
R[m]<-n-m-sum(R[1:m-1])
W<-runif(m, min = 0, max = 1)
for (i in 1:m){
V[i]<-W[i]^(1/(i+sum(R[(m-i+1):m])))
}
for (i in 1:m){
U[i]<- 1- prod(V[(m-i+1):m])
}
for (i in 1:m){
X[i]<- ((-1/theta[1])*log(1-U[i]))^(1/theta[2])
}
Then, I defined three functions 1- alpha.update() for updating alpha (Gibbs
step) 2- bettarg(), for the target distribution of beta 3- beta.update()
for updating beta using a Metropolis Hastings technique.
alpha.update=function(X, R, alpha.curr, beta.curr ,m, hyp)
{
o<-numeric(m)
for (i in 1:m) {
o[i]<- (1+R[i])*((X[i])^(beta.curr))
}
sh<-sum(o) + hyp[2] + (hyp[4]* beta.curr)
rg<-rgamma(1, shape= m+hyp[1]+hyp[3] , rate = sh )
return(rg)
}
bettarg<- function(X, R, alpha.curr, beta.curr ,m, hyp)
{
o<-numeric(m)
for (i in 1:m) {
o[i]<- (1+R[i])*((X[i])^( beta.curr))
}
bt<- beta.curr ^(m+hyp[3]-1) * prod((X)^( beta.curr -1)) *
exp(-1*alpha.curr *(sum(o) + (hyp[4]* beta.curr)))
return(bt)
}
beta.update<- function (X, R, alpha.curr, beta.curr ,m, hyp, cand.sd)
{
beta.cand<- rnorm(1, mean = beta.curr, sd = cand.sd)
AB<- bettarg (X, R, alpha.curr, beta.curr = beta.cand ,m, hyp)
CD<- bettarg (X, R, alpha.curr, beta.curr ,m, hyp)
accept.prob <- AB/CD
if (runif(1) <= accept.prob)
beta.cand
else beta.curr
}
Then I started a tiny chain using ten iterations but got this ERROR:
for (k in 1:JH)
{
alpha.curr<- alpha.update(X, R, alpha.curr, beta.curr ,m, hyp)
Bay.alpha[k]<-alpha.curr
beta.curr<- beta.update (X, R, alpha.curr, beta.curr ,m, hyp, cand.sd=2)
Bay.beta[k]<- beta.curr
Bay.Surv [k]<- exp (-1*Bay.alpha[k] * (t)^Bay.beta[k])
}
Error in if (runif(1) <= accept.prob) beta.cand else beta.curr :
missing value where TRUE/FALSE needed
In addition: Warning message:
In rgamma(1, shape = m + hyp[1] + hyp[3], rate = sh) : NAs produced
I ran the code several times for only one iteration but everything was fine
with no errors or warnings, so I don't know from where does the missing
value/ NA come from?
In addition, I want to calculate the acceptance rate for the
Metropolis-Hastings step, is this possible?
Any help would be 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.
[R] Metropolis-within-Gibbs algorithm
Hi all, I'm trying to write a code that performs the Metropolis-within-Gibbs algorithm, to draw values of a 2x1 parameter vector from a posterior distribution that doesn't have a well known form. So one of the parameters, theta1, has a well known full conditional distribution( for which the gibbs sampler can be used), but the other, theta2, doesn't have a well known full conditional (for which a random walk Metropolis-Hastings algorithm should be used). But theta1 depends on the generated value of theta 2. I know this code can be hand-written, but is there any package that can perform such update of the Metropolis-within-gibbs algorithm and provide me with acceptance rates or different scaling for the proposal distribution? I've checked the gibbs.met package but it used an independent proposal distribution. I went through the MCMCpack but couldn't find a function that performs what I want. So is there any other packages that you recommend or do I have to write my own functions? Thanks for helping, 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.
Re: [R] Alternatives for explicit for() loops
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 <jholt...@gmail.com> 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 <marammagdysa...@gmail.com>
> 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 thought 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(0,x[i],1)
>>}
>>LD[[m]]<-seq(0,(n-m-sum(x)),1)
>>LED<-expand.grid (LD)
>>LED<-as.matrix(LED)
>>store1<-numeric(nrow(LED))
>> for (j in 1:length(store1) )
>> {
>> incomb<-function(x,alpha,beta) {
>>
>> g<-((-1)^(sum(LED[j,])))*(gamma((1/beta)+1))*((alpha)^(-(1/beta)))
>> h <- choose(x, LED[j,-m])
>>ik<-prod(h)*choose((n-m-sum(x)),LED[j,m])
>> lm<-cumsum(LED[j,-m]) + (1:(m-1))
>> plm<-prod(lm)
>>gil<-g*ik/(plm)
>> hlm<-numeric(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<-pbapply(s,1,simpfun,n=6,m
Re: [R] Alternatives for explicit for() loops
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 thought 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(0,x[i],1)
}
LD[[m]]<-seq(0,(n-m-sum(x)),1)
LED<-expand.grid (LD)
LED<-as.matrix(LED)
store1<-numeric(nrow(LED))
for (j in 1:length(store1) )
{
incomb<-function(x,alpha,beta) {
g<-((-1)^(sum(LED[j,])))*(gamma((1/beta)+1))*((alpha)^(-(1/beta)))
h <- choose(x, LED[j,-m])
ik<-prod(h)*choose((n-m-sum(x)),LED[j,m])
lm<-cumsum(LED[j,-m]) + (1:(m-1))
plm<-prod(lm)
gil<-g*ik/(plm)
hlm<-numeric(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<-pbapply(s,1,simpfun,n=6,m=4,p=0.3,alpha=0.2,beta=2)
EXP<-sum(va)
Any help would be greatly appreciated.
Thanks a lot for your time.
Best Regards,
Maram Salem
On 2 November 2015 at 00:27, jim holtman <jholt...@gmail.com> wrote:
> Why are you recreating the incomb function within the loop instead of
> defining it outside the loop? Also you are referencing several variables
> that are global (e.g., m & j); you should be passing these in as parameters
> to the function.
>
>
> 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 1, 2015 at 7:31 AM, Maram SAlem <marammagdysa...@gmail.com>
> wrote:
>
>> Hi All,
>>
>> I'm writing a long code that takes long time to execute. So I used the
>> Rprof() function and found out that the function that takes about 80% of
>> the time is the incomb () fucntion (below), and this is most probably
>> because of the many explicit for() loops I'm using.
>>
>> n=18;m=4;p=0.3;alpha=0.2;beta=2
>> x=c(3,0,0)
>> 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<-numeric(nrow(LED))
>> h<- numeric(m-1)
>> lm<- numeric(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)
>
[R] Alternatives for explicit for() loops
Hi All,
I'm writing a long code that takes long time to execute. So I used the
Rprof() function and found out that the function that takes about 80% of
the time is the incomb () fucntion (below), and this is most probably
because of the many explicit for() loops I'm using.
n=18;m=4;p=0.3;alpha=0.2;beta=2
x=c(3,0,0)
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<-numeric(nrow(LED))
h<- numeric(m-1)
lm<- numeric(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<-numeric(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)
}
I'm trying to use alternatives (for ex. to vectorize things) to the
explicit for() loops, but things don't work out.
Any suggestions that can help me to speed up the execution of the incomb()
function are much appreciated.
Thanks a lot in advance.
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.
Re: [R] Time it takes to run a code
Thanks for helping Boris. Regards, Maram Salem On 25 October 2015 at 23:30, Boris Steipe <boris.ste...@utoronto.ca> wrote: > It may be useful for you to estimate the time complexity of your function: > try it with smaller input that takes short and noticeable time, see whether > the time increases linearly, quadratically, or exponentially with the > number of elements you process, then extrapolate to your full data set. > > > To see what your function is doing, you could > - add a print statement that tells you the progress every 1000 or 10,000 > items or so... > - add a progress bar > > https://stat.ethz.ch/R-manual/R-patched/library/utils/html/txtProgressBar.html > - install and use the pbapply package > https://cran.r-project.org/web/packages/pbapply/pbapply.pdf > > > You are right to make sure you are getting some feedback, it is easy to > make a mistake in function logic that will cause a function to fail to > terminate. > > > B. > > > > On Oct 25, 2015, at 5:05 PM, Maram SAlem <marammagdysa...@gmail.com> > wrote: > > > Hi All, > > > > I'm using a function, say func, and I want to apply it to all the rows > of a > > certain matrix. The problem is that my code kept on running for more than > > two days without giving any output. I've made some modifications.But is > > there a way to know the time needed to execute my code and reach an ouput > > but before running the code and not after it is run? I've tried > Sys.time() > > and system.time(), but still the console freezed for so long and I didn't > > reach anything. > > > > Any suggestions are much appreciated. > > Thanks in advance. > > > > 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.
[R] Avoiding for loops
Hi All,
I wonder if I can avoid the for() loop in any of the following loops.These
loops are a part of a larger code which I'm trying to accelerate.
n=6
m=4
x<-c(0,1,1)
1st loop
for (i in 1:m-1)
{
d[i]<- n- (sum(x[(1):(i)])) - i
}
e<- n*(prod(d))
2nd loop
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)
3rd loop
for (i in 1:(m-1))
{
h[i]<- choose(x[i],LED[j,i])
}
4th loop
for (i in 1:(m-1))
{
lm[i]<-(sum(LED[j,1:i])) + i
}
I appreciate if anyone has any suggestions or references.
Thanks in advance.
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.
Re: [R] Avoiding for loops
Thanks a lot Boris and Berend.
I'll consider the brackets ((m-1) in every loop). In addition, I'll read
more on profiling my code. In fact,I'm using the apply () in another part
of my code.
Thanks again for helping.
Maram Salem
On 25 October 2015 at 14:26, Berend Hasselman <b...@xs4all.nl> wrote:
>
> > On 25 Oct 2015, at 11:42, Maram SAlem <marammagdysa...@gmail.com> wrote:
> >
> > Hi All,
> >
> > I wonder if I can avoid the for() loop in any of the following
> loops.These
> > loops are a part of a larger code which I'm trying to accelerate.
> >
> > n=6
> > m=4
> > x<-c(0,1,1)
> >
> > 1st loop
> >
> > for (i in 1:m-1)
> > {
> > d[i]<- n- (sum(x[(1):(i)])) - i
> > }
> > e<- n*(prod(d))
> >
> >
> On the basis of the other loops I presume you mean
>
> for( in in 1:(m-1))
>
> in stead of what you wrote.
>
> You get the same result with
>
> d <- n - cumsum(x) - (1:(m-1))
>
>
> Berend
>
>
> > 2nd loop
> >
> > 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)
> >
> >
> > 3rd loop
> >
> > for (i in 1:(m-1))
> >
> > {
> >
> >h[i]<- choose(x[i],LED[j,i])
> >
> > }
> >
> >
> >
> > 4th loop
> >
> >
> > for (i in 1:(m-1))
> >
> > {
> >
> >lm[i]<-(sum(LED[j,1:i])) + i
> >
> > }
> >
> >
> > I appreciate if anyone has any suggestions or references.
> >
> >
> > Thanks in advance.
> >
> >
> > 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.
[R] Time it takes to run a code
Hi All, I'm using a function, say func, and I want to apply it to all the rows of a certain matrix. The problem is that my code kept on running for more than two days without giving any output. I've made some modifications.But is there a way to know the time needed to execute my code and reach an ouput but before running the code and not after it is run? I've tried Sys.time() and system.time(), but still the console freezed for so long and I didn't reach anything. Any suggestions are much appreciated. Thanks in advance. 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.
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
>>
>> [[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] 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))
> >
Re: [R] Error in rep.int() invalid 'times' value
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,]00
>
> [2,]10
>
> [3,]20
>
> [4,]01
>
> [5,]11
>
> [6,]02
>
> > f1(5,3)
>
> [,1] [,2]
>
> [1,]00
>
> [2,]10
>
> [3,]20
>
> [4,]01
>
> [5,]11
>
> [6,]02
>
> >
>
>
>
> 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
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 <wdun...@tibco.com> 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 <marammagdysa...@gmail.com>
> 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 <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
> >>>
> >>> [[alternative HTML version deleted]]
> >>>
> >>> __
> >>> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> >>> https://stat.ethz.ch/mailman/listinfo/r-help
> >>> PLEASE do read the posting guide
> >>> http://www.R-project.org/posting-guide.html
> >>> and provide commented, minimal, self-contained, reproducible code.
> >>
> >>
> >
>
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] Error in rep.int() invalid 'times' value
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.
Re: [R] using the apply() family to evaluate nested functions with common arguments
I'm sorry, I forgot to mention that in the last step I have to use the
choose() function for all the elements of x not only x[1],x[2], and x[3].
Any suggestions or recommendation for some reference or a package that can
help me sort that out.
Thanks.
Maram Salem
On 11 October 2015 at 18:52, Maram SAlem <marammagdysa...@gmail.com> wrote:
> Dear All,
>
> I'm trying to use the apply family to evaluate 2 nested functions whose
> arguments are somewhat overlapping. I'm doing this in order to evaluate
> nested multiple summations in a certain equation. First, I created the s
> matrix whose rows satisfy some logical condition.
>
> n=8
>
> m=4
>
> 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,]
>
>
> Then, I'm trying to evaluate a function called simpfun for each row of s,
> which could be easily done using the apply () function. The problem is that
> within the simpfun I need to evaluate another function, incomb(). This
> function has to be first evaluated for each row of the matrix LED, whose
> elements are sequences having the corresponding elements of each row of s
> as their upper limits.
>
>
> 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)
>
> incomb<-function(x,alpha,beta) {
>
>
> g<-((-1)^(sum(LED[1,])))*(gamma((1/beta)+1))*((alpha)^(-(1/beta)))
>
> fd<-
> choose(x[1],LED[1,1])*choose(x[2],LED[1,2])*choose(x[3],LED[1,3])
>
> return (g)
>
> }
>
>
>
> }
>
>
>
> where my x in the simpfun() refers to one of the rows of the s matrix, so
> that I'll be able to use something like
>
>
> va<-apply(s,1,simpfun,n=,m=,p=,alpha=,beta=)
>
>
>
> to apply it to each row of s.
>
>
>
> The problem now is that for each row of s (which is supposed to be my x in
> the simpfun()) ,I need to first apply the incomb() function for all the
> rows of LED.Thus, I need to modify what I wrote above in the body of the
> incomb() function in terms of something like y instead of LED[1,], so that
> I can again use the apply() function on all its rows first for one row of s
> ,say x, and then repeat this for all the other rows of s. I can't figure
> out how to do this while still having the rows of s, say x, as one of the
> arguments of the incomb() function for which I'm going to use the apply()
> function once more.
>
>
> I'm sorry if what I'm asking for is not that clear, but I'm trying to
> simplfy things as much as possible so that we don't go into tedious detalis.
>
>
> Thanks a lot in advance.
>
>
> 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.
[R] error with seq() used in for loop
Dear All,
I have a question concerning the seq() function arguments when used in a
for() loop.
I'll simplify this question to make it more clear. Suppose I have a
(5*3)matrix s (for ex.) and I need to write a function with for() loop, in
each step of the loop I need to generate a sequence whose upper limit is
the ith element of the first row of s, then put the resulting sequences in
a list. I used the following simple code (I've only included the first part
of the function)
> s<-matrix(c(1,0,1,0,1,1,0,0,1,0,0,0,0,0,1),nrow=5,byrow=TRUE)
> simpfun<- function (x,n,m,p,alpha,beta)
+ {
+ LD<-list()
+ for (i in 1:m-1)
+ {
+ LD[[i]]<-seq(0,x[i],1)
+ }
+ print(LD)
+ }
> mk<-simpfun(s[1,],n=6,m=4,p=0.3)
Error in seq.default(0, x[i], 1) : 'to' must be of length 1
Although x is supposed to be the vector
1 0 1
and thus x[1]=1, x[2]=0,x[3]=1.
So I don't get why the error "Error in seq.default(0, x[i], 1) : 'to' must
be of length 1" occurs in the first place.
Thanks for helping.
Maram
[[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] error with seq() used in for loop
Thanks Berend
On 11 October 2015 at 14:10, Berend Hasselman <b...@xs4all.nl> wrote:
>
> > On 11 Oct 2015, at 13:52, Maram SAlem <marammagdysa...@gmail.com> wrote:
> >
> > Dear All,
> > I have a question concerning the seq() function arguments when used in a
> > for() loop.
> > I'll simplify this question to make it more clear. Suppose I have a
> > (5*3)matrix s (for ex.) and I need to write a function with for() loop,
> in
> > each step of the loop I need to generate a sequence whose upper limit is
> > the ith element of the first row of s, then put the resulting sequences
> in
> > a list. I used the following simple code (I've only included the first
> part
> > of the function)
> >
> >
> >> s<-matrix(c(1,0,1,0,1,1,0,0,1,0,0,0,0,0,1),nrow=5,byrow=TRUE)
> >> simpfun<- function (x,n,m,p,alpha,beta)
> > + {
> > + LD<-list()
> > + for (i in 1:m-1)
> > + {
> > + LD[[i]]<-seq(0,x[i],1)
> > + }
> > + print(LD)
> > + }
> >> mk<-simpfun(s[1,],n=6,m=4,p=0.3)
> > Error in seq.default(0, x[i], 1) : 'to' must be of length 1
> >
> > Although x is supposed to be the vector
> > 1 0 1
> > and thus x[1]=1, x[2]=0,x[3]=1.
> > So I don't get why the error "Error in seq.default(0, x[i], 1) : 'to'
> must
> > be of length 1" occurs in the first place.
>
>
> Use
>
> for (i in 1:(m-1))
>
> in stead of what you put in the function.
> Read the section “Generating regular sequences” in the "An Introduction to
> R“ manual.
>
> Berend
>
>
>
[[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] error with seq() used in for loop
Thanks a lot for your help and patience Duncan.
It seems that my questions is really a trivial one but I never realized
that 1:4 means 0 1 2 3 4 , never knew it starts from 0.
On 11 October 2015 at 14:07, Duncan Murdoch <murdoch.dun...@gmail.com>
wrote:
> On 11/10/2015 7:52 AM, Maram SAlem wrote:
> > Dear All,
> > I have a question concerning the seq() function arguments when used in a
> > for() loop.
> > I'll simplify this question to make it more clear. Suppose I have a
> > (5*3)matrix s (for ex.) and I need to write a function with for() loop,
> in
> > each step of the loop I need to generate a sequence whose upper limit is
> > the ith element of the first row of s, then put the resulting sequences
> in
> > a list. I used the following simple code (I've only included the first
> part
> > of the function)
> >
> >
> >> s<-matrix(c(1,0,1,0,1,1,0,0,1,0,0,0,0,0,1),nrow=5,byrow=TRUE)
> >> simpfun<- function (x,n,m,p,alpha,beta)
> > + {
> > + LD<-list()
> > + for (i in 1:m-1)
> > + {
> > + LD[[i]]<-seq(0,x[i],1)
> > + }
> > + print(LD)
> > + }
> >> mk<-simpfun(s[1,],n=6,m=4,p=0.3)
> > Error in seq.default(0, x[i], 1) : 'to' must be of length 1
> >
> > Although x is supposed to be the vector
> > 1 0 1
> > and thus x[1]=1, x[2]=0,x[3]=1.
> > So I don't get why the error "Error in seq.default(0, x[i], 1) : 'to'
> must
> > be of length 1" occurs in the first place.
>
> The range of your loop is 1:m-1, where m is 4. That is
>
> > m <- 4
> > 1:m-1
> [1] 0 1 2 3
>
> and x[0] is length 0.
>
> I think you wanted 1:(m-1) (or even better, seq_len(m-1)) for your loop
> values.
>
> Duncan Murdoch
>
> >
> > Thanks for helping.
> >
> > Maram
> >
> > [[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] error with seq() used in for loop
I got you Berend.
Thanks again
On 11 October 2015 at 14:20, Berend Hasselman <b...@xs4all.nl> wrote:
>
> > On 11 Oct 2015, at 14:12, Maram SAlem <marammagdysa...@gmail.com> wrote:
> >
> > Thanks a lot for your help and patience Duncan.
> > It seems that my questions is really a trivial one but I never realized
> > that 1:4 means 0 1 2 3 4 , never knew it starts from 0.
> >
>
> It doesn’t mean what you mention.
>
> 1:4 is the sequence 1 2 3 4
>
> 1:4-1 is (1 2 3 4) -1 ==> 0 1 2 3
>
> The 1 is subtracted from all elements of the sequence.
>
> Berend
> >
> >
> > On 11 October 2015 at 14:07, Duncan Murdoch <murdoch.dun...@gmail.com>
> > wrote:
> >
> >> On 11/10/2015 7:52 AM, Maram SAlem wrote:
> >>> Dear All,
> >>> I have a question concerning the seq() function arguments when used in
> a
> >>> for() loop.
> >>> I'll simplify this question to make it more clear. Suppose I have a
> >>> (5*3)matrix s (for ex.) and I need to write a function with for()
> loop,
> >> in
> >>> each step of the loop I need to generate a sequence whose upper limit
> is
> >>> the ith element of the first row of s, then put the resulting sequences
> >> in
> >>> a list. I used the following simple code (I've only included the first
> >> part
> >>> of the function)
> >>>
> >>>
> >>>> s<-matrix(c(1,0,1,0,1,1,0,0,1,0,0,0,0,0,1),nrow=5,byrow=TRUE)
> >>>> simpfun<- function (x,n,m,p,alpha,beta)
> >>> + {
> >>> + LD<-list()
> >>> + for (i in 1:m-1)
> >>> + {
> >>> + LD[[i]]<-seq(0,x[i],1)
> >>> + }
> >>> + print(LD)
> >>> + }
> >>>> mk<-simpfun(s[1,],n=6,m=4,p=0.3)
> >>> Error in seq.default(0, x[i], 1) : 'to' must be of length 1
> >>>
> >>> Although x is supposed to be the vector
> >>> 1 0 1
> >>> and thus x[1]=1, x[2]=0,x[3]=1.
> >>> So I don't get why the error "Error in seq.default(0, x[i], 1) : 'to'
> >> must
> >>> be of length 1" occurs in the first place.
> >>
> >> The range of your loop is 1:m-1, where m is 4. That is
> >>
> >>> m <- 4
> >>> 1:m-1
> >> [1] 0 1 2 3
> >>
> >> and x[0] is length 0.
> >>
> >> I think you wanted 1:(m-1) (or even better, seq_len(m-1)) for your loop
> >> values.
> >>
> >> Duncan Murdoch
> >>
> >>>
> >>> Thanks for helping.
> >>>
> >>> Maram
> >>>
> >>> [[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.
>
>
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] using the apply() family to evaluate nested functions with common arguments
Dear All,
I'm trying to use the apply family to evaluate 2 nested functions whose
arguments are somewhat overlapping. I'm doing this in order to evaluate
nested multiple summations in a certain equation. First, I created the s
matrix whose rows satisfy some logical condition.
n=8
m=4
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,]
Then, I'm trying to evaluate a function called simpfun for each row of s,
which could be easily done using the apply () function. The problem is that
within the simpfun I need to evaluate another function, incomb(). This
function has to be first evaluated for each row of the matrix LED, whose
elements are sequences having the corresponding elements of each row of s
as their upper limits.
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)
incomb<-function(x,alpha,beta) {
g<-((-1)^(sum(LED[1,])))*(gamma((1/beta)+1))*((alpha)^(-(1/beta)))
fd<-
choose(x[1],LED[1,1])*choose(x[2],LED[1,2])*choose(x[3],LED[1,3])
return (g)
}
}
where my x in the simpfun() refers to one of the rows of the s matrix, so
that I'll be able to use something like
va<-apply(s,1,simpfun,n=,m=,p=,alpha=,beta=)
to apply it to each row of s.
The problem now is that for each row of s (which is supposed to be my x in
the simpfun()) ,I need to first apply the incomb() function for all the
rows of LED.Thus, I need to modify what I wrote above in the body of the
incomb() function in terms of something like y instead of LED[1,], so that
I can again use the apply() function on all its rows first for one row of s
,say x, and then repeat this for all the other rows of s. I can't figure
out how to do this while still having the rows of s, say x, as one of the
arguments of the incomb() function for which I'm going to use the apply()
function once more.
I'm sorry if what I'm asking for is not that clear, but I'm trying to
simplfy things as much as possible so that we don't go into tedious detalis.
Thanks a lot in advance.
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.
Re: [R] (subscript) logical subscript too long
Thanks Giorgio, I got it.
I managed to reach the matrix s whose rows represent all the possible
combinations. Here is the code:
> n=12
> m=7
> 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,]
The problem now is that the code works only for relatively small values of
n and m, but when I use, for ex., n=20 and m=9, I got this error
> n=20
> m=9
> 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))
*Error: cannot allocate vector of size 1.6 Gb*
*Any Suggestions please?*
*Thanks Again.*
*Maram*
On 30 September 2015 at 17:41, Giorgio Garziano <
giorgio.garzi...@ericsson.com> wrote:
> Be:
>
> log <- (rowSums(ED) <= (n - m))
>
>
> Compare the following two values:
>
>
>
>
>
> length(log)
>
>
>
>
>
> nrow(w)
>
>
>
>
>
> --
>
> GG
>
>
> [[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] (subscript) logical subscript too long
Thanks a lot Giorgio for your Help
Regards,
Maram
On 1 October 2015 at 14:27, Giorgio Garziano <giorgio.garzi...@ericsson.com>
wrote:
> If you are running a 32-bit Windows, there are following upper limits:
>
>
>
>
> https://cran.r-project.org/bin/windows/base/rw-FAQ.html#There-seems-to-be-a-limit-on-the-memory-it-uses_0021
>
>
>
>
>
> starts by:
>
>
>
> memory.limit(size=1920)
>
>
>
> and try increasing value of size as a parameter for memory.limit().
>
>
>
>
>
> I use Intel i5 Windows-7 64-bit 16GB RAM.
>
>
>
>
>
> GG
>
>
>
>
>
>
>
> *From:* Maram SAlem [mailto:marammagdysa...@gmail.com]
> *Sent:* giovedì 1 ottobre 2015 14:12
>
> *To:* Giorgio Garziano
> *Cc:* r-help@r-project.org
> *Subject:* Re: [R] (subscript) logical subscript too long
>
>
>
> Thanks a lot Giorgio, I used
>
>
>
> memory.limit(size=4096)
>
>
>
> but got
>
>
>
> don't be silly!: your machine has a 4Gb address limit
>
>
>
> I'm working on my Ph.D. thesis and I have a huge code of which this is
> just a very small part, so does this error mean that I need a new computer
> with extended capabilites to be able to execute my code?? I'm currently
> using intel core i3, windows 7
>
>
>
> Thanks for helping.
>
>
>
> Maram
>
>
>
>
>
> On 1 October 2015 at 13:37, Giorgio Garziano <
> giorgio.garzi...@ericsson.com> wrote:
>
> Check your memory size by:
>
>
>
> memory.limit()
>
>
>
> try to increase it by:
>
>
>
> memory.limit(size=4096)
>
>
>
>
>
>
>
> *From:* Maram SAlem [mailto:marammagdysa...@gmail.com]
> *Sent:* giovedì 1 ottobre 2015 13:22
> *To:* Giorgio Garziano
> *Cc:* r-help@r-project.org
> *Subject:* Re: [R] (subscript) logical subscript too long
>
>
>
> Thanks Giorgio, I got it.
>
>
>
> I managed to reach the matrix s whose rows represent all the possible
> combinations. Here is the code:
>
>
>
> > n=12
>
> > m=7
>
> > 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,]
>
>
>
> The problem now is that the code works only for relatively small values of
> n and m, but when I use, for ex., n=20 and m=9, I got this error
>
>
>
> > n=20
>
> > m=9
>
> > 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))
>
>
>
> *Error: cannot allocate vector of size 1.6 Gb*
>
>
>
> *Any Suggestions please?*
>
>
>
> *Thanks Again.*
>
> *Maram*
>
>
>
> On 30 September 2015 at 17:41, Giorgio Garziano <
> giorgio.garzi...@ericsson.com> wrote:
>
> Be:
>
> log <- (rowSums(ED) <= (n - m))
>
>
> Compare the following two values:
>
>
>
>
>
> length(log)
>
>
>
>
>
> nrow(w)
>
>
>
>
>
> --
>
> GG
>
>
> [[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] (subscript) logical subscript too long
Thanks a lot Giorgio, I used
memory.limit(size=4096)
but got
don't be silly!: your machine has a 4Gb address limit
I'm working on my Ph.D. thesis and I have a huge code of which this is just
a very small part, so does this error mean that I need a new computer with
extended capabilites to be able to execute my code?? I'm currently using
intel core i3, windows 7
Thanks for helping.
Maram
On 1 October 2015 at 13:37, Giorgio Garziano <giorgio.garzi...@ericsson.com>
wrote:
> Check your memory size by:
>
>
>
> memory.limit()
>
>
>
> try to increase it by:
>
>
>
> memory.limit(size=4096)
>
>
>
>
>
>
>
> *From:* Maram SAlem [mailto:marammagdysa...@gmail.com]
> *Sent:* giovedì 1 ottobre 2015 13:22
> *To:* Giorgio Garziano
> *Cc:* r-help@r-project.org
> *Subject:* Re: [R] (subscript) logical subscript too long
>
>
>
> Thanks Giorgio, I got it.
>
>
>
> I managed to reach the matrix s whose rows represent all the possible
> combinations. Here is the code:
>
>
>
> > n=12
>
> > m=7
>
> > 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,]
>
>
>
> The problem now is that the code works only for relatively small values of
> n and m, but when I use, for ex., n=20 and m=9, I got this error
>
>
>
> > n=20
>
> > m=9
>
> > 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))
>
>
>
> *Error: cannot allocate vector of size 1.6 Gb*
>
>
>
> *Any Suggestions please?*
>
>
>
> *Thanks Again.*
>
> *Maram*
>
>
>
> On 30 September 2015 at 17:41, Giorgio Garziano <
> giorgio.garzi...@ericsson.com> wrote:
>
> Be:
>
> log <- (rowSums(ED) <= (n - m))
>
>
> Compare the following two values:
>
>
>
>
>
> length(log)
>
>
>
>
>
> nrow(w)
>
>
>
>
>
> --
>
> GG
>
>
> [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
>
>
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] (subscript) logical subscript too long
Dear All,
I'm trying to write a function in the values of some numeric vectors
(d1,d2,...,d(m-1)). This function should be applied on some combinations
of the elements of the (m-1) d vectors that satisfy the condition of having
a sum less than or equal to (n-m). I've tried the following code, but got
an error:(subscript) logical subscript too long.
> n=20
> m=6
> 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)
> s<-w[rowSums(ED)<=(n-m),]
Error in w[rowSums(ED) <= (n - m), ] :
(subscript) logical subscript too long
Bearing in mind that the true values of n and m might be way larger than 20
and 6, consecutively.
Any Suggestions please?
Thanks.
Maram
[[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] writing an equation with multiple summation
My first problem is that the counters involved in the summation are
dependent, For instance for each r[i], 0<r[i]< m-(r[1] +r[2]+...+r[i-1])
so I tried using seq() and for loops to express all the possible
combinations but I couldn't.
On 28 September 2015 at 17:26, Maram SAlem <marammagdysa...@gmail.com>
wrote:
> Dear All,
>
> I'm trying to write and evaluate an equation which involves multiple
> summations but can't figure out how to do it.
>
> I've an numeric vector r
> r<-vector(mode = "numeric", length = m)
> and I have multiple summations (for ex.) of the form:
>
> [(sum from r[1]=0 to g(r[1])) (sum from r[2] =0 to g(r[2]))..(sum from
> r[m] to g(r[m]))] {the sum is over some complicated expression in
> r[1],r[2],.,r[m]},
>
> where g(r[i]) = m- (r[1] +r[2]+...+r[i-1])
> Any suggestions for some function or a package that can help me with this?
>
> Many Thanks,
> Maram
>
>
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] writing an equation with multiple summation
Dear All,
I'm trying to write and evaluate an equation which involves multiple
summations but can't figure out how to do it.
I've an numeric vector r
r<-vector(mode = "numeric", length = m)
and I have multiple summations (for ex.) of the form:
[(sum from r[1]=0 to g(r[1])) (sum from r[2] =0 to g(r[2]))..(sum from
r[m] to g(r[m]))] {the sum is over some complicated expression in
r[1],r[2],.,r[m]},
where g(r[i]) = m- (r[1] +r[2]+...+r[i-1])
Any suggestions for some function or a package that can help me with this?
Many Thanks,
Maram
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] Generating Progressively censored samples
Dear All, Is there some built-in function in R that can be used to generate progressively censored sample from a certain distribution, for example, the Weibull distribution? OR Do I have to write the code of the algorithm myself? Thanks for helping. 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.
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 arne.henning...@gmail.com 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 marammagdysa...@gmail.com:
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])/((theta[1]+theta[2])*T+(theta[1]*theta[2])))^2)+(
theta[3])*sum(((C+theta[1])/((theta[1]+theta[2])*C+(theta[1]*theta[2])))^2)
G[2,3]-(n/theta[2])+(-1)*sum((T+theta[1])/((theta[1]+theta[2])*T+(theta[1]*theta[2])))+(-1)*sum((C+theta[1])/((theta[1]+theta[2])*C+(theta[1]*theta[2])))
G[3,2]-G[2,3]
G[3,3]-((-1*r)/(theta[3])^2)
return(G)
}
and using numeric gradient and hessian matrix:
a - maxLik(loglik, start=c(40,50,2))
Warning messages:
1: In log(theta[3]) : NaNs produced
2: In log(theta[1] + theta[2
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 arne.henning...@gmail.com 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 marammagdysa...@gmail.com:
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])/((theta[1]+theta[2])*T+(theta[1]*theta[2])))^2)+(
theta[3])*sum(((C+theta[1])/((theta[1]+theta[2])*C+(theta[1]*theta[2])))^2)
G[2,3]-(n/theta[2])+(-1)*sum((T+theta[1])/((theta[1]+theta[2])*T+(theta[1]*theta[2])))+(-1)*sum((C+theta[1])/((theta[1]+theta[2])*C+(theta[1]*theta[2])))
G[3,2]-G[2,3]
G[3,3]-((-1*r)/(theta[3])^2)
return(G)
}
and using numeric gradient and hessian matrix:
a - maxLik(loglik, start=c(40,50,2))
Warning messages:
1: In log(theta[3]) : NaNs produced
2: In log(theta[1] + theta[2
[R] Warning message with maxLik()
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])/((theta[1]+theta[2])*T+(theta[1]*theta[2])))^2)+(
theta[3])*sum(((C+theta[1])/((theta[1]+theta[2])*C+(theta[1]*theta[2])))^2)
G[2,3]-(n/theta[2])+(-1)*sum((T+theta[1])/((theta[1]+theta[2])*T+(theta[1]*theta[2])))+(-1)*sum((C+theta[1])/((theta[1]+theta[2])*C+(theta[1]*theta[2])))
G[3,2]-G[2,3]
G[3,3]-((-1*r)/(theta[3])^2)
return(G)
}
and using numeric gradient and hessian matrix:
a - maxLik(loglik, start=c(40,50,2))
Warning messages:
1: In log(theta[3]) : NaNs produced
2: In log(theta[1] + theta[2]) : NaNs produced
3: In log(theta[1]) : NaNs produced
4: In log((T * (theta[1] + theta[2])) + (theta[1] * theta[2])) : NaNs
produced
5: In log((C * (theta[1] + theta[2])) + (theta[1] * theta[2])) : NaNs
produced
6: In log(theta[3]) : NaNs produced
7: In log(theta[1] + theta[2]) : NaNs produced
and so on…..
I don't know why I get these 50 warnings although:
1- The inputs of the log() function are strictly positive.
2- When I evaluated the log-likelihood fuction at the very begining it gave
me a number(which is -56.66) and not (NAN).
I've also tried to:
1- Reparamtrize my model using lamda(i)= log(theta(i)), for i=1,2,3, so
that it may solve the problem, but it didn't.
2- I've used the comparederivitive() function, and the analytic and numeric
gradients were so close.
Any help please?
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.
Re: [R] NaN produced from log() with positive input
Dear Arne, Will do. Thanks for helping. Maram Sent from my iPhone On Jul 12, 2015, at 8:23 AM, Arne Henningsen arne.henning...@gmail.com wrote: Dear Maram On 8 July 2015 at 17:52, Maram Salem marammagdysa...@gmx.com wrote: Dear Arne, On a second thought, as per your mail the warning messages occur each time, when maxLik() tries to calculate the logLik value for theta[1] = 0, theta[1] + theta[2] = 0, theta[3] = 0 or something similar. The component of the theta vector are all indeed strictly positive, and the initial values I used are c(40,50,2). and this means that theta[1] 0, theta[1] + theta[2] 0, and theta[3] 0 . These initial values are the parameter values that were used for generating the data (the C and T vectors). That's why I don't know why the warnings occur in the first place and why the estimates are far away from the initial values.a Any suggestions? - don't send your messages twice; - do what I suggested in my previous e-mail; - increase the number of observations; - check your log-likelihood function; - use an optimisation algorithm that does not use derivatives; - use numeric derivatives in the optimisation; - check your function for returning the derivatives of the log-likelihood function, e.g. with compareDerivatives(); - check the data generating process Best regards, Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Can't post to the list
Dear Rolf, I recieved this message exactly: The message's content type was not explicitly allowed Thanks for your concern, Maram Sent from my iPhone On Jul 9, 2015, at 11:36 AM, Rolf Turner r.tur...@auckland.ac.nz wrote: On 09/07/15 19:49, Maram Salem wrote: Sorry I mean NOT explicitly allowed You probably mean explicitly NOT allowed. The order of the words makes a substantial difference to the meaning. Sent from my iPhone On Jul 9, 2015, at 8:46 AM, Maram Salem marammagdysa...@gmx.com wrote: Dear all, I need help with posting to the list. Whenever i post any message, i get this reply : the message content was explicitly allowed. Any help please? Thanks a lot. You have just successfully posted two messages to the list. So think about the error message that you got. It told you that the content of your (rejected) message was disallowed. What was it about that content that would make the list disallow it? Then remove the offending content. cheers, Rolf Turner -- Technical Editor ANZJS Department of Statistics University of Auckland Phone: +64-9-373-7599 ext. 88276 __ 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] Can't post to the list
Sorry I mean NOT explicitly allowed Sent from my iPhone On Jul 9, 2015, at 8:46 AM, Maram Salem marammagdysa...@gmx.com wrote: Dear all, I need help with posting to the list. Whenever i post any message, i get this reply : the message content was explicitly allowed. Any help please? Thanks a lot. Maram 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] Can't post to the list
Dear all, I need help with posting to the list. Whenever i post any message, i get this reply : the message content was explicitly allowed. Any help please? Thanks a lot. Maram 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.
Re: [R] Can't post to the list
Thanks a lot John. Will check it out. Maram Sent from my iPhone On Jul 9, 2015, at 1:17 PM, John McKown john.archie.mck...@gmail.com wrote: On Thu, Jul 9, 2015 at 6:04 AM, Maram Salem marammagdysa...@gmx.com wrote: Dear Rolf, I recieved this message exactly: The message's content type was not explicitly allowed I'm just guessing here! But this list is plain text only. I wonder if your email client sometimes uses unsupported MIME types. Again, I'm no expert, but this might be if it included screen pictures in JPG format (Windows loves these from the Print Screen button). Or maybe BASE64 or UUENCODEd somehow. If the listserv software detect this, and postings with that type of data are not in the setup as explicitly allowed, then maybe. Here's a similar discussion which may be of interest: https://mailman-1.sys.kth.se/pipermail/gromacs.org_gmx-developers/2012-August/006176.html quote ... If you have received rejected mails saying The message's content type was not explicitly allowed it means that you have used in your mail special styles or formatting (e.g. bold typeface, colors etc.). ... /quote Thanks for your concern, Maram Sent from my iPhone -- Schrodinger's backup: The condition of any backup is unknown until a restore is attempted. Yoda of Borg, we are. Futile, resistance is, yes. Assimilated, you will be. He's about as useful as a wax frying pan. 10 to the 12th power microphones = 1 Megaphone Maranatha! John McKown [[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] Can't post to the list
Dear Rolf, There isn't any offending content. I think it's related to the styles used on the mail. Thanks for helping. Maram Sent from my iPhone On Jul 9, 2015, at 11:36 AM, Rolf Turner r.tur...@auckland.ac.nz wrote: On 09/07/15 19:49, Maram Salem wrote: Sorry I mean NOT explicitly allowed You probably mean explicitly NOT allowed. The order of the words makes a substantial difference to the meaning. Sent from my iPhone On Jul 9, 2015, at 8:46 AM, Maram Salem marammagdysa...@gmx.com wrote: Dear all, I need help with posting to the list. Whenever i post any message, i get this reply : the message content was explicitly allowed. Any help please? Thanks a lot. You have just successfully posted two messages to the list. So think about the error message that you got. It told you that the content of your (rejected) message was disallowed. What was it about that content that would make the list disallow it? Then remove the offending content. cheers, Rolf Turner -- Technical Editor ANZJS Department of Statistics University of Auckland Phone: +64-9-373-7599 ext. 88276 __ 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] NaN produced from log() with positive input
Dear Arne,
Sorry for posting my mail twice and thanks a lot for your help.
Best regards,
Maram
Sent from my iPhone
On Jul 7, 2015, at 9:55 PM, Arne Henningsen arne.henning...@gmail.com wrote:
Dear Maram
Please do NOT post your message twice!
The warning messages occur each time, when maxLik() tries to calculate
the logLik value for theta[1] = 0, theta[1] + theta[2] = 0, theta[3]
= 0 or something similar. According to the log-likelihood function,
it seems that the parameters theta[1], theta[2], and theta[3] must be
strictly positive. I suggest to re-parameterise your model so that the
estimated parameters can take any values between minus infinity and
infinity, e.g. by theta[1] - exp( param[1] ); theta[2] - exp(
param[2] ); theta[3] - exp( param[3] ) so that your estimated
parameter vector 'param' consists of log( theta[1] ), log( theta[2] ),
and log( theta[3] ). After the estimation, you can obtain the
estimated values of the thetas by exp( param[1] ), exp( param[2] ),
and exp( param[3] ) .
Best regards,
Arne
2015-07-06 2:29 GMT+02:00 Maram Salem marammagdysa...@gmx.com:
Dear All
I'm trying to find the maximum likelihood estimator of a certain
distribution based on the newton raphson method using maxLik package. I
wrote the log-likelihood , gradient, and hessian functionsusing the
following code.
#Step 1: Creating the theta vector
theta -vector(mode = numeric, length = 3)
# Step 2: Setting the values of r and n
r- 17
n -30
# Step 3: Creating the T vector
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)
# Step 4: Creating the C vector
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)
}
# Step 5: Creating the gradient vector and calculating its inputs
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)
}
# Step 6: Creating the G (Hessian) matrix and Calculating its inputs
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])/((theta[1]+theta[2])*T+(theta[1]*theta[2])))^2)+(
theta[3])*sum(((C+theta[1])/((theta[1]+theta[2])*C+(theta[1]*theta[2])))^2)
G[2,3]-(n/theta[2])+(-1)*sum((T+theta[1])/((theta[1]+theta[2])*T+(theta[1]*theta[2])))+(-1)*sum((C+theta[1])/((theta[1]+theta[2])*C+(theta[1]*theta[2])))
G[3,2]-G[2,3]
G[3,3]-((-1*r)/(theta[3])^2)
return(G)
}
mle-maxLik(loglik, grad = gradlik, hess = hesslik, start=c(40,50,2))
There were 50 or more warnings (use warnings() to see the first 50)
warnings ()
Warning messages:
1: In log(theta[3]) : NaNs produced
2: In log(theta[1] + theta[2]) : NaNs produced
3: In log(theta[1]) : NaNs produced
4: In log((T * (theta[1
[R] Can't get a reasonable the maximum likelihood estimator
Sent from my iPhone
Begin forwarded message:
From: Maram Salem marammagdysa...@gmx.com
Date: July 6, 2015 at 2:29:56 AM GMT+2
To: r-help@r-project.org r-help@r-project.org
Subject: [R] NaN produced from log() with positive input
Dear All
I'm trying to find the maximum likelihood estimator of a certain
distribution based on the newton raphson method using maxLik package. I wrote
the log-likelihood , gradient, and hessian functions using the following code.
#Step 1: Creating the theta vector
theta -vector(mode = numeric, length = 3)
# Step 2: Setting the values of r and n
r- 17
n -30
# Step 3: Creating the T vector
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)
# Step 4: Creating the C vector
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)
}
# Step 5: Creating the gradient vector and calculating its inputs
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)
}
# Step 6: Creating the G (Hessian) matrix and Calculating its inputs
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])/((theta[1]+theta[2])*T+(theta[1]*theta[2])))^2)+(
theta[3])*sum(((C+theta[1])/((theta[1]+theta[2])*C+(theta[1]*theta[2])))^2)
G[2,3]-(n/theta[2])+(-1)*sum((T+theta[1])/((theta[1]+theta[2])*T+(theta[1]*theta[2])))+(-1)*sum((C+theta[1])/((theta[1]+theta[2])*C+(theta[1]*theta[2])))
G[3,2]-G[2,3]
G[3,3]-((-1*r)/(theta[3])^2)
return(G)
}
mle-maxLik(loglik, grad = gradlik, hess = hesslik, start=c(40,50,2))
There were 50 or more warnings (use warnings() to see the first 50)
warnings ()
Warning messages:
1: In log(theta[3]) : NaNs produced
2: In log(theta[1] + theta[2]) : NaNs produced
3: In log(theta[1]) : NaNs produced
4: In log((T * (theta[1] + theta[2])) + (theta[1] * theta[2])) : NaNs produced
and so on ...
Although when I evaluate, for example, log(theta[3]) it gives me a number.
and the same applies for the other warnings.
Then when I used summary (mle), I got
Maximum Likelihood estimation
Newton-Raphson maximisation, 7 iterations
Return code 1: gradient close to zero
Log-Likelihood: -55.89012
3 free parameters
Estimates:
Estimate Std. error t value Pr( t)
[1,] 11.132Inf 0 1
[2,] 47.618Inf 0 1
[3,]1.293Inf 0 1
Where the estimates are far away from the starting values and they have
infinite standard errors. I think there is a problem with my gradlik or
hesslik functions, but I can't figure it out.
Any help?
Thank you in advance.
Maram
[R] NaN produced from log() with positive input
Dear All
I'm trying to find the maximum likelihood estimator of a certain distribution
based on the newton raphson method using maxLik package. I wrote the
log-likelihood , gradient, and hessian functionsusing the following code.
#Step 1: Creating the theta vector
theta -vector(mode = numeric, length = 3)
# Step 2: Setting the values of r and n
r- 17
n -30
# Step 3: Creating the T vector
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)
# Step 4: Creating the C vector
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)
}
# Step 5: Creating the gradient vector and calculating its inputs
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)
}
# Step 6: Creating the G (Hessian) matrix and Calculating its inputs
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])/((theta[1]+theta[2])*T+(theta[1]*theta[2])))^2)+(
theta[3])*sum(((C+theta[1])/((theta[1]+theta[2])*C+(theta[1]*theta[2])))^2)
G[2,3]-(n/theta[2])+(-1)*sum((T+theta[1])/((theta[1]+theta[2])*T+(theta[1]*theta[2])))+(-1)*sum((C+theta[1])/((theta[1]+theta[2])*C+(theta[1]*theta[2])))
G[3,2]-G[2,3]
G[3,3]-((-1*r)/(theta[3])^2)
return(G)
}
mle-maxLik(loglik, grad = gradlik, hess = hesslik, start=c(40,50,2))
There were 50 or more warnings (use warnings() to see the first 50)
warnings ()
Warning messages:
1: In log(theta[3]) : NaNs produced
2: In log(theta[1] + theta[2]) : NaNs produced
3: In log(theta[1]) : NaNs produced
4: In log((T * (theta[1] + theta[2])) + (theta[1] * theta[2])) : NaNs produced
and so on ...
Although when I evaluate, for example, log(theta[3]) it gives me a number. and
the same applies for the other warnings.
Then when I used summary (mle), I got
Maximum Likelihood estimation
Newton-Raphson maximisation, 7 iterations
Return code 1: gradient close to zero
Log-Likelihood: -55.89012
3 free parameters
Estimates:
Estimate Std. error t value Pr( t)
[1,] 11.132Inf 0 1
[2,] 47.618Inf 0 1
[3,]1.293Inf 0 1
Where the estimates are far away from the starting values and they have
infinite standard errors. I think there is a problem with my gradlik or hesslik
functions, but I can't figure it out.
Any help?
Thank you in advance.
Maram
[[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,
[R] histogram scott
Dear all, I want to use the histogtam as a density estimator, with the binwidths calculated using scott's formula which is binwidth = 3.49*ST.dev.*n^(-1/3) for the following data (30 data points) 12-9-3-6-1-23-21-7-18-16-15-4-19-22-20-2-3-18-8-10-1-7-5-4-11-12-3-9-19-7 so first,I' ve tried this manually, and substituted in the above formula and I got st.dev.=7.02745 and thus the binwidth=7.89313 But when I used hist with breaks = scott, that is h-hist(x,breaks=scott) I got the breaks in the histogram object = 0 10 20 30 that is, the binwidth used is equal to 10 not 7.89313?? I don't know why? shouldn't they be exactly the same?? Thanks a lot, Maram. [[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] solving a linear equation
Hi all, I've a linear equation of the form: 0.95=2 ( [3+ln(x/3)]^-13 + 4 [3+2ln(x/3)]^-13 + [3+4ln(x/3)]^-13 ) and I want to solve it for x, can I do this using R? Thanks in advance. Maram. [[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] Fw: solving a linear equation
Sorry, I didn't mean a linear equation, of course this equation is not linear, I meant an equation in one unknown. Hi all, I've a linear equation of the form: 0.95=2 ( [3+ln(x/3)]^-13 + 4 [3+2ln(x/3)]^-13 + [3+4ln(x/3)]^-13 ) and I want to solve it for x, can I do this using R? Thanks in advance. Maram. [[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] axis label
Hi all, I want the y-axis label to be ( in symbols) g(sigma given alpha) where given is the conditional sign. I've tried ylab=expression(g(sigma|alpha))) but it gave me g(|(sigma,alpha)) where the sigma and alpha are in greek but the conditional sign is misplaced (before the bracket) Any help would be appreciated. Maram. [[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] cdf
Dear all, I have the cdf of the following power fuction distribution: F(y)=(y/350)^a ,0y350, where a is some parameter with range a0. I want to use it as the argument of the discretize function of the actuar package. So I think I need to define this function to R so that if I entered a=1, I get the following F(y)=(y/350) and if I entered a=4.5, I get the following F(y) =(y/350)^4.5 ... and so on I've tried a-vector(mode=numeric,length=1) powercdf-function(a,y) (y/350)^a But when I typed: powercdf(10,y) instead of getting : (y/350)^10 (which is what I want) I got : object y not found ?? I want y to remain as it is, a continous variable, not for example seq(0,350). Thank you in advance. Maram [[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] finding the minimum value
Hi all, I'm using a certain procedure to calculate the value of some variable(Bayes risk),B. So I got the values B1, B2, , B1000, each under certain input values and using a long procedure. Now, I want to put the values I got in a nummerical vector and find their minimum value. I think c( ) should work.For example if I have only 10 values I could have used c(B1,B2,B3,B4,B5,B6,B7,B8,B9,B10) But how can I do this for 1000 values? I think the question is really trivial, but I tried to work on it and couldn't reach anythg. Thanks. Maram [[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] using histogram to find cdf
Dear all, How can I use the histogram density estimate (hist) to find the value of the cdf at a certain point? Thanks Maram [[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] converting a matrix into a fn
Dear All, I have a matrix m of the form m [,1] [,2] [1,] 9072.302 0.0004462366 [2,] 9086.811 0.0005628169 [3,] 9101.320 0.0007126347 [4,] 9115.830 0.0008986976 [5,] 9130.339 0.001126 [6,] 9144.848 0.0014018405 [7,] 9159.357 0.0017344229 [8,] 9173.866 0.0021363563 [9,] 9188.375 0.0026389996 [10,] 9202.884 0.0032406678 where the first col is my variable x and the second is f(x). I've obtained this from the result of a kernel density estimate so I don't know the explict form of f(x) . For some reason i want to integerate f(x) over a certain range using integrate( ). So i want R to consider the above matrix as a function where the first col is the argument x and the second is the realization f(x). I've tried as.function and vectorize but it didn't work out. Could u help me do it? Thanks Maram [[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] calculating probability
Hi all, I've a trivial question. If (q) is a continous variable,actually a vector of 1000 values. how to calculate the probability that q is greater than a specific value, i.e. P(q45)?? Thanks Maram [[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] kernel density estimates
Dear All, I have a variable q which is a vector of 1000 simulated positive values; that is I generated 1000 samples from the pareto distribution, from each sample I calculated the value of q ( a certain fn in the sample observations), and thus I was left with 1000 values of q and I don't know the distribution of q. Hence, I used the given code for kernel density estimation to estimate the density of q options(scipen=4) d - density(q, bw = nrd0,kernel=gaussian) d plot(d) But what I'm really intersed in is to estimate the probability that q is greater than a certain value , for ex.,P(q11000), using the kernel density estimate I obtained. Could u help me with a fn or some document to do this? Thank u so much Maram [[alternative HTML version deleted]] [[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] Fw: Hist kernel density estimates
Dear All, Here are the codes of a histogram a kernel density estimates that I used. For the hist estimate par(mex=1.3) dens-density(q) options(scipen=4) ylim-range(dens$y) h-hist(q,breaks=scott,freq=FALSE,probability=TRUE, + right=FALSE,xlim=c(9000,16000),ylim=ylim,main=Histogram of q(scott)) lines(dens) box() For the kernel estimateoptions(scipen=4) d - density(q, bw = nrd0,kernel=gaussian) d plot(d) In fact the variable q is a vector of 1000 simulated positive values; that is I generated 1000 samples from the pareto distribution, from each sample I calculated the value of q ( a certain fn in the sample observations), and thus I was left with 1000 values of q and I don't know the distribution of q. Hence, I used the above codes for histogram and kernel density estimation to estimate the density of q. But what I'm really intersed in is to estimate the probability that q is greater than a certain value , for ex.,P(q11000), using the density estimates I obtained. Could u help me with a fn or some document to do this? Thank u so much Maram [[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] Hist kernel density estimates
Dear All, Attached are the codes of a histogram a kernel density estimate and the output they produced. In fact the variable q is a vector of 1000 simulated values; that is I generated 1000 samples from the pareto distribution, from each sample I calculated the value of q ( a certain fn in the sample observations), and thus I was left with 1000 values of q and I don't know the distribution of q. Hence, I used the attached codes for histogram and kernel density estimation toestimate the density of q. But what I'm really intersed in is to estimate the probability that q is greater than a certain value , for ex.,P(q11000), using the density estimates I obtained. Could u help me with a fn or some document to do this? Thank u so much Maram __ 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] Fw: Hist kernel density estimates
For the hist estimate par(mex=1.3) dens-density(q) options(scipen=4) ylim-range(dens$y) h-hist(q,breaks=scott,freq=FALSE,probability=TRUE, + right=FALSE,xlim=c(9000,16000),ylim=ylim,main=Histogram of q(scott)) lines(dens) box() For the kernel estimateoptions(scipen=4) d - density(q, bw = nrd0,kernel=gaussian) d plot(d) In fact the variable q is a vector of 1000 simulated values; that is I generated 1000 samples from the pareto distribution, from each sample I calculated the value of q ( a certain fn in the sample observations), and thus I was left with 1000 values of q and I don't know the distribution of q. Hence, I used the attached codes for histogram and kernel density estimation toestimate the density of q. But what I'm really intersed in is to estimate the probability that q is greater than a certain value , for ex.,P(q11000), using the density estimates I obtained. Could u help me with a fn or some document to do this? Thank u so much Maram Dear All, Attached are the codes of a histogram a kernel density estimate and the output they produced. I'll copy the codes here in case there's something wrong with the attachement __ 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] axis scale
Dear All, I'm trying to plot a histogram (with the relative frequencies as the Y axis), But the scale of the y axis is given by 0e+00, 1e-04, 2e-04, 3e-04,. Now, I have 2 questions 1- Does (1e-04=0.01831563)? 2- If this true,how can i change the given scale to (0.01,0.03,0.05,0.07,0.09)? [[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] Fw: (no subject)
Dear group, Thank u so much 4 ur help. I've tried the link, http://finzi.psych.upenn.edu/R/library/quantreg/html/akj.html for adaptive kernel density estimation. But since I'm an R beginer and the topic of adaptive estimation is new for me, i still can't figure out some of the arguments of akj(x, z =, p =, h = -1, alpha = 0.5, kappa = 0.9, iker1 = 0) I've a vector of 1000 values (my X), but I don't know how to get the Z and what's Kappa? I'm sorry if the question is trivial but I hope u could recommend some refrence if u know one. Thank u so much again Maram From: John Kane jrkrid...@yahoo.ca Sent: Monday, June 29, 2009 10:35:49 PM Subject: Re: [R] (no subject) Perhaps? http://finzi.psych.upenn.edu/R/library/quantreg/html/akj.html Subject: [R] (no subject) To: r-help@r-project.org Received: Monday, June 29, 2009, 9:05 AM Hi group, I found a module for adaptive kernel density estimation for Stata users, but unfortunetly I don't have access to Stata, can I find a similar approach using R? [[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. __ Looking for the perfect gift? Give the gift of Flickr! [[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] (no subject)
Hi Group, I've a vector of 1000 numeric values for which I want to draw a histogram. I've read this vector into R with no variable name.I mean only the 1000 values, which makes V1 the name of the variable by default?? Then I tried hist(V1, breaks = Sturges, + freq = NULL, probability = !freq, + include.lowest = TRUE, right = TRUE, + density = NULL, angle = 45, col = NULL, border = NULL, + main = paste(Histogram of , V1name), + V1lim = range(breaks), ylim = NULL, + V1lab = V1name, ylab, + axes = TRUE, plot = TRUE, labels = FALSE, + nclass = NULL) It gave me this error: Error in hist(V1, breaks = Sturges, freq = NULL, probability = !freq, : object V1 not found I don't get what's wrong,( An R beginner)?? [[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] (no subject)
Hi group, I found a module for adaptive kernel density estimation for Stata users, but unfortunetly I don't have access to Stata, can I find a similar approach using R? [[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
Hi group, I found a module for adaptive kernel density estimation for Stata users, but unfortunetly I don't have access to Stata, can I find a similar approach using R? Thank u so much 4 ur time. [[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.
