[R] how do I get the legend?

http://en.wikipedia.org/wiki/Binomial_distribution Bino Hi there, How can I get the legend in the probability density plot in that http://en.wikipedia.org/wiki/File:Binomial_distribution_pmf.svg wiki page ? and the cumalative plot as well? # Binomial Dist

[R] test the goodness of it for negative binomial type 2

[code]library(MASS) x=c(rep(0,8096), rep(1,1629), rep(2,233), rep(3,38), rep(4,4) ) x.bar=round(mean(x),4) x.var=round(var(x),4) p.hat=round(x.bar/x.var,4) alpha.hat=round(x.bar*p.hat/(1-p.hat),4) fitdistr(x, "Negative Binomial") fitdistr(x, "Poisson")[/c

[R] TukeyHSD model thing

Hi, I am trying to reproduce a tukey test in R == x=c(145,40,40,120,180, 140,155,90,160,95, 195,150,205,110,160, 45,40,195,65,145, 195,230,115,235,225, 120,55,50,80,45 ) y2=c( rep(as.character(1),5), rep(as.c

[R] Fisher's LSD and tukey output thing

2-and how can I do that table in the attachement (page 8/9)? ( i dont know how to describe that thing, orginally produced in SAS) Q3-in the tukey output, can i ask to to show whether it's significant or not by indicating ***? (attachment page 6) Thanks! casperyc SAS output http://n4.nabbl

[R] barplot with factors problem

http://www.harding.edu/fmccown/R/#autosdatafile http://www.harding.edu/fmccown/R/#autosdatafile I am tring to get a barchat by factors, following the example in that link above. === x=c(145,40,40,120,180, 140,155,90,160,95, 195,150,205,110,160, 45

Re: [R] barplot with factors problem

http://n4.nabble.com/file/n1583733/100307070476876317b486a941.jpg I want to get a histogram by factors. Thanks. -- View this message in context: http://n4.nabble.com/barplot-with-factors-problem-tp1583671p1583733.html Sent from the R help mailing list archive at Nabble.com. _

[R] variance of discrete uniform distribution

Hi all, I am REALLY confused with the variance right now. for a discrete uniform distribution on [1,12] the mean is (1+12)/2=6.5 which is ok. y=1:12 mean(y) then var(y) gives me 13 1- on http://en.wikipedia.org/wiki/Uniform_distribution_%28discrete%29 wiki the variance is (12^2-1)/12=

Re: [R] variance of discrete uniform distribution

Hi Rolf Turner , God, it directed to the wrong page. I firstly find the formula in wiki, than tried to verify the answer in R, now, given that 143/12 ((n^2-1)/12 ) is the correct answer for a discrete uniform random variable, I am still not sure what R is calculating there? why it gives me 13?

Re: [R] test the goodness of it for negative binomial type 2

Hi Achim Zeileis-4, That's very helpful. Thanks! -- View this message in context: http://n4.nabble.com/test-the-goodness-of-it-for-negative-binomial-type-2-tp1575892p1585357.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-

Re: [R] curvedarrow (some graphics problem)

Hi there, This is now the code % library(grid) vp <- viewport( x = unit(0, "npc"), y = unit(0, "npc"), just = c("left", "bottom"), xscale = c(-1, 1) , yscale = c(-1, 1)) vp2 <- viewport( # probably not needed but I had trouble placing the xaxis x = unit(0,"npc"), y = unit(0.5, "npc"

[R] R loop help

Hi, I am tring to write a loop to compute this, == x1=c( rep(-1,4), rep(1,4) ) x2=c( rep(c(-1,-1,1,1),2) ) x3=c( rep(c(-1,1),4) ) x1*x2 x1*x3 x2*x3 suppose i have x1,x2,x3

[R] model reparameterization

== y=c(100,200,300,400,500) treatment=c(1,2,3,3,4) block=c(1,1,2,3,3) summary(lm(y~as.factor(treatment)+as.factor(block))) == The aim is to find a model that can estimate the comparison between treatment 1 with 2 and treatment 3 wit

Re: [R] R loop help

Hi there, That's exactly what I want. I have checked ?combn out, but I could get the following, suppose that I want ALL possible combinations of them, as this == apply( combn(paste('x', 1:4, sep =""), 2), 2, function(v) get(v[1])*get(v[2])

[R] writing function help

Hi, I have written a function, ( or 2 functions) becasue there are only tiny difference, One of them has lines == #c2 (price.mat) call option price c2=c( max(S.mat[1,3]-K,0), max(S.mat[2,3]-K,0), max(S.mat[3,3]-K,0) ) (some lines) list( "Stock Value"=round(S.mat,dp), "C

[R] writing function (plot problem）

Hi, === x=rnorm(20) y=rnorm(20) t=lm(y~x) plot(t) === you will get "click or hit enter to next page..." how do I write a function to archieve this ? say plot(x,y) then pause, wait plot(y,x) Thanks! casper -- View this message in context: http://n4

[R] writing function ( 'plot' and 'if') problem

=== myf=function(ds=1){ x=rnorm(10) y=rnorm(10) { #start of if if (ds==1) { list(x,y) } else (ds==2) { plot(x,y) } } # end of if } # end of function === Hi All, the problem i am having here is, that I want to be able to control the display, lf

Re: [R] writing function ( 'plot' and 'if') problem

what i am triyng to do is when ds=1 give me a list ds=2 plot Thanks casper -- View this message in context: http://n4.nabble.com/writing-function-plot-and-if-problem-tp1839091p1839125.html Sent from the R help mailing list archive at Nabble.com. _

[R] histogram breaks

=== Q2=c( + sample(10:19,8,T), + sample(20:24,15,T), + sample(25:29,25,T), + sample(30:39,18,T), + sample(40:49,12,T), + sample(50:64,7,T), + sample(65:89,5,T) + ) hist(Q2) can give me a histogram, however, how do i get a different 'b

[R] How to save as PDF when I used par(ask=T)

== y=c(9,9,17,11,7,8,15,5) treat=c("A","C","B","C","D","A","B","D") block=c(1,1,2,2,2,3,3,3) q1.mod=aov(y~as.factor(treat)+as.factor(block)) q1.mod summary(q1.mod) par(ask=T) plot(TukeyHSD(q1.mod))

[R] bootstrap help

Hi, I have 49 pairs in my data.frame 'data' x y 76 80 138 143 67 67 29 50 381 464 23 48 37 63 120 115 ... ... how do I get a bootstrap sample of size n=50? i have tried sample(data,size=50,replace=TRUE) and sample(data,replace=TRUE) both didnt see

Re: [R] bootstrap help

mated by the bar{x}/bar{y} ) i mis undertood that too. but now, i think what i actually need is 50 sets of this kind (49 rows) of tables from bootstrapping so that i can have 50 different 'bar{x}/bar{y}' s then to estimate the mean... Thanks. casper Ben Bolker wrote: > > ca

[R] confint for glm (general linear model)

Hi, I have a glm gives summary as follows, Estimate Std. Errorz valuePr(>|z|) (Intercept) -2.03693352 1.449574526 -1.405194 0.159963578 A0.01093048 0.006446256 1.695633 0.089955471 N0.41060119 0.224860819 1.826024 0.

Re: [R] bootstrap help

ize=nrow(data),replace=TRUE),])})) > > > r-help@r-project.org > > 2009/12/7 casperyc : >> >> Hi there, >> >> i think that's not what i was aiming for... >> >> i was aked to >> >> Generate 50 Bootstrap samples and corresponding estimates >>

