search for: theta

...o optimize for a set of parameters and want to set a constraint on only one parameter. Here is my function. What I want to do is estimate the parameters of a bivariate normal distribution where the correlation has to be between -1 and 1. Would you please advise how to revise it? ex=function(s,prob,theta1,theta,xa,xb,xc,xd,t,delta) { expo1= exp(s[3]*xa+s[4]*xb+s[5]*xc+s[6]*xd) expo2= exp(s[9]*xa+s[10]*xb+s[11]*xc+s[12]*xd) expo3= exp(s[15]*xa+s[16]*xb+s[17]*xc+s[18]*xd) expo4= exp(s[21]*xa+s[22]*xb+s[23]*xc+s[24]*xd) expo5= exp(s[27]*xa+s[28]*xb+s[29]*xc+s[30]*xd) nume1=prob[1]*(s[2]^-s[1]*s[1]*t...

...(-1) - 1 + x - x^2) s <- expression(x^1.3) s.x <- expression(1.3*(x^0.3)) # simulate process using the Milstein approximation X <- sde.sim(t0=0,T=100,X0=1,N=1000,drift=d,sigma=s,sigma.x=s.x,method="milstein") ## define the Lamperti transform function Transform <- function(t,x,theta) (1/(theta[5]*(theta[6]-1)))*(x^(1-theta[6])) ## define the diffusion function S <- function(t,x,theta) theta[5]*(x^theta[6]) ## define the transformed drift function and it's first six derivatives m0 <- function(t,x,theta) { b <- theta[5]*(theta[6]-1) p1 <- (theta[6]+1)/(theta[6]...

I got some R code that I don't understand. Question as comment in code //where is t comming from, what is phi inverse rAC <- function(name, n, d, theta){ #generic function for Archimedean copula simulation illegalpar <- switch(name, clayton = (theta < 0), gumbel = (theta < 1), frank = (theta < 0), BB9 = ((theta[1] < 1) | (theta[2] < 0)), GIG = ((theta[2] < 0) | (theta[3] < 0) | ((theta[1]>0) & (theta[3]==0)) | ((thet...

Hello, I want to have a title that will look something like: "Results for \theta=2.1", given that I have a variable theta=2.1, and \theta should show on the screen like the greek letter. I've tried a lot of things: theta <- 2.1 plot(1:10, main=expression(paste("Results for", theta, "=", eval(theta)))) or using bquote plot(1:10, main=paste(&quot...

dlogl <- -(n/theta)-sum((y/(theta)^2)*((1-exp(y/theta))/(1+exp(y/theta))) d2logl <- (n/theta^2) - sum((-2y/theta^3)*(1-exp(y/theta))/(1+exp(y/theta))) - sum(((2*y/theta^4)*exp(y/theta))/((1+exp(y/theta))^2)) returns the error message: Error: unexpected symbol in: "dlogl <- -(n/theta)-sum((y/(theta)^2)*((...

Hi, here comes my problem, say I have the following functions (example case) #------------------------------------------------------------ function1 <- function (x, theta) {a <- theta[1] ( 1 - exp(-theta[2]) ) * theta[3] ) b <- x * theta[1] / theta[3]^2 return( list( a = a, b = b )) } #----------------------------------------------------------- function2<-function (x, theta) {P <- function1(x, theta) c <- P$a * x * exp(-theta[2]) d <- P$b *...

Hello, I have written a function that demonstrates the CLT by generating samples following the exponential distribution, calculating the means, plotting the histogram, and drawing the limiting normal curve as an overlay. I have the title of each histogram state the sample size and rate (1/theta) for the exponential (the output is actually 4 histograms), but I can't get the greek letter theta to appear in the heading. Is this possible? I can get nontext for the xlab and ylab pretty easily. Thanks, Jason Below is my code. Currently it just writes out "theta." function (...

...te(numint, lower = -Inf, upper = Inf) : non-finite function value." My questions are How could I fix it? I tried to divide into several intervals and sum up, but I got same message. My code is ----------------------------------------------------------------------------------- alpha=1.3; theta0=c(0.5,0,sqrt(10),0.25,-.3,sqrt(0.05),.3,sqrt(0.05)) fc = function(theta){ numint = function(z) { (theta[1]*dnorm(z,mean=theta[2],sd=theta[3])+theta[4]*dnorm(z,mean=theta[5],sd=theta[6]) +(1-(theta[1]+theta[4]))*dnorm(z,mean=theta[7],sd=theta[8]))^alpha } inte = inte...

...e efficient code? # Begin script__________________________________________________ # Dichotomous scores for 100 respondents on 3 items with # probabilities of a correct response = .6, .4, and .7, # respectively x1 <- rbinom(100,1,.6) x2 <- rbinom(100,1,.4) x3 <- rbinom(100,1,.7) # 'theta.vec' is a vector holding 31 possible levels of theta # ranging from -3 to +3 in intervals of .2. theta.vec <- seq(-3,3,.2) theta <- sample(rep(theta.vec,5),100) x.mat <- (cbind(x1,x2,x3,theta)) rm(x1,x2,x3,theta) nc <- ncol(x.mat) ni <- nc - 1 nr <- nrow(x.mat) ntheta <-...

Hello, I have written a simple Metropolis-Hastings MCMC algorithm for a binomial parameter: MHastings = function(n,p0,d){ theta = c() theta[1] = p0 t =1 while(t<=n){ phi = log(theta[t]/(1-theta[t])) phisim = phi + rnorm(1,0,d) thetasim = exp(phisim)/(1+exp(phisim)) r = (thetasim)^4*(1-thetasim)^8/(theta[t]^4*(1-theta[t])^8) if(runif(1,0,1)<r){ theta[t+1] = thetasim } else { theta[t+1] = theta[t]...

...unnecessary braces, but I can't find any. False positive or am I missing something? If a false positive, is there any way to work around the warning? * checking Rd files against version 2 parser ... WARNING Warning: ./man/dbetabinom.Rd:32-34: Unnecessary braces at ?{p(x) = % (C(N,x)*Beta(N-x+theta*(1-p),x+theta*p))/% Beta(theta*(1-p),theta*p)}? ============= \details{ The beta-binomial distribution is the result of compounding a beta distribution of probabilities with a binomial sampling process. The density function is \deqn{p(x) = % \frac{C(N,x) \mbox{Beta}(N-x+\theta(1-p),x+\theta p)...

text() can draw text on a plot. Do we have a way/function to draw text with a box surround it? Thanks, -james

I need to wrap a loop inside a function and am having a small bit of difficulty getting the results I need. Below is a replicable example. # define functions pcm <- function(theta,d,score){ exp(rowSums(outer(theta,d[1:score],'-')))/ apply(exp(apply(outer(theta,d, '-'), 1, cumsum)), 2, sum) } foo <- function(theta,items, score){ like.mat <- matrix(numeric(length(items) * length(theta)), ncol = length(theta)) for(i in 1:length(item...

...it seems the mapply function in this case is slower (in terms of CPU speed) than the for loop. Here is an example. # data needed for example items <- list(item1 = c(0,1,2), item2 = c(0,1), item3 = c(0,1,2,3,4), item4 = c(0,1), item5=c(0,1,2,3,4), item6=c(0,1,2,3)) score <- c(2,1,3,1,3,2) theta <- c(-1,-.5,0,.5,1) # My old function using the for loop like.mat <- function(score, items, theta){ like.mat <- matrix(numeric(length(items) * length(theta)), ncol = length(theta)) for(i in 1:length(items)) like.mat[i, ] <- pcm(theta, items[[i]], score[[i]]) like.mat...

...(I'm running R 2.3.1, running takes ~1 min). library(survival); set.seed(10000) lambda <- 0.01 # Exp. hazard rate # Beta coefficients for Age,TC,HDLC,SBP,Diab,Smok beta <- c(0.0483,0.0064,-0.0270,0.0037,0.4284,0.5234) n <- 1000; Ngrp <- 2; # Nr patients, Nr frailty groups # Thetas for gamma-frailty thetaset <- c(1,2,10,100); Ntheta <- length(thetaset); # Define the simulated population age <-rnorm(n,48.6,11.7);tc<-rnorm(n,200,30) hdlc<-rnorm(n,47,6);sbp<-rnorm(n,135,6) rtmp <- runif(0,1,n=n); diab <- rep(0, n); diab[rtmp < 0.05223] <...

...ncing one small, albeit important problem. Here is part of my function below. This is a psychometric problem for dichotomously scored test items where x is a vector of 1s or 0s denoting whether the respondent answered the item correctly (1) or otherwise (0), b is a vector of item difficulties, and theta is an ability estimate for the individual. rasch <- function(b,theta){ 1 / ( 1 + exp(b - theta)) } The function rasch gives the probability of a correct response to item i conditional on theta, the individuals ability estimate myfun <- function(x, b, theta){ sum(rasch(b, theta)^x...

Dear R Users, When using zeroinfl() function to fit a Zero-Inflated Negative Binomial (ZINB) model to a dataset, the summary() gives an estimate of log(theta) and its standard error, z-value and Pr(>|z|) for the count component. Additionally, it also provided an estimate of Theta, which I believe is the exp(estimate of log(theta)). However, if I would like to have an standard error of Theta itself (not the SE.logtheta), how would I obtain or calcula...

...1]<-NA sample.deleted[delete.two,2]<-NA missing<-c(delete.one,delete.two) complete<-sample.deleted[!(is.na(sample.deleted[,1]) | is.na(sample.deleted[,2])),] deleted<-sample.deleted[missing,] Try to estimate the parameters of the data set less the deleted values: exact<-function(theta,complete,deleted){ one.only<-deleted[!(is.na(deleted[,1])),1] two.only<-deleted[!(is.na(deleted[,2])),2] u<-c(theta[1],theta[2]) sigma<-c(theta[3],theta[5],theta[5],theta[4]) dim(sigma)<-c(2,2) -sum(log(dmvnorm(x=complete,mean=u,sigma=sigma)))- sum(log...

Hello expeRts, I need generate symbolize the autocovariances matrix of a Gaussian ARMA(1,1), for derivate it and evaluate. I try this codes, but whitout sucess vacv<-NULL vacv[1]<-1-2*phi*theta-theta^2 vacv[2]<-(1-phi*theta)*(phi-theta) vacv[3:n]<-acv[2]*(phi^(1:(n-2))) facv<-list() for(i in 1:2) facv[[i]]<-deriv(y~vacv[i],c("phi","theta"),function(phi,theta){}) Erro en deriv.formula(y ~ vacv[i], c("phi", "theta"), function(phi, theta...

...found it useful for learning optim(), and you might find it similarly useful. I will be most grateful if you can guide me on how to do this better. Should one be using optim() or stats4::mle? set.seed(101) # For replicability # Setup problem X <- cbind(1, runif(100)) theta.true <- c(2,3,1) y <- X %*% theta.true[1:2] + rnorm(100) # OLS -- d <- summary(lm(y ~ X[,2])) theta.ols <- c(d$coefficients[,1], d$sigma) # Switch to log sigma as the free parameter theta.true[3] <- log(theta.true[3]) theta.ols[3] <- log(theta.ols[3]) # OLS likelihood function...