search for: sig2

...meters*? It?s easy enough to implement myself, but I?d like to avoid reinventing the wheel (with some bugs) if possible. Here?s a simple example of the result I?d like, using numerical differentiation: library(mvtnorm) library(numDeriv) f=function(pars, xx, yy) { mu=pars[1:2] sig1=pars[3] sig2=pars[4] rho=pars[5] sig=matrix(c(sig1^2,rho*sig1*sig2,rho*sig1*sig2,sig2^2),2) dmvnorm(cbind(x,y),mu,sig) } mu1=1 mu2=2 sig1=3 sig2=4 rho=.5 x=2 # or a x vector y=3 # or a y vector jacobian(f,c(mu1,mu2,sig1,sig2,rho),xx=x,yy=y) # (Can replace ?jacobian? with ?grad? if x and y have length 1.)...

..., which I can't check.) I believe the graphs they show, but when I code the function in R, according to the listed formula, I get very different graphs. Would someone please point out my error or tell me where it's already implemented in R? Here is my version: > sddist function(s,n) { sig2 <- n*s*s/(n-1) 2*(n/(2*sig2))^((n-1)/2) / gamma((n-1)/2) * exp(-n*s*s/(2*sig2)) * s^(n-2) } Version 2.3.1 (2006-06-01) on Windows XP SP2 Thanks for any help. David L. Reiner Rho Trading Securities, LLC Chicago? IL? 60605 312-362-4963

...m is to compare means of two independent samples that results in Beherens-Fisher posterior and in the model.file of WinBUGS "d2.bug" there is the following codes: model { for( i in 1 : 30 ) { x[i] ~ dnorm(mu1,sig1) } for( i in 1 : 30 ) { y[i] ~ dnorm(mu2,sig2) } mu1 ~ dnorm(50,1.0E-6) sig1 ~ dgamma(0.001,0.001) mu2 ~ dnorm(50,1.0E-6) sig2 ~ dgamma(0.001,0.001) lbda <- mu1 - mu2 } I??m a new user of WinBUGS and if someone detect error in the model codes too, I??m grateful. Thanks for help! Gilberto Matos. ----------...

...gma) mlebinorm<- function(startval, y)# startval = Initial Values to be passed , y= data { startval<-as.numeric(startval) oneside=matrix(c(0,0,0,0,-1,0,0,0,0,1),5,2) otherside = c(-0.9999999,-0.999999) n<- ncol(y) nf<- function(x) { mu1<-x[1] mu2<-x[2] sig1<-x[3] sig2<-x[4] rho<-x[5] nf<-(-n/2)*(log(sig1*sig2*(1-rho^2)))-(0.5/(1-rho^2)* (((y[1,1]-mu1)^2/sig1)- (2*rho*(y[1,1]-mu1)*(y[1,2]-mu2))/sqrt(sig1*sig2)+ (y[1,2]-mu2)^2/sig2)) for(i in 2:n) { nf<- nf - (0.5*(1-rho^2)* (((y[i,1]-mu1)^2/sig1)- (2*rho*(y[i,1]-mu1)*(y[i,2]-...

