[R] Two dimensional likelihood surface plot
Hi R users, I am trying to plot two dimensional posterior likelihood surface. I have a data like para1 para2 likehood ... ... ... ... I looked at contour plot but it needs a shorted values of parameters and a matrix of likelihood values. Is there any way to get the plot? or how can I change my likelihood values to a matrix for the function contour? Any suggestions are appreciated. GP University of Guelph Guelph, ON [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Two dimensional likelihood surface plot
Thanks David, what do you mean by organized data in regular manner? Gyanendra Pokharel University of Guelph Guelph, ON On Tue, Nov 11, 2014 at 3:53 PM, David Winsemius dwinsem...@comcast.net wrote: On Nov 11, 2014, at 11:17 AM, Gyanendra Pokharel wrote: Hi R users, I am trying to plot two dimensional posterior likelihood surface. I have a data like para1 para2 likehood ... ... ... ... I looked at contour plot but it needs a shorted values of parameters and a matrix of likelihood values. Is there any way to get the plot? or how can I change my likelihood values to a matrix for the function contour? If the data are organized in a regular manner, then this might succeed: with( df, contour( x=unique(para1), y=unique(para2) z= matrix( likehood, length(unique(para1), length(unique(para2) ) ) ) Any suggestions are appreciated. GP University of Guelph Guelph, ON [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. David Winsemius Alameda, CA, USA [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Two dimensional likelihood surface plot
Thanks a lot, David ! Thats what I wanted. it is greatly helpful. GP University of Guelph Guelph, ON On Tue, Nov 11, 2014 at 7:02 PM, David Winsemius dwinsem...@comcast.net wrote: On Nov 11, 2014, at 2:53 PM, Gyanendra Pokharel wrote: Thanks David, what do you mean by organized data in regular manner? This is what I meant by a regular manner: matrix( c( rep( seq(0,1, by=.1), 11), rep( seq(0,1, by=.1),each=11) , runif(121) ), 121,3) [,1] [,2][,3] [1,] 0.0 0.0 0.048946906 [2,] 0.1 0.0 0.332529489 [3,] 0.2 0.0 0.967099462 [4,] 0.3 0.0 0.565349269 [5,] 0.4 0.0 0.024230243 [6,] 0.5 0.0 0.421633329 [7,] 0.6 0.0 0.965847357 [8,] 0.7 0.0 0.719618276 [9,] 0.8 0.0 0.948675911 [10,] 0.9 0.0 0.180241643 [11,] 1.0 0.0 0.804828468 [12,] 0.0 0.1 0.713698501 [13,] 0.1 0.1 0.991003966 [14,] 0.2 0.1 0.936413540 [15,] 0.3 0.1 0.941731063 [16,] 0.4 0.1 0.373998953 [17,] 0.5 0.1 0.988915380 [18,] 0.6 0.1 0.500791201 [19,] 0.7 0.1 0.070137099 [20,] 0.8 0.1 0.968422057 [21,] 0.9 0.1 0.827396746 snipped the remain 100 lines with( df, contour( x=unique(V1), y=unique(V2), z= matrix( V3, length(unique(V1)), length(unique(V2) ) ) )) df - as.data.frame( matrix( c( rep( seq(0,1, by=.1), 11), rep( seq(0,1, by=.1),each=11) , runif(121) ), 121,3)) str(df) 'data.frame': 121 obs. of 3 variables: $ V1: num 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 ... $ V2: num 0 0 0 0 0 0 0 0 0 0 ... $ V3: num 0.628 0.661 0.163 0.57 0.527 ... My original untested suggestion had some missing parentheses and a missing comma: with( df, contour( x=unique(V1), y=unique(V2), z= matrix( V3, length(unique(V1)), length(unique(V2) ) ) )) Gyanendra Pokharel University of Guelph Guelph, ON On Tue, Nov 11, 2014 at 3:53 PM, David Winsemius dwinsem...@comcast.net wrote: On Nov 11, 2014, at 11:17 AM, Gyanendra Pokharel wrote: Hi R users, I am trying to plot two dimensional posterior likelihood surface. I have a data like para1 para2 likehood ... ... ... ... I looked at contour plot but it needs a shorted values of parameters and a matrix of likelihood values. Is there any way to get the plot? or how can I change my likelihood values to a matrix for the function contour? If the data are organized in a regular manner, then this might succeed: with( df, contour( x=unique(para1), y=unique(para2) z= matrix( likehood, length(unique(para1), length(unique(para2) ) ) ) Any suggestions are appreciated. GP University of Guelph Guelph, ON [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. David Winsemius Alameda, CA, USA David Winsemius Alameda, CA, USA [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Help to install package mlegp
Hi R users, I have successfully downloaded the package mlegp, but when I tried installing it says the following massage. package ‘mlegp’ successfully unpacked and MD5 sums checked Warning: cannot remove prior installation of package ‘mlegp’ library(mlegp) Error in library(mlegp) : there is no package called ‘mlegp’ Can some one suggest me whats going on? Gyanendra Pokharel University of Guelph Guelph, ON [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Help to install package mlegp
Thanks Uwe! I re-install R and re-install the package mlegp and I successfully installed it. Thanks again ! Gyan Gyanendra Pokharel University of Guelph Guelph, ON On Thu, Aug 21, 2014 at 6:36 PM, Uwe Ligges lig...@statistik.tu-dortmund.de wrote: On 22.08.2014 00:15, Gyanendra Pokharel wrote: Hi R users, I have successfully downloaded the package mlegp, but when I tried installing it says the following massage. package ‘mlegp’ successfully unpacked and MD5 sums checked Warning: cannot remove prior installation of package ‘mlegp’ library(mlegp) Error in library(mlegp) : there is no package called ‘mlegp’ Can some one suggest me whats going on? Start a fresh R session and reinstall without loading it in advance. Best, Uwe Ligges Gyanendra Pokharel University of Guelph Guelph, ON [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/ posting-guide.html and provide commented, minimal, self-contained, reproducible code. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Package gptk, no default option
Hi there, Could someone explain the input options in the function gpCreat(q,d,X,y,options) in package gptk? I tried looking at the program gpOptions.R but did not get the way to input the option Many thanks Gyan .. Gyanendra Pokharel University of Guelph Guelph, ON [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Naming columns
Import data: mydata - read.table(filename.txt) give following command for column name colnames(mydata) - c(V1,V2,V3,) Gyanendra Pokharel University of Guelph Guelph, ON On Mon, Aug 26, 2013 at 5:42 PM, Docbanks84 mban...@partners.org wrote: Hi , I just imported a large data set from notepad. I want to label the columns in R. I used 'import Dataset' to bring in my data set Now, I would like to label V1,V2,V3 etc?? Thanks -- View this message in context: http://r.789695.n4.nabble.com/Naming-columns-tp4674595.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] error message solution: cannot allocate vector of size 200Mb
Try in R 64 bit. Thanks Gyanendra Pokharel University of Guelph Guelph, ON On Thu, May 23, 2013 at 10:53 AM, Ray Cheung ray1...@gmail.com wrote: Dear All, I wrote a program using R 2.15.2 but this error message cannot allocate vector of size 200Mb appeared. I want to ask in general how to handle this situation. I try to run the same program on other computers. It is perfectly fine. Can anybody help? Thank you very much in advance. Best Regards, Ray [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Netcdf file in R
I can't see your attached file. Can you re-attache it? Thanks Gyanendra Pokharel University of Guelph Guelph, ON On Fri, Feb 22, 2013 at 7:58 AM, Anup khanal za...@hotmail.com wrote: Good afternoon, I am a new in R. I have to work with large climate data.I am not able to read the netcdf file. Can anyone try with this file attached ? Best Regards, ..Anup KhanalNorwegian Institute of science and Technology (NTNU)Trondheim, NorwayMob:(+47) 45174313 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Import multiple data frames and combine them using cbind
I am sorry David, I don't mean to give the answer of the impossible questions. Just you can say impossible. If there is no way to do what I explained, that's fine, we do have other alternatives what I wrote in the beginning. Thank you so much. On Wed, Dec 5, 2012 at 1:12 PM, David Winsemius dwinsem...@comcast.netwrote: beneath [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Import multiple data frames and combine them using cbind
Thanks Berend, your Idea is great, that,s what I was looking. Thanks again On Wed, Dec 5, 2012 at 2:06 PM, Berend Hasselman b...@xs4all.nl wrote: do.call(cbind,lapply(mydf, function(df) df[,1:2])) [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Import multiple data frames and combine them using cbind
Hi group,
I imported 16 data frames using the function list.files
temp - list.files(path=...)
myfiles = lapply(temp, read.table,sep = )
Now I have 16 data set imported in R window.
I want to combine them by row and tried some thing like (Here I am
considering only 20 columns)
for(i in 1:16){
data- cbind(myfiles[[i]][,1:20])
}
but it returns only first data set. I can combine them using
data - cbind(myfiles[[1]][,1:20],myfiles[[2]][1:20],...)
But I want in a loop so that I can make the efficient code.
Any kind of suggestion will be great for me.
Thanks
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Import multiple data frames and combine them using cbind
Thanks Jim but I need to use cbind otherwise I will not get the correct
information from the data. Also each data set has the dimension 200X365 but
I only need 200X20 each.
Thanks
On Tue, Dec 4, 2012 at 10:22 PM, jim holtman jholt...@gmail.com wrote:
try::
do.call(rbind, myfiles)
On Tue, Dec 4, 2012 at 9:37 PM, Gyanendra Pokharel
gyanendra.pokha...@gmail.com wrote:
Hi group,
I imported 16 data frames using the function list.files
temp - list.files(path=...)
myfiles = lapply(temp, read.table,sep = )
Now I have 16 data set imported in R window.
I want to combine them by row and tried some thing like (Here I am
considering only 20 columns)
for(i in 1:16){
data- cbind(myfiles[[i]][,1:20])
}
but it returns only first data set. I can combine them using
data - cbind(myfiles[[1]][,1:20],myfiles[[2]][1:20],...)
But I want in a loop so that I can make the efficient code.
Any kind of suggestion will be great for me.
Thanks
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
--
Jim Holtman
Data Munger Guru
What is the problem that you are trying to solve?
Tell me what you want to do, not how you want to do it.
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Import multiple data frames and combine them using cbind
Thanks Dennis,
Your code also produces same as Jim's but I am not looking that one, I need
to use cbind so that finally I will get the data frame of size 200X320
(i,e, 200 X(16X20)).
Thanks
On Tue, Dec 4, 2012 at 10:46 PM, Dennis Murphy djmu...@gmail.com wrote:
In addition to Jim's reply, had you used the plyr package, you could
have done this in one shot:
library(plyr)
DF - ldply(llply(temp, read.table), rbind)
The inner llply call is equivalent to your lapply() and the outer
ldply() is equivalent to Jim's do.call() code.
Dennis
On Tue, Dec 4, 2012 at 6:37 PM, Gyanendra Pokharel
gyanendra.pokha...@gmail.com wrote:
Hi group,
I imported 16 data frames using the function list.files
temp - list.files(path=...)
myfiles = lapply(temp, read.table,sep = )
Now I have 16 data set imported in R window.
I want to combine them by row and tried some thing like (Here I am
considering only 20 columns)
for(i in 1:16){
data- cbind(myfiles[[i]][,1:20])
}
but it returns only first data set. I can combine them using
data - cbind(myfiles[[1]][,1:20],myfiles[[2]][1:20],...)
But I want in a loop so that I can make the efficient code.
Any kind of suggestion will be great for me.
Thanks
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] random forest
Hi all, Can some one tell me the difference between the following two formulas? 1. epiG.rf -randomForest(gamma~.,data=data, na.action = na.fail,ntree = 300,xtest = NULL, ytest = NULL,replace = T, proximity =F) 2.epiG.rf -randomForest(gamma~.,data=data, na.action = na.fail,ntree = 300,xtest = NULL, ytest = NULL,replace = T, proximity =F) [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] random forest
Sorry, the previous was not right post. I want to know the difference between following to methods of random forest. 1. epiG.rf -randomForest(gamma~.,data=data, na.action = na.fail,ntree = 300,xtest = NULL, ytest = NULL,replace = T, proximity =F) 2. epiG.rf -randomForest(x = data,,y = data$gamma, na.action = na.fail,ntree = 300,xtest = NULL, ytest = NULL,replace = T, proximity =F) which one is the correct form of random forest? Thanls On Sun, Oct 21, 2012 at 9:59 PM, Gyanendra Pokharel gyanendra.pokha...@gmail.com wrote: Hi all, Can some one tell me the difference between the following two formulas? 1. epiG.rf -randomForest(gamma~.,data=data, na.action = na.fail,ntree = 300,xtest = NULL, ytest = NULL,replace = T, proximity =F) 2.epiG.rf -randomForest(gamma~.,data=data, na.action = na.fail,ntree = 300,xtest = NULL, ytest = NULL,replace = T, proximity =F) [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Random Forest for multiple categorical variables
Dear all, I have the following data set. V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 alpha beta 1111 111 111 11111alpha beta1 2122 122 12 2 12212alpha beta1 3133 133 13 3 13 313alpha beta1 4144 14414 4 144 14 alpha beta1 5155 15515 5 155 15 alpha beta1 6166166 16 6 166 16 alpha beta2 717717717 7 17 7 17 alpha beta2 8188 18 818 818 818alpha beta2 919919919 9 19 9 19alpha beta2 10 20 10 20 10 20 10 20 10 20 alpha beta2 I want to use the randomForest classification. If there is one categorical variable with different classes, we can use randomForest(resp~., data, ), but here I need to classify the data with two categorical variables. Any idea will be great. Thanks [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Problem with nls regression fit
Hi all, I got following problem in fitting the data. Any kind of suggestions are welcome beta - 3.5 d - seq(0.1,62.5,0.1) y - exp(-beta*d) y1 - y x - read.table(epidist.txt, header = TRUE) data.nls - as.data.frame(cbind(y1,x)) #attach(data.nls) nls.fit - nls(y1~dist,data.nls) Error in cll[[1L]] : object of type 'symbol' is not subsettable Best [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Bayesian 95% Credible interval
Hi all, I have the data from the posterior distribution for some parameter. I want to find the 95% credible interval. I think t.test(data) is only for the confidence interval. I did not fine function for the Bayesian credible interval. Could some one suggest me? Thanks [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Fwd: Importing a CSV file
-- Forwarded message -- From: Gyanendra Pokharel gyanendra.pokha...@gmail.com Date: Fri, Feb 10, 2012 at 11:13 AM Subject: Re: [R] Importing a CSV file To: sezgin ozcan szgoz...@gmail.com You save your file in your own folder either in desktop or in document where ever you want and serach the file using the command mydatafile - read.csv(file.choose()) then open the file Cheers On Thu, Feb 9, 2012 at 10:20 PM, sezgin ozcan szgoz...@gmail.com wrote: I have been trying to import a csv file to r. but I get the same message everytime. the message is Error in file(file, rt) : cannot open the connection In addition: Warning message: In file(file, rt) : cannot open file 'Users:/sezginozcan/Downloads/beer.data.csv': No such file or directory I use mac. I tried this command also a-read.table(clipboard,sep=\t,row.names=1,header=T) Error: unexpected input in a-read.table(clipboard,sep= I will appreciate if you help me before I get crazy. thank you __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Problem in while loop
Yes Michael, it works well and I got the result what I want but it totally
depends on how reliable result do I want. When I take very high rho (near
about 1) and very low psi, it takes very long time may be it gives us more
accurate result. But for lower rho and higher psi, it gives immediately,
and very poor result. Now as I asked before, how I can get the plot of a
versus s which determines the cooling curve with the global minimum of s
at a specific value of a.
On Tue, Dec 6, 2011 at 6:52 AM, R. Michael Weylandt
michael.weyla...@gmail.com wrote:
So I just ran your code verbatim with this one change and it finished
in less than 10 seconds. However, even without the change it doesn't
take more than 15 seconds: what exactly lead you to believe this was
an infinite loop?
Michael
On Tue, Dec 6, 2011 at 12:03 AM, R. Michael Weylandt
michael.weyla...@gmail.com michael.weyla...@gmail.com wrote:
Off the bat I'd suggest you vectorize loglikelihood as a simple one
liner:
sum(log(b^2 + (x-a)^2))
That alone will speed up your function many times over: I'll look at the
big
function in more detail tomorrow.
Michael
On Dec 5, 2011, at 10:37 PM, Gyanendra Pokharel
gyanendra.pokha...@gmail.com wrote:
Thanks Michael
Lets figure out the problem by using the following function. I found the
same problem in this code too.
loglikehood - function(a, b = 0.1, x = c(-4.2, -2.85, -2.3, -1.02, 0.7,
0.98, 2.72, 3.5))
{
s - 0
for(i in 1:length(x)){
s - s + log(b^2 + (x[i] - a)^2)
}
s
}
loglikelihood(0.1)
simann - function(T0 = 1, N = 500, rho = 0.9, x0 = 0, eps = 0.1, f){
moving - 1
count - 0
Temp - T0
x - x0
while(moving 0){
moving - 0
for(i in 1:N){
y - x + runif(1,-eps,eps)
alpha - min(1,exp((f(x) -f(y))/Temp))
if(runif(1)alpha){
moving - moving +1
x - y
}
}
Temp - Temp*rho
count - count + 1
}
fmin - f(x)
return(c(x,fmin,count))
}
simann(f = loglikelihood)
I hope we can analyze every problems from this function
On Mon, Dec 5, 2011 at 10:27 PM, R. Michael Weylandt
michael.weyla...@gmail.com michael.weyla...@gmail.com wrote:
It's not necessarily equivalent to your loglikelihood function but
since
that function wasn't provided I couldn't test it.
My broader point is this: you said the problem was that the loop ran
endlessly: I showed it does not run endlessly for at least one input so
at
least part of the problem lies in loglikelihood, which, being
unprovided, I
have trouble replicating.
My original guess still stands: it's either 1) a case of you getting
stuck
at probaccept = 1 or 2) so slow it just feels endless. Either way, my
prescription is the same: print()
Michael
On Dec 5, 2011, at 9:30 PM, Gyanendra Pokharel
gyanendra.pokha...@gmail.com wrote:
Yes, your function out- epiann(f = function(a,b)
log(dnorm(a)*dnorm(b))),
N = 10) works well.
But why you are changing the loglikelihood function to f = function(a,b)
log(dnorm(a)*dnorm(b))? how it is equivalent to loglikelihood? is there
any
mathematical relation? I also want to see the plot of aout and bout
versus
loglikelihood, and see the cooling rate in different temperature. how is
this possible?
On Mon, Dec 5, 2011 at 6:07 PM, R. Michael Weylandt
michael.weyla...@gmail.com wrote:
If you run
out- epiann(f = function(a,b) log(dnorm(a)*dnorm(b))), N = 10)
It takes less than 0.5 seconds so there's no problem I can see:
perhaps you want to look elsewhere to get better speed (like Rcpp or
general vectorization), or maybe your loglikihood is not what's
desired, but there's no problem with the loop.
Michael
On Mon, Dec 5, 2011 at 5:29 PM, Gyanendra Pokharel
gyanendra.pokha...@gmail.com wrote:
Yes, I checked the acceptprob, it is very high but in my view, the
while
loop is not stopping, so there is some thing wrong in the use of
while
loop.
When I removed the while loop, it returned some thing but not the
result
what I want. When i run the while loop separately, it never stops.
On Mon, Dec 5, 2011 at 5:18 PM, R. Michael Weylandt
michael.weyla...@gmail.com wrote:
Your code is not reproducible nor minimal, but why don't you put a
command print(acceptprob) in and see if you are getting reasonable
values. If these values are extremely low it shouldn't surprise you
that your loop takes a long time to run.
More generally, read up on the use of print() and browser() as
debugging
tools.
Michael
On Mon, Dec 5, 2011 at 3:47 PM, Gyanendra Pokharel
gyanendra.pokha...@gmail.com wrote:
I forgot to upload the R-code in last email, so heare is one
epiann - function(T0 = 1, N=1000, ainit=1, binit=1,rho = 0.99,
amean =
3,
bmean=1.6, avar =.1, bvar=.1, f){
moving - 1
count - 0
Temp - T0
Re: [R] Problem in while loop
Thanks Michael,
I am able to find the very nice plot of the function we discussed but still
have the problem of ploting the function loglikelihood(aout,bout) versus
aout posted in the initial massage.
Best
On Tue, Dec 6, 2011 at 10:40 AM, Gyanendra Pokharel
gyanendra.pokha...@gmail.com wrote:
Yes Michael, it works well and I got the result what I want but it totally
depends on how reliable result do I want. When I take very high rho (near
about 1) and very low psi, it takes very long time may be it gives us more
accurate result. But for lower rho and higher psi, it gives immediately,
and very poor result. Now as I asked before, how I can get the plot of a
versus s which determines the cooling curve with the global minimum of s
at a specific value of a.
On Tue, Dec 6, 2011 at 6:52 AM, R. Michael Weylandt
michael.weyla...@gmail.com wrote:
So I just ran your code verbatim with this one change and it finished
in less than 10 seconds. However, even without the change it doesn't
take more than 15 seconds: what exactly lead you to believe this was
an infinite loop?
Michael
On Tue, Dec 6, 2011 at 12:03 AM, R. Michael Weylandt
michael.weyla...@gmail.com michael.weyla...@gmail.com wrote:
Off the bat I'd suggest you vectorize loglikelihood as a simple one
liner:
sum(log(b^2 + (x-a)^2))
That alone will speed up your function many times over: I'll look at
the big
function in more detail tomorrow.
Michael
On Dec 5, 2011, at 10:37 PM, Gyanendra Pokharel
gyanendra.pokha...@gmail.com wrote:
Thanks Michael
Lets figure out the problem by using the following function. I found the
same problem in this code too.
loglikehood - function(a, b = 0.1, x = c(-4.2, -2.85, -2.3, -1.02, 0.7,
0.98, 2.72, 3.5))
{
s - 0
for(i in 1:length(x)){
s - s + log(b^2 + (x[i] - a)^2)
}
s
}
loglikelihood(0.1)
simann - function(T0 = 1, N = 500, rho = 0.9, x0 = 0, eps = 0.1, f){
moving - 1
count - 0
Temp - T0
x - x0
while(moving 0){
moving - 0
for(i in 1:N){
y - x + runif(1,-eps,eps)
alpha - min(1,exp((f(x) -f(y))/Temp))
if(runif(1)alpha){
moving - moving +1
x - y
}
}
Temp - Temp*rho
count - count + 1
}
fmin - f(x)
return(c(x,fmin,count))
}
simann(f = loglikelihood)
I hope we can analyze every problems from this function
On Mon, Dec 5, 2011 at 10:27 PM, R. Michael Weylandt
michael.weyla...@gmail.com michael.weyla...@gmail.com wrote:
It's not necessarily equivalent to your loglikelihood function but
since
that function wasn't provided I couldn't test it.
My broader point is this: you said the problem was that the loop ran
endlessly: I showed it does not run endlessly for at least one input
so at
least part of the problem lies in loglikelihood, which, being
unprovided, I
have trouble replicating.
My original guess still stands: it's either 1) a case of you getting
stuck
at probaccept = 1 or 2) so slow it just feels endless. Either way, my
prescription is the same: print()
Michael
On Dec 5, 2011, at 9:30 PM, Gyanendra Pokharel
gyanendra.pokha...@gmail.com wrote:
Yes, your function out- epiann(f = function(a,b)
log(dnorm(a)*dnorm(b))),
N = 10) works well.
But why you are changing the loglikelihood function to f =
function(a,b)
log(dnorm(a)*dnorm(b))? how it is equivalent to loglikelihood? is
there any
mathematical relation? I also want to see the plot of aout and bout
versus
loglikelihood, and see the cooling rate in different temperature. how
is
this possible?
On Mon, Dec 5, 2011 at 6:07 PM, R. Michael Weylandt
michael.weyla...@gmail.com wrote:
If you run
out- epiann(f = function(a,b) log(dnorm(a)*dnorm(b))), N = 10)
It takes less than 0.5 seconds so there's no problem I can see:
perhaps you want to look elsewhere to get better speed (like Rcpp or
general vectorization), or maybe your loglikihood is not what's
desired, but there's no problem with the loop.
Michael
On Mon, Dec 5, 2011 at 5:29 PM, Gyanendra Pokharel
gyanendra.pokha...@gmail.com wrote:
Yes, I checked the acceptprob, it is very high but in my view, the
while
loop is not stopping, so there is some thing wrong in the use of
while
loop.
When I removed the while loop, it returned some thing but not the
result
what I want. When i run the while loop separately, it never stops.
On Mon, Dec 5, 2011 at 5:18 PM, R. Michael Weylandt
michael.weyla...@gmail.com wrote:
Your code is not reproducible nor minimal, but why don't you put a
command print(acceptprob) in and see if you are getting reasonable
values. If these values are extremely low it shouldn't surprise you
that your loop takes a long time to run.
More generally, read up on the use of print() and browser() as
debugging
tools.
Michael
On Mon, Dec 5, 2011 at 3:47 PM, Gyanendra
[R] Problem in while loop
Hi all, I have the following code, When I run the code, it never terminate this is because of the while loop i am using. In general, if you need a loop for which you don't know in advance how many iterations there will be, you can use the `while' statement so here too i don't know the number how many iterations are there. So Can some one suggest me whats going on? I am using the Metropolis simulated annealing algorithm Best [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Problem in while loop
I forgot to upload the R-code in last email, so heare is one
epiann - function(T0 = 1, N=1000, ainit=1, binit=1,rho = 0.99, amean = 3,
bmean=1.6, avar =.1, bvar=.1, f){
moving - 1
count - 0
Temp - T0
aout - ainit
bout - binit
while(moving 0){
moving - 0
for (i in 1:N) {
aprop - rnorm (1,amean, avar)
bprop - rnorm (1,bmean, bvar)
if (aprop 0 bprop 0){
acceptprob - min(1,exp((f(aout, bout) -
f(aprop,bprop))/Temp))
u - runif(1)
if(uacceptprob){
moving - moving +1
aout - aprop
bout - bprop
}
else{aprob - aout
bprob - bout}
}
}
Temp - Temp*rho
count - count +1
}
fmin - f(aout,bout)
return(c(aout, bout,fmin, count) )
}
out- epiann(f = loglikelihood)
On Mon, Dec 5, 2011 at 3:46 PM, Gyanendra Pokharel
gyanendra.pokha...@gmail.com wrote:
Hi all,
I have the following code,
When I run the code, it never terminate this is because of the while loop
i am using. In general, if you need a loop for which you don't know in
advance how many iterations there will be, you can use the `while'
statement so here too i don't know the number how many iterations are
there. So Can some one suggest me whats going on?
I am using the Metropolis simulated annealing algorithm
Best
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Problem in while loop
Yes, I checked the acceptprob, it is very high but in my view, the while
loop is not stopping, so there is some thing wrong in the use of while
loop. When I removed the while loop, it returned some thing but not the
result what I want. When i run the while loop separately, it never stops.
On Mon, Dec 5, 2011 at 5:18 PM, R. Michael Weylandt
michael.weyla...@gmail.com wrote:
Your code is not reproducible nor minimal, but why don't you put a
command print(acceptprob) in and see if you are getting reasonable
values. If these values are extremely low it shouldn't surprise you
that your loop takes a long time to run.
More generally, read up on the use of print() and browser() as debugging
tools.
Michael
On Mon, Dec 5, 2011 at 3:47 PM, Gyanendra Pokharel
gyanendra.pokha...@gmail.com wrote:
I forgot to upload the R-code in last email, so heare is one
epiann - function(T0 = 1, N=1000, ainit=1, binit=1,rho = 0.99, amean =
3,
bmean=1.6, avar =.1, bvar=.1, f){
moving - 1
count - 0
Temp - T0
aout - ainit
bout - binit
while(moving 0){
moving - 0
for (i in 1:N) {
aprop - rnorm (1,amean, avar)
bprop - rnorm (1,bmean, bvar)
if (aprop 0 bprop 0){
acceptprob - min(1,exp((f(aout, bout) -
f(aprop,bprop))/Temp))
u - runif(1)
if(uacceptprob){
moving - moving +1
aout - aprop
bout - bprop
}
else{aprob - aout
bprob - bout}
}
}
Temp - Temp*rho
count - count +1
}
fmin - f(aout,bout)
return(c(aout, bout,fmin, count) )
}
out- epiann(f = loglikelihood)
On Mon, Dec 5, 2011 at 3:46 PM, Gyanendra Pokharel
gyanendra.pokha...@gmail.com wrote:
Hi all,
I have the following code,
When I run the code, it never terminate this is because of the while
loop
i am using. In general, if you need a loop for which you don't know in
advance how many iterations there will be, you can use the `while'
statement so here too i don't know the number how many iterations are
there. So Can some one suggest me whats going on?
I am using the Metropolis simulated annealing algorithm
Best
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Problem in while loop
Yes, your function out- epiann(f = function(a,b) log(dnorm(a)*dnorm(b))),
N = 10) works well.
But why you are changing the loglikelihood function to f = function(a,b)
log(dnorm(a)*dnorm(b))? how it is equivalent to loglikelihood? is there any
mathematical relation? I also want to see the plot of aout and bout versus
loglikelihood, and see the cooling rate in different temperature. how is
this possible?
On Mon, Dec 5, 2011 at 6:07 PM, R. Michael Weylandt
michael.weyla...@gmail.com wrote:
If you run
out- epiann(f = function(a,b) log(dnorm(a)*dnorm(b))), N = 10)
It takes less than 0.5 seconds so there's no problem I can see:
perhaps you want to look elsewhere to get better speed (like Rcpp or
general vectorization), or maybe your loglikihood is not what's
desired, but there's no problem with the loop.
Michael
On Mon, Dec 5, 2011 at 5:29 PM, Gyanendra Pokharel
gyanendra.pokha...@gmail.com wrote:
Yes, I checked the acceptprob, it is very high but in my view, the while
loop is not stopping, so there is some thing wrong in the use of while
loop.
When I removed the while loop, it returned some thing but not the result
what I want. When i run the while loop separately, it never stops.
On Mon, Dec 5, 2011 at 5:18 PM, R. Michael Weylandt
michael.weyla...@gmail.com wrote:
Your code is not reproducible nor minimal, but why don't you put a
command print(acceptprob) in and see if you are getting reasonable
values. If these values are extremely low it shouldn't surprise you
that your loop takes a long time to run.
More generally, read up on the use of print() and browser() as debugging
tools.
Michael
On Mon, Dec 5, 2011 at 3:47 PM, Gyanendra Pokharel
gyanendra.pokha...@gmail.com wrote:
I forgot to upload the R-code in last email, so heare is one
epiann - function(T0 = 1, N=1000, ainit=1, binit=1,rho = 0.99, amean
=
3,
bmean=1.6, avar =.1, bvar=.1, f){
moving - 1
count - 0
Temp - T0
aout - ainit
bout - binit
while(moving 0){
moving - 0
for (i in 1:N) {
aprop - rnorm (1,amean, avar)
bprop - rnorm (1,bmean, bvar)
if (aprop 0 bprop 0){
acceptprob - min(1,exp((f(aout, bout) -
f(aprop,bprop))/Temp))
u - runif(1)
if(uacceptprob){
moving - moving +1
aout - aprop
bout - bprop
}
else{aprob - aout
bprob - bout}
}
}
Temp - Temp*rho
count - count +1
}
fmin - f(aout,bout)
return(c(aout, bout,fmin, count) )
}
out- epiann(f = loglikelihood)
On Mon, Dec 5, 2011 at 3:46 PM, Gyanendra Pokharel
gyanendra.pokha...@gmail.com wrote:
Hi all,
I have the following code,
When I run the code, it never terminate this is because of the while
loop
i am using. In general, if you need a loop for which you don't know
in
advance how many iterations there will be, you can use the `while'
statement so here too i don't know the number how many iterations are
there. So Can some one suggest me whats going on?
I am using the Metropolis simulated annealing algorithm
Best
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.htmlhttp://www.r-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Problem in while loop
Thanks Michael
Lets figure out the problem by using the following function. I found the
same problem in this code too.
loglikehood - function(a, b = 0.1, x = c(-4.2, -2.85, -2.3, -1.02, 0.7,
0.98, 2.72, 3.5))
{
s - 0
for(i in 1:length(x)){
s - s + log(b^2 + (x[i] - a)^2)
}
s
}
loglikelihood(0.1)
simann - function(T0 = 1, N = 500, rho = 0.9, x0 = 0, eps = 0.1, f){
moving - 1
count - 0
Temp - T0
x - x0
while(moving 0){
moving - 0
for(i in 1:N){
y - x + runif(1,-eps,eps)
alpha - min(1,exp((f(x) -f(y))/Temp))
if(runif(1)alpha){
moving - moving +1
x - y
}
}
Temp - Temp*rho
count - count + 1
}
fmin - f(x)
return(c(x,fmin,count))
}
simann(f = loglikelihood)
I hope we can analyze every problems from this function
On Mon, Dec 5, 2011 at 10:27 PM, R. Michael Weylandt
michael.weyla...@gmail.com michael.weyla...@gmail.com wrote:
It's not necessarily equivalent to your loglikelihood function but
since that function wasn't provided I couldn't test it.
My broader point is this: you said the problem was that the loop ran
endlessly: I showed it does not run endlessly for at least one input so at
least part of the problem lies in loglikelihood, which, being unprovided, I
have trouble replicating.
My original guess still stands: it's either 1) a case of you getting stuck
at probaccept = 1 or 2) so slow it just feels endless. Either way, my
prescription is the same: print()
Michael
On Dec 5, 2011, at 9:30 PM, Gyanendra Pokharel
gyanendra.pokha...@gmail.com wrote:
Yes, your function out- epiann(f = function(a,b)
log(dnorm(a)*dnorm(b))), N = 10) works well.
But why you are changing the loglikelihood function to f = function(a,b)
log(dnorm(a)*dnorm(b))? how it is equivalent to loglikelihood? is there any
mathematical relation? I also want to see the plot of aout and bout versus
loglikelihood, and see the cooling rate in different temperature. how is
this possible?
On Mon, Dec 5, 2011 at 6:07 PM, R. Michael Weylandt
michael.weyla...@gmail.com wrote:
If you run
out- epiann(f = function(a,b) log(dnorm(a)*dnorm(b))), N = 10)
It takes less than 0.5 seconds so there's no problem I can see:
perhaps you want to look elsewhere to get better speed (like Rcpp or
general vectorization), or maybe your loglikihood is not what's
desired, but there's no problem with the loop.
Michael
On Mon, Dec 5, 2011 at 5:29 PM, Gyanendra Pokharel
gyanendra.pokha...@gmail.com wrote:
Yes, I checked the acceptprob, it is very high but in my view, the while
loop is not stopping, so there is some thing wrong in the use of while
loop.
When I removed the while loop, it returned some thing but not the result
what I want. When i run the while loop separately, it never stops.
On Mon, Dec 5, 2011 at 5:18 PM, R. Michael Weylandt
michael.weyla...@gmail.com wrote:
Your code is not reproducible nor minimal, but why don't you put a
command print(acceptprob) in and see if you are getting reasonable
values. If these values are extremely low it shouldn't surprise you
that your loop takes a long time to run.
More generally, read up on the use of print() and browser() as
debugging
tools.
Michael
On Mon, Dec 5, 2011 at 3:47 PM, Gyanendra Pokharel
gyanendra.pokha...@gmail.com wrote:
I forgot to upload the R-code in last email, so heare is one
epiann - function(T0 = 1, N=1000, ainit=1, binit=1,rho = 0.99,
amean =
3,
bmean=1.6, avar =.1, bvar=.1, f){
moving - 1
count - 0
Temp - T0
aout - ainit
bout - binit
while(moving 0){
moving - 0
for (i in 1:N) {
aprop - rnorm (1,amean, avar)
bprop - rnorm (1,bmean, bvar)
if (aprop 0 bprop 0){
acceptprob - min(1,exp((f(aout, bout) -
f(aprop,bprop))/Temp))
u - runif(1)
if(uacceptprob){
moving - moving +1
aout - aprop
bout - bprop
}
else{aprob - aout
bprob - bout}
}
}
Temp - Temp*rho
count - count +1
}
fmin - f(aout,bout)
return(c(aout, bout,fmin, count) )
}
out- epiann(f = loglikelihood)
On Mon, Dec 5, 2011 at 3:46 PM, Gyanendra Pokharel
gyanendra.pokha...@gmail.com wrote:
Hi all,
I have the following code,
When I run the code, it never terminate this is because of the while
loop
i am using. In general, if you need a loop for which you don't know
in
advance how many iterations there will be, you can use the `while'
statement so here too i don't know the number how many iterations
are
there. So Can some one suggest me whats going on?
I am using the Metropolis simulated annealing algorithm
Best
[[alternative HTML
Re: [R] Problem in while loop
I ma sorry, I miss typed the function, it should be loglikelihood instead.
On Mon, Dec 5, 2011 at 10:37 PM, Gyanendra Pokharel
gyanendra.pokha...@gmail.com wrote:
Thanks Michael
Lets figure out the problem by using the following function. I found the
same problem in this code too.
loglikehood - function(a, b = 0.1, x = c(-4.2, -2.85, -2.3, -1.02, 0.7,
0.98, 2.72, 3.5))
{
s - 0
for(i in 1:length(x)){
s - s + log(b^2 + (x[i] - a)^2)
}
s
}
loglikelihood(0.1)
simann - function(T0 = 1, N = 500, rho = 0.9, x0 = 0, eps = 0.1, f){
moving - 1
count - 0
Temp - T0
x - x0
while(moving 0){
moving - 0
for(i in 1:N){
y - x + runif(1,-eps,eps)
alpha - min(1,exp((f(x) -f(y))/Temp))
if(runif(1)alpha){
moving - moving +1
x - y
}
}
Temp - Temp*rho
count - count + 1
}
fmin - f(x)
return(c(x,fmin,count))
}
simann(f = loglikelihood)
I hope we can analyze every problems from this function
On Mon, Dec 5, 2011 at 10:27 PM, R. Michael Weylandt
michael.weyla...@gmail.com michael.weyla...@gmail.com wrote:
It's not necessarily equivalent to your loglikelihood function but
since that function wasn't provided I couldn't test it.
My broader point is this: you said the problem was that the loop ran
endlessly: I showed it does not run endlessly for at least one input so at
least part of the problem lies in loglikelihood, which, being unprovided, I
have trouble replicating.
My original guess still stands: it's either 1) a case of you getting
stuck at probaccept = 1 or 2) so slow it just feels endless. Either way,
my prescription is the same: print()
Michael
On Dec 5, 2011, at 9:30 PM, Gyanendra Pokharel
gyanendra.pokha...@gmail.com wrote:
Yes, your function out- epiann(f = function(a,b)
log(dnorm(a)*dnorm(b))), N = 10) works well.
But why you are changing the loglikelihood function to f = function(a,b)
log(dnorm(a)*dnorm(b))? how it is equivalent to loglikelihood? is there any
mathematical relation? I also want to see the plot of aout and bout versus
loglikelihood, and see the cooling rate in different temperature. how is
this possible?
On Mon, Dec 5, 2011 at 6:07 PM, R. Michael Weylandt
michael.weyla...@gmail.com wrote:
If you run
out- epiann(f = function(a,b) log(dnorm(a)*dnorm(b))), N = 10)
It takes less than 0.5 seconds so there's no problem I can see:
perhaps you want to look elsewhere to get better speed (like Rcpp or
general vectorization), or maybe your loglikihood is not what's
desired, but there's no problem with the loop.
Michael
On Mon, Dec 5, 2011 at 5:29 PM, Gyanendra Pokharel
gyanendra.pokha...@gmail.com wrote:
Yes, I checked the acceptprob, it is very high but in my view, the
while
loop is not stopping, so there is some thing wrong in the use of while
loop.
When I removed the while loop, it returned some thing but not the
result
what I want. When i run the while loop separately, it never stops.
On Mon, Dec 5, 2011 at 5:18 PM, R. Michael Weylandt
michael.weyla...@gmail.com wrote:
Your code is not reproducible nor minimal, but why don't you put a
command print(acceptprob) in and see if you are getting reasonable
values. If these values are extremely low it shouldn't surprise you
that your loop takes a long time to run.
More generally, read up on the use of print() and browser() as
debugging
tools.
Michael
On Mon, Dec 5, 2011 at 3:47 PM, Gyanendra Pokharel
gyanendra.pokha...@gmail.com wrote:
I forgot to upload the R-code in last email, so heare is one
epiann - function(T0 = 1, N=1000, ainit=1, binit=1,rho = 0.99,
amean =
3,
bmean=1.6, avar =.1, bvar=.1, f){
moving - 1
count - 0
Temp - T0
aout - ainit
bout - binit
while(moving 0){
moving - 0
for (i in 1:N) {
aprop - rnorm (1,amean, avar)
bprop - rnorm (1,bmean, bvar)
if (aprop 0 bprop 0){
acceptprob - min(1,exp((f(aout, bout) -
f(aprop,bprop))/Temp))
u - runif(1)
if(uacceptprob){
moving - moving +1
aout - aprop
bout - bprop
}
else{aprob - aout
bprob - bout}
}
}
Temp - Temp*rho
count - count +1
}
fmin - f(aout,bout)
return(c(aout, bout,fmin, count) )
}
out- epiann(f = loglikelihood)
On Mon, Dec 5, 2011 at 3:46 PM, Gyanendra Pokharel
gyanendra.pokha...@gmail.com wrote:
Hi all,
I have the following code,
When I run the code, it never terminate this is because of the
while
loop
i am using. In general, if you need a loop for which you don't
know in
advance how many iterations there will be, you can use the `while'
statement so here
Re: [R] Problem in log
The loglikelihood() looks ok and gives some value. But I am using this
function for the simulated annealing, generating the random samples from
uniform proposal density., The codes are as follows
epiannea - function(T0 = 1, N = 500,beta = 0.1,x0 = 0.1, rho = 0.90, eps =
0.1, loglikelihood, data){
moving - 1
count - 0
Temp - T0
alpha - x0
while(moving 0){
mmoving - 0
for(i in 1:N){
r - alpha + runif(1,-eps,eps)
if(r 0){
a - min(1,exp((loglikelihood(alpha,beta)-loglikelihood(r,beta))/Temp))
}
if(runif(1) a){
moving - moving +1
alpha - r
}
}
Temp - Temp*rho
count - count + 1
}
#plot(a,loglikelihood(alpha,beta), type = l)
fmin - loglikelihood(alpha, beta)
return(c(fmin, alpha, count))
}
epiannea(loglikelihood=loglikelihood)
Some times it shows the warnings that logp[i] produces NaNs, that means
p[i] is negative, but it should not be that because p[i] is the probabiliy
and can't be negative. Some times the code runs but never stop.
On Wed, Nov 30, 2011 at 7:34 AM, R. Michael Weylandt
michael.weyla...@gmail.com michael.weyla...@gmail.com wrote:
I'd suggest you do some leg-work and figure out why you are getting values
1. If your algorithm is motivated by some approximation then a min() or
pmin() *might* be the right fix, but if there are no approximations you may
need to start debugging properly to see why you are getting an out of
bounds value.
Since there's no random number generation involved, I'd hesitate to just
throw out the result without knowing its source. Also keep in mind the
limitations of floating point arithmetic if you expect alpha*d^beta to be
small.
Michael
On Nov 29, 2011, at 6:58 PM, Sarah Goslee sarah.gos...@gmail.com wrote:
On Tue, Nov 29, 2011 at 6:55 PM, Gyanendra Pokharel
gyanendra.pokha...@gmail.com wrote:
yes, log of negative number is undefined and R also do the same and
produces
NaNs. Here I want to reject the value of exp(-alpha*d^(-beta)) when
greater
than 1, and want to run the loop otherwise.
Thanks
Then just add another if() statement checking for that condition.
On Tue, Nov 29, 2011 at 6:48 PM, Sarah Goslee sarah.gos...@gmail.com
wrote:
Here p[i] - 1 - exp(-alpha*d^(-beta)) so, log(p[i]) produces NaNs
when exp(-alpha*d^(-beta)) is greater than 1. How can I remove
it.After
generating the out put we can omit it, but the problem is different.
Wait... you're complaining that you can't take the natural log of a
negative
number in R?
You can't do that anywhere. What do you expect to happen? The log of a
negative number IS NaN.
Sarah
On Tue, Nov 29, 2011 at 6:28 PM, Gyanendra Pokharel
gyanendra.pokha...@gmail.com wrote:
I have following code:
loglikelihood - function(alpha,beta= 0.1){
loglh-0
d-0
p-0
k-NULL
data-read.table(epidemic.txt,header = TRUE)
attach(data, warn.conflicts = F)
k -which(inftime==1)
d - (sqrt((x-x[k])^2+(y-y[k])^2))^(-beta)
p-1 - exp(-alpha*d)
for(i in 1:100){
if(i!=k){
if(inftime[i]==0){
loglh-loglh +log(1-p[i])
}
if(inftime[i]==2){
loglh-loglh + log(p[i])
}
}
}
return(loglh)
}
Here p[i] - 1 - exp(-alpha*d^(-beta))
so, log(p[i]) produces NaNs when exp(-alpha*d^(-beta)) is greater
than
1.
How can I remove it.After generating the out put we can omit it, but
the
problem is different.
On Tue, Nov 29, 2011 at 5:22 PM, Gyanendra Pokharel
gyanendra.pokha...@gmail.com wrote:
No, that,s not a problem Michael,
I have following code:
loglikelihood - function(alpha,beta= 0.1){
loglh-0
d-0
p-0
k-NULL
data-read.table(epidemic.txt,header = TRUE)
attach(data, warn.conflicts = F)
k -which(inftime==1)
d - (sqrt((x-x[k])^2+(y-y[k])^2))^(-beta)
p-1 - exp(-alpha*d)
for(i in 1:100){
if(i!=k){
if(inftime[i]==0){
loglh-loglh +log(1-p[i])
}
if(inftime[i]==2){
loglh-loglh + log(p[i])
}
}
}
return(loglh)
}
Here p[i] - 1 - exp(-alpha*d^(-beta))
so, log(p[i]) produces NaNs when exp(-alpha*d^(-beta)) is greater
than
1.
How can I remove it.After generating the out put we can omit it, but
the
problem is different.
--
Sarah Goslee
http://www.functionaldiversity.org
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.htmlhttp://www.r-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
[R] Problem in log
Hi all I have a function of log defined by y = log(1- exp(-a)), when exp(-a) is greater, 1, it produce NaN. How can I remove this in R? [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Problem in log
I have following code:
loglikelihood - function(alpha,beta= 0.1){
loglh-0
d-0
p-0
k-NULL
data-read.table(epidemic.txt,header = TRUE)
attach(data, warn.conflicts = F)
k -which(inftime==1)
d - (sqrt((x-x[k])^2+(y-y[k])^2))^(-beta)
p-1 - exp(-alpha*d)
for(i in 1:100){
if(i!=k){
if(inftime[i]==0){
loglh-loglh +log(1-p[i])
}
if(inftime[i]==2){
loglh-loglh + log(p[i])
}
}
}
return(loglh)
}
Here p[i] - 1 - exp(-alpha*d^(-beta))
so, log(p[i]) produces NaNs when exp(-alpha*d^(-beta)) is greater than 1.
How can I remove it.After generating the out put we can omit it, but the
problem is different.
On Tue, Nov 29, 2011 at 5:22 PM, Gyanendra Pokharel
gyanendra.pokha...@gmail.com wrote:
No, that,s not a problem Michael,
I have following code:
loglikelihood - function(alpha,beta= 0.1){
loglh-0
d-0
p-0
k-NULL
data-read.table(epidemic.txt,header = TRUE)
attach(data, warn.conflicts = F)
k -which(inftime==1)
d - (sqrt((x-x[k])^2+(y-y[k])^2))^(-beta)
p-1 - exp(-alpha*d)
for(i in 1:100){
if(i!=k){
if(inftime[i]==0){
loglh-loglh +log(1-p[i])
}
if(inftime[i]==2){
loglh-loglh + log(p[i])
}
}
}
return(loglh)
}
Here p[i] - 1 - exp(-alpha*d^(-beta))
so, log(p[i]) produces NaNs when exp(-alpha*d^(-beta)) is greater than 1.
How can I remove it.After generating the out put we can omit it, but the
problem is different.
On Tue, Nov 29, 2011 at 5:07 PM, R. Michael Weylandt
michael.weyla...@gmail.com wrote:
Do you mean remove the NaNs? Try na.omit() or complete.cases() or many
other options.
If you mean you want the complex log, try
log(as.complex(1-exp(-a)))
Michael
On Tue, Nov 29, 2011 at 5:02 PM, Gyanendra Pokharel
gyanendra.pokha...@gmail.com wrote:
Hi all I have a function of log defined by y = log(1- exp(-a)), when
exp(-a) is greater, 1, it produce NaN. How can I remove this in R?
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.htmlhttp://www.r-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Problem in log
yes, log of negative number is undefined and R also do the same and
produces NaNs. Here I want to reject the value of exp(-alpha*d^(-beta))
when greater than 1, and want to run the loop otherwise.
Thanks
On Tue, Nov 29, 2011 at 6:48 PM, Sarah Goslee sarah.gos...@gmail.comwrote:
Here p[i] - 1 - exp(-alpha*d^(-beta)) so, log(p[i]) produces NaNs
when exp(-alpha*d^(-beta)) is greater than 1. How can I remove it.After
generating the out put we can omit it, but the problem is different.
Wait... you're complaining that you can't take the natural log of a
negative
number in R?
You can't do that anywhere. What do you expect to happen? The log of a
negative number IS NaN.
Sarah
On Tue, Nov 29, 2011 at 6:28 PM, Gyanendra Pokharel
gyanendra.pokha...@gmail.com wrote:
I have following code:
loglikelihood - function(alpha,beta= 0.1){
loglh-0
d-0
p-0
k-NULL
data-read.table(epidemic.txt,header = TRUE)
attach(data, warn.conflicts = F)
k -which(inftime==1)
d - (sqrt((x-x[k])^2+(y-y[k])^2))^(-beta)
p-1 - exp(-alpha*d)
for(i in 1:100){
if(i!=k){
if(inftime[i]==0){
loglh-loglh +log(1-p[i])
}
if(inftime[i]==2){
loglh-loglh + log(p[i])
}
}
}
return(loglh)
}
Here p[i] - 1 - exp(-alpha*d^(-beta))
so, log(p[i]) produces NaNs when exp(-alpha*d^(-beta)) is greater than
1.
How can I remove it.After generating the out put we can omit it, but the
problem is different.
On Tue, Nov 29, 2011 at 5:22 PM, Gyanendra Pokharel
gyanendra.pokha...@gmail.com wrote:
No, that,s not a problem Michael,
I have following code:
loglikelihood - function(alpha,beta= 0.1){
loglh-0
d-0
p-0
k-NULL
data-read.table(epidemic.txt,header = TRUE)
attach(data, warn.conflicts = F)
k -which(inftime==1)
d - (sqrt((x-x[k])^2+(y-y[k])^2))^(-beta)
p-1 - exp(-alpha*d)
for(i in 1:100){
if(i!=k){
if(inftime[i]==0){
loglh-loglh +log(1-p[i])
}
if(inftime[i]==2){
loglh-loglh + log(p[i])
}
}
}
return(loglh)
}
Here p[i] - 1 - exp(-alpha*d^(-beta))
so, log(p[i]) produces NaNs when exp(-alpha*d^(-beta)) is greater than
1.
How can I remove it.After generating the out put we can omit it, but the
problem is different.
--
Sarah Goslee
http://www.functionaldiversity.org
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] Dataset fit for Wavelet regression
Hi all, Can any body suggest me the data set which has the length of power of two and fit for the wavelet smoothing for the regression? There are alot of datasets like ethanol in the package SemiPar that can be used for this smoothing but we have to extend them in the power of two. I have no idea of this. If any body give the technique for this, it will be great. Best [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] number of items to replace is not a multiple of replacement length
Hi all, following is my R -code and shows the error given below
n - 100
k-2
x1 -c(1, 3);x2 - c(2,5)
X - matrix(c(0,0), nrow = 2, ncol = n)
for(i in 1:k)
+ X[i, ] - mh1.epidemic(n,x1[i],x2[i])
Error in X[i, ] - mh1.epidemic(n,x1[i],x2[i]):
number of items to replace is not a multiple of replacement length
psi -t(apply(X, c(1,2), cumsum))
for(i in 1:nrow(psi)){
+ psi[i,] psi[i,]/(1:ncol(psi))
+ print(Gelman.Rubin(psi))
+
+ }
Can some one suggest me what this error means and how could I fix it?
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] Permutations
Hi all, why factorial(150) shows the error out of range in 'gammafn'? I have to calculate the number of subset formed by 150 samples taking 10 at a time. How is this possible? best [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Permutation matrix
Hi all, I have a set of elements (1, 1, 0,1,1,0,1,0,1,1) with ten elements. I have to construct the permutation matrix of this set which is of the size 10 by 2^10. Can some one help how is this possible? Is there is a particular function in R or I need to make function? Best [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Error in random walk Metroplis-hasting
Hi R community,
I have some data set and construct the likelihood as follows
likelihood - function(alpha,beta){
lh-1
d-0
p-0
k-NULL
data-read.table(epidemic.txt,header = TRUE)
attach(data, warn.conflicts = F)
k -which(inftime==1)
d - (sqrt((x-x[k])^2+(y-y[k])^2))^(-beta)
p-1 - exp(-alpha*d)
for(i in 1:100){
if(i!=k){
if(inftime[i]==0){
lh-lh*(1-p[i])
}
if(inftime[i]==2){
lh-lh*p[i]
}
}
}
return(lh)
}
Then, I want to simulate this by using random walk Metropolis- Hasting
algorithm with the single parameter update. i have the following R function
mh.epidemic -function(m,alpha, beta){
z- array(0,c(m+1, 2))
z[1,1] - alpha
z[1,2] - beta
for(t in 2:m){
u - runif(1)
y1 - rexp(1,z[t-1,1])
y2 -rexp(1,z[t-1,2])
z[t,1] - likelihood(y1,z[t-1,2])/likelihood(z[t-1,1],z[t-1,2])
a1 -min(1,z)
if(u=a1) z[t,] - y1 else
z[t,2] -z[t-1,1]
z[t,2]-likelihood(z[t,1], y2)/likelihood(z[t,1],z[t-1,2])
accept -min(1,z)
if(u=accept) z[t,] - y2 else
z[t,2] -z[t-1,1]
}
z
}
mh.epidemic(100, 1,2)
when I run it shows the following error:
Error in if (u = accept) z[t, ] - y2 else z[t, 2] - z[t - 1, 1] :
missing value where TRUE/FALSE needed
I know it is due to the NaN produce some where. So I tried by running the
code separately, as follows
m- 100
alpha -1
beta - 2
z- array(0,c(m+1, 2))
z[1,1] - alpha
z[1,2] - beta
for(t in 2:m){
+ u - runif(1)
+ y1 - rexp(1,z[t-1,1])
+ y2 -rexp(1,z[t-1,2])
+ z[t,1] - likelihood(y1,z[t-1,2])/likelihood(z[t-1,1],z[t-1,2])
+ a1 -min(1,z)
+ }
There were 50 or more warnings (use warnings() to see the first 50)
y1
[1] NaN
y2
[1] NaN
why, this simulation from exponential function is produce NaN?
If some one help me it will be great.
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] Help
Hi all, I have the mean vector mu- c(0,0) and variance sigma - c(10,10), now how to sample from the bivariate normal density in R? Can some one suggest me? I did not fine the function mvdnorm in R. Best Gyan [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Random-walk Metropolis-Hasting
Following is my code, can some one help on the error at the bottom?
mh-function(iterations,alpha,beta){
+ data-read.table(epidemic.txt,header = TRUE)
+ attach(data, warn.conflicts = F)
+ k-97
+ d - (sqrt((x-x[k])^2 + (y-y[k])^2))
+ p - 1-exp(-alpha*d^(-beta))
+ p.alpha-1 - exp(-3*d^(-beta))
+ p.beta - 1 - exp(alpha*d^(-2))
+ iterations-1000
+ mu.lambda - c(0,0);s.lambda - c(100,100)
+ prop.s - c(0.1,0.1)
+ lambda - matrix(nrow=iterations, ncol=2)
+ acc.prob-0
+ current.lambda - c(0,0)
+ for(t in 1:iterations){
+ prop.lambda - rnorm(2,current.lambda,prop.s)
+ a - p.beta/p.alpha
*(dnorm(prop.lambda,mu.lambda,s.lambda))*dnorm(current.lambda,mu.lambda,s.lambda)
+ accept - min(1,a)
+ u-runif(1)
+ if(u[t]=accept[t]){
+ current.lambda - prop.lambda
+ acc.prob - acc.prob +1
+ }
+ lambda[t,] - current.lambda
+ }
+ lambda
+ }
mh(1000,0,0)
Error in if (u[t] = accept[t]) { : missing value where TRUE/FALSE needed
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] Gibbs sampler
I have the following code,
gibbs -function(m,theta = 0.25, lambda =0.55, n =1){
alpha - 1.5
beta - 1.5
gamma - 1.5
x- array(0,c(m+1, 3))
x[1,1] - theta
x[1,2] - lambda
x[1,3]- n
for(t in 2:(m+1)){
x[t,1] - rbinom(1, x[t-1,3], x[t-1,1])
x[t,2]-rbeta(1, x[t-1,1] + alpha, x[t-1,3] - x[t-1,1] + beta)
x[t,3] - rpois(1,(1 - x[t-1,1])*gamma)
}
x
}
gibbs(100)
it returns only 1 or 2 values of theta, some times NaN, this may be if any
theta is greater than 1, which is used as the probability for the next
rbinom(), so returns NaN. Can some one suggest to solve this problem?
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] Error in drawing
I have got following error in drawing wavelet fitting. can some one help? library(faraway) data(lidar) newlidar-lidar[c(1:128),] library(wavethresh) wds - wd(newlidar$logratio) draw(wds) Error in plot.default(x = x, y = zwr, main = main, sub = sub, xlab = xlab, : formal argument type matched by multiple actual arguments [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] GAM
Hi R community! I am analyzing the data set motorins in the package faraway by using the generalized additive model. it shows the following error. Can some one suggest me the right way? library(faraway) data(motorins) motori - motorins[motorins$Zone==1,] library(mgcv) amgam - gam(log(Payment) ~ offset(log(Insured))+ s(as.numeric(Kilometres)) + s(Bonus) + Make + s(Claims),family = gaussian, data = motori) Error in smooth.construct.tp.smooth. spec(object, dk$data, dk$knots) : A term has fewer unique covariate combinations than specified maximum degrees of freedom summary(amgam) Error in summary(amgam) : object 'amgam' not found Gyan [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Correction in error
Hello R community, following is my code and it shows error, can some one
fix this error and explain why this occurs?
gibbs -function(m,n, theta = 0, lambda = 1){
alpha - 1.5
beta - 1.5
gamma - 1.5
x- array(0,c(m+1, 3))
x[1,1] - theta
x[1,2] - lambda
x[1,3]- n
for(t in 2:m+1){
x[t,1] - rbinom(x[t-1,3], 1, x[t-1,1])
x[t,2]-rbeta(m, x[t-1,1] + alpha, tx[t-1,3] - x[t-1,1] + beta)
x[t,3] - rpois(x[t-1,3] - x[t-1,1],(1 - x[t-1,2])*gamma)
}
x
}
gibbs(100, 10)
Gyn
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
