Try to search your question here, if you can't find : Ask Any Question Now ?

I am getting different results from my function for even same inputs

HomeCategory: stackoverflowI am getting different results from my function for even same inputs
craig asked 3 weeks ago

in the package below, a function is included. Some parameters are obtained with this function after running. And the code below is an example from package manual. Without any data, this code gives results by considering values of mean, var, covlag1 and pdr. So to try the function, I am opening 2 different R sessions, and copying the same code for each session. But unfortunately, I am getting different results.

My question is: Can someone just try this function by adding the package and copying the function. Because at first I was getting same results for different sessions, but now each session gives different results.

    # Equations of Random Parameter Bartlett-Lewis model
# modelled mean
library(EAS)
meanMBLRPM<-function(a,l,v,k,f,mx,h=1) {
x<-(h*l*mx*v*(1+k/f))/(a-1)
return(x)
}
# modelled variance
varMBLRPM<-function(a,l,v,k,f,mx,h=1) {
A<-(2*l*(1+k/f)*(mx^2)*(v^a))/((f^2)*((f^2)-1)*(a-1)*(a-2)*(a-3))
B<-(2*(f^2)-2+k*f)*(f^2)*((a-3)*h*(v^(2-a))-(v^(3-a))+((v+h)^(3-a)))
C<-k*(f*(a-3)*h*(v^(2-a))-(v^(3-a))+((v+f*h)^(3-a)))
D<-A*(B-C)
return(D)
}
# modelled covariance
covarMBLRPM<-function(a,l,v,k,f,mx,h=1,lag=1) {
A<-(l*(1+k/f)*(mx^2)*(v^a))/((f^2)*((f^2)-1)*(a-1)*(a-2)*(a-3))
B<-(2*(f^2)-2+k*f)*(f^2)*(((v+(lag+1)*h)^(3-a))-2*((v+lag*h)^(3-a))+((v+(lag-1)*h)^(3-a)))
C<-k*(((v+(lag+1)*h*f)^(3-a))-(2*((v+h*lag*f)^(3-a)))+((v+(lag-1)*h*f)^(3-a)))
D<-A*(B-C)
return(D)
}
# modelled probability dry
pdrMBLRPM<-function(a,l,v,k,f,h=1) {
mt<-((1+(f*(k+f))-(0.25*f*(k+f)*(k+4*f))+((f/72)*(k+f)*(4*(k^2)+27*k*f+72*(f^2))))*v)/(f*(a-1))
G00<-((1-k-f+1.5*k*f+(f^2)+0.5*(k^2))*v)/(f*(a-1))
A<-(f+(k*(v/(v+(k+f)*h))^(a-1)))/(f+k)
D<-exp(l*(-h-mt+G00*A))
return(D)
}
# Historical statistics (National Technical University of Athens rain gauge, Athens)
mean1 = 0.1226;var1 = 0.6323;cov1lag1 = 0.3271;pdr1 = 0.9183
mean6 = 0.7358;var6 = 10.1490;cov6lag1 = 4.0773;pdr6 = 0.8251
mean12 = 1.4705;var12 = 29.357;cov12lag1 =7.6865;pdr12 = 0.7476
mean24 = 2.9410;var24 = 76.667;cov24lag1 =10.2886;pdr24 = 0.6238
# Objective function
fopt <- function(x) {
a<-x[1];l<-x[2];v<-x[3];k<-x[4];f<-x[5];mx<-x[6]
w1=1;w2=1;w3=1;w4=1
S1 <- w1*((meanMBLRPM(a,l,v,k,f,mx,h=1)/mean1)-1)^(2)+
w2*((varMBLRPM(a,l,v,k,f,mx,h=1)/var1)-1)^(2)+
w3*((covarMBLRPM(a,l,v,k,f,mx,h=1,lag=1)/cov1lag1)-1)^(2)+
w4*((pdrMBLRPM(a,l,v,k,f,h=1)/pdr1)-1)^(2)
S6 <- w1*((meanMBLRPM(a,l,v,k,f,mx,h=6)/mean6)-1)^(2)+
w2*((varMBLRPM(a,l,v,k,f,mx,h=6)/var6)-1)^(2)+
w3*((covarMBLRPM(a,l,v,k,f,mx,h=6,lag=1)/cov6lag1)-1)^(2)+
w4*((pdrMBLRPM(a,l,v,k,f,h=6)/pdr6)-1)^(2)
S12 <- w1*((meanMBLRPM(a,l,v,k,f,mx,h=12)/mean12)-1)^(2)+
w2*((varMBLRPM(a,l,v,k,f,mx,h=12)/var12)-1)^(2)+
w3*((covarMBLRPM(a,l,v,k,f,mx,h=12,lag=1)/cov12lag1)-1)^(2)+
w4*((pdrMBLRPM(a,l,v,k,f,h=12)/pdr12)-1)^(2)
S24 <- w1*((meanMBLRPM(a,l,v,k,f,mx,h=24)/mean24)-1)^(2)+
w2*((varMBLRPM(a,l,v,k,f,mx,h=24)/var24)-1)^(2)+
w3*((covarMBLRPM(a,l,v,k,f,mx,h=24,lag=1)/cov24lag1)-1)^(2)+
w4*((pdrMBLRPM(a,l,v,k,f,h=24)/pdr24)-1)^(2)
S<-S1+S6+S12+S24
if(is.infinite(S)) {S<-10^8}
if(is.na(S)) {S<-10^8}
return(S)
}
# set the interior and exterior parameters bounds
xmin <- c(1.0001,0.001,0.001,0.001,0.001,0.001)
xmax <- c(5,0.1,5,1,1,20)
xlow <- c(1.0001,0.001,0.001,0.001,0.001,0.001)
xup <- c(15,0.1,20,20,1,50)
# apply the evalutionary annealing-simplex method
modecal <- eas(n=6,m=30,xmin,xmax,xlow,xup,fn=fopt,maxeval=1500,ftol=1.e-07,ratio=0.99,pmut=0.9,
beta=2,maxclimbs=5)
modecal
a<-modecal$bestpar[[1]];l<-modecal$bestpar[[2]];v<-modecal$bestpar[[3]]
k<-modecal$bestpar[[4]];f<-modecal$bestpar[[5]];mx<-modecal$bestpar[[6]]
par <- c(a=a,l=l,v=v,k=k,f=f,mx=mx)
par
# In order to use the derived parameters in the functions of HyetosR
# as well as in the classic version of Hyetos,please be sure that
# for parameters mx and sx the length units are millimeters (mm)
# and for parameters l, v, mx and sx the time units are days (d).
# For this reason, make the following unit conversions:
l<-l*24
v<-v/24
mx<-mx*24
# parameter set for implementation in HyetosR functions
par <- c(a=a,l=l,v=v,k=k,f=f,mx=mx)
par

required package

1 Answers
Best Answer
Jyoti answered 3 weeks ago
Your Answer

11 + 19 =

Popular Tags

WP Facebook Auto Publish Powered By : XYZScripts.com