Using nonlinear regression

2 messages · Mark Miller, Spencer Graves

Original

Mark Miller

Thu, Aug 11, 2005 1:51 AM #

Attached is a copy of my code, the data and the plots obtained by varying 
terms manually, I was told that nonlinear regression in R could find the 
values for me, but I am unable to figure how exactly I could implement this.  
Any help would be very greatly appreciated as I am completely stuck on this 
problem.

On Thursday 04 August 2005 13:40, you wrote:

______________________________________________
R-help at stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-
guide.html

-------------- next part --------------
A non-text attachment was scrubbed...
Name: cumIAT.ps
Type: application/postscript
Size: 47024 bytes
Desc: not available
Url : https://stat.ethz.ch/pipermail/r-help/attachments/20050811/3e2bdc83/cumIAT.ps
-------------- next part --------------
A non-text attachment was scrubbed...
Name: cumTFR.ps
Type: application/postscript
Size: 13981 bytes
Desc: not available
Url : https://stat.ethz.ch/pipermail/r-help/attachments/20050811/3e2bdc83/cumTFR.ps
-------------- next part --------------
A non-text attachment was scrubbed...
Name: cumTTR.ps
Type: application/postscript
Size: 35507 bytes
Desc: not available
Url : https://stat.ethz.ch/pipermail/r-help/attachments/20050811/3e2bdc83/cumTTR.ps
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: Rstuff.R
Url: https://stat.ethz.ch/pipermail/r-help/attachments/20050811/3e2bdc83/Rstuff.pl

4 days later

Mon, Aug 15, 2005 8:37 AM #

Do you want to estimate the parameters of a lognormal distribution or 
learn how to do nonlinear regression in R?

	  If the former, as far as I know, the best known method is maximum 
likelihood, for which the answer is to compute mean and standard 
deviations of the logs.  This assumes you are talking about the 
2-parameter lognormal.  I don't know the best method for a 3-parameter 
lognormal.  If that's what you want, PLEASE do read the posting guide! 
"http://www.R-project.org/posting-guide.html".  Doing so can increase 
the chances of getting a useful reply.

	  If you want examples of nonlinear regression, have you considered 
"nls" and "optim"?

	  spencer graves

Mark Miller wrote:

______________________________________________
R-help at stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-
guide.html


------------------------------------------------------------------------

rm(list=ls())

outPut = read.csv("dataOut2.csv")
arrive = outPut[1]
register = outPut[2]
complete = outPut[3]


IAT = 0
TTR = 0
TTC = 0
TFR = 0
cnt = 1

for(i in array(2:dim(arrive)[1]))
{
temp = outPut[i,3]-outPut[i,2]
if(temp > 0)
{
	IAT[cnt] = outPut[i,1]-outPut[i-1,1]
	TTR[cnt] = outPut[i,2]-outPut[i,1]
	TFR[cnt] = outPut[i,3]-outPut[i,2]
	cnt = cnt + 1
}
}

cumIAT = IAT/sum(IAT)
for(i in array(2:length(IAT)))
{
cumIAT[i] = cumIAT[i-1]+cumIAT[i]
}
cumIAT[1] = 0
plot(cumIAT,do.point=F)


TTR[cnt] = outPut[1,2]-outPut[1,1]
TFR[cnt] = outPut[1,3]-outPut[1,2]

# Plot for inter-arrival times
x = seq(0,30,0.01)
#postscript("cumIAT.ps")
plot(ecdf(IAT), do.point=FALSE)
lines(x, pexp(x,0.4))
dev.off()
# rexp(100,0.21)

x = seq(0,20,0.01)
postscript("cumTTR.ps")
plot(ecdf(TTR), do.point=FALSE)
lines(x, plnorm(x,1,0.7))
dev.off()
# rlnorm(100,1,0.7)

# Plot for Time to complete from registered
x = seq(0,30,0.01)
postscript("cumTFR.ps")
plot(ecdf(TFR), do.point=FALSE)
lines(x*600, pbeta(x,1.4,4.3))
dev.off()
# rbeta(100,1.6,5)*600

# Find the position with the leat time and hence the next avaliable ambulance
minimum = function(toFind)
{
min = 0;
pos = 0;
for(i in array(1:length(toFind)))
{
	if(i == 1)
	{
		min = toFind[i]
		pos = i
	}
	else
	{
		if(toFind[i] < min)
		{
			min = toFind[i]
			pos = i
		}
	}
} 

pos
}