Re: [R] confint for glm (general linear model)

does no one know this? -- View this message in context: http://n4.nabble.com/confint-for-glm-general-linear-model-tp954071p956658.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/m

Re: [R] confint for glm (general linear model)

ult # t-distribution > > Walmes Zeviani - Brazil > > > > casperyc wrote: >> >> Hi, >> >> I have a glm gives summary as follows, >> >>Estimate Std. Errorz valuePr(>|z|) >> (Intercept) -2.03693352

Re: [R] confint for glm (general linear model)

for an example, counts <- c(18,17,15,20,10,20,25,13,12) outcome <- gl(3,1,9); treatment <- gl(3,3) glm.D93 <- glm(counts ~ outcome + treatment, family=poisson()) confint(glm.D93) confint.default(glm.D93) # based on asymptotic normality to verify the confidence interval (confint.default(glm.D93

[R] curvedarrow (some graphics problem)

Hi there, I just wonder how to draw this kind of picture... http://www.nabble.com/file/p24158796/b.jpg http://www.nabble.com/file/p24158796/a.jpg and this is what i have done % library(shape) library(diagram) curve(sin(x),bty="n",-8,8,yaxt="n",ylab="",xaxt="n",type="n",xlab="") axis(1,la

Re: [R] curvedarrow (some graphics problem)

