[R] Smart way to sum a column in a matrix across all members in a list
Hi, I have a list that contains about 100 elements. Each element contains the following matrix structure a b 1 1 4 2 2 5 3 3 6 I'd like to sum column b across all elements. I'm thinking of using lapply or sapply but cannot get it to work. Any help would be really appreciated. Thank you. adschai __ 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] How to read and write an S4 class
Hi, I have an object in S4 class. Is there anyway that I can write to a file and read it back to R? Thank you. adschai __ 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] aggregating irregular time series
Hi, I have a couple of aggregation operations that I don't know how to accomplish. Let me give an example. I have the following irregular time series time x 10:00:00.02120 10:00:00.22420 10:00:01.00219 10:00:02:94820 1) For each entry time, I'd like to get sum of x for the next 2 seconds (excluding itself). Using the above example, the output should be timesumx 10:00:00.02139 10:00:00.22419 10:00:01.44220 10:00:02:948 0 2) For each i-th of x in the series, what's the first passage time to x[i]-1. I.e. the output should be time firstPassgeTime 10:00:00.0210.981 10:00:00.2240.778 10:00:01.442NA 10:00:02:948NA Is there any shortcut function that allows me to do the above? Thank you. adschai __ 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] POSIX time conversion doesn't display digit
Hi, I have the following string that I converted to a POSIXct: > a <- "2009-08-24 10:00:00.213" > a.p <- strptime(a,"%Y-%m-%d %H:%M:%OS") > as.double(as.POSIXct(a.p)) [1] 1251122400 I can't seem to get the decimal fraction of .213 out by casting to double or numeric. Is there anyway to make sure that POSIXct spits that out? Thank you. adschai __ 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] Aggregating irregular time series
Hi, I have a couple of aggregation operations that I don't know how to accomplish. Let me give an example. I have the following irregular time series time x 10:00:00.021 20 10:00:00.224 20 10:00:01.002 19 10:00:02:948 20 1) For each entry time, I'd like to get sum of x for the next 2 seconds (excluding itself). Using the above example, the output should be time sumx 10:00:00.021 39 10:00:00.224 19 10:00:01.442 20 10:00:02:948 0 2) For each i-th of x in the series, what's the first passage time to x[i]-1. I.e. the output should be time firstPassgeTime 10:00:00.021 0.981 10:00:00.224 0.778 10:00:01.442 NA 10:00:02:948 NA Is there any shortcut function that allows me to do the above? Thank you. adschai __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] [R-SIG-Finance] Aggregating irregular time series
Hi Gabor,
I appreciate your suggestion. I'm wondering if there's any faster
implementation that one could achieve. The dataset I have is about
400K rows and I have many of them. Just wondering if you have any
thoughts. Thanks.
Sakda
On Tue, Sep 8, 2009 at 6:57 AM, R_help Help wrote:
> -- Forwarded message --
> From: Gabor Grothendieck
> Date: Sun, Aug 30, 2009 at 11:08 PM
> Subject: Re: [R-SIG-Finance] Aggregating irregular time series
> To: R_help Help
> Cc: r-sig-fina...@stat.math.ethz.ch, r-help@r-project.org
>
>
> Try this for the first question:
>
> neighborapply <- function(z, width, FUN) {
> out <- z
> ix <- seq_along(z)
> jx <- findInterval(time(z) + width, time(z))
> out[] <- mapply(function(i, j) FUN(c(0, z[seq(i+1, length =
> j-i)])), ix, jx)
> out
> }
>
> # test - corrected :948 in last line
>
> library(zoo)
> library(chron)
>
> Lines <- "Time x
> 10:00:00.021 20
> 10:00:00.224 20
> 10:00:01.002 19
> 10:00:02.948 20"
>
> z <- read.zoo(textConnection(Lines), header = TRUE, FUN = times)
>
> neighborapply(z, times("00:00:02"), sum)
>
> # and here is an alternative neighborapply
> # using loops. The non-loop solution does
> # have the disadvantage that it does as
> # readily extend to other situations which
> # is why we add this second solution to
> # the first question.
>
> neighborapply <- function(z, width, FUN) {
> out <- z
> tt <- time(z)
> for(i in seq_along(tt)) {
> for(j in seq_along(tt)) {
> if (tt[j] - tt[i] > width) break
> }
> if (j == length(tt) && tt[j] - tt[i] <= width) j <- j+1
> out[i] <- FUN(c(0, z[seq(i+1, length = j-i-1)]))
> }
> out
> }
>
> The second question can be answered along the lines
> of the first by modifying neighborapply in the loop solution
> appropriately.
>
> On Sun, Aug 30, 2009 at 9:38 PM, R_help Help wrote:
>> Hi,
>>
>> I have a couple of aggregation operations that I don't know how to
>> accomplish. Let me give an example. I have the following irregular
>> time series
>>
>> time x
>> 10:00:00.021 20
>> 10:00:00.224 20
>> 10:00:01.002 19
>> 10:00:02:948 20
>>
>> 1) For each entry time, I'd like to get sum of x for the next 2
>> seconds (excluding itself). Using the above example, the output should
>> be
>>
>> time sumx
>> 10:00:00.021 39
>> 10:00:00.224 19
>> 10:00:01.442 20
>> 10:00:02:948 0
>>
>> 2) For each i-th of x in the series, what's the first passage time to
>> x[i]-1. I.e. the output should be
>>
>> time firstPassgeTime
>> 10:00:00.021 0.981
>> 10:00:00.224 0.778
>> 10:00:01.442 NA
>> 10:00:02:948 NA
>>
>> Is there any shortcut function that allows me to do the above? Thank you.
>>
>> adschai
>>
>> ___
>> r-sig-fina...@stat.math.ethz.ch mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-finance
>> -- Subscriber-posting only.
>> -- If you want to post, subscribe first.
>>
>
__
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] Creating Rcpp RcppDatetime from string with millisecond
Hi, I was trying to generate RcppDatetime from a string. The main problem for me is that my string contains millisecond. I saw that RcppDatetime takes in double of seconds since epoch. I try to generate a double using boost but has no success. I'm wondering if you have any sample snippet? The string looks like: "2010-03-18 15:50:51.232" Secondly, what time zone RcppDatetime is based on? Is it local machine time zone? And if I have a need to convert time zone, is there any tool at hand I can use? Or I have to resort to something like boost again? Thank you so much in advance. Rob __ 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] apply on xts
Hi,
I do not understand why after I called apply on a function that returns an
xts (getIdvAdjSeries) it returns a matrix whose columns are just numeric
value of time series in xts instead of a list of xts objects.
Basically, I called the following:
apply(matrix(tickers,ncol=1),1,FUN=getDivAdjSeries)
getDivAdjSeries <- function(ticker) {
seriesName <- paste(ticker,"Adjusted",sep=".");
command <- parse(text=paste(ticker,"[,'",seriesName,"']",sep=""));
s <- eval(command);
dimnames(s)[[2]] <- ticker;
command <- parse(text=paste(ticker,"@index",sep=""));
s <- xts(s,index=eval(command));
return(s);
}
This doesn't seem to work. Can anyone shed some light please? Thank you.
adschai
[[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] Can hist plot % of count instead of raw count?
I do not want to use histogram because I have more control over base graphs that R provide. I'm wondering if hist function provides a parameter that can set y axis as % of count instead of raw count? If not, I'd like to know how I can add text or line on histogram function? I don't like the idea of having to use raw coordinates to set text position. Thank you. adschai [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] apply on xts
Thanks David for reading for me. Anyhow, is there better way to do looping
on xts objects than using apply? Thanks again.
adschai
On Thu, Jun 25, 2009 at 12:39 AM, David Winsemius wrote:
>
> On Jun 24, 2009, at 9:26 PM, R_help Help wrote:
>
> Hi,
>>
>> I do not understand why after I called apply on a function that returns an
>> xts (getIdvAdjSeries) it returns a matrix whose columns are just numeric
>> value of time series in xts instead of a list of xts objects.
>>
>> Basically, I called the following:
>>
>> apply(matrix(tickers,ncol=1),1,FUN=getDivAdjSeries)
>>
>> getDivAdjSeries <- function(ticker) {
>> seriesName <- paste(ticker,"Adjusted",sep=".");
>> command <- parse(text=paste(ticker,"[,'",seriesName,"']",sep=""));
>> s <- eval(command);
>> dimnames(s)[[2]] <- ticker;
>> command <- parse(text=paste(ticker,"@index",sep=""));
>> s <- xts(s,index=eval(command));
>> return(s);
>> }
>>
>> This doesn't seem to work. Can anyone shed some light please? Thank you.
>>
>
> That would appear to be the expected behavior after reading the Value
> section of the help page for apply.
>
>
> David Winsemius, MD
> Heritage Laboratories
> West Hartford, CT
>
>
[[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] Defining class as a member of another class
Hi,
When I define a new class (through setClass), members defined in
representation argument doesn't seem to like a class. For example, if
I do the following:
setClass("NotWork",representation=(x="zoo"))
It seems to me that representation members will take in only primitive
type to R. Is there any way to stuff a class as a member in another
class? Thank you.
- adschai
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Defining class as a member of another class
Thank you Gabor. It works. I read a bit more and understood what you're doing.
One more question, I want to know more about .Data property of a
class. I know that it defines type of the class when asking typeof.
But I don't know much how to use them. Would you mind pointing me to
some reference and example? Thank you again!
ads
On Sat, Jun 27, 2009 at 10:29 PM, Gabor
Grothendieck wrote:
> Try this:
>
>> setClass("zoo")
> [1] "zoo"
>> setClass("Work",representation=(x="zoo"))
> [1] "Work"
>
>
> On Sat, Jun 27, 2009 at 10:01 PM, R_help Help wrote:
>> Hi,
>>
>> When I define a new class (through setClass), members defined in
>> representation argument doesn't seem to like a class. For example, if
>> I do the following:
>>
>> setClass("NotWork",representation=(x="zoo"))
>>
>> It seems to me that representation members will take in only primitive
>> type to R. Is there any way to stuff a class as a member in another
>> class? Thank you.
>>
>> - adschai
>>
>> __
>> 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] Adding RcppFrame to RcppResultSet causes segmentation fault
Hi,
I'm a bit puzzled. I uses exactly the same code in RcppExamples
package to try adding RcppFrame object to RcppResultSet. When running
it gives me segmentation fault problem. I'm using gcc 4.1.2 on redhat
64bit. I'm not sure if this is the cause of the problem. Any advice
would be greatly appreciated. Thank you.
Rob.
int numCol=4;
std::vector colNames(numCol);
colNames[0] = "alpha"; // column of strings
colNames[1] = "beta"; // column of reals
colNames[2] = "gamma"; // factor column
colNames[3] = "delta"; // column of Dates
RcppFrame frame(colNames);
// Third column will be a factor. In the current implementation the
// level names are copied to every factor value (and factors
// in the same column must have the same level names). The level names
// for a particular column will be factored out (pardon the pun) in
// a future release.
int numLevels = 2;
std::string *levelNames = new std::string[2];
levelNames[0] = std::string("pass"); // level 1
levelNames[1] = std::string("fail"); // level 2
// First row (this one determines column types).
std::vector row1(numCol);
row1[0].setStringValue("a");
row1[1].setDoubleValue(3.14);
row1[2].setFactorValue(levelNames, numLevels, 1);
row1[3].setDateValue(RcppDate(7,4,2006));
frame.addRow(row1);
// Second row.
std::vector row2(numCol);
row2[0].setStringValue("b");
row2[1].setDoubleValue(6.28);
row2[2].setFactorValue(levelNames, numLevels, 1);
row2[3].setDateValue(RcppDate(12,25,2006));
frame.addRow(row2);
RcppResultSet rs;
rs.add("PreDF", frame);
__
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] how to fit time varying coefficient regression model?
Hi - I read through dse package manual a bit. I'm not quite certain how I can use it to estimate a time varying coefficient regression model? I might pick up an inappropriate package. Any suggestion would be greatly appreciated. Thank you. rh __ 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] Evaluating/comparing dynamic linear model
Hi, I have two DLM model specifications (x[t] and y[t] are univariate): MODEL1: y[t] = b[t]x[t]+e[t], e[t] ~ N(0,v1^2) b[t] = b[t-1]+eta[t], eta[t] ~ N(0,w1^2) MODEL2: y[t] = a[t]+e[t], e[t] ~ N(0,v2^2) a[t] = a[t-1]+eta[t], eta[t] ~ N(0,w2^2) I run the filter through data recursively to obtain state variables for each model. However, how do I know if b[t]x[t] in MODEL1 is different from MODEL2? In other words, how do I know if x[t] makes a difference in explaining dynamic of y[t]? Another question is that how do I compare MODEL1 and MODEL2? From model specification point of view, how can one say that MODEL1 is better than MODEL2? Any suggestion/reference would be greatly appreciated. Thank you. ac __ 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] Plot log scale
Hi - a simple question. I know that on plot one can set log="xy" to set scale on both axes to log scale. I use this to plot autocorrelation of a slow decaying autocorrelated process. It has some points that are zero or slightly negative and cause this plot to return error. I'm wondering if there is anyway to force omitting these points? Thank you. rc __ 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] Implementation of gamma function for large number
Hi - I ran into a problem when the argument to gamma function is large. Says, if I do the following: > gamma(17000) [1] Inf Warning message: value out of range in 'gammafn' Is there anyway to get around this or other implementation? Thank you. -rc __ 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] Poisson dpois value is too small for double precision thus corrupts loglikelihood
Hi - I have a likelihood function that involves sums of two possions: L = a*dpois(Xi,theta1)*dpois(Yi,theta2)+b*(1-c)*a*dpois(Xi,theta1+theta3)*dpois(Yi,theta2) where a,b,c,theta1,theta2,theta3 are parameters to be estimated. (Xi,Yi) are observations. However, Xi and Yi are usually big (> 2). This causes dpois to returns 0 depending on values of theta1, theta2 and theta3. My first question is: is there anyway around this? My second question is regarding the likelihood function value. Since dpois returns 0, my log likelihood becomes infinite and causes problem when optimization. Is there any smart way to avoid this? Thank you, rc __ 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] Fast optimizer
Hi, I'm using optim with box constraints to MLE on about 100 data points. It goes quite slow even on 4GB machine. I'm wondering if R has any faster implementation? Also, if I'd like to impose equality/nonequality constraints on parameters, which package I should use? Any help would be appreciated. Thank you. rc __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Fast optimizer
Ok. I have the following likelihood function. L <- p*dpois(x,a)*dpois(y,b+c)+(1-p)*dpois(x,a+c)*dpois(y,b) where I have 100 points of (x,y) and parameters c(a,b,c,p) to estimate. Constraints are: 0 < p < 1 a,b,c > 0 c < a c < b I construct a loglikelihood function out of this. First ignoring the last two constraints, it takes optim with box constraint about 1-2 min to estimate this. I have to estimate the MLE on about 200 rolling windows. This will take very long. Is there any faster implementation? Secondly, I cannot incorporate the last two contraints using optim function. Thank you, rc On Thu, Oct 29, 2009 at 9:02 PM, Ravi Varadhan wrote: > > You have hardly given us any information for us to be able to help you. Give > us more information on your problem, and, if possible, a minimal, > self-contained example of what you are trying to do. > > Ravi. > > > Ravi Varadhan, Ph.D. > Assistant Professor, > Division of Geriatric Medicine and Gerontology > School of Medicine > Johns Hopkins University > > Ph. (410) 502-2619 > email: rvarad...@jhmi.edu > > > - Original Message - > From: R_help Help > Date: Thursday, October 29, 2009 8:24 pm > Subject: [R] Fast optimizer > To: r-help@r-project.org, r-sig-fina...@stat.math.ethz.ch > > >> Hi, >> >> I'm using optim with box constraints to MLE on about 100 data points. >> It goes quite slow even on 4GB machine. I'm wondering if R has any >> faster implementation? Also, if I'd like to impose >> equality/nonequality constraints on parameters, which package I should >> use? Any help would be appreciated. Thank you. >> >> rc >> >> __ >> R-help@r-project.org mailing list >> >> PLEASE do read the posting guide >> 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] Sequential MLE on time series with rolling window
Hi, Assuming I have a time series on which I will perform rolling-window MLE. In other words, if I stand at time t, I'm using points t-L+1 to t for my MLE estimate of parameters at time t (here L is my rolling window width). Next, at t+1, I'll do the same. My question is that is there anyway to avoid performing MLE each time like does the above. My impression is that rolling from point t to t+1, the likelihood function is equivalent to cutting out point t-L+1 and add back likelihood at point t+1. Is there any smart way to sequentially update the MLE instead of brute force calculation every time? Any suggestion or reference would be appreciated. Thank you. rc __ 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] Fitting ACD model
Hi - I'm wondering if there is any existing package in R that facilitates ACD (Autoregressive Conditional Duration) model. I googled. Apparently there were packages dynamo and fACD. However, they do not seem to be on any CRAN. I am wondering if anyone could point me to an existing package (if any). Thank you. -r __ 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] Need help installing xts package
Hi,
I am trying to install package xts. I am using OPENSUSE 11 64 bit. When I
invoke:
install.pacakges("xts",dependencies=T)
I received a lot of "ERROR: compilation failed for package" errors when R
tries to install xts dependencies. I do not understand this behavior.
I notice one thing. It complains.
ERROR: compilation failed for package 'Hmisc'
** Removing '/home/chibi/R/x86_64-unknown-linux-gnu-library/2.8/Hmisc'
* Installing *source* package 'Design' ...
** libs
WARNING: R include directory is empty -- perhaps need to install R-devel.rpm
or similar
Is it because I need R-devel.rpm? What is it? I search google around and
cannot find any explanation. Thank you.
adschai
[[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] Smart way to check if a column with certain name exists in a data.frame
Hi, I have a matrix or data.frame (df). When I would like to check if a column named "x" exists, it becomes my habit to do the following: if (length(which(colnames(df)=="x"))>0) Is there any smarter way to do this? Thank you. adschai [[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] lapply across using multiple columns
Hi, To set a simple an clear picture of what I'd like to do, here is an example. I would like to calculate mean of column A and B bucketed by factor in column X in a data.frame. lapply or aggregate operates on a single column. I found dapply function in some package which doesn't seem to exist any more. Is there anyway that I can accomplish something like this? Thank you so much in advance. adschai [[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] Mahalanobis distance with mixed data
Hi, I can see that there is a function for mahalanobis distance in stat package. But I'm wondering if there is any implementation that can take in a mix of ordinal data and numerical values together? If no, would you mind pointing me to some reference for a methodology to do this? Thank you in advance. Rob __ 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] Need help to understand integrate function
Hi,
I'm running into a wall when trying to use the integrate function. I
have the following setting:
powerLaw2 <- function(x,l1,l2,c0,t0) {
idx <- which(x <= 0);
if (length(idx) > 0) {
x[idx] <- 0;
}
xl <- (-l1+l2)*log(x/t0);
L <- log(c0)-l1*log(x)-log(1+exp(xl));
L <- exp(L);
return(L);
}
plCDF2 <- function(x,l1,l2,c0,t0) {
print(c(l1,l2,c0,t0));
cdf <-
integrate(f=powerLaw2,lower=x,upper=Inf,l1=l1,l2=l2,c0=c0,t0=t0,stop.on.error=T);
return(cdf$value);
}
I know that as long as my l1 > 1 and l2 < 1 the integration above
should converge. However, this doesn't seem to be the case when I call
integrate function with the following parameter:
Browse[1]>
integrate(f=powerLaw2,lower=1e-09,upper=Inf,l1=1.8980185,l2=-0.0804259,c0=1,t0=259.78,stop.on.error=T)
failed with message 'the integral is probably divergent'
Would you please shed some light? And how could I prevent this
situation? Thank you.
Rgds,
Rob
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Need help to understand integrate function
Hi Ravi - Thank you. I'm wondering how the magic happens here with
options(digits=10)? And the most important point, I do not quite
understand why divergent error could occur. In this case, the function
is analytically integrable. Any insight would be greatly appreciated
so that I learn something out of this. Thank you.
Rob
On Wed, Aug 11, 2010 at 11:46 PM, Ravi Varadhan wrote:
> Try this:
>
> options(digits=10)
>
>> integrate(f=powerLaw2, lower=0, upper=Inf, l1=1.8980185, l2=-0.0804259,
>> c0=1, t0=259.78, rel.tol=1.e-10)
> 0.01089019946 with absolute error < 3.7e-11
>
>
> Ravi.
>
>
> Ravi Varadhan, Ph.D.
> Assistant Professor,
> Division of Geriatric Medicine and Gerontology
> School of Medicine
> Johns Hopkins University
>
> Ph. (410) 502-2619
> email: rvarad...@jhmi.edu
>
>
> - Original Message -
> From: R_help Help
> Date: Wednesday, August 11, 2010 9:44 pm
> Subject: [R] Need help to understand integrate function
> To: r-help@r-project.org
>
>
>> Hi,
>>
>> I'm running into a wall when trying to use the integrate function. I
>> have the following setting:
>>
>> powerLaw2 <- function(x,l1,l2,c0,t0) {
>>
>> idx <- which(x <= 0);
>> if (length(idx) > 0) {
>> x[idx] <- 0;
>> }
>>
>> xl <- (-l1+l2)*log(x/t0);
>> L <- log(c0)-l1*log(x)-log(1+exp(xl));
>> L <- exp(L);
>> return(L);
>> }
>>
>> plCDF2 <- function(x,l1,l2,c0,t0) {
>>
>> print(c(l1,l2,c0,t0));
>> cdf <-
>> integrate(f=powerLaw2,lower=x,upper=Inf,l1=l1,l2=l2,c0=c0,t0=t0,stop.on.error=T);
>> return(cdf$value);
>>
>> }
>>
>> I know that as long as my l1 > 1 and l2 < 1 the integration above
>> should converge. However, this doesn't seem to be the case when I call
>> integrate function with the following parameter:
>>
>>
>> Browse[1]>
>> integrate(f=powerLaw2,lower=1e-09,upper=Inf,l1=1.8980185,l2=-0.0804259,c0=1,t0=259.78,stop.on.error=T)
>> failed with message 'the integral is probably divergent'
>>
>> Would you please shed some light? And how could I prevent this
>> situation? Thank you.
>>
>> Rgds,
>>
>> Rob
>>
>> __
>> r-h...@r-project.org mailing list
>>
>> PLEASE do read the posting guide
>> 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.
