similar to: Problem with numerical integration and optimization with BFGS

Hi, I'm trying to write a simple function like case1 <- function (m, cov, Q, R) { theta <- (acos(R/sqrt(Q^3))) beta <- (-2)*sqrt(Q)*cos(theta/3)+m[1]/3 rho1 <- (-2)*sqrt(Q)*cos((theta+2*pi)/3)+m[1]/3 rho2 <- (-2)*sqrt(Q)*cos((theta-2*pi)/3)+m[1]/3 stderrb <- deltamethod( ~(-2)*sqrt(Q)*cos(theta/3)+x1/3,m,cov) stderrr1 <- deltamethod(

I am looking for examples of code that demonstrates the fine tuning of the strip panels in lattice graphics and uses plotmath characters. The code for the graphic is as follows: xyplot(lagy ~ n | rho1 * rho2, data= data, layout=c(2,6), span = 1, xlab = "Sample Size", ylab = "Bias in the Coefficient for the Lag of X", type = "o") rho1 is a four level factor

[Apologies in advance if this is too "statistics" and not enough "R".] I've got an experiment with two sets of treatments. Each subject either received all treatments from set A or all treatments from set B. I can compute the N pairwise correlations for all treatments in either set using cor(). If I take the mean of these N pairwise correlations, I see that the effects

Hi All, I am doing econometric modeling of panel data (fixed effects). We currently use Eviews to do this, but I have discovered a bug in Eviews 7 and am exploring the use of R to build panel data models / longitudinal data models. I looked at the plm package but do not see how I can incorporate AR terms in the model using the plm package. I have an Eviews model with two AR terms, AR(1) and

Dear All: I am trying to fit a special case of a 2-banded Toeplitz correlation structure. A 2-banded Toeplitz has ones on the diagonal, a correlation, RHO1, on the first off-diagonal, and a correlation, RHO2, on the second off-diagonal, with zeros on all subsequent off-diagonals. After reading relevant sections in Mixed-Effects Models in S and S-PLUS (Pinheiro & Bates, 2000) and searching

Dear All: I am trying to fit a special case of a 2-banded Toeplitz correlation structure. A 2-banded Toeplitz has ones on the diagonal, a correlation, RHO1, on the first off-diagonal, and a correlation, RHO2, on the second off-diagonal, with zeros on all subsequent off-diagonals. After reading relevant sections in Mixed-Effects Models in S and S-PLUS (Pinheiro & Bates, 2000) and searching

Dear All: I am trying to fit a special case of a 2-banded Toeplitz correlation structure. A 2-banded Toeplitz has ones on the diagonal, a correlation, RHO1, on the first off-diagonal, and a correlation, RHO2, on the second off-diagonal, with zeros on all subsequent off-diagonals. After reading relevant sections in Mixed-Effects Models in S and S-PLUS (Pinheiro & Bates, 2000) and searching

Dear list, I have a problem on persp() x <- u1data #first coloum in attached data y <- u2data #second coloum in attached data f <- function(x,y){qgev(pnorm(rhoF*qnorm(pnorm((qnorm(y)-rho2*qnorm(x)/sqrt(1-rho2^2)))) +sqrt(1-rhoF^2)*qnorm(0.95)),-0.3935119, 0.4227890, 0.2701648)} z <- outer(x,y,f) persp(x,y,z) The R will display: "Error in persp.default(x, y,

Hello everybody, I recently wrote a "program" that to generate AR1 time series, here the code: #By Jomopo. Junio-2007, Leioa, Vizcaya #This program to create the AR1 syntetic series (one by one) #Where the mean is zero, but the variance of the serie AR1 and #the coef. of AR1 are be changed. If var serie AR1 = 1 then is standarized! #Final version for AR1 time series program #Mon Jul

Hello, I'm quite new with R and so I would like to know if there is a command to calculate an integral. In particular I simulated a bivariate normal distribution using these simple lines: rbivnorm <- function(n, # sample size mux, # expected value of x muy, # expected value of Y sigmax, # standard deviation of