baptiste auguie-4 wrote: > > Hi, > > Grid offers several functions to help drawing such graphs, > > see Paul Murrell's "Can R Draw Graphs? " (useR! 2006) > > I came up with this, as a quick example, > > vp <- viewport( > x = unit(0, "npc"), > y = unit(0, "npc"), > just = c("left", "bottom

Re: [R] curvedarrow (some graphics problem)

[quote]> grid.text("t", x=0.3, y=unit(-1, "line"), vp=vp2) Error in valid.units(units) : Invalid unit > grid.text("s", 0.8, unit(-1, "line"), gp=gpar(col="red"), vp=vp2) Error in valid.units(units) : Invalid unit[/quote] Hi, I have loaded the [grid] package. It seems that I have to load anoth

[R] legend help

i=1 for(rate in c(2,4) ){ for(shape in c(1,3,5) ){ curve(dgamma(x,rate,shape),xlim=c(0,3),ylab="",col=i,lty=i,add=T) i=i+1 } } How can I add some legend to represent these lines? i.e. the legend is displayed as col=1 lty=1 lambda=2 theta=1 col=2 lty=2 lambda=2 the

Re: [R] legend help

Yes, that is pretty much what I want. However, there was slightly a mistake. we need to use ''rate=rate[i]"" and "shape=shape[i]" because the default is == dgamma(x, shape, rate = 1, scale = 1/rate, log = FALSE) ==

Re: [R] legend help

This is very weird!! I know the 'lty=1' command and I tried yesterday , many times!! didnt work. However, today, it worked!!! Maybe the university's computer is stupid! Thanks! -- View this message in context: http://n4.nabble.com/legend-help-tp1461992p1470216.html Sent from the R help mailin

[R] dot plot by group

Hi all, I have the folloing data table %% TypeBATCH RESPONSE SHORT A 22 SHORT A 3 SHORT A 16 SHORT A 14 SHORT A 8 SHORT A 27 SHORT A 11 SHORT A 17 SHORT B 12 SHORT B 17 SHORT B 11 SHORT

Re: [R] dot plot by group

Hi Spector, Yes, that is exactly what I was aiming for. Thanks. Casper -- View this message in context: http://r.789695.n4.nabble.com/dot-plot-by-group-tp2990469p2990495.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-proj

Re: [R] dot plot by group

And now I just wonder why the ' bty='n' ' won't work? I did dotplot(BATCH~RESPONSE,data=d,subset=Type=='SHORT',bty='n') and tried other bty parameters, none is working Casper -- View this message in context: http://r.789695.n4.nabble.com/dot-plot-by-group-tp2990469p2990500.html Sent from t

Re: [R] dot plot by group

Thanks. It was 'bty=n' on a text graph, I was just having fun with it. It does not really matter that much. CasperYC -- View this message in context: http://r.789695.n4.nabble.com/dot-plot-by-group-tp2990469p3004323.html Sent from the R help mailing list archive at

[R] how do i plot this "hist"?

