Re: [Rd] solving nonlinear equations
Have you tried writing a function to compute SS = sum of squares deviations between the the left and right hand sides of your three equations, then using 'optim'? See also Venables and Ripley (2002) Modern Applied Statistics with S, 4th ed. (Springer). hope this helps. Spencer Graves p.s. I don't see how it's obvious that 'a=-c'. Kjetil Brinchmann Halvorsen wrote: > HAKAN DEMIRTAS wrote: >> I can't seem to get computationally stable estimates for the following >> system: >> >> Y=a+bX+cX^2+dX^3, where X~N(0,1). (Y is expressed as a linear combination of >> the first three powers of a standard normal variable.) Assuming that E(Y)=0 >> and Var(Y)=1, one can obtain the following equations after tedious algebraic >> calculations: >> >> 1) b^2+6bd+2c^2+15d^2=1 >> 2) 2c(b^2+24bd+105d^2+2)=E(Y^3) >> 3) 24[bd+c^2(1+b^2+28bd)+d^2(12+48bd+141c^2+225d^2)]=E(Y^4)-3 >> >> Obviously, a=-c. Suppose that distributional form of Y is given so we know >> E(Y^3) and E(Y^4). In other words, we have access to the third and fourth >> raw moments. How do we solve for these four coefficients? I reduced the >> number of unknowns/equations to two, and subsequently used a grid approach. >> It works well when I am close to the center of the support, but fails >> miserably at the tails. Any ideas? Hopefully, there is a nice R function >> that does this. >> >> Hakan Demirtas >> >> >> [[alternative HTML version deleted]] >> >> __ >> R-devel@r-project.org mailing list >> https://stat.ethz.ch/mailman/listinfo/r-devel >> > This is really a question for r-help not r-devel. > > I was about to say that this was a question for a symbolic algebra > system, but first tried in MuPAD 4.0, and left the machine alone. > returning after 2 hours MuPAD was still grinding and I had to kill it. > > Kjetil > > __ > R-devel@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-devel __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] unique.default() drops names (PR#9130)
Gregor Gorjanc <[EMAIL PROTECTED]> writes:
> Thank you for the reply! I appologize for not reading the latest
> documentation - there was no word about droping names in 2.3.1. However,
> I do wonder why simple fix (as shown in previous mail) is not OK.
I see value in unique() keeping names and from what I understand
the documentation could be changed to match ;-)
I don't know if there are good reasons for dropping names from
vectors.
Given that unique is very commonly used, I think the way to make such
a change is in the C code, not at the R level. So in that sense, I
think the patch you sent is not ideal. Below is a patch to
do_duplicated that keeps names. Lightly tested. No doc included. I
would consider more testing and doc if there was interest.
+ seth
diff --git a/src/main/unique.c b/src/main/unique.c
index a3c7a87..d8d31fa 100644
--- a/src/main/unique.c
+++ b/src/main/unique.c
@@ -382,7 +382,7 @@ SEXP duplicated(SEXP x)
*/
SEXP attribute_hidden do_duplicated(SEXP call, SEXP op, SEXP args, SEXP env)
{
-SEXP x, dup, ans;
+SEXP x, xnames, dup, ans, ansnames;
int i, k, n;
checkArity(op, args);
@@ -410,25 +410,38 @@ SEXP attribute_hidden do_duplicated(SEXP
k++;
PROTECT(dup);
-ans = allocVector(TYPEOF(x), k);
-UNPROTECT(1);
+PROTECT(ans = allocVector(TYPEOF(x), k));
+xnames = getAttrib(x, R_NamesSymbol);
+if (xnames != R_NilValue)
+ansnames = allocVector(STRSXP, k);
+else
+ansnames = R_NilValue;
+UNPROTECT(2);
k = 0;
switch (TYPEOF(x)) {
case LGLSXP:
case INTSXP:
for (i = 0; i < n; i++)
- if (LOGICAL(dup)[i] == 0)
+ if (LOGICAL(dup)[i] == 0) {
+if (ansnames != R_NilValue)
+SET_STRING_ELT(ansnames, k, STRING_ELT(xnames, i));
INTEGER(ans)[k++] = INTEGER(x)[i];
+}
break;
case REALSXP:
for (i = 0; i < n; i++)
- if (LOGICAL(dup)[i] == 0)
+ if (LOGICAL(dup)[i] == 0) {
+if (ansnames != R_NilValue)
+SET_STRING_ELT(ansnames, k, STRING_ELT(xnames, i));
REAL(ans)[k++] = REAL(x)[i];
+}
break;
case CPLXSXP:
for (i = 0; i < n; i++)
if (LOGICAL(dup)[i] == 0) {
+if (ansnames != R_NilValue)
+SET_STRING_ELT(ansnames, k, STRING_ELT(xnames, i));
COMPLEX(ans)[k].r = COMPLEX(x)[i].r;
COMPLEX(ans)[k].i = COMPLEX(x)[i].i;
k++;
@@ -436,22 +449,33 @@ SEXP attribute_hidden do_duplicated(SEXP
break;
case STRSXP:
for (i = 0; i < n; i++)
- if (LOGICAL(dup)[i] == 0)
+ if (LOGICAL(dup)[i] == 0) {
+if (ansnames != R_NilValue)
+SET_STRING_ELT(ansnames, k, STRING_ELT(xnames, i));
SET_STRING_ELT(ans, k++, STRING_ELT(x, i));
+}
break;
case VECSXP:
for (i = 0; i < n; i++)
- if (LOGICAL(dup)[i] == 0)
- SET_VECTOR_ELT(ans, k++, VECTOR_ELT(x, i));
+ if (LOGICAL(dup)[i] == 0) {
+if (ansnames != R_NilValue)
+SET_STRING_ELT(ansnames, k, STRING_ELT(xnames, i));
+SET_VECTOR_ELT(ans, k++, VECTOR_ELT(x, i));
+}
break;
case RAWSXP:
for (i = 0; i < n; i++)
- if (LOGICAL(dup)[i] == 0)
+ if (LOGICAL(dup)[i] == 0) {
+if (ansnames != R_NilValue)
+SET_STRING_ELT(ansnames, k, STRING_ELT(xnames, i));
RAW(ans)[k++] = RAW(x)[i];
+}
break;
default:
UNIMPLEMENTED_TYPE("duplicated", x);
}
+if (ansnames != R_NilValue)
+setAttrib(ans, R_NamesSymbol, ansnames);
return ans;
}
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
[Rd] BOD causes error in 2.4.0
Using "R version 2.4.0 Under development (unstable) (2006-08-08 r38825)"
on Windows XP and starting in a fresh session we get an error if we type BOD.
(There is no error in "Version 2.3.1 Patched (2006-06-04 r38279)".)
> BOD
Error in data.frame(Time = c("1", "2"), demand = c(" 8.3", "10.3"),
check.names = FALSE, :
row names contain missing values
In addition: Warning message:
corrupt data frame: columns will be truncated or padded with NAs in:
format.data.frame(x, digits = digits, na.encode = FALSE)
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] sort() generic? [Re: Suspicious behaviour of sort on POSIXct ..]
With all due respect, I don't think that the change to the sort function is "a
good change". New versions should be backward compatible whenever possible,
i.e. unless there's a very compelling reason to force large amounts of old code
to break. In this case, the sort function could have had an additional
"unclass" argument with a default value of FALSE. That way, old code wouldn't
break and new code could be written with unclass set to TRUE for speed. I hope
this will be the case for 2.3.4
___
Alexander Dannenberg
Geode Capital Management, LLC
53 State Street . Boston . MA . 02109
Tel 617.392.8123 Cell 617.372.6805
[EMAIL PROTECTED]
This e-mail, and any attachments hereto, are intended for us...{{dropped}}
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
[Rd] S4 Classes
Hello All, I'm trying to convince someone that they should transition a large project to use S4 instead of S3 classes. Does anyone have any good citations? Thanks! -- Dan Gerlanc -- Daniel Gerlanc Williams College '07 [[alternative HTML version deleted]] __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
[Rd] rownames() problem with 0-extent arrays (PR#9136)
[R 2.3.1 on Windows XP] Hello, The 'do.NULL=FALSE' option of 'rownames' and 'colnames' does not work as documented when the array has a 0-extent dimension: tmp <- matrix(1:15, ncol=3, dimnames=list(letters[1:5], LETTERS[1:3])) rownames(tmp) # Fine rownames(tmp[0,,drop=FALSE]) # NULL (okay, but see question below) rownames(tmp[0,,drop=FALSE], do.NULL=FALSE) # "row" The last line returns a character vector of length 1, "row", rather than character(0) as implied by the second paragraph of the details section of ?rownames: "If 'do.NULL' is 'FALSE', a character vector (of length 'NROW(x)' or 'NCOL(x)') is returned in any case, prepending 'prefix' to simple numbers, if there are no dimnames or the corresponding component of the dimnames is 'NULL'." -- A separate but related question (not a bug) has to do with the handling of dimnames for 0-extent arrays. If an array starts out with non-NULL dimnames and is subscripted down to have 0-extent, the corresponding component is set to NULL, rather than character(0): dimnames(tmp[0,,drop=FALSE]) # first component is NULL Why is that? Note that this is different from how data frames are handled: tmp.df <- data.frame(tmp) dimnames(tmp.df[0,]) # first component is character(0) The description in ?dimnames seems to allow character(0) components of dimnames, "The dimnames of a matrix or array can be 'NULL' or a list of the same length as 'dim(x)'. If a list, its components are either 'NULL' or a character vector the length of the appropriate dimension of 'x'." yet the implementation seems to resist them: tmp2 <- tmp[0,,drop=FALSE] dimnames(tmp2) <- list(character(0), LETTERS[1:3]) dimnames(tmp2) # still NULL So my question is whether it would be reasonable to change 'dimnames' and/or "dimnames<-" to use (or at least allow) character(0) rather than NULL for array dimensions of 0 extent. Rich Raubertas Merck & Co. __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] unique.default() drops names (PR#9130)
Seth Falcon fhcrc.org> writes: > I see value in unique() keeping names and from what I understand > the documentation could be changed to match > > I don't know if there are good reasons for dropping names from > vectors. > > Given that unique is very commonly used, I think the way to make such > a change is in the C code, not at the R level. So in that sense, I > think the patch you sent is not ideal. Below is a patch to > do_duplicated that keeps names. Lightly tested. No doc included. I > would consider more testing and doc if there was interest. Thank you Seth for your time on this. I also think that there is value in keeping names and I agree that C code is the best place to do the fix. Gregor __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] S4 Classes
Is there a capability that you would like for the package which could be achieved only if the package was transitioned to S4? If so, explain this to the author. If not, why ask them to change? Gordon >[Rd] S4 Classes >Daniel Gerlanc dgerlanc at gmail.com >Thu Aug 10 23:37:15 CEST 2006 > >Hello All, > >I'm trying to convince someone that they should transition a large project >to use S4 instead of S3 classes. Does anyone have any good citations? >Thanks! > >-- Dan Gerlanc > >-- >Daniel Gerlanc >Williams College '07 __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] S4 Classes
Gordon Smyth <[EMAIL PROTECTED]> writes: > Is there a capability that you would like for the package which could > be achieved only if the package was transitioned to S4? If so, > explain this to the author. If not, why ask them to change? 'achieved only if' sounds a bit strong. I've used S4 objects and found them useful for, among other things: * familiar object-oriented reasons (e.g., data structure and code reuse in derived classes and methods), * flexibility over S3 classes (e.g., slot and validity checking to at least partly relieve me of the checks, familiar to S3 programmers, at the top of each function to ensure that the data conforms to expectation) and methods (e.g., dispatch on other than the first argument), * interoperability with other packages, * object introspection, e.g., validating and updating serialized objects against current class definitions, exchanging complicated data structures with other languages, * the greater structure S4 classes and methods encourage in complicated programs, and * fun challenges to exploiting the unique (relative to, say, Java) features of S4 class structure and method dispatch. As Gordon says, though, package requirements should drive implementation decisions. Martin -- Bioconductor > Gordon > >>[Rd] S4 Classes >>Daniel Gerlanc dgerlanc at gmail.com >>Thu Aug 10 23:37:15 CEST 2006 >> >>Hello All, >> >>I'm trying to convince someone that they should transition a large project >>to use S4 instead of S3 classes. Does anyone have any good citations? >>Thanks! >> >>-- Dan Gerlanc >> >>-- >>Daniel Gerlanc >>Williams College '07 > > __ > R-devel@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-devel __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
[Rd] R's built-in Script Editor Bug (PR#9137)
Full_Name: Richard Cai
Version: R 2.3.1
OS: Windows XP SP2
Submission from: (NULL) (220.233.184.74)
My R scrpit code was interrupted when I was trying to load it into R's built-in
script editor.
R code for question 4
dat <- read.table("Drosophila.txt",header=T)
library(lattice)
trellis.device(color=F)
attach(dat)
## scatter plot of the data
data_plot <- xyplot(Weight~density)
print(data_plot)
dev.copy2eps(file="D:/CaiYu/UNE/Statistics/Stat
356/Assignment1/q4_data_plot.eps",horizontal=F)
# plot of group means against density ##
group_mean <- sapply(split(Weight,density),mean)
density_level <- as.numeric(levels(as.factor(density)))
group_plot <- xyplot(group_mean ~ density_level)
windows()
print(group_plot)
dev.copy2eps(file="D:/CaiYu/UNE/Statistics/Stat
356/Assignment1/q4_group_plot.eps",horizontal=F)
# fitted model ###
fitted_model <- lm(Weight ~density, data=dat)
print(anova(fitted_model))
print(summary(fitted_model)$coefficients)
residual_SS <- anova(fitted_model)$Sum[2]
df_residual_SS <- anova(fitted_model)$Df[2]
### One-way analysis of variance#
aov <- aov(Weight ~ as.factor(density), data=dat)
print(anova(fitted_model))
print(summary(aov))
pure_error_SS <- anova(aov)$Sum[2]
df_pure_error_SS <- anova(aov)$Df[2]
##lack of fittness ###
lack_of_fittness_SS <- residual_SS - pure_error_SS
df_lack_of_fittness_SS <- df_residual_SS - df_pure_error_SS
F_lack_of_fittness <-
(lack_of_fittness_SS/df_lack_of_fittness_SS)/(pure_error_SS/df_pure_error_SS)
P_value_lack_of_fittness <- pf(F_lack_of_fittness, df_lack_of_fittness_SS,
df_pure_error_SS, lower.tail=F)
print(P_value_lack_of_fittness)
##testing hypothesis beta1=0 ###
print(summary(fitted_model)$coefficients)
### prediction for density=15 #
new1 <- data.frame(density=15)
print(predict(lm(Weight~density), new1, se.fit = TRUE))
print(predict(lm(Weight~density), new1, interval="prediction"))
print(predict(lm(Weight~density), new1, interval="confidence"))
### prediction for density=40#
new2 <- data.frame(density=40)
print(predict(lm(Weight~density), new2, se.fit = TRUE))
print(predict(lm(Weight~density), new2, interval="prediction"))
print(predict(lm(Weight~density), new2, interval="confidence"))
# prediciton plot #
windows()
pred.w.plim <- predict(fitted_model, data.frame(density=seq(1,40,1)),
interval="prediction")
pred.w.clim <- predict(fitted_model, data.frame(density=seq(1,40,1)),
interval="confidence")
larval_density <- seq(1,40,1)
matplot(larval_density,cbind(pred.w.clim, pred.w.plim[,-1]),
lty=c(1,2,2,3,3), col=rep(1,5), type="l", ylab="predicted Weight")
dev.copy2eps(file="D:/CaiYu/UNE/Statistics/Stat
356/Assignment1/q4_prediction_plot.eps",horizontal=F)
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