ambsReq = 0
numAmbs = 0
numberAmbs = 0
avgWait = 1
numberAmbs2 = 0
avgWaitTime2 = 0
avgWaitTime = 0
counter = 0
counter2 = 1
cntO = 1

for(i in array(1:50))
{
while(avgWait > 0)
{
	counter = counter + 1
	numAmbs = numAmbs + 1
	numberAmbs[counter] = numAmbs
	numCalls = 1
	ambs = array(c(array(0,numAmbs)), dim=c(numAmbs,numCalls))
	waitTime = ambs
	totalTime = ambs
	currTime = 0
	timeTS = 0
	IotherAT = 0
	TotherTR = 0

	for(i in array(1:500))
	{
		#interAT = IAT(ceil(rand()*length(IAT)));
		interAT = rexp(1,0.21)
		#timeTR = TTR(ceil(rand()*length(TTR)));
		timeTR = rlnorm(1,1,0.7)
		#timeFR = TFR(ceil(rand()*length(TFR)));
		timeFR = rbeta(1,1.4,4.3)*600

		IotherAT[i] = interAT
		TotherTR[i] = timeTR

		currTime = currTime + interAT
	
		pos = minimum(totalTime)

		if(ambs[pos,numCalls] != 0)
		{
			numCalls = numCalls + 1
			ambs = array(c(ambs,array(0,numAmbs)), dim=c(numAmbs,numCalls))
			waitTime = array(c(waitTime,array(0,numAmbs)), dim=c(numAmbs,numCalls))
			if(totalTime[pos] > currTime)
				waitTime[pos,numCalls] = totalTime[pos] - currTime
			totalTime[pos] = totalTime[pos] + interAT + timeTR + timeFR
			ambs[pos,numCalls] = interAT + timeTR + timeFR
		}
		else
		{
			for(i in array(1:numCalls))
			{
				if(ambs[pos,i] == 0)
				{
					if(totalTime[pos] > currTime)
						waitTime[pos,i] = totalTime[pos] - currTime
					totalTime[pos] = totalTime[pos] + interAT + timeTR + timeFR
					ambs[pos,numCalls] = interAT + timeTR + timeFR
					break
				}
			}
		}
	}

	avgWait = sum(waitTime)/500
	avgWaitTime[counter] = avgWait
	if(numAmbs == 25)
	{
		avgWaitTime2[cntO] = avgWait
		numberAmbs2[cntO] = numAmbs
		cntO = cntO + 1
	}
}

postscript("timeAmbs.ps")
plot(numberAmbs,avgWaitTime,'lines') 
dev.off()

postscript("timeAmbs2.ps")
plot(numberAmbs2,avgWaitTime2,'lines') 
dev.off()

ambsReq[counter2] = numAmbs
numAmbs = 0
numberAmbs = 0
avgWait = 1
avgWaitTime = 0
counter = 0
counter2 = counter2 +1
cntO = 1
numberAmbs2 = 0
avgWaitTime2 = 1
}




TotherFR = 0
cnt = 1
for(i in array(1:dim(ambs)[1]))
{
for(j in array(1:dim(ambs)[2]))
{
	if(ambs[i,j] != 0)
	{
		TotherFR[cnt] = ambs[i,j]
		cnt = cnt +1
	}
}
}

postscript("dataIAT.ps")
hist(IAT,30)
dev.off()

postscript("simIAT.ps")
hist(IotherAT,30)
dev.off()

postscript("dataTTR.ps")
hist(TTR,30)
dev.off()

postscript("simTTR.ps")
hist(TotherTR,30)
dev.off()

postscript("dataTFR.ps")
hist(TFR,30)
dev.off()

postscript("simTFR.ps")
hist(TotherFR,30)
dev.off()


------------------------------------------------------------------------

______________________________________________
R-help at stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

Spencer Graves, PhD
Senior Development Engineer
PDF Solutions, Inc.
333 West San Carlos Street Suite 700
San Jose, CA 95110, USA

spencer.graves at pdf.com
www.pdf.com <http://www.pdf.com>
Tel:  408-938-4420
Fax: 408-280-7915