Hi all, I have the following data in abc.dat === 50 0 1 0 0 55 114 0 1 60 786 0 3 6522 324 2 3 7058 1035 1 7 7530 2568 034 80 9 293615 162 8527 216946 365 90

Re: [R] how do i plot this "hist"?

Hi all, Thank you both. The codes work perfectly. I now use barplot(t(x.m)[-1,], names.arg = x.m[,1]) to make the 'histogram' plot. Now how do I add a dentity line to it? like 'superpose' on the graph? say if I got a Normal with mean 100 and variance 3... And I can work out the mean and vari

Re: [R] how do i plot this "hist"?

Thanks！ now it's just a matter of 'superposing' the density onto the barplot I am looking at the example here, but it does not seem to applicable to barplot. casper -- View this message in context: http://r.789695.n4.nabble.com/how-do-i-plot-this-hist-tp3032796p3035110.html Sent from the R h

[R] par mfrow in "function" problem

Hi all, I defined the following # myhist=function(x){ hist(x,xlab="",main="") h=hist(x) xfit=seq(min(x),max(x),length=100) yfit=dnorm(xfit,mean(x),sd=sd(x)) yfit=yfit*diff(h$mids[1:2])*length(x) lines(xfit, yfit, col="bl

[R] what's wrong with this 'length' in function?

Hi all, I am having a trouble with this function I wrote ### p26=function(x,alpha){ # dummy variable j=1 ci=matrix(ncol=2,nrow=3) while (j<4){ if (j==2) {x=x+c(-1,1)*0.5}

[R] plot linear model problem

Hi all, Say I fit a linear model, and saved it as 'test.lm' Then if I use plot(test.lm) it gives me 4 graphs How do I ask for a 'subset' of it?? say just want the 1st graph, the residual vs fitted values, or the 1,3,4th graph? I think I can use plot(test.lm[c(1,3,4)]) before, but now, it's

Re: [R] plot linear model problem

Thank you both. casper -- View this message in context: http://r.789695.n4.nabble.com/plot-linear-model-problem-tp3045763p3045932.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/

[R] how exactly does 'identify' work?

Hi all, # test=data.frame(x=1:26,y=-23.5+0.45*(1:26)+rnorm(26)) rownames(test)=LETTERS[1:26] attach(test) #test test.lm=lm(y~x) plot(test.lm,2) identify(test.lm$res,,row.names(test)) # not working plot(x,y) identify(x,y,row.names(test)) # works fine ident

Re: [R] how exactly does 'identify' work?

Hi, I think the problem is 1 - when a linear model is fitted, ploting the qqnorm( test.lm$ res ) we dont 'know' what values are actually being used on the y-axis; and how do we refer to the ‘Index’ on the x-axis?? therefore, i dont know how to refer to the x and y coordinates in the identi

Re: [R] how exactly does 'identify' work?

yes, i tried to modify the "2L" part in plot.lm ### if (show[2L]) { ylim <- range(rs, na.rm = TRUE) ylim[2L] <- ylim[2L] + diff(ylim) * 0.075 qq <- qqnorm(rs, main = main, ylab = ylab23, ylim = ylim, ...) if (qqline)

Re: [R] how exactly does 'identify' work?

Omg! yes, it is working now! tmp <- qqnorm( resid(test.lm) ) What a simple nice trick!!! Actually, i wasnt looking for the 'i'th label, I was looking for the 'row.names' as label, like I stated in the 1st post. > identify(tmp, , row.names(test) ) is the label i have been trying to get. T

[R] how to add frequencies to barplot

Hi, I have count data x2=rep(c(0:3),c(13,80,60,27)) x2 0 1 2 3 13 80 60 27 I want to graph to be ploted as barplot(table(x2),density=4) how do I add relative frequency to it, like in hist(x2,labels=T) above the 'bar's Thanks. casper -- View this message in context: http://r.789

[R] how to divide each column in a matrix by its colSums?

