Re: [R] multicore and mclapply problem in calculation server
Hi Patrick,
Thanks for posting. The wrapper function works as it should. If you do not
specify return the function returns the last object declared, and this is
exactly what I intended. Doing some research on my own I found out that the
problem was here:
x <- c(1: length(linH))
stat.matrix <- do.call(rbind, mclapply(x, wrapper, mc.cores= 6))
Since linH = readLines(conH, n= N.perm) and N.perm was set to 999,
apparently mclapply sends _simultaneously_ the 999 jobs to the cpus
specified in mc.cores. Even if you use more than 6, as in the example, the
systems gets quickly overflowed. So the solution is to split length(linH)
into blocks equal to the number of cpus and the message to end-users is "be
careful with parallel processing functions in R".
Cheers
Juan Antonio Balbuena
El 19/12/2013 22:03, p_connolly escribió:
On Wed, 18-Dec-2013 at 09:28AM +0100, Juan Antonio Balbuena wrote:
|>
|>Hello
|>I am using package multicore for parallel computing in a Altix
UltraViolet
|>1000 server with 64 CPUs and 960 GB of RAM memory. Access is managed
by
|>means of a SGE queue system. This is the first time I am using
parallel
|>computing and my experience with supercomputers is quite limited. So
any
|>help will be much, much appreciated.
|>My experiment consists of a number of runs (N.runs) each involving a
number
|>of permutations (N. perms). (An excerpt of the code is included
below.) The
|>permutations are very time consuming and I am using mclapply to
distribute
|>the job among a given number of cpus (usually 12 to 24). The problem
is that
|>the system administrators notice that threads keep increasing as the
program
|>is executed to the point that they compromise the functioning of the
whole
|>system and have to abort the job.
So, does that mean you get nothing written to 1MH_30.tre? I'd be
surprised if you did get anything. Though there's lots of stuff
happening I know nothing about, I'm pretty sure there's an issue with
your wrapper() function.
[...]
|> wrapper <- function (x) { # THIS FUCNTION IS SUPPOSED TO
BE
|>PARALLELIZED (SEE BELOW)
|>treeH <- read.tree(text=linH[x])
|>treeP <- read.tree(text=linP[x])
|>mrcaL <- MRCALink.simul (treeH, treeP, HP)
|>stat.matrix[x,] <- three.stat(mrcaL)
|>}
Nothing is being returned, so your calls to rbind will have nothing to
put together. So what they end up trying to do, I've no idea. That
function probably needs a final line 'stat.matrix' before the '}'
One thing I discovered with mclapply is that to use the browser()
function, it's necessary to make what is x in your example of length 1
(and probably mc.cores to 1 also) before it's possible to examine
what's happening at various parts of the function being debugged.
Judicious use of cat() and Sys.time() to display what's happening at
various stages is also helpful.
HTH
--
~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.
___Patrick Connolly
{~._.~} Great minds discuss ideas
_( Y )_Middle minds discuss events
(:_~*~_:)Small minds discuss people
(_)-(_) . Anon
~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.
--
Dr. Juan A. Balbuena
Marine Zoology Unit
Cavanilles Institute of Biodiversity and Evolutionary Biology
University of
Valencia
[1]http://www.uv.es/~balbuena
P.O. Box 22085
[2]http://www.uv.es/cavanilles/zoomarin/index.htm
46071 Valencia, Spain
[3]http://cetus.uv.es/mullpardb/index.html
e-mail: [4]j.a.balbu...@uv.estel. +34 963 543 658fax +34 963 543 733
NOTE! For shipments by EXPRESS COURIER use the following street address:
C/ Catedrático José Beltrán 2, 46980 Paterna (Valencia), Spain.
References
1. http://www.uv.es/%7Ebalbuena
2. http://www.uv.es/cavanilles/zoomarin/index.htm
3. http://cetus.uv.es/mullpardb/index.html
4. mailto:j.a.balbu...@uv.es
__
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] ftable and data.frame
Hi, I used this command to produce a table: (tab1 = ftable(SEX, ESTCIV, Q1)) Q1 B L M N SEXOESTCIV F A 11 13 4 2 E 1 0 0 0 M A 5 0 3 1 E 0 0 0 0 but I need something like: SEXOESTCIVB L M N FA 11 13 4 2 FE 1 0 0 0 M A 5 0 3 1 M E 0 0 0 0 How can I get it? I need this format to use ordinal logistic regression and I have many tables. Thanks, Silvano. --- Este email está limpo de vírus e malwares porque a proteção do avast! Antivírus está ativa. [[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] Fitting particle size analysis data
Hi I made a simple spredsheet for PSD using Rosin Rammler equation and I am lazy to transform it to R. However for single purpose you can use nls. Reverse your cumulative values PSD$cum<-cumsum(PSD$ret) plot(PSD$size, PSD$cum) fit<-nls(cum~ exp(-((size/r)^gama))*100, data=PSD, start=c(r=80, gama=2)) summary(fit) Formula: cum ~ exp(-((size/r)^gama)) * 100 Parameters: Estimate Std. Error t value Pr(>|t|) r 88.9664 2.3360 38.09 2.35e-07 *** gama 2.5435 0.2244 11.33 9.36e-05 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 3.411 on 5 degrees of freedom Number of iterations to convergence: 7 Achieved convergence tolerance: 1.612e-06 lines(PSD$size, predict(fit)) Regards Petr > -Original Message- > From: r-help-boun...@r-project.org [mailto:r-help-bounces@r- > project.org] On Behalf Of Zorig Davaanyam > Sent: Friday, December 20, 2013 2:01 AM > To: r-help@r-project.org > Subject: [R] Fitting particle size analysis data > > Hi all, > > How do you fit a sieve analysis data to a statistical function? > I have many sieve analysis data of crushed rocks and I'd like to find > out which statistical distributions describe the particular particle > size distributions (PSD) the best. So basically I need to find fitted > parameters to statistical distributions (mostly weibull and truncated > lognormal). > Here is an example of particle size (in microns) versus percent weight > retained. > Sieve size Wt% Cumulative passing% > +250 0.1 99.9 > -250+1802.9 97 > -180+1259.5 87.5 > -125+90 21.266.3 > -90+6329.436.9 > -63+45 26 10.9 > -45 10.9 > > PSD<- > data.frame(size=c(250,180,125,90,63,45,0),retained=c(0.1,2.9,9.5,21.2,2 > 9.4,26,10.9),cumulative=c(99.9,97,87.5,66.3,36.9,10.9,0)) > > The above example is truncated to 350micron and I can't have particles > with minus dimension. Any help will be greatly appreciated. > > Thank you, > > Zorig > > [[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] Estimating and predicting using "segmented" Package
I don't know if this is still relevant for you, but probably for someone
else.
I ran into the same problem an found two dirty workarounds that can help
in many cases (but not all). Both can also be combined:
1. repeat the call to "segmented" with different seeds and different
initial estimates of the breakpoints (psi). Wrap this with "try" to catch
the error:
res= NULL
while (class(res)[1]!="segmented")
{
some_seed= runif(1) #you may also systematically vary this
some_psi = runif(2) #you may also systematically vary this.
#Better tailor this to the distribution of x values in
your_data
model_lm = lm(y~x, data=your_data)
res = try(segmented(obj=model_lm, seg.Z=~x, model=TRUE,
psi=list(x=some_psi),
control=seg.control(seed=some_seed)), silent=TRUE)
}
2. insert dummy points between your actual data to reduce the chance of
"only 1 datum in an interval". When you supply these dummy data with low
weights, their effect on the outcome should be low:
model_lm = lm(y~x, data=data_with_dummyentries,
weights=data_with_dummyentries$weight)
res=segmented(obj=model_lm, seg.Z=~x, model=TRUE, psi=list(x=c(10,20)))
Hope it helps,
Melwin
__
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,
excuse me
where did I mistake in these procedures?
(I want to fit these models to my data)
install.packages('drc',rep='http://cran.um.ac.ir')
install.packages('car',rep='http://cran.um.ac.ir')
install.packages('abind',rep='http://cran.um.ac.ir')
library(car)
library(abind)
library(alr3)
library(gtools)
library(magic)
library(plotrix)
library(drc)
library(lattice)
library(MASS)
library(nlme)
library(methods)
library(stats)
library(stats4)
dose<-c(0,1,2,3)
SRMKN<-scan()
100
64.82649842
50.15772871
9.936908517
a=data.frame(dose,SRMKN)
par(mfrow=c(3,3))
##fit the logistic model##
a.m1 <- drm(SRMKN ~ dose, data = a, fct = L.4())
summary(a.m1)
plot(a.m1)
logLik(a.m1)
ED(a.m1, 50, interval='delta')
##fitting the linear model##
a.m2<-lm(SRMKN ~ dose, data = a)
plot(a.m2)
##fitting the quadratic model##
a.m3<-lm(SRMKN ~ dose^2, data = a)
plot(a.m3)
[[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] Execute Excel-Macro using R
Hi, How to execute ms-Excel Macro(*.xlsm) using R function ? I tried but not get. There are method to call from R function from Excel macro, but i need Excel macro to execute from R. Is it possible ? I am using Excel-2010 macro. Thanks in advance, Antony. -- View this message in context: http://r.789695.n4.nabble.com/Execute-Excel-Macro-using-R-tp4682533.html Sent from the R help mailing list archive at Nabble.com. __ 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] Classification of polynomial regression: simple or multiple (conceptual doubt)
Dear list, I'm posting in the R-help list due to: - Not knowing a better place for it; - I would like to know the opinion of more specialized people. What is the best place to classify polynomial regressions (Y = bo + b1X + b2X^2 + ... + bnX^n): single or multiple linear regression? Regards, -- ///\\\///\\\///\\\///\\\///\\\///\\\///\\\///\\\ Jose Claudio Faria Estatistica UESC/DCET/Brasil joseclaudio.faria at gmail.com Telefones: 55(73)3680.5545 - UESC 55(73)9100.7351 - TIM 55(73)8817.6159 - OI ///\\\///\\\///\\\///\\\///\\\///\\\///\\\///\\\ __ 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] ftable and data.frame
What does your original data look like? I seems to me that it would be better to give us that information since you may not want to use ftable() at all. Ideally you should give us the original data or a sample of it. If it's confidential just replace the actual values with fake data. The structure of the data is more important than the values. Please use dput() (see ?dput for information) to supply the data. Just do : dput(mydata) and then paste the result into the email. John Kane Kingston ON Canada > -Original Message- > From: silv...@uel.br > Sent: Fri, 20 Dec 2013 11:01:26 -0200 > To: r-help@r-project.org > Subject: [R] ftable and data.frame > > Hi, > > I used this command to produce a table: > > (tab1 = ftable(SEX, ESTCIV, Q1)) > >Q1 B L M N > SEXOESTCIV > F A 11 13 4 2 > E 1 0 0 0 > M A 5 0 3 1 > E 0 0 0 0 > > but I need something like: > > > SEXOESTCIVB L M N > FA 11 13 4 2 > FE 1 0 0 0 > M A 5 0 3 1 > M E 0 0 0 0 > > How can I get it? > > I need this format to use ordinal logistic regression and I have many > tables. > > Thanks, > > Silvano. > > --- > Este email esta limpo de vmrus e malwares porque a protegco do avast! > Antivmrus esta ativa. > > > [[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. FREE 3D MARINE AQUARIUM SCREENSAVER - Watch dolphins, sharks & orcas on your desktop! __ 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] ftable and data.frame
Hi, Try: library(reshape2) dcast(as.data.frame(tab1), SEX+ESTCIV~Q1,value.var="Freq") ##not tested. A.K. On Friday, December 20, 2013 8:03 AM, "silv...@uel.br" wrote: Hi, I used this command to produce a table: (tab1 = ftable(SEX, ESTCIV, Q1)) Q1 B L M N SEXO ESTCIV F A 11 13 4 2 E 1 0 0 0 M A 5 0 3 1 E 0 0 0 0 but I need something like: SEXO ESTCIV B L M N F A 11 13 4 2 F E 1 0 0 0 M A 5 0 3 1 M E 0 0 0 0 How can I get it? I need this format to use ordinal logistic regression and I have many tables. Thanks, Silvano. --- Este email está limpo de vírus e malwares porque a proteção do avast! Antivírus está ativa. [[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.
Re: [R] Classification of polynomial regression: simple or multiple (conceptual doubt)
Dear Dr. José Faria, I think that the best category to put polynomial regressions is single regressions. Although, in polynomial regressions there are more then one term as in multiple regressions this is an adjustment consequence, not a design consequence. So, to me this is sufficient to justify the category. This is my point of view. Best regards, Gilson Carvalho 2013/12/20 Jose Claudio Faria > Dear list, > > I'm posting in the R-help list due to: > - Not knowing a better place for it; > - I would like to know the opinion of more specialized people. > > What is the best place to classify polynomial regressions (Y = bo + > b1X + b2X^2 + ... + bnX^n): single or multiple linear regression? > > Regards, > -- > ///\\\///\\\///\\\///\\\///\\\///\\\///\\\///\\\ > Jose Claudio Faria > Estatistica > UESC/DCET/Brasil > joseclaudio.faria at gmail.com > Telefones: > 55(73)3680.5545 - UESC > 55(73)9100.7351 - TIM > 55(73)8817.6159 - OI > ///\\\///\\\///\\\///\\\///\\\///\\\///\\\///\\\ > > __ > 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. > -- Gilson Correia de Carvalho Doutor em Ecologia *(Doctor in Ecology)* Mestre em Ecologia e Biomonitoramento* (Master in Ecology and Biomonitoring) * Bacharel em Ciências Biológicas *(Bachelor in Biological Sciences)* CV Lattes: http://lattes.cnpq.br/8361386734266580 Professor Assistente *(Assistant Professor)* Departamento de Biointeração *(Department of Biointeraction)* Instituto de Ciências da Saúde *(Health Sciences Institute)* Universidade Federal da Bahia *(Federal University of Bahia)* - Diretor Técnico *(Technical Director)* Holos Soluções Ambientais Ltda - Skype: bio_gilson GTalk: biogilson [[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] Using cbind to merge different variables
HI, May be this helps: set.seed(45) dat1 <- as.data.frame(matrix(sample(0:1,100*5,replace=TRUE),ncol=5)) dat1$Newvar <- 1*(!!rowSums(dat1)) A.K. Hello. I have a problem combining a number of variables. I have five columns with binary variables with the values 0 and 1. I would like to combine them into just one binary variable with 1 whenever just one of the other has value one, and 0 if none of them have value one. How can I do that? I tried to use cbind function, but for some reason I get results 1 even though all of the varuables included = 0 and for some rows I get 2. I tried changing the deparse.level but that doesn't seem to be the problem. I hope someone can help. Kind regards __ 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] Using assign with mapply
Trying to understand environments is not for the faint of heart. If lists do what you want then I would stick with a list and not worry about the environments. Most of the time that you deal with environments everything happens automatically behind the scenes and you don't need to worry about the details. If you want to learn more about environments there here is a start. Many of the ways of working with environments is the same as working with lists (you can access and assign with `[[` etc.) and there are `as.` functions to convert between them. The biggest difference is that environments use references instead of copies (powerful but dangerous). If you do something like mylist2 <- mylist mylist2$x <- 1 then mylist 2 will be a copy of mylist 1 and the value of x within mylist2 will be created or modified (but mylist will remain unchanged). However if you do the same thing with an environment then a copy is not made and the `x` variable in the original environment will be created or changed. This can be a big benefit when you have a large data object that needs to be passed to multiple functions, with an environment the data will never be copied, with a list or other object you may end up with multiple copies (though R is really good at not making copies when it does not need to, but sometimes you still end up with more copies than needed when R cannot tell if a copy is needed or not), but this is dangerous in that if you make any changes then the original is changed as well (the regular mechanism of making copies on changes protects the original). Environments also use hashing for name look ups which can speed things up when you have a large number of variables that you are accessing by name (but most usual cases are quick enough that you will not notice when using lists or other objects). Hope that helps, On Thu, Dec 19, 2013 at 10:17 AM, Julio Sergio Santana wrote: > Greg Snow <538280 gmail.com> writes: > >> >> The take home message that you should be learning from your struggles >> is to "Not Use The 'assign' Function!" and "Do Not Use Global >> Variables Like This". >> >> R has lists (and environments) that make working with objects that are >> associated with each other much simpler and fits better with the >> functional programming style of R. >> > > Thanks, Greg! > > Yours is a very smart solution to the problem I posed. > > By the way, what I'm trying to do is reading from a file a set of user given > parameters, in two paired columns: parameter-name, value; and then, managing > these parameters inside my R program. Now I do understand a bit more about > lists. What about environments? Are they similar to lists, and when, and how > are they created? > > Best regrards, > > -Sergio. > > __ > 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. -- Gregory (Greg) L. Snow Ph.D. 538...@gmail.com __ 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] hist() : is there a way to change the border width?
Thanks a lot. It works! On Thursday, December 19, 2013 10:53 PM, Jim Lemon wrote: On 12/20/2013 08:19 AM, capricy gao wrote: > I have played around with it and found that the only color could be changed. > But I really would like to change the width... > Hi Capricy, Try this on the first example for hist: hist(islands) par(lwd=3) hist(islands) Jim [[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] A function which is a sum of other functions...
My first thought was to use Reduce, but I think for this case that is a bit of overkill. You can have a vector or list of functions and just use sapply/lapply on the list of functions then sum the result. A quick example: > funs <- c(sin,cos,tan) > sapply( funs, function(f) f(pi/6) ) [1] 0.500 0.8660254 0.5773503 > sum(sapply( funs, function(f) f(pi/6) )) [1] 1.943376 Just wrap the above in a function with whatever options you want to use. If you need the functions to return vectors (of the same length) then you can still use sapply, but use rowSums, colSums, or apply on the result instead of sum. On Thu, Dec 19, 2013 at 12:05 PM, Onur Uncu wrote: > > Dear R Users > > I have a list of functions. Each function in the list is a function of single > variable. I would like to create a function (of one variable) which > represents the sum of all the functions in the list. So, if the functions in > my list are f1(x),..,f5(x) then I would like a new function > f(x)=f1(x)+f2(x)+...f5(x) > > Appreciate any suggestions on how to do this. > > I need the above f(x) function because I would like to minimise it with > respect to x using the nlm function. > > Thanks. > __ > 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. -- Gregory (Greg) L. Snow Ph.D. 538...@gmail.com __ 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] Inconsistent computation of an integral
Thanks William.
I was convinced by pmax, until I played with the following example today. I
tried to reproduce a computation from a paper. Here is the code:
P<- function(x) {
ab<-100*exp((0.0435-0.0269-0.5*0.3315^2)*4.3122+x*sqrt(4.3122)*0.3315)
return(ab)
}
integrand<- function(x){
cd<-
-1/2*((0.094+1.1465)*exp(0.0435*4.3122)+0.29/100*exp(0.0269*4.3122)*P(x)+0.89/100*max(P(x)-88.254,0))^(-2)*pnorm(x)
return(cd)
}
Solution<- integrate(integrand, lower=-Inf, upper=Inf)
Solution
The above code gives: -0.1366377
The paper reports: -0.1349
If I use pmax instead of max, I get: -0.1965606, which is much worse.
It may look like small discrepancies, but it makes a big difference to me.
I am still very puzzled by these discrepancies.
2013/12/19 William Dunlap
> I think you want to use pmax(x-50, 0), which returns a vector
> the length of x, instead of max(x-50,0), which returns a scalar.
>
> Bill Dunlap
> Spotfire, TIBCO Software
> wdunlap tibco.com
>
>
> > -Original Message-
> > From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org]
> On Behalf
> > Of Aurélien Philippot
> > Sent: Thursday, December 19, 2013 10:38 AM
> > To: R-help@r-project.org
> > Subject: [R] Inconsistent computation of an integral
> >
> > Dear R experts,
> > I computed the same integral in two different ways, and find different
> > values in R.
> > The difference is due to the max function that is part of the integrand.
> In
> > the first case, I keep it as such, in the second case, I split it in two
> > depending on the values of the variable of integration.
> >
> > 1) First computation
> >
> > # Function g
> > g<-
> > function(x){1/(x*0.20*sqrt(10)*sqrt(2*pi))*exp(-0.5*((log(x/50)-
> > 0.1*10)/(0.20*sqrt(10)))^2)}
> >
> > ### Function f1
> > f1<- function(x) {1/(500+10*x+1*max(x-50,0))}
> >
> > integrand1<- function(x) {
> > out<- f1(x)*g(x)
> > return(out)
> > }
> >
> > i2<- integrate(integrand1, lower=0, upper=Inf )$value
> >
> > It gives me: i2= 3.819418e-08
> >
> > 2) Second computation
> > I break the max function in two, depending on the values of the variable
> of
> > integration x (and I use the same density g as before):
> >
> > f11<- function(x) {1/(500+10*x)}
> > f12<- function(x) {1/(500+10*x+1*(x-50))}
> >
> >
> > integrand11<- function(x) {
> > out<- f11(x)*g(x)
> > return(out)
> > }
> >
> > integrand12<- function(x) {
> > out<- f12(x)*g(x)
> > return(out)
> > }
> >
> >
> > i21<- integrate(integrand11, lower=0, upper=50 )$value
> > +integrate(integrand12, lower=50, upper=Inf)$value
> >
> > I get i21=5.239735e-08
> >
> > The difference makes a huge difference for the computations I do. Does
> > anyone know where it comes from?
> > Thanks in advance!
> >
> > [[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.
>
[[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] Inconsistent computation of an integral
You should try to match up values of integrand() instead of the value of its
integral
on a single integral if you want to see if the function is being computed
correctly.
integrate(f) calls f(x) where x is a vector of values (typically 21 values).
It expects
that f() is vectorized: that f(x[i]) == f(x)[i] for i in seq_along(x). When
you put max(x)
into your function it messes that up:
> integrand(0.5)
[1] -0.07377961
> integrand(1.0)
[1] -0.05514443
> integrand(c(0.5, 1.0)) # expect the above two numbers in output vector
[1] -0.05106583 -0.05514443
Bill Dunlap
Spotfire, TIBCO Software
wdunlap tibco.com
From: Aurélien Philippot [mailto:aurelien.philip...@gmail.com]
Sent: Friday, December 20, 2013 8:59 AM
To: William Dunlap; R-help@r-project.org
Subject: Re: [R] Inconsistent computation of an integral
Thanks William.
I was convinced by pmax, until I played with the following example today. I
tried to reproduce a computation from a paper. Here is the code:
P<- function(x) {
ab<-100*exp((0.0435-0.0269-0.5*0.3315^2)*4.3122+x*sqrt(4.3122)*0.3315)
return(ab)
}
integrand<- function(x){
cd<-
-1/2*((0.094+1.1465)*exp(0.0435*4.3122)+0.29/100*exp(0.0269*4.3122)*P(x)+0.89/100*max(P(x)-88.254,0))^(-2)*pnorm(x)
return(cd)
}
Solution<- integrate(integrand, lower=-Inf, upper=Inf)
Solution
The above code gives: -0.1366377
The paper reports: -0.1349
If I use pmax instead of max, I get: -0.1965606, which is much worse.
It may look like small discrepancies, but it makes a big difference to me. I am
still very puzzled by these discrepancies.
2013/12/19 William Dunlap mailto:wdun...@tibco.com>>
I think you want to use pmax(x-50, 0), which returns a vector
the length of x, instead of max(x-50,0), which returns a scalar.
Bill Dunlap
Spotfire, TIBCO Software
wdunlap tibco.com
> -Original Message-
> From: r-help-boun...@r-project.org
> [mailto:r-help-boun...@r-project.org] On
> Behalf
> Of Aurélien Philippot
> Sent: Thursday, December 19, 2013 10:38 AM
> To: R-help@r-project.org
> Subject: [R] Inconsistent computation of an integral
>
> Dear R experts,
> I computed the same integral in two different ways, and find different
> values in R.
> The difference is due to the max function that is part of the integrand. In
> the first case, I keep it as such, in the second case, I split it in two
> depending on the values of the variable of integration.
>
> 1) First computation
>
> # Function g
> g<-
> function(x){1/(x*0.20*sqrt(10)*sqrt(2*pi))*exp(-0.5*((log(x/50)-
> 0.1*10)/(0.20*sqrt(10)))^2)}
>
> ### Function f1
> f1<- function(x) {1/(500+10*x+1*max(x-50,0))}
>
> integrand1<- function(x) {
> out<- f1(x)*g(x)
> return(out)
> }
>
> i2<- integrate(integrand1, lower=0, upper=Inf )$value
>
> It gives me: i2= 3.819418e-08
>
> 2) Second computation
> I break the max function in two, depending on the values of the variable of
> integration x (and I use the same density g as before):
>
> f11<- function(x) {1/(500+10*x)}
> f12<- function(x) {1/(500+10*x+1*(x-50))}
>
>
> integrand11<- function(x) {
> out<- f11(x)*g(x)
> return(out)
> }
>
> integrand12<- function(x) {
> out<- f12(x)*g(x)
> return(out)
> }
>
>
> i21<- integrate(integrand11, lower=0, upper=50 )$value
> +integrate(integrand12, lower=50, upper=Inf)$value
>
> I get i21=5.239735e-08
>
> The difference makes a huge difference for the computations I do. Does
> anyone know where it comes from?
> Thanks in advance!
>
> [[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.
[[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] A function which is a sum of other functions...
Makes sense. Simple and works. Thank you Greg. Sent from my iPhone > On 20 Dec 2013, at 16:15, Greg Snow <538...@gmail.com> wrote: > > My first thought was to use Reduce, but I think for this case that is > a bit of overkill. You can have a vector or list of functions and > just use sapply/lapply on the list of functions then sum the result. > A quick example: > >> funs <- c(sin,cos,tan) >> sapply( funs, function(f) f(pi/6) ) > [1] 0.500 0.8660254 0.5773503 >> sum(sapply( funs, function(f) f(pi/6) )) > [1] 1.943376 > > Just wrap the above in a function with whatever options you want to > use. If you need the functions to return vectors (of the same length) > then you can still use sapply, but use rowSums, colSums, or apply on > the result instead of sum. > >> On Thu, Dec 19, 2013 at 12:05 PM, Onur Uncu wrote: >> >> Dear R Users >> >> I have a list of functions. Each function in the list is a function of >> single variable. I would like to create a function (of one variable) which >> represents the sum of all the functions in the list. So, if the functions in >> my list are f1(x),..,f5(x) then I would like a new function >> f(x)=f1(x)+f2(x)+...f5(x) >> >> Appreciate any suggestions on how to do this. >> >> I need the above f(x) function because I would like to minimise it with >> respect to x using the nlm function. >> >> Thanks. >> __ >> 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. > > > > -- > Gregory (Greg) L. Snow Ph.D. > 538...@gmail.com __ 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] by class...
I used the by() function on a data.frame to get sums of the data grouped by 2 factors. The function worked however the output is in a class called 'by'. Not familiar with this class. How can I turn the output into a nice table where columns represent values of factor1, row represent values of factor2 and the entries in the table are the sums that were calculated using the by function? I did some web search which suggested using do.call(rbind, datframe_object) but this command gave the following error: "Second argument must be a list"... Thank you. __ 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] ftable and data.frame
On Dec 20, 2013, at 5:01 AM, wrote: > Hi, > > I used this command to produce a table: It's not actually a 'table'. > > (tab1 = ftable(SEX, ESTCIV, Q1)) is.table(tab1) # will return FALSE > > Q1 B L M N > SEXOESTCIV >F A 11 13 4 2 > E 1 0 0 0 >M A 5 0 3 1 > E 0 0 0 0 > > but I need something like: > > > SEXOESTCIVB L M N > FA 11 13 4 2 > FE 1 0 0 0 > M A 5 0 3 1 > M E 0 0 0 0 > > How can I get it? You need to describe the purpose of this effort. If is for display, the answer will be to look at the code for `print.ftable` (which eventually passes its arguments to `format.ftable`. If it is for creating something other than an ftable (such a data.frame or character-matrix), then you will probably need to coerce the ftable to a matrix to get the entry values and then extract the row labels from the 'tab1'-object with: attr(tab1, "row.vars") and `rep` them approriately. Or you could work with the output from format(ftable(tab1))[ ,1:2]. cbind( format(ftable(Titanic), quote=FALSE)[ ,1:3], format(ftable(Titanic), quote=FALSE)[ ,5:6]) [,1][,2] [,3][,4] [,5] [1,] " " " " " " " No" "Yes" [2,] "Class" "Sex " "Age " " " " " [3,] "1st " "Male " "Child" " 0" " 5" [4,] " " " " "Adult" "118" " 57" [5,] " " "Female" "Child" " 0" " 1" [6,] " " " " "Adult" " 4" "140" [7,] "2nd " "Male " "Child" " 0" " 11" [8,] " " " " "Adult" "154" " 14" [9,] " " "Female" "Child" " 0" " 13" [10,] " " " " "Adult" " 13" " 80" [11,] "3rd " "Male " "Child" " 35" " 13" [12,] " " " " "Adult" "387" " 75" [13,] " " "Female" "Child" " 17" " 14" [14,] " " " " "Adult" " 89" " 76" [15,] "Crew " "Male " "Child" " 0" " 0" [16,] " " " " "Adult" "670" "192" [17,] " " "Female" "Child" " 0" " 0" [18,] " " " " "Adult" " 3" " 20" > > I need this format to use ordinal logistic regression and I have many tables. That makes me think using ftable is the completely misguided approach. You should describe in more detail you plans to use ordinal logistic regression. I'm also wondering if you have used attach(). It would make more sense to keep these variables in a dataframe. Almost all regression functions expect input in the form of a normalized dataframe rather than in separate variables of aggregated counts. > (tab1 = ftable(SEX, ESTCIV, Q1)) > > Thanks, > > Silvano. > > --- > Este email está limpo de vírus e malwares porque a proteção do avast! > Antivírus está ativa. > > > [[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. David Winsemius Alameda, CA, USA __ 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] by class...
On Dec 20, 2013, at 9:38 AM, Onur Uncu wrote: > I used the by() function on a data.frame to get sums of the data grouped by 2 > factors. The function worked however the output is in a class called 'by'. > Not familiar with this class. How can I turn the output into a nice table > where columns represent values of factor1, row represent values of factor2 > and the entries in the table are the sums that were calculated using the by > function? > > I did some web search which suggested using do.call(rbind, datframe_object) > but this command gave the following error: > "Second argument must be a list"... The output of by() is always a list, so I would have expected: do.call(rbind, by.object) to have retruned a different error message that what you suggest. > > > Thank you. > __ > 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. David Winsemius Alameda, CA, USA __ 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] The difference between SAS and R
SAS uses words. R uses symbols. See this: http://xkcd.com/1306/ (Yes, I know IML uses plenty of symbols. It's just supposed to be funny. And somewhat true.) -- Kevin Wright [[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] by class...
Hi, Try: as.table(by(warpbreaks[,1],warpbreaks[,-1],sum)) #or to convert to data.frame as.data.frame(as.table(by(warpbreaks[,1],warpbreaks[,-1],mean))) A.K. On Friday, December 20, 2013 12:39 PM, Onur Uncu wrote: I used the by() function on a data.frame to get sums of the data grouped by 2 factors. The function worked however the output is in a class called 'by'. Not familiar with this class. How can I turn the output into a nice table where columns represent values of factor1, row represent values of factor2 and the entries in the table are the sums that were calculated using the by function? I did some web search which suggested using do.call(rbind, datframe_object) but this command gave the following error: "Second argument must be a list"... Thank you. __ 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.
Re: [R] by class...
Hi, You can also try: library(reshape2) dcast(as.data.frame(as.table(by(warpbreaks[,1],warpbreaks[,-1],sum))),wool~tension, value.var="Freq") A.K. On , arun wrote: Hi, Try: as.table(by(warpbreaks[,1],warpbreaks[,-1],sum)) #or to convert to data.frame as.data.frame(as.table(by(warpbreaks[,1],warpbreaks[,-1],mean))) A.K. On Friday, December 20, 2013 12:39 PM, Onur Uncu wrote: I used the by() function on a data.frame to get sums of the data grouped by 2 factors. The function worked however the output is in a class called 'by'. Not familiar with this class. How can I turn the output into a nice table where columns represent values of factor1, row represent values of factor2 and the entries in the table are the sums that were calculated using the by function? I did some web search which suggested using do.call(rbind, datframe_object) but this command gave the following error: "Second argument must be a list"... Thank you. __ 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.
Re: [R] by class...
On Dec 20, 2013, at 9:45 AM, David Winsemius wrote: > > On Dec 20, 2013, at 9:38 AM, Onur Uncu wrote: > >> I used the by() function on a data.frame to get sums of the data grouped by >> 2 factors. The function worked however the output is in a class called 'by'. >> Not familiar with this class. How can I turn the output into a nice table >> where columns represent values of factor1, row represent values of factor2 >> and the entries in the table are the sums that were calculated using the by >> function? >> >> I did some web search which suggested using do.call(rbind, datframe_object) >> but this command gave the following error: >> "Second argument must be a list"... > > The output of by() is always a list, so I would have expected: Or maybe my memories are false. ?by "This is always a list if simplify is false, otherwise a list or array (see tapply)." So maybe you need simplify=FALSE > > do.call(rbind, by.object) to have retruned a different error message that > what you suggest. > > >> >> >> Thank you. >> __ >> 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. > > David Winsemius > Alameda, CA, USA > > __ > 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. David Winsemius Alameda, CA, USA __ 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] t test for multiple classes
dear all,
I have a dataset composed by different classes: I have two main classes,
"active (Act.)" and "latent (Lat.)", furher subdivided in the gene
subclasses "IP10" and "MIG" and then into 7 stimulus sub-subclasses "ESAT6"
to "PHA".
Would it be possible to test in a direct way the statistical differences
between the response variable ("ratio") of the classes active and latent
values for each of the sub-subclasses and stratified by the subclass gene?
that is: how can i find differences for, let's say, the ESAT6 ratios
between active and latent for the gene IP10, then for the ESAT6 active
versus latent for the MIG gene; then reiterate for the CFP10 for IP10,
CFP10 for MIG and so on and so forth.
the data looks like:
TB gene stimulus ratio
Act. IP10 ESAT6 0.752020287
Act. IP10 CFP10 0.893628674
Act. IP10 Rv3615c 0.958108603
Act. IP10 Rv2654 0.319492309
Act. IP10 Rv3879 0.074381059
Act. IP10 Rv3873 0.874629325
Act. IP10 PHA 0.308093623
Act. MIG ESAT6 0.398817499
Act. MIG CFP10 0.022455057
Act. MIG Rv3615c 0.186013561
Act. MIG Rv2654 0.040591894
Act. MIG Rv3879 0.716143409
Act. MIG Rv3873 0.271241311
Act. MIG PHA 0.154207839
Lat. IP10 ESAT6 0.579278956
Lat. IP10 CFP10 0.553393777
Lat. IP10 Rv3615c 0.918602425
Lat. IP10 Rv2654 0.56821826
Lat. IP10 Rv3879 0.429991707
Lat. IP10 Rv3873 0.938468864
Lat. IP10 PHA 0.621095162
Lat. MIG ESAT6 0.486344894
Lat. MIG CFP10 0.288455361
Lat. MIG Rv3615c 0.602348528
Lat. MIG Rv2654 0.717395188
Lat. MIG Rv3879 0.178745118
Lat. MIG Rv3873 0.330885875
Lat. MIG PHA 0.327117342
Many thanks for any possible help.
Best regards,
Luigi
[[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] by class...
You seem to be falling prey to a common misconception that "R" is some monolithic tool, when in fact it is a herd of cats. The "by" function, from the "base" package, returns a list of results returned by your function. One approach to making a data frame out of that is to use the simplify2array function, transpose, and convert to data frame. However, there are many other ways to aggregate data as well. Packages plyr, sqldf, data.table all have strengths and weaknesses, but you must always keep in mind that they are contributed packages and you have to refer to the package name if you ask here about using functions from them. --- Jeff NewmillerThe . . Go Live... DCN:Basics: ##.#. ##.#. Live Go... Live: OO#.. Dead: OO#.. Playing Research Engineer (Solar/BatteriesO.O#. #.O#. with /Software/Embedded Controllers) .OO#. .OO#. rocks...1k --- Sent from my phone. Please excuse my brevity. Onur Uncu wrote: >I used the by() function on a data.frame to get sums of the data >grouped by 2 factors. The function worked however the output is in a >class called 'by'. Not familiar with this class. How can I turn the >output into a nice table where columns represent values of factor1, row >represent values of factor2 and the entries in the table are the sums >that were calculated using the by function? > >I did some web search which suggested using do.call(rbind, >datframe_object) but this command gave the following error: >"Second argument must be a list"... > > >Thank you. >__ >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.
Re: [R] Strange subvector output --> x[n] != x[1:n][n]
You've found an interesting corner of Circle 1 of 'The R Inferno'. http://www.burns-stat.com/documents/books/the-r-inferno/ The issue is that your 'n' in the final case is slightly less than 11. So: > n [1] 11 > as.integer(n) [1] 10 > 1:n [1] 1 2 3 4 5 6 7 8 9 10 11 The mystery to me is why `:` thinks it is doing an integer sequence but ends in 11 rather than 10. Most people are likely to think the mystery is why as.integer(n) is 10. The reason is that coercion to integer is truncation (except if the number is really close to the integer farther from 0). Why that and not round? Well, just because. (Actually probably speed back in the day when it could matter.) Pat On 20/12/2013 04:33, Gewart wrote: Hi, Can anyone explain what is going on...!? For a vector "x=seq(min,max,0.01)", when generating sub-vector "a" based on a starting value "st", things go as expected as long as "st" is not too close to the beginning of "x". For example, if x starts at -5 and increments by 0.01, whenever I try to generate the sub-vector "a" (as below) with a starting value of 0.49 or less it does not generate the expected output: The initial value of "a" is wrong. Thanks in advance for any clarity you can shed. Gary ...(please see two versions of code below) #THIS WORKS...(st > -4.9) min = -5; max = 1; x=seq(min,max,0.01) st= -4.8 ; end= 0 a=x[((st-min)/0.01+1):((end-min)/0.01+1)] n=(st-min)/0.01+1 #compare a[1:10]; c(x[n:(n+9)]) #test... n x[1:n]; x[n] ### x[n]== x[1:n][n] ; As expected ## # BUT THIS IS WEIRD!!...(st <= -4.9) st= -4.90 ; end= 0 ### -> BUG in generation of a!! a=x[((st-min)/0.01+1):((end-min)/0.01+1)]; n=(st-min)/0.01+1 #compare a[1:10]; c(x[n:(n+9)]) #test n x[1:n]; x[n] ### NOW x[n] != x[1:n][n] !!?? What is going on!? -- View this message in context: http://r.789695.n4.nabble.com/Strange-subvector-output-x-n-x-1-n-n-tp4682526.html Sent from the R help mailing list archive at Nabble.com. __ 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. -- Patrick Burns pbu...@pburns.seanet.com twitter: @burnsstat @portfolioprobe http://www.portfolioprobe.com/blog http://www.burns-stat.com (home of: 'Impatient R' 'The R Inferno' 'Tao Te Programming') __ 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] Grab multinomial coefficients
Hello! I am using function multinom: library(nnet) library(MASS) mnl<-multinom(myDV~. ,data=mydata,na.action="na.omit", MaxNWts = 2000, maxit = 1000) I can see the resulting coefficients: print(mnl) But how could I grab them? I am not finding them when I do: str(mnl) Thank you! -- Dimitri Liakhovitski [[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] Grab multinomial coefficients
Hello, Try coef(mnl) Hope this helps, Rui Barradas Em 20-12-2013 20:21, Dimitri Liakhovitski escreveu: Hello! I am using function multinom: library(nnet) library(MASS) mnl<-multinom(myDV~. ,data=mydata,na.action="na.omit", MaxNWts = 2000, maxit = 1000) I can see the resulting coefficients: print(mnl) But how could I grab them? I am not finding them when I do: str(mnl) Thank you! __ 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] Grab multinomial coefficients
That's it! Thanks a lot! On Fri, Dec 20, 2013 at 4:18 PM, Rui Barradas wrote: > Hello, > > Try > > coef(mnl) > > Hope this helps, > > Rui Barradas > > Em 20-12-2013 20:21, Dimitri Liakhovitski escreveu: > > Hello! >> >> I am using function multinom: >> >> library(nnet) >> library(MASS) >> >> mnl<-multinom(myDV~. ,data=mydata,na.action="na.omit", MaxNWts = 2000, >> maxit = 1000) >> >> >> I can see the resulting coefficients: >> >> print(mnl) >> >> But how could I grab them? I am not finding them when I do: >> str(mnl) >> >> Thank you! >> >> > -- Dimitri Liakhovitski [[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.
