From: Spencer Graves <spencer.graves_at_pdf.com>

Date: Tue 16 Aug 2005 - 01:37:52 EST

>>It might be good to have an example of your problem.

*>>
*

*>>On Aug 4, 2005, at 5:57 AM, Mark Miller wrote:
*

*>>
*

*>>>Hi, I have been trying to figure out how to use the nonlinear
*

*>>>regression to
*

*>>>fit the cumulative lognormal distribution to a number of data
*

*>>>points I have
*

*>>>but I am a new R user and I cant quite decipher the notes on nonlinear
*

*>>>regression. Any help in this regard will be greatly appreciated,
*

*>>>my email
*

*>>>address is mmiller@nassp.uct.ac.za
*

*>>>
*

*>>>______________________________________________
*

*>>>R-help@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@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
*

Date: Tue 16 Aug 2005 - 01:37:52 EST

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:

> 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: >

>>It might be good to have an example of your problem.

-- Spencer Graves, PhD Senior Development Engineer PDF Solutions, Inc. 333 West San Carlos Street Suite 700 San Jose, CA 95110, USA spencer.graves@pdf.com www.pdf.com <http://www.pdf.com> Tel: 408-938-4420 Fax: 408-280-7915 ______________________________________________ R-help@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.htmlReceived on Tue Aug 16 01:53:55 2005

*
This archive was generated by hypermail 2.1.8
: Sun 23 Oct 2005 - 15:20:12 EST
*