Hi, I have a matrix, say m=matrix(c( 983,679,134, 383,416,84, 2892,2625,570 ),nrow=3 ) i can find its row/col sum by rowSums(m) colSums(m) How do I divide each row/column by its rowSum/colSums and still return in the matrix form? (i.e. the new rowSums/colSums =

Re: [R] how to divide each column in a matrix by its colSums?

In that case, there are values >1, which is clearly not what I wanted. Thanks. I think I should use prop.table -- View this message in context: http://r.789695.n4.nabble.com/how-to-divide-each-column-in-a-matrix-by-its-colSums-tp3062739p3062752.html Sent from the R help mailing list archive

Re: [R] how to divide each column in a matrix by its colSums?

I am using > prop.table(m,1) and > prop.table(m,2) my aim. which I think is the most 'easy' way. Thanks. Casper -- View this message in context: http://r.789695.n4.nabble.com/how-to-divide-each-column-in-a-matrix-by-its-colSums-tp3062739p3062797.html Sent from the R help mailing list archive

[R] how to add these "axis" label？

Hi All, How do I add these axis labels? ### p=seq(0,1,length.out=500) p=p[-c(1,length(p))] g1=log(p/(1-p)) g2=qnorm(p) g3=log(-log(1-p)) g4=-log(-log(p)) plot(p,g1, 'n',ylim=c(-5,5),las=1, bty='n', xaxt='n',yaxt='n', xla

Re: [R] how to add these

Thomas Stewart wrote: > > It isn't clear to me what you want to do. Do you want the axes to show? > Do > you want labels for the lines? Do you want a legend? What is your > desired > output? > > -tgs > > On Wed, Dec 8, 2010 at 2:42 PM, casperyc wrote:

Re: [R] how to add these

thanks! problem solved! casper -- View this message in context: http://r.789695.n4.nabble.com/how-to-add-these-tp3078908p3079340.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/

[R] How to print colorful R output??

Hi All, I wonder if there is a way to print the R output with COLOR? Not the color plots, but the outputs in the console. Thank. casper -- View this message in context: http://r.789695.n4.nabble.com/How-to-print-colorful-R-output-tp3082750p3082750.html Sent from the R help mailing list archiv

Re: [R] How to print colorful R output??

Hi All, My aim is actually not that complicated as you guys understand. What I want is this, when I print it by clicking File--> Print... It gaves me a black white output. But what I want is 'red', for all the codes i typed in, 'blue', for the R output, just like the console. Thanks! (I

[R] use of 'apply' for 'hist'

Hi all, ## dof=c(1,2,4,8,16,32) Q5=matrix(rt(100,dof),100,6,T,dimnames=list(NULL,dof)) par(mfrow=c(2,6)) apply(Q5,2,hist) myf=function(x){ qqnorm(x);qqline(x) } apply(Q5,2,myf) ## These looks ok. However, I would lik

[R] how to use paste in function for expression?

Hi all, ## dof=c(1,2,4,8,16,32) Q5=matrix(rt(100,dof),100,6,T,dimnames=list(NULL,dof)) par(mfcol=c(2,6)) for (i in 1:6){ dof2=dof[i] hist(Q5[,i],main=paste("t[",dof2,"]",sep="")) qqnorm(Q5[,i]) qqline(Q5[,i]) } ###

Re: [R] use of 'apply' for 'hist'

Thanks Dieter. Casper -- View this message in context: http://r.789695.n4.nabble.com/use-of-apply-for-hist-tp3093811p3093938.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailm

Re: [R] use of 'apply' for 'hist'

HI Sarah, I will just use a loop then. I think my colnames are fine. Thanks! casper -- View this message in context: http://r.789695.n4.nabble.com/use-of-apply-for-hist-tp3093811p3093937.html Sent from the R help mailing list archive at Nabble.com. _

[R] how to use expression as function arguements?

Hi all, # integ=function(f,n){ # computes vector y = f(x) x=runif(1) y=f hatI=mean(y) hatI } # example of use integ(x^2+x+100,100) # it returns an error says no obj 'x' how do I '

[R] how to see what's wrong with a self written function?

Hi all, I am writing a simple function to implement regularfalsi (secant) method. ### regulafalsi=function(f,x0,x1){ x=c() x[1]=x1 i=1 while ( f(x[i])!=0 ) { i=i+1 if (i==2) {

Re: [R] how to see what's wrong with a self written function?

Hi Petr Savicky, It is very interesting and a good way to plot, so that helps me to visualized the method, also easier to see what's wrong. The "readline" is new to me, and I should learn how to use it in the future. Thanks! casper -- View this message in context: http://r.789695.n4.nabble.c

[R] Binary response GLM Question

Hi all, I have a problem with binary response data in GLM fitting. The problem is that the "y" take only 1 or 0, and if I use logit link, it is the log of the odds ratio, which is p/(1-p). In my situation, think "y" is "p", so sometimes the odds is 0, sometimes it is "1/0", which is (should be) un

Re: [R] Binary response GLM Question

Hi Josh, Thank you for your reply. The reason I thank "Y" (0 and 1) here as p is because I think each observation is just a bernulli trial, so in this case the binomial n=1. And yet R still fits it (with the logit link) . I know the expression for the logit link, so I assumed I can take y here as

[R] what controls the limits of x axis in hclust

Hi all, I followed the example in http://www.statmethods.net/advstats/cluster.html and apply it to one of my own dataset, I got this tiny problem with the boreders, http://r.789695.n4.nabble.com/file/n3340238/temp.png The red rectangle are 'outside' the plot, so I want to know how do I cont

[R] how to quote "factors" in a function?

### code ### x=sample(LETTERS[1:26],100,T) prob=function(y){ count=sum(x==y) total=length(x) count/total } ### end ### How do I quote the letters A,B,C,D ect, in order to make the "prob" function return the weights of the desired Letter? Thanks! --

[R] function input as variable name (deparse/quote/paste) ??

Hi all Say I have a function: myname=function(dat,x=5,y=6){ res<<-x+y-dat } for various input such as myname(dat1) myname(dat2) myname(dat3) myname(dat4) myname(dat5) how should I modify the 'res' line, to have new informative variable name correspondingly, such as dat1.res dat2.res dat3.

Re: [R] function input as variable name (deparse/quote/paste) ??

Sorry if I wasn't stating what I really wanted or it was a bit confusing. Basically, there are MANY datasets to run suing the same function I have written a function to analyze it and returns a LIST of useful out put in the variable 'res' (to the workspace). I also created another script run.r s

Re: [R] function input as variable name (deparse/quote/paste) ??

Thank you everyone for your reply. Like I said in my original post, this is just a demonstrative example of my 'big' self written script. My 'big' function take several inputs, of which the first 1 is the dataset and returns a LIST variable 'res<<-list()' to the workspace with many information.

[R] R CMD Rd2pdf [options] files

Hi all, Sorry if I seem a bit pissed because I am! Basically, I have written some script and I wanted to make a "documentation" for it using .Rd format. I followed the example here in http://stat.ethz.ch/R-manual/R-de

[R] Quicker way to apply values to a function

Hi all, myint=function(mu,sigma){ integrate(function(x) dnorm(x,mu,sigma)/(1+exp(-x)),-Inf,Inf)$value } mymu=seq(-3,3,length(1000)) mysigma=seq(0,1,length(500))[-1] k=1 v=c() for (j in 1:length(mymu)) { for (i in 1:length(mysigma)) { v[k]=myint(mymu[j],mysigma[i])

[R] Vectorize (scalar) function

Hi all, myint=function(mu,sigma){ integrate(function(x) dnorm(x,mu,sigma)/(1+exp(-x)),-Inf,Inf)$value } x=seq(0,50,length=3000) x=x[-1] plot(x,myint(4,x)) # not working yet I think I have to 'Vectorize' it somehow? What's a scalar function? and a primitive function? Thanks. casper

[R] show and produce PDF file with pdf() and dev.off( ) in function

Hi all, I know how to use pdf() and dev.off() to produce and save a graph. However, when I put them in a function say myplot(x=1:20){ pdf("xplot.pdf") plot(x) dev.off() } the function work. But is there a way show the graph in R as well as saving it to the workspace? Thanks. casper --

Re: [R] Quicker way to apply values to a function

Hi Petr, Thanks for confirming that the integral is bounded. I was thinking about the same thing. However, this requires that 'sigma' is positive. The actual problem occurred in my optimization routine, where i have set the parameter sigma=exp(para), where para is the logit of a uniform random

[R] R numerical integration

Hi all, Is there any other packages to do numerical integration other than the default 'integrate'? Basically, I am integrating: integrate(function(x) dnorm(x,mu,sigma)/(1+exp(-x)),-Inf,Inf)$value The integration is ok provided sigma is >0. However, when mu=-1.645074 and sigma=17535.26 It sto

Re: [R] R numerical integration

The quadinf command in library pracma still fails when mu=-2.986731 with sigma=53415.18. While Maple gives me an estimate of 0.5001701024. Maple: (for those who are interested) myf:=(mu,sigma)-> evalf(Int(exp(-(x-mu)^2/2/sigma^2)/sigma/sqrt(2*Pi)/(1+exp(-x)

[R] Create Model Object (setClass?setMethod?)

Hi all, I have a self written likelihood as a model and functions to optimize and get fitted values, confidence intervals ect. I wonder if there is a way to define a 'class', or a 'model' (or a certain object)? so that I can use 'summary' to produce a summary like it does for a lm object. Also,

[R] how to avoid this : underflow occurred in 'lgammacor'

Hi all, I am constructing a likelihood involving the following lbeta(j + a, k - j + b) where j,k are constants and a and b are parameters (>0). While doing the optimization, the error sometimes occurs, In lbeta(j + a, k - j + b) : underflow occurred in 'lgammacor' Is there a way to avoid it?

[R] Numerical integration again

Hi all, Here is an integration function require(pracma) # for 'quadinf' myint=function(j) { quadinf(function(x) (1/(1+exp(-x)))^j*(1-1/(1+exp(-x)))^(k-j)*dnorm(x,mu,casigma),-Inf,Inf) } in any optimization routine. It works fine most of the time but failed with some particular sets of va

[R] something weird in integration (pracma library)

Hi, library(pracma) k=20 mu=4.5 casigma=17000 myint=function(j) { quadinf(function(x) (1/(1+exp(-x)))^j*(1-1/(1+exp(-x)))^(k-j)*dnorm(x,mu,casigma),-Inf,Inf) } sapply(0:k,myint) ##

[R] problem with Gauss Hermite ( x and w )

Hi all, I am using the 'gaussHermite' function from the 'pracma' library CODES ### library(pracma) cc=gaussHermite(10) cc$x^2 cc$x^5 cc$x^4 CODES ### as far so good. However, it does NOT work for any NON integer values, say CODES ##

Re: [R] problem with Gauss Hermite ( x and w )

Hi, I know what complex number are, but I am not sure what you meant by that? ##CODES### > 2.5^(-2.4) [1] 0.1109032 > -2.5^(-2.4) [1] -0.1109032 ##CODES### works fine. Negative powers mean they take the reciprocal and as far as I am concerned, real^real

Re: [R] problem with Gauss Hermite ( x and w )

Rui Barradas wrote > > > real^real is not necessarily real. > > The most well known example is (-1)^0.5 = imaginary unit. > Damn, can't believe it! It's a silly mistake! Now that something wonders me is that when applying the Gaussian Hermit, sum w f(x_i) What happens when f(x_i) does not

[R] triangular matrices input/output

Hi, Is there any package that deals with triangular matrices? Say ways of inputting an upper (lower) triangular matrix? Or convert a vector of length 6 to an upper (lower) triangular matrix (by row/column)? Thanks! - ## PhD candidate in Statistics Big R Fan Big LEGO Fan