Hi, Using indicator variables I have been able to fit and run the code for fitting a bivariate mixed model using unstructured covariance matrix The code is lme.fit1<- lme(one.var~-1+indic1+indic2+I(indic1*d.time)+I(indic2*d.time), random =~ -1+indic1+indic2|m.unit, weights = varIdent(~1|indic1) ,data = new.data) My variables are one.var :- the two response variables stacked one after

Hi, I have a bivariate longitudinal dataset. As an example say, i have the data frame with column names var1 var2 Unit time trt (trt represents the treatment) Now suppose I want to fit a joint model of the form for the *i* th unit var1jk = alpha1 + beta1*timejk + gamma1* trtjk + delta1* timejk:trtjk + error1jk var2 = alpha2 + beta2*timejk + gamma2* trtjk + delta2* timejk:trtjk +

I have gone over the examples and can't figure this out: rho<-.77 text(x=.05,y=.5,paste(expression(rho),rho)) I was hoping to get this to print a Greek rho with 0.77 beside it. Instead I get: rho 0.77 (i.e. Roman lettering) The help on expression() is quite opaque so I don't understand how it works. Thanks for any help. Bill Simpson

As part of carrying out a complicated maximum likelihood estimation, I am trying to learn to program likelihoods in R. I started with a simple probit model but am unable to get the code to work. Any help or suggestions are most welcome. I give my code below: ************************************ mlogl <- function(mu, y, X) { n <- nrow(X) zeta <- X%*%mu llik <- 0 for (i in 1:n) { if

Dear all, I ran into problems with the function "optim" when I tried to do an mle estimation of a simple lognormal regression. Some warning message poped up saying NANs have been produced in the optimization process. But I could not figure out which part of my code has caused this. I wonder if anybody would help. The code is in the following and the data is in the attachment. da <-

Dear List, I am not sure how should i optimize a log-likelihood numerically: Here is a Text book example from Statistical Inference by George Casella, 2nd Edition Casella and Berger, Roger L. Berger (2002, pp. 355, ex. 7.4 # 7.2.b): data = x = c(20.0, 23.9, 20.9, 23.8, 25.0, 24.0, 21.7, 23.8, 22.8, 23.1, 23.1, 23.5, 23.0, 23.0) n <- length(x) # likelihood from a 2 parameter Gamma(alpha,

I'm looking at what would be involved in converting from MeetMe to ConfBridge and there seems to be a lot of missing administrative things, but I hope I'm just missing it. We all know about the missing realtime linkage. That's a major nuisance, but can be worked around. More serious is that the CLI command to display users in a ConfBridge don't show the caller ID information, so

I feel like a complete dolt, as I know this question has been asked by others on a fairly regular basis, but I'm going in circles trying to get the following to work: id.prob<-function (tt) { library(mvtnorm) #============================ Makeham<-function(tt) { a2=0.030386513 a3=0.006688287 b3=0.039047537 t<-tt-20 h.t<-a2+a3*exp(b3*t) S.t<-exp(-a2*t+a3/b3*(1-exp(b3*t)))

Hi: I am using R 1.7.0 on Windows. I am having trouble getting "outer" to work on one of my functions. Here is a simple example illustrating my problem: > b1 <- c(1.2,2.3) > b2 <- c(0.5,0.6) > x <- c(3e+01, 1e+02, 3e+02, 5e+02, 1e+03, 1e+04, 1e+05, 1e+06) > y <- c(2,4,2,5,2,3,1,1) > n <- c(5,8,3,6,2,3,1,1) > outer(b1,b2,FUN=bpllkd,x,y,n)

R Users: I am trying to estimate a model of fertility behaviour using birth history data with maximum likelihood. My code works but is extremely slow (because of several for loops and my programming inefficiencies); when I use the genetic algorithm to optimize the likelihood function, it takes several days to complete (on a machine with Intel Core 2 processor [2.66GHz] and 2.99 GB RAM). Computing