[R] Area() artefacts??
Hello, everybody
I run the following program, and depending on the size of eps I get
different results.
With eps=1e-05, the program calculates wrong values for x=65:67 and
others. The program runs fine with eps=1e-07. Why is this so?
Also, I am using area() instead of integrate() because I cannot make
integrate to work, especially with imaginary numbers. Maybe someone
can show me how to use integrate in this particular code? THanks in
advance!
#Computes the p.m.f. via (10.53) of the number of i.i.d. Ber(p) trials
#required until m consecutive successes occur.
#Requires MASS package
#==#
consecpmf - function(xvec, m, p, eps=1e-05){
library(MASS)
f-numeric()
for(j in seq(xvec)){
x - xvec[j]
f[j] - area(fun, -pi, pi, limit=1000, eps=eps, x, m, p)
}
f-Re(f)
round(f,4)
}
fun - function(t,x,m,p){
I - exp(-1i*t*x)*cf(t,m,p)/(2*pi)
I
}
cf - function(t,m,p){
q - 1-p
if(m==1) {g - p*exp(1i*t)/(1-q*exp(1i*t))} else
{kk - exp(1i*t)*cf(t, m-1, p); g - (p*kk)/(1-q*kk)}
g
}
#===TESTING#
consecpmf(seq(0,200),3,0.5)
__
R-help@stat.math.ethz.ch 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 split() several times in a row?
Hi, fellow R users. I have a question about sapply and split combination. I have a big dataframe (4 observations, 21 variables). First variable (factor) is date and it is in format 8.29.97, that is, I have monthly data. Second variable (also factor) has levels 1 to 6 (fractiles 1 to 5 and missing value with code 6). The other 19 variables are numeric. For each month I have several hunder observations of 19 numeric and 1 factor. I am normalizing the numeric variables by dividing val1 by val2, where: val1: (for each month, for each numeric variable) difference between mean of ith numeric variable in fractile 1, and mean of ith numeric variable in fractile 5. val2: (for each month, for each numeric variable) standard deviation for ith numeric variable. Basically, as far as I understand, I need to use split() function several times. To calculate val1 I need to use split() twice - first to split by month and then split by fractile. Is this even possible to do (since after first application of split() I get a list)?? Is there a smart way to perform this normalization computation? My knowledge of R is not so advanced, but I need to know an efficient way to perform calculations of this kind. Would really appreciate some help from experienced R users! Regards, S -- Laziness is nothing more than the habit of resting before you get tired. - Jules Renard (writer) Experience is one thing you can't get for nothing. - Oscar Wilde (writer) When you are finished changing, you're finished. - Benjamin Franklin (Diplomat) __ R-help@stat.math.ethz.ch 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] Additional args to fun in integrate() not found?
Hello, fellow Rdicts,
I have the code for the program below. I need to integrate a function
of x and p. I use integrate to integrate over x and pass p as
an additional argument. p is specified and given default value in
the argument list. Still, integrate() cannot read p, unless I
explicitly insert a numeric value in the integrate() argument list.
And when I do that, I get the right result, but still some warnings.
Please, help me with these problems:
1) why is p not recognized?
2) what are these warning messages?
PROGRAM CODE:
---
#THIS LIBRARY IS NEEDED FOR THE INCOMPLETE GAMMA FUNCTION
library(zipfR)
#
gedCDF = function(yvec, p=2, mu=0, scale=1, numint=0)
{
#-
#Setting k to sqrt(2) and the GED with p=2 coincides with standard normal.
#Set k=1 and GED with p=1 coincides with Laplace.
k-sqrt(2)
#k-1
scale-scale*k
zvec-(yvec-mu)/scale
cdf-matrix(0, length(zvec),1)
for(i in 1:length(zvec))
{
z-zvec[i]
if(numint==0)
{
if(z=0)
{
t-0.5*(1-1/gamma(1/p)*Igamma((1/p),(-z)^p,lower=TRUE))
}
else
{
t-1-(0.5*(1-1/gamma(1/p)*Igamma((1/p),(z)^p,lower=TRUE )))
}
}
else
{
t-integrate(geddenstandard, -35, z, subdivisions=1000,
rel.tol=100*.Machine$double.eps, abs.tol=rel.tol,
stop.on.error=TRUE, keep.xy=FALSE, aux=NULL,p)
}
cdf[i]-t
}
cdf
}
#-
geddenstandard = function(z,p)
{
f-p/(2*gamma(1/p))*exp(-abs(z)^p)
}
---
If I run with this definition I get the following error message and abort:
gedCDF(c(1,2,3,4,5), numint=1)
Error in eval(expr, envir, enclos) : ..1 used in an incorrect context,
no ... to look in
If I replace p in integrate() with 2, I get correct answers, but
still some warning messages:
gedCDF(c(1,2,3,4,5), numint=1)
[[1]]
[1] 0.8413447
[[2]]
[1] 0.9772499
[[3]]
[1] 0.9986501
[[4]]
[1] 0.683
[[5]]
[1] 0.997
Warning messages:
1: number of items to replace is not a multiple of replacement length
2: number of items to replace is not a multiple of replacement length
3: number of items to replace is not a multiple of replacement length
4: number of items to replace is not a multiple of replacement length
5: number of items to replace is not a multiple of replacement length
---
Do I get these warnings because I define cdf as a matrix and the
output-cdf is a list?
Please, help me with these!
Email to my gmail account, please: [EMAIL PROTECTED]
THanks in advance
Sergey
__
R-help@stat.math.ethz.ch 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] Additional: to integrate()
Hi again, people I found that there is some problem with rel.tol agrument in the integrate() function. It is not getting found. Maybe because it is a call to .Machine$double.eps? If I specify rel.tol as an argument in gedCDF, everything works, except for those warning messages. Strange. -- Laziness is nothing more than the habit of resting before you get tired. - Jules Renard (writer) Experience is one thing you can't get for nothing. - Oscar Wilde (writer) When you are finished changing, you're finished. - Benjamin Franklin (President) __ R-help@stat.math.ethz.ch 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] Automated figure production
Hello, everybody
Two questions:
1) I am new to maillists, and particularly to r-help maillist. Where
specifically do I go online to see my question and answers to it??? If
I use the searchable archives, or archives by Robert King, I see my
question but not the answers, though I know that at least one persion
posted the answer to r-help. Why do not I see the answers?
2)
I need to produce 104 figures. In each graphic window I want to plot 8
figures. How do I automate the process so that it opens 13 separate
graphic windows in R Gui and plots the figures? (I then paste them in
Word, one by one). Also, how do I save all 104 figures to one PDF
file?
Here is the code to produce the figures, and at this point I change
the j-index in the for-loop by hand, produce 8 figures, copy them to
Word, then change j to 9:16 etc...
#HOW TO PRODUCE ALL 104 GRAPHS AT ONCE
old.par-par(no.readonly=TRUE)
par(mfrow=c(4,2))
for(j in 1:8) {
plot(Cleaned[zz[,j],4], type=b, main=long.names[j], xlab=Month,
ylab=MFR, col=blue)
abline(h=0, col=red)
}
par(old.par)
Thanks for your help!
Sergey
--
Laziness is nothing more than the habit of resting before you get tired.
- Jules Renard (writer)
Experience is one thing you can't get for nothing.
- Oscar Wilde (writer)
When you are finished changing, you're finished.
- Benjamin Franklin (President)
__
R-help@stat.math.ethz.ch 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] Combining tapply() and cor.test()?
Hello, fellow R-users. Let me describe the setup first. I have a data.frame, a sample of which is reported below: Company.Name Periods Returns MFR.Factor 350 Wartsila Oyj A 1996-07-31 6.82 0.02 351Custodia Holding AG 1996-07-31 4.15-0.02 352 Wartsila Oyj 1996-07-31 7.73 0.09 353 GEA Group AG 1996-07-3110.12 0.04 354LEGRAND ORD 1996-07-31 -7.46-0.20 355 Mayr-Melnhof Karton AG 1996-07-31 4.71-0.05 356GEVAERT NPV 1996-08-30 NA NA 357NOKIA K FMA2.50 1996-08-30 7.65 0.03 358 Altadis S.A. 1996-08-30 7.65 0.55 359 Metrovacesa S.A. 1996-08-30 4.55-0.17 360 Oce N.V. 1996-08-309.43 0.23 The variable Periods is a date object, shows the month. Variables Returns and MFR.Factor are numeric. For each month the number of Returns and MFR.Factors varies, sometimes it is 350, sometimes 320 etc. What I need is to use cor.test(Returns, MFR.Factor,...) for each month, and produce a dataframe with columns: Period, cor.estimate, p.value. The simplest way would be with tapply() using variable Period as a factor, but tapply() only applies FUN to just one cell. What is the most painless way to achieve my objective? Thank you in advance for your help! Best, Sergey __ R-help@stat.math.ethz.ch 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.