...problem below: R message is like below: ------- Error in garch.gjr.d0(param, y, X, Z, iterate = FALSE)$res : $ operator is invalid for atomic vectors plot(Ch.mic,Ch.mic.rep,xlim=c(0,6),ylim=c(0,6)) ----------------------- Coding -------- res<-garch.gjr.d0(a$estimate,y,X,X,iterate=FALSE)$res sig2<-garch.gjr.d0(a$estimate,y,X,X,iterate=FALSE)$sig2 test.mic <- function(param,y,X,Z) { n <- length(y) res <- garch.gjr.d0(param,y,X,Z,iterate=FALSE)$res s2 <- garch.gjr.d0(param,y,X,Z,iterate=FALSE)$sig2 s<-sqrt(s2) r <- res u <- r[-1] r2<- (...

...scale' to use when finding, for example, the MLEs of a bivariate normal distribution using optim. Here's code for this example: ----------------------------------------------- library(MASS) # needed mvrnorm library(mvtnorm) # needed for dmvnorm set.sed(20080516) n=1000 mu1=3 mu2=5 sig1=7 sig2=20 rho=.5 sigmat=matrix(c(sig1^2,sig1*sig2*rho,sig1*sig2*rho,sig2^2),2) xy=mvrnorm(n,c(mu1,mu2),sigmat) # n = 1000 observations from this # distribution. obj=function(par,xy) # The function to maximize. { mu=par[1:2] sigmat=matrix(c(par[3]^2,par[3]*par[4]*par[5],par...

Hi, R-listers I tried to plot several graphs in a sigle x-y coordinate settings, like the following: |(y) s | ****** s | ***** s | sssssssssssssssssss |_______________________________(x) where "*" and "s" denote two diffrent plots. However, when I used plot(data1); % data1 is the data points of "*"

.... Error in integrate(slice, 0, Inf, , , , , , , a, b) * delta : non-numeric argument to binary operator [ By the way, I'm trying to evaluate a two-dimensional integral by slicing it up into one-dimensional bits which I will loop over to evaluate.] Here's the code: mu <- 0.3 sig2 <- 0.01 alpha <- 20 beta <- 15.75 rho <- 0.1 k <- 1/(rho^2.5*gamma(rho)^2*sqrt(2*pi*sig2)) slice <- function(w,a,b) { g <- w^(1/rho) g1 <- w1^(1/rho) h <- g1^a*g^b E <- -(y-rho*mu -g1/alpha + g/beta)^2/(2*sig2*rho) k*h*exp(E- g1 - g) } hi <- 10 delta &lt...

...0.00 13.6 218.1 36.8 -999 65.1 -999 -999 18.2 0.2 0506010005 0.00 12.9 214.3 36.8 -999 65.5 -999 -999 16.9 0.2 0506010010 0.00 14.4 215.7 36.9 -999 65.4 -999 -999 20.4 0.2 *dataM.txt* doy yr mon day hr hgt1 hgt2 hgt3 co21 co22 co23 sig1 sig2 sig3 dif flag 244.02083 2005 09 01 00 2.5 5.8 9.1 -999.99 -999.99 -999.99 -999.99 -999.99 -999.99 -999.99 PRE 244.0625 2005 09 01 01 2.5 5.8 9.1 -999.99 -999.99 -999.99 -999.99 -999.99 -999.99 -999.99 PRE 244.10417 2005 09 01 02 2.5 5.8 9.1 -999.99 -999.99 -999.99 -999.99 -999.99 -999.99 -999.9...

...to his post.  Someone could help me? ## commands BBtraj <- list() for (i in 1:(nrow(loc@data)-1)) { BBtraj[[i]] <- kernelbb(as.ltraj(loc@coords[i:(i+1),], date=loc@data$time[i:(i+1)], id = as.character(loc@data$id[i:(i+1)]), typeII = TRUE), sig1=as.numeric(as.character(loc@data$sig1[i])), sig2= 5, grid = 1000) }   Rémi Lesmerises, biol. M.Sc., Candidat Ph.D. en Biologie Université du Québec à Rimouski remilesmerises@yahoo.ca [[alternative HTML version deleted]]

...t;-array(0,3) # array for diagnostics + hat<-array(0,n) # array for diagnostics Vp[1,1]<-1.0 M$gcv.ubre<-1.0; direct.mesh<-100 # number of points for overall s.p. initial direct search @@ -448,7 +449,8 @@ as.integer(n),as.integer(q),as.integer(C.r),as.double(M$sig2),as.double(Vp), as.double(edf),as.double(M$conv.tol),as.integer(M$max.half),as.double(ddiag) , as.integer(idiag),as.double(sdiag),as.integer(direct.mesh),as.double(M$min.e df), - as.double(M$gcv.ubre),as.double(M$target.edf),PACKAGE="mgcv") + as.double(M$gcv.ubre),as.d...

Dear List I am using R with mcgv package to model spatial variation in density estimates of dorcas gazelle in Sinai. I have 59 points of data and 4 explanatory variables(distance from mountain edge, camel presence, Latitude & Longitude). I want to test the model fir via bootstraping. I have used the jacknife bootstraping but it have the limitation of allowing only 58 trials. I tried to use the

All, Sorry for bombarding you with GAM related questions, but... I know a partial residual option in plot.gam() is on Simon Wood's todo list, but since I'm in the midst of a project and not yet having acquired sufficient R knowledge to code something usable myself I'll have to put my trust in you. Anybody got some code lying around for doing this? Or if someone can supply me with

...3) # and this is what I did with the help of Dennis Murphy last week: myframe3$day <- as.Date(myframe3$myframestime) hrfun <- function(d) { dxy <- d[, c("Long", "Lat")] d_ltraj <- as.ltraj(dxy, d$myframestime, id = d$AnimalID) dLiker <- liker (d_ltraj, sig2 = 18, rangesig1= c(1, 100) ) dBB <- kernelbb(d_ltraj, sig1=dLiker$Ernie$sig1, sig2=18) data.frame(AnimalID = unique(d$AnimalID), day = unique(d$day), h = kernel.area(dBB, unout=c("km2"))) } ddply(myframe3, .(day), hrfun) # i.e. I calculated the size of the homerange for Animal...

...pValue of the calculated statistic with respect of the distribution of this statistic under the null hypothesis. 2-how is it possible to test an overdispersion in the poisson model? for choosing a best model I need this mesure of dispersion. Should I build a glm(cases~expected,family=quasipoisson)$sig2 or is it possible directly in DClusters? Thank you, for advance for answers. Best regards, Erik A. ============================================ Erik-André SAULEAU SEAIM Hôpital du Hasenrain BP 1070 87, Avenue de Altkirch 68051 MULHOUSE Cédex Tel: 03-89-64-79-95 Fax: 03-89-64-79-71 Mél: sauleau...

...var(y) e[1] <- y[1]-mu0 l<-0 for(t in 2:T) { e[t] <- y[t]-mu0 s2[t] <- s2[t]+ a0+a1*e[t-1]^2+b1*s2[t-1] + ifelse(e[t-1]<0, a2*e[t-1]^2, 0) l <- l -0.5*log(2*pi*s2[t])-0.5*e[t]^2/s2[t] if(s2[t]>20000) return(20000) } if(iterate) return(-l) else return(list(loglik=l,sig2=s2,res=e/sqrt(s2))) } res<-garch.gjr.d0(a$estimate,y,X,X,iterate=FALSE)$res sig2<-garch.gjr.d0(a$estimate,y,X,X,iterate=FALSE)$sig2 plot(ts(y)) lines(-2*sqrt(sig2),lty=2) lines(2*sqrt(sig2),lty=2) sims.matrix <- rbind(sims[-(1:n.iter/2),1,],sims[-(1:n.iter/2),2,],sims[-(1:n.iter/2),3,]...

Dear All- I have a dataset that is comprised of the following: doy yr mon day hr hgt1 hgt2 hgt3 co21 co22 co23 sig1 sig2 sig3 dif flag 244.02083 2005 09 01 00 2.6 9.5 17.8 375.665 373.737 373.227 3.698 1.107 0.963 -0.509 PRE 244.0625 2005 09 01 01 2.6 9.5 17.8 393.66 384.773 379.466 15.336 11.033 5.76 -5.307 PRE 244.10417 2005 09 01 02 2.6 9.5 17.8 411.162 397.866 387.755 6.835 5.61 6.728 -10.112 PRE 244.14583 20...

I have estimated a multiple nonparametric regression using the loess command in R. I have also estimated an additive version of the model using the gam function. Is there a way of using the output of these two models to test the restrictions imposed by the additive model?

Hello, I tried to construct my very first loop today and completly failed :-( Maybe someone can help me? I have a dataframe somewhat like this one: myframe <- data.frame (Timestamp=c("24.09.2012 09:00", "24.09.2012 10:00", "24.09.2012 11:00", "25.09.2012 09:00", "25.09.2012 10:00", "25.09.2012

Hi, I'm in dire need of a fast fourier transformation for me stupid biologist, i.e. I have a heartbeat signal and would like to decompose it into pure sin waves, getting three vectors, one containing the frequencies of the sin waves, one the magnitudes and one the phases (that's what I get from my data acquisition software's FFT function). I'd be very much obliged, if someone