search for: sigma2

The background: I'm trying to fit a Poisson-lognormal distrbutuion to some data. This is a way of modelling species abundances: N ~ Pois(lam) log(lam) ~ N(mu, sigma2) The number of individuals are Poisson distributed with an abundance drawn from a log-normal distrbution. To fit this to data, I need to integrate out lam. In principle, I can do it this way: PLN1 <- function(lam, Count, mu, sigma2) { dpois(Count, exp(lam), log=F)*dnorm(LL, mu, sqrt(sigm...

Dear R developers, sessionInfo() below Please have a look at the following two versions of the same function: 1. Intended behavior: > Su1 = function(u=100, l=u, mu=0.53, sigma2=4.3^2) + { + print(c(u, l, mu)) # here, l is set to u?s value + u = u/sqrt(sigma2) + l = l/sqrt(sigma2) + mu = mu/sqrt(sigma2) + print(c(u, l, mu)) + } > > Su1() [1] 100.00 100.00 0.53 [1] 23.2558140 23.2558140 0.1232558 In the first version, both u and l are correctly divided...

Another way to avoid the problem is to not redefine variables that are arguments. E.g., > Su3 <- function(u=100, l=u, mu=0.53, sigma2=4.3^2, verbose) { if (verbose) { print(c(u, l, mu)) } uNormalized <- u/sqrt(sigma2) lNormalized <- l/sqrt(sigma2) muNormalized <- mu/sqrt(sigma2) c(uNormalized, lNormalized, muNormalized) } > Su3(verbose=TRUE) [1] 100.00 100.00 0.53 [1] 23.2558140 23....

...This is not really an R question but a statistical one. If someone could either give me the brief explanation or point me to a reference that might help, I'd appreciate it. I want to estimate the mean of a log-normal distribution, given the (log scale normal) parameters mu and sigma squared (sigma2). I understood this should simply be: exp(mu + sigma2) ... but I the following code gives me something strange: R <- 10000000 mu <- -400 sigma2 <- 200 tmp <- rlnorm(R, mu, sqrt(sigma2)) # a sample from the desired log-normal distribution muh <- mean(log(tmp)) sigma2h <- var(log...

...elps, Rui Barradas Citando Matthias Gondan <matthias-gondan at gmx.de>: > Dear R developers, > > sessionInfo() below > > Please have a look at the following two versions of the same function: > > 1. Intended behavior: > >> Su1 = function(u=100, l=u, mu=0.53, sigma2=4.3^2) > + { > + print(c(u, l, mu)) # here, l is set to u?s value > + u = u/sqrt(sigma2) > + l = l/sqrt(sigma2) > + mu = mu/sqrt(sigma2) > + print(c(u, l, mu)) > + } >> >> Su1() > [1] 100.00 100.00 0.53 > [1] 23.2558140 23.2558140 0.1232558 > &...

...ore u is assigned a new value. Of course, this would require the R interpreter to track all these dependencies in both ways incl. more complicated ones in which L might depend on more than just u. In the future, I?ll avoid dependencies between parameters. Su4 <- function(u=100, l=100, mu=0.53, sigma2=4.3^2) # instead of l=u And maybe also ?in-place? changes of values? Best regards, Matthias Von: William Dunlap Gesendet: Samstag, 2. September 2017 19:41 An: Rui Barradas Cc: Matthias Gondan; r-help at r-project.org Betreff: Re: [R] Strange lazy evaluation of default arguments Another way to...

Dear R users, I have two nonlinear equations, f1(x1,x2)=0 and f2(x1,x2)=0. I try to use optim command by minimize f1^2+f2^2 to find x1 and x2. I found the optimal solution changes when I change initial values. How to solve this? BTW, I also try to use grid searching. But I have no information on ranges of x1 and x2, respectively. Any suggestion to solve this question? Thanks, Kate

Hi. I am trying to analyze with BRugs the Box-Tiao variance components example in WinBUGS. The output from BRugs, mean sd MC_error val2.5pc median val97.5pc start sample sigma2.btw 681.9 1161 10.89 0.7016 253.8 4232 25001 100000 sigma2.with 4266.0 1246 4.92 2480.0000 4057.0 7262 25001 100000 doesn't match the output from WinBUGS, node mean sd MC error 2.5% median 97.5% start sample sigma2.btw 2929.0 2191.0 12.0...

...y, in Stata my coefficient estimates are identical in both the OLS estimation by -reg- and the maximum likelihood estimation using Stata's ml lf command. In R my coefficient estimates differ slightly between lm() and my own maximum likelihood estimation. Secondly, the estimates for sigma2 are very different between R and Stata. 3.14 in R compared to 1.78 in Stata. I have copied my maximum likelihood program below. It is copied from a great intro to MLE in R by Macro Steenbergen http://artsci.wustl.edu/~jmonogan/computing/r/MLE_in_R.pdf Any comments are welcome. In particular...

Dear all, I am stuck on the following problem with integrate(). I have been out of luck using RSiteSearch().. My function is g2<-function(b,theta,xi,yi,sigma2){ xi<-cbind(1,xi) eta<-drop(xi%*%theta) num<-exp((eta + rep(b,length(eta)))*yi) den<- 1 + exp(eta + rep(b,length(eta))) result=(num/den)*exp((-b^2)/sigma2)/sqrt(2*pi*sigma2) r=prod(result) return(r) } And I am interested in evalua...

...ject.org > Subject: [R] Strange lazy evaluation of default arguments > > Dear R developers, > > sessionInfo() below > > Please have a look at the following two versions of the same function: > > 1. Intended behavior: > > > Su1 = function(u=100, l=u, mu=0.53, sigma2=4.3^2) > + { > + print(c(u, l, mu)) # here, l is set to u?s value > + u = u/sqrt(sigma2) > + l = l/sqrt(sigma2) > + mu = mu/sqrt(sigma2) > + print(c(u, l, mu)) > + } > > > > Su1() > [1] 100.00 100.00 0.53 > [1] 23.2558140 23.2558140 0.1232558 >...

...loess() and return the mimimizing loess fit for a specified criterion. loess.aic <- function (x) { # extract values from loess object if (!(inherits(x,"loess"))) stop("Error: argument must be a loess object") span <- x$pars$span n <- x$n traceL <- x$trace.hat sigma2 <- sum( x$residuals^2 ) / (n-1) delta1 <- x$one.delta delta2 <- x$two.delta enp <- x$enp aicc <- log(sigma2) + 1 + 2* (2*(traceL+1)) / (n-traceL-2) aicc1<- n*log(sigma2) + n* ( (delta1/(delta2*(n+enp)))/(delta1^2/delta2)-2 ) gcv <- n*sigma2 / (n-traceL)^2 result &lt...

Dear UseRs, I wrote the following function to use MLE. --------------------------------------------- mlog <- function(theta, nx = 1, nz = 1, dt){ beta <- matrix(theta[1:(nx+1)], ncol = 1) delta <- matrix(theta[(nx+2):(nx+nz+1)], ncol = 1) sigma2 <- theta[nx+nz+2] gamma <- theta[nx+nz+3] y <- as.matrix(dt[, 1], ncol = 1) x <- as.matrix(data.frame(1, as.matrix(dt[, 2:(nx+1)], ncol = 2))) z <- as.matrix(dt[, (nx+2):(nx+nz+1)], ncol = nz) d <- z %*% delta / (gamma * sigma2)^.5 mustar <- (1-gamma) * z %*%...

Hi, A question. When I run gelman.diag and gelman.plot with mcmc lists obtained from MCMCregress, the results are following. > post.R <- MCMCregress(Size~Age+Status, data = data, burnin = 5000, mcmc = 100000, + thin = 10, verbose = FALSE, beta.start = NA, sigma2.start = NA, + b0 = 0, B0 = 0, nu = 0.001, delta = 0.001) > post1.R <- MCMCregress(Size~Age+Status, data = data, burnin = 5000, mcmc = 100000, + thin = 10, verbose = FALSE, beta.start = -3, sigma2.start = NA, + b0 = 0, B0 = 0, nu = 0.001, delta = 0.001) > post2.R &...

...help at stat.math.ethz.ch> References: <4405D98F.6030709 at dssm.unipa.it> vito muggeo wrote: > Dear all, > I am stuck on the following problem with integrate(). I have been out of > luck using RSiteSearch().. > > My function is > > g2<-function(b,theta,xi,yi,sigma2){ > xi<-cbind(1,xi) > eta<-drop(xi%*%theta) > num<-exp((eta + rep(b,length(eta)))*yi) > den<- 1 + exp(eta + rep(b,length(eta))) > result=(num/den)*exp((-b^2)/sigma2)/sqrt(2*pi*sigma2) > r=prod(result) > return(r) >...

...nd Ogutu 2002), to test for a trend over time, using multiple sites: y[ij]=mu+b[j]+a[i]+w[j]*(beta +t[i])+c[ij] where: y[ij]= a response variable at site i and year j mu = fixed intercept Beta=fixed slope w[j]=constant representing the jth year (covariate) b[j]=random effect of jth year, iid N(0,sigma2[b]) a[i]=random effect of the ith site, iid N(0, sigma2[a]) t[i]=random effect of ith site, iid N(0, sigma2[t]) c[ij]=random error associated with ith site and jth year I would like to assume that an unstructured relationship applies to a[i] and t[i] (i.e., I would like to assume that the random e...

...therefore the log is -Inf. > > Where is my mistake? > > My code: > > # EM algorithm manually > # dat is the data > > > # initial values > pi1<-0.5 > pi2<-0.5 > mu1<--0.01 > mu2<-0.01 > sigma1<-0.01 > sigma2<-0.02 > loglik[1]<-0 > loglik[2]<-sum(pi1*(log(pi1)+ > log(dnorm(dat,mu1,sigma1))))+sum(pi2*(log(pi2)+log(dnorm(dat,mu2,sigma2)))) > > > > tau1<-0 > tau2<-0 > k<-1 > > # loop > while(abs(loglik[k+1]-loglik[k])>= 0....

...in R. The equation I want to estimate is: p(t)=a+b*p(t-1)+error Using STATA I get 0.92 for a, and 0.73 for b. Code that I use in R is: p<-matrix(data$p) # price at time t lp<-cbind(1,data$lp) # price at time t-1 mle <- function(theta) { sigma2<-theta[1] b<- theta[-1] n<-length(p) e<-p-lp%*%b logl<- -(n/2)*log(sigma2)-((t(e)%*%e)/(2*sigma2)) return(-logl) } out <- optim(c(0,0,0),mle, method = "L-BFGS-B", lower = c(0, -Inf, -Inf), upper = c(Inf, Inf, Inf)) The "result" I get is:...

I am sorry to bug you, I am having this error whenever I want to run random effects regression in software R: Error in if (sigma2$id < 0) stop(paste("the estimated variance of the", : missing value where TRUE/FALSE needed. Please help me look into it. [[alternative HTML version deleted]]

Hi, I am trying to fit a function of the form: y = A0 + A1 * exp( -0.5* ( (X - Mu1) / Sigma1 )^2 ) - A2 * exp ( -0.5* ( (X-Mu2)/Sigma2 )^2 ) i.e. a mean term (A0) + a difference between two gaussians. The constraints are A1,A2 >0, Sigma1,Sigma2>0, and usually Sigma2>Sigma1. The plot looks like a "Mexican Hat". I had trouble (poor fits) fitting this function to toy data in Matlab and now I am playing with...